# Galaxy Cluster Statistics in Modified Gravity Cosmologies

Marcell Howard  Arthur Kosowsky  Department of Physics and Astronomy,
University of Pittsburgh, Pittsburgh, PA 15260 USA  Georgios Valogiannis
Department of Physics, Harvard University, Cambridge, MA 02138, USA

###### Abstract

If the accelerated expansion of the universe is due to a modification of
general relativity at late times, it is likely that the growth of structure on
large scales would also display deviations from the standard cosmology. We
investigate the statistics of the distribution of galaxy cluster-sized halos
as a probe of gravity. We analyze the output of several matched N-body
simulations with the same initial conditions and expansion histories but using
both DGP and $f(R)$ gravity with various parameters. From each simulation we
extract the cluster mass function, power spectrum, and mean pairwise velocity
at redshifts $1$, $0.3$, and $0$. All three statistics display systematic
differences between gravity theories. The mean pairwise velocity provides an
important consistency test for any posited departure from general relativity
suggested by measurements of the power spectrum and cluster mass function.
Upcoming microwave background experiments, including Simons Observatory and
CMB-S4, will detect tens of thousands of galaxy clusters via the thermal
Sunyaev-Zeldovich effect, probe their masses with lensing of the microwave
background, and potentially measure velocities using the transverse lensing
effect or the kinematic Sunyaev-Zeldovich effect. These cluster measurements
promise to be a substantial probe of modified gravity.

##  I Introduction

The most common explanation for the current accelerating expansion of the
universe is some undetected stress-energy component with negative pressure and
dominant energy density. However, no compelling ideas for this component have
emerged. A simple dimensional analysis argument for the natural size of any
vacuum energy component combines the fundamental constants $c$, $\hbar$, and
$G$ to obtain an energy density which is too large by the infamous 120 orders
of magnitude. Dynamical components such as a scalar field face similar
problems. Within the framework of general relativity, accelerating cosmic
expansion today requires new physics (”Dark Energy”) which dominates the
energy budget of the universe with a characteristic energy scale corresponding
roughly to room temperatures – a scale which is constrained by innumerable
experimental results.

This conceptual difficulty has led to consideration of the only other logical
alternative: that general relativity (GR) is not the correct theory to
describe cosmological dynamics (for a comprehensive review, see Clifton et al.
(2012)). general relativity has passed numerous stringent observational tests,
but all of these are on length scales corresponding to the size of our solar
system or smaller. It is therefore a possibility that a fundamental theory of
gravity mimics general relativity on small scales, while having a different
form on scales much larger than the solar system. Constructing such theories
is challenging, because any modification must be made small enough to be
consistent with current measured values of the parameterized post-Newtonian
(PPN) parameters Will (2014). The propagation speed of gravitational waves
being that of light to within one part in $10^{16}$ Abbott et al. (2016, 2017)
also rules out many possible alternate gravity theories Langlois et al.
(2018).

Gravity has two cosmological effects: it determines the expansion history of
the universe given its energy-momentum content, and it drives the growth of
primordial perturbations into cosmological structures. The combination of
these effects, along with initial conditions visible in the anisotropies of
the cosmic microwave background, determines the detailed distribution of
density perturbations at late times on large scales. Any modification of GR on
large scales therefore will, to some extent, affect the distribution of matter
on large scales.

Here we investigate the impact of differing gravitational theories on the
largest bound masses in the universe, corresponding to galaxy clusters. Galaxy
clusters are rare objects which form in locations corresponding to rare peaks
in the initial mass distribution. The mass of a given individual halo is
determined largely by gravitational growth moderated by the cosmic expansion.
At a finer level of detail, tidal interactions and mergers with other nearby
halos also contribute to the present halo mass. (Radiation feedback from the
formation of stars and active galactic nuclei can also have an effect, but the
largest galaxy clusters have a deep enough gravitational potential that the
effect appears to be minor Bocquet et al. (2015); we ignore all baryonic
effects in this work.) Both the growth factor and the expansion history can be
altered in modified gravity compared to general relativity, as can merger
histories and tidal effects.

While it is not possible to point to properties of an individual cluster as
reflecting modified gravity, the statistical distribution of observable
cluster properties reflects the underlying theory of gravity. In this work, we
analyze a number of matched cosmological N-body simulations with identical
initial conditions and differing gravitational physics to quantify the
statistical differences in galaxy cluster populations due to the nature of
gravity. We first demonstrate the close correspondence between the largest-
mass halos in all simulations, which arises due to the identical initial
conditions. We then extract the distribution of large halo masses, the halo
spatial power spectrum, and the mean pairwise halo velocity for each
simulation and compare with those from the standard cosmology with general
relativity. Systematic differences ranging from a few percent to tens of
percent are evident in the simulations considered.

Observationally, galaxy cluster statistics are a promising route to
constraints on gravity. In recent years, cosmic microwave background
experiments including the Planck Satellite Ade et al. (2016), the Atacama
Cosmology Telescope Ade et al. (2019), and the South Pole Telescope Bleem et
al. (2020) have mapped large sky regions at angular resolution and sensitivity
sufficient to detect thousands of galaxy clusters via the thermal Sunyaev-
Zeldovich (tSZ) effect. The tSZ effect, which creates a spectral distortion in
the microwave radiation passing through the galaxy cluster’s ionized gas, is
essentially independent of cluster redshift and gives a relatively simple
selection function for clusters Bleem et al. (2020); Hilton et al. (2021).
Once identified, cluster positions can be determined in three dimensions
through follow-up galaxy redshift observations; cluster masses can be measured
using optical weak gravitational lensing, such as Murata et al. (2019), or
statistically using lensing of the microwave background Madhavacheril et al.
(2015); and cluster velocities can potentially be obtained from the kinematic
Sunyaev-Zeldovich effect Calafut et al. (2021); Schaan et al. (2021) or
potentially from the transverse lensing effect Hotinli et al. (2019). As we
anticipate a new generation of microwave experiments with steadily increasing
sensitivity, quantifying the possibilities for constraining modified gravity
is of particular interest. These results can also be used to calibrate
approximate calculations of structure growth in modified gravity theories such
as the mass peak-patch formalism Stein et al. (2018, 2020).

The paper is organized as follows: in section 2, we briefly review the Hu-
Sawacki Hu and Sawicki (2007) and (normal-)DGP Dvali et al. (2000) modified
gravity models and their effects on cosmological scales. Section 3 summarizes
the cosmological simulations on which our analysis is based, and establishes
the correspondence of massive halos between the different simulations.
Sections 4, 5, and 6 present the cluster halo mass function, cluster power
spectrum, and halo mean pairwise velocity, and compare between the different
models. A summary and outlook section discusses extending these results to
more general parameterizations of modified gravity, the potential of upcoming
data to measure these cluster statistics well enough to constrain modified
gravity, and the relation of cluster statistics to other probes of large-scale
structure.

##  II Modified Gravity Models

###  II.1 Hu-Swacki $f(R)$ Gravity

This work focuses on scalar-tensor theories of gravity, with extra spin-zero
degrees of freedom added to the Einstein-Hilbert gravitational action. In
particular, we consider two definite models. The first is $f(R)$-gravity,
where the Ricci curvature scalar in the action is replaced with a general
function of the curvature scalar. In particular, the Hu-Sawicki model Hu and
Sawicki (2007) is given by (in units of $c=1$)

$S=\frac{1}{2\kappa^{2}}\int\differential[4]{x}\sqrt{-g}(R+f(R))+S_{m},$ (1)

where $S_{m}$ is the action for any additional matter fields, the reduced
four-dimensional Planck mass is $\kappa^{-1}\equiv\frac{1}{\sqrt{8\pi
G}}\approx$2.43\text{\times}{10}^{18}\text{\,}\mathrm{GeV}$$, and

$f(R)\equiv-\mu^{2}\frac{c_{1}(R/\mu^{2})^{n}}{1+c_{2}(R/\mu^{2})^{n}},$ (2)

where $\mu^{2}\equiv\Omega_{m}H_{0}^{2}$ with $\Omega_{m}$ being the present
time fractional energy density for matter and $H_{0}$ is the value of the
Hubble parameter today. In its original inception, the model was defined by
three parameters $c_{1}$, $c_{2}$, and $n$. At large curvature such that
$\absolutevalue{\mu^{2}/R}\ll 1$, we get

$f(R)\approx-\frac{c_{1}}{c_{2}}\mu^{2}+\frac{c_{1}}{c_{2}^{2}}\mu^{2}\quantity(\frac{\mu^{2}}{R})^{n}.$
(3)

We take $n=1$, a simple special case which captures the qualitative behavior
of the model. To reproduce accelerated expansion, the ratio $c_{1}/c_{2}$ is
set to recover the effective cosmological constant in the high-curvature
limit:

$f(R)\approx-6\Omega_{\Lambda}H_{0}^{2}+|f_{R0}|\frac{\bar{R}^{2}}{R},$ (4)

where $\bar{R}$ is the background curvature value evaluated in $\Lambda$CDM
today

$\bar{R}=3\Omega_{m}H_{0}^{2}\quantity(1+4\frac{\Omega_{\Lambda}}{\Omega_{m}}),$
(5)

$f_{R0}$ is given by Valogiannis and Bean (2019)

${f}_{R0}\equiv\evaluated{\derivative{f}{R}}_{R={\bar{R}}}=-\frac{c_{1}}{c_{2}^{2}}\quantity(\frac{\Omega_{m}}{3(\Omega_{m}+4\Omega_{\Lambda})})^{2},$
(6)

and $f_{R}\equiv\derivative{f}{R}$ functions as a new scalar degree of freedom
termed the ”scalaron”. Notice in the limit of $c_{1}/c_{2}^{2}\rightarrow 0$,
we recover the standard $\Lambda$CDM cosmology, meaning we can parameterize
deviations from GR solely by varying the magnitude of $f_{R0}$.

We consider three specific models by choosing
$\absolutevalue{{f}_{R0}}=\\{10^{-6},10^{-5},10^{-4}\\}$, which will be
referred to as F6, F5, and F4 respectively (which we have ordered by
increasing deviation from GR). The Hu-Sawicki theory is an example of a class
of modified gravity theories that possesses a chameleon mechanism or a
mechanism which screens the effects of modified gravity in regions where the
density exceeds a characteristic threshold. This is equivalent to imparting a
spatial dependence on the mass of the scalaron i.e.

$m^{2}(x)=\frac{1}{3}\quantity(\frac{1+f_{R}}{f_{RR}}-R)\approx\frac{1}{3f_{RR}}\equiv\frac{1}{3}\absolutevalue{\derivative[2]{f}{R}}^{-1},$
(7)

where we took the high curvature limit. As a result, a fifth force with a
Yukawa-like potential is introduced which acts on all massive particles Alam
et al. (2021). The modified Poisson’s Equation (in the quasi-static limit) is

$2\nabla^{2}\Phi=\kappa^{2}a^{2}\variation{\rho_{m}}+\nabla^{2}f_{R},$ (8)

where $\Phi$ is the gravitational potential in the conformal Newtonian gauge,
$\variation{\rho}_{m}$ is the perturbation in the matter energy density and in
the $f(R)\rightarrow 0$ limit, we get the usual unperturbed Poisson’s
Equation. The perturbed field satisfies the constraint

$-3\frac{\nabla^{2}f_{R}}{a^{2}}=\kappa^{2}\variation{\rho}_{m}+\variation{R},$
(9)

where $\variation{R}$ is the perturbed Ricci scalar and the dynamics of this
scalar results in the modified Friedman Equation

$H^{2}+\frac{1}{6}f-(\dot{H}+H^{2})f_{R}+6H(\ddot{H}+4H\dot{H})f_{RR}=\frac{\mu^{2}}{3},$
(10)

where overhead dots denote derivatives with respect to time. The new terms
involving $f$ and its derivatives are negligible compared to the Hubble
parameter $H^{2}$, which results in the background expansion to be
approximately the same as it is in $\Lambda$CDM.

###  II.2 DGP Gravity

The second model we consider is the normal DGP (or nDGP as its more commonly
known) gravity Dvali et al. (2000); Lombriser et al. (2009) which posits that
our four dimensional spacetime is embedded within a five dimensional brane.
The action for the theory is

$S=\int\differential[4]{x}\sqrt{-g}\quantity(\frac{1}{2\kappa^{2}}R+\mathcal{L}_{m})+\frac{1}{2\kappa^{2}r_{c}}\int\differential[5]{x}\sqrt{-g_{5}}R_{5},$
(11)

where $r_{c}$ (called the Vainshtein radius) is the ratio between the 5D and
4D Planck mass Lombriser et al. (2009) which acts as a free parameter that
determines the length-scale for which we recover standard $\Lambda$CDM
cosmology, and $R_{5}$ and $g_{5}$ are the five-dimensional curvature scalar
and metric determinant respectively. The equations of motion are

$2\nabla^{2}\Phi=\kappa^{2}a^{2}\variation{\rho}_{m}+\nabla^{2}\varphi,$ (12)
$\nabla^{2}\varphi+\frac{r_{c}^{2}}{3a^{2}\beta(a)}\quantity[(\nabla^{2}\varphi)^{2}-(\nabla_{i}\nabla_{j}\varphi)^{2}]=\frac{\kappa^{2}a^{2}}{3\beta(a)}\variation{\rho}_{m},$
(13)

where $\Phi$ is the gravitational potential in the conformal Newtonian gauge
and $\beta(a)$ is

$\beta(a)=1+2Hr_{c}\quantity(1+\frac{\dot{H}}{3H^{2}}).$ (14)

In the limit $\varphi\rightarrow 0$ (equivalent to taking the
$r_{c}\rightarrow\infty$ limit in $\beta(a)$), we recover the standard
Poisson’s Equation for the Newtonian gravitational potential. The modified
Friedman equation is given by Valogiannis and Bean (2019)

$H^{2}+\frac{1}{r_{c}}H=\frac{\kappa^{2}}{3}\rho_{\rm tot},$ (15)

where $\rho_{tot}$ contains the contributions from radiation, matter, and the
cosmological constant. In the $r_{c}\rightarrow\infty$ limit, we recover the
usual Friedman Equation. The extra term has the effect of slowing the
expansion and thus we must put in by hand the cosmological constant in order
to replicate late-time acceleration Lombriser et al. (2009). We set
$H_{0}r_{c}=1$ or $H_{0}r_{c}=5$ and refer to these two models as N1 and N5
respectively (here N1 is the stronger deviations from GR and N5 is the weaker
deviation). nDGP is an example of a theory that exhibits the Vainshtein
mechanism i.e. the mechanism in which the effects of modified gravity are
screened at scales smaller than a characteristic length scale (the Vainshtein
radius).

We emphasize here that because modifications to GR are considered small, they
are treated as perturbations on top of the background expansion.

Figure 1: The number of halos in that ratio bin versus The ratio of the mass
between a modified gravity halo and its corresponding matched pair in GR (at
$z=1$) under the constraint that the MG halo mass is at least 20% of the GR
halo’s mass. Figure 2: On the horizontal axis is the distance that the
modified gravity halo has traveled from its matched halo GR pair. The vertical
axis is the number of halos in a given separation bin. Figure 3: The magnitude
of the relative velocity between a modified gravity halo and its matched GR
pair for the same set as halos as in Figure 1.

##  III The _ELEPHANT_ Simulations and _ECOSMOG_ & _ECOSMOG-V_

We extract statistical properties of cluster-mass halos from the set of
_ELEPHANT_ (Extended LEnsing PHysics using ANaltyic ray Tracing) simulations
of cosmic structure growth with modified gravity Cautun et al. (2018). These
cosmological N-body simulations use the _ECOSMOG_ (Efficient COde for
Simulating MOdified Gravity) Li et al. (2012) for $f(R)$ gravity, and the
_ECOSMOG-V_ Barreira et al. (2015) code for nDGP gravity. Both of these codes
are modifications of the N-body and hydrodynamic simulation code _RAMSES_
Teyssier (2002). The _ELEPHANT_ simulations have been used previously for Alam
et al. (2021). The _ECOSMOG_ algorithm is adapted specifically to solve the
equations of motion for an $f(R)$ theory, particularly theories that
incorporate a chameleon mechanism whereas _ECOSMOG-V_ is particularly well
suited for MG theories that have a Vainshtein mechanism (the V in _ECOSMO-V_
stands for Vainstein).

These codes differ from the standard cosmological N-body code in two ways: (1)
they evolve density perturbations using a modified Poisson equation reflecting
the modified gravity effective force law and (2) they incorporate a screening
mechanism to mitigate the effect of modified gravity on small scales.

Each simulation is contained in a box with side length $L_{\rm Box}=1024$
Mpc/h and periodic boundary conditions, with $N=1024^{3}$ particles of mass
$7.78\times 10^{10}h^{-1}M_{\odot}$. Each simulation began with the same
initial conditions, given by a specific realization of a Gaussian random
density field with a power spectrum consistent with the cosmic microwave
background anisotropies at early times. Density perturbations are assumed to
be adiabatic with equal fractional density perturbations in all components at
every point in space. The density power spectrum at the starting redshift
$z_{i}=49$ can be obtained from standard codes such as _CLASS_ Blas et al.
(2011), assuming evolution of the initial density field at recombination using
the usual linear growth factor from general relativity. In other words, the
simulations assume that any modified gravity effects are only relevant at
redshifts below $z=z_{i}$, and general relativity is a good approximation at
earlier times. We make use of code outputs at redshifts $z\in\\{1,0.3,0\\}$.
We also made use of only one realization (Box 1) out of the five that were
available.

Halo catalogs were compiled from the simulations using the halo finder
_ROCKSTAR_ (Robust Overdensity Calculation using K-Space Topologically
Adaptive Refinement) Behroozi et al. (2012). _ROCKSTAR_ is a 6D phase-space
algorithm to locate dark matter halos, their substructure, and non-spherical
distortions. Halo masses are defined as

$M_{200c}=\frac{4\pi}{3}\rho_{c}R^{3}_{200c},$ (16)

i.e. the amount of mass contained within a sphere of radius $R_{200c}$
enclosing 200 times the critical density. The lowest mass objects in these
catalogues are $\order{10^{11}}M_{\odot}/h$.

The cosmological parameters used for these simulations were the best-fit WMAP9
parameters Hinshaw et al. (2013): $\Omega_{b}=0.046$, $\Omega_{\rm
CDM}=0.235$, $\Omega_{\Lambda}=0.719$, $h=0.697$, $n_{s}=0.971$, and
$\sigma_{8}=0.820$. The _ROCKSTAR_ catalogue includes masses, co-moving
coordinates, and peculiar velocities which we use in this analysis.

##  IV Halo Matching

In order to characterize the differences between GR and a modified gravity
model, we first find corresponding objects in the two simulations. First, we
collect all GR halos with a mass of $M>10^{14}M_{\odot}/h$ ($N=18138$). We
focus on this mass range because upcoming high-resolution CMB experiments will
plausibly detect individual clusters in this mass range via the thermal SZ
effect. The most massive MG halo within a co-moving distance of 2 Mpc/h from
the GR halo is considered the match of the GR halo, if the MG halo mass is
within a factor of 5. In the simulations, every GR halo had at least one such
matching in each MG simulation. A small fraction of the GR halos (759 out of
18138, or 4.18%) had more than one matching halo in at least one of the MG
simulations. In these cases, determining the MG match in one or more
simulations could prove ambiguous, so for simplicity we simply exclude these
GR halos from our analysis.

All of the matched pairs have a separation less than 2 Mpc/h. This makes
sense, because typically a cluster will have a characteristic peculiar
velocity of a few hundred km/s (say 300 km/s), and over a cosmological time
scale of a few Gyr (say 10 Gyr), it would move by 3 Mpc. The difference in
peculiar velocities between MG and GR will generally be significantly smaller
than the GR peculiar velocity, provided that the modification to gravity is
small.

To compare the matched pairs, we looked at three different statistics: the
ratio between a halo mass as produced within a MG model to a halo produced
within GR, the magnitude of the relative position (which we call the
separation) between a MG halo and a GR halo, and the magnitude of the relative
velocity between the MG halo and the GR halo.

Figure 1 shows the distribution of the ratio between the mass of a MG halo and
its matched pair GR halo at a redshift of $z=1$. While each distribution is
peaked near a mass ratio of 1, each MG simulation has some matched halos with
substantial mass differences, up to a factor of 5. Since the initial
conditions are the same for each model, the discrepant masses reflect physical
differences in the evolution of the halos. These likely include differing
amounts of subhalos being accreted, and enhanced or depressed tidal
disruption. Figure 2 is the distribution of halo separations between a MG halo
and its GR matched pair. MG models that more strongly deviate from GR have a
larger spread of pair separations. For models further from GR, the separation
distribution also becomes somewhat asymmetric around the mean. Figure 3
displays the distribution of the relative velocities between a MG halo and its
matched pair GR halo, which appear generally consistent with the spread of
displacements. All of these distributions tend to be slightly skewed towards
larger differences, because the modifications of gravity considered all tend
to increase the strength of gravity at late times.

##  V Halo Mass Function

The simplest statistic characterizing the cluster halo distribution is the
halo mass function, counting the number of halos in various mass bins. Figure
4 shows the cluster mass function for halos with masses
$M_{200c}>10^{14}M_{\odot}/h$ at three different redshifts $z=1$, $z=0.3$, and
$z=0$, while Fig. 5 replots these data showing the ratio of the number of
clusters in each modified theory mass bin to the number in the corresponding
GR mass bin. Poisson error bars for the number of clusters in each mass bin
are shown; these characterize how well the given simulations determine the
mass function in each bin. At any redshift, higher-mass bins contain fewer
clusters and thus have larger uncertainties in the signal. The number of halos
in the largest mass bins increases significantly between $z=1$ and $z=0.3$,
but the growth of massive halos is considerably reduced between $z=0.3$ and
$z=0$ due to the accelerating expansion of the universe suppressing halo
growth compared to the matter-dominated era.

The models that have the largest increase in gravitational strength from GR
(i.e. F4 and N1) produce larger numbers of massive halos, as expected. Similar
results were shown in Ref. Alam et al. (2021) using the same simulations. This
effect appears consistent across redshifts. Halos at the lower mass limit of
$10^{14}M_{\odot}/h$ show significantly less departure from GR than higher-
mass halos, as is clear in Fig. 5; this is evidence that the chameleon or
Vainshtein screening mechanisms in these models are having significant impact
at the lower end of the cluster mass range.

Cluster halo mass correlates with the easily measured thermal SZ distortion,
by which the majority of clusters have been detected. The thermal SZ has a
relatively simple selection function which is approximately mass limited.
Probing small differences in the mass function will likely require direct mass
estimates from either optical weak lensing or CMB lensing. Of the statistics
considered here, the mass function has the largest signal but possibly the
most complicated systematics to account for.

Figure 4: The halo mass function for each gravity model at each of the given
redshifts. Figure 5: Ratio of the number of halos in each mass bin for a
modified gravity model and general relativity.

##  VI Power Spectrum

Figure 6: The power spectrum for each gravity model in order of decreasing
redshift plotted over halos with $M_{200c}\geq 10^{14}M_{\odot}/h$, for
redshifts $z=1$ (top), $z=0.3$ (middle), $z=0$ (bottom). All spectra were
computed in the publicly available code _nbodykit_ Hand et al. (2018).

Figure 7: The ratio of power spectra between a modified gravity model and
general relativity for the same mass range as prior, for redshifts $z=1$
(top), $z=0.3$ (middle), $z=0$ (bottom).

Next we consider the spatial matter power spectrum of the cluster halo
distribution. All spectra were computed with the publicly available code
_nbodykit_ Hand et al. (2018). Figure 6 shows the power spectra for our models
at each of the considered redshifts. The overall shape of the spectrum remains
consistent between the different models and redshift, varying mostly in
amplitude. Figure 7 replots the ratio of each MG power spectrum to the
corresponding GR spectrum, for each different redshift.

Figure 7 shows that the power spectrum of each model is broadly consistent
with an overall change in amplitude compared to GR, but little change in
shape. The stronger $f(R)$ models display an overall reduction in power
spectrum amplitude at all redshifts. The nDGP models, in contrast, show a
suppressed power spectrum at high redshift but return close to GR at $z=0$.
Additionally, the nDGP models show a larger power-spectrum scatter at high $k$
than the $f(R)$ models. However, while there was little change in shape were
we to calculate the fractional difference, $1-P_{MG}(k)/P_{GR}(k)$, with
catalogues containing more galaxies in the chosen mass range, which would thus
give less noisy curves, some change in the shape would be found.

Our results, with a cluster mass cutoff, are similar to those in Figure 4 of
Alam et al. (2021) (which uses a smaller mass cutoff). Additionally, because
large halos tend to cluster more, the amplitude is greater than that of the
full matter power spectrum. Even though MG generically boosts the strength of
gravity, we see a reduction of power in the case of the $f(R)$ gravity models.
While this result is unintuitive, it is consistent with the explanation in
Ref. Hernández-Aguayo et al. (2018): (1) Since the enhanced strength of
gravity is environmentally sensitive, whenever the chameleon screening is
efficient we see more halos form in low-density regions compared to GR. This
reduced the difference in the density of large halos between high and low-
density regions, leading to less power on the scales of large overdensities
and voids. (2) Stronger gravity overall leads initially to more halos at a
given mass, but that results in a higher merger rate. The combination of these
also will modify the power spectrum, but the direction depends on the detailed
balance between halo formation and merger rate. However, both of these effects
are due to the imposition of a halo mass cutoff. If we were to consider the
full matter distribution, we would recover the enhanced matter power spectrum
expected if the strength of gravity has been boosted at all scales (cf. Fig. 1
of Ref. Alam et al. (2021)).

This leads us to the conclusion that the relationship between the cluster halo
power spectrum and the underlying gravity model is complex, and care must be
taken in the interpretation of measurements.

##  VII Mean Pairwise Velocity

Figure 8: Comparison of the mean pairwise velocity as a function of co-moving
separation, in separation bins of 10 Mpc/h, for all gravity models, for
redshifts $z=1$ (left), $z=0.3$ (middle), $z=0$ (right). The standard error on
the mean is displayed on each point. Figure 9: The difference between the mean
pairwise velocity for a modified gravity model and GR, in separation bins of
10 Mpc, for redshifts $z=1$ (left), $z=0.3$ (middle), $z=0$ (right).The
standard error for each point combines the GR and MG errors in quadrature.

We also consider the cluster mean pairwise peculiar velocity, as future
microwave background maps may provide a path to a clean measurement of this
quantity. The mean pairwise peculiar velocity (MPV) $v(r)$ Davis and Peebles
(1977) is the mean of the relative velocity component along the line of
separation between each possible pair of cluster halos,

$v(r)=\langle\hat{\mathbf{r}}_{ij}\cdot(\overrightarrow{\mathbf{v}}_{i}-\overrightarrow{\mathbf{v}}_{j})\rangle,$
(17)

where
$r=\absolutevalue{\overrightarrow{\mathbf{r}}_{i}-\overrightarrow{\mathbf{r}}_{j}}$
is the magnitude of the co-moving separation between the i-th and j-th
cluster, $\hat{\mathbf{r}}_{ij}$ given by

$\hat{\mathbf{r}}_{ij}=\frac{\overrightarrow{\mathbf{r}}_{i}-\overrightarrow{\mathbf{r}}_{j}}{|\overrightarrow{\mathbf{r}}_{i}-\overrightarrow{\mathbf{r}}_{j}|},$
(18)

is the unit vector parallel to the line that connects the i-th and j-th halo,
and $\overrightarrow{\mathbf{v}}_{i}$ is the velocity of the i-th halo and
$\langle\cdots\rangle$ is mean over all possible pairs $i,j$ of clusters. If
the distribution of cluster velocities were random, we would anticipate the
mean to be zero (clusters would be just as likely to fall to move toward each
other as they are to move outward). However, since gravity is an attractive
force, we expect to find, on average, a small negative mean velocity for all
pair of clusters at a given separation. This makes MPV a promising tool to
test GR on scales of tens of Mpc because it directly probes the effect of of
gravity on these scales.

Figure 8 shows the MPV at each redshift as a function of the co-moving
separation between two halos. Normal errors on the sums over pairs are shown.
The velocities are all negative, as expected. The models which boost the
strength of gravity at small separations tend to show more negative
velocities. This can more clearly be seen in Fig. 9 which shows the difference
between the MPV in each gravity model and GR. The various models show a range
of behaviors. The two nDGP models show consistently higher MPV than general
relativity at all separations, while the $f(R)$ models have a less consistent
behavior with separation. This is likely due to screening in $f(R)$ depending
on the local density, inducing more complex velocity variations.

Another notable feature is that the MPV decreases with decreasing redshift,
due to the effect of late-time acceleration on large-scale structure
formation. The differences between models also becomes less pronounced at late
times, due to the damping of the velocities by the cosmic acceleration common
to all the models. Compared to the power spectrum, the regularity of the mean
pairwise velocity and its distinctive evolution with redshift make it a
promising diagnostic for any departures from general relativity.

##  VIII Summary, Conclusion, and Outlook

We have analyzed the output of several modified gravity simulations of large-
scale structure to extract statistical properties of galaxy clusters. Starting
from identical initial conditions, the simulations host corresponding galaxy
clusters to the standard gravity case, with some scatter in cluster mass and
displacements in cluster positions and velocities which we quantify. We then
extract the cluster mass function, power spectrum, and mean pairwise velocity
for each simulation. For the simulations we consider, the differences in these
statistics range from a few percent to tens of percent. If gravity is
modified, differences in these cluster statistics will be correlated,
providing leverage in distinguishing modified gravity from systematic errors.

The models considered, of the DGP or $f(R)$ form, are a range of
representative modified gravity models. While a cosmological gravity theory
may be different in detail, we expect cluster statistics to be roughly within
the range of models considered here. Greater effects would likely have already
been spotted, while it is plausible that a modification of gravity substantial
enough to generate accelerating expansion will have some non-negligible effect
on galaxy cluster scales. Any cosmological modification of gravity must be
shielded from significant effects at solar system scales where departures from
general relativity are strongly constrained observationally (e.g., Vikram et
al. (2013)).

Modified gravity constraints on scales of galaxy clusters or larger scales can
provide a powerful complement to the usual post-Newtonian constraints on solar
system scales. Here we have evaluated three galaxy cluster statistics – mass
function, spatial power spectrum, and mean pairwise velocity – which show a
range of variations from a few percent to tens of percent for the models
considered. If a modification of gravity is correct, these three signals are
correlated, which gives further leverage in detecting small departures from
general relativity.

Galaxy clusters are particularly interesting candidates for probing subtleties
of structure growth, for two reasons: current and upcoming experiments compile
large cluster catalogs with well-understood completeness, and cluster catalogs
compiled through their thermal Sunyaev-Zeldovich distortions of the cosmic
microwave background detect all clusters along the line of sight regardless of
redshift. Arcminute-scale temperature decrements at microwave frequencies
below the 220 GHz null are an unambiguous cluster signature. The latest
arcminute-resolution microwave maps over large sky areas from the Atacama
Cosmology Telescope yield a catalog of over 4000 galaxy clusters with masses
above $5\times 10^{14}M_{\odot}$ Hilton et al. (2021), and the next data
release will likely double this number to give the largest cluster catalog to
date. Cluster SZ signals detected at greater than $5\sigma$ statistical
significance have been demonstrated to give essentially a complete mass-
limited cluster catalog over a given sky region Hasselfield et al. (2013) for
clusters at $z>0.2$; the cluster mass limit therefore decreases as map
sensitivity improves. Our cluster sample is anticipated to increase to at
least 16,000 clusters with masses down to $2\times 10^{14}M_{\odot}$ detected
at the same significance with the upcoming Simons Observatory Ade et al.
(2019), and likely many more down to lower cluster mass limits with CMB-S4.

Masses for clusters can be obtained through the correlation of cluster mass
with the thermal SZ signal, although galactic emission from within clusters
complicates this route Dicker et al. (2021). Gravitational lensing masses will
ultimately be available, either through weak lensing in large optical surveys
(see, e.g., Murata et al. (2019)) or lensing of the microwave background
Madhavacheril et al. (2015, 2020). Either of these techniques can be used to
compile the cluster mass function. The cluster power spectrum and the mean
pairwise velocity as a function of separation require positions in three-
dimensional space, which can be obtained from cluster redshifts, either
photometric or spectroscopic. The Dark Energy Survey and the upcoming Rubin
Telescope will provide photometric redshifts for sky regions largely
overlapping with ACT, Simons Observatory, and CMB-S4, while current large
spectroscopic survey instruments have demonstrated the capability to obtain
great numbers of spectroscopic redshifts.

Cluster velocities along the line of sight induce the kinematic Sunyaev-
Zeldovich temperature shift in the microwave background; this signal has been
detected directly via the mean pairwise momentum statistic Calafut et al.
(2021); Bernardis et al. (2017); Hand et al. (2012); and P. A. R. Ade et al.
(2016); Soergel et al. (2016) and by cross-correlating with a cosmic velocity
field inferred from the large-scale galaxy distribution Schaan et al. (2016).
Extracting cluster velocities from the kSZ signal is complicated by the signal
being proportional to the bulk gas momentum rather than velocity; any claimed
signature of modified gravity might also arise as a systematic error in
cluster gas modeling.

A cleaner alternative is the transverse velocity lensing signal Birkinshaw and
Gull (1983); Gurvits and Mitrofanov (1986), a dipole CMB temperature
distortion which is proportional to the cluster transverse velocity and the
total cluster mass. Combined with cluster mass estimates from lensing, this
signal holds the promise of a clean velocity probe. Its amplitude is small,
and so far undetected, although a first detection is likely with coming
experiments Hotinli et al. (2019); Yasini et al. (2019).

Given these basic observables, estimators for the galaxy cluster statistics
considered in this paper can be formulated, accounting for various noise
sources which must be mitigated. SZ-selected galaxy clusters avoid some of the
more difficult selection effects associated with many astronomical catalogs.
While precise masses of individual clusters from microwave background lensing
are far off, statistical mass determinations will rapidly improve, and optical
weak lensing masses from LSST (whose sky coverage is an excellent match for
ACT, Simons Observatory, and CMB-S4) will leverage significantly greater depth
and image quality compared to current optical surveys. Transverse lensing of
the microwave background to determine cluster velocities likely has
significantly lower astrophysical systematics compared to other velocity
probes such as redshift-space distortions or kinematic SZ. The ultimate
limiting systematic effects in the galaxy cluster statistics considered here
are not well understood; detailed studies will be required to determine the
eventual ability of galaxy cluster statistics to constrain modifications of
gravity.

The analysis presented here demonstrates that if one galaxy cluster statistic
departs from the predictions of standard cosmological models due to a
modification of gravity, other statistics will also display potentially
measurable differences. This provides a way to distinguish subtle differences
in structure growth due to modified gravity from possible systematic errors in
these challenging measurements. Correlations between the statistics computed
here in generic modifications of gravity will be considered elsewhere. Other
constraints on modified gravity from galaxy clustering and cluster abundances
can be found in Ref. Liu et al. (2021).

An alternate approach to probing structure formation is to construct
estimators of the total baryon velocity or density by cross-correlating
Fourier modes of the the entire distributions of galaxies and microwave
background signals (for example, Ref. Hotinli et al. (2021)) rather than using
only the galaxy clusters as tracers of the mass distribution. These types of
estimators should be equivalent to using all of the halos down to some mass
limit below that of galaxy clusters. This technique has the advantage of
simplicity and of using all available data; but it also picks out spurious
signals from other unrelated signals in the data which must be subtracted
through modeling. It remains to be seen whether the increase in data is an
advantage over using only galaxy clusters given the increased difficulty in
accounting for systematic errors. Ultimately, estimating structure growth
using multiple different techniques with different systematics is the best
route to robust conclusions about subtle signals such as those from modified
gravity.

The results of our halo-matching exercise shows that a peak-patch description
of halos in modified gravity should be equally good as for general relativity.
We have shown that large galaxy clusters are in one-to-one correspondence in
the modified gravity and general relativity simulations with identical initial
conditions. Individual cluster masses differ due to systematic differences in
linear growth and random differences in merger histories. This suggests the
possibility of constructing rapid approximate calculations of galaxy cluster
locations, masses, and velocities as has been done for standard general
relativity with the Websky code Stein et al. (2020). This peak-patch code
needs to be modified for a different linear growth function, and masses for
the resulting cluster-sized halos determined by abundance matching with the
kind of modified gravity simulations used in this paper instead of
conventional gravity simulations. The abundance matching procedure naturally
incorporates astrophysical differences in clusters between general relativity
and modified gravity such as tidal effects or different frequency of minor
mergers. The linear growth factor can generally be determined (analytically or
numerically) for a given gravity theory. So in principle it is straightforward
to create Websky-style sky simulations for other gravity theories which are
close to general relativity. We will address this possibility in more detail
elsewhere.

One other interesting extension of this work would be to use a general
parameterization of modified gravity, such as the Parameterized Post-Friedman
formulation of modified gravity Baker et al. (2013). Combining a general
framework such as this, plus Websky-style simulations of large-scale structure
and corresponding microwave background signatures, will provide an efficient
method of constructing useful mocks for confronting upcoming cosmological data
with modified gravity models. Given the pressing nature of the fundamental
physics questions surrounding the accelerating expansion of the universe, we
look forward to the near future when new probes of gravity on cosmological
scales become available.

## Acknowledgements

We thank Baojiu Li and Wojciech Hellwing on behalf of Cautun et al. (2018) and
Hellwing et al. (2017) for making available the _ELEPHANT_ and nDGP
simulations, respectively.

## References

  * Clifton et al. (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physics Reports 513, 1 (2012), URL https://doi.org/10.1016%2Fj.physrep.2012.01.001.
  * Will (2014) C. M. Will, Living Reviews in Relativity 17 (2014), URL https://doi.org/10.12942%2Flrr-2014-4.
  * Abbott et al. (2016) B. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al., Physical Review Letters 116 (2016), URL https://doi.org/10.1103%2Fphysrevlett.116.061102.
  * Abbott et al. (2017) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, et al., The Astrophysical Journal 848, L13 (2017), URL https://doi.org/10.3847/2041-8213/aa920c.
  * Langlois et al. (2018) D. Langlois, R. Saito, D. Yamauchi, and K. Noui, Physical Review D 97 (2018), URL https://doi.org/10.1103%2Fphysrevd.97.061501.
  * Bocquet et al. (2015) S. Bocquet, A. Saro, K. Dolag, and J. J. Mohr, Monthly Notices of the Royal Astronomical Society 456, 2361 (2015), URL https://doi.org/10.1093%2Fmnras%2Fstv2657.
  * Ade et al. (2016) P. A. R. Ade, N. Aghanim, F. Argüeso, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., Astronomy & Astrophysics 594, A26 (2016), URL https://doi.org/10.1051%2F0004-6361%2F201526914.
  * Ade et al. (2019) P. Ade, J. Aguirre, Z. Ahmed, S. Aiola, A. Ali, D. Alonso, M. A. Alvarez, K. Arnold, P. Ashton, J. Austermann, et al., Journal of Cosmology and Astroparticle Physics 2019, 056 (2019), URL https://doi.org/10.1088%2F1475-7516%2F2019%2F02%2F056.
  * Bleem et al. (2020) L. E. Bleem, S. Bocquet, B. Stalder, M. D. Gladders, P. A. R. Ade, S. W. Allen, A. J. Anderson, J. Annis, M. L. N. Ashby, J. E. Austermann, et al., The Astrophysical Journal Supplement Series 247, 25 (2020), URL https://doi.org/10.3847%2F1538-4365%2Fab6993.
  * Hilton et al. (2021) M. Hilton, C. Sifón, S. Naess, M. Madhavacheril, M. Oguri, E. Rozo, E. Rykoff, T. M. C. Abbott, S. Adhikari, M. Aguena, et al., The Astrophysical Journal Supplement Series 253, 3 (2021), URL https://doi.org/10.3847%2F1538-4365%2Fabd023.
  * Murata et al. (2019) R. Murata, M. Oguri, T. Nishimichi, M. Takada, R. Mandelbaum, S. More, M. Shirasaki, A. J. Nishizawa, and K. Osato, Publications of the Astronomical Society of Japan 71 (2019), URL https://doi.org/10.1093%2Fpasj%2Fpsz092.
  * Madhavacheril et al. (2015) M. Madhavacheril, N. Sehgal, R. Allison, N. Battaglia, J. R. Bond, E. Calabrese, J. Caligiuri, K. Coughlin, D. Crichton, R. Datta, et al., Physical Review Letters 114 (2015), URL https://doi.org/10.1103%2Fphysrevlett.114.151302.
  * Calafut et al. (2021) V. Calafut, P. Gallardo, E. Vavagiakis, S. Amodeo, S. Aiola, J. Austermann, N. Battaglia, E. Battistelli, J. Beall, R. Bean, et al., Physical Review D 104 (2021), URL https://doi.org/10.1103%2Fphysrevd.104.043502.
  * Schaan et al. (2021) E. Schaan, S. Ferraro, S. Amodeo, N. Battaglia, S. Aiola, J. E. Austermann, J. A. Beall, R. Bean, D. T. Becker, R. J. Bond, et al., Physical Review D 103 (2021), URL https://doi.org/10.1103%2Fphysrevd.103.063513.
  * Hotinli et al. (2019) S. C. Hotinli, J. Meyers, N. Dalal, A. H. Jaffe, M. C. Johnson, J. B. Mertens, M. Münchmeyer, K. M. Smith, and A. van Engelen, Physical Review Letters 123 (2019), URL https://doi.org/10.1103%2Fphysrevlett.123.061301.
  * Stein et al. (2018) G. Stein, M. A. Alvarez, and J. R. Bond, Monthly Notices of the Royal Astronomical Society 483, 2236 (2018), URL https://doi.org/10.1093%2Fmnras%2Fsty3226.
  * Stein et al. (2020) G. Stein, M. A. Alvarez, J. R. Bond, A. van Engelen, and N. Battaglia, Journal of Cosmology and Astroparticle Physics 2020, 012 (2020), URL https://doi.org/10.1088%2F1475-7516%2F2020%2F10%2F012.
  * Hu and Sawicki (2007) W. Hu and I. Sawicki, Physical Review D 76 (2007), URL https://doi.org/10.1103%2Fphysrevd.76.064004.
  * Dvali et al. (2000) G. Dvali, G. Gabadadze, and M. Porrati, Physics Letters B 485, 208 (2000), URL https://doi.org/10.1016%2Fs0370-2693%2800%2900669-9.
  * Valogiannis and Bean (2019) G. Valogiannis and R. Bean, Physical Review D 99 (2019), URL https://doi.org/10.1103%2Fphysrevd.99.063526.
  * Alam et al. (2021) S. Alam, C. Arnold, A. Aviles, R. Bean, Y.-C. Cai, M. Cautun, J. L. Cervantes-Cota, C. Cuesta-Lazaro, N. C. Devi, A. Eggemeier, et al., Journal of Cosmology and Astroparticle Physics 2021, 050 (2021), URL https://doi.org/10.1088%2F1475-7516%2F2021%2F11%2F050.
  * Lombriser et al. (2009) L. Lombriser, W. Hu, W. Fang, and U. Seljak, Physical Review D 80 (2009), URL https://doi.org/10.1103%2Fphysrevd.80.063536.
  * Cautun et al. (2018) M. Cautun, E. Paillas, Y.-C. Cai, S. Bose, J. Armijo, B. Li, and N. Padilla, Monthly Notices of the Royal Astronomical Society 476, 3195 (2018), URL https://doi.org/10.1093%2Fmnras%2Fsty463.
  * Li et al. (2012) B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama, Journal of Cosmology and Astroparticle Physics 2012, 051 (2012), URL https://doi.org/10.1088%2F1475-7516%2F2012%2F01%2F051.
  * Barreira et al. (2015) A. Barreira, S. Bose, and B. Li, Journal of Cosmology and Astroparticle Physics 2015, 059 (2015), URL https://doi.org/10.1088%2F1475-7516%2F2015%2F12%2F059.
  * Teyssier (2002) R. Teyssier, Astronomy & Astrophysics 385, 337 (2002), URL https://doi.org/10.1051%2F0004-6361%3A20011817.
  * Blas et al. (2011) D. Blas, J. Lesgourgues, and T. Tram, Journal of Cosmology and Astroparticle Physics 2011, 034 (2011), URL https://doi.org/10.1088%2F1475-7516%2F2011%2F07%2F034.
  * Behroozi et al. (2012) P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, The Astrophysical Journal 762, 109 (2012), URL https://doi.org/10.1088%2F0004-637x%2F762%2F2%2F109.
  * Hinshaw et al. (2013) G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, J. Dunkley, M. R. Nolta, M. Halpern, R. S. Hill, N. Odegard, et al., The Astrophysical Journal Supplement Series 208, 19 (2013), URL https://doi.org/10.1088%2F0067-0049%2F208%2F2%2F19.
  * Hand et al. (2018) N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak, and Z. Slepian, The Astronomical Journal 156, 160 (2018), URL https://doi.org/10.3847%2F1538-3881%2Faadae0.
  * Hernández-Aguayo et al. (2018) C. Hernández-Aguayo, C. M. Baugh, and B. Li, Monthly Notices of the Royal Astronomical Society 479, 4824 (2018), URL https://doi.org/10.1093%2Fmnras%2Fsty1822.
  * Davis and Peebles (1977) M. Davis and P. J. E. Peebles, The Astrophysical Journal Supplement Series 34, 425 (1977).
  * Vikram et al. (2013) V. Vikram, A. Cabré, B. Jain, and J. VanderPlas, Journal of Cosmology and Astroparticle Physics 2013, 020 (2013), URL https://doi.org/10.1088%2F1475-7516%2F2013%2F08%2F020.
  * Hasselfield et al. (2013) M. Hasselfield, M. Hilton, T. A. Marriage, G. E. Addison, L. F. Barrientos, N. Battaglia, E. S. Battistelli, J. R. Bond, D. Crichton, S. Das, et al., Journal of Cosmology and Astroparticle Physics 2013, 008 (2013), URL https://doi.org/10.1088%2F1475-7516%2F2013%2F07%2F008.
  * Dicker et al. (2021) S. R. Dicker, E. S. Battistelli, T. Bhandarkar, M. J. Devlin, S. M. Duff, G. Hilton, M. Hilton, A. D. Hincks, J. Hubmayr, K. Huffenberger, et al., Monthly Notices of the Royal Astronomical Society 508, 2600 (2021), URL https://doi.org/10.1093%2Fmnras%2Fstab2679.
  * Madhavacheril et al. (2020) M. S. Madhavacheril, C. Sifón, N. Battaglia, S. Aiola, S. Amodeo, J. E. Austermann, J. A. Beall, D. T. Becker, J. R. Bond, E. Calabrese, et al., The Astrophysical Journal 903, L13 (2020), URL https://doi.org/10.3847%2F2041-8213%2Fabbccb.
  * Bernardis et al. (2017) F. D. Bernardis, S. Aiola, E. Vavagiakis, N. Battaglia, M. Niemack, J. Beall, D. Becker, J. Bond, E. Calabrese, H. Cho, et al., Journal of Cosmology and Astroparticle Physics 2017, 008 (2017), URL https://doi.org/10.1088%2F1475-7516%2F2017%2F03%2F008.
  * Hand et al. (2012) N. Hand, G. E. Addison, E. Aubourg, N. Battaglia, E. S. Battistelli, D. Bizyaev, J. R. Bond, H. Brewington, J. Brinkmann, B. R. Brown, et al., Physical Review Letters 109 (2012), URL https://doi.org/10.1103%2Fphysrevlett.109.041101.
  * and P. A. R. Ade et al. (2016) and P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, E. Aubourg, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., Astronomy & Astrophysics 586, A140 (2016), URL https://doi.org/10.1051%2F0004-6361%2F201526328.
  * Soergel et al. (2016) B. Soergel, S. Flender, K. T. Story, L. Bleem, T. Giannantonio, G. Efstathiou, E. Rykoff, B. A. Benson, T. Crawford, S. Dodelson, et al., Monthly Notices of the Royal Astronomical Society 461, 3172 (2016), URL https://doi.org/10.1093%2Fmnras%2Fstw1455.
  * Schaan et al. (2016) E. Schaan, S. Ferraro, M. Vargas-Magañ a, K. M. Smith, S. Ho, S. Aiola, N. Battaglia, J. R. Bond, F. D. Bernardis, E. Calabrese, et al., Physical Review D 93 (2016), URL https://doi.org/10.1103%2Fphysrevd.93.082002.
  * Birkinshaw and Gull (1983) M. Birkinshaw and S. F. Gull, Nature 302, 315 (1983).
  * Gurvits and Mitrofanov (1986) L. I. Gurvits and I. G. Mitrofanov, Nature 324, 349 (1986).
  * Yasini et al. (2019) S. Yasini, N. Mirzatuny, and E. Pierpaoli, The Astrophysical Journal 873, L23 (2019), URL https://doi.org/10.3847%2F2041-8213%2Fab0bfe.
  * Liu et al. (2021) R. Liu, G. Valogiannis, N. Battaglia, and R. Bean, Physical Review D 104 (2021), URL https://doi.org/10.1103%2Fphysrevd.104.103519.
  * Hotinli et al. (2021) S. C. Hotinli, M. C. Johnson, and J. Meyers, Physical Review D 103 (2021), URL https://doi.org/10.1103%2Fphysrevd.103.043536.
  * Baker et al. (2013) T. Baker, P. G. Ferreira, and C. Skordis, Physical Review D 87 (2013), URL https://doi.org/10.1103%2Fphysrevd.87.024015.
  * Hellwing et al. (2017) W. A. Hellwing, K. Koyama, B. Bose, and G.-B. Zhao, Physical Review D 96 (2017), URL https://doi.org/10.1103%2Fphysrevd.96.023515.


11institutetext: Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, France
11email<EMAIL_ADDRESS>22institutetext: Osservatorio Astrofisico di
Torino, Istituto Nazionale di Astrofisica (INAF), I-10025 Pino Torinese, Italy

# The chemical signature of the Galactic spiral arms
revealed by Gaia DR3

E. Poggio 1122 A. Recio-Blanco 11 P. A. Palicio 11 P. Re Fiorentin 22 P. de
Laverny 11 R. Drimmel 22
G. Kordopatis 11 M. G. Lattanzi 22 M. Schultheis 11 A. Spagna 22 E. Spitoni 11

(Received XXXX; accepted YYY)

Taking advantage of the recent _Gaia_ Data Release 3 (DR3), we map chemical
inhomogeneities in the Milky Way’s disc out to a distance of $\sim$ 4 kpc from
the Sun, using different samples of bright giant stars. The samples are
selected using effective temperatures and surface gravities from the GSP-Spec
module, and are expected to trace stellar populations of different typical
age. The cool (old) giants exhibit a relatively smooth radial metallicity
gradient with an azimuthal dependence. Binning in Galactic azimuth $\phi$, the
slope gradually varies from $d$[M/H]$/dR\sim-0.054$ dex kpc-1 at
$\phi\sim-20^{\circ}$ to $\sim-0.035$ dex kpc-1 at $\phi\sim 20^{\circ}$. On
the other hand, the relatively hotter (and younger) stars present remarkable
inhomogeneities, apparent as three (possibly four) metal-rich elongated
features in correspondence of the spiral arms’ locations in the Galactic disc.
When projected onto Galactic radius, those features manifest themselves as
statistically significant bumps on top of the observed radial metallicity
gradients with amplitudes up to $\sim 0.05-0.1$ dex, making the assumption of
a linear radial decrease not applicable to this sample. The strong correlation
between the spiral structure of the Galaxy and the observed chemical pattern
in the young sample indicates that the spiral arms might be at the origin for
the detected chemical inhomogeneities. In this scenario, the spiral arms would
leave in the younger stars a strong signature, which progressively disappears
when cooler (and older) giants are considered.

##  1 Introduction

Figure 1: _Chemical inhomogenities in the Galactic disc for three different
samples._ _Left panel_ : Selection of Sample A, B and C in the Kiel diagram.
_Right panels_ : Mean metallicity in heliocentric coordinates for Sample A, B,
and C (from left to right, respectively) in the Galactic plane. The position
of the Sun is shown by the white cross in (X,Y)=(0 kpc, 0 kpc); the Galactic
center is to the right; Galactic rotation is going clockwise.

It has been known for several decades that the disc of the Milky Way contains
large-scale non-axisymmetric features, including the spiral arms and a central
bar (see e.g. Georgelin & Georgelin, 1976; Okuda et al., 1977; Shen & Zheng,
2020; Gaia Collaboration, Drimmel et al., 2022, and references threin), which
can play an important role in the disc’s evolution, and are expected to leave
their signature in the observed properties of the stars.

Several models have tried to explain the behaviour of spiral arms in disc
galaxies, although their dynamical nature and physical origin still remain
unknown. Spiral arms have been proposed to induce large-scale shock waves that
trigger the gravitational collapse of gas clouds, enhancing star formation
(Roberts, 1969). They can also induce radial migration of stars (Lynden-Bell &
Kalnajs, 1972; Sellwood & Binney, 2002; Schönrich & Binney, 2009a, b; Roškar
et al., 2012), and trap or scatter stars close to orbital resonances
(Contopoulos & Grosbol, 1986). As a consequence, they can generate variations
in the mean chemical composition of stars in the arm and inter-arm regions
(Minchev et al., 2012; Kordopatis et al., 2015; Grand et al., 2016; Khoperskov
et al., 2018; Spitoni et al., 2019; Khoperskov & Gerhard, 2021).

Mapping the spatial distribution of metals within galaxies therefore
represents a key aspect to study the processes of chemical enrichment and
mixing in the interstellar medium. An extensive review of the chemical
enrichment of galaxies across the cosmic epochs can be found in Maiolino &
Mannucci (2019) (see also references therein). In external galaxies, several
works have found evidence of azimuthal variations in the mean chemical
composition of stars (e.g. Sánchez-Menguiano et al., 2016; Ho et al., 2017;
Vogt et al., 2017). Recently, using data from the PHANGS–MUSE survey (Emsellem
et al., 2022), Williams et al. (2022) found that most of the star-forming
galaxies considered in their sample showed significant 2D metallicity
variations.

In the Milky Way, chemical inhomogenities were explored by several studies.
Balser et al. (2011) analyzed the oxygen distribution of HII regions in the
Galactic disc, divided their sample into three Galactic azimuth bins, and
found that the corresponding radial gradients were significantly different.
Chemical inhomogeneities in iron and oxygen were also found in Cepheids
(Pedicelli et al., 2009; Genovali et al., 2014; Kovtyukh et al., 2022). Large
metallicity variations were also found in the interstellar medium (De Cia et
al., 2021). By combining data from HII regions, Cepheids, B stars and red
supergiants, Davies et al. (2009) presented evidence for large azimuthal
variations ($\sim$ 0.4 dex) in oxygen, magnesium and silicon on a scale of few
kpc in the inner Galactic disc. On the other hand, using APOGEE red-clump
stars, Bovy et al. (2014) constrained azimuthal variations in the median
metallicity to be $\leq$ 0.02 dex in the region covered by their sample. Using
the large spectroscopic survey RAVE and the Geneva Copenhagen Survey, Antoja
et al. (2017) analyzed stars in a cylinder of 0.5 kpc radius centered on the
Sun, and found asymmetric metallicity patterns in velocity space.

Recently, _Gaia_ Data Release 3 (hereafter DR3, Gaia Collaboration, Vallenari
et al., 2022) published the largest stellar catalog with chemical abundances,
atmospheric parameters and radial velocities ever created. Radial Velocity
Spectrometer (RVS) data were parameterised by the General Stellar
Parameteriser - spectroscopy (_GSP-Spec_) module (Recio-Blanco et al., 2022),
delivering chemo-physical parameters for 5.6 million stars over the entire
sky. Gaia Collaboration, Recio-Blanco et al. (2022) observed azimuthal
variations in red giant branch stars, but didn’t present quantitative results
or in-depth discussions on their origin.

Thanks to the large number of chemical and astrometric measurements, as well
as _Gaia_ ’s all-sky sampling, it is now possible to map chemical
inhomogeneities in the Milky Way’s disc as never before, with the goal of
detecting the possible signature left by the Galactic spiral arms in the
stellar metallicity.

The Letter is structured as follows: in Section 2 we describe the datasets; in
Section 3 we present our results; in Section 4 we discuss our findings and
future perspectives.

##  2 Data selection

In this study, we use the main stellar atmospheric parameters (effective
temperature $T_{\rm eff}$, surface gravity log($g$), global metallicity [M/H])
derived from _Gaia_ RVS spectra by the _GSP-Spec_ module, to select tracers of
the disc population and map large-scale chemical inhomogeneities in our
Galaxy. The quality constraints on spectroscopic and astrometric measurements
(as well as a Toomre kinematical selection) are detailed in Appendix A.

To serve our purpose, we focus on giant stars, as they are intrinsically
bright, and therefore allow us to sample a relatively large volume of the
Galactic disc. Also, taking advantage of the large number of high-quality
chemical measurements for these stars111As an indication, for stars with
log($g$)$<$ 2.5 and $6000<$ $T_{\rm eff}$$<3500$ K, the metallicity
uncertainty is $<$ 0.1 dex for 1.4 million stars, $<$ 0.05 dex for $\sim$680
000 stars, and $<$ 0.01 dex for $\sim$43 000 stars., a robust statistical
analysis can be performed, thanks to high signal-to-noise spectra.

Figure 1 (left panel) shows the portion of the Kiel diagrams populated by the
giant stars in our selected sample. Here we select stars brighter than the red
clump (which is located at approximately log($g$)$\simeq$ 2.3 and $T_{\rm
eff}$$\simeq$ 4750 K), to avoid selection function effects due to the
superposition of different stellar populations (see Section 4 of Gaia
Collaboration, Recio-Blanco et al., 2022). Indeed, such effects can cause
artefacts when investigating radial gradients, which might be erroneously
interpreted as signatures of the spiral arms. Selected areas in Fig. 1 (left
panel) define three different subsamples (details in Appendix A), which cover
different portions of the red giant branch. Additionally, we apply a vertical
cut $|Z|=|d\sin{b}|<0.75$ kpc, where $b$ is the galactic latitude and $d$ is
the heliocentric distance according to the geometric model by Bailer-Jones et
al. (2021). The adopted samples are here labelled as A, B, and C, and contain
respectively 19 340, 151 139, and 346 570 stars.

The portion of the Kiel diagram covered by our three samples is expected to be
populated by stars of different ages. Both theoretical considerations (based
on isochrones) and empirical evidence (based on spatial distribution and
stellar kinematics) indicate that Sample A is expected to be typically younger
than the other two samples (more details in Appendix B). For instance, the
comparison between the selected areas on the Kiel diagram and the PARSEC
isochrones (Bressan et al., 2012; Chen et al., 2014, 2015; Tang et al., 2014;
Pastorelli et al., 2019) for metallicities between $[M/H]=$ 0 dex and
$[M/H]=$-0.5 dex suggests that Sample A should mostly contain stars younger
than 100-300 Myr (taking into account typical uncertainties on $T_{\rm eff}$
and log($g$)), while Sample C should be dominated by much older stars, i.e.
$\gtrsim$ 1 Gyr (see Fig. 4). On the other hand, Sample B is expected to be
intermediate between the two, containing a mixture of young and old stellar
populations. As will be shown in the following, the different content of the
three samples will determine remarkable differences in their observed
metallicity maps.

##  3 Results

Figure 1 (right panels) shows the maps of the mean metallicity
$\langle\,$[M/H]$\,\rangle$ in the Galactic plane for Samples A, B and C,
using a smoothing Gaussian kernel with a bandwidth of 175 pc (see Appendix C).
As discussed in the following, the three maps exhibit both similarities and
differences. On a large scale, stars are typically more metal-rich toward the
inner parts of the Galaxy for all three samples. This is in agreement with
previous observations of radial metallicity gradients in the Galactic disc
(e.g. Genovali et al., 2014; Gaia Collaboration, Recio-Blanco et al., 2022),
and also expected from the inside-out formation scenario (Bird et al., 2013;
Matteucci, 2021). However, the maps also present smaller-scale metallicity
inhomogenehities. The relative importance, extent, and shape of the observed
inhomogeneities vary for the three different samples. Sample A exhibits three,
possibly four, metal-rich elongated features, that diagonally cross the
portion of the XY-plane covered by our dataset. One feature extends
approximately from (X,Y)=($-$1.5 kpc, 3 kpc) to ($-$3 kpc, $-$1 kpc). A second
feature stretches from (X,Y)=(1 kpc, 3.5 kpc) to (X,Y)=($-$1.5 kpc, $-$1 kpc),
and then continues almost vertically out to (X,Y)=($-$1.5 kpc, $-$3 kpc). A
third feature (which possibly contains multiple segments) extends from
(X,Y)=(2.5 kpc, 1 kpc) to (X,Y)=(1 kpc, $-$4 kpc). This pattern becomes
progressively less evident as we move from Sample A to Sample B and C. In
Sample C, the elongated features disappear almost completely, and simply
result in an observed asymmetry about the Sun of stars at $X\simeq$ 0 kpc,
typically more metal rich at $Y>0$ kpc than at $Y<0$ kpc. While such asymmetry
was already noted by Gaia Collaboration, Recio-Blanco et al. (2022) in their
RGB sample, the scenario presented in this contribution indicates that it
might represent a remnant of the spiral arms chemical signature.

Figure 2:  _Correlation between the observed metal-rich features in Sample A
and the segments of the nearest spiral arms in the Galaxy._ _Upper left_
panel: Map of the mean metallicity in the Galactic plane for Sample A shown in
Fig. 1, but now compared to the position of the spiral arms from Upper Main
Sequence stars (UMS) in Poggio et al. (2021), shown as grey shaded areas. The
labels show the position of the Perseus arm, the Local (Orion) arm and the
Sagittarius-Carina arm based on the UMS map. _Upper right panel_ : bivariate
kernel density estimation for Sample A using an Epanechnikov kernel with a
smoothing length of 250 pc. The contours are the same as the upper left panel.
_Lower left panel_ : Same as Upper panels, but now showing the metallicity
excess. Regions with positive (negative) metallicity excess are more metal-
rich (-poor) than the average (see text for a detailed explanation). _Lower
right panel_ : Using the map shown in the lower left panel, the overdensity of
the spiral arms vs. the metallicity excess is shown for each pixel (black
dots). When the spiral arm overdensity is $>0$, the region can be considered
belonging to a spiral arm, whereas regions with $<0$ can be considered
interarm regions. Black contours show the distribution of the pixels in the
spiral arms overdensity vs. metallicity excess space. The red line shows a
linear fit to the points.

Figure 2 (upper left panel) shows again the mean metallicity for Sample A, but
now compared to the position of the spiral arms in the Galaxy as mapped by
Upper Main Sequence UMS stars in Poggio et al. (2021), shown as grey shaded
areas and black contours, indicating (from left to right) the segments of: the
Perseus arm, the Local (Orion) arm, and the Sagittarius-Carina (Sag-Car) arm
(the latter contour possibly containing the Scutum arm).

In the upper right panel of Fig. 2, we can see the spatial distribution of
sources in Sample A, obtained with an Epanechnikov kernel density estimator
with a bandwidth of 250 pc. White contours show again the position of the
nearest spiral arms from UMS stars. As we can see, Sample A traces well the
position of the spiral arms, confirming that it contains typically young
stars.

To better analyse and map chemical inhomogeneities in the Galactic disc, we
define a new variable, called Metallicity Excess, defined as
$\langle\,\rm{[M/H]}\,\rangle_{\rm{loc}}-\langle\,\rm{[M/H]}\,\rangle_{\rm{large}}$,
where $\langle\,[\rm{M/H}]\,\rangle_{\rm{loc}}$ and
$\langle\,[\rm{M/H}]\,\rangle_{\rm{large}}$ represent the mean metallicity
smoothed on a local or large scale, respectively. The smoothing is performed
as described in Appendix C. The lower left panel of Fig. 2 shows the obtained
Metallicity Excess for Sample A if we adopt for the local smoothing the same
value chosen for the mean metallicity in the upper left panel of Fig. 2 (i.e.
$\langle\,[\rm{M/H}]\,\rangle_{\rm{loc}}=\langle\,[\rm{M/H}]\,\rangle$), and
for the large-scale metallicity a scale-length of 5 times the value adopted
for the local smoothing. The red (blue) regions in the Metallicity Excess map
should be interpreted as places where the stars are more metal-rich (-poor)
than the average.

We therefore study the correlation between the observed metal-rich regions and
the geometry of the spiral arms by dissecting the XY map (e.g. lower left
panel of Fig. 2) in pixels of 0.3 kpc width. For each pixel, we consider the
spiral arms overdensity (calculated as described in Poggio et al., 2021) and
the Metallicity Excess. The bottom right panel of Fig. 2 shows the correlation
between these two quantities (every black dot corresponds to a pixel in the XY
map). The red line shows a linear fit to the black dots. As we can see,
regions located in correspondence of the spiral arms (i.e. with overdensity
larger than 0) tend to be more metal-rich than the average (i.e. with a
Metallicity Excess larger than 0). To quantify the correlation between these
two variables, we calculate Kendall’s correlation coefficient $\tau_{B}=0.38$,
which can be considered as an indication of strong positive correlation222A
correlation is usually considered very weak for values less than 0.10, weak
between 0.10 and 0.19, moderate between 0.20 and 0.29, strong for 0.30 or
above (Botsch, 2011).. Therefore, the elongated metal-rich features observed
in Sample A appear to be statistically correlated with the position of the
spiral arms in the Galaxy.

It should be noted, however, that there are some exceptions. For instance, the
enhancement of metallicity at approximately (X,Y)=( 0 kpc, $-$1 kpc) does not
seem to correspond to an overdensity in the spatial distribution. The
comparison between our map and the data from Hottier et al. (2021) reveals
that it might correspond to the Vela Molecular Ridge. This region can be seen
as a dust enhancement in the maps of Lallement et al. (2019); Vergely et al.
(2022), and possibly prevents us from seeing an overdensity in the spatial
distribution of stars (for both Sample A and UMS stars, see the upper right
panel in Fig. 2).

Furthermore, it is interesting to note that not all the regions in the spiral
arms exhibit the same chemical behaviour: the metallicity excess seems to be
typically high in the Sag-Car arm and in the upper part (Y¿0) of the Local
arm, but seems to be milder in the lower part of the Local arm (Y¡0) and the
Perseus arm (Y $\lesssim$ 1 kpc).

Figure 3 shows the impact of the above-described metal-rich features on the
observed metallicity gradient in the Galactic disc. To explore the metallicity
variations in the disc, we dissect the Galactic plane into slices of 10∘ width
in Galactic azimuth $\phi$. For each slice, we then show the median
metallicity for Sample A and Sample C. Sample B presents an intermediate
behaviour between Sample A and C (here not shown to preserve the clarity of
the Figure). As we can see, Sample A presents some peaks of higher
metallicity, superimposed on a global decrease as a function of $R$, here
calculated assuming a distance to the Galactic center $R_{\odot}=8.249$ kpc
(Gravity Collaboration et al., 2021). Based on bootstrap uncertainties, the
peaks deviate more than 3-sigma from a linear decrease as a function of
Galactic radius, indicating that the radial metallicity gradient is dominated
by chemical undulations, that reach maximum deviations of $\sim$ 0.05-0.1 dex.
In this context, the observed peaks can be interpreted as the projection on
the radial direction of the metal-rich features observed in the XY maps. For
Sample C, on the other hand, a relatively smooth decrease of the median
metallicity as a function of $R$ is apparent. The slopes of the radial
metallicity gradients for different azimuthal bins can be found in Table 1. As
we can see, the radial metallicity gradient gradually varies depending on
azimuth, and becomes gradually steeper for $\phi<0^{\circ}$.

Additional tests to verify the robustness of our results can be found in
Appendix E.

Figure 3:  _Impact of the observed chemical inhomogeneities on the Galactic gradients._ _Upper panel_ : Map of the mean metallicity in the Galactic plane for Sample A, with overlaid azimuthal slices of 5∘, from $\phi=-20^{\circ}$ to $\phi=20^{\circ}$. _Lower panel_ : median metallicity as a function of Galactocentric radius $R$ for different slices in Galactic azimut $\phi$ for Sample A (filled dots) and Sample C (open dots). Regions that might be associated with the Perseus (red, letters ‘a’–‘e’), Local (green, ‘f’–‘n’), and Sagittarius/Scutum (blue, ‘o’–‘x’) arms are shown as colored regions both in the upper and lower panels.  Table 1: Radial gradient for Sample C in different azimuthal bins, as shown in Fig. 3, assuming a linear model. Only stars with $R>6.5$ kpc are here considered, to avoid selection effects in the inner regions due to extinction, as well as $R<11.5$ kpc, to avoid low statistics regions. The corresponding table for Sample A and B is not presented here due to the spiral arms-related metallicity variations, which are invalidating the assumption of a linear model.  $\phi$ (deg) | Slope (dex kpc-1) | Intercept (dex)
---|---|---
$-$20 $\pm$ 5 |  $-$0.054 $\pm$ 0.001 | 0.15 $\pm$ 0.01
$-$15 $\pm$ 5 |  $-$0.052 $\pm$ 0.002 | 0.14 $\pm$ 0.01
$-$10 $\pm$ 5 |  $-$0.047 $\pm$ 0.002 | 0.09 $\pm$ 0.02
$-$5 $\pm$ 5 |  $-$0.046 $\pm$ 0.002 | 0.09 $\pm$ 0.01
0 $\pm$ 5 |  $-$0.047 $\pm$ 0.001 | 0.09 $\pm$ 0.01
5 $\pm$ 5 |  $-$0.046 $\pm$ 0.001 | 0.09 $\pm$ 0.01
10 $\pm$ 5 |  $-$0.043 $\pm$ 0.001 | 0.06 $\pm$ 0.01
15 $\pm$ 5 |  $-$0.039 $\pm$ 0.001 | 0.04 $\pm$ 0.01
20 $\pm$ 5 |  $-$0.035 $\pm$ 0.002 | 0.01 $\pm$ 0.01

##  4 Discussion and Conclusion

Taking advantage of the new _Gaia_ DR3 data, we have mapped the mean
metallicity in the Galactic disc out to 3-4 kpc of the Sun for three different
samples of bright giant stars. We found evidence of chemical inhomogeneities,
which manifest themselves as statistically significant undulations on top of
the expected radial metallicity gradients. The appearance and magnitude of the
detected inhomogeneities are different for the three considered samples, being
more structured and pronounced for the sample containing hotter (and a larger
relative fraction of young) stars. In this sample, three (possibly four)
extended metal rich features are detected, which appear to be statistically
correlated with the position of the spiral arms (independently derived from
previous works). Such connection points to the spiral arms as a plausible
origin for the observed metal-rich features.

It is interesting to compare our maps of the mean metallicity in the Milky
Way’s disc to the results based on external spiral galaxies. For the galaxy
HCG91c, Vogt et al. (2017) found that the chemical enrichment of the
interstellar medium proceeds preferentially along spiral structures, and less
efficiently across them, in agreement with our results. Williams et al. (2022)
considered a sample of 19 galaxies, for which detected significant chemical
variations in most cases, but found no clear signs of enrichment along the
spiral arms. On the other hand, Sánchez-Menguiano et al. (2020) detected the
presence of more metal-rich HII regions in the spiral arms with respect to the
corresponding interarm regions for a large subsample of galaxies, in agreement
with the maps presented here for the Milky Way.

While the observed inhomogeneities are more pronounced and structured for the
youngest stars, azimuthal variations in the old stellar populations also
deserve some discussion. The oldest stars considered here (Sample C) present a
metallicity gradient that changes with Galactic azimuth, with a slope that
varies from $\sim-$0.054 dex kpc-1 to $\sim-$0.035 dex kpc-1. To our
knowledge, this is the first time that such azimuthal dependence is shown for
old stellar populations, both in the Milky Way and in external galaxies. The
slopes obtained here for Sample C using different azimuthal bins are in good
agreement (considering the uncertainties) with those from the giant stars at
all azimuths in the Galactic plane in Gaia Collaboration, Recio-Blanco et al.
(2022) (see their Table 1), which have a slope of $-$0.055 $\pm$ 0.007 dex
kpc-1 above the plane, and $-$0.057 $\pm$ 0.007 dex kpc-1 below the plane.
Previous works reported slopes for thin disc (i.e. low $\alpha-$sequence)
field stars ranging from $\sim-$0.053 dex kpc-1 to $\sim-$0.068 dex kpc-1
(Bergemann et al., 2014; Recio-Blanco et al., 2014; Hayden et al., 2015;
Anders et al., 2017).

Recently, Hawkins (2022) mapped azimuthal variations in the Galactic disc,
using both _Gaia_ DR3 data and a sample OBAF-type stars from the LAMOST survey
published by Xiang et al. (2022). Using the LAMOST sample, the author found
azimuthal structure on top of the radial gradient, which doesn’t necessarily
follow the spiral arms. On the other hand, using _Gaia_ DR3, Hawkins (2022)
found good agreement with the results presented in this _Letter_ , finding
that giant stars near the spiral arms typically contain more metals than those
in the interarm regions. The author suggests that the discrepancy between the
two samples can be due to the fact that either LAMOST data is too local to see
the overall structure, or that azimuthal structure depends strongly on age,
selection function and spectral type of the tracer population.

Several works available in the literature can offer a theoretical framework
for the obtained observational results. Using high-resolution $N$ -body
simulations, Khoperskov et al. (2018) showed that kinematically hot and cold
stellar populations in the Galactic disc react in a different way to a spiral
arm perturbation, naturally leading to azimuthal variations in the mean
metallicity of stars in the simulated disc (based on the assumption that
younger stars typically tend to have higher metallicity and smaller random
motions than the older stellar components, see for example Holmberg et al.,
2007). Figure 4 in Khoperskov et al. (2018) shows the azimuthal metallicity
variations at different times of the simulation for three different initial
metallicity profiles. They found that the spiral arms locii typically tend to
be more metal-rich than the interarm-regions. Their predictions are in good
agreement with the observations presented in this work.

In the future, a detailed modelling will be crucial to fully explain the
observed chemical inhomogeneities in the Galactic disc. Chemo-dynamical
evolution models should be adopted to interpret the impact of the spiral arms
on chemical variations in the Galactic disc, and to understand how different
mechanisms can influence the observed present-day chemical abundance patterns
(e.g. Roškar et al., 2012; Grand et al., 2016; Minchev et al., 2018; Spitoni
et al., 2019, 2022a, 2022b; Carr et al., 2022). Future works comparing
observations from _Gaia_ DR3 and ground-based surveys to sophisticated
theoretical models will be able to show further insights on the evolutionary
processes that shaped the present-day appearance of the Milky Way’s disc.

###### Acknowledgements.

The authors thank the referee for comments and suggestions that improved the
overall quality of the paper. This work has made use of data from the European
Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed
by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.
esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by
national institutions, in particular the institutions participating in the
Gaia Multilateral Agreement. EP acknowledges support by the Centre National
d’études Spatiales (CNES). ARB and ES received funding from the European
Union’s Horizon 2020 research and innovation program under SPACE-H2020 grant
agreement number 101004214 (EXPLORE project).

## References

  * Anders et al. (2017) Anders, F., Chiappini, C., Minchev, I., et al. 2017, A&A, 600, A70
  * Antoja et al. (2017) Antoja, T., Kordopatis, G., Helmi, A., et al. 2017, A&A, 601, A59
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Balser et al. (2011) Balser, D. S., Rood, R. T., Bania, T. M., & Anderson, L. D. 2011, ApJ, 738, 27
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527
  * Bensby et al. (2014) Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71
  * Bergemann et al. (2014) Bergemann, M., Ruchti, G. R., Serenelli, A., et al. 2014, A&A, 565, A89
  * Bird et al. (2013) Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43
  * Botsch (2011) Botsch, R. 2011, Scopes and Methods of Political Science
  * Bovy et al. (2014) Bovy, J., Nidever, D. L., Rix, H.-W., et al. 2014, ApJ, 790, 127
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Carr et al. (2022) Carr, C., Johnston, K. V., Laporte, C. F. P., & Ness, M. K. 2022, arXiv e-prints, arXiv:2201.04133
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Chen et al. (2014) Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525
  * Contopoulos & Grosbol (1986) Contopoulos, G. & Grosbol, P. 1986, A&A, 155, 11
  * Davies et al. (2009) Davies, B., Origlia, L., Kudritzki, R.-P., et al. 2009, ApJ, 696, 2014
  * De Cia et al. (2021) De Cia, A., Jenkins, E. B., Fox, A. J., et al. 2021, Nature, 597, 206
  * Emsellem et al. (2022) Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191
  * Gaia Collaboration, Creevey et al. (2022) Gaia Collaboration, Creevey, O. L., Sarro, L. M., Lobel, A., et al. 2022, arXiv e-prints, arXiv:2206.05870
  * Gaia Collaboration, Drimmel et al. (2022) Gaia Collaboration, Drimmel, R., Romero-Gomez, M., Chemin, L., et al. 2022, arXiv e-prints, arXiv:2206.06207
  * Gaia Collaboration, Recio-Blanco et al. (2022) Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., de Laverny, P., et al. 2022, arXiv e-prints, arXiv:2206.05534
  * Gaia Collaboration, Vallenari et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., Prusti, T., & et al. 2022, A&A in press
  * Genovali et al. (2014) Genovali, K., Lemasle, B., Bono, G., et al. 2014, A&A, 566, A37
  * Georgelin & Georgelin (1976) Georgelin, Y. M. & Georgelin, Y. P. 1976, A&A, 49, 57
  * Grand et al. (2016) Grand, R. J. J., Springel, V., Kawata, D., et al. 2016, MNRAS, 460, L94
  * Gravity Collaboration et al. (2021) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2021, A&A, 654, A22
  * Hawkins (2022) Hawkins, K. 2022, arXiv e-prints, arXiv:2207.04542
  * Hayden et al. (2015) Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132
  * Ho et al. (2017) Ho, I. T., Seibert, M., Meidt, S. E., et al. 2017, ApJ, 846, 39
  * Holmberg et al. (2007) Holmberg, J., Nordström, B., & Andersen, J. 2007, A&A, 475, 519
  * Hottier et al. (2021) Hottier, C., Babusiaux, C., & Arenou, F. 2021, A&A, 655, A68
  * Khoperskov et al. (2018) Khoperskov, S., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 611, L2
  * Khoperskov & Gerhard (2021) Khoperskov, S. & Gerhard, O. 2021, arXiv e-prints, arXiv:2111.15211
  * Kordopatis et al. (2015) Kordopatis, G., Binney, J., Gilmore, G., et al. 2015, MNRAS, 447, 3526
  * Kordopatis et al. (2022) Kordopatis, G., Schultheis, M., McMillan, P. J., et al. 2022, arXiv e-prints, arXiv:2206.07937
  * Kovtyukh et al. (2022) Kovtyukh, V., Lemasle, B., Bono, G., et al. 2022, MNRAS, 510, 1894
  * Lallement et al. (2019) Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
  * Lindegren (2018) Lindegren, L. 2018, Tech. rep., GAIA-C3-TN-LU-LL-124, http://www.rssd. esa.int/doc-fetch.php?id=3757412
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell, D. & Kalnajs, A. J. 1972, MNRAS, 157, 1
  * Maiolino & Mannucci (2019) Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3
  * Matteucci (2021) Matteucci, F. 2021, A&A Rev., 29, 5
  * Minchev et al. (2018) Minchev, I., Anders, F., Recio-Blanco, A., et al. 2018, MNRAS, 481, 1645
  * Minchev et al. (2012) Minchev, I., Famaey, B., Quillen, A. C., et al. 2012, A&A, 548, A126
  * Okuda et al. (1977) Okuda, H., Maihara, T., Oda, N., & Sugiyama, T. 1977, Nature, 265, 515
  * Pastorelli et al. (2019) Pastorelli, G., Marigo, P., Girardi, L., et al. 2019, MNRAS, 485, 5666
  * Pedicelli et al. (2009) Pedicelli, S., Bono, G., Lemasle, B., et al. 2009, A&A, 504, 81
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104
  * Poggio et al. (2018) Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21
  * Re Fiorentin et al. (2019) Re Fiorentin, P., Lattanzi, M. G., & Spagna, A. 2019, MNRAS, 484, L69
  * Re Fiorentin et al. (2021) Re Fiorentin, P., Spagna, A., Lattanzi, M. G., & Cignoni, M. 2021, ApJ, 907, L16
  * Recio-Blanco et al. (2014) Recio-Blanco, A., de Laverny, P., Kordopatis, G., et al. 2014, A&A, 567, A5
  * Recio-Blanco et al. (2022) Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2022
  * Roberts (1969) Roberts, W. W. 1969, ApJ, 158, 123
  * Roškar et al. (2012) Roškar, R., Debattista, V. P., Quinn, T. R., & Wadsley, J. 2012, MNRAS, 426, 2089
  * Sánchez-Menguiano et al. (2016) Sánchez-Menguiano, L., Sánchez, S. F., Kawata, D., et al. 2016, ApJ, 830, L40
  * Sánchez-Menguiano et al. (2020) Sánchez-Menguiano, L., Sánchez, S. F., Pérez, I., et al. 2020, MNRAS, 492, 4149
  * Schönrich & Binney (2009a) Schönrich, R. & Binney, J. 2009a, MNRAS, 396, 203
  * Schönrich & Binney (2009b) Schönrich, R. & Binney, J. 2009b, MNRAS, 399, 1145
  * Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829
  * Sellwood & Binney (2002) Sellwood, J. A. & Binney, J. J. 2002, MNRAS, 336, 785
  * Shen & Zheng (2020) Shen, J. & Zheng, X.-W. 2020, Research in Astronomy and Astrophysics, 20, 159
  * Spitoni et al. (2022a) Spitoni, E., Aguirre Børsen-Koch, V., Verma, K., & Stokholm, A. 2022a, A&A, 663, A174
  * Spitoni et al. (2019) Spitoni, E., Cescutti, G., Minchev, I., et al. 2019, A&A, 628, A38
  * Spitoni et al. (2022b) Spitoni, E., Recio-Blanco, A., de Laverny, P., et al. 2022b, arXiv e-prints, arXiv:2206.12436
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287
  * Vergely et al. (2022) Vergely, J. L., Lallement, R., & Cox, N. L. J. 2022, arXiv e-prints, arXiv:2205.09087
  * Vogt et al. (2017) Vogt, F. P. A., Pérez, E., Dopita, M. A., Verdes-Montenegro, L., & Borthakur, S. 2017, A&A, 601, A61
  * Williams et al. (2022) Williams, T. G., Kreckel, K., Belfiore, F., et al. 2022, MNRAS, 509, 1303
  * Xiang et al. (2022) Xiang, M., Rix, H.-W., Ting, Y.-S., et al. 2022, A&A, 662, A66
  * Zari et al. (2021) Zari, E., Rix, H. W., Frankel, N., et al. 2021, A&A, 650, A112

##  Appendix A Data quality cuts and selection

We select Gaia DR3 sources with atmospheric parameters (available through the
`gaiadr3.astrophysical_parameters` table), radial velocities, and 5-parameters
astrometric solution having quality index `ruwe` $<1.4$ (Lindegren 2018). In
addition, _Gaia_ duplicated sources are excluded from the sample. As for the
_GSP-Spec_ parameters, we reject objects based on their metallicity
uncertainty and _GSP-Spec_ quality flags (Recio-Blanco et al. 2022) reporting
on the degree of biases from line broadening and radial velocity errors
affecting [M/H], flux noise, extrapolation on stellar atmospheric parameters,
as well as KM-type stars.

The definition of the working sample is as follows:

     [M/H]_unc<0.5 & vbroadM<2 & vradM<2 &
     fluxNoise<4 & extrapol<3 & KMtypestars<2
     astrometric_params_solved = 31 & ruwe < 1.4
     & duplicated_source == false

From the Gaia archive the following query provides the data set employed:

    SELECT g.*, ap.teff_gspspec,
    ap.logg_gspspec, ap.mh_gspspec
    FROM gaiadr3.gaia_source AS g
    INNER JOIN gaiadr3.astrophysical_parameters AS ap
    ON g.source_id = ap.source_id
    WHERE (
    ((mh_gspspec_upper-mh_gspspec_lower)<0.25)
    AND ((flags_gspspec LIKE ‘__0%’)
    OR (flags_gspspec LIKE ‘__1%’))
    AND ((flags_gspspec LIKE ‘_____0%’)
    OR (flags_gspspec LIKE ‘_____1%’))
    AND ((flags_gspspec LIKE ‘______0%’)
    OR (flags_gspspec LIKE ‘______1%’)
    OR (flags_gspspec LIKE ‘______2%’)
    OR (flags_gspspec LIKE ‘______3%’))
    AND ((flags_gspspec LIKE ‘_______0%’)
    OR (flags_gspspec LIKE ‘_______1%’)
    OR (flags_gspspec LIKE ‘_______2%’))
    AND ((flags_gspspec LIKE ‘____________0%’)
    OR (flags_gspspec LIKE ‘____________1%’))
    ) AND g.astrometric_params_solved = 31
    AND g.ruwe < 1.4

Additionally, objects without $d_{Geo}$ distances were discarded.

For the obtained sample, we calculate the Galactocentric coordinates
$R,\phi,Z$ and the corresponding velocities $V_{R},V_{\phi},V_{Z}$ following
the same conventions and parameters adopted in Gaia Collaboration, Recio-
Blanco et al. (2022). We apply a cut —Z—¡0.75 kpc, as we are mainly interested
in selecting disc stars. To this end, we also apply a kinematical cut based on
the Toomre diagram (see for example Bensby et al. 2003, 2014; Re Fiorentin et
al. 2019, 2021; Gaia Collaboration, Recio-Blanco et al. 2022; Gaia
Collaboration, Creevey et al. 2022) in order to remove possible halo
contaminants:

$\sqrt{V_{R}^{2}+V_{Z}^{2}+(V_{\phi}-V_{0})^{2}}<210\,\,\,\mbox{km\leavevmode\nobreak\
s}^{-1},$ (1)

where $V_{0}=238.5$ km s-1 is the velocity of the local standard of rest at
the Sun’s position, based on Gravity Collaboration et al. (2021) and Schönrich
et al. (2010).

Finally, we select different populations of RGB stars in the log($g$)-$T_{\rm
eff}$ plane of the Kiel diagram. This is accomplished by selecting stars in
the three boxes shown in the left panel of Fig. 1, that is:

    logg < 1.5 & logg > 0.5 &
    (logg > (coeff*teff + interc_left)) &
    (logg < (coeff*teff + interc_right))

where `coeff`= 0.00192 $\rm{dex\,K}^{-1}$ is the adopted slope, chosen to
follow the natural inclination of the RGB branch (as done also for the
selection of the massive sample in Gaia Collaboration, Recio-Blanco et al.
2022), whereas the intercept of the two inclined lines, delimiting the
selected regions on the left and right side, are: `interc_left` = -8.3 +
$\Delta$ dex and `interc_right` = -7.3 + $\Delta$ dex, where $\Delta$ = 0, 0.5
and 1 for Sample A, B and C, respectively.

##  Appendix B Sample characterization

Figure 4:  PARSEC isochrones for different metallicites and ages, compared to
the regions of the Kiel diagram selected here for Sample A, B and C.  Figure
5:  Distribution of azimuthal velocities $V_{\phi}$ for Sample A, B, and C.
The histogram is normalized so that the area under the histogram integrates to
1 for all three samples.

Here we performed some tests to further explore the content of Sample A, B and
C. First of all, we compared our selected areas of the Kiel diagram with the
prediction from the PARSEC isochrones for different metallicities (Bressan et
al. 2012; Chen et al. 2014, 2015; Tang et al. 2014; Pastorelli et al. 2019).
Figure 4 shows the isochrones for 10 Gyr, 1 Gyr, 500 Myr, 100 Myr, 30 and 20
Myr, for both solar metallicity [M/H]=0 dex and [M/H]=-0.5 dex. It should be
noted that effective temperatures from _GSP-Spec_ have been found to be in
agreement with those from the isochrones (see Section 10.2 of Recio-Blanco et
al. 2022), and no important bias as been found with the literature. As we can
see in Fig 4, in the region covered by Sample A, only isochrones for 100 Myr
or younger are present. Of course, in real data, we do expect some
contamination from older stars, given that uncertainties will tend to blur the
distribution. Nevertheless, the fact that Sample B and C contain stars
typically older than those in Sample A is supported by the location of the
isochrones in the Kiel diagram.

Further indications can be also found empirically. Although ages can be also
inferred for stars (e.g. Kordopatis et al. 2022), here we simply use stellar
kinematics as a proxy for the typical age of the sample. For Sample A, B and
C, respectively, the velocity dispersions are:
$\sigma_{V_{Z}}=(19.5,21.8,23.6)$ km s-1, $\sigma_{V_{R}}=(28.6,35.1,38.9)$ km
s-1, $\sigma_{V_{\phi}}=(26.6,26.9,29.0)$ km s-1. Figure 5 shows the
distribution of the azimuthal velocities $V_{\phi}$. As we can see, the
distribution of stars in Sample A presents a prominent peak at high azimuthal
velocity ($\sim$ 230-240 km s-1). On the contrary, Sample B and Sample C
present a broader distribution, with a large fraction of stars at low
$V_{\phi}$, as would be expected as a consequence of the asymmetric drift for
a kinematically hot stellar population.

##  Appendix C Bivariate smoothing

For a given position (X,Y) in the Galactic plane, we calculate the local
smoothed metallicity:

$\langle\,{\rm[M/H]}\,\rangle_{\rm
loc}\,(X,Y)=\frac{\sum^{N}_{i}{\rm[M/H]}_{i}\,K\left(\frac{X-X_{i}}{h_{X}}\right)\,K\left(\frac{Y-Y_{i}}{h_{Y}}\right)}{\sum^{N}_{i}K\left(\frac{X-X_{i}}{h_{X}}\right)\,K\left(\frac{Y-Y_{i}}{h_{Y}}\right)}\quad,$
(2)

where $N$ is the total number of stars in the sample, $X_{i}$ and $Y_{i}$ are
the X- and Y-coordinates of the $i$-th star of the sample, $K$ is a Gaussian
kernel

$K\biggl{(}\frac{X-X_{i}}{h_{X}}\biggr{)}=\frac{1}{\sqrt{2\pi}\,h_{X}}\exp{\Biggl{(}-\frac{1}{2}\Biggl{(}\frac{(X-X_{i})}{h_{X}}\Biggr{)}^{2}\Biggr{)}}$
(3)

with a scale-length $h_{X}$ (and similarly for the Y-coordinate). In this
work, local smoothed metallicity (see Fig. 1, right panels) is shown using a
scale-length of 0.175 kpc. As discussed in the main text, the metallicity
excess is calculated by taking the difference of the local metallicity (i.e.
Equation 2, assuming a scale-length of 0.175 kpc) and the metallicity on a
large scale (i.e. again using Equation 2, but now adopting a scale-length 5
times larger than the one used for the local metallicity).

##  Appendix D Azimuthal variations

Figure 6 shows a projection of the observed metal-rich spiral features when
the Galactic plane is cut into Galactocentric rings. For each ring, we show
the azimuthal variations for both Sample A and Sample C, calculated as the
median metallicity as a function of azimuth, after subtracting the median
metallicity of the sample in each ring. We note that the azimuthal variations
depend on the adopted sample, being more prominent for Sample A; typically,
they do not exceed 0.08 dex. The amplitude and location of the observed peaks
depend on how the rings are selected with respect to the spiral arms’
geometry.

Figure 6:  _Impact of the observed chemical inhomogeneities on azimuthal
variations._ _Upper panel_ : Map of the mean metallicity in the Galactic plane
for Sample A, with overlaid rings of constant Galactocentric radius at
$R=6,7,8,9,10$ and $11$ kpc. _Lower panels_ : median metallicity for sample A
and C (filled/open dots, respectively) as a function of Galactic azimuth
$\phi$ for different rings in $R$ (as specified in each panel), after the
median metallicity of the stars for each ring is subtracted.

##  Appendix E Additional tests

Here we present some additional tests that we perform to verify the robustness
of our results.

First, we checked the quality of the astrometric measurements in our selected
stars. For Sample A, 99.7 % of the stars have a parallax signal-to-noise
greater than 5. For those stars, as a further check, we constructed a map
using distances inferred by the inverse of the parallax, obtaining consistent
results.

As an additional test, we selected only the stars satisfying the criteria on
_Gaia_ DR3 astrophysical parameters listed for the high-quality sample in Gaia
Collaboration, Recio-Blanco et al. (2022). In this case, the chemical
signature of the spiral arms is still present, but the entire map is shifted
toward more metal-rich values, as expected (given that more metal-rich stars
are more likely to be included in the high-quality sample).

To test the possible contribution from Asymptotic Giant Branch (AGB) stars, we
cross-matched our three samples with the long-period variables catalog
available in _Gaia_ DR3 `gaiadr3.vari_long_period_variable`. We found that,
for all three samples, the percentage of AGB stars from that catalog is less
than 0.2 %.

Our maps are also quite robust against the selection of the regions in the
Kiel diagram. Indeed, a change in the shape of the selection leave the maps
almost unchanged - although, of course, a significant shift in $T_{\rm eff}$
(or log($g$)) can make an impact, as discussed in the main text.

Finally, we test other spiral arms maps. The spiral arms contours used in the
main text are from Poggio et al. (2021), which relies on distances calculated
as described in Poggio et al. (2018). As a test, we also used the spiral arms
overdensity map from Gaia Collaboration, Drimmel et al. (2022), where the
$d_{Geo}$ distances are adopted from Bailer-Jones et al. (2021), as well as
the one based the sample from Zari et al. (2021) (which used astro-kinematic
distances), always obtaining consistent results.