# Selecting energy-momentum trace dependent gravity theories with LSS

Jonas Pinheiro da Silva<EMAIL_ADDRESS>PPGCosmo, CCE, Universidade
Federal do Espírito Santo (UFES)
Av. Fernando Ferrari, 540, CEP 29.075-910, Vitória, ES, Brasil.  Hermano
Velten<EMAIL_ADDRESS>Departamento de Física, Universidade Federal
de Ouro Preto (UFOP), Campus Universitário Morro do Cruzeiro, 35.402-136, Ouro
Preto, Brazil

###### Abstract

We study scalar cosmological perturbations in $f(R,T)$ modified gravity
theories, where $T$ represents the trace of the energy-momentum tensor. We
provide detailed equations for the matter energy density contrast. We solve
them numerically to facilitate a comparison with available large-scale
structure (LSS) formation observational data on $f\sigma_{8}$, while also
addressing the $S_{8}$ tension. We identify $f(R,T)$ models that either lead
to growth enhancement or suppression. Since recent results in the literature
indicate a preference for the latter feature, this type of analysis is quite
useful for selecting viable modifications of gravity. The studied class of
such $f(R,T)$ models are either ruled out or severely restricted.

##  I Introduction

A number of considerable efforts have been devoted to assessing the validity
of $f(R,T)$ theories, where ($T\equiv g^{\mu\nu}T_{\mu\nu}$) means the trace
of the energy-momentum tensor, with special emphasis on the late time
cosmological evolution. Since in such approach there is no associated
component to the dark energy fluid, the pressureless matter, presumably a Cold
Dark Matter (CDM) component, influenced by this modified gravity theory is
responsible for driving the accelerated expansion of the universe
Velten:2017hhf ; Moraes:2019hgx ; Jeakel:2023xlp ; Myrzakulov:2023lcl ;
Iosifidis:2021kqo ; Fortunato:2023ypc ; Mishra:2024uwq ; Bhattacharjee:2019oim
; Bouali:2023fid . However, it is essential to recognize that recent research
has imposed substantial limitations on the practical usefulness of $f(R,T)$
approaches at the solar system level Bertini:2023pmp . On the other hand, in
the cosmological scenario, the results highlighted in Ref. Jeakel:2023xlp  are
of particular interest since they demonstrate that though $f(R,T)$ theories
can be consistent on late cosmological times, when a radiation component is
added to the total cosmic energy budget the ability of these theories to
maintain coherence fails.

The formulation of $f(R,T)$ theories from its Lagrangian to the field
equations level is still a matter of debate. Even the standard practise to
move the extra modified gravity terms from the left hand side to the right
hand side of the field equations forming an effective energy-momentum tensor
is disputed Fisher:2019ekh ; PhysRevD.101.108501 ; Fisher:2020zwx . Also,
issues about the definition of the coupling term between the trace of the
energy-momentum tensor and geometry are intrinsic to the formulation of
$f(R,T)$ theories. In particular, there is the claim that trace free
components do not couple to geometry as discussed in Ref. Shabani:2014xvi .
The present work is not dedicated to an exhaustive analysis of this aspect.
Instead, we direct our focus to the formulation of the dynamical equations
that dictate the behavior of cosmological matter scalar perturbations within
the context of $f(R,T)$ theories. Therefore, in our approach, by neglecting
the influence of radiation the issue about the coupling term is irrelevant. In
fact, concerning only the growth of matter density fluctuations in an
expanding universe, the cosmological matter dominated phase plays a crucial
role in the final large scale structure patterns. According to the standard
cosmological model, baryonic perturbations start to follow the already evolved
CDM fluctuations after the last scattering surface when the background
dynamics is dominated by a Einstein-de Sitter (matter like) expansion.
Therefore, in order to analyse the matter (dark plus baryonic) clustering
features in the observed LSS the radiative contribution is irrelevant for our
purposes.

We expand previous analysis performed in Ref.Alvarenga:2013syu  by considering
now the non-conservative aspect of the continuity equation, an issue avoided
in such reference. We also obtain first order dynamical equations for the
scalar perturbations which are slightly different from those obtained in
Alvarenga:2013syu . The relevance of this formulation emerges as it offers a
more complete and precise exploration of the observable manifestations that
arise during the process of structure formation in a cosmological universe
dominated by matter under the influence of energy-momentum trace dependent
modified gravity theories. Given the growth of the density contrast, an
inherent feature of a matter-dominated universe, our investigation turns to a
deep exploration of the free model parameter region where the $f(R,T)$ models
remain theoretically sound. Therefore, the main goal of this study lies in
producing results that shed light on the status of these theories concerning
the structure formation process.

In the next section II we provide a brief analytical development for obtaining
the field equations in the $f(R,T)$ family of theories. We also present the
first order expressions for the field equations which shall be used later.
Subsequently, in section III we revisit the cosmological scalar perturbations
equations in $f(R,T)$ theories. In particular, we obtain solutions for a
perfect fluid cosmological component. We then solve such equations numerically
and promote a comparison with currant available observational data. This
analysis will be very important to select viable $f(R,T)$ models. The final
section is devoted to our conclusions.

##  II General framework

The general action for the $f(R,T)$ theories is written according to the
general action Harko:2011kv ,

$S\left(g^{\mu\nu}R_{\mu\nu},g^{\mu\nu}T_{\mu\nu},\psi_{m}\right)=\int
d^{4}x\sqrt{-g}\left[L_{G}\left(R,T\right)+L_{m}(g_{\mu\nu},\psi_{m})\right],$
(1)

where the Lagrangian density for the gravitational sector,
$L_{G}=f(R,T)/\kappa^{2}$, depends on the curvature scalar ($R\equiv
g^{\mu\nu}R_{\mu\nu}$) and on the energy-momentum tensor trace $T$. The
coupling constant is defined exactly as in General Relativity (GR) i.e.,
$\kappa^{2}=8\pi G$. All matter fields are denoted by $\psi_{m}$. In turn, the
energy-momentum tensor is obtained as usual

$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta\left(L_{m}\sqrt{-g}\right)}{\delta
g^{\mu\nu}}.$ (2)

The field equations are obtained after variation of action (1). In this
process, one has to obtain the expression for $\delta f(R,T)$. In order to
focus on the contribution of the energy-momentum tensor let us now follow the
arguments presented in Ref. Shabani:2014xvi . The total energy-momentum tensor
(2) can be split into relativistic ($L^{(r)}$) and non-relativistic
($L^{(m)}$) components

$L_{m}\rightarrow L^{(m)}+L^{(r)},$ (3)

then, Eq. (2) is rewritten as

$T_{\mu\nu}=T^{(m)}_{\mu\nu}+T^{(r)}_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}L^{(m)}\right)}{\delta
g^{\mu\nu}}-\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}L^{(r)}\right)}{\delta
g^{\mu\nu}}.$ (4)

After varying the action (1) with respect to the metric, the contribution from
the energy-momentum tensor is encoded in the term,

$\displaystyle\delta f(T)=f_{T}\delta
T=f_{T}\delta\left(g^{\mu\nu}T_{\mu\nu}\right).$ (5)

Since $g^{\mu\nu}T^{(r)}_{\mu\nu}=0$, then

$\frac{\delta f(T)}{\delta
g^{\mu\nu}}=f_{T}\left(T^{(m)}_{\mu\nu}+g^{\mu\nu}\frac{\delta
T^{(m)}_{\mu\nu}}{\delta g^{\mu\nu}}\right).$ (6)

This means that only the non-relativistic matter couples to the geometric
sector in $f(R,T)$ gravity. This is a quite important issue in this sort of
theory.

Before we proceed, a cautionary remark is necessary. In a previous work by the
same authors Jeakel:2023xlp , we analysed the FLRW late-time background
expansion of $f(R,T)$ cosmologies. Regarding the late-time phenomenology, we
found that $f(R,T)$ cosmologies have a narrow free model parameter space in
agreement with current observations. However, in such reference, when we
attempted to analyse the complete cosmological scenario, including the early-
time aspects, the inclusion of radiation has been made by adding the
relativistic pressure at the field equations level i.e., we replaced the total
pressure $p$ by the sum of non-relativistic and relativistic components
$p_{m}+p_{r}$. Consequently, we found the inviability of $f(R,T)$ cosmologies
when considering the radiation fluid, as this component would no longer scale
according to $\rho_{r}\propto a^{-4}$ as usual, due to the non-conservation of
the energy-momentum tensor in $f(R,T)$. In summary, the conclusions presented
in Jeakel:2023xlp  are irrelevant for the purposes of the present work, as we
will not take into account relativistic components into account here.

The action (1) leads to the following field equations111The coupling between
matter and curvature implies in the existence of a matter tensor
$\Theta_{\mu\nu}$ (see Harko:2011kv  for details). In this case, we write
$\Theta_{\mu\nu}^{m}=g^{\mu\nu}\frac{\delta T^{m}_{\mu\nu}}{\delta
g^{\mu\nu}}$. Due to the perfect fluid structure, one has two possibilities
for the matter Lagrangian: $L_{m}=\rho$ or $L_{m}=-p$. A detailed discussion
about the choice of the Lagrangian and its implications on the dynamics is
found in Refs.Sotiriou:2008it ; Bertolami:2008ab ; Faraoni:2009rk .

$f_{R}R^{\mu}_{\nu}-\frac{f(R,T)}{2}\delta^{\mu}_{\nu}-\Delta^{\mu}_{\nu}f_{R}=\kappa^{2}T^{\mu}_{\nu}+f_{T}T^{\mu(m)}_{\nu}.$
(7)

The term $\Delta^{\mu}_{\nu}$, is defined as a combination of operators

$\Delta^{\mu}_{\nu}\equiv
g^{\mu\alpha}\nabla_{\alpha}\nabla_{\nu}-\delta^{\mu}_{\nu}g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}.$
(8)

The terms proportional to $f_{T}$ in equation (7) are the geometry-matter
coupling terms that arises from the variation of the trace of the energy-
momentum tensor. It is worth noting that for a pressureless fluid ($p_{m}=0$),
mimicking a matter dominated universe, the right hand side of (7) is written
simple as $(\kappa^{2}+f_{T})T^{\mu(m)}_{\nu}$ i.e., the combination
$(\kappa^{2}+f_{T})$ becomes the effective coupling between the energy-
momentum tensor and geometry.

A complete description of $f(R,T)$ theories in the metric and metric-affine
formalism can be found, respectively, in Refs.Harko:2011kv ;
Barrientos:2018cnx .

An important feature of the $f(R,T)$ class of theories is its non-conservative
nature. Generally, the presence of matter and (or) coupling with second order
terms in the $R$-curvature in the Lagrangian density of the gravitational
sector implies in such non-conservation. Via the Bianchi identity
$\nabla_{\mu}G^{\mu}_{\nu}=0$ (the restricted case where $f(R,T)$ depends
linearly on $R$ plus a trace dependence) and propriety’s from Levi-Civita
connection, $\nabla_{\mu}g^{\mu\nu}=0$, one finds that equation (7) turns into

$\nabla_{\mu}T^{\mu(m)}_{\nu}=-\frac{1}{(\kappa^{2}+f_{T})}\left(\nabla_{\mu}f_{T}T^{\mu(m)}_{\nu}+f_{T}\frac{1}{2}\nabla_{\mu}T\right).$
(9)

Once more, conservative $f(R,T)$ models can be conceived by forcing the right
hand side of the above expression to vanish. For pressureless fluids this
happens when $f(R,T)=R+T^{1/2}$. This is the case studied in Ref.
Alvarenga:2013syu .

Hereafter, we will avoid using the index $m$ for non-relativistic matter.
Then, all equations below apply to a universe fully dominated by non-
relativistic matter.

###  II.1 Perturbed field equations

Since we will be interested in the perturbative dynamics of $f(R,T)$ theories,
let us present the first order version of the field equation (7). It reads

$\displaystyle\delta G^{\mu}_{\nu}=\kappa^{2}\delta
T^{\mu}_{\nu}+f_{\bar{T}\bar{T}}\delta
T\left(\bar{T}^{\mu}_{\nu}+\bar{p}\delta^{\mu}_{\nu}\right)+f_{\bar{T}}\left(\delta
T^{\mu}_{\nu}+\delta p\delta^{\mu}_{\nu}+\frac{1}{2}\delta
T\delta^{\mu}_{\nu}\right).$ (10)

From now on all non-perturbed (background) terms carry an overbar. In order to
arrive to the above equation we have assumed the minimally coupled scenario
$f(R,T)=R+f(T)$. In this case $f_{R}=\partial f(R,T)/\partial R=1$.
Correspondingly, the derivative $f_{T}=\partial f(R,T)/\partial T$ is written
at first order as $\delta f_{T}=f_{\bar{T}\bar{T}}\delta T$.

From (9), the general form of the first order continuity equation is written
as

$\displaystyle\left(\partial_{\mu}\delta
T^{\mu}_{\nu}+\delta\Gamma^{\mu}_{\mu\lambda}\bar{T}^{\lambda}_{\nu}+\bar{\Gamma}^{\mu}_{\mu\lambda}\delta
T^{\lambda}_{\nu}-\delta\Gamma^{\lambda}_{\mu\nu}\bar{T}^{\mu}_{\lambda}-\bar{\Gamma}^{\lambda}_{\mu\nu}\delta
T^{\mu}_{\lambda}\right)(\kappa^{2}+f_{\bar{T}})=-\biggl{[}f_{\bar{T}\bar{T}}\nabla_{\mu}\bar{T}^{\mu}_{\nu}\delta
T+f_{\bar{T}\bar{T}}\partial_{\mu}\bar{T}\left(\delta T^{\mu}_{\nu}+\delta
p\delta^{\mu}_{\nu}\right)$
$\displaystyle+\left(f_{\bar{T}\bar{T}\bar{T}}\partial_{\mu}\bar{T}\delta
T+f_{\bar{T}\bar{T}}\partial_{\mu}\delta
T\right)\left(\bar{T}^{\mu}_{\nu}+\bar{p}\delta^{\mu}_{\nu}\right)+f_{\bar{T}\bar{T}}\delta
T\left(\frac{1}{2}\partial_{\mu}\bar{T}+\partial_{\mu}\bar{p}\delta^{\mu}_{\nu}\right)+f_{\bar{T}}\left(\frac{1}{2}\partial_{\mu}\delta
T+\partial_{\mu}\delta p\delta^{\mu}_{\nu}\right)\biggr{]}.$ (11)

This will be used in the next section.

##  III Scalar perturbations for a perfect fluid

Cosmological perturbation theory is designed to indicating spacetime
anisotropies and inhomogeneities that are not considered in the Friedmann-
Lemaître-Robertson-Walker (FLRW) metric. Thus, we can define the FLRW (flat)
background metric, $\bar{g}_{\mu\nu}$, as

$\bar{g}_{\mu\nu}\equiv g_{\mu\nu}-\delta g_{\mu\nu},$ (12)

where in this case, $g_{\mu\nu}$ denotes the metric of the FLRW-perturbed
spacetime with the contribution of the perturbations $\delta g_{\mu\nu}$. In
terms of the conformal time $\tau$, the FLRW metric reads

$d\bar{s}^{2}=a(\tau)^{2}\left(d\tau^{2}-\delta_{ij}dx^{i}dx^{j}\right),$ (13)

Via (12) this leads us to the FLRW-perturbed metric

$ds^{2}=a(\tau)^{2}\eta_{\mu\nu}+\delta g_{\mu\nu},$ (14)

where $\eta_{\mu\nu}$ is the Minkowski metric. As we are only interested in
the first order scalar perturbation structure, let us fix the gauge

$\delta g_{\mu\nu}=2a(\tau)^{2}\left(\begin{array}[]{cc}\phi&0\\\
0&\psi\delta_{ij}\\\ \end{array}\right),$ (15)

hence, the line element (14) becomes

$ds^{2}=a(\tau)^{2}\left[(1+2\phi)d\tau^{2}-(1-2\psi)\delta_{ij}dx^{i}dx^{j}\right].$
(16)

With the metric signature used here we have the same line element used as used
e.g., in Mukhanov:1988jd ; Alvarenga:2013syu . In literature the line element
(16) is called of conformal Newtonian gauge or longitudinal gauge. The
components of perturbed energy-momentum tensor for cosmological perfect fluids
in this gauge is given by

$\delta T^{0}_{0}=\delta\rho;\hskip 14.22636pt\delta T^{i}_{j}=-\delta
p\delta^{i}_{j};\hskip 14.22636pt\text{and}\hskip 14.22636pt\delta
T^{0}_{i}=-\delta T^{i}_{0}=-(1+\bar{\omega})\bar{\rho}v_{i},$ (17)

being $\bar{\rho}$ and $\bar{p}$ the background contribution of the fluid’s
energy density and pressure, respectively, connected by the equation of state
parameter $\bar{\omega}=\bar{p}/\bar{\rho}$. The four velocity is defined as
$v_{\mu}=(0,v_{i})$. Developing the right side of the (10) one finds

$\displaystyle\delta G^{\mu}_{\nu}=\left(\kappa^{2}+f_{\bar{T}}\right)\delta
T^{\mu}_{\nu}+\left[f_{\bar{T}\bar{T}}(1-3c^{2}_{s})(1+\bar{\omega})\bar{\rho}v^{\mu}v_{\nu}+\frac{f_{\bar{T}}}{2}(1-c^{2}_{s})\delta^{\mu}_{\nu}\right]\delta\rho.$
(18)

For the sake of generality, let us consider an equation of state
$\bar{\omega}\equiv\bar{\omega}(\bar{\rho})$. Therefore, an adiabatic process
is defined by the relation

$\delta
p=\left(\bar{\omega}+\bar{\rho}\frac{d\bar{\omega}}{d\bar{\rho}}\right)\delta\rho=c^{2}_{s}\delta\rho,$
(19)

where $c^{2}_{s}$ is the square sound speed. If one adopts a constant equation
of state parameter $\bar{\omega}$, then $\bar{\omega}=c^{2}_{s}$. After such
definitions and trivial calculations, the $\\{00\\}$, $\\{0i\\}$ and
$\\{ij\\}$ components of (18) are, respectively

$\nabla^{2}\psi-3\mathcal{H}\left(\psi^{\prime}+\mathcal{H}\phi\right)=\frac{a^{2}(\eta)}{2}\kappa^{2}\bar{\rho}\left\\{1+\frac{1}{\kappa^{2}}\left[f_{\bar{T}\bar{T}}(1-3c^{2}_{s})(1+c^{2}_{s})\bar{\rho}+\frac{f_{\bar{T}}}{2}(3-c^{2}_{s})\right]\right\\}\delta,$
(20)
$(\phi^{\prime}+\mathcal{H}\phi)_{,i}=\frac{a^{2}(\eta)}{2}\kappa^{2}\bar{\rho}(1+\bar{\omega})\left(1+\frac{f_{\bar{T}}}{\kappa^{2}}\right)v_{i},$
(21)
$\phi^{\prime\prime}+3\mathcal{H}\phi^{\prime}+(2\mathcal{H}^{\prime}+\mathcal{H}^{2})\phi=\frac{a^{2}(\eta)}{2}\kappa^{2}\bar{\rho}\left[c^{2}_{s}-\frac{f_{\bar{T}}}{2\kappa^{2}}(1-3c^{2}_{s})\right]\delta,\hskip
14.22636pt\text{(if $i=j$)},$ (22)

and

$\left(\phi-\psi\right)_{,ij}=0\hskip 14.22636pt\text{(if $i\neq j$)};$ (23)

being the matter density contrast defined as
$\delta\equiv\delta\rho/\bar{\rho}$. Note that for perfect (shear free) fluids
we only have one degree of freedom associated with the scalar perturbed metric
potentials (the result $\phi=\psi$ in (23) arises naturally).

Elaborating further the continuity equation (II.1) we find for the
$\nabla_{\mu}T_{0}^{\mu}$ component, in Fourier space
($\nabla^{2}\rightarrow-k^{2}$, being $k$ the wavenumber),

$\displaystyle(\delta\rho)^{\prime}+3\mathcal{H}(1+c_{s}^{2})\delta\rho-(1+\bar{\omega})(k^{2}v+3\phi^{\prime})\bar{\rho}=\frac{(1-3c_{s}^{2})^{2}\bar{\rho}^{\prime}\bar{\rho}\delta}{(\kappa^{2}+f_{\bar{T}})^{2}}\left[(1+\bar{\omega})f_{\bar{T}\bar{T}}\bar{\rho}+\frac{f_{\bar{T}}}{2}+\frac{f_{\bar{T}}c_{s}^{2}}{1-3c_{s}^{2}}\right]$
$\displaystyle-\frac{(1-3c_{s}^{2})\bar{\rho}\delta^{\prime}}{\kappa^{2}+f_{\bar{T}}}\biggl{[}f_{\bar{T}\bar{T}}\bar{\rho}+\frac{1}{2}f_{\bar{T}}+\bar{\omega}f_{\bar{T}\bar{T}}\bar{\rho}+\frac{c_{s}^{2}f_{\bar{T}}}{1-3c_{s}^{2}}\biggr{]}-\frac{(1-3c_{s}^{2})^{2}\bar{\rho}^{\prime}\delta}{\kappa^{2}+f_{\bar{T}}}\biggl{\\{}f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}^{2}(1+\bar{\omega})+\frac{1}{2}f_{\bar{T}\bar{T}}\bar{\rho}+$
$\displaystyle\frac{1}{1-3c_{s}^{2}}\biggl{[}2f_{\bar{T}\bar{T}}\bar{\rho}(1+\bar{\omega})+\frac{1}{2}f_{\bar{T}}+c_{s}^{2}f_{\bar{T}\bar{T}}\bar{\rho}+\frac{c_{s}^{2}f_{\bar{T}}}{1-3c_{s}^{2}}\biggr{]}\biggr{\\}}.$
(24)

At the same time, the momentum component $\nabla_{\mu}T^{\mu}_{i}$ reads,

$\displaystyle-(1+\bar{\omega})\left[(\bar{\rho}v)^{\prime}+\bar{\rho}\phi+4\bar{\rho}\mathcal{H}v\right]=-(1+\bar{\omega})(1-3c_{s}^{2})f_{\bar{T}\bar{T}}\bar{\rho}\bar{\rho}^{\prime}v+\frac{(1-3c_{s}^{2})}{2}\\{[(1-3c_{s}^{2})f_{\bar{T}\bar{T}}\bar{\rho}+f_{\bar{T}}]\bar{\rho}^{\prime}\delta+f_{\bar{T}}\bar{\rho}\delta^{\prime}\\}+$
$\displaystyle+\bar{\omega}(1-3c_{s}^{2})[f_{\bar{T}\bar{T}}\delta\rho\bar{\rho}^{\prime}+f_{\bar{T}\bar{T}\bar{T}}\delta\rho\bar{\rho}^{\prime}+f_{\bar{T}\bar{T}}\bar{\rho}(\delta\rho)^{\prime}].$
(25)

Now, for the sake of simplicity, let us define the following combination

$\bar{\xi}=1+\frac{1}{\kappa^{2}+f_{\bar{T}}}\left(f_{\bar{T}\bar{T}}\bar{\rho}+\frac{1}{2}f_{\bar{T}}\right).$
(26)

By eliminating the $T$ dependence on $f(R,T)$ one gets $\bar{\xi}=1$. This
happens in either RG or $f(R)$ theories. This definition will appear below in
several equations. The relevance of this definition is understood when
analysing the conservation equation (9). In fact, for the background FLRW
metric and considering a pressureless fluid i.e., $\bar{w}=c_{s}^{2}=0$, the
continuity equation (9) is written as

$\bar{\xi}\bar{\rho}^{\prime}+3\mathcal{H}\bar{\rho}=0.$ (27)

Back to the perturbed level, adopting a pressureless fluid and eliminating the
$\bar{\rho}^{\prime}$ terms with help of (27), the Eqs. (III) and (III) turn
into, respectively

$\bar{\xi}\delta^{\prime}-(k^{2}v+3\phi^{\prime})=-\frac{3\mathcal{H}\bar{\rho}\delta}{\kappa^{2}+f_{\bar{T}}}\left[f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\left(f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\right)\right],$
(28)

and

$v^{\prime}+\phi+\mathcal{H}v\left[4-\frac{3}{\bar{\xi}}\left(1+\frac{f_{\bar{T}\bar{T}}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\right)\right]=\frac{1}{2(\kappa^{2}+f_{\bar{T}})}\left[f_{\bar{T}}\delta^{\prime}-\frac{3\mathcal{H}\delta}{\bar{\xi}}\left(f_{\bar{T}\bar{T}}\bar{\rho}+f_{\bar{T}}\right)\right].$
(29)

By combining the time derivative of Eq. (28) with (29) we obtain a general
equation for the evolution of the matter energy density contrast $\delta$. It
reads

$\displaystyle\biggl{\\{}\delta^{\prime\prime}+\delta^{\prime}\mathcal{H}\biggl{[}4-\frac{3}{\bar{\xi}}\biggl{(}1+\frac{f_{\bar{T}\bar{T}}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\biggr{)}\biggr{]}\biggr{\\}}\bar{\xi}+\frac{3}{\kappa^{2}+f_{\bar{T}}}\biggl{[}f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\left(f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\right)\biggr{]}\biggl{[}\mathcal{H}^{\prime}\bar{\rho}\delta+\delta^{\prime}\mathcal{H}\bar{\rho}+$
$\displaystyle+\frac{3\mathcal{H}^{2}\bar{\rho}}{\bar{\xi}}\delta\biggl{(}1-\frac{f_{\bar{T}\bar{T}}}{\kappa^{2}+f_{\bar{T}}}\biggr{)}\biggr{]}+3\mathcal{H}\biggl{[}4-\frac{3}{\bar{\xi}}\biggl{(}1+\frac{f_{\bar{T}\bar{T}}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\biggr{)}\biggr{]}\biggl{\\{}\frac{\mathcal{H}\bar{\rho}\delta}{\kappa^{2}+f_{T}}\biggl{[}f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\left(f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\right)\biggr{]}-\phi^{\prime}\biggr{\\}}-$
$\displaystyle-3\phi^{\prime\prime}+\delta^{\prime}\xi^{\prime}+\frac{3\mathcal{H}\bar{\rho}\delta}{\kappa^{2}+f_{\bar{T}}}\biggl{[}-\frac{3\mathcal{H}\bar{\rho}}{\bar{\xi}}f_{\bar{T}\bar{T}\bar{T}}+\frac{\bar{\xi}^{\prime}}{\bar{\xi^{2}}}\biggl{(}\frac{5}{2}f_{\bar{T}\bar{T}}+f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}\biggr{)}+\frac{3\mathcal{H}\bar{\rho}}{\bar{\xi}^{2}}\biggl{(}f_{\bar{T}\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{7}{2}f_{\bar{T}\bar{T}\bar{T}}\biggr{)}\biggr{]}+$
$\displaystyle+k^{2}\biggl{\\{}\phi-\frac{1}{2(\kappa^{2}+f_{\bar{T}})}\biggl{[}f_{\bar{T}}\delta^{\prime}-\frac{3\mathcal{H}\delta}{\bar{\xi}}(f_{\bar{T}\bar{T}}\bar{\rho}+f_{\bar{T}})\biggr{]}\biggr{\\}}=0.$
(30)

General relativity is recovered by making $f(T)=0$ and the $\Lambda$CDM limit
is achieved with a specific choice of free model parameters.

Equation (III) is the general equation for the matter energy density contrast
in $f(R,T)$ models. In order to solve it and, since we will consider scales
that are within the sound horizon, we apply to it the quasi-static
approximation adopting that time derivatives of $\phi$ are subdominant with
respect to spatial derivatives and we will also rewrite $\phi$ in (III) in
terms of the remaining scalar sector quantities using (20) and (21) (see
Maggiore:2018sht ; Brandenberger:2003vk  for a review). This leaves us finally
with a more compact format

$\displaystyle\delta^{\prime\prime}+\mathcal{F}\mathcal{H}\delta^{\prime}-\frac{3}{2}\mathcal{H}^{2}\mathcal{G}\delta=0,$
(31)

where we have defined the effective Hubble friction term $\mathcal{F}$ and the
effective gravitational coupling $\mathcal{G}$, respectively:

$\displaystyle\mathcal{F}\equiv\frac{1}{\bar{\xi}}\biggl{\\{}\biggl{[}4-\frac{3}{\bar{\xi}}\left(1+\frac{f_{\bar{T}\bar{T}}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\right)\biggr{]}\bar{\xi}+\frac{3}{\kappa^{2}+f_{\bar{T}}}\biggl{[}f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\left(f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\right)\biggr{]}\bar{\rho}+\frac{\bar{\xi}^{\prime}}{\mathcal{H}}-\frac{k^{2}f_{\bar{T}}}{2\mathcal{H}(\kappa^{2}+f_{\bar{T}})}\biggr{\\}};$
(32)

and

$\displaystyle\mathcal{G}\equiv-\frac{2}{3\mathcal{H}^{2}\bar{\xi}}\biggl{\\{}\frac{3}{\kappa^{2}+f_{\bar{T}}}\biggl{[}f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\biggl{(}f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\biggr{)}\biggr{]}\biggl{[}\mathcal{H}^{\prime}\bar{\rho}+\frac{3\mathcal{H}^{2}\bar{\rho}}{\bar{\xi}}\biggl{(}1-\frac{f_{\bar{T}\bar{T}}}{\kappa^{2}+f_{\bar{T}}}\biggr{)}\biggr{]}+\frac{3\mathcal{H}^{2}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\biggl{[}4-$
$\displaystyle-\frac{3}{\bar{\xi}}\left(1+\frac{f_{\bar{T}\bar{T}}\bar{\rho}}{\kappa^{2}+f_{\bar{T}}}\right)\biggr{]}\biggl{[}f_{\bar{T}\bar{T}}-\frac{1}{\bar{\xi}}\biggl{(}f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{5}{2}f_{\bar{T}\bar{T}}\biggr{)}\biggr{]}+\frac{3\mathcal{H}\bar{\rho}}{\kappa^{2}+f_{T}}\biggl{[}-\frac{3\mathcal{H}\bar{\rho}}{\bar{\xi}}f_{\bar{T}\bar{T}\bar{T}}+\frac{\bar{\xi}^{\prime}}{\bar{\xi}^{2}}\biggl{(}\frac{5}{2}f_{\bar{T}\bar{T}}+f_{\bar{T}\bar{T}\bar{T}}\bar{\rho}\biggr{)}+$
$\displaystyle+\frac{3\mathcal{H}\bar{\rho}}{\bar{\xi}^{2}}\biggl{(}f_{\bar{T}\bar{T}\bar{T}\bar{T}}\bar{\rho}+\frac{7}{2}f_{\bar{T}\bar{T}\bar{T}}\biggr{)}\biggr{]}-\frac{3k^{2}\mathcal{H}}{2(\kappa^{2}+f_{\bar{T}})\bar{\xi}}\left(f_{\bar{T}\bar{T}}\bar{\rho}+f_{\bar{T}}\right)-\frac{3}{2}\mathcal{H}^{2}\left[1+\frac{1}{\kappa^{2}}\left(f_{\bar{T}\bar{T}}\bar{\rho}+\frac{3}{2}f_{\bar{T}}\right)\right]\biggr{\\}}.$
(33)

In (31), by recovering GR with $f(R,T)=R$, then $\mathcal{F}=1$ and
$\mathcal{G}=1$ as expected. Our task now is to assess how the chosen $f(R,T)$
models impact the evolution of the linear scalar perturbations.

##  IV Numerical results and data analysis

We proceed now solving Eq. (31) for competitive $f(R,T)$ models. We will test
two different approaches i.e., the exponential one

$f(R,T)=R+\alpha e^{\beta T},$ (34)

and the polynomial one

$f(R,T)=R+\gamma_{n}T^{n}.$ (35)

These models cover all existing proposed $f(R,T)$ models in literature.
Hereafter, we will study them separately i.e., we either apply the exponential
one or the polynomial one. We treat them as separate cases. Each model has the
GR limit when either $\alpha$ or $\gamma_{n}$ vanish, respectively. In the
exponential model, if $\alpha$ has a non-vanishing value and the trace
dependence disappears with $\beta=0$, then the $\Lambda$CDM model is achieved.
The later limit is also obtained in the polynomial model with $\gamma_{n}\neq
0$ and $n=0$.

For each model there are two free parameters. One more than the flat
$\Lambda$CDM model. Both parameters will be redefined as dimensionless
quantities (denoted by an overbar) according to the definitions

$\bar{\alpha}=\frac{\alpha}{\kappa^{2}\rho_{0}};\hskip
28.45274pt\bar{\beta}=\beta\rho_{0};\hskip 28.45274pt\text{and}\hskip
28.45274pt\bar{\gamma}_{n}=\frac{\gamma_{n}\rho_{0}^{n-1}}{\kappa^{2}}.$ (36)

But what are the preferred free model parameters according to current
cosmological data? This issue has been addressed recently in Ref.
Jeakel:2023xlp . This reference has used the gas mass fraction in galaxy
clusters, the background cosmic chronometers data and the age of the universe
to test these two models. This analysis has shown that the late time dynamics
can be fairly described by $f(R,T)$ gravity but only for a narrow range of the
free model parameters. For both models the accepted parameter region is the
one in which the trace dependence is quite weak ($\bar{\beta}$ and $n$ values
close to zero), making these models very close to the $\Lambda$CDM cosmology.
We will explore latter $\bar{\beta}$ and $n$ values within the upper limit
established by Ref. Jeakel:2023xlp . However, in the exponential model the
accepted $\bar{\alpha}$ and $\bar{\beta}$ values extend further than the
region close to the $\Lambda$CDM model. For example, combinations like
$\\{\bar{\alpha},\bar{\beta}\\}=\\{1.2,0.3\\},\\{1.1,0.45\\}$ and
$\\{1,0.6\\}$ are allowed by available background observational data. We ask
the reader to check Figs. $1$ and $2$ of Ref. Jeakel:2023xlp  for such
results.

Given the expressions presented in (32) and (III) the analytical treatment of
(31) is not viable. We have then to solve this equation numerically for both
models (34)-(35).

We adopt $k=0.1hMpc^{-1}$ for the wavenumber appearing in (32) since it
corresponds scales inside the horizon until today. In all of our results the
initial conditions are taken at the matter-radiation equality and they are set
with the help of the CAMB code Lewis:1999bs . Specifically, we fix initial
conditions for the density contrast of order $\delta(z_{eq})\sim 10^{-4}$
where $z_{eq}=3400$ refers to the matter-radiation equality in the concordance
$\Lambda$CDM cosmology.

After calculating the evolution of the matter density contrast $\delta$, we
can compare the perturbative dynamics of $f(R,T)$ models with recent available
data from LSS surveys. In this task, the relevant quantity is the combination
$f\sigma_{8}$ where $f$ is the growth rate given by

$f\equiv\frac{d\ln{\delta}}{d\ln{a}},$ (37)

while $\sigma_{8}^{2}$ is the variance of the perturbed cosmic density field
smoothed within a sphere of radius $8h^{-1}Mpc$ calculated in the linear
regime as

$\sigma_{8}(z)=\frac{\delta(z)}{\delta(z=0)}\sigma_{8}(z=0).$ (38)

As a function of the redshift this combination reads

$f\sigma_{8}(z)=-(1+z)\frac{\sigma_{8}(z=0)}{\delta(z=0)}\frac{d}{dz}\delta(z).$
(39)

One can also calculate the so called $S_{8}$ parameter defined according to

$S_{8}=\sigma_{8}(z)\left(\frac{\Omega_{m}}{0.3}\right)^{1/2}.$ (40)

Current attempts to measure this parameter have revealed an apparent
cosmological tension. Inferred $S_{8}$ values from the Cosmic Microwave
Background (CMB) data obtained by the PLANCK mission are in disagreement with
estimations from galaxy surveys like eBOSS (Extended Baryon Oscillation
Spectroscopic Survey) eBOSS:2020yzd , DES (Dark Energy Survey) DES:2017myr  or
KiDS (Kilo-Degree Survey) Heymans:2020gsg . These surveys provide a lower
value for $S_{8}$ with relevant statistical confidence level than the one
inferred by CMB data. Since $S_{8}$ quantifies the amplitude of matter
fluctuations, there is the possibility to have the following interpretation:
galactic surveys are indicating a suppressed clustering growth. This same
trend has also been recently pointed out in Ref. Nguyen:2023fip  via the
direct fitting of available $f\sigma_{8}$ data. These results will be used
below the select $f(R,T)$ models that have such clustering growth suppression.

We show in Fig. (1) the evolution of the $f\sigma_{8}$ combination as a
function of the redshift in the studied $f(R,T)$ models. The red curve
corresponds to the flat $\Lambda$CDM model with $\Omega_{\Lambda}=0.7$. The
curves for the exponential $f(R,T)$ model are displayed in the black lines.
The polynomial one in green. For both models we fix the parameters
$\bar{\alpha}=\bar{\gamma}_{n}=1.4$ and vary either $\bar{\beta}$ (the
exponential model) or $n$ (the polynomial model). The values used for
$\bar{\beta}$ and $n$ are shown in the figure. Available observational data
shown in gray is obtained from Refs. Benisty:2020kdt ; Said:2020epb ;
Kazantzidis:2018rnb ; Sakr:2023bms . The inset of this figure amplifies the
low redshift region. The qualitative behavior seen here has a general trend.
The polynomial model (green lines) remains quite similar to the $\Lambda$CDM
model in the past but it has an amplified growth at low redshift. This feature
is better seen in the inset of Fig. (1). The exponential models present the
oscillatory feature at redshifts of order $10$ to $20$, depending on the
$\bar{\beta}$ values. The solid black line, corresponding to
$\bar{\beta}=0.001$, is identical to the red one ($\Lambda$CDM) at low
redshfits. However, even with such small value $\bar{\beta}=0.001$ there is an
oscillatory feature at $z\sim 20$.

Figure 1: Evolution of $f\sigma_{8}$ as a function of the redshift. The data
points are taken from Refs. Benisty:2020kdt ; Said:2020epb ;
Kazantzidis:2018rnb ; Sakr:2023bms . The black (green) lines refer to
different $\bar{\beta}(n)$ values for the exponential (polynomial) model
fixing $\bar{\alpha}(\bar{\gamma}_{n})=1.4$. The inset amplifies the low
redshift region.

For a better visualization of how the different modified gravity model
parameters influence the evolution of the density contrast given by Eq. (31)
let us plot the functions $\bar{\xi}$ (26), $\mathcal{F}$ (32) and
$\mathcal{G}$ (III). They are displayed respectively in Fig. (2), in the top
panel of (3) and in the botton panel of (3).

Figure 2: Evolution of the parameter $\bar{\xi}$ as a function of the
redshift. Black lines represent the exponential model where we have fixed
$\bar{\alpha}=1.4$. For the polynomial models (green lines) we have fixed
$\bar{\gamma}_{n}=1.4$. The red line is flat $\Lambda$CDM limit with
$\Omega_{\Lambda}=0.7$.

The behavior of $\bar{\xi}$ seen in (2) influences the conservation law (27).
For the different $n$ values used in the polynomial model, the resulting
curves are almost indistinguishable from GR. Again, the adopted values for the
parameter $n$ are based on the results presented in Ref. Jeakel:2023xlp  i.e.,
only $n\lesssim 0.1$ are observationally allowed. Conversely, the exponential
model has a quite different evolution. In the exponential model, by keeping
$\bar{\alpha}=1.4$, only $\bar{\beta}\lesssim 0.1$ are allowed from the
background analysis. With such range of values for $\bar{\alpha}$ and
$\bar{\beta}$, the resulting behavior of $\bar{\xi}$ is shown in the black
lines of (2). The quantity $\bar{\xi}$ increases for high redshifts. This
makes the non-conservation aspect more relevant during the early stage of the
universe. Apart from the numerical visualization provided by this figure, it
is also possible to have an analytical hint about this behavior. By
substituting the $f(R,T)$ models (34) and (35) into (26) one finds that in the
high density (past) limit the $\bar{\xi}$ quantity tends to a constant value
in the polynomial model. This constant value differs from the unity
proportionally to $n$. However, for the exponential model, $\bar{\xi}$ is
proportional to the matter density parameter $\Omega$ which is a monotonically
increasing function towards high redshifts. This explains the behavior
observed in Fig. (2).

The difference between the polynomial and the exponential models is also clear
when analysing the curves for $\mathcal{F}$ and $\mathcal{G}$ in Fig. (3). The
curves for the exponential models bounces at high redshift values. The smaller
the $\bar{\beta}$ value, later this bounce happens and the amplitude of
oscillation becomes smaller. We have checked this bounce is related to the
change in the sign of the derivative of the $\bar{\xi}$ function as shown in
Fig. (2). A visual inspection of both figures reveals this feature.

Figure 3: Evolution of the term effective Hubble friction term $\mathcal{F}$
and the funtion $\mathcal{G}$ as a function of the redshift.

We have so far investigated the free model parameters in which the studied
$f(R,T)$ models are very close to the $\Lambda$CDM at background level.
However, concerning the exponential model, even such small $\bar{\beta}$
values are able to leave pathological imprints like oscillations on the
$f\sigma_{8}$ observable. We remind now that according to Ref. Jeakel:2023xlp
the observationally allowed free space parameters in the exponential model can
be extended to regions beyond the $\Lambda$CDM limit. As mentioned before, the
combinations like
$\\{\bar{\alpha},\bar{\beta}\\}=\\{1.2,0.3\\},\\{1.1,0.45\\}$ and
$\\{1,0.6\\}$ are in agreement with background cosmological data. We test
therefore the perturbative dynamics with this set of free model parameter
values in the exponential model. We show in Fig. (4) the $f\sigma_{8}$ curves
for these parameter values (black lines) in comparison to the flat
$\Lambda$CDM model with $\Omega_{\Lambda}=0.7$ (the red curve) obtained in the
limit $\bar{\beta}=0$ and $\bar{\alpha}=1.4$. Both black curves have a
slightly larger $f\sigma_{8}$ amplitude at high redshifts but they transit to
a suppressed amplitude at small redshifts, in the region where observational
data is available.

Figure 4: Evolution of $f\sigma_{8}$ as a function of the redshift. The data
points are taken from Refs. Benisty:2020kdt ; Said:2020epb ;
Kazantzidis:2018rnb ; Sakr:2023bms . The black lines refer to different
combinations of $\bar{\alpha}$ and $\bar{\beta}$ parameters for the
exponential model. The inset amplifies the low redshift region.

Apart from the analysis using the $f\sigma_{8}$ observable, let us also use
equations for the gravitational potential (22), derived in the previous
section, in order to compute the Integrated Sachs-Wolfe (ISW) effect of both
$f(R,T)$ models and compare it with the prediction of a flat $\Lambda$CDM
model with $\Omega_{m0}=0.3$. The $H_{0}$ parameter plays no role in this
analysis. We calculate the relative amplifications ($Q$) of the ISW effect as

$Q=\frac{\left(\frac{\Delta T}{T}\right)^{f(R,T)}_{ISW}}{\left(\frac{\Delta
T}{T}\right)^{\Lambda CDM}_{ISW}}-1.$ (41)

If $Q>0$ $(<0)$ the $f(R,T)$ models produce more (less) temperature variation
to the CMB photons via the ISW effect than the adopted $\Lambda$CDM model.
Such an analysis of the ISW effect has been proposed in Dent:2008ek  and used
in Velten:2011bg ; Velten:2012uv . Due to our previous results, we present
only the results for the exponential model. With the same parameter values as
used in (4) the $Q$ values for
$\\{\bar{\alpha},\bar{\beta}\\}=\\{1.2,0.3\\},\\{1.1,0.45\\}$ and
$\\{1,0.6\\}$ are, respectively, $32.9\%,29.9\%,$ and $24.0\%$. This means
that all tested values enhance the ISW signal in comparison to the
$\Lambda$CDM model. Such level of enhancement can not be probed by current
cross-correlation studies aiming to detect the ISW effect (see, for example,
Bahr-Kalus:2022yny ).

##  V Conclusions

We have studied $f(R,T)$ modified gravity models, being $T$ the trace of the
energy-momentum tensor. We have focused on the analysis of cosmological linear
scalar perturbations. The equation for the evolution of the matter density
contrast $\delta$ has been presented in the compact form (31). We remark our
equation is different from the one previously found in Ref. Alvarenga:2013syu
even adopting the same approximation, the quasi-static one. We have then
studied the perturbative behavior of the exponential $f(R,T)$ model (34) and
the polynomial one (35).

In order to assess the viability of such models we have used the $f\sigma_{8}$
observable. The polynomial model presents an excess of clustering amplitude at
low redshifts leading to high $f\sigma_{8}$ values (green lines in Fig. (1)).
This behavior is challenged by the recent findings of Ref. Nguyen:2023fip
which indicates that a suppressed clustering is preferred by data. Therefore,
we can conclude the polynomial model is not a competitive $f(R,T)$ model in
light of $f\sigma_{8}$ data.

Concerning the exponential model, our first analysis focused on the free
parameter space near the flat $\Lambda$CDM model i.e., the limit $\beta\sim 0$
and $\bar{\alpha}=1.4$. The latter value corresponds to $\Omega_{\Lambda}=0.7$
when $\bar{\beta}$ vanishes. With the parameter values near the $\Lambda$CDM
limit we found oscillations in the $f\sigma_{8}$ evolution (black lines of
Fig. (1)). This pathological behavior also compromises the viability of such
parameter values. We have then moved to another region of the free model
parameters, by testing the set
$\\{\bar{\alpha},\bar{\beta}\\}=\\{1.2,0.3\\},\\{1.1,0.45\\}$ and
$\\{1,0.6\\}$, also allowed at background level as shown in Ref.
Jeakel:2023xlp . The corresponding $f\sigma_{8}$ curves for this set of
parameter values, as shown in Fig. (4), are in good agreement with
observational data. The calculated $S_{8}$ values are, respectively,
$0.81,0.75$ and $0.75$, alleviating the $S_{8}$ tension.

In conclusion, we have extended in this work the analysis performed in Ref.
Jeakel:2023xlp  to the scalar linear perturbations. Now, we are finally able
to exclude the polynomial $f(R,T)$ model as a competitive candidate to explain
the late time accelerated cosmological expansion. On the other hand, the
exponential model has passed this test, but its allowed free model parameter
space has been severely constrained.

## Acknowledgments

Partial support from Brazilian funding agencies PROPPi/UFOP, CAPES, CNPq,
FAPEMIG and FAPES is acknowledged. We thank Sunny Vagnozzi for useful
comments.

## References

  * (1) H. Velten and T. R. P. Caramês, Phys. Rev. D 95 (2017) no.12, 123536, doi:10.1103/PhysRevD.95.123536, [arXiv:1702.07710 [gr-qc]].
  * (2) P. H. R. S. Moraes, P. K. Sahoo and S. K. J. Pacif, Gen. Rel. Grav. 52 (2020) no.4, 32, doi:10.1007/s10714-020-02681-3, [arXiv:1905.00417 [gr-qc]].
  * (3) A. P. Jeakel, J. Pinheiro da Silva and H. Velten, Phys. Dark Univ. 43, 101401 (2024), doi:10.1016/j.dark.2023.101401, [arXiv:2303.15208v3 [astro-ph.CO]]
  * (4) N. Myrzakulov, M. Koussour, A. H. A. Alfedeel and E. I. Hassan, Chin. Phys. C 47, no.11, 115107 (2023), doi:10.1088/1674-1137/acf2fa, [arXiv:2308.09913 [astro-ph.CO]].
  * (5) D. Iosifidis, N. Myrzakulov and R. Myrzakulov, Universe 7 (2021) no.8, 262 doi:10.3390/universe7080262, [arXiv:2106.05083 [gr-qc]].
  * (6) J. A. S. Fortunato, P. H. R. S. Moraes, J. G. d. Júnior and E. Brito, Eur. Phys. J. C 84, no.2, 198 (2024), doi:10.1140/epjc/s10052-024-12544-9, [arXiv:2305.01325 [gr-qc]].
  * (7) R. K. Mishra and N. Jain, Int. J. Theor. Phys. 63, no.1, 29 (2024), doi:10.1007/s10773-024-05553-7
  * (8) S. Bhattacharjee and P. K. Sahoo, Grav. Cosmol. 26, no.3, 281-284 (2020), doi:10.1134/S0202289320030032, [arXiv:1908.06759 [gr-qc]].
  * (9) A. Bouali, H. Chaudhary, T. Harko, F. S. N. Lobo, T. Ouali and M. A. S. Pinto, Mon. Not. Roy. Astron. Soc. 526, 4192-4208 (2023), doi:10.1093/mnras/stad2998, [arXiv:2309.15497 [gr-qc]].
  * (10) N. R. Bertini and H. Velten, Phys. Rev. D 107 (2023) no.12, 124005, doi:10.1103/PhysRevD.107.124005, [arXiv:2303.09699 [gr-qc]].
  * (11) S. B. Fisher and E. D. Carlson, Phys. Rev. D 100, no.6, 064059 (2019), doi:10.1103/PhysRevD.100.064059, [arXiv:1908.05306 [gr-qc]].
  * (12) Harko, T. & Moraes, P. Comment on “Reexamining f(R,T) gravity”. Phys. Rev. D. 101, 108501 (2020,5), doi:10.1103/PhysRevD.101.108501.
  * (13) S. B. Fisher and E. D. Carlson, Phys. Rev. D 101, no.10, 108502 (2020), doi:10.1103/PhysRevD.101.108502, [arXiv:2004.02934 [gr-qc]].
  * (14) H. Shabani and M. Farhoudi, Phys. Rev. D 90, no.4, 044031 (2014), doi:10.1103/PhysRevD.90.044031, [arXiv:1407.6187 [gr-qc]].
  * (15) F. G. Alvarenga, A. de la Cruz-Dombriz, M. J. S. Houndjo, M. E. Rodrigues and D. Sáez-Gómez, Phys. Rev. D 87 (2013) no.10, 103526 [erratum: Phys. Rev. D 87 (2013) no.12, 129905], doi:10.1103/PhysRevD.87.103526, [arXiv:1302.1866 [gr-qc]].
  * (16) T. Harko, F. S. N. Lobo, S. Nojiri and S. D. Odintsov, Phys. Rev. D 84 (2011), 024020, doi:10.1103/PhysRevD.84.024020, [arXiv:1104.2669 [gr-qc]].
  * (17) E. Barrientos, F. S. N. Lobo, S. Mendoza, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 97 (2018) no.10, 104041, doi:10.1103/PhysRevD.97.104041, [arXiv:1803.05525 [gr-qc]].
  * (18) T. P. Sotiriou and V. Faraoni, Class. Quant. Grav. 25, 205002 (2008), doi:10.1088/0264-9381/25/20/205002, [arXiv:0805.1249 [gr-qc]].
  * (19) O. Bertolami, F. S. N. Lobo and J. Paramos, Phys. Rev. D 78, 064036 (2008), doi:10.1103/PhysRevD.78.064036, [arXiv:0806.4434 [gr-qc]].
  * (20) V. Faraoni, Phys. Rev. D 80, 124040 (2009), doi:10.1103/PhysRevD.80.124040, [arXiv:0912.1249 [astro-ph.GA]].
  * (21) V. F. Mukhanov, Sov. Phys. JETP 67 (1988), 1297-1302.
  * (22) R. H. Brandenberger, Lect. Notes Phys. 646 (2004), 127-167, doi:10.1007/978-3-540-40918-2_5, [arXiv:hep-th/0306071 [hep-th]].
  * (23) M. Maggiore, Oxford University Press, 2018, ISBN 978-0-19-857089-9
  * (24) A. Lewis, A. Challinor and A. Lasenby, Astrophys. J. 538, 473-476 (2000), doi:10.1086/309179, [arXiv:astro-ph/9911177 [astro-ph]].
  * (25) S. Alam et al. [eBOSS], Phys. Rev. D 103 (2021) no.8, 083533, doi:10.1103/PhysRevD.103.083533, [arXiv:2007.08991 [astro-ph.CO]].
  * (26) T. M. C. Abbott et al. [DES], Phys. Rev. D 98 (2018) no.4, 043526, doi:10.1103/PhysRevD.98.043526, [arXiv:1708.01530 [astro-ph.CO]].
  * (27) C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hildebrandt, B. Joachimi, K. Kuijken, C. A. Lin, A. G. Sánchez and J. L. van den Busch, et al. Astron. Astrophys. 646 (2021), A140, doi:10.1051/0004-6361/202039063, [arXiv:2007.15632 [astro-ph.CO]].
  * (28) N. M. Nguyen, D. Huterer and Y. Wen, Phys. Rev. Lett. 131 (2023) no.11, 111001, doi:10.1103/PhysRevLett.131.111001, [arXiv:2302.01331 [astro-ph.CO]].
  * (29) D. Benisty, Phys. Dark Univ. 31 (2021), 100766, doi:10.1016/j.dark.2020.100766, [arXiv:2005.03751 [astro-ph.CO]].
  * (30) K. Said, M. Colless, C. Magoulas, J. R. Lucey and M. J. Hudson, Mon. Not. Roy. Astron. Soc. 497 (2020) no.1, 1275-1293, doi:10.1093/mnras/staa2032, [arXiv:2007.04993 [astro-ph.CO]].
  * (31) L. Kazantzidis and L. Perivolaropoulos, Phys. Rev. D 97 (2018) no.10, 103503, doi:10.1103/PhysRevD.97.103503, [arXiv:1803.01337 [astro-ph.CO]].
  * (32) Z. Sakr, Universe 9 (2023) no.8, 366, doi:10.3390/universe9080366, [arXiv:2305.02863 [astro-ph.CO]].
  * (33) J. B. Dent, S. Dutta and T. J. Weiler, Phys. Rev. D 79 (2009), 023502, doi:10.1103/PhysRevD.79.023502, [arXiv:0806.3760 [astro-ph]].
  * (34) H. Velten and D. J. Schwarz, JCAP 09 (2011), 016, doi:10.1088/1475-7516/2011/09/016, [arXiv:1107.1143 [astro-ph.CO]].
  * (35) H. Velten and D. Schwarz, Phys. Rev. D 86 (2012), 083501, doi:10.1103/PhysRevD.86.083501, [arXiv:1206.0986 [astro-ph.CO]].
  * (36) B. Bahr-Kalus, D. Parkinson, J. Asorey, S. Camera, C. Hale and F. Qin, Mon. Not. Roy. Astron. Soc. 517 (2022) no.3, 3785-3803, doi:10.1093/mnras/stac2040, [arXiv:2204.13436 [astro-ph.CO]] .


# First-order coherence of light emission from inhomogeneously broadened
mesoscopic ensembles

A. Delteil, V. Blondot, S. Buil, J.-P. Hermier  Université Paris-Saclay, UVSQ,
CNRS, GEMaC, 78000, Versailles, France
<EMAIL_ADDRESS>
###### Abstract

Inhomogeneous broadening is well known to hinder individual characteristics of
emitters, supplanting the single-particle properties by their broader
probability distribution. Here, we present an analysis of the emission spectra
of mesoscopic ensembles of inhomogeneously distributed emitters below the
thermodynamic limit (101 – 104 emitters). Based on a simple analytical model
and an extensive numerical analysis, we show that the number and individual
linewidths of the emitters can be directly estimated from the ensemble
autocorrelation function in spite of an inhomogeneously broadened emission. As
an application, we analyze the photoluminescence of colloidal nanocrystal
aggregates embedded in a gold shell. Our general method can be applied to a
wide range of mesoscopic many-body systems and could provide new insights into
their first-order coherence properties.

##  I I. Introduction

Spectroscopy of emitter ensembles distinguishes between two main regimes,
depending on the way the environment affects the individual components. In the
homogeneous broadening regime, the spectral shape of the ensemble emission is
identical to the single-particle spectrum. The field autocorrelation function
then decays as the inverse of the individual linewidth $\delta E$. On the
other hand, in the inhomogeneously distributed regime, the spectral shape is
solely given by the probability distribution of the center frequencies. This
is typically the case for Doppler-broadened atomic gases Patel63 , ensembles
of size-distributed semiconductor quantum dots Gammon96 ; Empedocles99AM ,
color centers with varying mechanical or electromagnetic environments
Rosenzweig18 ; Lindner18 , excitons in quantum wells Houdre96 , spin-dephased
Raman transitions Hemmer01 ; Sun16 , and so on. In such cases, the
autocorrelation function is therefore given by the ensemble linewidth and
rapidly decays as its inverse $T_{2}^{*}$ Sun16 ; Chen13 . In some cases,
techniques such as spectral hole burning Palinginis03  or dynamical decoupling
Abella66 ; Press10  can provide access to the underlying homogeneous dephasing
time.

Here, we investigate mesoscopic ensembles of emitters that are inhomogeneously
distributed, yet are insufficiently numerous to lead to a fully
inhomogeneously dephased autocorrelation function. We show analytically and
numerically that, in contrast to the thermodynamic limit where the coherence
is fully lost in a typical time $T_{2}^{*}\sim\hbar/\Delta E_{\mathrm{inh}}$,
in the case of mesoscopic ensembles the loss is only partial and some
information can be retrieved from the self-normalized first-order
autocorrelation function. The reminiscent exponential tail indeed retains
information on both the number and intrinsic linewidth of the individual
emitters that compose the ensemble. Using numerical simulations, we propose a
simple method to experimentally access these quantities. We emphasize that
there is no need for prerequisite information about the number of individual
emitters, their intensity, their individual properties or their inhomogeneous
distribution.

As an illustration, we then apply our method to experimental data from
mesoscopic ($N\sim 10^{1}$–$10^{4}$) ensembles of nanocrystals (NCs) of
inhomogeneous linewidth $\sim 40$ meV, exceeding by about one order of
magnitude the individual linewidths, to infer an estimation of the individual
emitter properties. We compare our results to geometrical estimations of the
NC population of 22 NC aggregates of varying sizes as well as to the emission
brightness. This analysis allows us to deduce the fraction of photoactive
emitters. Finally, we implement a stochastic optimization algorithm to show
that the emission spectra are well reproduced using the extracted parameters.

##  II II. Model and simulations

###  II.1 A. Framework and model

For the sake of simplicity, we consider ensembles of $N$ individual emitters
of identical Lorentzian emission spectrum (at the exception of their center
emission wavelength), and we term their individual coherence time $T_{2}$ –
even though part of the individual linewidth could contain time-averaged
broadening such as spectral diffusion. We denote $E_{i}=E_{0}+\Delta E_{i}$
($i=1,...,N$) their center energies, with $E_{0}$ the center energy of the
ensemble and $\Delta E_{i}$ the individual deviations from $E_{0}$. $\Delta
E_{i}$ are therefore random variables that are identically distributed
according to the thermodynamic limit inhomogeneous spectrum of typical width
$\hbar/T_{2}^{*}$. The $N$-particle power spectrum reads:

$S(E)\propto\sum_{i=1}^{N}\mathcal{L}(E-E_{i})=\sum_{i=1}^{N}\mathcal{L}(E)*\delta(E-E_{i})$
(1)

where $\mathcal{L}(E)=\dfrac{1}{(\hbar/T_{2})^{2}+(E-E_{0})^{2}}$ is the
Lorentzian distribution. The first-order autocorrelation function
$g^{(1)}(\tau)$ can be obtained by Fourier transform of the power spectrum:

$g^{(1)}(\tau)=Ce^{-\frac{|\tau|}{T_{2}}}e^{i\frac{E_{0}}{\hbar}\tau}\sum_{i=1}^{N}e^{i\frac{\Delta
E_{i}}{\hbar}\tau}$

where $C$ is a normalization constant. Let us focus on its modulus
$\left|g^{(1)}(\tau)\right|=Ce^{-\frac{|\tau|}{T_{2}}}M_{N}$, where we
introduce the random variable $M_{N}=\left|\sum_{i=1}^{N}e^{i\frac{\Delta
E_{i}}{\hbar}\tau}\right|$. We consider two opposite limits: when $\tau=0$, we
simply obtain $\left|g^{(1)}(\tau)\right|=CN$ owing to the vanishing phase
argument in $M_{N}$. On the other hand, when $\tau\gg T_{2}^{*}$, the standard
deviation of the random variable $X=\frac{\Delta E_{i}}{\hbar}\tau$ is much
greater than 1 and therefore, thanks to the $2\pi$ periodicity of the complex
exponential, $X$ can be approximated by a uniformly distributed variable on
the $[0,2\pi)$ interval. The calculation of $M_{N}$ is then analogous to a
two-dimensional random walk in the complex plane with fixed step size but
random direction, and was treated by Lord Rayleigh in 1880 Rayleigh80 . It can
be easily derived using the central limit theorem and its probability
distribution function when $N\gg 1$ is a Rayleigh distribution
$P(x)=2x/N\times e^{-x^{2}/N}$ with expectation value
$E[M_{N}]=\sqrt{N\pi}/2$. This leads us to the following expression for the
expectation value of the self-normalized first-order correlation function:

$E\left[\left|\dfrac{g^{(1)}(\tau)}{g^{(1)}(0)}\right|\right]=\dfrac{\sqrt{\pi}\exp\left(-\frac{|\tau|}{T_{2}}\right)}{2\sqrt{N}}$
(2)

This result shows that the first-order autocorrelation function contains
information about the individuals (namely $N$ and $T_{2}$) in its long-time
tail. Interestingly, Eq. 2 is independent of the distribution function of the
emitters center energies. However, this information is contained in the
intermediate timescales $0<\tau\lesssim T_{2}^{*}$ where the autocorrelation
quickly converges towards the Fourier transform of the emitter spectral
distribution. In addition, it is noteworthy that Eq. 2 does not involve any
additional normalization or calibration procedure based for instance on the
intensity of the emitters, since the autocorrelation is simply normalized by
its value at zero delay.

In practice, in a given physical system, $N$ and $\Delta E_{i}$ are fixed.
Therefore, it is not possible to repeat a random draw to estimate the
expectation value of Eq. 2. However, in the next section we will show that a
fitting procedure of the autocorrelation can yield the desired estimation
since, when $\tau$ increases, the phase argument in the exponential terms of
$M_{N}$ continuously varies, thereby randomly probing the whole support of
$X$.

###  II.2 B. Numerical simulations

We generate example spectra of varying emitter numbers based on Eq. 1. With no
loss of generality, we choose parameters that are close to those of the
emitters experimentally studied in Sec. III: $E_{0}=1.96$ eV, $\Delta
E_{\mathrm{inh}}=40$ meV. We use a Gaussian distribution of the center
energies $\Delta E_{i}$ Kagan96 ; Lipta07 . Figure 1 displays six examples
with $N$ varying between 100 and 1000. Fluctuations of the spectrum envelope
can be observed, with randomly varying widths and amplitudes – preventing to
directly extract individual properties with a fitting procedure. In the
following, we show that a fit of the autocorrelation function allows to access
characteristics of the individual spectra.

Figure 1: Six examples of simulated spectra with the number of emitters
varying from $N=100$ to $N=1000$.

We use fast Fourier transform (FFT) to calculate the autocorrelation functions
from the spectra. Figure 2 shows typical examples of
$A(\tau)=|{g^{(1)}(\tau)}/{g^{(1)}(0)}|$ calculated from ensembles of various
$N$ and $\delta E$. Around $\tau=0$, $A(\tau)$ exhibits a Gaussian shape of
width $\sim 1/\Delta E_{\mathrm{inh}}$ that originates from the inhomogeneous
distribution, and that is identical for all the ensembles, independently of
$N$ and $\delta E$. At longer $\tau$, a fluctuating exponential tail is
observed. Its amplitude depends on $N$ [fig. 2(a)] and its decay time depends
on $T_{2}$ [Fig. 2(b)]. As plotted in orange solid lines on fig. 2(a) and (b),
these tails are well fitted by exponential functions $Ae^{-\tau/\tau_{0}}$
using log-linear fitting. The least mean square fitting procedure ensures that
the fluctuations are averaged out if the fitting interval is longer than the
typical timescale of the fluctuations (which is of order $T_{2}^{*}$), such
that the fitting parameters constitute estimators of $N$ and $T_{2}$. The
theoretical curves calculated using Eq. 2 are plotted on the same graphs
(brown curves), showing the excellent agreement between the theoretical
expectation value and the fit results.

Figure 2: (a) Autocorrelation of ensembles with varying size. (b)
Autocorrelation of ensembles with varying individual linewidths. The orange
lines are fits to the data for $\tau>T_{2}^{*}$, and the brown lines are
calculated with Eq. 2.

We repeat this fitting procedure on a large number of numerically generated
signals, with $N$ varying between $10$ and $10^{4}$, and $\delta E$ varying
between 1 and 10 meV. The fitting results are shown on Fig. 3. The extracted
parameters follow the expected dependences, with moderate uncertainty ($\Delta
A/A\approx 15$ % for the amplitude, and $\Delta\tau_{0}/\tau_{0}\approx 10$ %
for the decay time), showing that our method allows to extract $N$ and $T_{2}$
from the first-order autocorrelation function calculated from the emission
spectrum.

In the next section, we illustrate an implementation of our method on
experimental data.

Figure 3: Fitting parameters extracted from $\sim 3\times 10^{4}$ numerically
generated datasets. (a) Relative amplitude, plotted as a function of $N$, in
log-log scale. Inset: same data, plotted in linear scale as a function of
$\sqrt{N}$. (b) Decay time, plotted as a function of $\delta E$, in log-log
scale. Inset: Same data, plotted in linear scale as a function of $1/\delta
E$. The error bars represent the standard deviation of the extracted
parameters.

##  III III. Experimental data

In this section, we focus on the luminescence spectra of mesoscopic aggregates
of CdSe/CdS/CdZnS NCs Blondot20  embedded within a silica and a gold shells,
termed golden supraparticles (GSPs). The synthesis of these emitters is
detailed in Blondot22 . The gold shell enhances the light emission through the
Purcell effect, thus ensuring that the role of Förster resonant energy
transfer can be neglected Blondot20 ; Blondot22 . The GSPs considered in this
section comprise between $10^{2}$ and $10^{4}$ NCs. Their photoluminescence
spectra were measured in a confocal microscope at 4 K under non-resonant laser
excitation. Figure 4 shows six representative spectra from GSPs of various
sizes. As with the simulated spectra of Fig. 1, the envelopes exhibit shape
irregularities that originate from the random center wavelength of a finite
number of emitters. The overall distribution is not completely Gaussian but
slightly bimodal, with a main subpopulation centered around 630 nm, and a
secondary redder subpopulation centered around 640 nm.. As we evidenced in
Sec. II, this has no impact on the analysis, which is independent of the
lineshape of the inhomogeneous distribution.

Figure 4: Six typical photoluminescence spectra of individual GSPs of various
sizes.

In the same way as in Sec. II, we calculate the autocorrelation using FFT. The
results are shown on Fig. 5 (blue curves) for the six spectra from Fig. 4. The
autocorrelation functions exhibit a long-time tail that can be well fitted by
an exponential function with a constant offset using log-linear least mean
square optimization (orange curves in Fig. 5).

Figure 5: Zoom-in of the autocorrelation function (blue curves) and fit
results (orange curve) for the same emitters as in Fig. 4. The gray curves
represent the parts of the data that are not fitted.

We applied the same procedure to 22 GSPs. Figure 6 shows the results of the
estimation of $N$ and $\delta E$ on these emitters. The error bars are
slightly higher than those shown in Fig. 3 due to the shorter time range on
which the exponential tail is fitted. Concerning $\delta E$, we obtain values
in the 2 to 5 meV range, in agreement with previous studies on individual NCs
at low temperature Empedocles99 ; Coolen08 . These values also consistently
match the narrowest features observed on the spectra. Concerning $N$, to gain
more insight, we now compare our conclusions to the number $N_{\mathrm{geom}}$
estimated from the measured volume of the GSPs.

Figure 6: $N$ and $\delta E$ as deduced from the autocorrelation fit for the
22 GSPs. The error bars are calculated from the uncertainties on the fitting
parameters.

To estimate the number of NCs in the GSPs, we measured the diameters of the 22
GSPs using a scanning electron microscope (SEM) and an atomic force microscope
(AFM). Knowing the center-to-center distance of individual NCs (8.5 nm), the
number of NCs in a GSP can be calculated by assuming a compact random packing
of the NCs Zaccone22 . Figure 7(a) shows a comparison of $N$ and
$N_{\mathrm{geom}}$. While these two numbers are sizably different for some
GSPs, we note that we always obtain $N\lesssim N_{\mathrm{geom}}$. This
observation shows that some GSPs contain a large fraction of optically
inactive NCs, which could be due to either damage during the synthesis or to
photobleaching of part of the NC population of these GSPs. The fraction
$F=N/N_{\mathrm{geom}}$ of optically active NCs is plotted in Fig. 7(b) as a
function of $N_{\mathrm{geom}}$. To further confirm that the observed
discrepancy between $N$ and $N_{\mathrm{geom}}$ is due to optically inactive
NCs, we plot the count rate measured at a fixed laser power as a function of
$N$ [Fig. 7(c)] and of $N_{\mathrm{geom}}$ [figure 7(d)], together with a
linear fit. The count rate exhibits a higher correlation with $N$ than
$N_{\mathrm{geom}}$, with a twice higher covariance, thus providing further
evidence of our assumption.

Figure 7: (a) Estimated number $N$ plotted against the geometrically estimated
number $N_{\mathrm{geom}}$. (b) Fraction of optically active NCs plotted
against $N_{\mathrm{geom}}$ (c) Count rate as a function of $N$. (d) Count
rate as a function of $N_{\mathrm{geom}}$.

To provide additional evidence of the relevance of the extracted parameters,
we implemented a fitting algorithm where $N$ and $\delta E$ are fixed, given
by the values extracted from the autocorrelation function, while the
individual center wavelengths constitute the parameter vector to optimize.
This algorithm is based on simulated annealing, a variant of random wall
climbing where the parameter vector to optimize is allowed to scatter
backwards, with a probability that decreases as a function of a decreasing
parameter $T$. The results of this optimization are shown in Fig. 8. The
excellent agreement between simulated and experimental spectra indicates that
the spectra can be perfectly reproduced using the parameters extracted from
our autocorrelation analysis. We note that, although the number of free
parameters is large (the center frequencies of all the emitters), a good
agreement cannot be obtained if $N$ or $\delta E$ are not well chosen.

Figure 8: Blue curves: Spectra from Fig. 4. Orange curves: Simulated annealing
optimized spectra.

##  IV IV. Conclusion

We analyzed spectral properties of mesoscopic ensembles in a regime that lies
halfway between homogeneous ($T_{2}$-limited) and ensemble-averaged
($T_{2}^{*}$-limited) regimes. In this mesoscopic regime, the spectral
distribution is inhomogeneously broadened, yet retains information about
individual components. This information can be retrieved from the self-
normalized first-order correlation obtained from the emission spectrum by
Fourier transform, as pointed out by our simple analytical model. We performed
a numerical analysis showing that an exponential fitting of the long-time
components of the self-normalized autocorrelation allows to reliably access
the number and individual linewidths of the constituents. We applied this
method to experimental data from colloidal II-VI NC aggregates, allowing to
estimate the number and linewidths of the optically active NCs. The extracted
parameters provide a better agreement with the count rate and optical spectra
than simple counting from geometrical parameters. Our method is very general
and could be applied to a wide range of mesoscopic systems, for which it could
constitute a useful characterization tool. We therefore expect our work to
shine a new light on the coherence properties of composite physical systems in
the intermediate dephasing regime.

##  V Acknowledgments

The authors thank J.-J. Greffet and C. Arnold for fruitful discussions. The
authors also acknowledge A. Bogicevic, T. Pons and N. Lequeux for GSP
synthesis. This work is supported by the French Agence Nationale de la
Recherche (ANR), under grant ANR-17-CE24-0046 (project GYN).

##  VI Data availability

The data generated in this study are available at
https://doi.org/10.5281/zenodo.7034845.

## References

  * (1) C. K. N. Patel, Determination of Atomic Temperature and Doppler Broadening in a Gaseous Discharge with Population Inversion, Phys. Rev. 131, 1582 (1963).
  * (2) D. Gammon, E. S. Snow, B. V. Shanabrook, D. S. Katzer, and D. Park, Fine Structure Splitting in the Optical Spectra of Single GaAs Quantum Dots, Phys. Rev. Lett. 76, 3005 (1996).
  * (3) S. A. Empedocles, R. Neuhauser, K. Shimizu, and M. G. Bawendi, Photoluminescence from Single Semiconductor Nanostructures, Advanced Materials 11, 1243 (1999).
  * (4) Y. Rosenzweig, Y. Schlussel, and R. Folman, Probing the origins of inhomogeneous broadening in nitrogen-vacancy centers with Doppler-free-type spectroscopy, Phys. Rev. B 98, 014112 (2018).
  * (5) S. Lindner, A. Bommer, A. Muzha, A. Krueger, L. Gines, S. Manda, O. Williams, E. Londero, A. Gali, and C. Becher, Strongly inhomogeneous distribution of spectral properties of silicon-vacancy color centers in nanodiamonds, New J. Phys. 20 115002 (2018).
  * (6) R. Houdré, R. P. Stanley, and M. Ilegems, Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: resolution of a homogeneous linewidth in an inhomogeneously broadened system. Phys. Rev. A 53, 2711 (1996).
  * (7) P. R. Hemmer, A. V. Turukhin, M. S. Shahriar, and J. A. Musser, Raman-excited spin coherences in nitrogen-vacancy color centers in diamond, Optics Letters 26, 361 (2001).
  * (8) Z. Sun, A. Delteil, S. Faelt and A. Imamoğlu, Measurement of spin coherence using Raman scattering, Phys. Rev. B 93, 241302 (2016).
  * (9) O. Chen et al., Compact high-quality CdSe–CdS core–shell nanocrystals with narrow emission linewidths and suppressed blinking, Nature Mater. 12, 445 (2013).
  * (10) P. Palinginis, S. Tavenner, M. Lonergan and H. Wang, Spectral hole burning and zero phonon linewidth in semiconductor nanocrystals, Phys. Rev. B 67, 201307 (2003).
  * (11) I. D. Abella, N. A. Kurnit, and S. R. Hartmann, Photon Echoes, Phys. Rev. 141, 391 (1966).
  * (12) D. Press, K. De Greve, P. L. McMahon, T. D. Ladd, B. Friess, C. Schneider, M. Kamp, S. Höfling, A. Forchel, and Y. Yamamoto, Ultrafast optical spin echo in a single quantum dot, Nature Photonics 4, 367 (2010).
  * (13) Lord Rayleigh, On the resultant of a large number of vibrations of the same pitch and of arbitrary phase, Phil. Mag. S. 10, 73 (1880).
  * (14) C. R. Kagan, C. B. Murray, and M. G. Bawendi, Long-range resonance transfer of electronic excitations in close-packed CdSe quantum-dot solids, Phys. Rev. B 54, 8633 (1996).
  * (15) T. J. Liptay, L. F. Marshall, P. S. Rao, R. J. Ram, and M. G. Bawendi, Anomalous Stokes shift in CdSe nanocrystals, Phys. Rev. B 76, 155314 (2007).
  * (16) V. Blondot, A. Bogicevic, A. Coste1, C. Arnold, S. Buil, X. Quélin, T. Pons, N. Lequeux, and J.-P. Hermier, Fluorescence properties of self assembled colloidal supraparticles from CdSe/CdS/ZnS nanocrystals, New J. Phys. 22 113026 (2020).
  * (17) V. Blondot, C. Arnold, A. Delteil, A. Bogicevic, T. Pons, N. Lequeux, D. Gérard, S. Buil, J.-P. Hermier, Fluorescence decay enhancement and FRET inhibition in self-assembled hybrid gold CdSe/CdS/CdZnS colloidal nanocrystals supraparticles, arXiv:2207.03135
  * (18) S. A. Empedocles and M. G. Bawendi, Influence of Spectral Diffusion on the Line Shapes of Single CdSe Nanocrystallite Quantum Dots, J. Phys. Chem. B 103, 1826 (1999).
  * (19) L. Coolen, X. Brokmann, P. Spinicelli, and J.-P. Hermier, Emission Characterization of a Single CdSe-ZnS Nanocrystal with High Temporal and Spectral Resolution by Photon-Correlation Fourier Spectroscopy, Phys. Rev. Lett. 100, 027403 (2008).
  * (20) A. Zaccone, Explicit Analytical Solution for Random Close Packing in $d=2$ and $d=3$, Phys. Rev. Lett. 128, 028002 (2022).