# Joint reconstructions of growth and expansion histories from stage-IV
surveys with minimal assumptions. II. Modified gravity and massive neutrinos.

Rodrigo Calderón*  —<EMAIL_ADDRESS>Korea Astronomy and Space Science
Institute, Daejeon 34055, South Korea  Laboratoire Charles Coulomb, Université
de Montpellier & CNRS, 34095 Montpellier, France  Benjamin L’Huillier*  —
<EMAIL_ADDRESS>Department of Physics and Astronomy, Sejong University,
Seoul 05006, South Korea  David Polarski*  —<EMAIL_ADDRESS>Laboratoire Charles Coulomb, Université de Montpellier & CNRS, 34095
Montpellier, France  Arman Shafieloo*  —<EMAIL_ADDRESS>Korea Astronomy
and Space Science Institute, Daejeon 34055, South Korea  University of Science
and Technology, Daejeon 34113, South Korea  Alexei A. Starobinsky*  —
<EMAIL_ADDRESS>L. D. Landau Institute for Theoretical Physics RAS,
Chernogolovka, Moscow region 142432, Russia  Bogolyubov Laboratory of
Theoretical Physics, Joint Institute for Nuclear Research, Dubna 141980,
Russia

###### Abstract

Based on a formalism introduced in our previous work, we reconstruct the
phenomenological function $G_{\rm eff}(z)$ describing deviations from general
relativity (GR) in a model-independent manner. In this alternative approach,
we model $\mu\equiv G_{\mathrm{eff}}/G$ as a Gaussian process and use
forecasted growth-rate measurements from a stage-IV survey to reconstruct its
shape for two different toy-models. We follow a two-step procedure: (i) we
first reconstruct the background expansion history from supernovae (SNe) and
baryon acoustic oscillation (BAO) measurements; (ii) we then use it to obtain
the growth history $f\sigma_{8}$, that we fit to redshift-space distortions
(RSD) measurements to reconstruct $G_{\mathrm{eff}}$. We find that upcoming
surveys such as the Dark Energy Spectroscopic Instrument (DESI) might be
capable of detecting deviations from GR, provided the dark energy behavior is
accurately determined. We might even be able to constrain the transition
redshift from $G\to G_{\mathrm{eff}}$ for some particular models. We further
assess the impact of massive neutrinos on the reconstructions of
$G_{\mathrm{eff}}$ (or $\mu$) assuming the expansion history is given, and
only the neutrino mass is free to vary. Given the tight constraints on the
neutrino mass, and for the profiles we considered in this work, we recover
numerically that the effect of such massive neutrinos do not alter our
conclusions. Finally, we stress that incorrectly assuming a $\Lambda$CDM
expansion history leads to a degraded reconstruction of $\mu$, and/or a non-
negligible bias in the ($\Omega_{\text{m},0}$,$\sigma_{8,0}$)-plane.

##  I Introduction

Addressing the late-time accelerated phase of expansion of the Universe
remains a major challenge for fundamental physics [1, 2]. Though most
observations to date are in agreement with the standard (concordance) model of
cosmology ($\Lambda$CDM), alternative explanations for dark energy (DE)—other
than a cosmological constant $\Lambda$—are still up for debate (see e.g. [3]).
In particular, modifying the laws of gravity (beyond Einstein’s GR) at large-
scales remains a tantalizing possibility [4, 5]. Besides the exact nature of
the dark energy (DE) component and its (effective) equation of state,
additional modifications come with the properties of the relativistic degrees
of freedom, notably the neutrino sector. Interestingly, despite the wide class
of modified-gravity (MG) scenarios explored in the last decades, observations
seem to suggest that GR remains our best description of gravitational
interactions, where dark energy is in the form of a cosmological constant in
the Einstein field equations. For example, the detection of GW 170817,
together with its electromagnetic counterpart GRB 170817A [6], implies that
gravitational waves travel at the speed of light—ruling out a large subclass
of Horndeski models predicting a tensor speed $c_{T}\neq c$ at the present
epoch [7]. Hence the detection of gravitational waves (GW) has added stringent
constraints on modified gravity models in addition to local constraints. Note
that a viable cosmic expansion history can give additional strong constraints,
for example, on $f(R)$ models [8].111Viable cosmological models of the present
Universe in $f(R)$ gravity satisfying these constraints were independently
constructed soon after that paper in [9, 10, 11]. At the phenomenological
level, most modified theories of gravity predict a time (and possibly scale)
dependent _effective_ gravitational coupling $G_{\rm eff}(z)~{}$[12, 13]
entering the equation for the growth of perturbations. Thus, detecting a
deviation from Newton’s constant would be a smoking gun for physics beyond
$\Lambda$CDM and even beyond GR.

Let us present now the basic formalism of our approach, starting with the
background. We consider here spatially flat Friedmann-Lemaître-Robertson-
Walker universes with

$h^{2}(z)\equiv
H^{2}/H_{0}^{2}=\Omega_{\text{m},0}(1+z)^{3}+(1-\Omega_{\text{m},0})f_{\mathrm{DE}}(z)~{},$
(1)

where $f_{\mathrm{DE}}=\rho_{\rm DE}(z)/\rho_{\rm DE}(z=0)$. While the second
term in (1) becomes generically subdominant in the past for viable
cosmologies, this has to be enforced explicitly at high redshifts (where no
data are available) once we use Gaussian processes in order to reconstruct
$h(z)$ [14]. We stress further that the parameter $\Omega_{\text{m},0}$ refers
to _clustered_ dustlike matter only. The second term of (1) is more general
than the compact notation suggests, see the discussion given in [14]. We turn
now to the perturbations. We use the following conventions and notations [12]
(see also e.g. [15]) in the conformal Newtonian gauge, where the perturbed
FLRW metric is described by ($c=1$)

$\displaystyle\mathrm{d}s^{2}=-(1+2\phi)\mathrm{d}t^{2}+(1-2\psi)a^{2}\mathrm{d}\bm{x}^{2},$
(2)

where $\phi$ and $\psi$ are the Bardeen potentials. Phenomenologically, on
subhorizon scales, in many modified gravity models the departure from the
standard perburbations growth in GR is encoded in the modified Poisson
equation [12] (see also e.g. [16, 15, 17])

$\displaystyle\nabla^{2}\phi=4\pi G_{\text{eff}}(a,\bm{k})~{}\rho\equiv 4\pi
G\mu(a,\bm{k})~{}\rho~{}.$ (3)

GR corresponds obviously to $\mu\equiv 1$. The relation between the Bardeen
potentials is expressed as follows

$\displaystyle\phi\equiv\eta(a,\bm{k})~{}\psi~{},$ (4)

the two potentials are generically unequal in these models. The subhorizon
modes are essentially affected by $\mu$ as is explicit from Eq. (5) given
below, while super horizon modes are affected by both $\mu$ and $\eta$ [16].
In this work, given the datasets considered, we restrict our attention to
$\mu$ (see e.g. [18, 19] for constraints on $\eta$). In what follows, we will
use $G_{\mathrm{eff}}$ and $\mu$ interchangeably, since $\mu$ is just
$G_{\mathrm{eff}}$ in units of $G$. The growth of dustlike subhorizon matter
perturbations in the quasi-static Approximation (QSA) is then governed by [12]

$\ddot{\delta}+2H\dot{\delta}=4\pi G\,\rho\,\delta\,\mu(z,\bm{k}),$ (5)

where $\delta\equiv\delta\rho/\rho$ is the density contrast of dustlike
matter. For modes of cosmological interest, the $k$-dependence of $\mu$ is
often mild and can be neglected in a first approach [20, 21, 22]—see e.g. [23,
24, 25, 26] for current and future constraints on the scaledependence of
$\mu$. Note that this is certainly the case for the unscreened scalar-tensor
model considered in [12]. We will restrict ourselves here to phenomenological
models where $\mu$ or $G_{\text{eff}}$ is scale independent.

The above equation can be rewritten in terms of the growth factor
$f\equiv\delta^{\prime}/\delta$, to give

$\displaystyle
f^{\prime}+\left(f+2+\frac{h^{\prime}}{h}\right)f-\frac{3}{2}\Omega_{\text{m}}(z)\mu(z)$
$\displaystyle=0~{},$ (6)

where a prime stands for derivative with respect to $N\equiv\ln{a}$. From an
observational standpoint, redshift space distortions (RSD) provide us with
growth rate measurements of the quantity

$f\sigma_{8}\equiv\frac{\sigma_{8,0}}{\delta_{0}}f\delta=\frac{\sigma_{8,0}}{\delta_{0}}\delta^{\prime},\;\;\;\text{
with }\;\;\;\delta_{0}=\delta(z=0)~{}.$ (7)

We remind that the quantities $\Omega_{i}$ appearing in (1) and (6) are
defined in the standard way as in GR with the help of Newton’s constant $G$.

In this work, we will use the synergy between geometrical background probes
(type Ia supernovae [SN] and baryon acoustic oscillations [BAO]) and growth
measurements from RSD to constrain the phenomenological function $\mu(z)$
describing the departures from GR. While current analysis pipelines rely on
various assumptions (namely, $\Lambda$+GR) when extracting the cosmological
information from large-scale structure, in particular the BAO and RSD
measurements, we expect that our results will remain essentially unaffected
when such effects are taken into account.

The paper is organized as follows. We start by describing in detail the
methodology and the data used in Sec. II. In Sec. III, we apply the method to
simulated RSD data generated with $\mu\neq 1$ in both idealistic and realistic
scenarios and further discuss the implications of the results. We also comment
on the effects of incorrectly assuming a $\Lambda$CDM expansion history on the
reconstructions in Sec. III.3. In Sec. IV, we consider separately the
inclusion of massive neutrinos.

##  II Method and Data

###  II.1 Models and Mock Data

For the data, we generate mock $f\sigma_{8}$ measurements for a (stage-IV)
DESI-like survey following Tables 2.3–2.7 in [27] (covering $14$ K $\deg^{2}$)
and for different behaviors of $G_{\mathrm{eff}}$ that we aim to reconstruct.
Namely, we consider an $f(R)$-inspired _bumplike_ profile (which we refer to
simply as “bump”) and a smooth _steplike_ transition (“dip” hereafter) in the
recent past towards the weak gravity regime ($G_{\mathrm{eff}}<G$), see e.g.
[28, 29].222Indeed, both such profiles can occur in viable cosmological models
in $f(R)$ gravity, see [30] in particular, especially in the case of
oscillations around phantom divide [31]. These two profiles are treated purely
phenomenologically here,333Ref. [32] presented a concrete MG model with
similar profiles for $\mu(z)$ considered in this work (or rather their
reflections along the $\mu=1$ axis), which could simultaneously ease the
$H_{0}$ and $\sigma_{8}$ tensions. indeed viable $f(R)$ theories are actually
screened and allow $G_{\rm eff,0}$ to deviate from $G$ today. Nonetheless, due
to the $k$-dependence of $\mu$ which we do not discuss here, cosmic scales
smaller than some critical scale would experience a boost in their growth in
the recent past.

In the case of the dip, we consider it mainly to assess whether such profiles
can be accurately reconstructed using our model-independent approach. Note in
this context that a decreasing $\mu$ is impossible in massless scalar-tensor
models [33]. To summarize, these hybrid profiles allow us to test our
reconstruction independently of any theoretical prior.

Figure 1: Schematic representation of our analysis pipeline. We sample the
(log) of the hyperparameters ($\sigma_{f},\ell_{f},z_{c}$) to obtain samples
of $\mu(z)$, which we use – together with a given expansion history $h(z)$ and
a value of $\sigma_{8}$ – to compute the growth history $f\sigma_{8}(z)$ that
we compare to DESI-like RSD measurements to “reconstruct” $\mu(z)$. The
expansion history is reconstructed in a model-independent fashion, using the
methodology introduced in [14].

The behaviors of the phenomenological functions $\mu^{\rm fid}(z)$ used to
generate the data are depicted by the dashed lines in the upper panel of Fig.
2, while the corresponding growth $f\sigma_{8}(z)$ evolutions are shown in the
lower panel. We also make use of stage-IV SN+BAO data to determine the
background expansion history $h(z)$ without relying on a specific parametric
model, as explained in Sec. III.2. The fiducial background used to generate
the data is a Chevallier-Polarski-Linder (CPL) model [34, 35], extensively
discussed in [14] with

$\theta^{\mathrm{fid}}=\\{\Omega_{\mathrm{m,0}}^{\mathrm{fid}}=0.28,w_{0}^{\mathrm{fid}}=-0.95,w_{a}^{\mathrm{fid}}=0.3,\\\
h_{\mathrm{fid}}=0.7,\sigma_{8,0}^{\mathrm{fid}}=0.81\\},$ (8)

where
$H_{0}=100\,h\,$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\,\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$.
More details on the background-only (SN+BAO) mock data can also be found in
[14]. Already at this stage, let us note that modified theories of gravity can
lead to a modified Chandrasekhar mass (with $m_{\rm ch}\sim
G_{\mathrm{eff}}^{-3/2}$ [36]), relevant for SNeIa analyses, which can affect
the absolute magnitude [e.g. $\Delta M=\frac{15}{4}\log\mu(z)$ in scalar-
tensor theories444Note however that this theoretical correction can be even
smaller, if the stretch correction is taken into account [linderprivate] [37,
38])]] and hence the distance measurements obtained from such standard candles
[39, 40, 41]. This effect has even been proposed as a possible explanation for
the mismatch between early and late-time measurements of the Hubble constant
$H_{0}$, see e.g. [42, 43, 44, 45, 46, 47, 48]. However, for our purposes, we
neglect these effects and assume the $h(z)$ measurements obtained from SNe are
independent of $\mu$ in the current analysis. The inclusion of these effects
for a specific model might be the subject of future works.

###  II.2 The method

To explore possible modifications of gravity at late times, we model
$G_{\mathrm{eff}}(z)$ as a Gaussian process555We do not delve into the details
of Gaussian process modeling here, instead we refer the reader to our previous
work [14] and the excellent review [49] for more. Note that in this work,
unlike common notations in the GP literature, $\mu$ refers to the
phenomenological function appearing in (3), and the mean of the GP is denoted
by $\bar{f}$. (GP) centered around Newton’s constant $G$, such that

$\mu(z;\sigma_{f},\ell_{f},z_{c})=\begin{cases}\mathcal{GP}(\bar{f}(z)=1,k(\sigma_{f},\ell_{f})),&\text{for}\
z<z_{c}\\\ 1,&\text{for}\ z\geq z_{c}\end{cases}$ (9)

so that we recover GR at large-$z$. We “pretrain” our GP with the following
theoretical priors:

$\displaystyle\mu(0)$ $\displaystyle=1\pm 0.05,$ (10a)
$\displaystyle\mu(z_{c})$ $\displaystyle=1\pm 10^{-10},$ (10b)
$\displaystyle\mu^{\prime}(0)$ $\displaystyle=\mu^{\prime}(z_{c})=0\pm
10^{-10}.$ (10c)

These conditions allow us to smoothly recover $G_{\mathrm{eff}}=G$ above a
certain $z_{c}$ and at $z=0$, while exploring possible departures from GR at
intermediate redshifts $0.1<z<10$ (see e.g. [50, 51, 52, 46, 19, 53, 54, 55]
for other approaches). Recovering $G_{\mathrm{eff}}(0)=G$ is not strictly
necessary (see our discussion at the beginning of this section), but from a
technical point of view it can help guide our reconstructions at very low $z$
where we are volume limited and uncertainties become quite large. Furthermore,
when dealing with real data, we do not know the true behavior of $\mu$, and
whether the underlying model is screened or not, hence the two representative
behaviors at $z=0$ chosen for our profiles. It is comforting to find that the
first condition does not alter the reconstruction of the second profile around
$z=0$ as illustrated by the blue curves in Fig. 2.

We use a squared exponential kernel given by

$k(z,z^{\prime};\sigma_{f},\ell_{f})=\sigma_{f}^{2}\,e^{-(z-z^{\prime})^{2}/2\ell_{f}^{2}},$
(11)

where $\sigma_{f}$ and $\ell_{f}$ determine the amplitude and typical length
scale of the correlations, respectively [49].

In a Bayesian spirit, we give flat uninformative and wide priors to the
cosmological and (hyper)parameters, listed in Table 1. We sample the parameter
space using Markov chain Monte Carlo (MCMC) methods, as implemented in emcee
[56, 57]. At each step in the MCMC, we draw a sample of $\mu(N=\ln{a})\equiv
G_{\mathrm{eff}}/G\sim\mathcal{GP}(1,K)$, characterized by
$(\sigma_{f},\ell_{f},z_{c})$, and solve the growth equation for a given value
of $\sigma_{8,0}$ and expansion $h(z)$, to obtain a solution
$f\sigma_{8}(z)$—see the diagram in Fig. 1—prior to any comparison with the
data. Note that the parameter $z_{c}$ also enters prior to any computation of
the likelihood and is irrespective of the data points. In other words, this
can be seen as _forward-modeling_ , rather than training the GP with the data
in the usual sense. Thus, we rely on the maximization of the following
likelihood function

$\ln{\mathcal{L}}=-\frac{1}{2}~{}\chi^{2}=-\frac{1}{2}~{}\bm{r}^{\intercal}\cdot\mathbf{C}^{-1}\cdot\bm{r},$
(12)

where $\bm{r}=f\sigma_{8}^{\rm GP}-f\sigma_{8}^{\rm data}$ is the residual
vector and $\mathbf{C}^{-1}$ is the inverse of the covariance matrix. The
growth history $f\sigma_{8}^{\rm GP}$ is obtained by solving the Eq. (6) for
each “pre-trained” sample of $\mu(z)$ drawn from Eq. (9). Those samples of
$\mu(z)$ retracing a similar shape to $\mu^{\rm fid}$ will yield a better fit
to growth data, and thus will be statistically favored in the long run.
Averaging over a large number of realizations gives the median shape of
$\mu(z)$ and $95\%\,(2\sigma)$ confidence intervals around it. This is along
the lines of what was done in [14] to reconstruct $f_{\mathrm{DE}}$, but this
time we also include conditions on the derivatives of the GP, to smoothly
recover the form in Eq. (9), following the formalism described in Appendix A.

Parameter | $\sigma_{8,0}$ | $\log_{10}\sigma_{f}$ | $\log_{10}\ell_{f}$ | $\log_{10}z_{c}$
---|---|---|---|---
Prior | $[0.5,1.2]$ | $[-3,0.5]$ | $[-1,0.2]$ | $[-1,1]$
Table 1: Uniform priors for the parameters used in the MCMC analyses.

##  III Results and Discussions

Figure 2: Reconstructions of $G_{\mathrm{eff}}$ in the idealistic case where
the background $h(z)$ and amplitude of fluctuations $\sigma_{8,0}$ are
perfectly known. Solid lines and shaded regions correspond to the median,
$68\%$ and $95\%$ confidence intervals around it, respectively. Dashed lines
correspond to the fiducial cosmologies generating the DESI-like (RSD) data.
The redshift $z_{c}$ of the transition to GR, as well as the hyperparameters
$\sigma_{f}$ and $\ell_{f}$ appearing in (9) are nonetheless free parameters
to be determined by the data. Both of these reconstructions detect deviations
from GR ($\mu=1$) at more than $2\sigma$ for $z\sim 1$.

###  III.1 Ideal case: Background is perfectly known

We first consider the idealistic case where the background expansion history
is perfectly known. In other words, we fix $\Omega_{\text{m},0}$ and
$\sigma_{8,0}$ to their fiducial values and further assume that the dark
energy evolution is known $f_{\mathrm{DE}}(z)=f_{\mathrm{DE}}^{\rm fid}$.
Although this is far from being a realistic scenario, it allows us to test our
method and quantify the uncertainties purely coming from the modifications of
gravity, encoded in $G_{\mathrm{eff}}$.

The posterior distributions for $\mu(z)$ assuming perfect knowledge of $h(z)$
and $\sigma_{8,0}$ are shown in Fig. 2. If the background (and the amplitude
of fluctuations $\sigma_{8,0}$) are perfectly known, the RSD data alone are
enough to perform an accurate (within $2\sigma$) reconstruction of the
underlying theory of gravity, _i.e._ $G_{\mathrm{eff}}(z)$. In the next
subsection, we take a more realistic approach, where only minimal assumptions
on the background are made666We only assume a flat FLRW universe, and that the
Hubble rate is a sum of a matter term and an “effective” DE component [14] and
$h(z)$ is purely determined from the data.

###  III.2 Realistic case: $\sigma_{8,0}$ free—$\Omega_{\text{m},0}$ and
$f_{\mathrm{DE}}(z)$ determined by SN+BAO

In this section, instead of assuming a parametric form for $h(z)$, we use the
reconstructed expansion history as determined by SN+BAO data. In practice,
this amounts to obtaining an expansion history $h(z)$ from the samples of
$f_{\mathrm{DE}}$ and calculating angular and luminosity distances which are
then fitted to the data, as explained in [14]. The degeneracies between
$\sigma_{8,0},\Omega_{\text{m},0}$ and $G_{\mathrm{eff}}$ make it very hard to
say something about the underlying theory of gravity, given the quality of the
data and, in particular, when all parameters are free to vary. To circumvent
this issue, we assume a single expansion history, as determined solely by the
data. More specifically, the expansion history $h(z)$, along with the value of
$\Omega_{\text{m},0}$–needed for solving the growth equation (6)–is the median
of all the realizations drawn from the Markov SN+BAO chains,777The posterior
distributions correspond to the blue contours shown in Fig. 6 of Calderón  _et
al._ [14]. obtained in [14]. Indeed, it was shown in [14] that our method is
able to capture a large class of DE models, even those where the contribution
from DE is not negligible at high-$z$. Our reconstruction of $h(z)$ is
accurate to $\lesssim 1\%$ across the entire redshift range of interest–see
Fig. 3. The amplitude of the fluctuations, $\sigma_{8,0}$, now becomes a free
parameter, and we sample the full parameter space
$\theta=\\{\sigma_{8,0},\log_{10}\sigma_{f},\log_{10}\ell_{f},\log_{10}z_{c}\\}$
in the range given by Table 1. In Fig. 4, we show the reconstructions when
using the median of $h(z)$ and median $\Omega_{\text{m},0}$ from the SN+BAO
chains. As expected, the uncertainties in the reconstructions increase with
respect to those in Fig.2, as $\sigma_{8,0}$ is now a free parameter that is
somewhat degenerate with $G_{\mathrm{eff}}$, allowing for more flexibility in
the samples of $G_{\mathrm{eff}}$ drawn at each step in MCMC.

Figure 3: _Top:_ Reconstruction of the DE evolution $f_{\mathrm{DE}}(z)$.
_Bottom:_ Relative (percentage) errors in the background reconstructions from
forecasted SN+BAO measurements. The orange line corresponds to the true
fiducial background in (8), while gray lines depict the reconstructed median,
68 and 95% confidence levels around it. Dashed-black line correspond to
$\Lambda$CDM’s best-fit ($f_{\mathrm{DE}}=1$, $\Omega_{\text{m},0}=0.3103$) to
SN+BAO data.

The advantage of taking this approach is that we do not make any assumption on
the evolution of DE, and we are able to effectively reconstruct any expansion
history directly from the data, by reconstructing $f_{\mathrm{DE}}(z)$.
Moreover, this disentangles the uncertainties coming from the growth evolution
$f\sigma_{8}(z)$ and those coming from the background expansion $h(z)$. This
also allows us to point down a value for $\Omega_{\text{m},0}$, which is of
course anticorrelated with $\sigma_{8,0}$, which is in turn anticorrelated
with $G_{\mathrm{eff}}$. Thus, allowing for more constraining power on the
quantity of interest $\mu(z)$ from RSD alone. The two-dimensional posteriors
of the quantity $\mu$ at two different redshifts $z=0$ and $z=1.4$ are shown
in Fig. 5. At $z=1.4$, where most of the constraining power of RSD
measurements lies, the bumplike posteriors in red exclude GR ($\mu=1$, in
dashed) at $>2\sigma$, while the posteriors for the diplike profile in blue
are marginally consistent with GR at the $2\sigma$ level. At low redshift,
because of the large uncertainties in $f\sigma_{8}$, the posteriors are much
broader and provide a $\sim 20\%$ constraint on $\mu(z=0)$. We note that the
study of peculiar velocities using SNIa from ZTF and LSST can potentially
improve the measurements of the growth at very low-$z$ by a factor of 2 with
respect to DESI [58]—see also [59] for other interesting constraints using
gravitational waves and galaxies’ peculiar velocities. Interestingly, because
the redshift $z_{c}$ in (9) of the transition from $G\to G_{\mathrm{eff}}$ is
a free parameter, our method allows us to constrain when the departures from
GR start taking place (see Fig. 8 and the discussions in Appendix B). For the
particular profiles considered in this work, the posterior distribution of
$z_{c}$ is quite peaked, and we have a “detection” of a transition from $G\to
G_{\mathrm{eff}}$ in both cases, as seen from Fig. 8. The corresponding
constraints on $z_{c}$ are given in Table 2

Figure 4: Realistic case where $\sigma_{8,0}$ is allowed to vary, and the
background $h(z)$ is determined by SN+BAO (gray lines in Fig. 3). The fiducial
cosmologies used to generate the $f\sigma_{8}(z)$ measurements are shown by
the dashed lines. Despite having larger error confidence intervals with
respect to the idealistic case in Fig. 2, both of these reconstructions are
still able to rule out GR at more than $2\sigma$ at $z\sim 1$. Figure 5:
Marginalized posterior distributions of the relevant cosmological parameters,
when using our model-independent reconstructions of $h(z)$, shown in gray in
Fig. 3 [where the unknown function $f_{\mathrm{DE}}(z)$ is reconstructed in a
fully model-independent way and $\Omega_{\text{m},0}$ is fixed to the median
of all possible values obtained from the SN+BAO chains— _cf._ Sec. III.2].

###  III.3 Incorrectly assuming a $\Lambda$CDM background

Figure 6: Reconstructions of $\mu(z)$ when assuming the best-fit
$\Lambda$CDM’s expansion history, with $\Omega_{\text{m},0}=0.31$. Incorrectly
assuming a $\Lambda$CDM background leads to biased determinations of
$\sigma_{8,0}$ and a degraded reconstruction of $\mu(z)$, despite being
perfectly consistent with $f\sigma_{8}(z)$, as can be seen from the lower
panel.

Cosmological observations suggest that dark energy is in the form of a
cosmological constant $\Lambda$. Because of its simplicity and agreement with
observations, it remains the standard model of cosmology today. Thus, most
cosmological analyses are done within the $\Lambda$CDM framework, which might
lead to biased reconstructions if DE is _not constant_ , as for the fiducial
cosmology considered here. In this section, we explore the effects of
incorrectly assuming a $\Lambda$CDM background expansion history in the
reconstructions of $\mu(z)$. In other words, we fit a $\Lambda$CDM model to
the SN+BAO mock data described before and find the corresponding best-fit
value for $\Omega_{\text{m},0}$ (and thus
$\Omega_{\Lambda,0}=1-\Omega_{\text{m},0}$). We remind the reader that the
mock data are generated from a time-evolving CPL dark energy model, given by
Eq. (8). We then use this expansion history to solve for the perturbations and
reconstruct $\mu(z)$, as explained in the previous sections. The black dashed
lines in Fig. 3 show the best-fit $\Lambda$CDM expansion history (with
$\Omega_{\text{m},0}^{\Lambda\rm CDM,bf}=0.3103^{+0.0025}_{-0.0024}$),
compared to the fiducial one with $\Omega_{\text{m},0}^{\rm fid}=0.28$ in
orange (hence, representing a $\sim 12\sigma$ bias in the fractional matter
density). Despite having almost identical $H(z)$, the differences in the DE
evolution $f_{\mathrm{DE}}(z)$ and biased $\Omega_{\text{m},0}$ translate into
a degraded reconstruction of $G_{\mathrm{eff}}$, shown in Fig. 6—to be
compared with Fig. 4. We also find that the inferred value of $\sigma_{8,0}$
can be biased $\sigma_{8,0}\sim 0.78$ vs. $\sigma_{8,0}^{\rm fid}=0.81$ (which
corresponds to a $\sim 1.2\sigma$ bias in the inferred amplitude of
fluctuations) for the case of the dip (in blue)—see Table 2. As understood
from our previous work [14], from the background-only (SN+BAO) standpoint, the
lack of DE at high-$z$ is compensated by higher values of
$\Omega_{\text{m},0}$, which translates into lower values of $\sigma_{8,0}$
(or lower $G_{\mathrm{eff}}<G_{\mathrm{eff}}^{\rm fid}$) to maintain the
agreement with growth-rate measurements of $f\sigma_{8}(z)$. This is a perfect
example of what might happen if one incorrectly assumes DE is constant, the
background expansion history might be consistent with the geometrical probes
(SN+BAO), but a tension might appear in the amplitude of fluctuations
$\sigma_{8,0}$ inferred from LSS observables. Despite the bias in the
cosmological parameters $\Omega_{\text{m},0}$ and $\sigma_{8,0}$—and for the
specific cases of $\mu(z)$ considered here—the reconstructions are still able
to capture the main trends in $\mu(z)$.

Finally, let us note that for the steplike transition in blue, the reason why
the reconstructions deviate somehow from the fiducial $\mu^{\rm fid}(z)$ (in
dashed) at very low-$z$ is because of our theoretical prior
$G_{\mathrm{eff}}(z=0)\simeq G$, which tends to draw our GP samples back to 1.
We stress that this prior does not need to be imposed, as we do not
necessarily have $G_{\mathrm{eff}}(z=0)\simeq G$ in most MG theories. We have
in mind here theories without screening mechanisms that do require
$G_{\mathrm{eff}}\simeq G$ today to satisfy local constraints, e.g. [60].
Despite this prior, because of the large uncertainties in RSD measurements at
$z\sim 0$, our reconstructions are still able to capture (within $2\sigma$)
the true fiducial $\mu^{\rm fid}_{\rm Dip}$.

##  IV Effect of massive neutrinos

In this section, we consider universes containing massive neutrinos. We want
to investigate how well our reconstruction of $\mu$ fares in their presence.
It is well known that free-streaming species with nonzero mass (here massive
neutrinos) lead to a suppression of gravitational clustering on scales below a
characteristic scale, corresponding to their free-streaming length. Hence,
while massive neutrinos contribute to the universe’s expansion in the same way
as usual dustlike matter (corresponding to $\Omega_{m}$), they are absent from
the driving term in the matter perturbations growth. Hence we have in front of
us a situation where the parameter $\Omega_{{\rm m},0}$ does not represent all
dustlike components at low redshifts. Indeed, one cannot distinguish massive
neutrinos from dustlike matter purely from geometric probes at low $z$. In
this case, the splitting in (1), while sensible theoretically, is somewhat
ambiguous regarding expansion data if we have no additional information on
$\Omega_{{\rm m},0}$ or $\Omega_{\nu,0}$. This ambiguity however gets broken
once we consider the perturbations growth. In a first step, we assume the
presence of massive neutrinos and we work with equation (14) below [instead of
(1)]. So, while we reconstruct $\mu$ as a Gaussian process, we assume the
background expansion is known up to the two parameters $\Omega_{{\rm m},0}$
and $m_{\nu}$. Here however, we have only one free parameter left. Indeed, in
this section we fix the present relative energy density
$\Omega_{\text{m},0}^{\text{tot}}$ of all components which behave like dust at
low $z$, namely,

$\displaystyle\Omega_{\text{m},0}^{\text{tot}}\equiv\Omega_{\text{m},0}+\Omega_{\nu,0}=\Omega_{\text{cdm},0}+\Omega_{\text{b},0}+\Omega_{\nu,0}~{},$
(13)

where $\Omega_{\text{cdm},0}$, $\Omega_{\text{b},0}$, and $\Omega_{\nu,0}$ are
the present relative densities of cold dark matter, baryons, and massive
neutrinos respectively. Note that the couple of parameters $\left(\Omega_{{\rm
m},0},m_{\nu}\right)$ and $\left(\Omega_{{\rm m},0}^{\rm tot},m_{\nu}\right)$
carry the same information.

We assume now that $h^{2}(z)$ is given by

$\displaystyle h^{2}(z)$
$\displaystyle=~{}\Omega_{\text{m},0}~{}(1+z)^{3}+\Omega_{\Lambda,0}$ (14)
$\displaystyle+\Omega_{\gamma,0}~{}(1+z)^{4}\left(1+0.2271\,\frac{N_{\text{eff}}}{3}\,\sum_{i}f_{\nu}\left(\frac{m_{\nu_{i}}}{T_{\nu}}\right)\right),$

where $f_{\nu}(y)\simeq(1+(Ay)^{p})^{1/p}$ is a fit provided in Ref. [61],
with $A=\frac{180\zeta(3)}{7\pi^{4}}$ and $p=1.83$; where $\zeta$ is Riemann’s
$\zeta$-function. This fitting function $f_{\nu}$ describes the evolution from
the relativistic behavior when $m_{\nu}\ll T_{\nu}$ ($T_{\nu}\sim a^{-1}$) to
the nonrelativistic regime when we have eventually $m_{\nu}\gg T_{\nu}$. Like
in (1), the first term appearing in (14) corresponds to the fractional amount
of matter that clusters. In order to test our reconstruction in the presence
of massive neutrinos, it is more relevant to consider universes sharing
identical $\Omega_{{\rm m},0}^{{\rm tot}}$ rather than identical $\Omega_{{\rm
m},0}$, but with different $\Omega_{{\rm m},0}$, or equivalently different
neutrino masses $m_{\nu}$. Clearly, the parameters $\Omega_{{\rm m,0}}^{\rm
tot}$ and $m_{\nu}$, completely define the background expansion (14).

The driving term in the perturbations growth equation depends on the
combination $\mu~{}\Omega_{m}$. Hence for modified gravity and in the presence
of massive neutrinos, this combination is modified at low redshifts as follows

$\displaystyle G\Omega_{\rm m}^{{\rm tot}}\to G_{\mathrm{eff}}\Omega_{\rm m}$
$\displaystyle=G\Omega_{\rm m}^{{\rm
tot}}~{}\mu\left(1-\frac{\Omega_{\nu}}{\Omega_{\rm m}^{{\rm tot}}}\right)$
$\displaystyle\approx
0.965\frac{m_{\nu}}{$0.5\text{\,}\mathrm{e}\mathrm{V}$}h_{70}^{-2}\mu~{}G\Omega_{\rm
m}^{{\rm tot}},$ (15)

where we evidently have $\Omega_{\rm m}^{{\rm tot}}=\Omega_{\rm m}$ in the
absence of massive neutrinos, and
$h_{70}=H_{0}/$70\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\,{Mpc}^{-1}$$.
For the values we take here, the change comes essentially from modified
gravity.

Here, we forecast the future surveys’ potential to reconstruct the coupling
strength $\mu(z)$ in the presence of massive neutrinos and purely from RSD
measurements of $f\sigma_{8}(z)$. As before, we generate mock data from a
fiducial model; this time we choose a ($\Lambda$CDM$\nu$) cosmology containing
2 massless and 1 massive neutrinos, with $m_{\nu}^{\rm fid}=0.5\;\rm eV$.
Although this mass is larger than what is currently allowed by cosmological
observations888Cosmological constraints are indirect and somewhat model
dependent, unlike ground-based experiments. [62, 63], it is still within the
allowed mass range probed by terrestrial experiments, which constrain
$m_{\nu}^{2}\equiv\Sigma_{i}\left\lvert
U_{ei}\right\rvert^{2}m_{i}^{2}=0.26^{+0.34}_{-0.34}\;\rm eV^{2}$, yielding an
upper bound on the electron (anti)-neutrino mass $m_{\nu}<0.8\;\rm eV$ at
$90\%$ confidence level [64]999Note that masses of usual and sterile neutrinos
$m_{\nu}\sim 1\,\rm eV$ are well possible in viable $f(R)$ cosmological models
[65, 66].. The rest of the cosmological parameters are fixed to Planck’s best-
fit values. Because of the growth suppression from such a massive neutrino,
the normalization of the matter power spectrum $P_{\rm m}(k,z=0)$,
characterized by $\sigma_{8,0}$, is now $\sigma_{8,0}^{\rm fid}\simeq 0.73$,
lower than in the previous sections (where $\sigma_{8,0}$ was fixed to
$\sigma_{8,0}^{\rm fid}=0.81$).
In what follows, we assume that this normalization ($\sigma_{8,0}^{\rm
fid}=0.73$, as obtained for $\mu=1$) is the same for all profiles of
$G_{\mathrm{eff}}$. Although the actual normalization of the $P_{\rm
m}(k,z=0)$ would indeed depend on the theory of gravity, we generate mock data
for different profiles of $\mu$ from the same value of $\sigma_{8,0}$. We
stress that this choice is arbitrary, as we are dealing with simulated data
and we are interested in assessing whether the theory of gravity $\mu(z)$ and
$\sigma_{8,0}$ are accurately recovered by our model-independent
reconstructions, which do not know anything about the underlying theory that
generates the data.

We then sample the parameters
$\theta=\\{\sigma_{8,0},m_{\nu},\log_{10}\sigma_{f},\log_{10}\ell_{f},\log_{10}z_{c}\\}$,
with $m_{\nu}\in[0,1]$ to see the impact of a varying neutrino mass on the
reconstructions of $\mu(z)$. The posterior distributions for the relevant
cosmological parameters are shown on Fig. 7. Although we sample $m_{\nu}$, we
show the posterior distributions for the derived parameter
$\Omega_{\text{m},0}$, as it corresponds to the driving term for the growth in
the right-hand side of Eq.(6) and the actual neutrino mass is unconstrained.
The value of $\sigma_{8,0}$ is anticorrelated with the reconstructions of
$\mu$, mainly seen in the $(\sigma_{8,0},\mu(z=1.4))$-plane. Large deviations
from GR, up to $\mu(z=1.4)\sim 1.8$ can be achieved, provided that the
amplitude of fluctuations $\sigma_{8,0}$ is low ($\sigma_{8,0}\sim 0.65$). A
slight (negative) correlation between $\Omega_{\text{m},0}$ and $\sigma_{8,0}$
is also obtained, as expected. The enhanced suppression of growth (due to
larger mass $m_{\nu}$, hence smaller
$\Omega_{\text{m},0}=\Omega_{\text{m},0}^{\rm tot}-\Omega_{\nu,0}$) needs to
be compensated by larger values of $\sigma_{8,0}$, to maintain the agreement
with $f\sigma_{8}$ measurements. Despite these correlations, the
reconstructions of $\mu(z)$ remain accurate, and does not seem to be affected
by a varying neutrino mass (other than increasing the uncertainties in the
reconstructions, due to an additional free parameter). The fiducial value for
$\sigma_{8,0}$, shown as a dashed vertical line in Fig. 7, is also accurately
recovered.

Finally, let us note that we separately tested our reconstructions in the
presence of massive neutrinos without assuming the functional form of $h(z)$,
given by (14) but using instead the (reconstructed) _effective_
$f_{\mathrm{DE}}(z)$ in Eq. (1), which captures the effect of relativistic
species [14]. Our conclusions remain unaltered, but no information on the
neutrino mass can be obtained.

Figure 7: Marginalized posterior distributions for the parameters in the
presence of massive neutrinos. This figure is the same as Fig. 5, but this
time assuming the background is known $h(z)$ (up to 1 free parameter
$m_{\nu}$) given by Eq. (14), including relativistic species and when the
neutrino mass is free to vary (_cf._ Sec. IV).

##  V Conclusions

In a companion paper Calderón  _et al._ [14], we jointly reconstructed the
growth and expansion histories inside GR directly from the data and using
minimal assumptions. We showed that our framework is able to capture a wide
variety of behaviors in the DE component. In this work, we extend our
methodology to include possible modifications of gravity at late times, as
encoded by the function $G_{\mathrm{eff}}(z)$ appearing in the (modified)
Poisson equation. We illustrate the efficiency of our method in reconstructing
different theories of gravity by reconstructing two phenomenological _shapes_
of $\mu(z)\equiv G_{\mathrm{eff}}/G$. As an example, we consider a “bump” and
a smooth transition (“dip”) towards the weak gravity regime in the recent
past. We used the reconstructed $h(z)$ from background-only data, as obtained
in [14] in order to fit $f\sigma_{8}(z)$ to RSD mock data, thereby
constraining $\mu(z)$ using minimal assumptions. We also explore the effects
of incorrectly assuming a $\Lambda$CDM background. In both cases, the fiducial
$\mu(z)$ is within the $1\sigma$ confidence intervals of our reconstructions,
if the background is accurately determined, and within $2\sigma$ if we
incorrectly assume the $\Lambda$CDM’s best-fit $h(z)$. Finally, we explored
the impact of massive neutrinos on the reconstructions of $\mu(z)$. To
summarize, let us list a few important results.

  * •

If the background is _given_ (Fig. 2), or _accurately reconstructed_ from
SN+BAO (Fig. 3), our reconstructions of $G_{\mathrm{eff}}(z)$ are able to
distinguish both fiducial $\mu$-profiles from GR at $\gtrsim 2\sigma$ (see
Figs. 2 and 4).

  * •

Incorrectly assuming a $\Lambda$CDM expansion (with the best-fit
$\Omega_{\text{m},0}$ to background probes) can lead to _biased/degraded
reconstructions_ (red-shaded regions in Fig. 6) and/or _biased estimations_ of
the amplitude of fluctuations $\sigma_{8,0}$ (see Table 2). This is despite
the perfect agreement with $f\sigma_{8}(z)$ measurements, as shown in the
lower panel of Fig. 6.

  * •

The posterior distributions for the hyperparameters _clearly_ show the _need_
for a deviation from the mean $\bar{f}=1$, _i.e._ GR is _not_ a good
description of the data. This is understood because the marginalized contours
in Fig. 8 suggest the data are not consistent with vanishing values of
$\sigma_{f}$, i.e. the posterior does not extend to $\sigma_{f}\to 0$, and
therefore _require_ deviations from the considered mean function.
Interestingly, the redshift of the transition $z_{c}$ is also not compatible
with small values of $z_{c}$, and we have a “detection” on when this
transition from $G\to G_{\mathrm{eff}}$ happens; seen as a clear bump in Fig.
8.

In this work, we used forecasted (stage-IV) SN+BAO data to reconstruct the DE
evolution $f_{\mathrm{DE}}(z)$—which determines the expansion history
$h(z)$—and separately reconstructed $\mu(z)$ using DESI-like $f\sigma_{8}(z)$
measurements for two different toy models of $G_{\mathrm{eff}}$. We expect our
methodology to hold for essentially any (viable) form of $G_{\mathrm{eff}}$.
We showed that for both profiles considered in this work, the reconstructions
are able to detect the deviations from GR at $\gtrsim 2\sigma$ in the redshift
range $0.5\lesssim z\lesssim 1.5$ where DESI’s (RSD) constraining power lies.
The inclusion of external data sets, such as the (modified) luminosity
distance of gravitational waves $d_{L}^{\rm GW}(z)$ [67] or the Integrated
Sachs-Wolfe effect (ISW) seen in the temperature anisotropies of the cosmic
microwave background (CMB) in cross-correlation with LSS surveys would provide
interesting (model-independent) constraints on the allowed deviations from GR
[68]. Moreover, we note that the effect of massive neutrinos would be tracked
more accurately if we allow for a scale-dependent growth. We leave such
extensions for future work.

## Acknowledgements

We thank Eric Linder for comments on the draft. B.L. acknowledges the support
of the National Research Foundation of Korea (NRF-2019R1I1A1A01063740 and
NRF-2022R1F1A1076338) and the support of the Korea Institute for Advanced
Study (KIAS) grant funded by the government of Korea. A.S. would like to
acknowledge the support by National Research Foundation of Korea
NRF2021M3F7A1082053, and the support of the Korea Institute for Advanced Study
(KIAS) grant funded by the government of Korea. A.A.S. was partly supported by
the Project No. 0033-2019-0005 of the Russian Ministry of Science and Higher
Education.

##  Appendix A GAUSSIAN PROCESS WITH OBSERVATIONS ON THE DERIVATIVES

Figure 8: Marginalized posterior distributions for the relevant parameters from the RSD chains. The background expansion history used in the analysis is fixed to the median $h(z)$ obtained from the SN+BAO chains, shown as a gray line in Fig. 3. Model | $\sigma_{8,0}$ | $\log_{10}z_{c}$ | $\ln{\mathcal{L}}$
---|---|---|---
Bump | $0.814^{+0.022}_{-0.032}$ | $0.524^{+0.342}_{-0.092}$ | $-10.15^{+0.78}_{-1.48}$
Dip | $0.806^{+0.017}_{-0.013}$ | $0.79^{+0.15}_{-0.40}$ | $-10.36^{+0.96}_{-1.46}$
Bump $(\Lambda\rm CDM)$ | $0.812^{+0.013}_{-0.021}$ | $0.86^{+0.10}_{-0.42}$ | $-10.35^{+0.94}_{-1.97}$
Dip $(\Lambda\rm CDM)$ | $0.784^{+0.023}_{-0.015}$ | $-0.17^{+0.20}_{-0.34}$ | $-10.5^{+1.0}_{-1.3}$
Table 2: Marginalized constraints (median and $68\%$ confidence levels) on the
relevant parameters. We compare the effects of incorrectly assuming
$\Lambda$CDM (bottom-two rows) vs the reconstructed $h(z)$ from SN+BAO in the
upper-two rows.In the last column, we report the log-likelihood values
$\ln\mathcal{L}\equiv-\frac{1}{2}\chi^{2}$.

In this section, we describe a less common use of Gaussian process when we
also observe the derivative of the function $f$ to be reconstructed [49, 69].
We note that in this section, $f$ denotes a general function, not the growth
rate. In our case, $f=\mu(z)$. In addition to observations of $y$, we also
“observe” $y^{\prime}=f^{\prime}(x)+\varepsilon$, where,

$\varepsilon\sim\mathcal{N}(0,C_{y^{\prime}}),$ (16)

is a Gaussian noise and $C_{y^{\prime}}$ is the covariance of the derivatives.

We further assume that $y$ and $y^{\prime}$ are uncorrelated. Therefore, the
vector

$\begin{bmatrix}y\\\ y^{\prime}\\\ f\\\ f^{\prime}\end{bmatrix}$ (17)

is jointly Gaussian, and the posterior predictive distribution can be
calculated using

$\displaystyle\begin{bmatrix}f\\\
f^{\prime}\end{bmatrix}|y,y^{\prime},X,X_{*}\sim\mathcal{N}\left(\begin{bmatrix}\bar{f}\\\
\bar{f}^{\prime}\end{bmatrix};\begin{bmatrix}A-CB^{-1}C^{T}\end{bmatrix}\right),$
(18)

where the mean is

$\displaystyle\begin{bmatrix}\bar{f}\\\ \bar{f}^{\prime}\end{bmatrix}$
$\displaystyle=CB^{-1}\begin{bmatrix}y-\mu_{y}\\\
y^{\prime}-\mu_{y^{\prime}}\end{bmatrix},$ (19)

and the covariance matrix is given by

$\displaystyle A$ $\displaystyle=\begin{bmatrix}K_{**}&K_{**}^{01}\\\
K_{**}^{10}&K_{**}^{11}\\\ \end{bmatrix}\in\mathbb{M}_{2n_{*},2n_{*}},$ (20a)
$\displaystyle B$ $\displaystyle=\begin{bmatrix}K+C_{y}&K^{01}\\\
K^{10}&K^{11}+C_{y^{\prime}}\end{bmatrix}\in\mathbb{M}_{n+n^{\prime}},$ (20b)
$\displaystyle C^{T}$ $\displaystyle=\begin{bmatrix}K_{*}&K_{*}^{01}\\\
K_{*}^{10}&K_{*}^{11}\end{bmatrix}\in\mathbb{M}_{n+n^{\prime},2n_{*}},$ (20c)

where

$\displaystyle K$ $\displaystyle=k(X,X),$ (21a) $\displaystyle K_{*}$
$\displaystyle=k(X,X_{*}),$ (21b) $\displaystyle K_{**}$
$\displaystyle=k(X_{*},X_{*}),$ (21c)

and for any matrix $X$,

$\displaystyle X^{i,j}$ $\displaystyle=\frac{\partial^{i+j}X}{\partial
X^{i}\partial X^{j}}.$ (22)

The subscript $*$ denote the set of points $X$ where the observations are
done.101010In our analysis, the observations are done in redshift, so that
$X_{*}=z_{*}$. However, this formalism is general and can be applied to any
input variable $X$. This formalism allows us to impose theoretical priors on
the samples of $\mu(z)$ and its derivative $\mu^{\prime}(z)$ to smoothly
recover the expected GR behavior at early times [see Eq. (9)].

##  Appendix B Distribution of the hyperparameters

Inspecting the posterior distributions of the hyperparameters, shown in Fig.
8, can yield additional information on the $G_{\mathrm{eff}}$ reconstructions
and put interesting constraints on the departures from GR. First, let us note
that the inferred value of $\sigma_{8,0}$ is unbiased in both cases, when the
evolution of DE $f_{\mathrm{DE}}(z)$ is reconstructed using our model
independent approach [14]. This is not the case when one (incorrectly) assumes
a $\Lambda$CDM expansion history (see Table 2). Second, both the bump and dip
reconstructions seem to require a deviation from the mean function
$\bar{f}=\mu=1$ (_i.e._ GR), as the posteriors of $\log_{10}\sigma_{f}$ are
not compatible with $\sigma_{f}\to 0$. This suggest that GR is not a good
description of the growth $f\sigma_{8}(z)$ history and that the data requires
extra flexibility, as encoded by the GP kernel in Eq. (11). Lastly, the
posteriors of $\log_{10}z_{c}$ seem to peak at the redshift $z_{c}\sim 3$
where the departures from GR actually takes place (depicted by the vertical
dashed line in Fig. 8).

## References

  * Weinberg [1989] S. Weinberg, The Cosmological Constant Problem, Rev. Mod. Phys. 61, 1 (1989).
  * Sahni and Starobinsky [2000] V. Sahni and A. A. Starobinsky, The Case for a positive cosmological Lambda term, Int. J. Mod. Phys. D 9, 373 (2000), arXiv:astro-ph/9904398 .
  * Perivolaropoulos and Skara [2021] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update (2021), arXiv:2105.05208 [astro-ph.CO] .
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified gravity and cosmology, Physics Reports 513, 1–189 (2012).
  * Tsujikawa [2010] S. Tsujikawa, Modified gravity models of dark energy, Lect. Notes Phys. 800, 99 (2010), arXiv:1101.0191 .
  * 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, and et al., Gravitational waves and gamma-rays from a binary neutron star merger: GW 170817 and GRB 170817A, The Astrophysical Journal 848, L13 (2017).
  * Baker  _et al._ [2017] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Strong Constraints on Cosmological Gravity from GW170817 and GRB 170817A, Phys. Rev. Lett. 119, 251301 (2017), arXiv:1710.06394 .
  * Amendola  _et al._ [2007] L. Amendola, D. Polarski, and S. Tsujikawa, Are $f(R)$ dark energy models cosmologically viable?, Phys. Rev. Lett. 98, 131302 (2007).
  * Hu and Sawicki [2007] W. Hu and I. Sawicki, Models of f(R) cosmic acceleration that evade solar system tests, Phys. Rev. D 76, 064004 (2007), arXiv:0705.1158 .
  * Appleby and Battye [2007] S. A. Appleby and R. A. Battye, Do consistent $F(R)$ models mimic General Relativity plus $\Lambda$?, Phys. Lett. B 654, 7 (2007), arXiv:0705.3199 [astro-ph] .
  * Starobinsky [2007] A. A. Starobinsky, Disappearing cosmological constant in f(R) gravity, JETP Lett. 86, 157 (2007), arXiv:0706.2041 [astro-ph] .
  * Boisseau  _et al._ [2000] B. Boisseau, G. Esposito-Farèse, D. Polarski, and A. A. Starobinsky, Reconstruction of a Scalar-Tensor Theory of Gravity in an Accelerating Universe, Phys. Rev. Lett. 85, 2236 (2000), arXiv:gr-qc/0001066 .
  * Felice  _et al._ [2011] A. D. Felice, T. Kobayashi, and S. Tsujikawa, Effective gravitational couplings for cosmological perturbations in the most general scalar–tensor theories with second-order field equations, Physics Letters B 706, 123 (2011).
  * Calderón  _et al._ [2022] R. Calderón, B. L’Huillier, D. Polarski, A. Shafieloo, and A. A. Starobinsky, Joint reconstructions of growth and expansion histories from stage-IV surveys with minimal assumptions: Dark energy beyond $\Lambda$, Phys. Rev. D 106, 083513 (2022), arXiv:2206.13820 .
  * Simpson  _et al._ [2013] F. Simpson, C. Heymans, D. Parkinson, C. Blake, M. Kilbinger, J. Benjamin, T. Erben, H. Hildebrandt, H. Hoekstra, T. D. Kitching, Y. Mellier, L. Miller, L. Van Waerbeke, J. Coupon, L. Fu, J. Harnois-Déraps, M. J. Hudson, K. Kuijken, B. Rowe, T. Schrabback, E. Semboloni, S. Vafaei, and M. Velander, CFHTLenS: testing the laws of gravity with tomographic weak lensing and redshift-space distortions, MNRAS 429, 2249 (2013), arXiv:1212.3339 .
  * Pogosian  _et al._ [2010] L. Pogosian, A. Silvestri, K. Koyama, and G.-B. Zhao, How to optimally parametrize deviations from general relativity in the evolution of cosmological perturbations, Phys. Rev. D 81, 104023 (2010), arXiv:1002.2382 .
  * Pogosian and Silvestri [2016] L. Pogosian and A. Silvestri, What can cosmology tell us about gravity? Constraining Horndeski gravity with $\Sigma$ and $\mu$, Phys. Rev. D 94, 104014 (2016).
  * Marta Pinho  _et al._ [2018] A. Marta Pinho, S. Casas, and L. Amendola, Model-independent reconstruction of the linear anisotropic stress $\eta$, J. Cosmology Astropart. Phys 2018, 027 (2018), arXiv:1805.00027 .
  * Alestas  _et al._ [2022] G. Alestas, L. Kazantzidis, and S. Nesseris, Machine learning constraints on deviations from general relativity from the large scale structure of the Universe (2022), arXiv:2209.12799 .
  * Hojjati  _et al._ [2012] A. Hojjati, L. Pogosian, A. Silvestri, and S. Talbot, Practical solutions for perturbed $f(R)$ gravity, Phys. Rev. D 86, 123503 (2012).
  * Hojjati  _et al._ [2014] A. Hojjati, L. Pogosian, A. Silvestri, and G.-B. Zhao, Observable physical modes of modified gravity, Phys. Rev. D 89, 083505 (2014).
  * Bose  _et al._ [2015] S. Bose, W. A. Hellwing, and B. Li, Testing the quasi-static approximation in $f(R)$ gravity simulations, J. Cosmology Astropart. Phys 2015, 034 (2015).
  * Baker  _et al._ [2014] T. Baker, P. G. Ferreira, C. D. Leonard, and M. Motta, New gravitational scales in cosmological surveys, Phys. Rev. D 90, 124030 (2014), arXiv:1409.8284 .
  * Johnson  _et al._ [2016] A. Johnson, C. Blake, J. Dossett, J. Koda, D. Parkinson, and S. Joudaki, Searching for modified gravity: scale and redshift dependent constraints from galaxy peculiar velocities, MNRAS 458, 2725 (2016), arXiv:1504.06885 .
  * Garcia-Quintero  _et al._ [2020] C. Garcia-Quintero, M. Ishak, and O. Ning, Current constraints on deviations from general relativity using binning in redshift and scale, J. Cosmology Astropart. Phys 2020, 018 (2020).
  * Denissenya and Linder [2022] M. Denissenya and E. V. Linder, Constraining scale dependent growth with redshift surveys, J. Cosmology Astropart. Phys 2022, 029 (2022), arXiv:2208.10508 .
  * Aghamousa  _et al._ [2016] A. Aghamousa _et al._ (DESI), The DESI Experiment Part I: Science,Targeting, and Survey Design (2016), arXiv:1611.00036 .
  * Kazantzidis and Perivolaropoulos [2021] L. Kazantzidis and L. Perivolaropoulos, $\sigma$8 Tension. Is Gravity Getting Weaker at Low z? Observational Evidence and Theoretical Implications, in  _Modified Gravity and Cosmology; An Update by the CANTATA Network_, edited by E. N. Saridakis, R. Lazkoz, V. Salzano, P. V. Moniz, S. Capozziello, J. Beltrán Jiménez, M. De Laurentis, and G. J. Olmo (2021) pp. 507–537, arXiv:1907.03176 .
  * Gannouji  _et al._ [2021] R. Gannouji, L. Perivolaropoulos, D. Polarski, and F. Skara, Weak gravity on a $\Lambda$CDM background (2021), arXiv:2011.01517 .
  * Motohashi  _et al._ [2010] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Phantom boundary crossing and anomalous growth index of fluctuations in viable f(R) models of cosmic acceleration, Prog. Theor. Phys. 123, 887 (2010), arXiv:1002.1141 [astro-ph.CO] .
  * Motohashi  _et al._ [2011] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Future Oscillations around Phantom Divide in f(R) Gravity, JCAP 06, 006, arXiv:1101.0744 [astro-ph.CO] .
  * Saridakis [2023] E. N. Saridakis, Solving both $H_{0}$ and $\sigma_{8}$ tensions in $f(T)$ gravity (2023) arXiv:2301.06881 [gr-qc] .
  * Gannouji  _et al._ [2018] R. Gannouji, L. Kazantzidis, L. Perivolaropoulos, and D. Polarski, Consistency of modified gravity with a decreasing $G_{\rm eff}(z)$ in a $\Lambda$CDM background, Phys. Rev. D 98, 104044 (2018), arXiv:1809.07034 .
  * Chevallier and Polarski [2001] M. Chevallier and D. Polarski, Accelerating Universes with Scaling Dark Matter, Int. J. Mod. Phys. D 10, 213 (2001), gr-qc/0009008 .
  * Linder [2003] E. V. Linder, Exploring the Expansion History of the Universe, Physical Review Letters 90, 091301 (2003), astro-ph/0208512 .
  * Amendola  _et al._ [1999] L. Amendola, S. Corasaniti, and F. Occhionero, Time variability of the gravitational constant and Type Ia supernovae (1999), arXiv:astro-ph/9907222 .
  * Gaztañaga  _et al._ [2001] E. Gaztañaga, E. García-Berro, J. Isern, E. Bravo, and I. Domínguez, Bounds on the possible evolution of the gravitational constant from cosmological type-Ia supernovae, Phys. Rev. D 65, 023506 (2001), arXiv:astro-ph/0109299 .
  * Gannouji  _et al._ [2006] R. Gannouji, D. Polarski, A. Ranquet, and A. A. Starobinsky, Scalar tensor models of normal and phantom dark energy, J. Cosmology Astropart. Phys 2006, 016 (2006), arXiv:astro-ph/0606287 .
  * Wright and Li [2018] B. S. Wright and B. Li, Type Ia supernovae, standardizable candles, and gravity, Phys. Rev. D 97, 083505 (2018), arXiv:1710.07018 .
  * Zhao  _et al._ [2018] W. Zhao, B. S. Wright, and B. Li, Constraining the time variation of Newton’s constant G with gravitational-wave standard sirens and supernovae, J. Cosmology Astropart. Phys 2018, 052 (2018), arXiv:1804.03066 .
  * Ballardini and Finelli [2022] M. Ballardini and F. Finelli, Type Ia supernovae data with scalar-tensor gravity, Phys. Rev. D 106, 063531 (2022), arXiv:2112.15126 .
  * Benevento  _et al._ [2020] G. Benevento, W. Hu, and M. Raveri, Can Late Dark Energy Transitions Raise the Hubble constant?, Phys. Rev. D 101, 103517 (2020), arXiv:2002.11707 [astro-ph.CO] .
  * Braglia  _et al._ [2020] M. Braglia, M. Ballardini, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, K. Koyama, and D. Paoletti, Larger value for $H_{0}$ by an evolving gravitational constant, Phys. Rev. D 102, 023529 (2020), arXiv:2004.11161 .
  * Ballardini  _et al._ [2020] M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A. Starobinsky, and C. Umiltà, Scalar-tensor theories of gravity, neutrino physics, and the $H_{0}$ tension, JCAP 10, 044, arXiv:2004.14349 [astro-ph.CO] .
  * Alestas  _et al._ [2021] G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, w -M phantom transition at zt¡0.1 as a resolution of the Hubble tension, Phys. Rev. D 103, 083517 (2021), arXiv:2012.13932 .
  * Pogosian  _et al._ [2021] L. Pogosian, M. Raveri, K. Koyama, M. Martinelli, A. Silvestri, and G.-B. Zhao, Imprints of cosmological tensions in reconstructed gravity (2021), arXiv:2107.12992 [astro-ph.CO] .
  * Perivolaropoulos and Skara [2022] L. Perivolaropoulos and F. Skara, A reanalysis of the latest SH0ES data for $H_{0}$: Effects of new degrees of freedom on the Hubble tension (2022), arXiv:2208.11169 [astro-ph.CO] .
  * Alestas  _et al._ [2022] G. Alestas, L. Perivolaropoulos, and K. Tanidis, Constraining a late time transition of $G_{\rm eff}$ using low-z galaxy survey data (2022), arXiv:2201.05846 [astro-ph.CO] .
  * Rasmussen and Williams [2006] C. Rasmussen and C. Williams,  _Gaussian Processes for Machine Learning_, Adaptative computation and machine learning series (University Press Group Limited, 2006).
  * Nesseris  _et al._ [2011] S. Nesseris, C. Blake, T. Davis, and D. Parkinson, The WiggleZ Dark Energy Survey: constraining the evolution of Newton’s constant using the growth rate of structure, J. Cosmology Astropart. Phys 2011, 037 (2011), arXiv:1107.3659 [astro-ph.CO] .
  * Espejo  _et al._ [2019] J. Espejo, S. Peirone, M. Raveri, K. Koyama, L. Pogosian, and A. Silvestri, Phenomenology of Large Scale Structure in scalar-tensor theories: joint prior covariance of $w_{\textrm{DE}}$, $\Sigma$ and $\mu$ in Horndeski, Phys. Rev. D 99, 023512 (2019), arXiv:1809.01121 [astro-ph.CO] .
  * Raveri  _et al._ [2021] M. Raveri, L. Pogosian, K. Koyama, M. Martinelli, A. Silvestri, G.-B. Zhao, J. Li, S. Peirone, and A. Zucca, A joint reconstruction of dark energy and modified growth evolution (2021), arXiv:2107.12990 [astro-ph.CO] .
  * Ruiz-Zapatero  _et al._ [2022a] J. Ruiz-Zapatero, D. Alonso, P. G. Ferreira, and C. Garcia-Garcia, Impact of the Universe’s expansion rate on constraints on modified growth of structure, Phys. Rev. D 106, 083523 (2022a), arXiv:2207.09896 [astro-ph.CO] .
  * Ruiz-Zapatero  _et al._ [2022b] J. Ruiz-Zapatero, C. García-García, D. Alonso, P. G. Ferreira, and R. D. P. Grumitt, Model-independent constraints on $\Omega$m and H(z) from the link between geometry and growth, MNRAS 512, 1967 (2022b), arXiv:2201.07025 [astro-ph.CO] .
  * Heisenberg  _et al._ [2022] L. Heisenberg, H. Villarrubia-Rojo, and J. Zosso, Can late-time extensions solve the $H_{0}$ and $\sigma_{8}$ tensions? (2022), arXiv:2202.01202 [astro-ph.CO] .
  * Goodman and Weare [2010] J. Goodman and J. Weare, Ensemble samplers with affine invariance, Communications in Applied Mathematics and Computational Science 5, 65 (2010).
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publications of the Astronomical Society of the Pacific 125, 306–312 (2013).
  * Graziani  _et al._ [2020] R. Graziani, M. Rigault, N. Regnault, P. Gris, A. Möller, P. Antilogus, P. Astier, M. Betoule, S. Bongard, M. Briday, J. Cohen-Tanugi, Y. Copin, H. M. Courtois, D. Fouchez, E. Gangler, D. Guinet, A. J. Hawken, Y. L. Kim, P. F. Léget, J. Neveu, P. Ntelis, P. Rosnet, and E. Nuss, Peculiar velocity cosmology with type ia supernovae (2020).
  * Palmese and Kim [2021] A. Palmese and A. G. Kim, Probing gravity and growth of structure with gravitational waves and galaxies’ peculiar velocity, Phys. Rev. D 103, 103507 (2021), arXiv:2005.04325 .
  * Ballardini  _et al._ [2022] M. Ballardini, F. Finelli, and D. Sapone, Cosmological constraints on the gravitational constant, J. Cosmology Astropart. Phys 2022, 004 (2022), arXiv:2111.09168 .
  * Komatsu  _et al._ [2011] E. Komatsu _et al._ , Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation, ApJS 192, 18 (2011), arXiv:1001.4538 .
  * Alam  _et al._ [2021] S. Alam _et al._ (eBOSS), Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory, Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 .
  * Aker  _et al._ [2022] M. Aker _et al._ (KATRIN), Direct neutrino-mass measurement with sub-electronvolt sensitivity, Nature Phys. 18, 160 (2022), arXiv:2105.08533 .
  * Motohashi  _et al._ [2013] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Cosmology Based on f(R) Gravity Admits 1 eV Sterile Neutrinos, Phys. Rev. Lett. 110, 121302 (2013), arXiv:1203.6828 [astro-ph.CO] .
  * Chudaykin  _et al._ [2015] A. S. Chudaykin, D. S. Gorbunov, A. A. Starobinsky, and R. A. Burenin, Cosmology based on f(R) gravity with O(1) eV sterile neutrino, JCAP 05, 004, arXiv:1412.5239 [astro-ph.CO] .
  * Belgacem  _et al._ [2018] E. Belgacem, Y. Dirian, S. Foffa, and M. Maggiore, Modified gravitational-wave propagation and standard sirens, Phys. Rev. D 98, 023510 (2018), arXiv:1805.08731 .
  * Kreisch and Komatsu [2018] C. Kreisch and E. Komatsu, Cosmological constraints on horndeski gravity in light of GW170817, Journal of Cosmology and Astroparticle Physics 2018 (12), 030.
  * Shafieloo  _et al._ [2012] A. Shafieloo, A. G. Kim, and E. V. Linder, Gaussian process cosmography, Phys. Rev. D 85, 123530 (2012), arXiv:1204.2272 .


# Universal contact and collective excitations of a strongly interacting Fermi
gas

Yun Li  Sandro Stringari  Dipartimento di Fisica, Università di Trento and
INO-CNR BEC Center, I-38123 Povo, Italy

###### Abstract

We study the relationship between Tan’s contact parameter and the macroscopic
dynamic properties of an ultracold trapped gas, such as the frequencies of the
collective oscillations and the propagation of sound in one-dimensional (1D)
configurations. We find that the value of the contact, extracted from the most
recent low-temperature measurements of the equation of state near unitarity,
reproduces with accuracy the experimental values of the collective frequencies
of the radial breathing mode at the lowest temperatures. The available
experiment results for the 1D sound velocities near unitarity are also
investigated.

###### pacs:

03.75.Kk, 03.75.Ss, 05.30.Fk, 67.85.-d

##  I Introduction

Significant experimental and theoretical work has been devoted in recent years
to understand the universal properties of interacting Fermi gases along the
BEC-BCS crossover (for a review, see, for example, Giorgini2008 ). More
recently, Tan has introduced a new concept for investigating universality
based on the so-called contact parameter, which relates the short-range
features of these systems to their thermodynamic properties Tan2008 ;
Werner2009 . The universality of the Tan’s relations has been proven in a
series of experiments based on the measurement of the molecular fraction
Partridge2005 , the momentum distribution, the RF spectroscopic rate at high
frequencies, the adiabatic sweep and virial theorems Stewart2010 , the spin
structure factor Kuhnle2010  and the equation of state Navon2010 . The
temperature dependence of the contact parameter has been the object of recent
theoretical Hu2011  and experimental Kuhnle2011  work.

In this paper we discuss the relationship between Tan’s contact parameter and
the frequencies of the collective oscillations of a harmonically trapped Fermi
gas near unitarity. We also investigate the relationship with the sound
velocity in highly elongated configurations. The study of the collective
oscillations along the BEC-BCS crossover in terms of the contact parameter has
been already addressed in DelloStritto2010  where upper bounds to the
collective frequencies were calculated using a sum rule approach. However, the
sum rule method developed in this work significantly overestimates the
hydrodynamic frequencies of trapped Fermi gases, being consequently
ineffective for a useful quantitative comparison with experimental data in the
superfluid hydrodynamic regime of low temperatures.

The present approach is based on a perturbative solution of the hydrodynamic
equations of superfluids near unitarity Bulgac2005 . This allows for an exact
analytic relationship between Tan’contact parameter and the deviations of the
frequencies of the collective oscillations as well as of the one-dimensional
(1D) sound velocity from their values at unitarity. The high precision
achievable in the frequency measurements is in particular expected to provide
an alternative accurate determination of the contact parameter and to further
confirm the universality of this physical quantity.

##  II Contact, equation of state and hydrodynamic equations

We start from the following definition of the contact parameter, based on the
so-called adiabatic sweep theorem Tan2008 :

$\left[\frac{dE}{d(1/a)}\right]_{N}=-\frac{\hbar^{2}\mathcal{I}}{4\pi M},$ (1)

where $E$ is the total energy of the system, $\mathcal{I}$ is the contact
parameter, $M$ is the atomic mass, and $a$ is the $s$-wave scattering length.
By using the local density approximation (LDA), the total energy can be
calculated as $E=\int d^{3}r(\epsilon+nV_{\text{ext}})$, where
$V_{\text{ext}}$ is the trapping potential and $\epsilon$ is the energy
density of a uniform gas. The equilibrium profile, in the LDA, is determined
by the equation

$\mu(n)+V_{\text{ext}}=\bar{\mu},$ (2)

where

$\mu(n)=\frac{\partial\epsilon(n,a)}{\partial n}$ (3)

is the chemical potential of uniform matter, providing the equation of state
of the gas, while $\bar{\mu}$ is the chemical potential of the trapped system,
fixed by the normalization condition. The derivative of the total energy with
respect to $1/a$ in Eq.(1) can then be conveniently written as

$\left[\frac{dE}{d(1/a)}\right]_{N}=\int
d^{3}r\,\left[\frac{\partial\epsilon(n,a)}{\partial(1/a)}\right]_{n},$ (4)

exploiting the link between the contact parameter and the equation of state.

On the other hand the chemical potential of uniform matter $\mu(n)$ is a
crucial ingredient of the hydrodynamic equations of superfluids. At zero
temperature, these equations actually read Pitaevskii2003 :

$\displaystyle\frac{\partial n}{\partial t}+\nabla\left(\mathbf{v}n\right)=0,$
(5) $\displaystyle\frac{\partial\mathbf{v}}{\partial
t}+\frac{1}{M}\nabla\left[\frac{1}{2}M\mathbf{v}^{2}+\mu(n)+V_{\text{ext}}\right]=0,$

so that an insightful question is to understand the link between the contact
parameter and the behavior of the collective modes emerging from the solutions
of Eq.(5). In the following we will discuss the problem near unitarity, where
we can expand the chemical potential in the form Tan2008 ; Bulgac2005

$\displaystyle\mu(n)=\frac{\hbar^{2}}{2M}\xi\left(3\pi^{2}n\right)^{2/3}-\frac{2\hbar^{2}}{5M}\frac{\zeta}{a}\left(3\pi^{2}n\right)^{1/3}$
(6)

with $\xi$ and $\zeta$ being two universal dimensionless parameters accounting
for the effects of the interactions. The first term in (6) exhibits the same
density dependence as the ideal Fermi gas with the renormalization factor
$\xi$. The second term (first-order correction in $1/a$) is directly related
to Tan’s parameter calculated at unitarity. Indeed, using Eqs.(1-4), one
easily finds that the contact, for a harmonically trapped system, is given by

$\frac{\mathcal{I}}{Nk_{F}}=\frac{512}{175}\zeta,$ (7)

where we have defined the Fermi wave vector $k_{F}=[3\pi^{2}n(0)]^{1/3}$
depending on the density in the center of the trap.

The small-amplitude oscillations of the gas can be studied by solving the
linearized hydrodynamic equations (5),

$-\omega^{2}\delta
n=\frac{1}{M}\nabla\cdot\left\\{n\nabla\left[\frac{\partial\mu(n)}{\partial
n}\delta n\right]\right\\},$ (8)

where $\delta n$ is the amplitude of the density oscillations around the
equilibrium value $n$. At unitarity $\mu(n)\propto n^{2/3}$, and the solutions
of Eq.(8), in the presence of harmonic trapping with axial symmetry, exhibit
the analytic form

$\omega^{2}(\lambda)=\frac{\omega_{\perp}^{2}}{3}\left[\left(4\lambda^{2}+5\right)\pm\sqrt{16\lambda^{4}-32\lambda^{2}+25}\,\right],$
(9)

with $\lambda=\omega_{z}/\omega_{\perp}$, holding for the lowest $m=0$ modes,
and include the coupling between the monopole and quadrupole oscillations
caused by the non spherical shape of the potential (here $m$ is the $z$
component of angular momentum carried by the excitation). Notice that,
remarkably, result (9) does not depend on the parameter $\xi$ characterizing
the equation of state at unitarity. This relation for the collective
frequencies was actually first obtained in contexts different from the unitary
Fermi gas, such as the ideal Bose gas above $T_{c}$ Kagan1997 ; Griffin1997
and the ideal Fermi gas Amoruso1999  in the hydrodynamic regime. In both cases
the equation of state [$\mu(n,s)$ actually exhibits the same $n^{2/3}$
dependence for fixed entropy per particle $s$] and the hydrodynamic equations
yield the same dispersion law (9) for the scaling solutions of coupled
quadrupole monopole type in the presence of harmonic trapping. The prediction
(9) for the collective frequencies was checked experimentally at unitarity
providing a direct confirmation of the universality exhibited by the unitary
Fermi gas Kinast2004 ; Altmeyer2007 . The same oscillations were investigated
out of unitarity along the BEC-BCS regime, confirming the predictions of
theory Stringari2004 ; Astrakharchik2005  and in particular the fine details
of the equation of state accounted for by quantum Monte Carlo simulations.

##  III Collective frequencies and sound velocity shifts near unitarity

In the following we calculate the deviations of the collective frequencies
from the unitary value (9) holding for small values of the dimensionless
parameter $1/k_{F}a$. To this purpose we solve the hydrodynamic equations
using a perturbative procedure, which is generally applicable to any equation
of state having the form $\mu(n)=\mu_{0}(n)+\mu_{1}(n)$, where
$\mu_{1}(n)\ll\mu_{0}(n)$ represents the first-order correction. The density
profile $n(\mathbf{r})=n_{0}(\mathbf{r})+n_{1}(\mathbf{r})$, including the
first-order correction, can be obtained starting from the expansion of
$\mu(n)$ around the zero-order ground-state density profile
$n_{0}(\mathbf{r})$ and using the equilibrium condition (2) in LDA. This gives

$n_{1}=\left[\delta\bar{\mu}-\mu_{1}(n_{0})\right]\left/\frac{\partial\mu_{0}(n_{0})}{\partial
n_{0}}\right.,$ (10)

where $\delta\bar{\mu}$ is the first-order correction to the chemical
potential $\bar{\mu}$. Solving the linearized hydrodynamic equation (8)
perturbatively, we find the following expression for the frequency shift:

$\frac{\delta\omega}{\omega}=-\frac{\int
d^{3}r\left(\nabla^{2}f_{0}^{\ast}\right)\left[n_{1}-n_{0}\left(\partial
n_{1}/\partial n_{0}\right)\right]f_{0}}{2\omega^{2}M\int
d^{3}r\,f_{0}^{\ast}\;\delta n_{0}}$ (11)

for the collective oscillations caused by the perturbation $\mu_{1}(n)$ in the
chemical potential. In Eq.(11) $f_{0}=[\partial\mu_{0}(n_{0})/\partial
n_{0}]\delta n_{0}$ is the zero-order eigenfunction of (8), and $n_{1}$ is
given by (10). According to Eq.(11), one always has $\delta\omega=0$ for the
surface modes satisfying the condition $\nabla^{2}f_{0}=0$, as expected, due
to their independence of the equation of state. For compression modes one
instead expects a correction due to the changes in the equation of state. In
general, one can show that the first term in (10) proportional to
$\delta\bar{\mu}$ gives no contribution to the frequency shift and will be
consequently neglected in the following Note_02 .

Result (11) is valid for both Bose and Fermi systems. In the case of weakly
interacting Bose-Einstein condensed gas, it allows for the calculation of the
frequency shifts caused by the Lee-Huang-Yang corrections in the equation of
state. In this case, the corresponding density correction, neglecting the term
proportional to $\delta\bar{\mu}$, can be written as
$n_{1}=-32(n_{0}a)^{3/2}/3\sqrt{\pi}$, and using Eq.(11), one finds the
frequency shift of the compression mode as derived in Pitaevskii1998 . For the
Fermi gas near unitarity, we instead employ the expansion of the equation of
state (6) around the density profile $n_{0}(\mathbf{r})$ calculated at
unitarity. Ignoring also in this case the irrelevant term proportional to
$\delta\bar{\mu}$, we find

$n(\mathbf{r})=n_{0}(\mathbf{r})+n_{1}(\mathbf{r})=n_{0}(\mathbf{r})-\frac{3\beta}{2\alpha}n_{0}^{2/3}(\mathbf{r}),$
(12)

where $\alpha=\hbar^{2}\xi(3\pi^{2})^{2/3}/2M$,
$\beta=-2\hbar^{2}\zeta(3\pi^{2})^{1/3}/(5Ma)$, and
$n_{0}=[(\bar{\mu}_{0}-V_{\text{ext}})/\alpha]^{3/2}$, with $\bar{\mu}_{0}$
being the chemical potential evaluated for a trapped system at unitarity.
After some straightforward algebra, one finds the following expression for the
frequency shift of the collective oscillations near unitarity:

$\frac{\delta\omega}{\omega}=\frac{\beta}{6\omega^{2}M}\frac{\int
d^{3}r\left(\nabla^{2}f_{0}^{\ast}\right)n_{0}^{2/3}\,f_{0}}{\int
d^{3}r\,f_{0}^{\ast}n_{0}^{1/3}\,f_{0}}.$ (13)

One can check that this result is equivalent to the result of Eq.(24) and (25)
in Bulgac2005  where a similar expansion was carried out near unitarity.

Figure 1: (Color online) Functions $\eta_{+}$ (solid blue) and $\eta_{-}$
(dashed red) relative to the higher and lower $m=0$ modes as a function of the
deformation parameter $\lambda=\omega_{z}/\omega_{\perp}$.

The eigenfunctions for the $m=0$ modes (9) have the form $f_{0}\sim
a+br_{\perp}^{2}+cz^{2}$, where

$\displaystyle\frac{a}{b}$
$\displaystyle=-\frac{\bar{\mu}_{0}[4\lambda^{2}+3(\omega/\omega_{\perp})^{2}-10]}{6M\omega_{z}^{2}},$
(14) $\displaystyle\frac{c}{b}$
$\displaystyle=\frac{3(\omega/\omega_{\perp})^{2}-10}{2}.$

After some length but straightforward algebra, one finally obtains the
following result for the frequency shift:

$\frac{\delta\omega}{\omega}=\left[\frac{128\zeta}{525\pi\xi}\,\eta_{\pm}(\lambda)\right]\left(k_{F}a\right)^{-1}=\left[\frac{\mathcal{I}/Nk^{0}_{F}}{12\pi\xi^{1/2}}\,\eta_{\pm}(\lambda)\right]\left(k^{0}_{F}a\right)^{-1},$
(15)

where

$\eta_{\pm}(\lambda)=\frac{1}{2}\pm\frac{3}{\phantom{\hat{1}}2\sqrt{16\lambda^{4}-32\lambda^{2}+25}\phantom{\hat{1}}},$
(16)

with the index $\pm$ referring to the higher $(+)$ and lower $(-)$ frequencies
of (9). In the second equality of (15), we have used relation (7) for the
contact parameter calculated for an harmonically trapped atomic cloud and we
have expressed the Fermi momentum $k_{F}$ in terms of the ideal Fermi gas wave
vector $k^{0}_{F}=(24N)^{1/6}a_{ho}^{-1}$, with $a_{ho}$ being the geometrical
average of the harmonic oscillator lengths. For the same total number of
atoms, the density $n(0)$ in the center of the trap for the unitary gas is
$\xi^{-3/4}$ times larger than for the ideal gas, yielding
$k_{F}=k^{0}_{F}\xi^{-1/4}$.

Eq.(15) represents the main result of the present paper. It relates Tan’s
contact, a central quantity for the universality relations holding in
interacting systems, with the low-energy macroscopic dynamics of the system,
namely, the frequencies of the collective oscillations. These equations can be
used either to predict theoretically the frequency shifts, once the
dimensionless parameters $\xi$ and $\zeta$ or Tan’s contact are known, or to
determine experimental constraints on the value of $\mathcal{I}$. In Fig.1 we
plot $\eta_{\pm}$ as a function of the deformation parameter $\lambda$. For a
spherical trap ($\lambda=1$) one obtains $\eta_{+}=1$ for the monopole mode
and $\eta_{-}=0$, confirming, as already anticipated, that there is no
frequency shift for the surface quadrupole mode. When $\lambda\neq 1$, the two
modes are coupled. In the limit of both spherical traps and highly elongated
traps ($\lambda\ll 1$) it reproduces the results of Bulgac2005 .

Figure 2: (Color online) Frequency of the radial compression mode for an
elongated Fermi gas in the units of radial frequency ($\omega_{\perp}$)
Altmeyer2007 . The deformation of the trap $\lambda=0.038$. The dashed red
curve refers to the equation of state based on Monte Carlo simulations. Open
and solid circles correspond to experiment measurements. The black straight
line shows the slope of the frequency shifts around unitarity obtained from
Eq.(15) using $\xi=0.41$ and $\zeta=0.93$ as extracted from the direct
measurement of the equation of state Navon2010  (corresponding to
$\mathcal{I}/Nk^{0}_{F}\simeq 3.4$). The green (light gray) straight line
instead corresponds to the values $\xi=0.41$ and $\zeta=0.74$ yielding
$\mathcal{I}/Nk^{0}_{F}\simeq 2.7$ (see text).

In Fig.2 we show the prediction of Eq.(15) in the trapping conditions of the
experiment of Altmeyer2007 , using the values $\xi=0.41$, $\zeta=0.93$
extracted from the direct measurement of the equation of state Navon2010
carried out at the lowest temperatures and yielding the value
$\mathcal{I}/Nk_{F}^{0}\simeq 3.4$ for the contact parameter. These values
differ by a few percent from the most recent theoretical predictions based on
Monte Carlo simulations at $T=0$ (see, for example, Gandolfi2011  and
references therein). The predicted slope $\delta\omega/\omega\sim
0.11/(k^{0}_{F}a)$ (black line) of the collective frequencies of the radial
breathing mode at unitarity turns out to be in very good agreement with
experiments ($\delta\omega/\omega\sim 0.12/(k_{F}^{0}a)$ Note_01 ).

In order to appreciate the sensitivity of the slope to the choice of the
values of $\xi$ and $\zeta$, in Fig.2, we also show the predictions [green
(light gray) line] for the frequency shifts using the values $\xi=0.41$ and
$\zeta=0.74$, yielding the smaller value $\mathcal{I}/Nk^{0}_{F}\simeq 2.7$
for the contact. This value is closer to the measurement of the contact
carried out in Stewart2010  and Kuhnle2011  at slightly higher values of
temperature. The resulting slope $\delta\omega/\omega\sim 0.09/(k_{F}^{0}a)$
provides a worse description of the experimental data for the collective
frequencies.

The collective oscillations discussed above represent the discretized values
of the usual sound waves described by hydrodynamics. It is actually useful to
calculate also the changes of the sound velocity of a uniform sample near
unitarity in terms of the dimensionless parameters $\xi$ and $\zeta$ or,
equivalently, Tan’s contact parameter. For bulk Fermi gases, one finds $\delta
c/c=-\zeta/(5\xi k_{F}a)$.

Figure 3: (Color online) Normalized 1D sound velocity $c_{1}/v_{F}$ vs the
interaction parameter $1/k^{0}_{F}a$ in the unitary regime, where $v_{F}=\hbar
k_{F}^{0}/M$ is the Fermi velocity for ideal Fermi gas and
$k_{F}^{0}=k_{F}\xi^{1/4}$ is the corresponding Fermi wave vector. Circles
with error bars correspond to the experiment measurements Joseph2007 . The
dashed red curve is a quadratic fitting Thomas_private , and the black
straight line is the slope of $\delta c_{1}/c_{1}$ calculated from Eq.(18)
using the values $\xi=0.41$ and $\zeta=0.93$ extracted from Navon2010 .

In the case of cylindrical geometry with radial harmonic trapping the sound
velocities can be also calculated starting from the 1D hydrodynamic expression
Capuzzi2006

$c_{1}=\left\\{\left.\frac{1}{M}\int d^{2}r_{\perp}\,n\right/\int
d^{2}r_{\perp}\,\left[\frac{\partial\mu(n)}{\partial
n}\right]^{-1}\right\\}^{1/2},$ (17)

where the the density profile along the radial direction should be evaluated
in LDA. The above results hold for sound waves characterized by wavelengths
significantly larger than the radial size of the gas. Carrying out a
perturbative development around unitarity, similar to the one employed above
for the calculation of the excitation frequencies, we obtain the result

$\frac{\delta c_{1}}{c_{1}}=-\frac{3\zeta}{20\xi}\left(k_{F}a\right)^{-1},$
(18)

where $c_{1}=(\xi/5)^{1/2}(\hbar k_{F}/M)$ is the 1D sound velocity depending
on the density in the center of the trap via $k_{F}$. Notice that $c_{1}$
differs from the sound velocity calculated at the central value of the density
by the factor $\sqrt{3/5}$. In Fig.3 we show the experimental values of
$c_{1}$ in the unitary regime measured in Joseph2007  together with the slope
evaluated from Eq.(18) using for $\xi$ and $\zeta$ the $T\simeq 0$ values
obtained from the experiment Navon2010 . A quadratic fit is applied to the
experiment data (dashed red curve) Thomas_private . The figure shows that the
value of the slope at unitarity [black curve, $\delta
c_{1}/c_{1}\sim-0.27/(k_{F}^{0}a)$] overestimates the experimental linear
changes in a visible way, suggesting that these experimental data were carried
out at relatively higher temperatures, such that they cannot be accurately
reproduced by employing the $T=0$ values of the contact parameter. Another
source of disagreement might be due to the fact that the conditions of
applicability of the 1D expression (17) for the sound velocity are not fully
satisfied in the experiment of Joseph2007 .

##  IV Conclusion

In conclusion our analysis reveals consistency between the experimental
results for the contact, obtained through the measurement of the equation of
state carried out at $T\simeq 0$ in Navon2010 , and the behavior of the
collective frequencies carried out in Altmeyer2007  at the lowest
temperatures. It would be interesting to extend our analysis of the frequency
shifts to finite temperature. The analysis could be simplified by the fact
that the scaling modes of monopole and quadrupole types have a universal
behavior at unitarity and their frequencies, calculated in the hydrodynamic
regime, are independent of temperature. Furthermore they are not coupled to
second sound. The proper calculation of the resulting slope at finite
temperature and the explicit connection with Tan’s contact parameter at finite
temperature will be the object of a future work.

###### Acknowledgements.

Useful discussions with R. Combescot and L. P. Pitaevskii are acknowledged.
This work has been supported by ERC through the QGBE grant.

## References

  * (1) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys., 80, 1215 (2008).
  * (2) S. Tan, Ann. Phys. 323, 2952 (2008); ibid 323, 2971 (2008); 323, 2987 (2008).
  * (3) E. Braaten and L. Platter, Phys. Rev. Lett. 100, 205301 (2008); R. Combescot, F. Alzetto, and X. Leyronas, Phys Rev. A 79, 053640 (2009); F. Werner, L. Tarruell, and Y. Castin, Eur. Phy. J. B 68, 401 (2009).
  * (4) G. B. Partridge, K. E. Strecker, R. I. Kamar, M. W. Jack, and R. G. Hulet, Phys. Rev. Lett., 95, 020404 (2005).
  * (5) J. T. Stewart, J. P. Gaebler, T. E. Drake, and D. S. Jin, Phys. Rev. Lett., 104, 235301 (2010).
  * (6) E. D. Kuhnle, H. Hu, X.-J. Liu, P. Dyke, M. Mark, P. D. Drummond, P. Hannaford, and C. J. Vale, Phys. Rev. Lett., 105, 070402 (2010).
  * (7) N. Navon, S. Nascimbène, F. Chevy, and C. Salomon, Science, 328, 729 (2010).
  * (8) H. Hui, X. J. Liu and P. D. Drummond, New Journal of Physics, 13, 035007 (2011).
  * (9) E. D. Kuhnle, S. Hoinka, and P. Dyke, H. Hu, P. Hannaford, and C. J. Vale, Phys. Rev. Lett., 106, 170402 (2011).
  * (10) M. DelloStritto, and T. M. De Silva, arXiv:1012.2329.
  * (11) A. Bulgac, and G. F. Bertsch, Phys. Rev. Lett., 94, 070401 (2005).
  * (12) L. P. Pitaevskii, and S. Stringari, 2003, _Bose-Einstein Condensation_ (Clarendon, Oxford).
  * (13) Y. Kagan, E. L. Surkov, and G. V. Shlyapnikov, Phys. Rev. A, 55, R18 (1997).
  * (14) A. Griffin, W. C. Wu, and S. Stringari, Phys. Rev. Lett., 78, 1838 (1997).
  * (15) M. Amoruso, I. Meccoli, A. Minguzzi, and M. P. Tosi, Eur. Phys. J. D, 7, 441 (1999).
  * (16) J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E. Thomas, Phys. Rev. Lett., 92, 150402 (2004).
  * (17) A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen, M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm, Phys. Rev. Lett., 98, 040401 (2007).
  * (18) S. Stringari, Europhys. Lett., 65, 749 (2004).
  * (19) G. E. Astrakharchik, R. Combescot, X. Leyronas, and S. Stringari, Phys. Rev. Lett., 95, 030404 (2005).
  * (20) This follows from the fact that the unperturbed frequencies [see for example Eq.(9)] do not depend on the value of the chemical potential $\bar{\mu}$.
  * (21) L. Pitaevskii, and S. Stringari, Phys. Rev. Lett., 81, 4541 (1998).
  * (22) S. Gandolfi, K. E. Schmidt, J. Carlson, Phys. Rev. A 83, 041601 (2011).
  * (23) R. Grimm (private communication).
  * (24) P. Capuzzi, P. Vignolo, F. Federici, and M. P. Tosi, Phys. Rev. A, 73, 021603 (2006).
  * (25) J. Joseph, B. Clancy, L. Luo, J. Kinast, A. Turlapov, and J. E. Thomas, Phys. Rev. Lett., 98, 170401 (2007).
  * (26) J. Joseph and J. E. Thomas (private communication).


[a]Evan Wickenden

# Approaching the Chiral and Continuum Limit of Large-N QCD

Thomas DeGrand

###### Abstract

We present preliminary results from our calculation of the low energy
constants (LECs) of the chiral effective theory for 3, 4 and 5 color QCD with
$N_{f}=2$ dynamical fermion flavors. We simulate with clover fermions over a
range of lattice couplings and quark masses. We observe the expected $N_{c}$
scaling for the LECs appropriate to the condensate $B$ and pseudoscalar decay
constant $F$. The range of quark masses over which leading order chiral
perturbation theory describes the data grows as $N_{c}$ rises.

##  1 Introduction

The large $N_{c}$ limit provides interesting insight into nonperturbative
features of QCD [1, 2, 3]. The lattice offers a means to test its predictions.
The objective of our study is to simulate across $N_{c}$ with $N_{f}=2$
flavors of dynamical fermions, measure the low energy constants (LECs) of
chiral perturbation theory ($\chi$PT), and observe their $N_{c}$ scaling. We
aim to extrapolate our results to $a\to 0$ and $N_{c}\to\infty$.

The results shown here are preliminary. We observe the expected $N_{c}$
scaling for the leading LECs: $B\sim N_{c}^{0}$ and $F\sim N_{c}^{1/2}$. A
feature we observe is that the region of parameter space where $\chi$PT
describes the data is larger at larger $N_{c}$. This is due to the fact that
the chiral expansion may be written in terms of the parameter
$\xi=\frac{m_{\pi}^{2}}{8\pi^{2}f_{\pi}^{2}}$, which, because $f_{\pi}^{2}\sim
N_{c}$, is naturally suppressed at large $N_{c}$.

A large $N_{c}$ study with $N_{f}=2$ flavors in the chiral and continuum limit
is absent from the literature. A number of other approaches to the large
$N_{c}$ limit have been studied. Hernandez et al. considered the $N_{f}=4$
case with $3-6$ colors at a single lattice spacing [4]. Bali et al. studied
spectroscopy with 2-7 and 17 colors in the quenched approximation [5]. Using
Eguchi-Kawai reduction, according to which finite-volume effects vanish as
$N_{c}\to\infty$, Gonzalez-Arroyo and collaborators have performed studies at
extremely large $N_{c}$ [6]. Our goal is to make comparisons with these and
other studies.

##  2 Wilson Chiral Perturbation Theory

$\chi$PT has a set of LECs which characterize the quark mass dependence of
quantities such as $m_{\pi}^{2}$, $f_{\pi}$. Predictions for these observables
through NNLO are present in the literature. They are typically presented using
two different expansions: the $x$ and $\xi$ expansions, where

$\displaystyle x\equiv\frac{2Bm_{q}}{4\pi^{2}F^{2}},$
$\displaystyle\xi\equiv\frac{m_{\pi}^{2}}{4\pi f_{\pi}^{2}}$ (1)

Using $\xi$ rather than $x$ amounts to a resummation of the chiral expansion.
Careful studies in the $\text{SU}(3)$ case with $N_{f}=2(+1)$ fermions
demonstrate that the $\xi$ expansion has better convergence properties, with
the NNLO formulae describing lattice data to larger fermion or pseudoscalar
masses [7, 8, 9]. A theoretical justification for this behavior is that the
NNLO chiral logs in the $\xi$ expansion have smaller numerical coefficients
than they do in the $x$ expansion. Following these results, we work primarily
with the $\xi$ expansion. We plan to do fits with the $x$ expansion as a
consistency check.

QCD on the lattice has chiral symmetry explicitly broken by both nonzero quark
mass and nonzero lattice spacing. Wilson Chiral Perturbation Theory
(W$\chi$PT) enhances chiral perturbation theory at the level of the effective
Lagrangian by introducing operators with an explicit dependence on lattice
spacing. These operators arise from the Symanzik action for QCD via a spurion
analysis. Using the power counting scheme $\delta\sim p^{2}\sim m_{q}\sim a$,
Bar, Rupak and Shoresh derive the effective chiral Lagrangian for Wilson-type
fermions through $\mathcal{O}(\delta^{2})$ [10]. The NLO W$\chi$PT predictions
for $m_{\pi}^{2}/m_{q}$ and $f_{\pi}$ are:

$\displaystyle\frac{m_{\pi}^{2}}{m_{q}}$
$\displaystyle=2B(1+\frac{1}{2}\xi\log\xi)+L_{M}\xi+W_{q}a+\tilde{W}_{qq}\frac{a^{2}}{\xi}$
$\displaystyle f_{\pi}$ $\displaystyle=F(1-\xi\log\xi)+L_{f}\xi+W_{F}a.$ (2)

The $W$’s are the correction terms due to the use of Wilson fermions.

We use the Wilson flow parameter $t_{0}$ to set the scale. $t_{0}$ depends on
the bare parameters in the simulation, and thus on measured observables. The
dependence of the gradient flow scale $t_{0}$ on $m_{PS}^{2}$ was described by
Bar and Golterman [11]. Using a $t_{0}$ which is mass-dependent (as was done
by Ref. [12]) does not alter the chiral logarithms, but it does affect the
$L_{i}$’s in Eq. 2. We need a mass-independent definition of a lattice
spacing, which we obtain by interpolating $m_{\pi}^{2}\sqrt{t_{0}}$ versus
$t_{0}$ to $0.15$ for our data at a given $\beta$.

We match scales across $N_{c}$ in the commonly - used way, taking

$t_{0}^{2}\left\langle E(t_{0})\right\rangle=C(N_{c})$ (3)

with

$C(N_{c})=C\left(\frac{3}{8}\frac{N_{c}^{2}-1}{N_{c}}\right),$ (4)

and $C=0.3$ the usual value used in $SU(3)$. Our data span a range of
$t_{0}=1$ to 3 or so, corresponding (with an $SU(3)$ value of
$\sqrt{t_{0}}\sim 0.15$ fm) to $a=0.09-0.15$ fm.

##  3 Simulation and Analysis Methodology

We use the HMC algorithm with the Wilson gauge action with nHYP links and
clover fermions [13]. We have mostly simulated on $16^{3}\times 32$ lattices,
with a few $24^{3}\times 32$ volumes used as cross checks. For each data
point, we collected between 400 and several thousand trajectories, keeping
every 10th configuration for measurements. A jackknife analysis is used to
obtain errors on and correlations between lattice observables. By increasing
the number of deletions and monitoring how the errors grow, we can roughly
estimate how autocorrelated the data are, and find them to be minimally so. We
discard at least the first hundred trajectories for thermalization, and
further test for thermalization by walking out the minimum trajectory and
watching for drifting observables. Most of the data points included in our
analysis have $m_{\pi}L>4$; a few are slightly below this threshold. The
Python package gvar is used to automatically track errors and correlations,
alongside lsqfit for correlated nonlinear least squares fitting [14, 15].

Pion masses, decay constants and Axial Ward Identity quark masses are all
determined in the usual way of fitting correlators to hyperbolic functions. As
a way to reduce systematic errors associated with choosing fit ranges, we
employ Neil and Jay’s Bayesian model-averaging method [16]. This weights each
fit proportionally to $e^{-\chi^{2}/2+N}$, where $N$ is the number of points
included in the fit, which balances a small chi-squared against a larger fit
range.

##  4 Chiral Fits

In Fig. 1 we present plots of $m_{\pi}^{2}\sqrt{t_{0}}/m_{q}$ and
$f_{\pi}\sqrt{t_{0}}$ versus $\xi$ for 3, 4, and 5 colors. We fit these to
standard NLO $\chi$PT at each lattice spacing (i.e. neglecting the Wilson
terms), and in Fig. 2 plot the $a$-dependence we obtain in $B$, $F$, $L_{M}$,
and $L_{F}$.

Figure 1: Left: pion decay constant scaled by $\sqrt{t_{0}}\sqrt{3/N_{c}}$ for
3,4,5 colors. Right: squared pion mass divided by AWI quark mass (again scaled
by $\sqrt{t_{0}}$). The lines and colored bands show the results of NLO
$\chi$PT fits to the individual $\beta$ values, neglecting the $W_{i}$’s in
Eq. 2.

Figure 2: Preliminary plots of $B$, $F$, $L_{M}$, $L_{F}$ (all in units of
$\sqrt{t_{0}}$) versus $a^{2}/t_{0}$. $F$ is scaled by $\sqrt{3/N_{c}}$. Blue
points are $\text{SU}(3)$, orange points $\text{SU}(4)$, green points
$\text{SU}(5)$.

##  5 Analysis

The plots for $F$ and $L_{F}$ show small $a$ dependence and appear to
extrapolate to a common value of $\sqrt{t_{0}}\sqrt{3/N_{c}}F\sim 0.07-0.08$
or $\sqrt{3/N_{c}}F\sim 90-100$ MeV (recall we are working in a convention
where $f_{\pi}\sim 131$ MeV), or $\sqrt{3/N_{c}}F\sim 60$\- 70 MeV in the “93
MeV” $f_{\pi}$ convention. This is a reasonable number, according to the
tables in Ref. [7]. $\sqrt{t_{0}}B$ has larger $a-$ dependence, but a value of
$\sqrt{t_{0}}B\sim 2.0-2.4$ gives $\Sigma^{1/3}\sim 250$ MeV, again not an
unreasonable value.

The above plots also convey a novel feature of this project: the chiral limit
becomes easier at larger $N_{c}$, in that a constant range of $\xi$
corresponds to a much wider range of $m_{q}$. Again, this is simply due to the
fact that $\xi\sim m_{\pi}^{2}/f_{\pi}^{2}\sim m_{q}/N_{c}$. This is
nevertheless very useful, since it means that at larger $N_{c}$, one can
obtain good chiral data while staying further away from the computational
challenges of very light dynamical quarks. $\xi$ also governs the size of
finite volume effects, so these are also reduced at larger $N_{c}$. In Fig. 2,
one can see the LECs are much better determined for $5$ colors than $4$ or $3$
(and they appear to have less lattice spacing dependence), despite arising
from fewer net trajectories. At three colors, on a $16^{3}\times 32$ lattice,
we are squeezed between finite volume effects and keeping $\xi$ small enough
for good fits to NLO chiral formulae. In fact, the bottleneck to our project
is $N_{c}=3$, which is of course much less interesting than $N_{c}>3$ due to
the availability of much better data by other collaborations. We only need it
to cross check our results. For this reason we have found it necessary to
begin running on $24^{3}\times 32$ lattices.

##  6 Conclusion

In this project we aim to make continuum predictions for the low energy
constants of chiral perturbation theory with $N_{f}=2$ dynamical fermions. We
aim to assess how well large $N_{c}$ predictions describe the physical
$\text{SU}(3)$ case and to compare our results with other approaches to the
large $N_{c}$ limit, namely studies with different fermion content—different
$N_{f}$’s and/or different representations.

We have so far only used the next-to-leading order prediction of
$\text{SU}(N_{f})$ Wilson chiral perturbation theory, with the standard power
counting scheme $\delta\sim p^{2}\sim m_{q}\sim a$. In the large $N_{c}$
limit, the axial anomaly vanishes and the $\eta$ meson recovers its status as
a Goldstone boson, so $U(N_{f})$ $\chi$PT becomes the correct low-energy
description of QCD. Considering alternative power counting schemes, NNLO
predicitons, and $U(N_{f})$ $\chi$PT are directions for future study.

Data collection and analysis continue. The goal remains first
$\text{SU}(2)\times\text{SU}(2)$ fits at each $N_{c}$, then global $U(2)\times
U(2)$ fits incorporating all data across $N_{c}$.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics under Award Number DE-
SC-0010005. Some of the computations for this work were also carried out with
resources provided by the USQCD Collaboration, which is funded by the Office
of Science of the U.S. Department of Energy using the resources of the Fermi
National Accelerator Laboratory (Fermilab), a U.S. Department of Energy,
Office of Science, HEP User Facility. Fermilab is managed by Fermi Research
Alliance, LLC (FRA), acting under Contract No. DE- AC02-07CH11359. Our
computer code is based on the publicly available package of the MILC
collaboration [17]. The version we use was originally developed by Y. Shamir
and B. Svetitsky.

## References

  * [1] G. ’t Hooft, Nucl. Phys. B 72, 461 (1974) doi:10.1016/0550-3213(74)90154-0
  * [2] G. ’t Hooft, Nucl. Phys. B 75, 461 (1974). doi:10.1016/0550-3213(74)90088-1
  * [3] E. Witten, Annals Phys. 128, 363 (1980) doi:10.1016/0003-4916(80)90325-5
  * [4] P. Hernández, C. Pena and F. Romero-López, Eur. Phys. J. C 79, no.10, 865 (2019) doi:10.1140/epjc/s10052-019-7395-y [arXiv:1907.11511 [hep-lat]].
  * [5] G. S. Bali, F. Bursa, L. Castagnini, S. Collins, L. Del Debbio, B. Lucini and M. Panero, JHEP 1306, 071 (2013) doi:10.1007/JHEP06(2013)071 [arXiv:1304.4437 [hep-lat]].
  * [6] M. G. Pérez, A. González-Arroyo and M. Okawa, JHEP 04, 230 (2021) doi:10.1007/JHEP04(2021)230 [arXiv:2011.13061 [hep-lat]].
  * [7] S. Aoki et al. [Flavour Lattice Averaging Group], Eur. Phys. J. C 80, no.2, 113 (2020) doi:10.1140/epjc/s10052-019-7354-7 [arXiv:1902.08191 [hep-lat]].
  * [8] J. Noaki et al. [JLQCD and TWQCD], Phys. Rev. Lett. 101, 202004 (2008) doi:10.1103/PhysRevLett.101.202004 [arXiv:0806.0894 [hep-lat]].
  * [9] S. Dürr et al. [Budapest-Marseille-Wuppertal], Phys. Rev. D 90, no.11, 114504 (2014) doi:10.1103/PhysRevD.90.114504 [arXiv:1310.3626 [hep-lat]].
  * [10] O. Bar, G. Rupak and N. Shoresh, Phys. Rev. D 70, 034508 (2004) doi:10.1103/PhysRevD.70.034508 [arXiv:hep-lat/0306021 [hep-lat]].
  * [11] O. Bar and M. Golterman, Phys. Rev. D 89, no.3, 034505 (2014) [erratum: Phys. Rev. D 89, no.9, 099905 (2014)] doi:10.1103/PhysRevD.89.034505 [arXiv:1312.4999 [hep-lat]].
  * [12] V. Ayyar, T. DeGrand, M. Golterman, D. C. Hackett, W. I. Jay, E. T. Neil, Y. Shamir and B. Svetitsky, Phys. Rev. D 97, no.7, 074505 (2018) doi:10.1103/PhysRevD.97.074505 [arXiv:1710.00806 [hep-lat]].
  * [13] A. Hasenfratz and F. Knechtli, Phys. Rev. D 64, 034504 (2001) doi:10.1103/PhysRevD.64.034504 [arXiv:hep-lat/0103029 [hep-lat]].
  * [14] P. Lepage, C. Gohlke, & D. Hacket, (2021). gplepage/gvar: gvar version 11.9.4 (v11.9.4). Zenodo. https://doi.org/10.5281/zenodo.5479009
  * [15] P. Lepage, C. Gohlke. (2021). gplepage/lsqfit: lsqfit version 12.0.1 (v12.0.1). Zenodo. https://doi.org/10.5281/zenodo.5512582
  * [16] W. I. Jay and E. T. Neil, Phys. Rev. D 103, 114502 (2021) doi:10.1103/PhysRevD.103.114502 [arXiv:2008.01069 [stat.ME]].
  * [17] https://github.com/milc-qcd/milc_qcd/