Anomaly calculation with groupby leaves seasonal cycle

What happened?

When calculating EOFS with sea surface temperature (SST) anomaly, We see the seasonal cycle still in the data.
Tested with xarray v2023.03.0 and v202304.02, both have the same problem (seasonal cycle still in the data).
So the anomaly calculated with groupby is not correct.

Attached is the results of using xarray v2023.04.02

What did you expect to happen?

The EOFs will not have seasonal cycle. Attached is the results from xarray v0.20.2

Minimal Complete Verifiable Example

import xarray as xr
import numpy as np
import pandas as pd
from eofs.xarray import Eof

url = 'http://kage.ldeo.columbia.edu:81/SOURCES/.LOCAL/.sst.mon.mean.nc/.sst/dods'
ds = xr.open_dataset(url)

gb = ds.sst.groupby('time.month')
ds_anom = (gb - gb.mean(dim='time')).compute()
solver = Eof(ds_anom)
pcs = solver.pcs(npcs=4,pcscaling=1)
eofs = solver.eofsAsCorrelation(neofs=4)

eofs.plot(x='lon',y='lat',col='mode',col_wrap=4,add_colorbar=1)
pcs.plot(x='time',col='mode',col_wrap=4);

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

No response

Anything else we need to know?

No response

Environment

The issue names no repository file or test. Start by running the provided Python MVCE and inspecting the groupby anomaly calculation and its result across the reported xarray versions. Done means the anomaly no longer retains the seasonal cycle and the behavior is covered by a reproducible test.