# Cosmological coupling of local gravitational systems

Mariano Cadoni  Mirko Pitzalis  Davi C. Rodrigues  Andrea P. Sanna

###### Abstract

We investigate the cosmological coupling of spherical, local astrophysical
systems. We derive a general formula quantifying the cosmological coupling of
the Misner-Sharp mass of these objects. We show that, in the weak-field limit,
the cosmological coupling is only allowed if there are pressure anisotropies.
We also apply our results to galaxies, modelling them with the Navarro-Frenk-
White and Einasto profiles. We show that the galactic mass can be coupled to
the cosmological dynamics and examine its dependence on the scale factor of
the universe.

##  1 Introduction

The possibility of a dynamical coupling between local astrophysical objects
like, e.g., black holes or stars, and the large-scale cosmological dynamics in
a General Relativity (GR) framework has a quite long and partially
contradictory history (see, e.g., Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36]). In the last few years, there have been rather
interesting developments on the issue, both from the theoretical and
observational point of view.

From the theoretical side, two different approaches have emerged. One employs
perturbations and an averaging procedure [37, 38, 39] finding that the
coupling can be traced back to pressures within local compact astrophysical
objects. This effect manifests itself as a power-law relation between the
local object masses and the scale factor $a$, i.e., $M(a)\propto a^{k}$ [37],
with $-3\leq k\leq 3$. These results and the consequent claim that black holes
might be the source of dark energy [40] have been criticised from different
points of view. First, the framework adopted could be flawed from the
beginning [41]. Second, the huge separation of scales between local objects
and the cosmological dynamics should make such coupling physically implausible
[42, 43]. Third, the equation of state of matter inside a black hole could
hardly mimic dark energy [44, 45, 46].

The other approach tackles the cosmological coupling problem using a solid GR
framework [47, 48], along the lines of the first work on the subject by
McVittie [2]. It is based on a general cosmological embedding of local objects
and gives a general description of the coupling between local inhomogeneities
and the cosmological background. The origin of the cosmological coupling is
here traced back to a curvature term in Einstein’s equations, giving the
precise prediction $M(a)\propto a^{k}$ with $k=1$. Moreover, it was shown that
the cosmological coupling is not generic for compact objects, but limited to
specific ones characterised by a mismatch between their ADM and Misner-Sharp
(MS) masses [48]. A striking consequence of this result is that, while usual
singular black holes do not couple to the cosmological background, non-
singular ones do. Thus, an observational evidence of a cosmological coupling
of astrophysical black holes with coupling exponent $k=1$ would be the smoking
gun of their non-singular nature [48].

From the observational side, the situation is quite intricate. The first test
of the cosmological mass shift formula was performed in Ref. [40] using a
sample of supermassive black holes at the centre of elliptical galaxies at
different redshifts. This set of data showed a preference for $k\sim 3$ [40,
47], which eventually motivated the speculative connection between dark energy
and black holes. Other observational constraints on the scaling parameter $k$
are rather controversial, as they seem to strongly depend on the astrophysical
sample employed in the analysis [49, 50, 51, 52, 53].

Leaving aside these issues on the observational side, the main goal of this
paper is to generalize the results of Refs. [47, 48] to generic, local, static
and spherically-symmetric, astrophysical objects, in particular those in the
weak-field regime. Indeed, a limitation of the results of Refs. [47, 48] is
that the cosmological mass shift formula derived there applies only to
strongly-coupled, compact objects, characterized by a linear relation between
their mass $M$ and radius $L$. Therefore, it is of interest to generalize
these results to generic local objects, without any assumption on their
compactness. Another issue is related to understanding how specific is the
embedding of astrophysical bodies in a Friedmann-Lemaître-Robertson-Walker
(FLRW) cosmological background. Specifically, one would like to understand
whether generic local gravitational systems, especially Newtonian ones, allow
for cosmological embedding. Answering this question is crucial to shed light
on the possible cosmological coupling of astrophysical systems of interest,
such as galaxies, galaxy clusters or binary systems, which are known to be in
the weak-field regime. Addressing the above issues is further complicated by
the fact that the existence of cosmologically-embedded solutions typically
requires an anisotropic fluid as local source. Thus, the weak-field limit of
gravitational systems characterized by this kind of sources needs to be kept
under control.

In this paper, we address these issues building on the results of Refs. [47,
48]. We first generalize the mass formula of Ref. [47] to generic spherically-
symmetric, local astrophysical objects (section 2). Then, we delve into the
weak-field limit of these solutions sourced by an anisotropic fluid and
discuss their cosmological embedding in a FLRW background (section 3).
Consequently, we apply our results to the cosmological coupling of galaxies
(section 4). Lastly, we draw our conclusions in section 5.

##  2 Cosmological coupling of generic local astrophysical objects

In Refs. [47, 48], a mass formula describing the cosmological coupling of
strongly-coupled, ultra-compact astrophysical objects was derived. A crucial
assumption in the derivation was the linear scaling of the mass of the
spherically-symmetric object with its radius [48]. This is surely the case for
black holes or other compact objects, but not for other astrophysical bodies,
such as stars, galaxies, or galaxy clusters. From general arguments, we expect
the mass of any local, spherically-symmetric solution of GR, allowing for a
cosmological embedding, to depend on the redshift.

In this section, we will generalize the mass formula of Refs. [47, 48] to
generic, spherically-symmetric objects that allow for a cosmological
embedding. Throughout this paper, we will consider the cosmological background
as fixed, i.e., fully determined by the usual FLRW cosmology. This is because
we are only interested in the cosmological embedding of local structures in a
regime sufficiently far away from the transition scale to the large-scale
homogeneity and isotropy. Therefore, throughout the paper, we will assume that
such embedding is permitted. This assumption is rather strong, as it is not
allowed by every static, spherically-symmetric solution of GR [54].
Additionally, we assume that the object is locally static on the time-scale of
the cosmological evolution. In other words, the only time-dependence comes
from its cosmological coupling. This implies the exclusion of any other local
dynamical process, such as astrophysical accretion mechanisms.

Following Refs. [47, 48], the spacetime is described by the metric

$\text{d}s^{2}=a^{2}(\eta)\left[-e^{\alpha(\eta,r)}\text{d}\eta^{2}+e^{\beta(\eta,r)}\text{d}r^{2}+r^{2}\text{d}\Omega^{2}\right]\,,$
(2.1)

where $\eta$ is the conformal time.

We consider the general situation in which the local object is sourced by an
anisotropic fluid, characterized by radial and tangential pressures
$p_{\parallel}$, $p_{\perp}$, respectively. The usual isotropic fluid is
recovered as a particular case by setting $p_{\parallel}=p_{\perp}\equiv p$.

The resulting Einstein’s equations allow for two branches of solutions (for
details, see Ref. [47]). For our purposes, the interesting solutions are those
characterized by $\dot{\alpha}=0$. In this case, the field equations and
stress-energy tensor conservation take the form [47, 48]

$\displaystyle
e^{-\beta(r,\eta)}=e^{-\beta_{0}(r)}\left(\frac{a(\eta)}{a(\eta_{0})}\right)^{k(r)}\,;$
(2.2a)
$\displaystyle\frac{\dot{a}^{2}}{a^{2}}\left(3-r\alpha^{\prime}\right)e^{-\alpha}+\frac{1-e^{-\beta}+r\beta^{\prime}e^{-\beta}}{r^{2}}=8\pi
Ga^{2}\rho\,;$ (2.2b)
$\displaystyle\frac{e^{-\beta}+re^{-\beta}\alpha^{\prime}-1}{r^{2}}+e^{-\alpha}\left(-2\frac{\ddot{a}}{a}+\frac{\dot{a}^{2}}{a^{2}}\right)=8\pi
Ga^{2}p_{\parallel}\,;$ (2.2c) $\displaystyle
p^{\prime}_{\parallel}+\frac{\alpha^{\prime}}{2}\left(\rho+p_{\parallel}\right)+\frac{2}{r}\left(p_{\parallel}-p_{\perp}\right)=0\,,$
(2.2d)

where

$k(r)\equiv r\alpha^{\prime}(r)\,,$ (2.3)

$\eta_{0}$ is a reference time, whereas $\beta_{0}(r)$ is the metric function
evaluated at the latter. $\beta_{0}(r)$ represents more than a simple initial
condition. We are assuming that local profiles of the objects remain the same
throughout the cosmological evolution. Physically, this is consistent with our
assumption of neglecting any other local, time-dependent dynamics of the
object. The latter is then assumed to be static. The only permitted time-
dependence comes from the cosmological expansion, represented by the
factorization in eq. 2.2a. Thus, $\beta_{0}(r)$ gives the static profile of
the local solution at every value of the cosmological time. Notice that the
factorized form (2.2a) is compatible with eqs. 2.2b, 2.2c and 2.2d since the
system is not closed. The latter equations have to be considered as the
definitions of the $\eta$\- and $r$-dependent density and pressures. We still
have the freedom to choose the reference time $\eta_{0}$ and the normalization
of the scale factor at that time. Throughout this paper, we will choose this
reference time as the present one and adopt the standard scale factor
normalization $a(\eta_{0})=a_{0}=1$.

Before solving the relevant time-dependent equations, we consider a simpler
problem. At $\eta_{0}$, assuming that $\dot{a}$ and $\ddot{a}$ are
sufficiently small, eqs. 2.2b, 2.2c and 2.2d give

$\displaystyle\frac{1-e^{-\beta_{0}}+r\beta_{0}^{\prime}e^{-\beta_{0}}}{r^{2}}=8\pi
G\,\rho(r)\,;$ (2.4a)
$\displaystyle\frac{e^{-\beta_{0}}+re^{-\beta_{0}}\alpha^{\prime}-1}{r^{2}}=8\pi
G\,p_{\parallel}(r)\,;$ (2.4b) $\displaystyle
p^{\prime}_{\parallel}+\frac{\alpha^{\prime}}{2}\left(\rho+p_{\parallel}\right)+\frac{2}{r}\left(p_{\parallel}-p_{\perp}\right)=0\,.$
(2.4c)

In the above, $\rho(r)$, $p_{\parallel}(r)$ are the time-independent density
and pressure profiles of the local object at the present time.

One can easily integrate eq. 2.4a to obtain

$e^{-\beta_{0}(r)}=1-\frac{2Gm(r)}{r}\,,\qquad
m(r)=4\pi\,\int^{r}\text{d}\tilde{r}\,\tilde{r}^{2}\rho(\tilde{r})\,.$ (2.5)

Then one can solve eqs. 2.4b and 2.4c for $\alpha(r)$.

Once $\beta_{0}$ is determined, the full cosmologically-embedded solution is
obtained by using eq. 2.2a to get $\beta(\eta,r)$. The time-dependent
pressures and density $p_{\parallel}(\eta,r)$, $\rho(\eta,r)$, sourcing the
cosmologically coupled solution, are then given by eqs. 2.2b and 2.2c.

As a matter of fact, the existence of the time-dependent solution is generally
only guaranteed if the tangential pressure $p_{\perp}$ is not fixed, i.e.,
given by a free function. Conversely, if $p_{\perp}(r)$ is fixed, it may
happen that the conservation equation (2.2d) gives an independent equation,
the full system becomes over-constrained and the cosmologically-embedded
solution does not exist. We will see in the next section that this is crucial
for the existence of cosmologically coupled solutions describing local
gravitational objects in the weak-field regime.

The MS mass of the cosmologically-coupled solution can be inferred from the
integration of eq. 2.2b over the proper volume (see [47, 48]). Neglecting the
cosmological density term, which is assumed to be negligible with respect to
the density of the local object, we get the mass shift between some arbitrary
initial $\eta_{i}$ and final $\eta_{f}$ cosmological times

$M(\eta_{f},L)=M(\eta_{i},L)\frac{a_{f}}{a_{i}}\frac{1-e^{-\beta_{0}(L)}a_{f}^{k(L)}}{1-e^{-\beta_{0}(L)}a_{i}^{k(L)}}\,,$
(2.6)

where $L$ is the radius (or size) of the local object. For black holes, this
is identified as the event horizon radius. For horizonless objects, $L$ is
defined as the radius of the sphere that contains $99\%$ of the object’s mass
[47].

To ease the notation, in the following, we use $k(L)\equiv k_{L}$ and
$\beta_{0}(L)\equiv\beta_{0}$. Equation 2.6 represents the most general
expression describing the cosmological redshifting of the mass of generic,
comologically-coupled local objects.

The mass formula obtained in [47, 48] can be recovered applying eq. 2.6 to
very compact objects, for which the initial mass satisfies $M\propto a_{i}L$.
If the object is endowed with an event horizon, we get a pure linear scaling
of $M(a)$

$M(\eta_{f},L)=M(\eta_{i},L)\frac{a_{f}}{a_{i}}\,,$ (2.7)

which is also a particular case of Croker’s scaling law [37, 38, 55, 39, 56].
Indeed, let $r_{h}$ be the horizon radius, thus in the limit $L\to r_{h}$, one
finds $e^{-\beta_{0}(L)}a^{k_{L}}\to 0$. The same behavior of eq. 2.7 is also
recovered from eq. 2.6 whenever $k_{L}=0$ and regardless of the presence of
horizons.

For generic objects allowing for cosmological embedding, instead, eq. 2.6
holds in its generality, giving also a sublinear behavior. In this case, the
linear behavior (2.7) holds only approximately near $a_{f}\sim a_{i}$. Indeed,
by expanding eq. 2.6 near this point, we get eq. 2.7 as leading term.

###  2.1 Quasi-asymptotic flatness and behavior of the mass formula

At fixed time and at distances much greater than the object radius $L$, but
much smaller than the Hubble radius, the static gravitational field of any
local, isolate, gravitational system must decay to zero similarly to the
Newtonian behavior. This simple physical feature constrains the behavior of
the mass formula (2.6), resulting in $\beta_{0}>0$. We write the mass $M$ as a
function of $a$, for $0\leq a\leq 1$,

$M(a)=\frac{La}{2G}\left(1-e^{-\beta_{0}}a^{k_{L}}\right)\,.$ (2.8)

$M(a)$ can behave in two possible ways:

  1. 1.

If $k_{L}$ is close to zero, it increases monotonically from $0$ to $M(1)$;

  2. 2.

If $k_{L}$ is large, it starts increasing linearly near $a=0$, thus reaching a
maximum at $a_{M}$, and then it goes down to $M(1)$.

The presence of the maximum can be detected by looking for the zeroes of the
first derivative of $M(a)$. One finds

$k_{L}\ln a_{M}=-\ln\left[e^{-\beta_{0}}(k_{L}+1)\right]\,.$ (2.9)

For $a<1$, this equation admits a solution only for

$e^{-\beta_{0}}\geq\frac{1}{k_{L}+1}\,.$ (2.10)

Whenever this inequality is satisfied, case $2$ is realized.

The above analysis shows the non-trivial evolution of the quasi-local mass and
its dependence on the cosmological background. Considering the physical origin
of the cosmological coupling of local structures, this behavior can be traced
back to the local curvature and these results might, therefore, be interpreted
as follows. On the one hand, if eq. 2.10 is satisfied, the local structure
becomes weakly coupled at $r=L$, leading to a decrease in mass when $a$
increases. On the other hand, if eq. 2.10 is violated, the physical system is
still in a strong-coupling regime, in a similar manner to the black-hole case,
implying the increase of $M$ with $a$.

##  3 Cosmological embedding in the weak-field limit

In the previous section, we derived the MS mass for a generic local,
spherically-symmetric object, assuming that a GR solution describing its
cosmological embedding exists. The existence of the latter is a quite involved
issue, as there are no general theorems ensuring that a static, spherically-
symmetric solution of GR can be promoted to the form (2.1). Nevertheless,
cosmological embedded solutions are known in a number of cases, including the
Schwarzschild black hole [2] and regular black holes [47, 57].

Most of the objects considered so far are gravitationally strongly-coupled,
i.e., they cannot be described using the weak-field approximation at scales of
the order of their radius $L$. In this section, we will discuss the
cosmological embedding of local objects in the weak-field regime. Among others
we have galaxies and galaxy clusters.

The investigation of the cosmological embedding in the weak-field limit is
also important from an alternative standpoint. At great distances $r\gg L$,
even strongly-coupled objects allow for a weak-field description in GR.
Therefore, the cosmological embedding in the weak-field approximation is
expected to provide insights about known results on the cosmological embedding
of black holes [48].

We start our investigation by first considering the weak-field limit of the
field equations at constant time (2.4).

The weak-field limit alone is not sufficient for recovering the Newtonian
limit within GR. In particular, the latter requires the pressures (or any
anisotropic stress) to be sufficiently small, so as not to influence the
metric solution directly. If the pressures are not sufficiently small to be
compatible with the Newtonian limit, the potential associated to $g_{00}$ does
not need to satisfy a Poisson equation. Therefore, we deal with three cases of
weak-field limit:

  1. 1.

Newtonian: the weak-field gravitational potential $\phi$ satisfies the
standard Poisson equation and the pressures are negligible (Newtonian fluid
limit), i.e., $p_{\parallel}(r)=p_{\perp}(r)=0$;

  2. 2.

Non-Newtonian, which includes two subcases:
$2.a\,$ Poissonian: $\phi$ still satisfies the Poisson equation, but some
components of the pressure are not zero;
$2.b\,$ Non-Poissonian: $\phi$ does not satisfy the Poisson equation in its
standard form.

We will show below that any isotropic Poissonian weak-field limit is
necessarily Newtonian. Notice also that case $2.a$ implies that the far-field
generated by a gravitational system decays as $1/r$ for $r\to\infty$, but
anisotropies in the pressures are necessary. This is, for instance, the case
of regular black holes with de Sitter core discussed in [47, 57].
Schematically, we will show that

$\displaystyle\mbox{Weak field limit
}\xrightarrow{p_{\parallel}\approx-2p_{\perp}}\mbox{ Poissonian
}\xrightarrow{p_{\parallel}\approx p_{\perp}\approx 0}\mbox{ Newtonian };\,$
(3.1a) $\displaystyle\mbox{Weak field limit
}\xrightarrow{p_{\parallel}\neq-2p_{\perp}}\mbox{ non-Poissonian
}.\quad\quad\quad\quad\quad\quad\quad\,\,\,\,\,$ (3.1b)

To investigate the weak-field limit in the static case, we perturb the metric
functions around the Minkowski background

$\displaystyle e^{\alpha}$ $\displaystyle\simeq
1+2\phi(r)+\mathcal{O}(\phi^{2})$ (3.2a) $\displaystyle e^{-\beta_{0}}$
$\displaystyle\simeq 1+2\psi(r)+\mathcal{O}(\psi^{2})$ (3.2b)

where $\phi(r)$ and $\psi(r)$ are two potentials (the former is the Newtonian
one). We expand all quantities up to the linear order $\mathcal{O}(\phi)$ and
$\mathcal{O}(\psi)$. Therefore

$\displaystyle\alpha$ $\displaystyle=\ln(1+2\phi)\simeq
2\phi+\mathcal{O}(\phi^{2})\,;$ (3.3a) $\displaystyle\beta_{0}$
$\displaystyle=-\ln(1+2\psi)\simeq-2\psi+\mathcal{O}(\psi^{2})\,.$ (3.3b)

At linear order, eqs. 2.4a and 2.4b become

$\displaystyle 4\pi G\rho$
$\displaystyle\simeq-\frac{1}{r^{2}}\partial_{r}\left(\psi r\right)\,;$ (3.4a)
$\displaystyle 4\pi Gp_{\parallel}$
$\displaystyle\simeq\left(\frac{\psi}{r^{2}}+\frac{\phi^{\prime}}{r}\right)\,.$
(3.4b)

The conservation equation (2.4c) finally reads

$\left(r^{2}p_{\parallel}\right)^{\prime}\simeq\frac{1}{4\pi
G}\left(\psi^{\prime}+r\phi^{\prime\prime}+\phi^{\prime}\right)\simeq
2rp_{\perp}\,,$ (3.5)

where the higher-order terms, originating from
$\alpha^{\prime}(\rho+p_{\parallel})$, have been neglected. Let us now discuss
separately the different weak-field regimes.

###  3.1 Newtonian and Poissonian weak-field limits

The Poisson equation holds in its usual form in both the Poissonian and
Newtonian weak-field limits:

$\nabla^{2}\phi=4\pi G\rho\,.$ (3.6)

As a result, we can classify the weak-field limit either as Newtonian or
Poissonian depending on the pressures $p_{\parallel}(r)$, $p_{\perp}(r)$.
Comparing eqs. 3.4a and 3.6 yields

$\nabla^{2}\phi=\frac{1}{r^{2}}\partial_{r}(r^{2}\phi^{\prime})=-\frac{1}{r^{2}}\partial_{r}\left(\psi
r\right)\,,$ (3.7)

which implies the following relation between the potentials

$\psi=\frac{c_{1}}{r}-r\phi^{\prime}\,,$ (3.8)

with $c_{1}$ an integration constant. Accordingly, the pressure
$p_{\parallel}$ is

${}p_{\parallel}=\frac{c_{1}}{4\pi Gr^{3}}\,.$ (3.9)

From eq. 3.5 we finally get

$p_{\perp}\simeq-\frac{c_{1}}{8\pi Gr^{3}}=-\frac{1}{2}p_{\parallel}\,.$
(3.10)

Depending on the value of the integration constant $c_{1}$, we can have two
subcases.

####  3.1.1 Case $c_{1}=0$: Newtonian weak-field limit

From eqs. 3.9 and 3.10, we see that this case corresponds to a Newtonian
fluid, in which the pressures are negligible, i.e.,
$p_{\parallel}=p_{\perp}=0$.

The remaining equations are

$\displaystyle\psi$ $\displaystyle=-r\phi^{\prime}\,;$ (3.11a)
$\displaystyle\rho$ $\displaystyle=\frac{1}{4\pi
Gr^{2}}\left(r^{2}\phi^{\prime}\right)^{\prime}\,.$ (3.11b)

To close the system, the density profile must be specified. We note that, when
$\rho=0$ (the point-like source case), $\phi=-\psi=-c_{2}/r$ (with $c_{2}$
another integration constant), which is the standard Newtonian contribution.

####  3.1.2 Case $c_{1}\neq 0$: Poissonian weak-field limit

In this case, the pressures are not zero and we have

$\displaystyle\psi$ $\displaystyle=\frac{c_{1}}{r}-r\phi^{\prime}\,;$ (3.12a)
$\displaystyle\rho$ $\displaystyle=\frac{1}{4\pi
Gr^{2}}\left(r^{2}\phi^{\prime}\right)^{\prime}\,,$ (3.12b)

whereas the pressure components must satisfy

$p_{\parallel}+2p_{\perp}=0\,.$ (3.13)

Also here the system must be closed by specifying a density profile $\rho(r)$.
An interesting example of a Poissonian weak-field limit, which will also be
considered when we will discuss the case of galaxies, is recovered by choosing
a logarithmic gravitational potential

$\phi=A\ln\left(\frac{r}{\delta}\right)\,,$ (3.14)

where $A$ and $\delta$ are some constants. Equation 3.12a gives

$\psi=-\left(A-\frac{c_{1}}{r}\right)\,,$ (3.15)

while, from the Poisson equation, the density profile reads as

$\rho=\frac{A}{4\pi Gr^{2}}\,.$ (3.16)

From eq. 3.10 the pressures are

$p_{\perp}=-\frac{1}{2}p_{\parallel}=\frac{c_{1}}{4\pi Gr^{3}}\,.$ (3.17)

Summarizing, for a generic anisotropic fluid in the weak-field limit, if we
impose the Poisson equation, we can have: $(1)$ the usual Newtonian limit,
when $p_{\perp}=p_{\parallel}=0$; $(2)$ a non-Newtonian, but still Poissonian,
weak-field limit, in which both pressures do not vanish and satisfy eq. 3.13.

It is important to stress that eq. 3.13 forbids any isotropic Poissonian weak-
field limit different from the Newtonian one $p_{\perp}=p_{\parallel}=0$.

###  3.2 Non-Poissonian weak-field limit

If we do not impose the validity of the Poisson equation in the form of eq.
3.6, the weak-field approximation gives

$\nabla^{2}\phi=4\pi G\left(\rho+p_{\parallel}+2p_{\perp}\right)\,,$ (3.18)

i.e., the Poisson equation is sourced by the active mass and we have a generic
non-Poissonian weak-field limit. Notice that the Poissonian and Newtonian
weak-field limits can be obtained as particular cases of eq. 3.18, by setting
$p_{\parallel}+2p_{\perp}=0$ or $p_{\parallel}=p_{\perp}=0$, respectively.

###  3.3 Cosmological embedding

Having under control the weak-field limit of static, spherically-symmetric
solutions, we are ready to embed these solutions into a fixed FLRW
cosmological background. In order to do this, we will use a perturbative
expansion of the $r-$dependent part of the metric functions $\beta$ and
$\alpha$ into the system (2.2).

This approach relies on the possibility of separating the space- and time-
dependent components of the metric functions. Actually, in the case under
consideration, this is indeed possible since the metric function $\alpha$ does
not depend on time, whereas the time dependence of $\beta$ occurs in the
simple functional form (2.2a).

Notice also that, in our perturbative approach, we will completely neglect the
gauge fixing issues. This approach is expected to be consistent as long as we
are only interested in the existence of cosmologically embedded solutions.

By virtue of eq. 2.2a, the first order in the perturbations gives

$\begin{split}e^{-\beta(r,\eta)}&\simeq 1+2r\phi^{\prime}\,\ln a+2\psi\,,\\\
\beta(r,\eta)&\simeq-2\psi-2r\phi^{\prime}\,\ln a\,,\\\
\beta^{\prime}&\simeq-2\psi^{\prime}-2r\phi^{\prime\prime}\,\ln
a-2\phi^{\prime}\,\ln a\,,\\\ e^{-\alpha}&\simeq 1-2\phi\,.\end{split}$ (3.19)

By inserting these expansions into eqs. 2.2b and 2.2c, we obtain

$\displaystyle\frac{\dot{a}^{2}}{a^{2}}\left(3-2r\phi^{\prime}-6\phi\right)-\frac{2\psi+4r\phi^{\prime}\ln
a+2r\psi^{\prime}+2r^{2}\phi^{\prime\prime}\ln a}{r^{2}}\simeq 8\pi
Ga^{2}\rho\,,$ (3.20a) $\displaystyle\frac{2r\phi^{\prime}\ln
a+2\psi+2r\phi^{\prime}}{r^{2}}+(1-2\phi)\left(\frac{\dot{a}^{2}}{a^{2}}-2\frac{\ddot{a}}{a}\right)\simeq
8\pi Ga^{2}p_{\parallel}\,.$ (3.20b)

Let us first consider a static, isotropic, Newtonian source, for which
$p_{\parallel}(r)\simeq p_{\perp}(r)\simeq 0$. The cosmological uplifting of
this local solution can be performed in two different ways: $1.)$ by assuming
a fully isotropic source at all scales, i.e.,
$p_{\parallel}(r,\eta)=p_{\perp}(r,\eta)$ everywhere; $2.)$ by allowing for
$p_{\parallel}(r,\eta)\neq p_{\perp}(r,\eta)$. This second possibility
introduces the difficulty of explaining the physical mechanism producing
anisotropies from isotropic local solutions, which then average again on large
cosmological scales to isotropic configurations, giving rise to the usual FLRW
model. Currently, we are not interested in describing these transition
phenomena. In this section, we will only consider the fully isotropic case
$1.)$. In the next section we will discuss solutions with small anisotropies
in the pressures.

Using eq. 3.11a, eqs. 3.20a and 3.20b become

$\displaystyle\frac{\dot{a}^{2}}{a^{2}}\left(3-2r\phi^{\prime}-6\phi\right)+\frac{2\phi^{\prime}}{r}-\frac{4\phi^{\prime}}{r}\ln
a+\frac{2}{r}(\phi^{\prime}+r\phi^{\prime\prime})-2\phi^{\prime\prime}\ln
a\simeq 8\pi Ga^{2}\rho\,;$ (3.21a) $\displaystyle\frac{2\phi^{\prime}}{r}\ln
a+\left(1-2\phi\right)\left(\frac{\dot{a}^{2}}{a^{2}}-2\frac{\ddot{a}}{a}\right)\simeq
8\pi Ga^{2}p_{\parallel}\,.$ (3.21b)

We now plug these into the conservation equation (2.4c), which, in our fully
isotropic case, gives

$\frac{\phi^{\prime\prime}}{r\phi^{\prime}}-\frac{1}{r^{2}}=-\frac{1}{\ln
a}\left(\frac{\dot{a}^{2}}{a^{2}}+\frac{\ddot{a}}{a}\right)\,,$ (3.22)

where we have neglected the higher-order terms and we have assumed
$\phi\neq\text{constant}$. The above is characterized by a complete separation
between the $r$\- and $\eta$-dependent terms. Generally, such equation would
have to be solved by setting both terms equal to the same constant. However,
this would determine the form of the scale factor. Since we are working at
fixed, but fully arbitrary, FLRW cosmological background, the only possibility
for the above equation to be consistent is that $\phi=\text{constant}$, which
however corresponds to the Minkowski solution. This shows that a gravitational
system, described by the usual Newtonian weak-field limit with a fully
isotropic source, does not allow for a cosmological embedding in FLRW
cosmologies, unless the spherical object induces a local anisotropy in the
cosmological fluid.

The derivation above only holds for fully isotropic, Newtonian systems. In
particular, it does not apply to Poissonian or non-Poissonian objects. Indeed,
in these two cases, eq. 2.2d does not give an independent equation for $\phi$,
$\psi$, $\rho$, but it must be considered as an additional equation defining
the tangential pressure $p_{\perp}$. Thus, generally, non-Newtonian systems in
the weak-field regime will allow for a cosmological embedding.

Notice that the above statement remains true for any weak-field non-Newtonian
system, even if we impose for the local system an equation of state
$p_{\parallel}=p_{\parallel}(\rho(r))$ or a given profile for $p_{\perp}(r)$
(this is the case of galaxies, which will be discussed in the next section).
In fact, in this case, $\phi$ and $\psi$ are independent quantities and eqs.
3.20a, 3.20b and 2.2d simply define the (time-dependent) quantities
$\rho(\eta,r)$, $p_{\parallel}(\eta,r)$ and $p_{\perp}(\eta,r)$, respectively.

In order to elucidate these points, let us apply the same method used for the
Newtonian system to the non-Newtonian weak-field solution given by eqs. 3.14,
3.15, 3.16 and 3.17. Using eqs. 3.14 and 3.15 we get from eq. 2.2d,

$-\frac{c_{1}}{r^{2}}a^{2}-2r\left[-1+\ln\left(\frac{r}{\delta}\right)\right]\dot{a}^{2}+ra\left[-1+4\ln\left(\frac{r}{\delta}\right)\right]\ddot{a}=8\pi
Ga^{4}\,r\,p_{\perp}\,.$ (3.23)

This equation has always a solution which determines $p_{\perp}$.

####  3.3.1 Anisotropic Newtonian systems

In the section above, we have considered the two cases of a fully isotropic,
Newtonian system and that of a non-Newtonian, anisotropic system. However, one
can also consider local systems in the weak-field regime characterized by an
isotropic density profile and small anisotropies in the pressures. From a
physical point of view, these small anisotropies could naturally be induced by
the backreaction between the spherical body and the nearby cosmological fluid.
In order to describe the cosmological coupling of these local systems, we will
decompose the density and the pressure profiles of the source as follows (see
also Ref. [54])

$\displaystyle\rho(r,\,\eta)$
$\displaystyle=\frac{\rho(r)}{a(\eta)^{2}}+\rho_{\mathrm{c}o}(\eta)+\rho^{+}(r,\eta)\,;$
(3.24a) $\displaystyle p_{\parallel}(r,\,\eta)$
$\displaystyle=\frac{p_{\parallel}(r)}{a(\eta)^{2}}+p_{\mathrm{c}o}(\eta)+p_{\parallel}^{+}(r,\eta)\,;$
(3.24b) $\displaystyle p_{\perp}(r,\,\eta)$
$\displaystyle=\frac{p_{\perp}(r)}{a(\eta)^{2}}+p_{\mathrm{c}o}(\eta)+p_{\perp}^{+}(r,\eta)\,.$
(3.24c)

In the above, $\rho(r)$, $p_{\parallel}(r)$ and $p_{\perp}(r)$ are given by
the field equations in the static and weak-field limits, namely by eqs. 3.4a,
3.4b and 3.18. $\rho_{\mathrm{c}o}$ and $p_{\mathrm{c}o}$ represent, instead,
the purely cosmological background contributions, satisfying the Friedmann
equations

$\displaystyle a^{2}\rho_{\mathrm{c}o}$ $\displaystyle=\frac{3}{8\pi
G}\left(\frac{\dot{a}}{a}\right)^{2}\,;$ (3.25a) $\displaystyle
a^{2}p_{\mathrm{c}o}$ $\displaystyle=\frac{1}{8\pi
G}\left[\left(\frac{\dot{a}}{a}\right)^{2}-2\frac{\ddot{a}}{a}\right]\,.$
(3.25b)

Finally, the cross contributions (denoted with the superscript $+$) are the
ones which allow for small anisotropies induced on the object when
$p_{\parallel}^{+}\neq p_{\perp}^{+}$.

The existence of the cosmologically embedded solutions can be shown by
deriving the explicit form of $\rho^{+}$, $p_{\parallel}^{+}$,
$p_{\perp}^{+}$. Assuming, for simplicity, the cosmological fluid to be dust,
i.e., $p_{\mathrm{c}o}=0$, from eq. 2.2c, together with eqs. 3.25b and 3.4b,
we obtain

$p_{\parallel}^{+}(r,\eta)=\frac{\phi^{\prime}}{4\pi
Gra^{2}}\ln\left(\frac{a}{a_{0}}\right)\,.$ (3.26)

From eq. 2.2b, instead, we have

$\rho^{+}(r,\eta)=-\frac{\nabla^{2}\phi}{4\pi
Ga^{2}}\ln\left(\frac{a}{a_{0}}\right)-\frac{2}{3}\frac{\rho_{\mathrm{c}o}}{r^{2}}\left(r^{3}\phi\right)^{\prime}\,.$
(3.27)

The energy-momentum conservation, to lowest order, gives

$p_{\perp}^{+}(r,\eta)=\frac{1}{8\pi
Gra^{2}}\left\\{\left[1+\ln\left(\frac{a}{a_{0}}\right)\right](\phi^{\prime}+r\phi^{\prime\prime})+\psi^{\prime}\right\\}+\frac{1}{2}r\phi^{\prime}\rho_{\mathrm{c}o}\,.$
(3.28)

If the static, local object has isotropic pressure, i.e.,
$p_{\parallel}(r)=p_{\perp}(r)\equiv p(r)$, from eq. 2.4c, it can only be
stable if $p(r)$ is much smaller than $\rho(r)$, such that
$p(r)\sim\alpha^{\prime}\rho$ (thus implying the Newtonian limit). Following
the prescription above, such object can only be cosmologically embedded and
continues to be stable (in the sense that $p(r)$ and $\rho(r)$ are kept time-
independent quantities) as long as we assume the presence of small, dynamical
anisotropies in the pressure, i.e., if
$p_{\parallel}^{+}(r,\eta)\not=p^{+}_{\perp}(r,\eta)$. This result is
perfectly consistent with the absence of a cosmological embedding of a purely-
isotropic source at all scales, outlined in the previous subsection.

Instead, if the static central object is characterized by different pressure
components, both of which are subleading with respect to the corresponding
density, we recover again eq. 3.11a. In this limit, eq. 3.27 can be expressed
as

$\rho^{+}(r,\eta)=\frac{\rho(r)}{a^{2}(\eta)}\ln\left(\frac{a_{0}}{a(\eta)}\right)-\frac{2}{3}\frac{\rho_{\mathrm{c}o}(\eta)}{r^{2}}\left(r^{3}\phi\right)^{\prime}\,.$
(3.29)

Assuming that $\rho_{\mathrm{c}o}$ is negligible within the spherical body,
$a_{0}$ fixes the instant in which the additional density correction is zero
($\rho^{+}=0$). As $a(\eta)$ increases, $\rho^{+}$ always decreases and can
become negative. If the $\rho_{\mathrm{c}o}$-term is relevant in the weak-
field limit, $(r^{3}\phi)^{\prime}$ is expected to be always negative (since
$\phi$ would be negative and would approach zero slower than $1/r$). Since
$\rho_{\mathrm{c}o}\propto a^{-3}$, at large redshift this term dominates.
Therefore, qualitatively, the energy density of such object is expected to
increase up to some value of $a$, and then to start decreasing. The turning
point should depend on the cosmological background and the object internal
structure. This is in qualitative agreement with our general results for the
behaviour of the mass $M(a)$ obtained in section 2.1. It will be further
confirmed in section 4.2.

Summing up, in the weak-field limit, a spherical object with fully isotropic
pressure ($p_{\perp}(r,\eta)=p_{\parallel}(r),\eta$), immersed in an isotropic
fluid, cannot be cosmologically coupled. However, the presence of the
spherical body can induce anisotropies in the pressures, i.e.,
$p^{+}_{\parallel}(r,\eta)\not=p^{+}_{\perp}(r,\eta)$. We have explicitly
shown that non-negligible anisotropies are a necessary ingredient to describe,
in a manner consistent with GR, the embedding of local objects into the
cosmological background. At the moment, however, we are unable to describe the
dynamical mechanism leading to such configurations. If this is the case, there
exists a non-trivial cosmological coupling and the MS mass of local objects
will evolve cosmologically.

##  4 Cosmological coupling of galaxies and galaxy clusters

Galaxies are understood as bound gravitational systems in the weak-field
regime. Here, we consider that the dark-matter halo can be approximated by a
spherical matter distribution, as in the standard picture, but which has a
non-negligible relativistic pressure. In the standard picture, dark matter is
made of particles that only interact through gravity (self-interaction being
negligible). In the latter case, the effective pressure of the fluid comes
directly from the velocity dispersion of the particles. It is impossible for
such non-interacting particles to describe a fluid with relativistic pressure
since this would imply relativistic particle velocities, which would far
exceed the escape velocity of the galaxy. Self-interacting dark matter does
not have this limit.

An important application of our results concerns galaxies (elliptical galaxies
in particular) and galaxy clusters. Indeed, these are local objects in the
weak-field regime and can be described, in a first approximation, as static
and spherically-symmetric. Moreover, their rotational curves and the local
velocity dispersion indicate that they can be described as gravitational
systems sourced by an anisotropic fluid [58, 59, 60, 61, 62]. Therefore, in a
GR framework, the local dynamics of galaxy and clusters can be described by
eqs. 2.4a, 2.4b and 2.4c. In all generality, we expect that they couple with
the large-scale cosmological dynamics along the lines discussed in the
previous sections.

To close the system of equations (2.4a), (2.4b), (2.4c), we need information
about their density $\rho(r)$ and pressure ($p_{\perp}$, $p_{\parallel}$)
profiles. Several density profiles for galaxies, devised to reproduce the
observed rotational curves, dark-matter halos and/or coming from $N-$body
simulations, have been proposed in the literature. Here, we will consider
three cases: the simple $1/r^{2}$, the Navarro-Frenk-White (NFW) and the
Einasto profiles.

Providing a description of pressure profiles is much more complicated, as it
is very difficult to relate them to direct observations. Galactic pressure
profiles are phenomenologically related to the velocity dispersion in the
galaxy. A way to codify this information is either through $(1)$ an equation
of state $p_{\parallel}=p_{\parallel}(\rho)$ or $(2)$ an explicit profile
$p_{\perp}(r)$. The simplest assumption one can start from is that tangential
pressures can be neglected:

$p_{\perp}=0\,.$ (4.1)

Within the standard approach to cold dark matter, the velocity dispersion of
the latter is expected to be anisotropic. Hence, this can be seen as a first
approximation [63, 64, 65, 66, 67]. If one considers a two-fluid model for the
galaxy, describing baryonic and dark matter separately, this choice is also
quite natural (see [62] and the discussion therein). By fixing $p_{\perp}=0$,
we are also fixing the profile of $p_{\parallel}$ in the weak-field limit.
Indeed, eq. 3.5 implies that $p_{\parallel}\propto 1/r^{2}$.

For a given density profile $\rho(r)$, the time-independent metric functions
$e^{\alpha}$, $e^{\beta_{0}}$ can be derived using eqs. 2.5 and 2.4 and
imposing $p_{\perp}=0$. While $\beta_{0}(r)$ can be easily obtained by direct
integration of eq. 2.5, $\alpha(r)$, instead, is obtained by adding together
eqs. 2.4a and 2.4b, by differentiating eq. 2.4b and finally by substituting
the results into eq. 2.4c. One gets a Riccati equation for $Z\equiv
r\alpha^{\prime}$, namely

$Z^{\prime}+\frac{Z^{2}}{2r}-\frac{Z\beta_{0}^{\prime}}{2}=\beta_{0}^{\prime}\,.$
(4.2)

Notice that the function $Z(r)$ coincides with the function $k(r)$ (2.3),
which in turn enters the cosmological mass shift formula (2.6).

One can also transform the Riccati equation into a linear second order ODE by
defining $Z=2r(u^{\prime}/u)$

$u^{\prime\prime}+\frac{u^{\prime}}{r}-\frac{\beta_{0}^{\prime}}{2}u^{\prime}-\frac{\beta_{0}^{\prime}}{2r}u=0\,.$
(4.3)

Typically, eq. 4.2 or eq. 4.3 cannot be solved in closed form; however, they
can be solved perturbatively or numerically. Let us discuss separately the
$1/r^{2}$, NFW and Einasto profiles. First, it is instructive to consider the
case of the weak-field limit with a constant equation-of-state parameter.

###  4.1 Weak-field limit with constant equation-of-state parameter

As a first application of the derived cosmological coupling with
$p_{\perp}=0$, we consider the following equation of state

$p_{\parallel}=w\rho\,,$ (4.4)

with $w$ a constant. From eq. 3.5, which assumes the weak-field limit, one
finds $p_{\parallel}\propto 1/r^{2}$. Hence, we write the dark-matter density
as

$\rho(r)=\frac{\rho_{0}}{r^{2}}\,.$ (4.5)

This profile was previously considered in Ref. [62]. It also describes the
large radii behavior of the pseudo-isothermal profile [68].

Figure 1: The MS mass ratio ($M_{f}/M_{i}$) as a function of $a_{f}$, as given
by eq. 4.10. Two different initial redshift values ($z_{i}$) are considered
(10 or 50). For each $z_{i}$ case, two regions are shown, one in gray for
$-1<w<1$ and the other in blue ($z_{i}=50$) or red ($z_{i}=10$) for the
restricted range $0<w<1/3$. The dotted line is the line with $M_{f}/M_{i}=1$.

From the field equation (2.4b), in the weak-field limit we get,

$\beta_{0}(r)=r\alpha^{\prime}(r)-8\pi Gw\rho_{0}\,.$ (4.6)

Substituting this solution into the time-time field equation (2.4a) we have,

$\nabla^{2}\alpha=8\pi G(1+w)\rho(r)\,.$ (4.7)

For $w\approx 0$, $2\alpha$ is the Newtonian potential and one recovers the
Newtonian limit. The case $w=0$ implies $p_{\parallel}=p_{\perp}=0$, which is
not possible since the spherical object would not be stable. For $w>0$, it is
in general possible to find stability with $p_{\parallel}\not=p_{\perp}$,
which allows for cosmological coupling.

The $\alpha(r)$ solution, from eq. 4.7 reads,

$\alpha(r)=8\pi G(1+w)\rho_{0}\ln\frac{r}{r_{0}}\,,$ (4.8)

where $r_{0}$ is an arbitrary constant. In the weak-field limit, as expected,
the time-time metric component can be redefined by an additive constant
without changing the physical results.

With the $\alpha(r)$ solution, it is possible to explicitly compute $k_{L}$
and $\beta_{0}(L)$, which read

$\displaystyle k_{L}=L\alpha^{\prime}(L)=8\pi G(1+w)\rho_{0}\,,$
$\displaystyle\beta_{0}(L)=L\alpha^{\prime}(L)-8\pi Gw\rho_{0}=8\pi
G\rho_{0}\,.$ (4.9)

We point out that both $\beta_{0}$ and $k_{L}$ do not depend on $L$ in this
case. Thus, the MS mass reads

$\displaystyle M(\eta_{f},L)$
$\displaystyle=M(\eta_{i},L)\frac{a_{f}}{a_{i}}\frac{1-e^{-\beta_{0}(L)+k_{L}\ln
a_{f}}}{1-e^{-\beta_{0}(L)+k_{L}\ln a_{i}}}$
$\displaystyle=M(\eta_{i},L)\frac{a_{f}}{a_{i}}\frac{1-e^{8\pi
G\rho_{0}(-1+(1+w)\ln a_{f})}}{1-e^{8\pi G\rho_{0}(-1+(1+w)\ln a_{i})}}\,$
(4.10) $\displaystyle\approx M(\eta_{i},L)\frac{a_{f}}{a_{i}}\frac{1-(1+w)\ln
a_{f}}{1-(1+w)\ln a_{i}}\,.$

The above expression is meaningful only when there is an adequate pressure
value to support the spherical object, so that the case $w=0$, as an exact
expression, is not a valid case; however, small and positive $w$ values are
possible. We add that there are exotic compact objects with $w\approx-1$ that
can be stable and, in this case, eq. 4.10 implies a linear increase of the
mass. Figure 1 shows the evolution of $M_{f}/M_{i}$ with different conditions.

###  4.2 $1/r^{2}$ density profile

In this section, we consider three cases related to this halo profile. First
we consider the weak-field limit, then we consider two classes of solutions
that go beyond the weak-field limit.

####  4.2.1 Exact case

This profile is not able to correctly describe the observed phenomenology near
the galactic center, since $\rho$ diverges as $1/r^{2}$ for $r\to 0$.
Therefore, we will employ it here given its simplicity, as it allows us to
analytically solve our equations. This will enable us to test the cosmological
coupling of non-Newtonian systems in the weak-field limit in a controlled
manner.

Instead of solving the Riccati equation (4.2) or the associated second-order
equation (4.3), in this case $\alpha$ can be found by recasting eqs. 2.4b and
2.4c into the following forms

$\alpha^{\prime}=-\frac{2}{r^{2}}\frac{\frac{\text{d}}{\text{d}r}\left(r^{2}p_{\parallel}\right)}{\rho+p_{\parallel}},\quad-\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}\left(r^{2}p_{\parallel}\right)(r-2GM)=\left(8\pi
Gr^{2}p_{\parallel}+\frac{2GM}{r}\right)(\rho+p_{\parallel})\,.$ (4.11)

We can now use the density profile (4.5) into eq. 2.5 to get the mass function
$m(r)$ and $\beta_{0}$

$m(r)=4\pi\sigma r+m_{B}\,,\qquad e^{-\beta_{0}(r)}=1-\frac{2Gm_{B}}{r}-8\pi
G\sigma\,.$ (4.12)

$m_{B}$ is an integration constant that might be interpreted as the baryonic
matter. This contribution can be neglected as long as we consider large
galactic scales over $\sim 20\ \text{kpc}$ (where rotation curves starts
flattening, according to observations).

A quick inspection of the right equation in eq. 4.11 reveals that there are
two families of solutions.

####  4.2.2 First family: $p_{\parallel}=-\rho$

In this case, $p_{\parallel}=-\rho=-\sigma/r^{2}$ automatically satisfies the
$p_{\parallel}$ equation in eq. 4.11. Therefore, plugging this result into eq.
2.2c gives

$r\alpha^{\prime}(1-8\pi
G\sigma)=0\,\qquad\Rightarrow\qquad\alpha=\text{constant}\,,$ (4.13)

where the constant can be set equal to zero without loss of generality (we can
always perform a time rescaling to absorb it).

####  4.2.3 Second family

The general solution of eq. 4.11 for $p_{\parallel}\neq-\rho$ together with
(4.12) reads as

$p_{\parallel}=-\frac{\sigma}{r^{2}}+\frac{1-8\pi G\sigma}{4\pi
Gr^{2}\ln\left(r/\ell\right)}\,,\qquad
e^{\alpha}=B\,\ln\left(\frac{r}{\ell}\right)\,,$ (4.14)

where $\ell$, $A$ are integration constants. Additionally, in order for
$e^{\alpha}$ to be positive, we must require $r>\ell$, which corresponds to
the regime of validity of this solution (for example, this can be
phenomenologically identified as the regime in which the baryonic matter is
negligible and the rotational curves begin to flatten out).

We see that the solutions of the first family, characterized by
$p_{\parallel}=-\rho$ and a negative pressure, dominate for $r\to\infty$,
i.e., in the transition to the cosmological regime. Conversely, at smaller
(galactic) scales, the second positive term for the pressure dominates.
Physically, this means that, at galactic scales, the hydrostatic equilibrium
is reached when the dark-matter halo balances the gravitational attraction
through positive radial pressure. On the other hand, at great distances, at
the transition to the cosmological regime, the hydrostatic equilibrium is not
reached in the way we are used to. It is a local equilibrium in which both
sides of the hydrostatic equilibrium equation (right equation in eq. 4.11)
vanish separately.

In the generic case, our static solution describes a spacetime with a conical
singularity [62]. In the simplest case given in eq. 4.13, the metric is

$\text{d}s=-\text{d}\eta^{2}+\frac{\text{d}r^{2}}{(1-8\pi
G\sigma)}+r^{2}\text{d}\Omega^{2}\,.$ (4.15)

The solution is physically acceptable in the weak-field limit when the deficit
angle is very small

$\sigma\ll\frac{1}{8\pi G}.$ (4.16)

In this limit, the spacetime can be well approximated by the Minkowski one.

By virializing the galactic motion, one gets the asymptotic constant velocity
$v^{2}=2G\sigma$, and $\sigma\approx ra_{0}/G$, where $a_{0}$ is the typical
acceleration scale [69, 70, 71, 72]. From observations we know that
$a_{0}\approx H=\frac{1}{\hat{L}}$, where $\hat{L}$ is the cosmological
horizon. From condition (4.16), it follows that $r\ll\hat{L}$, which is always
satisfied at galactic scales.

####  4.2.4 Cosmological coupling of galaxies with $1/r^{2}$ density profile

Figure 2: Plot of eq. 4.18 as a function of the scale factor $a$, for
different values of $L/\ell$. The solid curves are obtained by evaluating eq.
4.18 at redshift $z=50$, representing the epoch of formation of old galaxies.
The dashed curves are obtained by evaluating eq. 4.18 at $z=10$, representing
the epoch of formation of young galaxies. We also plot for comparison the
linear scaling $M(a)\propto a$ (dash dotted black line). We set $8\pi
G\sigma\sim 10^{-6}$.

From the results of sections 2 and 3, we can discuss the cosmological coupling
of galaxies and evaluate the mass shift (2.6) for the two families of
solutions found above.

#### First family

In this case, working in the weak-field limit for which eq. 4.16 holds, we can
neglect the deficit angle. Thus, the solution reduces essentially to the
Minkowski space and the density and pressures are zero. The weak-field limit
is Newtonian, so there is no cosmological coupling of the solution.

#### Second family

In this case, the pressure $p_{\parallel}(r)$ and the density are not zero,
the weak-field limit is non-Newtonian, hence the cosmological embedding of the
solution is possible. From eq. 4.14, we have

$k(r)=r\alpha^{\prime}(r)=\frac{2}{\ln\left(r/\ell\right)}\,,$ (4.17)

which should be evaluated at a given galactic radius $r=L$, representing the
size of the galaxy. Regardless of the value of $L$, we see that $k(r)>0$ since
$r>\ell$ (see eq. 4.14). Therefore, the cosmological mass shift (2.6) reads

$\begin{split}\frac{M(a)}{M(a_{i})}=\frac{a}{a_{i}}\frac{1-(1-8\pi
G\sigma)a^{k_{L}}}{1-(1-8\pi G\sigma)a_{i}^{k_{L}}}\,.\end{split}$ (4.18)

$8\pi G\sigma$ is a small parameter but must be kept different from zero. It
acts as a kind of regulator for the the mass $M(a)$, as it prevents the latter
from going to zero at the present time, $a=1$. Additionally, it allows to
recover the linear scaling $M(a)\propto a$ in the limit $k_{L}\to 0$.

Since $e^{-\beta_{0}}$ is very close to $1$, we expect the inequality (2.10)
to hold. Indeed, there is a maximum at

$a_{M}=[1+k_{L}-8\pi
G\sigma(1+k_{L})]^{-1/k_{L}}\simeq(1+k_{L})^{-1/k_{L}}\,,$ (4.19)

where we have neglected the deficit angle. By varying the values of the
galactic radius $L$ and selecting (small) values of the parameter $8\pi
G\sigma$, we show some curves for $M(a)$ in fig. 2.

$M(a)$ behaves in an unphysical way. The mass grows linearly near its
formation and keeps growing in a sub-linear way until it reaches a maximum;
then it starts decreasing. This behavior is a consequence of the validity of
the inequality (2.10) and gives a further indication that such a simple
$r^{-2}$ density profile is not appropriate for a physically realistic
description of galactic dynamics.

###  4.3 Navarro-Frenk-White density profile

The NFW density profile is employed to fit dark-matter profiles derived from
$N-$body simulations. It gives a more realistic description of galactic dark-
matter halos than the rough $r^{-2}$ profile. It reads as

$\rho(r)=\frac{\rho_{c}}{\frac{r}{r_{c}}\left(1+\frac{r}{r_{c}}\right)^{2}}\,,$
(4.20)

where $\rho_{c}$ and $r_{c}$ are two parameters characterizing the model.

Although the $1/r$ behavior of the NFW density profile near the galactic
center represents an improvement with respect to the $1/r^{2}$-divergence of
the model analyzed in the previous section, it still does not provide a
satisfactory description of halos in the bulge region (see Ref. [73]).

$r_{c}$ is the typical length-scale at which dark-matter effects begin to
become relevant, whereas $\rho_{c}$ is related to the critical density of the
universe $\rho_{\text{critic}}=3H^{2}/8\pi G$, where $H$ is the Hubble
constant. Both $\rho_{c}$ and $r_{c}$ are defined in terms of a reference
scale, which is the virial radius, commonly known as $r_{200}$. It physically
represents the radius within which the halo density corresponds approximately
to $200$ times the critical density of the universe. $r_{200}$ is inferred
from the mass of the halo
$M_{200}=200\rho_{\text{critical}}(4\pi/3)r_{200}^{3}$.

The ratio of $r_{200}$ and $r_{c}$ also defines a relevant parameter
$\mathcal{C}\equiv r_{200}/r_{c}$, called the concentration, which can
actually be measured or inferred from numerical simulations. Considering
another parameter $\delta_{c}$ (the overdensity parameter), the NFW profile
can be rewritten as follows

$\frac{\rho(r)}{\rho_{\text{critical}}}=\frac{\delta_{c}}{\frac{r}{r_{c}}\left(1+\frac{r}{r_{c}}\right)^{2}}\,.$
(4.21)

By requiring that, at $r_{200}$, $\rho=200\rho_{\text{critical}}$, we also
obtain the expression of $\delta_{c}$

$\delta_{c}=\frac{200}{3}\frac{\mathcal{C}^{3}}{\log\left(1+\mathcal{C}\right)-\frac{\mathcal{C}}{1+\mathcal{C}}}\,,$
(4.22)

which also defines $\rho_{c}$ as

$\rho_{c}=\delta_{c}\,\rho_{\text{critical}}\,.$ (4.23)

Summarizing, the quantities measured/inferred from numerical simulations are
$\rho_{\text{critical}}$, $r_{200}$, $\delta_{c}$ and $\mathcal{C}$. In terms
of them, we can express the parameters characterizing the NFW profile (4.20):

$r_{c}=\frac{r_{200}}{\mathcal{C}}\,,\qquad\rho_{c}=\delta_{c}\,\rho_{\text{critical}}\,.$
(4.24)

Using $r_{c}$, $\rho_{c}$, we can construct a dimensionless
parameter111Reinstating the speed of light $c$, it should read as $8\pi
G\rho_{c}r_{c}^{2}/c^{2}$.,

$\mathcal{K}\equiv 8\pi G\rho_{c}r_{c}^{2}\,,$ (4.25)

which will be used below to describe the cosmological mass shift of galaxies.

Some values of the parameters $M_{200}$, $r_{200}$, $\mathcal{C}$, $r_{c}$,
$\delta_{c}$, $\mathcal{K}$ are given in table 1 for a sample of galaxies
taken from [74].

$M_{200}[10^{12}\,M_{\odot}]$ | $r_{200}[\text{kpc}]$ | $\mathcal{C}^{-1}=r_{c}/r_{200}$ | $r_{c}[\text{kpc}]$ | $\delta_{c}$ | $\mathcal{K}$
---|---|---|---|---|---
$0.319$ | $177$ | $0.052$ | $9.204$ | $230538$ | $3.46\cdot 10^{-7}$
$2.425$ | $348$ | $0.06$ | $20.88$ | $160060$ | $1.24\cdot 10^{-6}$
$2.552$ | $354$ | $0.124$ | $43.896$ | $26596.3$ | $9.09\cdot 10^{-7}$
$29.67$ | $802$ | $0.131$ | $105.062$ | $23322.9$ | $4.57\cdot 10^{-6}$
$102.68$ | $1213$ | $0.110$ | $133.43$ | $35504.7$ | $1.12\cdot 10^{-5}$
$1109.9$ | $2682$ | $0.151$ | $404.982$ | $16659.5$ | $4.85\cdot 10^{-5}$
$1931.5$ | $3226$ | $0.188$ | $606.488$ | $10014.7$ | $6.53\cdot 10^{-5}$
$3009.7$ | $3740$ | $0.143$ | $534.82$ | $18940.5$ | $9.61\cdot 10^{-5}$
Table 1: Table of astrophysical parameters related to the NFW profile,
according to the data contained in the original paper [74].

Using the NFW density profile (4.20) into eq. 2.5, one easily finds the mass
function $m(r)$ and the metric function $\beta_{0}$

$m(r)=\frac{{\cal
K}r_{c}}{2}\left[\ln\left(1+\frac{r}{r_{c}}\right)-\frac{\frac{r}{r_{c}}}{1+\frac{r}{r_{c}}}\right]\,,\qquad
e^{-\beta_{0}}=1-\frac{{\cal
K}r_{c}}{r}\left[\ln\left(1+\frac{r}{r_{c}}\right)-\frac{\frac{r}{r_{c}}}{1+\frac{r}{r_{c}}}\right],$
(4.26)

where ${\cal K}$ is given by eq. 4.25 and we have fixed the integration
constant such that $m(0)=0$.

In order to find the second metric function $\alpha$, we have either to solve
the Riccati equation (4.2) or the associated second order equation (4.3). This
cannot be done in closed analytic form for the present case. However, to
compute the cosmological mass shift (2.6), we need to evaluate the solution at
the galactic radius $r=L$. Being $L\gg r_{c}$, we just need approximate
solutions for $r/r_{c}\gg 1$. Expanding around $r/r_{c}\to\infty$ we get

$e^{-\beta_{0}}=1-{\cal K}\frac{r_{c}}{r}\ln\frac{r}{r_{c}}+{\cal
O}\left(\frac{r_{c}}{r}\right)\,.$ (4.27)

Using eq. 4.27 into the Riccati equation (4.2), we obtain the approximate
solution for $Z=r\alpha^{\prime}$

$Z(r)=k(r)={\cal K}\,r_{c}\,\frac{\ln\frac{r}{r_{c}}}{r}+{\cal
O}\left(\frac{r_{c}}{r}\right)\,.$ (4.28)

Finally, this gives the exponent $k(L)\equiv k_{L}$ in this approximation

$k_{L}={\cal K}\,\frac{r_{c}}{L}\ln\frac{L}{r_{c}}\,.$ (4.29)

####  4.3.1 Cosmological mass shift

The static spherically-symmetric solution given by eqs. 4.26, 4.27 and 4.28
has a non-Newtonian weak-field limit. According to the results of section 3,
it, thus, allows for a cosmological embedding in a FLRW cosmology. We thus
calculate the cosmological mass shift for the galaxy. By virtue of eqs. 4.27,
4.29 and 2.6, it gives

$M(\eta_{f},L)=M(\eta_{i},L)\frac{a_{f}}{a_{i}}\frac{1-\left[1-k\right]a_{f}^{k_{L}}}{1-\left[1-k_{L}\right]a_{i}^{k_{L}}}\,,$
(4.30)

whereas eq. 2.8 gives

$M(a)=\frac{La}{2G}\left[1-(1-k_{L})e^{k_{L}\ln a}\right]\,,$ (4.31)

with $k_{L}$ given by eq. 4.29.

The behavior of the mass can be easily inferred by inserting eqs. 4.27 and
4.29 into eq. 2.10. One finds that the inequality is always violated, thus the
mass increases monotonically with the scale factor.

Introducing the dimensionless ratio $y=L/r_{c}$, the parameter $k_{L}$ given
by eq. 4.29 can also be written as $k_{L}=\mathcal{K}(\ln y/y)$. Even if the
scale $L$ is not really determined (but presumably it should be of the order
of $r_{200}$), we should expect $y\gg 1$ to hold. Therefore, the function
$y^{-1}\ln y$ is bounded from above by $\text{e}^{-1}\simeq 0.367879$ since it
has a maximum at $\bar{y}=\text{e}$. We show the plot of $k_{L}/\mathcal{K}$
in fig. 3.

Figure 3: Plot of $k_{L}/\mathcal{K}$ as a function of $L/r_{c}$ for the NFW
profile.

In the best case scenario, namely for very massive galaxies (see the last
lines in table 1) and $y\sim\text{e}$ (maximum of fig. 3), we should expect a
$k_{L}$ of the order $10^{-5}$. Actually this estimate is not far from being
true, if we identify $L$ with $r_{200}$. Indeed, if we take, e.g., the last
line in table 1, we have that $y=r_{200}/r_{c}\sim 6.99301$. In this case we
have therefore $k_{L}\simeq 2.7\cdot 10^{-5}$. This also explains why the
inequality (2.10) cannot be satisfied. Indeed, eq. 4.31 has a maximum at
$a_{M}=(k_{L}+k_{L}^{2})^{-1/k_{L}}$, which for very small $k_{L}$ goes far
beyond $a=1$.

Figure 4: Plot of eq. 4.30 as a function of the scale factor $a$ for different
values of $k_{L}$. The solid curves are obtained by evaluating eq. 4.30 at
redshift $z=50$, representing the epoch of formation of old galaxies. The
dashed curves are obtained by evaluating eq. 4.30 at $z=10$, representing the
epoch of formation of young galaxies. We also plot for comparison the linear
scaling $M(a)\propto a$ (dash dotted black line).

In fig. 4 we provide the plots of the function $M(a)$ for different values of
the parameter $k_{L}$. In the physical range $a_{i}<a<1$, the mass grows
monotonically with $a$. The growth starts linearly near $a_{i}$, but becomes
almost immediately strongly sublinear, due to the smallness of $k_{L}$. The
larger values of $k_{L}$ used in the figure are unrealistic. They were adopted
just to get a qualitative idea of the behavior of the mass shift when it
becomes relatively large.

###  4.4 Einasto profiles

Another density profile that can be used to model dark matter halos in
galaxies, which can be roughly described as spherically-symmetric, is the
Einasto profile

$\rho=\rho_{0}\,e^{-\left(r/r_{c}\right)^{1/n}}\,,$ (4.32)

where $n$ is a positive index and $r_{c}$ a characteristic length scale of the
profile. $n=1$ is typically used to describe galactic disks, whereas $n=4$ is
used for elliptical galaxies. In general, we can consider that the mass $M$ of
the halo is contained within a radius $R$, so that the mass function (2.5) can
be written in the more convenient way

$m(r)=M\frac{\Gamma\left(3n\right)-\Gamma\left[3n,\left(r/r_{c}\right)^{1/n}\right]}{\Gamma\left(3n\right)-\Gamma\left[3n,\left(R/r_{c}\right)^{1/n}\right]}\,,$
(4.33)

where $\Gamma(a,x)$ is the incomplete Gamma function.

We will use the same strategy for both cases $n=1$ and $n=4$: first solve the
Riccati equation (4.2) approximately near the origin $r\sim 0$, then the
numerical result will be used inside eq. 2.6 to infer the behavior of the
cosmological mass shift. The first step will then be employed as a boundary
condition at the point $r=\epsilon$ to determine the numerical solution of eq.
4.2 up to a certain distance from the center of the halo. Similarly to the NFW
profile, also here the GR solution possesses a non-Newtonian weak-field limit,
allowing for its cosmological embedding. Let us discuss separately the $n=1$
and $n=4$ cases.

Figure 5: Plot of the numerical solution of the Riccati equation (4.2) with
the Einasto profile (4.32) with indexes $n=1$ (blue curve) and $n=4$ (red
curve).

####  4.4.1 $n=1$

Using eq. 4.33 into eq. 2.5 and expanding near the galactic center, i.e.,
$r=0$, the metric function $\beta$ reads

$e^{-\beta_{0}}\simeq
1+\frac{2MGr^{2}}{3r_{c}^{3}\left[-2+\Gamma\left(3,\frac{R}{r_{c}}\right)\right]}+\mathcal{O}(r^{3})\,,\qquad\beta_{0}\simeq-\frac{2GM\,r^{2}}{3r_{c}^{3}\left[-2+\Gamma\left(3,\frac{R}{r_{c}}\right)\right]}\,.$
(4.34)

It is also easy to derive the corresponding approximate solution of the
Riccati equation (4.2),

$Z(r)\simeq-\frac{4GM\,r}{3r_{c}^{3}\left[-2+\Gamma\left(3,\frac{R}{r_{c}}\right)\right]}\,.$
(4.35)

We have all the ingredients to solve eq. 4.2 numerically. It is of interest to
look into the qualitative behavior of the cosmological mass shift (2.6). We
will consider selected values of the parameters:

$R=100\,M\,;\qquad r_{c}=10\,M\,;\qquad M=1\,;\qquad\epsilon=10^{-6}\,.$
(4.36)

The result of the numerical integration is plotted in fig. 5. Notice the
similarity of the curve in fig. 5 with the behavior of $k_{L}(y)$ of the NFW
profile.

This result is then used to compute the mass shift (2.6). The result is
plotted in fig. 6.

Figure 6: Plot of eq. 4.33 as a function of $a$, for the Einasto profile with
$n=1$. The solid curves are obtained by evaluating eq. 4.33 at redshift
$z=50$, representing the epoch of formation of old galaxies. The dashed curves
are obtained by evaluating eq. 4.33 at $z=10$, representing the epoch of
formation of young galaxies. We also plot for comparison the linear scaling
$M(a)\propto a$ (dash dotted black line).

The mass $M(a)$ grows monotonically in the physical range $a_{i}<a<1$,
similarly to what happens for the NFW profile. However, here the curves are
not only strongly sublinear, but tend to flatten at values of $a$ close to 1.

####  4.4.2 $n=4$

We essentially repeat the same steps used for $n=1$, taking now $n=4$. The
expansion of the metric function $\beta$ near $r=0$ now reads

$e^{-\beta_{0}}\simeq
1+\frac{GM\,r^{2}}{6r_{c}^{3}\left[-39916800+\Gamma\left(12,\frac{R^{1/4}}{r_{c}^{1/4}}\right)\right]}+\mathcal{O}(r^{3})\,,\qquad\beta_{0}\simeq-\frac{GM\,r^{2}}{6r_{c}^{3}\left[-39916800+\Gamma\left(12,\frac{R^{1/4}}{r_{c}^{1/4}}\right)\right]}\,.$
(4.37)

The near-$r=0$ approximate solution of eq. 4.2 for $Z$ now reads

$Z\simeq-\frac{GM\,r^{2}}{6r_{c}^{3}\left[-39916800+\Gamma\left(12,\frac{R^{1/4}}{r_{c}^{1/4}}\right)\right]}\,.$
(4.38)

For the numerical integration, we adopt the same parameters as in eq. 4.36.
The result is plotted in fig. 5.

Finally, the mass shift (2.6) is plotted in fig. 7.

Figure 7: Plot of eq. 4.33 as a function of $a$, for the Einasto profile with
$n=4$. The solid curves are obtained by evaluating eq. 4.33 at redshift
$z=50$, representing the epoch of formation of old galaxies. The dashed curves
are obtained by evaluating eq. 4.33 at $z=10$, representing the epoch of
formation of young galaxies. We also plot for comparison the linear scaling
$M(a)\propto a$ (dash dotted black line).

The behavior of $M(a)$ is again quite similar to the NFW and $n=1$ Einasto
profiles. It is characterized by a monotonic and strongly sublinear growth and
flattening.

##  5 Conclusions

In this paper we have derived a general formula for the MS mass giving the
cosmological coupling of local, static, spherically-symmetric, cosmologically-
embedded astrophysical objects. This formula generalizes the one given in
Refs. [47, 48] for black holes or other ultra-compact objects. We have also
proven general statements concerning the cosmological embedding of weak-field
local systems sourced by anisotropic fluids. We have shown that anisotropies
in the pressures, which could be induced by the influence of the cosmological
background, are essential to allow for a cosmological coupling. Although
adopting simple and not completely realistic models to describe the halo
density profiles, by applying our results to galaxies, we have shown that
their mass can be cosmologically coupled. Our results are rather solid and
general. Generically, the MS mass formula holds true once the cosmological
embedding is possible. The latter follows from a quite general feature, namely
the existence of a non-Newtonian weak-field limit. Moreover, our results are
consistent with what is already known about the cosmological coupling of black
holes and other strongly-coupled objects, for which the anisotropy in the
stress-energy tensor seems to be a necessary condition to allow the
cosmological coupling [47, 48].

Let us briefly discuss some drawbacks of our approach and related open issues.
Apart from black holes, black-hole mimickers and galaxies/galaxy clusters, one
would like to test the cosmological coupling also with other local objects in
the weak-field regime, like stars with anisotropies (see, e.g., Refs. [75, 76,
77, 78, 79] and references therein) or binary systems. However, this may not
be so easy not only because of observational difficulties, but also because of
an inherent limitation of our approach. Our mass formula (2.6) holds for
static (eternal) local objects, i.e., objects with no local time-dependent
dynamics. The only allowed time-dependence is the large-scale, cosmological
one. This is a rather strong limitation. It excludes, for instance, binary
systems in which the period of rotation changes over time due gravitational-
wave emission. Incidentally, the inclusion of local time-dependent dynamics,
together with the cosmological one in the system (2.2), seems to be a
formidable task, unless the time-scale of the local dynamics is several orders
of magnitude less than the cosmological one.

From the observational point of view, we expect the detection of the
cosmological coupling of galaxies or other objects in the weak-field regime to
be challenging. This is due to the uncertainties in measuring galactic masses,
which are well above 10%. In addition, it is quite difficult to find
homogeneous samples at different redshifts. Therefore, the search for other
physical parameters, which are sensitive to the cosmological mass shift, is of
paramount importance. They might provide us with indirect proof of the
cosmological coupling of objects in the weak-field regime.

## Acknowledgements

We thank Luca Amendola, Paolo Serra and Andrea Tarchi for fruitful
discussions. DCR thanks Heidelberg University, Germany, and Universidade
Federal do Rio de Janeiro, Brazil, for hospitality and support. He also
acknowledges partial support from Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq-Brazil) and Fundação de Amparo à Pesquisa e
Inovação do Espírito Santo (FAPES-Brazil).

## References

  * Gao and Li [2023] S.-J. Gao and X.-D. Li, Astrophys. J. 956, 128 (2023), arXiv:2307.10708 [astro-ph.HE] .
  * McVittie [1933] G. C. McVittie, Mon. Not. Roy. Astron. Soc. 93, 325 (1933).
  * Einstein and Straus [a] A. Einstein and E. G. Straus, Rev. Mod. Phys. 17 (1945), 120-124 (a), 10.1103/RevModPhys.17.120.
  * Einstein and Straus [b] A. Einstein and E. G. Straus, Rev. Mod. Phys. 18 (1946), 148-149 (b), 10.1103/RevModPhys.18.148.
  * Pachner [1963] J. Pachner, Phys. Rev. 132, 1837 (1963).
  * Dicke and Peebles [1964] R. H. Dicke and P. J. E. Peebles, Physical Review Letters 12, 435 (1964).
  * Vaidya and Yashodhara [1968] P. C. Vaidya and P. S. Yashodhara, Tensor (Japan) 19, 191 (1968).
  * D’Eath [1975] P. D. D’Eath, Phys. Rev. D 11, 1387 (1975).
  * Gautreau [1984] R. Gautreau, Phys. Rev. D 29, 198 (1984).
  * Cooperstock  _et al._ [1998] F. I. Cooperstock, V. Faraoni, and D. N. Vollick, Astrophys. J. 503, 61 (1998), arXiv:astro-ph/9803097 .
  * Nayak  _et al._ [2001] K. R. Nayak, M. A. H. MacCallum, and C. V. Vishveshwara, Phys. Rev. D 63, 024020 (2001), arXiv:gr-qc/0006040 .
  * Baker [2000] G. A. Baker, (2000), arXiv:astro-ph/0003152 .
  * Bolen  _et al._ [2001] B. Bolen, L. Bombelli, and R. Puzio, Class. Quant. Grav. 18, 1173 (2001), arXiv:gr-qc/0009018 .
  * Dominguez and Gaite [2001] A. Dominguez and J. C. Gaite, EPL 55, 458 (2001), arXiv:astro-ph/0106199 .
  * Ellis [2002] G. F. R. Ellis, Int. J. Mod. Phys. A 17, 2667 (2002), arXiv:gr-qc/0102017 .
  * Gao and Zhang [2004] C. J. Gao and S. N. Zhang, Phys. Lett. B 595, 28 (2004), arXiv:gr-qc/0407045 .
  * Sheehan and Kriss [2004] D. P. Sheehan and V. G. Kriss, (2004), arXiv:astro-ph/0411299 .
  * Nesseris and Perivolaropoulos [2004] S. Nesseris and L. Perivolaropoulos, Phys. Rev. D 70, 123529 (2004), arXiv:astro-ph/0410309 .
  * Sultana and Dyer [2005] J. Sultana and C. C. Dyer, Gen. Rel. Grav. 37, 1347 (2005).
  * Li and Wang [2007] Z.-H. Li and A. Wang, Mod. Phys. Lett. A 22, 1663 (2007), arXiv:astro-ph/0607554 .
  * Adkins  _et al._ [2007] G. S. Adkins, J. McDonnell, and R. N. Fell, Phys. Rev. D 75, 064011 (2007), arXiv:gr-qc/0612146 .
  * McClure and Dyer [2006] M. L. McClure and C. C. Dyer, Class. Quant. Grav. 23, 1971 (2006).
  * Sereno and Jetzer [2007] M. Sereno and P. Jetzer, Phys. Rev. D 75, 064031 (2007), arXiv:astro-ph/0703121 .
  * Faraoni and Jacques [2007] V. Faraoni and A. Jacques, Phys. Rev. D 76, 063510 (2007), arXiv:0707.1350 [gr-qc] .
  * Balaguera-Antolinez and Nowakowski [2007] A. Balaguera-Antolinez and M. Nowakowski, Class. Quant. Grav. 24, 2677 (2007), arXiv:0704.1871 [gr-qc] .
  * Mashhoon  _et al._ [2007] B. Mashhoon, N. Mobed, and D. Singh, Class. Quant. Grav. 24, 5031 (2007), arXiv:0705.1312 [gr-qc] .
  * Carrera and Giulini [2010] M. Carrera and D. Giulini, Rev. Mod. Phys. 82, 169 (2010), arXiv:0810.2712 [gr-qc] .
  * Gao  _et al._ [2011] C. Gao, X. Chen, Y.-G. Shen, and V. Faraoni, Phys. Rev. D 84, 104047 (2011), arXiv:1110.6708 [gr-qc] .
  * Faraoni  _et al._ [2014] V. Faraoni, A. F. Zambrano Moreno, and A. Prain, Phys. Rev. D 89, 103514 (2014), arXiv:1404.3929 [gr-qc] .
  * Kopeikin [2015] S. Kopeikin, Eur. Phys. J. Plus 130, 11 (2015), arXiv:1407.6667 [astro-ph.CO] .
  * Faraoni  _et al._ [2015] V. Faraoni, M. Lapierre-Léonard, and A. Prain, JCAP 10, 013 (2015), arXiv:1508.01725 [gr-qc] .
  * Mello  _et al._ [2017] M. M. C. Mello, A. Maciel, and V. T. Zanchin, Phys. Rev. D 95, 084031 (2017), arXiv:1611.05077 [gr-qc] .
  * Faraoni [2018] V. Faraoni, Universe 4, 109 (2018), arXiv:1810.04667 [gr-qc] .
  * Guariento  _et al._ [2019] D. C. Guariento, A. Maciel, M. M. C. Mello, and V. T. Zanchin, Phys. Rev. D 100, 104050 (2019), arXiv:1908.04961 [gr-qc] .
  * Spengler  _et al._ [2022] F. Spengler, A. Belenchia, D. Rätzel, and D. Braun, Class. Quant. Grav. 39, 055005 (2022), arXiv:2109.03280 [gr-qc] .
  * Agatsuma [2022] K. Agatsuma, Phys. Dark Univ. 38, 101134 (2022), arXiv:2211.15668 [gr-qc] .
  * Croker and Weiner [2019] K. S. Croker and J. L. Weiner, Astrophys. J. 882, 19 (2019), arXiv:2107.06643 [gr-qc] .
  * Croker  _et al._ [2020a] K. S. Croker, K. A. Nishimura, and D. Farrah, Astrophys. J. 889 (2020) no.2, 115 889, 115 (2020a).
  * Croker  _et al._ [2020b] K. S. Croker, J. Runburg, and D. Farrah, Astrophys. J. 900, 57 (2020b).
  * Farrah  _et al._ [2023] D. Farrah, K. S. Croker, M. Zevin, G. Tarlé, V. Faraoni, S. Petty, J. Afonso, N. Fernandez, K. A. Nishimura, C. Pearson, L. Wang, D. L. Clements, A. Efstathiou, E. Hatziminaoglou, M. Lacy, C. McPartland, L. K. Pitchford, N. Sakai, and J. Weiner, The Astrophysical Journal Letters 944, L31 (2023).
  * Mistele [2023] T. Mistele, Res. Notes AAS 7, 101 (2023), arXiv:2304.09817 [gr-qc] .
  * Wang and Wang [2023] Y. Wang and Z. Wang, (2023), arXiv:2304.01059 [gr-qc] .
  * Gaur and Visser [2023] R. Gaur and M. Visser, (2023), arXiv:2308.07374 [gr-qc] .
  * Parnovsky [2023] S. L. Parnovsky, (2023), arXiv:2302.13333 [gr-qc] .
  * Avelino [2023] P. P. Avelino, JCAP 08, 005 (2023), arXiv:2303.06630 [gr-qc] .
  * Dahal  _et al._ [2023] P. K. Dahal, S. Maharana, F. Simovic, I. Soranidis, and D. R. Terno, (2023), arXiv:2312.16804 [gr-qc] .
  * Cadoni  _et al._ [2023] M. Cadoni, A. P. Sanna, M. Pitzalis, B. Banerjee, R. Murgia, N. Hazra, and M. Branchesi, JCAP 11, 007 (2023), arXiv:2306.11588 [gr-qc] .
  * Cadoni  _et al._ [2024] M. Cadoni, R. Murgia, M. Pitzalis, and A. P. Sanna, JCAP 03, 026 (2024), arXiv:2309.16444 [gr-qc] .
  * Rodriguez [2023] C. L. Rodriguez, Astrophys. J. Lett. 947, L12 (2023), arXiv:2302.12386 [astro-ph.CO] .
  * Andrae and El-Badry [2023] R. Andrae and K. El-Badry, Astron. Astrophys. 673, L10 (2023), arXiv:2305.01307 [astro-ph.CO] .
  * Lei  _et al._ [2024] L. Lei, L. Zu, G.-W. Yuan, Z.-Q. Shen, Y.-Y. Wang, Y.-Z. Wang, Z.-B. Su, W.-k. Ren, S.-P. Tang, and H. Zhou, Sci. China Phys. Mech. Astron. 67, 229811 (2024), arXiv:2305.03408 [astro-ph.CO] .
  * Amendola  _et al._ [2024] L. Amendola, D. C. Rodrigues, S. Kumar, and M. Quartin, Mon. Not. Roy. Astron. Soc. 528, 2377 (2024), arXiv:2307.02474 [astro-ph.CO] .
  * Lacy  _et al._ [2024] M. Lacy, A. Engholm, D. Farrah, and K. Ejercito, Astrophys. J. Lett. 961, L33 (2024), arXiv:2312.12344 [astro-ph.CO] .
  * Cadoni and Sanna [2021a] M. Cadoni and A. P. Sanna, Int. J. Mod. Phys. A 36, 2150156 (2021a), arXiv:2012.08335 [gr-qc] .
  * Croker  _et al._ [2021] K. S. Croker, M. J. Zevin, D. Farrah, K. A. Nishimura, and G. Tarle, Astrophys. J. Lett. 921, L22 (2021), arXiv:2109.08146 [gr-qc] .
  * Croker  _et al._ [2019] K. Croker, K. Nishimura, and D. Farrah,  (2019), 10.3847/1538-4357/ab5aff, arXiv:1904.03781 [astro-ph.CO] .
  * Cadoni  _et al._ [2022] M. Cadoni, M. Oi, and A. P. Sanna, Phys. Rev. D 106, 024030 (2022), arXiv:2204.09444 [gr-qc] .
  * Cadoni  _et al._ [2018] M. Cadoni, R. Casadio, A. Giusti, W. Mück, and M. Tuveri, Phys. Lett. B 776, 242 (2018), arXiv:1707.09945 [gr-qc] .
  * Tuveri and Cadoni [2019a] M. Tuveri and M. Cadoni, (2019a), arXiv:1904.08209 [gr-qc] .
  * Tuveri and Cadoni [2019b] M. Tuveri and M. Cadoni, Phys. Rev. D 100, 024029 (2019b), arXiv:1904.11835 [gr-qc] .
  * Cadoni  _et al._ [2020] M. Cadoni, A. P. Sanna, and M. Tuveri, Phys. Rev. D 102, 023514 (2020), arXiv:2002.06988 [gr-qc] .
  * Cadoni and Sanna [2021b] M. Cadoni and A. P. Sanna, Class. Quant. Grav. 38, 135004 (2021b), arXiv:2101.07642 [gr-qc] .
  * van der Marel [1993] R. P. van der Marel, 45, 79 (1993).
  * Fukushige  _et al._ [2004] T. Fukushige, A. Kawai, and J. Makino, Astrophys. J. 606, 625 (2004), arXiv:astro-ph/0306203 .
  * Gebhardt  _et al._ [2000] K. Gebhardt, R. Bender, G. Bower, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, C. Grillmair, L. C. Ho, J. Kormendy, T. R. Lauer, J. Magorrian, J. Pinkney, D. Richstone, and S. Tremaine, The Astrophysical Journal 539, L13 (2000).
  * Hernquist [1990] L. Hernquist, Astrophys. J. 356, 359 (1990).
  * Tremaine  _et al._ [2002] S. Tremaine, K. Gebhardt, R. Bender, G. Bower, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, C. Grillmair, L. C. Ho, J. Kormendy, T. R. Lauer, J. Magorrian, J. Pinkney, and D. Richstone, The Astrophysical Journal 574, 740 (2002).
  * Begeman [1989] K. G. Begeman, Astron. & Astrop. 223, 47 (1989).
  * Milgrom [1983a] M. Milgrom, Astrophys. J. 270, 365 (1983a).
  * Milgrom [1983b] M. Milgrom, Astrophys. J. 270, 371 (1983b).
  * de Martino  _et al._ [2020] I. de Martino, S. S. Chakrabarty, V. Cesare, A. Gallo, L. Ostorero, and A. Diaferio, Universe 6, 107 (2020), arXiv:2007.15539 [astro-ph.CO] .
  * McGaugh [2020] S. McGaugh, Galaxies 8, 35 (2020), arXiv:2004.14402 [astro-ph.GA] .
  * Salucci [2019] P. Salucci, Astron. Astrophys. Rev. 27, 2 (2019), arXiv:1811.08843 [astro-ph.GA] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 .
  * Herrera and Santos [1997] L. Herrera and N. O. Santos, Phys. Rept. 286, 53 (1997).
  * Mak and Harko [2003] M. K. Mak and T. Harko, Proc. Roy. Soc. Lond. A 459, 393 (2003), arXiv:gr-qc/0110103 .
  * Cardoso and Pani [2019] V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Raposo  _et al._ [2019] G. Raposo, P. Pani, M. Bezares, C. Palenzuela, and V. Cardoso, Phys. Rev. D 99, 104072 (2019), arXiv:1811.07917 [gr-qc] .
  * Becerra  _et al._ [2024] L. M. Becerra, E. A. Becerra-Vergara, and F. D. Lora-Clavijo, Phys. Rev. D 109, 043025 (2024), arXiv:2401.10311 [astro-ph.HE] .