# Cassini and extra force constraints to nonminimally coupled gravity with
screening mechanism

Riccardo March<EMAIL_ADDRESS>Istituto per le Applicazioni del Calcolo,
CNR, Via dei Taurini 19, 00185 Roma, Italy  INFN - Laboratori Nazionali di
Frascati (LNF), Via E. Fermi 40, Frascati 00044 Roma, Italy  Orfeu Bertolami
<EMAIL_ADDRESS>Departamento de Física e Astronomia,
Faculdade de Ciências da Universidade do Porto, Rua do Campo Alegre 687,
4169-007 , Porto, Portugal  Centro de Física das Universidades do Minho e do
Porto, Rua do Campo Alegre 687, 4169-007, Porto, Portugal  Marco Muccino
<EMAIL_ADDRESS>INFN - Laboratori Nazionali di Frascati (LNF), Via
E. Fermi 40, Frascati 00044 Roma, Italy  Cláudio Gomes<EMAIL_ADDRESS>Departamento das Ciências da Física, Química e Engenharia, Faculdade de
Ciências e Tecnologia da Universidade dos Açores,
Campus de Ponta Delgada, Rua da Mãe de Deus 9500-321 Ponta Delgada, Portugal
Centro de Física das Universidades do Minho e do Porto, Rua do Campo Alegre
687, 4169-007, Porto, Portugal  Simone Dell’Agnello
<EMAIL_ADDRESS>INFN - Laboratori Nazionali di Frascati (LNF),
Via E. Fermi 40, Frascati 00044 Roma, Italy

(today)

###### Abstract

We consider a nonminimally coupled curvature-matter gravity theory at the
Solar System scale. Both a fifth force of Yukawa type and a further non-
Newtonian extra force that arises from the nonminimal coupling are present in
the solar interior and in the solar atmosphere up to interplanetary space. The
extra force depends on the spatial gradient of space-time curvature $R$. The
conditions under which the effects of such forces can be screened by the
chameleon mechanism and be made consistent with Cassini measurement of PPN
parameter $\gamma$ are examined. This consistency analysis requires a specific
study of Sun’s dynamical contribution to the arising forces at all its layers.

##  I Introduction

Cosmological and astrophysical phenomena, such as the accelerated expansion of
the Universe and the flattening of the galaxy rotation curves, can be
explained by resorting to dark energy and dark matter, respectively.
Nevertheless, such phenomena could arise, in principle, from a modification of
General Relativity (GR) at astrophysical (galactic and extragalactic) and
cosmological scales. Among such modifications, $f(R)$ gravity Capoz-1 ;
Carroll ; Capoz-2 ; DeFTs  involves the replacement of the Ricci curvature
scalar $R$, in the Einstein-Hilbert action, by a nonlinear function $f(R)$.

A further modification of GR is nonminimally (NMC) coupled gravity, where the
Einstein-Hilbert action is replaced with a more general form involving two
functions of curvature $f^{1}(R)$ and $f^{2}(R)$ BBHL . The function
$f^{1}(R)$ has a role analogous to $f(R)$ gravity theory, and the function
$f^{2}(R)$ multiplies the matter Lagrangian density giving rise to a
nonminimal coupling between geometry and matter. This possibility has been
extensively studied in the context of dark matter drkmattgal , dark energy
curraccel , inflation c , energy density fluctuations d , gravitational waves
e , cosmic virial theorem f , Jeans’ instability and star formation in Bok
globules g . This model has also been examined with the Newton-Schrodinger
approach h ; i .

In a previous paper MPBD  the case of functions $f^{1}(R),f^{2}(R)$ analytic
at $R=0$ was considered, and constraints to the resulting NMC gravity model
have been computed through perturbations to perihelion precession by using
data from observations of Mercury’s orbit.

It turns out that NMC gravity modifies the gravitational attraction by
introducing both a fifth force of the Yukawa type and an extra force which
depends on the spatial gradient of the Ricci scalar $R$. While the Yukawa
force is typical also of $f(R)$ gravity, the existence of the extra force is
specific of NMC gravity BBHL ; BLP , and it is an effect of the nonminimal
coupling that induces a non-vanishing covariant derivative of the energy-
momentum tensor. The arising Yukawa contribution can give origin to static
solutions even though in the absence of pressure i .

In Ref. MBMBD  constraints to the NMC gravity model with analytic
$f^{1},f^{2}$ functions have been computed by using the results of a
geophysical experiment which looked for deviations from Newton’s inverse
square law in the ocean Zum . It turns out that the presence of the extra
force in a fluid such as seawater imposes more stringent constraints on the
NMC gravity model than the observation of both Mercury’s perihelion precession
and lunar geodetic precession. Hence for the NMC gravity model with analytic
functions, because of the extra force, Solar System constraints are weaker
than geophysical constraints.

In order to look for meaningful Solar System constraints to NMC gravity, in
the present paper we consider the case of a function $f^{2}(R)$ which contains
a term proportional to $R^{\alpha}$, with $\alpha<0$, so that $f^{2}(R)$ is
not analytic at $R=0$. The resulting model has been used in Ref. drkmattgal
to predict the flattening of the galaxy rotation curves, and to predict the
current accelerated expansion of the Universe curraccel .

Since $f^{2}(R)$ is not analytic, the method based on the $1/c$ expansion used
in Ref. MPBD  does not work for the present model. A completely different
nonlinear method has to be used. It turns out that in Solar System the above
NMC model exhibits a screening mechanism, which is the NMC version of the so
called chameleon mechanism KW , which makes it possible to compute Solar
System constraints.

In the present paper we adapt to NMC gravity the nonlinear computations made
in Ref. HS  for the chameleon mechanism in the gravity field of the Sun,
assuming spherical symmetry. The nonminimal coupling complicates the
computation with respect to the case of $f(R)$ gravity considered in Ref. HS ,
though the solution turns out to have essentially the same general qualitative
properties. We compute an analytic approximation of the chameleon solution of
the field equations and we find constraints on the parameters of the NMC
gravity model from the Cassini measurement of PPN parameter $\gamma$ Cassini .
In order to satisfy the Cassini constraint the chameleon solution turns out to
be close to GR inside a _screening radius_ $r_{s}$ that has to be large
enough, particularly, $r_{s}$ either lies inside the solar convection zone,
close to the top of the zone, or it is larger. Deviations from GR are sourced
by the fraction of solar mass, including solar atmosphere, contained in the
region with radii $r>r_{s}$, so that if $r_{s}$ lies in the convection zone
then such deviations are essentially sourced by a _thin shell_ of mass in the
convection zone, which is a typical property of the chameleon mechanism KW .

Moreover, by computing the equations of hydrodynamics in the solar interior
and atmosphere, both modeled as a perfect fluid, we find the expressions of
the fifth force of Yukawa type and of the non-Newtonian extra force which is a
consequence of the nonminimal coupling. The shape of both forces is affected
by the screening mechanism and depends on the dynamical contribution of the
various Sun’s layers. Then we compute the effect of the NMC extra force that
is expected to be relevant where the mass density gradient is large: for
instance, close to the Sun’s edge in the solar atmosphere, particularly in the
chromosphere and in the transition region between the chromosphere and the
inner solar corona. Further constraints on the gravity model are then computed
by resorting to spectroscopic observations of the solar atmosphere.

Before performing our computation, let us point out that in an interesting
recent work Fisher-Carlson  it has been argued that for NMC models where the
coupling function, $f^{2}(R)$, is dominated by negative powers of the scalar
curvature, one should expect sizeable effects of small curvatures at very
dense nuclear physics situations, which demands for a considerable
suppression. This is clearly a quite special situation, as negative scalar
curvature effects were designed to have implications in astrophysical and
cosmological situations. Negative powers models also lead to pathological
situations for black holes for which $R=0$. In fact, the $f^{2}(R)$ function
of the NMC curvature-matter models were not devised to have a single universal
term, but to be a sum of terms that would be relevant at specific scales.
Thus, the point made by of Ref. Fisher-Carlson  must be regarded as the
statement that at nuclear physics scales the NMC effects are suppressed. In
what concerns the calculation that will be carried out in the present work,
this issue can be addressed through the proposed screening mechanism which
would yield a plausible scenario for the desired suppression of the extra
force inside a nucleus. Indeed, the nucleus has to be unscreened, hence the
screening radius is zero, due to its small mass. Thus, with the expression of
the extra force outside the screening radius computed in this paper, the
expression for the density gradient is multiplied by a negative power of
density. In other words, the screening mechanism would suppress the extra
force inside the nucleus where the density is large.

The paper is organized as follows. In Section II the NMC gravity model is
specified. In Section III we compute an approximate solution of the field
equations inside and around the Sun, namely the chameleon solution. In Section
IV we compute the equations of hydrodynamics of a perfect fluid (under the
assumption of spherical symmetry) and we find the expressions of the fifth
force and extra force in the fluid; particularly, we find an approximation of
the extra force in regions of the solar atmosphere where the density gradient
is large. In Section V we compute the constraints from Cassini measurement of
$\gamma$. In Section VI we compute the effect of the extra force in the
chromosphere-corona transition region an we look for constraints from
spectroscopic measurements. Conclusions are drawn in Section 7. A model of
mass density profile, both for the solar interior and the atmosphere, is
reported in an appendix and used to find analytical order of magnitude
estimates of the constraints on the parameters of the gravity model.

##  II Nonminimally coupled gravity

In the present work we consider gravitational theories with an action
functional of the form BBHL ,

$S=\int\left[\frac{1}{2}f^{1}(R)+[1+f^{2}(R)]\mathcal{L}_{m}\right]\sqrt{-g}\,d^{4}x,$
(1)

where $f^{i}(R)$ (with $i=1,2$) are functions of the Ricci scalar curvature
$R$, $\mathcal{L}_{m}$ is the Lagrangian density of matter, and $g$ is the
metric determinant. The standard Einstein-Hilbert action of GR is recovered by
taking

$f^{1}(R)=\frac{c^{4}}{8\pi G}R,\qquad f^{2}(R)=0,$ (2)

where $G$ is Newton’s gravitational constant.

The variation of the action functional with respect to the metric $g_{\mu\nu}$
yields the field equations:

$\displaystyle\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)R_{\mu\nu}-\frac{1}{2}f^{1}g_{\mu\nu}$
(3)
$\displaystyle=\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)+\left(1+f^{2}\right)T_{\mu\nu},$

where $f^{i}_{R}\equiv df^{i}/dR$. The trace of the field equations is given
by

$\displaystyle\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)R-2f^{1}+3\square
f^{1}_{R}+6\square\left(f^{2}_{R}\mathcal{L}_{m}\right)$
$\displaystyle=\left(1+f^{2}\right)T,$ (4)

where $T$ is the trace of the energy-momentum tensor $T_{\mu\nu}$.

A relevant feature of NMC gravity is that the energy-momentum tensor of matter
is not covariantly conserved, indeed, applying the Bianchi identities to Eq.
(3), one finds that

$\nabla_{\mu}T^{\mu\nu}=\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu\nu}\mathcal{L}_{m}-T^{\mu\nu})\nabla_{\mu}R,$
(5)

a result that, as discussed thoroughly in Refs. multiscalar ; Sotiriou1 ,
cannot be “gauged away” by a convenient conformal transformation, but is
instead a distinctive feature of the model under scrutiny.

This property, when applied to the hydrostatic equilibrium of the Sun, will
play an important role in constraining NMC gravity.

###  II.1 Metric and energy-momentum tensors

We use the following notation for indices of tensors: Greek letters denote
space-time indices ranging from 0 to 3, whereas Latin letters denote spatial
indices ranging from 1 to 3. The signature of the metric tensor is
$(-,+,+,+)$.

The Sun is modelled as a stationary distribution of matter with spherical
symmetry. Then the metric which describes the spacetime in the gravitational
field of the Sun, and in the solar neighbourhood of the galaxy, has the
general spherically symmetric isotropic form around a source centred at $r=0$:

$\displaystyle ds^{2}=$ $\displaystyle-$
$\displaystyle\left[1-2\Phi(r)+2\Psi(r)\right]c^{2}dt^{2}$ (6)
$\displaystyle+$
$\displaystyle\left[1+2\Phi(r)\right](dr^{2}+r^{2}d\Omega^{2}),$

where the potentials $\Phi$ and $\Psi$ are perturbations of the Minkowski
metric such that $\left|\Phi(r)\right|\ll 1$ and $\left|\Psi(r)\right|\ll 1$.
For the purpose of the present paper the functions $\Phi$ and $\Psi$ will be
computed at order ${\cal O}(1/c^{2})$. These considerations are consistent
with previous assumptions to tackle the hydrostatic equilibrium in the
relativistic limit stellequil .

The components of the energy-momentum tensor in spherical coordinates, to the
relevant order for our computations and in the case of spherical symmetry and
radial motion, are given by (Ref. Wi , Chapter 4.1):

$\displaystyle T^{tt}$ $\displaystyle=$ $\displaystyle\rho c^{2}+{\cal
O}\left(1\right),\qquad T^{tr}=\rho cv+{\cal O}\left(\frac{1}{c}\right),$ (7)
$\displaystyle T^{rr}$ $\displaystyle=$ $\displaystyle\rho v^{2}+p+{\cal
O}\left(\frac{1}{c^{2}}\right),$ (8) $\displaystyle T^{\theta\theta}$
$\displaystyle=$ $\displaystyle\frac{p}{r^{2}}+{\cal
O}\left(\frac{1}{c^{2}}\right),$ (9) $\displaystyle T^{\varphi\varphi}$
$\displaystyle=$ $\displaystyle\frac{p}{r^{2}\sin^{2}\theta}+{\cal
O}\left(\frac{1}{c^{2}}\right),$ (10) $\displaystyle T^{t\theta}$
$\displaystyle=$ $\displaystyle
T^{t\varphi}=T^{r\theta}=T^{r\varphi}=T^{\theta\varphi}=0,$ (11)

where matter (Sun’s interior and solar atmosphere) is considered as a perfect
fluid with matter density $\rho$, velocity $v$, and pressure $p$.

The Sun is assumed in hydrostatic equilibrium with the exception of the outer
solar corona where the dynamical equilibrium of a steady atmosphere is
considered in order to take into account solar wind Stix . Density and
pressure are radial functions $\rho=\rho(r),p=p(r)$ and velocity $v=v(r)$ in
the outer corona is also radial.

The trace of the energy-momentum tensor is

$T=-\rho c^{2}+{\cal O}\left(1\right).$ (12)

In the present paper we use $\mathcal{L}_{m}=-\rho c^{2}+{\cal O}(1)$ for the
Lagrangian density of matter BLP .

###  II.2 Choice of functions $f^{1}(R)$ and $f^{2}(R)$

Part of computations will be made for general functions $f^{1}(R)$ and
$f^{2}(R)$, while the constraints to NMC gravity will be exploited for the
following specific choice of functions:

$f^{1}(R)=\frac{c^{4}}{8\pi G}R,\qquad
f^{2}(R)=q_{1}R+q_{2}R^{\alpha},\quad\alpha<1,$ (13)

where the function $f^{1}(R)$ corresponds to GR and $q_{1},q_{2}$ and $\alpha$
are real numbers that have to be considered as parameters of the NMC model of
gravity.

The functions (13), with $q_{1}=0$ and negative values of the exponent
$\alpha$, have been used in Ref. drkmattgal  to model the rotation curves of
galaxies, and in Ref. curraccel  to model the current accelerated expansion of
the Universe. The case $q_{2}=0$ and $q_{1}\neq 0$ has been used in Ref.
stellequil  to model stellar equilibrium with computation of solar
observables. If both terms with $q_{1}$ and $q_{2}$ are present, then in Refs.
drkmattgal ; curraccel  the authors conjecture, for $\alpha<0$, that the term
with coefficient $q_{2}$ dominates at large distances and low densities, while
the term with $q_{1}$ dominates inside astrophysical objects with high
densities such as the Sun.

In the present paper we consider the simultaneous presence of both terms and
we find, for negative $\alpha$, that the term in $f^{2}(R)$ with coefficient
$q_{1}$ dominates in the Sun’s interior, while the parameter $q_{2}$ becomes
important in the constraint by Cassini measurement. The results are in
agreement with the above conjecture.

##  III Approximate solution of the field equations

We begin by adapting to NMC gravity the method of solution applied in HS  to
$f(R)$ gravity. Since the metric potentials $\Phi$ and $\Psi$ are small
perturbations we neglect the higher order terms that include products of
potentials or their derivatives, and cross-products between their derivatives
and the potentials. By computing the Ricci tensor and curvature it follows
that the functions $\Phi$ and $\Psi$ satisfy the following equations HS :

$\displaystyle\nabla^{2}(\Phi+\Psi)=-\frac{R}{2},$ (14)
$\displaystyle\nabla^{2}\Psi=-\frac{1}{2}\left(R^{0}_{0}+\frac{R}{2}\right).$
(15)

We introduce the scalar field $\eta$ which is a function of curvature $R$ also
explicitly depending on the radial coordinate $r$:

$\eta=\eta(r,R)=f^{1}_{R}-2f^{2}_{R}\,\rho(r)c^{2}.$ (16)

Aiming to consider the stationary case we have

$\square\eta=\left(1-2\Phi\right)\frac{d^{2}\eta}{dr^{2}}+\left(\frac{2}{r}-\frac{4}{r}\Phi+\frac{d\Psi}{dr}\right)\frac{d\eta}{dr},$
(17)

and neglecting the products of $\Phi$ and $d\Psi/dr$ by the derivatives of
$\eta$, which are second order quantities, we get
$\square\eta\approx\nabla^{2}\eta$.

In this approximation the time-time component of the field equations (3) is
given by

$\eta R_{00}+\frac{1}{2}\,f^{1}=\nabla^{2}\eta+(1+f^{2})\rho c^{2},$ (18)

and the trace (II) of the field equations becomes

$\eta R+3\nabla^{2}\eta-2f^{1}=-(1+f^{2})\rho c^{2}.$ (19)

Following Ref. HS  and combining the time-time component with the trace of the
field equations we obtain

$R^{0}_{0}=\frac{1}{6\eta}\left[2\eta R-f^{1}-4(1+f^{2})\rho c^{2}\right],$
(20)

then Eq. (15) for the potential $\Psi$ becomes

$\nabla^{2}\Psi=-\frac{1}{12\eta}\left[5\eta R-f^{1}-4(1+f^{2})\rho
c^{2}\right].$ (21)

In the sequel the functions $f^{1}(R)$ and $f^{2}(R)$ are required to satisfy
the following conditions in the gravitational field of the Sun that will have
to be verified a posteriori:

$\left|\frac{8\pi G}{c^{4}}\,\frac{f^{1}}{R}-1\right|\ll
1,\qquad\left|f^{2}\right|\ll 1,$ (22)

and the following condition on the derivatives of $f^{1}$ and $f^{2}$ with
respect to $R$,

$\left|\frac{8\pi G}{c^{4}}\eta-1\right|\ll 1.$ (23)

The conditions (22) mean that the Lagrangian density in Eq. (1) is a small
perturbation of the Lagrangian of GR, while the condition (23), in the case
$f^{2}(R)=0$, becomes a condition used in Ref. HS  for $f(R)$ gravity theory.
Using such conditions the equations (14) and (21) for $\Phi$ and $\Psi$ are
approximately given by:

$\displaystyle\nabla^{2}\Phi$ $\displaystyle=$ $\displaystyle-\frac{4\pi
G}{c^{2}}\rho+\frac{1}{6}\left(\frac{8\pi G}{c^{2}}\rho-R\right),$ (24)
$\displaystyle\nabla^{2}\Psi$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(\frac{8\pi G}{c^{2}}\rho-R\right),$ (25)

and the trace (19) of the field equations is approximated by

$\nabla^{2}\eta=\frac{c^{4}}{24\pi G}\,R-\frac{1}{3}\rho c^{2}.$ (26)

Note that the above equations are formally the same as the one found in Ref.
HS  for $f(R)$ gravity, with the difference that the scalar field $\eta$,
defined in (16), depends explicitly on $r$ through the multiplication by
$\rho(r)$ due to the nonminimal coupling. Such dependence on $\rho(r)$ will be
exploited in the sequel.

###  III.1 High-curvature solution

Following the approach commonly used for this kind of problem HS ; KW , we
introduce a potential function $V=V(r,\eta)$ and an effective potential
$V_{\rm eff}=V-\rho\eta c^{2}/3$ such that

$\nabla^{2}\eta=\frac{\partial V_{\rm eff}}{\partial\eta},\qquad\frac{\partial
V}{\partial\eta}=\frac{c^{4}}{24\pi G}\,\omega(\eta,\rho),$ (27)

where the function $\omega(\eta,\rho)$ is obtained by solving the equation
(16) with respect to $R$. Here we assume that such a solution exists and it is
unique: particularly, this property is satisfied for the specific choice of
functions $f^{1}(R),f^{2}(R)$ defined in Eq. (13). Again, the difference with
respect to $f(R)$ gravity consists in the explicit dependence of $\partial
V/\partial\eta$ on $\rho(r)$ due to the nonminimal coupling.

The effective potential has an extremum which corresponds to the GR solution

$R=\omega(\eta,\rho)=\frac{8\pi G}{c^{2}}\rho,$ (28)

which will be called the high-curvature solution as in Ref. HS . We require
that such an extremum is a minimum (see Refs. HS ; KW  and the discussion in
the sequel), which yields the condition

$\frac{\partial^{2}V_{\rm eff}}{\partial\eta^{2}}=\frac{c^{4}}{24\pi
G}\,\frac{1}{\eta_{R}}\geq 0,$ (29)

with $\eta_{R}=f^{1}_{RR}-2f^{2}_{RR}\rho c^{2}$ and $R=\omega(\eta,\rho)$,
the double subscript in $f^{i}_{RR}$ denoting second derivative with respect
to $R$.

Case of specific choice of $f^{1},f^{2}$. The minimum condition (29) is
equivalent to $\eta_{R}\geq 0$ and, for the specific choice (13) of functions
of curvature, the condition becomes

$\alpha(\alpha-1)q_{2}[\omega(\eta,\rho)]^{\alpha-2}\leq 0.$ (30)

In the following we assume $q_{2}\neq 0$, $\alpha\neq 0$ and $\alpha<1$, and
the curvature $R=\omega(\eta,\rho)$ positive in the Solar System, so that the
minimum condition requires

$\begin{cases}0<\alpha<1&\implies q_{2}>0,\\\ \alpha<0&\implies
q_{2}<0.\end{cases}$ (31)

####  III.1.1 Consistency condition in the Sun’s interior

We require the solution $\eta$ of the trace equation to be a perturbation of
the GR solution in part of the Sun’s interior, where we then look for an
approximation of the high-curvature solution. Such a requirement must be met
in order to satisfy the constraint from Cassini measurement of PPN parameter
$\gamma$. Moreover, the finiteness of $\nabla^{2}\eta$ at the origin (Sun’s
center) imposes the boundary condition

$\frac{d\eta}{dr}=0,\quad\mbox{at }r=0.$ (32)

The high-curvature solution (28) is an exact solution of the equation (26) for
$\eta$ only if $\nabla^{2}\eta=0$, hence if $\eta$ is a harmonic function.
Under spherical symmetry, the only harmonic function which satisfies the
boundary condition (32) is a constant and, by definition (16) of the function
$\eta$, it follows that the Sun’s density would satisfy the condition

$f^{1}_{R}-2f^{2}_{R}\,\rho(r)c^{2}=constant.$ (33)

For instance, for the choice (13) of functions $f^{1}(R),f^{2}(R)$ it follows
that the density $\rho(r)$ must also be constant, which is not the case for
the Sun’s interior. Hence the high-curvature solution can only be an
approximate solution of Eq. (26), and a consistency condition for such a
solution is then

$\left|\nabla^{2}\left(\eta(r,R=8\pi
G\rho(r)/c^{2})\right)\right|\ll\frac{1}{3}\rho c^{2}.$ (34)

Computing the Laplacian of $\eta$ according to (16), and setting at the
minimum of $V_{\rm eff}$,

$\frac{\partial^{2}V_{\rm
eff}}{\partial\eta^{2}}=\frac{1}{\lambda^{2}}>0,\quad\mbox{for}\quad\omega(\eta,\rho)=\frac{8\pi
G}{c^{2}}\rho,$ (35)

where $\lambda=\lambda(\rho)>0$ has dimension of length and depends on
density, the consistency condition (34) for the high-curvature solution reads

$\left|(\lambda^{2}-6f^{2}_{R})\nabla^{2}\rho+\frac{8\pi
G}{c^{2}}\left(\frac{d\lambda^{2}}{dR}-12f^{2}_{RR}\right)\left(\frac{d\rho}{dr}\right)^{2}\right|\ll\rho,$
(36)

with $R=8\pi G\rho/c^{2}$.

Case of specific choice of $f^{1},f^{2}$. For the specific choice (13) of
functions of curvature we have

$\lambda^{2}=6q_{2}\alpha(1-\alpha)\left(\frac{8\pi
G}{c^{2}}\rho\right)^{\alpha-1},$ (37)

and the consistency condition Eq. (36) reads

$\left|\lambda^{2}\left[\frac{\alpha}{1-\alpha}\nabla^{2}\rho-\frac{\alpha+1}{\rho}\left(\frac{d\rho}{dr}\right)^{2}\right]+6q_{1}\nabla^{2}\rho\right|\ll\rho.$
(38)

If $\alpha<0$ and the order of magnitude of $|\alpha|$ is unity, and if the
following conditions hold separately:

$\begin{cases}\lambda^{2}|\nabla^{2}\rho|&\ll\rho,\\\
\lambda|d\rho/dr|&\ll\rho,\\\ |q_{1}\nabla^{2}\rho|&\ll\rho,\end{cases}$ (39)

then the consistency condition is satisfied. If $f^{2}(R)=0$ NMC gravity
reduces to $f(R)$ gravity, and in this case the first two conditions of (39)
correspond to the consistency condition found in Ref. HS .

The explicit expression of the solution for $\eta$ in the Sun’s interior will
be found in Section III.1.3.

####  III.1.2 Solution for $\eta$ in the outskirts of the Solar System

The outer coronal atmosphere of the Sun escapes supersonically into
interstellar space giving rise to the solar wind Stix . We approximate the
galactic mass density in the solar neighbourhood of the Milky Way with the
constant value $\rho_{g}\approx 10^{-24}\,{\rm g}/{\rm cm}^{3}$. We denote by
$r_{g}$ a distance from the Sun’s center such that mass density is dominated
by the galactic density component for $r>r_{g}$. We choose $r_{g}$ at the
heliopause, the boundary between the solar wind and the interstellar medium,
corresponding to a heliocentric radial distance of about $120\mbox{ AU
}=2.58\times 10^{4}R_{\odot}$, where $R_{\odot}$ is the Sun’s radius. By using
the density model of the outer solar corona in Section .4.5 of Appendix B, we
find for the electron density of the solar wind, at distance $r_{g}$, the
value $n_{e}\approx 3.24\times 10^{-4}\,{\rm cm}^{-3}$ not far from that
measured by the Voyager spacecrafts Voyag . Across the heliopause a large
(factor of 20 to 50) density increase takes place so that mass density reaches
the galactic value Voyag .

We will assume that in the solar neighbourhood of the galaxy, for $r>r_{g}$,
the spacetime curvature $R$ is approximately given by the GR solution,
$R_{g}=8\pi G\rho_{g}/c^{2}$. Hence, for the gravitational field of the galaxy
the high-curvature solution holds in the solar vicinity for $r>r_{g}$, and the
field $\eta$ approximately minimizes the effective potential $V_{\rm eff}$.
The meaning of such an assumption will be discussed later.

If the condition $R=R_{g}$ is exactly satisfied for $r>r_{g}$, then Eq. (26)
implies that $\eta$ is a harmonic function so that, under spherical symmetry,
$\eta(r)=a_{1}+a_{2}/r$, with $a_{1},a_{2}$ suitable constants, and by
definition (16) of $\eta$,

$\eta(r)=f^{1}_{R_{g}}-2f^{2}_{R_{g}}\,\rho_{g}c^{2}=constant,\quad\mbox{for
}r>r_{g},$ (40)

and $a_{2}=0$. Then, integrating Eq. (26) over spheres of radii $r\geq r_{g}$
centred at $r=0$, and using the divergence theorem, we find

$\int_{0}^{r_{g}}\left(R-\frac{8\pi G}{c^{2}}\,\rho\right)r^{2}dr=0.$ (41)

However, we will find at the end of Section III.2.4 that the presence of the
gravitational field of the Sun makes the integral in (41) strictly negative
(see also HS ), so that equality (41) is not satisfied and the solution with
$R=R_{g}$ for $r>r_{g}$ is not consistent. Then the GR solution can only be an
approximate solution, $R\approx R_{g}$, of Eq. (26) for $r>r_{g}$.

Applying the divergence theorem to Eq. (26) over the sphere with radius
$r_{g}$, we have that the integral in Eq. (41) is negative if and only if

$\frac{d\eta}{dr}(r=r_{g})<0.$ (42)

If we now denote by $\eta_{g}$ the minimizer of the effective potential
$V_{\rm eff}(\eta)$ corresponding to $\rho(r)=\rho_{g}$, since $\eta$
approximately minimizes $V_{\rm eff}$ in the solar neighbourhood of the
galaxy, then Eq. (27), using Eq. (35), becomes

$\nabla^{2}\eta\approx\frac{1}{\lambda_{g}^{2}}(\eta-\eta_{g}),\qquad\mbox{for
}r>r_{g},$ (43)

where $\lambda_{g}=\lambda(\rho_{g})$. The computations in the present paper
will be made under the condition $\lambda_{g}\gg r_{g}$ which will permit us
to find analytic estimates of the results.

We require the gravitational field of the Sun to become negligible in
comparison with the galactic field at large distances from the Sun’s center,
so that we impose the boundary condition

$\eta(r)\approx\eta_{g},\qquad\mbox{for }\frac{r}{\lambda_{g}}\gg 1.$ (44)

The solution of equation (43) with such a boundary condition is the Yukawa
profile

$\eta(r)=C\frac{e^{-r/\lambda_{g}}}{r}+\eta_{g},\qquad r>r_{g},$ (45)

with $C$ constant to be determined by matching the Yukawa profile with the
solution for $r<r_{g}$ that will be computed in the next section. The constant
$C$ measures the deviation of $R$ from $R_{g}$ in the solar neighborhood of
the galaxy. Note that if the extremum of the effective potential were a
maximum, then the Yukawa profile would be replaced by a damped (according to
$1/r$) trigonometric one.

Using $\lambda_{g}\gg r_{g}$ the Yukawa profile is approximated by

$\eta(r)\approx\frac{C}{r}+\eta_{g},\qquad r_{g}<r\ll\lambda_{g}.$ (46)

Note that the condition (42) imposed by the presence of the Sun requires
$C>0$.

Case of specific choice of $f^{1},f^{2}$. For the specific choice (13) of
functions of curvature the minimizer $\eta_{g}$ of the effective potential
$V_{\rm eff}(\eta)$, corresponding to $\rho(r)=\rho_{g}$, is given by

$\eta_{g}=\frac{c^{4}}{8\pi G}\left\\{1-\frac{16\pi
G}{c^{2}}\rho_{g}\left[q_{1}+\alpha q_{2}\left(\frac{8\pi
G}{c^{2}}\rho_{g}\right)^{\alpha-1}\right]\right\\}.$ (47)

Using Eq. (37) we see that for $\alpha<1$ the quantity $\lambda(\rho)$
increases as density decreases, and this is a typical property of the
chameleon mechanism KW . Particularly, the Yukawa range
$\lambda_{g}=\lambda(\rho_{g})$ is an upper bound for $\lambda$ in the Solar
System.

####  III.1.3 Solution for $\eta$ in the Sun’s interior

In the region of the Sun’s interior where the high-curvature solution
approximately holds, according to definition (16), $\eta$ is given by

$\eta(r,R)\approx\eta(r,8\pi G\rho(r)/c^{2}).$ (48)

Case of specific choice of $f^{1},f^{2}$. For the specific choice (13) of
functions of curvature, using (16) and the expression (37) of $\lambda^{2}$,
we have

$\displaystyle\eta$ $\displaystyle\approx$ $\displaystyle\frac{c^{4}}{8\pi
G}-2\left[q_{1}+q_{2}\alpha\left(\frac{8\pi
G}{c^{2}}\rho\right)^{\alpha-1}\right]\rho c^{2}$ (49) $\displaystyle=$
$\displaystyle\frac{c^{4}}{8\pi
G}-2\left[q_{1}+\frac{\lambda^{2}(\rho)}{6(1-\alpha)}\right]\rho c^{2},$

where

$\lambda^{2}(\rho)=\lambda_{g}^{2}\left(\frac{\rho}{\rho_{g}}\right)^{\alpha-1}.$
(50)

Hence for $\alpha<0$ the length $\lambda(\rho)$ becomes very small in the part
of the solar interior where the high-curvature solution holds, although
$\lambda_{g}\gg r_{g}$. For instance, if $\alpha=-1$ and
$\lambda_{g}=10^{3}r_{g}\sim 10^{5}$ AU, at the bottom of the solar convection
zone where $\rho\approx 1.65\times 10^{-1}\,{\rm g}/{\rm cm}^{3}$, we have
$\lambda(\rho)\approx 10^{-7}\,{\rm m}$. At the top of the convection zone
where $\rho\approx 2.73\times 10^{-7}\,{\rm g}/{\rm cm}^{3}$, we have
$\lambda(\rho)\approx 7\times 10^{-2}\,{\rm m}$, hence completely negligible
quantities in both cases (the values of density are computed by using the
density profile in Section .3.1 of Appendix B).

Hence in the case $\alpha<0$ and $|\alpha|$ of order of magnitude unity,
neglecting the contribution of the term proportional to $\lambda^{2}$, in the
high-curvature region we have

$\eta\approx\frac{c^{4}}{8\pi G}-2q_{1}\rho c^{2},$ (51)

so that the linear term $q_{1}R$ in the NMC function $f^{2}(R)$ is dominant in
the Sun’s interior stellequil , where density is large enough. Then, using the
trace equation (26), the curvature $R$ is given by

$R\approx\frac{8\pi G}{c^{2}}\,\rho-48q_{1}\frac{\pi
G}{c^{2}}\,\nabla^{2}\rho,$ (52)

so that, after neglecting $\lambda(\rho)$, the consistency condition (38)
becomes

$6|q_{1}|\,\frac{|\nabla^{2}\rho|}{\rho}\ll 1,$ (53)

which corresponds to the third condition of Eq. (39). By using the model of
mass density profile in the radiative part of the solar interior reported in
Section .3.2 of Appendix B, the maximum of the quantity
$|\nabla^{2}\rho|/\rho$ in the radiative interior is achieved at $r\approx
0.22R_{\odot}$ and it is of order of $10^{2}/R_{\odot}^{2}$, from which it
follows the constraint

$|q_{1}|\ll 10^{-2}R_{\odot}^{2}\qquad\mbox{in the radiative interior.}$ (54)

In the case of the solar convection zone, using the model of mass density
profile reported in Section .3.1 of Appendix B, the maximum of the quantity
$|\nabla^{2}\rho|/\rho$ in the convection zone is achieved at the top of the
zone and it is of order of $10^{7}/R_{\odot}^{2}$, from which it follows the
constraint

$|q_{1}|\ll 10^{-7}R_{\odot}^{2}\qquad\mbox{in the convection zone.}$ (55)

Hence the upper bound on $|q_{1}|$ in the Sun’s interior becomes more
stringent as mass density decreases.

###  III.2 Screening mechanism

In the outer zone of the solar interior (the outer part of the convective
zone), in the solar atmosphere, and in the interplanetary space, the solution
$\eta$ interpolates between the value (48) in the inner Sun’s interior, and
the value (46) in the outskirts of the Solar System, respectively. In order to
compute such interpolating function we adapt to NMC gravity the chameleon
mechanism developed in Ref. KW  and further analysed in Ref. DS .

According to Ref. KW , in the inner zone of the Sun’s interior the solution
$\eta(r)$ remains close to the minimizer of the effective potential $V_{\rm
eff}$ for $r<r_{s}$, where $r_{s}$ is a critical radius, called the screening
radius, that has to be determined. Hence the high-curvature solution holds for
$r<r_{s}$, GR is approximately satisfied, the consistency condition (36) has
also to be satisfied, and the solution for $\eta$ is given by Eq. (48).

In order to compute $\eta$ for $r>r_{s}$, first we integrate $\nabla^{2}\eta$
over the spherical shell of radii $r_{s}$ and $r>r_{s}$, and we use the
divergence theorem:

$4\pi\int_{r_{s}}^{r}\nabla^{2}\eta(r^{\prime})^{2}dr^{\prime}=4\pi\left[\frac{d\eta}{dr}(r)r^{2}-\frac{d\eta}{dr}(r_{s})r_{s}^{2}\right],$
(56)

from which, solving with respect to $d\eta/dr$ and integrating, it follows

$\displaystyle\eta(r)$ $\displaystyle=$
$\displaystyle\eta(r_{s})+\frac{d\eta}{dr}(r_{s})\left(r_{s}-\frac{r_{s}^{2}}{r}\right)$
(57) $\displaystyle+$
$\displaystyle\int_{r_{s}}^{r}\frac{1}{(r^{\prime})^{2}}\int_{r_{s}}^{r^{\prime}}\nabla^{2}\eta(r^{\prime\prime})^{2}dr^{\prime\prime}dr^{\prime},$

then, integrating by parts, using Eq. (27) and the expression of the effective
potential $V_{\rm eff}$, it follows that the function $\eta$ satisfies the
integral equation

$\eta(r)=\eta(r_{s})+\frac{d\eta}{dr}(r_{s})\left(r_{s}-\frac{r_{s}^{2}}{r}\right)+\frac{1}{r}\int_{r_{s}}^{r}\left[\frac{1}{3}\,\rho(r^{\prime})c^{2}-\frac{\partial
V}{\partial\eta}\right]\left(r^{\prime}\right)^{2}dr^{\prime}-\int_{r_{s}}^{r}\left[\frac{1}{3}\,\rho(r^{\prime})c^{2}-\frac{\partial
V}{\partial\eta}\right]r^{\prime}dr^{\prime},$ (58)

where $\eta(r_{s})$ and $d\eta/dr(r_{s})$ are evaluated by using Eq. (48).

Now we proceed to estimate the expression in square brackets in the
integrands. For $r>r_{s}$, in the outer zone of the Sun’s interior and in the
inner solar atmosphere, where mass density is significantly larger than the
galactic density $\rho_{g}$, we require the potential $V(\eta,\rho)$ to
satisfy the condition (see also Ref. KW )

$\left|\frac{\partial V}{\partial\eta}(\eta,\rho)\right|\ll\frac{1}{3}\,\rho
c^{2},$ (59)

as soon as $\eta$ is displaced enough from the minimizer of $V_{\rm eff}$. The
explicit dependence of $\partial V/\partial\eta$ on density $\rho$ is a
distinctive feature of the application of the chameleon mechanism to NMC
gravity with respect to $f(R)$ gravity.

Using Eq. (27), condition (59) reads

$\frac{c^{2}}{8\pi G}\,\omega(\eta,\rho)\ll\rho.$ (60)

This inequality will be verified a posteriori in the case of the specific
choice (13) of functions of curvature. Given that the proof is a bit
technical, it is given in Section .1 of Appendix A.

Since the function $\omega(\eta,\rho)$ is equal to curvature $R$, then $R$ is
small in comparison with the GR solution and we say that Eq. (60) is the
condition for a local low-curvature solution as in Ref. HS . Hence, for
$r>r_{s}$ the solution locally deviates from GR, while deviations from GR are
screened for $r<r_{s}$. We also say that the Sun is screened for $r<r_{s}$.

In the outer solar atmosphere and interplanetary space, where mass density
becomes smaller and gradually approaches the galactic density $\rho_{g}$, we
expand the derivative of the potential around the minimizer $\eta_{g}$:

$\frac{\partial V}{\partial\eta}(\eta,\rho)\approx\frac{\partial
V}{\partial\eta}(\eta_{g},\rho)+\frac{\partial^{2}V}{\partial\eta^{2}}(\eta_{g},\rho)(\eta-\eta_{g}).$
(61)

Case of specific choice of $f^{1},f^{2}$. For the specific choice (13) of
functions of curvature, the function $\omega(\eta,\rho)$, that has to be used
in condition (60), is given by

$\displaystyle\omega(\eta,\rho)$ $\displaystyle=$
$\displaystyle\left(\frac{16\pi G}{c^{2}}\,\alpha
q_{2}\rho\right)^{1/(1-\alpha)}$ (62) $\displaystyle\times$
$\displaystyle\left(1-\frac{16\pi G}{c^{2}}\,q_{1}\rho-\frac{8\pi
G}{c^{4}}\,\eta\right)^{1/(\alpha-1)}.$

####  III.2.1 Solution for $\eta$ in the case of specific functions
$f^{1},f^{2}$

In the sequel we consider the specific choice (13) of functions of curvature.
We assume $\alpha<0$ and $|\alpha|$ of order of magnitude unity, moreover, we
require the consistency condition (38) to be satisfied for $r<r_{s}$ in the
Sun’s interior.

Using now (13) and the definition (16) of $\eta$ we have

$\eta=\frac{c^{4}}{8\pi G}-2\left(q_{1}+\alpha q_{2}R^{\alpha-1}\right)\rho
c^{2}.$ (63)

In Section III.1.3 we have argued that the linear term $q_{1}R$ in the NMC
function $f^{2}(R)$ is dominant in the Sun’s interior, where density is large
enough. Conversely, in regions with low mass density, such as galactic and
interplanetary space and even the solar atmosphere, the effect of the term
$q_{2}R^{\alpha}$ becomes dominant for $\alpha<0$ as it will be shown in the
sequel (see also Ref. stellequil  and the discussion in Ref. curraccel ).
Then, because of the smallness of the spacetime curvature $R$ in the outer
solar atmosphere and interplanetary space, we assume

$|q_{1}|\ll|q_{2}|R^{\alpha-1},$ (64)

so that in this region of space $\eta$ is given by

$\eta\approx\frac{c^{4}}{8\pi G}-2\alpha q_{2}R^{\alpha-1}\rho c^{2}.$ (65)

The validity of assumption (64) will be verified a posteriori by resorting to
the constraint from Cassini measurement. Since the proof is a bit technical it
is given in Section .2 of Appendix A. Solving Eq. (65) with respect to
curvature $R=\omega(\eta,\rho)$ we find

$\omega(\eta,\rho)\approx\left(\frac{16\pi G}{c^{2}}\alpha
q_{2}\rho\right)^{1/(1-\alpha)}\left(1-\frac{8\pi
G}{c^{4}}\eta\right)^{1/(\alpha-1)},$ (66)

from which, using (27) we obtain the property

$\frac{\partial
V}{\partial\eta}(\eta,\rho)\approx\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\frac{\partial
V}{\partial\eta}(\eta,\rho_{g}).$ (67)

Taking now into account that at density $\rho_{g}$ (in the solar vicinity of
the galaxy) the field $\eta$ approximately minimizes the effective potential
$V_{\rm eff}$, so that

$\frac{\partial
V}{\partial\eta}(\eta_{g},\rho_{g})\approx\frac{1}{3}\,\rho_{g}c^{2},$ (68)

we can compute the approximation (61) of the derivative of the potential:

$\displaystyle\frac{\partial V}{\partial\eta}(\eta,\rho)$
$\displaystyle\approx$
$\displaystyle\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\left[\frac{\partial
V}{\partial\eta}(\eta_{g},\rho_{g})+\frac{\partial^{2}V}{\partial\eta^{2}}(\eta_{g},\rho_{g})(\eta-\eta_{g})\right]\approx\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\left[\frac{1}{3}\,\rho_{g}c^{2}+\frac{1}{\lambda_{g}^{2}}(\eta-\eta_{g})\right]$
(69) $\displaystyle\approx$
$\displaystyle\frac{1}{3}\,\rho_{g}c^{2}\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)},$

where we have taken into account that $\lambda_{g}$ is assumed large.

We now proceed to solve the equation for $\eta$ by evaluating the integrals in
Eq. (58). For $r>r_{s}$, in the outer zone of the Sun’s interior and in the
inner solar atmosphere, using Eq. (59) we have

$\frac{1}{3}\,\rho c^{2}-\frac{\partial
V}{\partial\eta}\approx\frac{1}{3}\,\rho c^{2}.$ (70)

In the outer solar atmosphere and interplanetary space, using (69) we have

$\frac{1}{3}\,\rho c^{2}-\frac{\partial
V}{\partial\eta}\approx\frac{1}{3}\,\rho
c^{2}\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right].$
(71)

Since we have

$\alpha<0\implies 0<-\frac{\alpha}{1-\alpha}<1,$ (72)

for $|\alpha|$ not too small, the second term inside the square bracket in Eq.
(71) is negligible in comparison to 1 in the outer zone of the Sun’s interior
and in the inner solar atmosphere, where mass density is significantly larger
than the galactic density $\rho_{g}$, so that Eq. (71) is valid with a good
approximation for any $r>r_{s}$.

We are ready to write the solution for $\eta$ in terms of the screening
radius. Using Eq. (49), for $r<r_{s}$ we have the solution

$\eta=\frac{c^{4}}{8\pi
G}-2\left[q_{1}+\frac{\lambda^{2}(\rho)}{6(1-\alpha)}\right]\rho c^{2},$ (73)

which satisfies the boundary condition (32) by using the expression (37) of
$\lambda^{2}(\rho)$ and the density model of the Sun’s radiative interior in
Section .3.2 of Appendix B, which yields $d\rho/dr=0$ at $r=0$.

Then we introduce the effective mass $M_{\rm eff}(r)$ defined for $r>r_{s}$ as
follows:

$M_{\rm
eff}(r)=4\pi\int_{r_{s}}^{r}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right](r^{\prime})^{2}dr^{\prime}.$
(74)

Substituting now the expression (71) in Eq. (58), and using Eqs. (73-74), we
obtain the solution $\eta$ for $r_{s}<r<r_{g}$:

$\displaystyle\eta(r)$ $\displaystyle=$ $\displaystyle\frac{c^{4}}{8\pi
G}-2c^{2}\left\\{\left[q_{1}+\frac{\lambda_{s}^{2}}{6(1-\alpha)}\right]\rho_{s}+\left[q_{1}+\frac{\alpha\lambda_{s}^{2}}{6(1-\alpha)}\right]\rho_{s}^{\prime}\left(r_{s}-\frac{r_{s}^{2}}{r}\right)\right\\}+\frac{c^{2}}{12\pi}\,\frac{M_{\rm
eff}(r)}{r}$ (75) $\displaystyle-$
$\displaystyle\frac{c^{2}}{3}\int_{r_{s}}^{r}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]r^{\prime}dr^{\prime},$

where $\rho_{s}=\rho(r_{s})$, $\rho_{s}^{\prime}=d\rho/dr(r_{s})$ and
$\lambda_{s}=\lambda(\rho_{s})$. Here we have kept the terms involving
$\lambda_{s}$, notwithstanding that $\lambda_{s}$ is negligibly small (see
Section III.1.3), because such terms will be necessary for the verification a
posteriori of inequality (60), which is a crucial property of a chameleon
solution KW . Since this is the only reason to keep the terms with
$\lambda_{s}$ in the expression of $\eta$, in the following formulae we
neglect such terms.

The solution (73) and (75) is continuous with its derivative at $r=r_{s}$,
moreover, $\eta$ has to be continuous with its derivative at $r=r_{g}$.
Imposing the continuity of the derivative in $r_{g}$, and using the expression
(46) of $\eta$ for $r>r_{g}$, the constant $C$ in Eq. (46) is determined:

$C=\frac{c^{2}}{12\pi}\,M_{\rm
eff}(r_{g})+2c^{2}q_{1}\rho_{s}^{\prime}r_{s}^{2}.$ (76)

Imposing the continuity of $\eta$ in $r_{g}$ the integral identity follows:

$\displaystyle\frac{c^{2}}{8\pi G}$ $\displaystyle-$ $\displaystyle
2q_{1}(\rho_{s}+\rho_{s}^{\prime}r_{s})-\frac{\eta_{g}}{c^{2}}$ (77)
$\displaystyle=$
$\displaystyle\frac{1}{3}\int_{r_{s}}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]rdr.$

The solution $\eta$ is completely determined once the screening radius $r_{s}$
is determined. Then, using the expression (47) of the minimizer $\eta_{g}$ we
have

$\frac{\eta_{g}}{c^{2}}=\frac{c^{2}}{8\pi G}-2q_{1}\rho_{g}-2\alpha
q_{2}\left(\frac{8\pi G}{c^{2}}\right)^{\alpha-1}\rho_{g}^{\alpha}.$ (78)

Substituting in Eq. (77), and taking into account that
$\rho_{g}\ll\rho(r_{s})$, we find an integral equation which determines the
screening radius $r_{s}$ and completes the solution for $\eta$:

$\displaystyle\frac{1}{6}\int_{r_{s}}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]rdr$
$\displaystyle=$ $\displaystyle-q_{1}(\rho_{s}+\rho_{s}^{\prime}r_{s})+\alpha
q_{2}\left(\frac{8\pi G}{c^{2}}\right)^{\alpha-1}\rho_{g}^{\alpha}.$

This is the NMC version of the integral equation found in Ref. DS  for the
chameleon mechanism. Eventually, the solution for $\eta$ is given by formula
(75) for $r_{s}<r<r_{g}$.

####  III.2.2 Verification of inequalities

The solution for $\eta$ has been computed by assuming some inequalities,
necessary in order to find an analytic approximation of the solution, that
have to be verified a posteriori. Thus, we now show that the computed solution
satisfies some of these inequalities (the inequalities that require a
technical proof are verified in Appendix A).

We see immediately that the expressions (73) and (75), which give the solution
for $r\leq r_{g}$, satisfy

$\left|\frac{8\pi G}{c^{4}}\eta-1\right|={\cal
O}\left(\frac{1}{c^{2}}\right)\ll 1,$ (80)

which yields inequality (23). Moreover, using formula (76) for the constant
$C$, and formula (47) for the minimizer $\eta_{g}$ together with the integral
equation (III.2.1), we see that also the solution (46) satisfies the same
inequality for $r>r_{g}$. Hence the computed solution $\eta$ satisfies
inequality (23) for any $r$.

Using now the expression (52) of curvature, for $r<r_{s}$ we have $R={\cal
O}(1/c^{2})$. Then, the integral equation (III.2.1) implies

$|q_{2}|=\rho_{g}^{-\alpha}\cdot{\cal
O}\left(\frac{1}{c^{2-2\alpha}}\right)\ll 1,$ (81)

being $\alpha<0$, from which, using formula (62) for curvature
$R=\omega(\eta,\rho)$ and property (80) for $\eta$, we have $R={\cal
O}(1/c^{2})$ for $r_{s}<r<r_{g}$. Eventually, since $R\approx R_{g}$ for
$r>r_{g}$ this property is satisfied for any $r$.

Now we observe that the quantity $|q_{1}|R$ is largest at the Sun’s centre
and, using inequality (54), we have

$|q_{1}|R(0)\ll 10^{-2}R_{\odot}^{2}\frac{8\pi G}{c^{2}}\,\rho(0)\ll 10^{-5},$
(82)

so that, from definition (13) of $f^{1}$ and $f^{2}$, using property (81) it
follows

$\left|f^{2}(R)\right|=\left|q_{1}R+q_{2}R^{\alpha}\right|<|q_{1}|R(0)+\rho_{g}^{-\alpha}\cdot{\cal
O}\left(\frac{1}{c^{2}}\right)\ll 1,$ (83)

which yields the second of inequalities (22), the first being trivial.

####  III.2.3 Solution for the potentials $\Phi$ and $\Psi$

We have assumed that the high-curvature solution holds for the gravitational
field of the galaxy in the solar neighbourhood for $r>r_{g}$, so that GR is
approximately satisfied. This assumption implies that the Milky Way is
screened within a distance of about $8\,{\rm kpc}$ from its center, where the
Solar System is approximately located. Such a screening condition may impose
additional constraints on the NMC gravity model whose assessment requires the
solution for the gravitational field of the Milky Way, possibly taking also
into account the effect of the other galaxies in the local group, however that
will be the object of future research.

In what follows we denote by $U$ the Newtonian potential of the mass
distribution with density $\rho$,

$U=G\int\frac{\rho(\mathbf{y})}{|\mathbf{x}-\mathbf{y}|}d^{3}y,$ (84)

which satisfies the Poisson equation $\nabla^{2}U=-4\pi G\rho$.

Using Eqs. (24-25) it follows that the potential $\Psi$ of the metric is
related to the deviation from GR, then we impose the following boundary
conditions in the galaxy at large distances from the Sun’s center, where GR is
satisfied by our assumptions:

$\Phi(r)\approx\frac{1}{c^{2}}U(r),\quad\Psi(r)\approx 0,\qquad\mbox{for
}\frac{r}{\lambda_{g}}\gg 1.$ (85)

Combining equations (25) and (26) for $\Psi$ and $\eta$ we have

$\nabla^{2}\left(\Psi+\frac{8\pi G}{c^{4}}\,\eta\right)=0,$ (86)

which implies in the case of spherical symmetry:

$\Psi+\frac{8\pi G}{c^{4}}\,\eta=a_{1}+\frac{a_{2}}{r},$ (87)

with $a_{1},a_{2}$ constants to be determined. Requiring $\Psi$ to be not
singular at $r=0$ imposes $a_{2}=0$.

Using the boundary conditions (44) and (85) for $\eta$ and $\Psi$,
respectively, it then follows

$a_{1}=\frac{8\pi G}{c^{4}}\,\eta_{g},$ (88)

from which we find

$\Psi(r)=-\frac{8\pi G}{c^{4}}\left[\eta(r)-\eta_{g}\right].$ (89)

The solution for $\Psi$ then follows immediately from the solution for $\eta$:
using Eqs. (73), (78) and (64) with $R=R_{g}$, we find for $r<r_{s}$,

$\Psi(r)=16\frac{\pi G}{c^{2}}\left[q_{1}\rho(r)-\alpha
q_{2}R_{g}^{\alpha-1}\rho_{g}\right],$ (90)

using Eqs. (75) and (78-III.2.1), and neglecting $\lambda_{s}$, we find for
$r_{s}<r<r_{g}$,

$\Psi(r)=-16\frac{\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}\frac{r_{s}^{2}}{r}-\frac{2G}{3c^{2}}\,\frac{M_{\rm
eff}(r)}{r}-\frac{8\pi
G}{3c^{2}}\int_{r}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]r^{\prime}dr^{\prime},$
(91)

and, using Eq. (46) and formula (76), we find for $r>r_{g}$,

$\Psi(r)=-\frac{2G}{3c^{2}}\,\frac{M_{\rm eff}(r_{g})}{r}-16\frac{\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}\frac{r_{s}^{2}}{r}.$ (92)

Combining now equations (24) and (25) for $\Phi$ and $\Psi$ we have

$\nabla^{2}\left(\Phi-\frac{U}{c^{2}}-\frac{1}{2}\,\Psi\right)=0.$ (93)

The solution $\Phi$ of this equation which is not singular at $r=0$ and
satisfies the boundary conditions (85) is given by

$\Phi=\frac{U}{c^{2}}+\frac{1}{2}\,\Psi.$ (94)

The solutions found for $\Phi$ and $\Psi$ define the space-time metric (6).

####  III.2.4 Screening radius in the convection zone

In this section we consider the case of the screening radius in the solar
convection zone which is relevant for the constraint from the Cassini
measurement. Using the results of Section III.1.3, we neglect $\lambda$ in the
convection zone so that the consistency condition (38) becomes inequality (53)
which we write in the form

$|q_{1}|\,\frac{|\nabla^{2}\rho|}{\rho}<\varepsilon\ll 1.$ (95)

For $r>R_{\odot}$ we consider the term

$-\frac{2G}{rc^{2}}\left[8\pi
q_{1}\rho_{s}^{\prime}r_{s}^{2}+\frac{1}{3}\,M_{\rm eff}(r)\right]$ (96)

in the expression (91) of $\Psi(r)$. We approximate the effective mass (74)
for $r>R_{\odot}$, and the integral in Eq. (III.2.1), with

$\displaystyle M_{\rm eff}(r)\approx M_{\rm eff}(R_{\odot})\approx
4\pi\int_{r_{s}}^{R_{\odot}}\rho(r)r^{2}dr,$ (97)
$\displaystyle\int_{r_{s}}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]rdr\approx\int_{r_{s}}^{R_{\odot}}\rho(r)rdr,$
(98)

where we have neglected the contribution from the solar atmosphere outside the
photosphere, and we have taken into account that $\rho_{g}\ll\rho(r)$ for
$r\leq R_{\odot}$ and $\alpha<0$.

Then, using inequality (95) and the expressions of $\rho$ and $\nabla^{2}\rho$
in the convection zone reported in Section .3.1 of Appendix B, we find that

$\displaystyle 8\pi\left|q_{1}\rho_{s}^{\prime}\right|r_{s}^{2}$
$\displaystyle<$ $\displaystyle\frac{20}{3}\varepsilon M_{\rm
eff}(R_{\odot}),$ (99)
$\displaystyle|q_{1}|\cdot|\rho_{s}+\rho_{s}^{\prime}r_{s}|$ $\displaystyle<$
$\displaystyle 20\varepsilon\,\frac{1}{6}\int_{r_{s}}^{R_{\odot}}\rho(r)rdr,$
(100)

for any $r_{s}$ in the convection zone. Hence, for $\varepsilon\leq 1/20$ the
quantity (96) is negative and, using (91-92), it follows $\Psi(r)<0$ for any
$r>R_{\odot}$.

Moreover, for such values of $\varepsilon$, using Eq. (76), it follows that
the constant $C$ in Eq. (46) is positive, in agreement with the statement
immediately following Eq. (46). Hence, it follows that condition (42) is
satisfied, so that the integral in Eq. (41) is negative as anticipated in
Section III.1.2.

Eventually, with the above approximations the solution (75) for $\eta$ in the
interval $(R_{\odot},r_{g})$ takes a simpler form, which, also neglecting
$\lambda_{s}$, we write as

$\displaystyle 1-\frac{8\pi G}{c^{4}}\eta$ $\displaystyle\approx$
$\displaystyle 2\alpha q_{2}\left(\frac{8\pi
G}{c^{2}}\rho_{g}\right)^{\alpha}-\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}\frac{r_{s}^{2}}{r}$ (101) $\displaystyle-$
$\displaystyle\frac{2G}{3c^{2}}\,\frac{M_{\rm eff}(R_{\odot})}{r}.$

This approximation will be used in the sequel so to compute in the solar
atmosphere the extra force due to the nonminimal coupling.

##  IV Fifth force and extra force

In this section we compute the equations of hydrodynamics of a perfect fluid
by resorting to the covariant divergence of the energy-momentum tensor BBHL ,
as given by Eq. (5) that we repeat for convenience:

$\nabla_{\mu}T^{\mu\nu}=\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu\nu}\mathcal{L}_{m}-T^{\mu\nu})\nabla_{\mu}R.$
(102)

First we compute the $\nu=t$ component of this equation. Using the components
of the energy-momentum tensor given by Eqs. (7)-(11), the left-hand side of
Eq. (102) yields

$\frac{1}{c}\,\nabla_{\mu}T^{\mu t}=\frac{\partial\rho}{\partial
t}+\frac{\partial}{\partial r}(\rho v)+\frac{2}{r}\,\rho v+{\cal
O}\left(\frac{1}{c^{2}}\right).$ (103)

In order to compute the right-hand side of Eq. (102), first we observe that,
using definition (13) of $f^{2}$ and property (81) of $q_{2}$, we have

$f^{2}_{R}=q_{1}+\alpha q_{2}R^{\alpha-1}=q_{1}+\rho_{g}^{-\alpha}\cdot{\cal
O}(1)={\cal O}(1),$ (104)

with respect to $1/c^{2}$. Then, taking into account that $f^{2}$ satisfies
the second of inequalities (22), the evaluation of the right-hand side of Eq.
(102) yields

$\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu t}\mathcal{L}_{m}-T^{\mu t})\frac{\partial
R}{\partial x^{\mu}}={\cal O}\left(\frac{1}{c}\right).$ (105)

Neglecting terms of order ${\cal O}(1/c^{2})$ the continuity equation, in the
case of spherical symmetry and radial motion, then in the nonrelativistic
limit, it follows as usual:

$\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial r}(\rho
v)+\frac{2}{r}\,\rho v=0.$ (106)

The NMC term on the right-hand side of Eq. (102) gives a distinctive
contribution to the spatial part of this equation that now we compute. For
$\nu=r$, using the components (7)-(11) of the energy-momentum tensor, the
left-hand side yields

$\displaystyle\nabla_{\mu}T^{\mu r}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial t}(\rho v)+\frac{\partial}{\partial
r}(\rho v^{2})-\rho c^{2}\left(\frac{d\Phi}{dr}-\frac{d\Psi}{dr}\right)$ (107)
$\displaystyle+$ $\displaystyle\frac{2}{r}\,\rho v^{2}+\frac{\partial
p}{\partial r}+{\cal O}\left(\frac{1}{c^{2}}\right).$

Using now the continuity equation Eq. (106) at order ${\cal O}(1)$, and the
solution (94) for the metric potential $\Phi$, we get

$\nabla_{\mu}T^{\mu r}=\rho\frac{dv}{dt}+\frac{\partial p}{\partial
r}-\rho\frac{dU}{dr}+\frac{1}{2}\rho c^{2}\frac{d\Psi}{dr}+{\cal
O}\left(\frac{1}{c^{2}}\right),$ (108)

where $d/dt=\partial/\partial t+v\partial/\partial r$ is the time derivative
following the fluid.

The right-hand side of Eq. (102) yields

$\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu r}\mathcal{L}_{m}-T^{\mu r})\frac{\partial
R}{\partial x^{\mu}}=-c^{2}f^{2}_{R}\,\rho\frac{dR}{dr}+{\cal
O}\left(\frac{1}{c^{2}}\right).$ (109)

Eventually, the components $\nu=\theta$ and $\nu=\varphi$ yield terms of order
${\cal O}(1/c^{2})$.

Combining equations (108) and (109), and neglecting terms of order ${\cal
O}(1/c^{2})$, we obtain the radial equation of NMC hydrodynamics for a perfect
fluid in the nonrelativistic limit, which in the stationary case reads

$\rho\frac{dv}{dt}=\rho\frac{dU}{dr}-\frac{dp}{dr}-\frac{1}{2}\rho
c^{2}\frac{d\Psi}{dr}-c^{2}f^{2}_{R}\,\rho\frac{dR}{dr}.$ (110)

We observe the presence of two additional terms in comparison with Eulerian
equations of Newtonian hydrodynamics:

  * (i)

a fifth force density proportional to the gradient of the metric potential
$\Psi$;

  * (ii)

an extra force density proportional to the product of $f^{2}_{R}$ by the
gradient of curvature $R$.

The extra force density in (ii) has been extensively discussed in Ref. BBHL ,
and for relativistic perfect fluids in Ref. BLP . While the fifth force is
typical of $f(R)$ gravity theory, the extra force is specific of NMC gravity.
Now we derive the explicit expressions of such force densities corresponding
to the specific choice (13) of functions of curvature.

###  IV.1 Forces inside the screening radius

If $r<r_{s}$, using the solution (90) for $\Psi$, we find for the fifth force:

$F_{\rm f}=-8\pi Gq_{1}\rho\frac{d\rho}{dr}.$ (111)

In order to compute the extra force we need $f^{2}_{R}$ which, using the
consistency condition (53), is approximated by means of the GR value of
curvature:

$f^{2}_{R}\approx q_{1}+\alpha q_{2}\left(\frac{8\pi
G}{c^{2}}\rho\right)^{\alpha-1}=q_{1}+\frac{\lambda^{2}(\rho)}{6(1-\alpha)}\approx
q_{1},$ (112)

being $\lambda^{2}(\rho)$ negligible inside the screening radius. Using the
solution (52) for the curvature inside the screening radius, we have for the
derivative:

$\frac{dR}{dr}=\frac{8\pi G}{c^{2}}\frac{d\rho}{dr}-48q_{1}\frac{\pi
G}{c^{2}}\frac{d\nabla^{2}\rho}{dr},$ (113)

and neglecting the term with $q_{1}^{2}$ (see also Ref. stellequil  for
analogous computations), we find for the extra force

$F_{\rm e}=-8\pi Gq_{1}\rho\frac{d\rho}{dr},$ (114)

which shows that the fifth force and the extra force are equal inside the
screening radius.

###  IV.2 Forces outside the screening radius

If $r_{s}<r<r_{g}$, using the solution (91) for $\Psi$, we find for the fifth
force:

$F_{\rm f}=-8\pi
Gq_{1}\rho_{s}^{\prime}\rho(r)\left(\frac{r_{s}}{r}\right)^{2}-\frac{G}{3}\rho(r)\frac{M_{\rm
eff}(r)}{r^{2}}.$ (115)

The extra force for $r_{s}<r<r_{g}$ is obtained by computing the derivative of
the curvature $R=\omega(\eta,\rho)$, given by formula (62), and using the
solution (75) for $\eta$ where we neglect the terms with $\lambda_{s}$:

$\displaystyle F_{\rm e}$ $\displaystyle=$
$\displaystyle-\frac{c^{4}}{(1-\alpha)16\pi G}\left(\frac{16\pi
G}{c^{2}}\alpha q_{2}\right)^{\frac{1}{1-\alpha}}\left(1-\frac{8\pi
G}{c^{4}}\eta\right)\left(1-\frac{16\pi G}{c^{2}}q_{1}\rho-\frac{8\pi
G}{c^{4}}\eta\right)^{\frac{\alpha-2}{1-\alpha}}$ (116) $\displaystyle\times$
$\displaystyle\left[\frac{d\rho}{dr}\left(1-\frac{8\pi
G}{c^{4}}\eta\right)-\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}\rho\left(\frac{r_{s}}{r}\right)^{2}-\frac{2G}{3c^{2}}\rho\frac{M_{\rm
eff}(r)}{r^{2}}\right]\rho^{\frac{\alpha}{1-\alpha}}.$

In the case of screening radius in the convection zone, using the
approximation (101) of the solution for $\eta$, we can further approximate the
term between square brackets in the expression of $F_{\rm e}$ as follows:

$\displaystyle-$ $\displaystyle\frac{2G}{rc^{2}}\left(8\pi
q_{1}\rho_{s}^{\prime}r_{s}^{2}+\frac{1}{3}M_{\rm
eff}(R_{\odot})\right)\left(\frac{d\rho}{dr}+\frac{\rho}{r}\right)$ (117)
$\displaystyle+$ $\displaystyle 2\alpha q_{2}R_{g}^{\alpha}\frac{d\rho}{dr}.$

Let us now consider $r$ varying in the solar atmosphere, particularly in
regions close to the Sun’s edge where the density gradient is large. Using the
model of mass density profile in the chromosphere and in the chromosphere-
corona transition region, reported in Section .4.2 and Section .4.3 of
Appendix B, respectively, we observe that in such regions the density gradient
is large enough to have

$\left|\frac{d\rho}{dr}\right|\gg\frac{\rho}{r}.$ (118)

Consequently, the extra force is approximated by

$\displaystyle F_{\rm e}$ $\displaystyle\approx$
$\displaystyle-\frac{c^{4}}{(1-\alpha)16\pi G}\left(\frac{16\pi
G}{c^{2}}\alpha q_{2}\right)^{\frac{1}{1-\alpha}}\left(1-\frac{8\pi
G}{c^{4}}\eta\right)^{2}$ (119) $\displaystyle\times$
$\displaystyle\left(1-\frac{16\pi G}{c^{2}}q_{1}\rho-\frac{8\pi
G}{c^{4}}\eta\right)^{\frac{\alpha-2}{1-\alpha}}\rho^{\frac{\alpha}{1-\alpha}}\,\frac{d\rho}{dr},$

where the expression (101) has to be used for function $\eta$.

In Section .2 of Appendix A, values of $r_{s}$ in the convection zone, $r_{0}$
in the inner region of solar atmosphere, and $\alpha$, are given such that the
term with $q_{1}$ can be neglected for any $r>r_{0}$:

$\frac{16\pi G}{c^{2}}|q_{1}|\rho(r)\ll 1-\frac{8\pi G}{c^{4}}\eta(r).$ (120)

Now we observe that for $r$ in the solar atmosphere and close enough to the
Sun’s edge the quantity $1/r$ is slowly varying in comparison with density
$\rho(r)$ and its derivative, so that we may approximate the expression (101)
of $\eta(r)$ with $\eta(R_{\odot})$. Then we obtain the further approximation
of the extra force:

$F_{\rm
e}\approx\mathcal{F}(\alpha,q_{1},q_{2})\rho^{\frac{\alpha}{1-\alpha}}\,\frac{d\rho}{dr},$
(121)

where the coefficient

$\displaystyle\mathcal{F}(\alpha,q_{1},q_{2})$ $\displaystyle=$
$\displaystyle-\frac{c^{4}}{(1-\alpha)16\pi G}\left(\frac{16\pi
G}{c^{2}}\alpha q_{2}\right)^{\frac{1}{1-\alpha}}$ (122) $\displaystyle\times$
$\displaystyle\left[1-\frac{8\pi
G}{c^{4}}\eta(R_{\odot})\right]^{\frac{\alpha}{\alpha-1}}$

is independent of $r$. This last approximation of the extra force can be used
in regions of the solar atmosphere with high density gradient, such as the
chromosphere and the chromosphere-corona transition region.

##  V Constraint from Cassini measurement

In this section we assume the screening radius in the convection zone and
$\varepsilon$ in inequality (95) such that $\varepsilon\leq 1/20$ so that
$\Psi(r)<0$ for $r>R_{\odot}$. Indeed, we will find that the screening radius
that saturates the Cassini constraint lies in the convection zone. Using the
metric (6) we have for PPN parameter $\gamma$:

$1-\gamma=-\frac{\Psi}{\Phi-\Psi}.$ (123)

Since $\Psi(r)<0$, using formulae (84) and (94) we have $1-\gamma(r)>0$ for
$r>R_{\odot}$.

The most stringent bound on $\gamma$ is given by the Cassini measurement
Cassini :

$\gamma-1=(2.1\pm 2.3)\times 10^{-5},$ (124)

which implies

$0<1-\gamma<2\times 10^{-6}.$ (125)

Using then (94) and (123), and neglecting $1-\gamma$ with respect to $1$
because of the Cassini bound (124), we get

$(1-\gamma)\frac{U}{c^{2}}=-\Psi\left(1-\frac{1-\gamma}{2}\right)\approx-\Psi,$
(126)

from which, using $U(r)=GM_{\odot}/r$ for the Newtonian potential for
$r>R_{\odot}$ in the interplanetary space (neglecting the contribution of the
solar atmosphere), where $M_{\odot}$ is the Sun’s mass, we obtain

$1-\gamma\approx-\frac{c^{2}}{GM_{\odot}}\,r\Psi.$ (127)

Using now the solution (91) for $\Psi$, we have

$\frac{d}{dr}\left(r\Psi(r)\right)<0,\qquad\mbox{for }r_{s}<r<r_{g},$ (128)

from which it follows that $1-\gamma(r)$ is an increasing function of $r$ in
the interplanetary space. For $r>r_{g}$, where the galactic potential $U_{g}$
dominates over the Newtonian potential of the Sun (both approaching zero at
infinity), the potential profile of the galaxy varies slowly enough in such a
way that, for or purposes, it can be considered constant over the solar
neighbourhood. Then, using (92) and (126) with $U=U_{g}$, it follows that
$1-\gamma(r)$ is a decreasing function of $r$ in this region. Such a
qualitative behaviour of $1-\gamma(r)$ is in agreement with the result found
numerically in Ref. HS  for $f(R)$ gravity.

Inserting now the expression (91) for $\Psi$ in formula (127) for $1-\gamma$,
and using the upper bound (125) from the Cassini measurement, we obtain the
following inequality:

$\frac{16\pi}{M_{\odot}}q_{1}\rho_{s}^{\prime}r_{s}^{2}+\frac{2}{3}\,\frac{M_{\rm
eff}(r)}{M_{\odot}}+\frac{8\pi}{3}\,\frac{r}{M_{\odot}}\int_{r}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]r^{\prime}dr^{\prime}<2\times
10^{-6},$ (129)

for $r$ varying in the interplanetary space. The integral is extended over the
solar atmosphere and, using the corresponding values of mass density reported
in Appendix B, it turns out that the integral term is negligible in comparison
with $10^{-6}$ so that, taking into account the approximation (97) of the
effective mass, the previous inequality is approximated by

$\frac{16\pi}{M_{\odot}}q_{1}\rho_{s}^{\prime}r_{s}^{2}+\frac{2}{3}\,\frac{M_{\rm
eff}(R_{\odot})}{M_{\odot}}<2\times 10^{-6}.$ (130)

Since the effective mass $M_{\rm eff}(R_{\odot})$ depends on $r_{s}$, and the
screening radius is determined by the integral equation (III.2.1) which
contains the NMC gravity parameters $\alpha,q_{1},q_{2}$, then inequality
(130) is a constraint on NMC gravity parameters.

We observe that if $q_{1}=0$ inequality (130) holds true also for $f(R)$
gravity (see related results in Ref. HS ).

###  V.1 Computation of the screening radius and of the effective mass

The mass of the Sun’s convection zone is about $0.02M_{\odot}$ Christ  and the
mass of the solar atmosphere is about $10^{-10}M_{\odot}$. If $\varepsilon$ in
inequality (95) is small enough ($\varepsilon\leq 10^{-2}$ is sufficient),
then using inequalities (99) and (130), it follows that the screening radius
that saturates the Cassini constraint, i.e. that makes the inequality (130) an
equality, lies in the convection zone.

Moreover, for such values of $\varepsilon$ inequality (130) is satisfied only
if $R_{\odot}-r_{s}\ll R_{\odot}$, so that the contribution to the effective
mass only comes from a thin shell KW  of radii $(r_{s},R_{\odot})$ in the
upper part of the convection zone and in the photosphere.

In this section we estimate the screening radius in the convection zone, as a
function of parameters $\alpha,q_{1},q_{2}$, by means of an approximate
solution of the integral equation (III.2.1), obtained by resorting to the
solar density model reported in Appendix B.

In order to solve the integral equation (III.2.1) we use the approximation
(98). We denote $r_{p}$ the radius at the base of the photosphere, which
coincides with the top of the convection zone. If $r_{s}$ lies in the
convection zone, since the integration is extended over the thin shell, then
we approximate the integral on the right hand side of (98), considered as a
function of $r_{s}$, by means of the second order Taylor expansion around
$r_{p}$, obtaining

$\frac{1}{6}\int_{r_{s}}^{R_{\odot}}\rho(r)rdr\approx I_{\rm ph}+P(r_{s}),$
(131)

where $I_{\rm ph}=4.19\times 10^{10}\,{\rm g\,cm^{2}}$ is the contribution of
the photosphere computed by using the mass density profile reported in Section
.4.1 of Appendix B, and

$P(r_{s})=\frac{1}{12}\left[\rho_{p}\left(r_{p}^{2}-r_{s}^{2}\right)-\rho^{\prime}_{p}r_{p}\left(r_{p}-r_{s}\right)^{2}\right],$
(132)

where $\rho_{p}=\rho(r_{p})$, $\rho_{p}^{\prime}=d\rho/dr(r_{p})$. The
quantities $\rho_{p}$ and $\rho^{\prime}_{p}$ are computed by using the
density profile of the convection zone reported in Section .3.1 of Appendix B,
and the derivative is computed from the side of the convection zone.

The term $\rho_{s}+\rho_{s}^{\prime}r_{s}$ in the integral equation (III.2.1)
is approximated in the thin shell by means of the quadratic approximation of
density inside the convection zone starting at $r_{p}$:

$\displaystyle\rho_{s}+\rho_{s}^{\prime}r_{s}$ $\displaystyle\approx$
$\displaystyle Q(r_{s})=\rho_{p}-\rho_{p}^{\prime}(r_{p}-2r_{s})$ (133)
$\displaystyle+$
$\displaystyle\frac{1}{2}\rho_{p}^{\prime\prime}(r_{p}-r_{s})(r_{p}-3r_{s}).$

Substituting such second order approximations in the integral equation
(III.2.1), the screening radius $r_{s}$ then solves the following quadratic
equation:

$P(r_{s})+q_{1}Q(r_{s})+I_{\rm ph}=\alpha q_{2}\left(\frac{8\pi
G}{c^{2}}\right)^{\alpha-1}\rho_{g}^{\alpha}.$ (134)

The solution $r_{s}=r_{s}(\alpha,q_{1},q_{2})$ has to be selected out of the
two roots of the quadratic equation and then substituted into inequality
(130), particularly into the expression of the effective mass. The
contribution of the thin shell to the effective mass is computed by means of
an analogous second order Taylor approximation of the integral in (97):

$\displaystyle M_{\rm eff}(R_{\odot})$ $\displaystyle\approx$ $\displaystyle
4\pi
r_{p}\left(r_{p}-r_{s}\right)\left[\rho_{p}r_{s}-\frac{1}{2}\rho^{\prime}_{p}r_{p}\left(r_{p}-r_{s}\right)\right]$
(135) $\displaystyle+$ $\displaystyle M_{\rm ph},$

where $M_{\rm ph}=2.2\times 10^{23}\,{\rm g}$ is the mass of the photosphere
computed by using the mass density profile reported in Section .4.1 of
Appendix B.

Eqs. (130), (134) and (135) will be used in the next section to find the
Cassini constraint on NMC gravity parameters $\alpha,q_{1},q_{2}$.

###  V.2 Constraints on NMC gravity parameters

The constraint from Cassini measurement determines an admissible region in the
three-dimensional parameter space with coordinates $\alpha,q_{1},q_{2}$. We
represent the admissible region by means of two-dimensional exclusion plots
obtained using sections with planes $\alpha=constant$ and $q_{1}=constant$.

In order to avoid either too small or too large numbers we replace parameters
$q_{1},q_{2}$ with the following rescaled, dimensionless parameters:

$\widetilde{q}_{1}=\frac{q_{1}}{R_{\odot}^{2}},\qquad\widetilde{q}_{2}=q_{2}R_{g}^{\alpha},$
(136)

with $R_{g}=8\pi G\rho_{g}/c^{2}$. With this substitution, the quadratic
equation (134) which determines the screening radius $r_{s}$ becomes

$\frac{8\pi
G}{c^{2}}\left[P(r_{s})+\widetilde{q}_{1}R_{\odot}^{2}Q(r_{s})+I_{\rm
ph}\right]=\alpha\widetilde{q}_{2}.$ (137)

The admissible root, which has to satisfy $r_{s}\leq r_{p}$, expressed in the
form $r_{s}=r_{s}(\alpha,\widetilde{q}_{1},\widetilde{q}_{2})$, is substituted
in the approximation (135) of the effective mass, then the admissible values
of parameters $\alpha,\widetilde{q}_{1},\widetilde{q}_{2}$ satisfy inequality
(130). The admissible region in parameter space is restricted by means of the
condition $\lambda_{g}\gg r_{g}$ introduced in Section III.1.2 which,
expressed in terms of parameters $\alpha,\widetilde{q}_{2}$ and using Eq. (37)
with $\rho=\rho_{g}$, becomes

$\left[\frac{3}{4\pi}\alpha(1-\alpha)\frac{c^{2}}{G\rho_{g}}\widetilde{q}_{2}\right]^{1/2}>10^{2}r_{g},$
(138)

where we have required $\lambda_{g}>10^{2}r_{g}$. The admissible region is
further restricted by the consistency condition, inequality (95):

$|q_{1}|<\varepsilon\,\frac{\rho(r)}{|\nabla^{2}\rho(r)|}\quad\mbox{for any
}r<r_{s}.$ (139)

Using the expression of $\nabla^{2}\rho$ in Section .3.1 of Appendix B, the
ratio $\rho/|\nabla^{2}\rho|$ is a decreasing function of $r$ in the
convection zone, so that, for given $\varepsilon$, the upper bound on
$|q_{1}|$ is smaller for $r=r_{s}$. Moreover, for given values of $q_{1}$ and
$\varepsilon$ the consistency condition yields an upper bound
$r_{s}<r_{s}^{\ast}=r_{s}^{\ast}(q_{1},\varepsilon)$ on the screening radius,
so that, using the integral equation (III.2.1), the approximation (98), and
inequality (100), we get the following lower bound on
$\alpha\widetilde{q}_{2}$ for $\varepsilon\leq 1/20$:

$\alpha\widetilde{q}_{2}>\frac{4\pi
G}{3c^{2}}\left[(1-20\varepsilon)\int_{r_{s}^{\ast}}^{R_{\odot}}\rho(r)rdr\right].$
(140)

Our results are graphically reported in Figures 2-7: admissible regions for
parameters are plotted in white, while the excluded regions are plotted in
grey. Fig. 2 shows the section of the admissible region in three-dimensional
parameter space with the plane $\widetilde{q}_{1}=0$; since both $\alpha$ and
$\widetilde{q}_{2}$ are negative the admissible region is plotted in the
quarter plane with coordinates $(|\alpha|,|\widetilde{q}_{2}|)$. Fig. 2 shows
the sections with the plane $\widetilde{q}_{1}=10^{-7}$ corresponding to
values $\varepsilon=10^{-2}$ and $\varepsilon=10^{-3}$. Using inequality (139)
the value $\widetilde{q}_{1}=10^{-7}$ corresponds to the upper bound
$r_{s}^{\ast}\approx 0.99R_{\odot}$ for $\varepsilon=10^{-2}$. Fig. 2 shows
that the lower bound (140) is greater than the one resulting from inequality
$\lambda_{g}>10^{2}r_{g}$.

Figure 1: Cassini constraint on the parameter quarter plane
$|\alpha|,|\widetilde{q}_{2}|$ for $\widetilde{q}_{1}=0$. The solid line
yields the upper bound on $|\widetilde{q}_{2}|$ from Cassini measurement, the
dotted line yields the lower bound on $|\widetilde{q}_{2}|$ from inequality
$\lambda_{g}>10^{2}r_{g}$.

Figure 2: Cassini constraint on the parameter quarter plane
$|\alpha|,|\widetilde{q}_{2}|$ for $\widetilde{q}_{1}=10^{-7}$. The dotted
colored lines yield the lower bound on $|\widetilde{q}_{2}|$ from inequality
(140). If $\varepsilon=10^{-2}$ the zone plotted in light grey between the
dotted colored lines is admissible.

Figures 4 and 4 show the sections of the admissible region in three-
dimensional parameter space with the planes $\alpha=-1$ and $\alpha=-10$,
respectively. The admissible region is plotted in the half plane with
coordinates $(\widetilde{q}_{1},|\widetilde{q}_{2}|)$. In the case
$\widetilde{q}_{1}=0$, the value $\alpha=-1$ has been used in Ref. drkmattgal
to model the rotation curves of galaxies, and the value $\alpha=-10$ has been
used in Ref. curraccel  to model the current accelerated expansion of the
Universe. The asymmetry of the admissible region with respect to the axis
$\widetilde{q}_{1}=0$ is due to the impact of the sign of $q_{1}$ in the
solution of the integral equation (III.2.1).

Figure 3: Cassini constraint on the parameter half plane
$\widetilde{q}_{1},|\widetilde{q}_{2}|$ for $\alpha=-1$. The solid line is the
Cassini bound, the dotted black line is the bound from inequality
$\lambda_{g}>10^{2}r_{g}$, the dotted colored lines are the bounds from
inequality (140). If $\varepsilon=10^{-2}$ the zone plotted in light grey
between the dotted colored lines and the Cassini bound is admissible.

Figure 4: Cassini constraint on the parameter half plane
$\widetilde{q}_{1},|\widetilde{q}_{2}|$ for $\alpha=-10$. The solid line is
the Cassini bound, the dotted black line is the bound from inequality
$\lambda_{g}>10^{2}r_{g}$, the dotted colored lines are the bounds from
inequality (140). If $\varepsilon=10^{-2}$ the zone plotted in light grey
between the dotted colored lines and the Cassini bound is admissible.

Figure 5: Allowed values of $\widetilde{q}_{1}$ as function of the
dimensionless screening radius $\xi=r_{s}/R_{\odot}$, for different values of
the parameter $\varepsilon$ from the consistency condition, Eq. (95). The
excluded regions are colored in grey.

Figure 5 shows the allowed region for $\widetilde{q}_{1}$ as function of the
dimensionless screening radius $\xi=r_{s}/R_{\odot}$. That is computed by
intersecting the region allowed by the inequality of the Cassini constraint,
Eq. (130), with the region allowed by the consistency condition, Eq. (95), for
the values $\varepsilon=10^{-2}$ and $\varepsilon=10^{-3}$. Since
$\rho_{s}^{\prime}<0$ the Cassini constraint yields a lower bound on
$\widetilde{q}_{1}$ while the consistency condition yields both an upper
(positive) bound on $\widetilde{q}_{1}$ and a lower (negative) bound. The
intersection point between the curve of the Cassini constraint and the curve
of the upper bound from the consistency condition, defines a lower bound on
the screening radius. Such a lower bound $\xi_{\rm low}$ is the value of the
screening radius that saturates the Cassini constraint: values of $r_{s}$
below $\xi_{\rm low}R_{\odot}$ are excluded. Figure 5 shows that the lower
bound on $\xi$ weakly depends on $\varepsilon$ and its value is $\xi_{\rm
low}\approx 0.98$ which corresponds to a maximum radial thickness of the thin
shell of about $0.02R_{\odot}$.

Figure 6: Allowed values of $|\widetilde{q}_{2}|$ for $\alpha=-1$ as function
of the dimensionless screening radius $\xi=r_{s}/R_{\odot}$ for different
values of the parameter $\varepsilon$ from the consistency condition, Eq.
(95). The excluded regions are colored in grey.

Figure 7: Allowed values of $|\widetilde{q}_{2}|$ for $\alpha=-10$ as function
of the dimensionless screening radius $\xi=r_{s}/R_{\odot}$ for different
values of the parameter $\varepsilon$ from the consistency condition, Eq.
(95). The excluded regions are colored in grey.

Figures 7 and 7 show the allowed region for $|\widetilde{q}_{2}|$ as function
of $\xi$ for a fixed value of $\alpha$. That is computed by solving the
integral equation (III.2.1) with respect to $q_{2}$, then allowing $q_{1}$ to
vary in its admissible region. This has been done for $\alpha=-1$ in Fig. (7),
and $\alpha=-10$ in Fig. 7; in both cases the values $\varepsilon=10^{-2}$ and
$\varepsilon=10^{-3}$ have been used. The constant $\widetilde{C}_{\alpha}$
corresponds to the lower bound on $\widetilde{q}_{2}$ from inequality (138),
and both figures show that its value is much smaller than the lower bound
computed by means of the integral equation. Eventually, we observe that
$|\widetilde{q}_{2}|$ decreases monotonically as the screening radius
increases.

###  V.3 Verification of the consistency condition

Once the constraint from Cassini measurement has been satisfied, we have to
check that the consistency condition (38) is verified in the solar interior
for $r\leq r_{s}$. We verify such a condition in the convection zone: using
the model of density profile $\rho(r)$ in this zone, reported in Section .3.1
of Appendix B, we have

$\nabla^{2}\rho=\frac{n_{c}-1}{n_{c}\rho}\left(\frac{d\rho}{dr}\right)^{2},\qquad
n_{c}=2.33,$ (141)

from which, using Eq. (50) for $\lambda^{2}$ and taking into account that
$\nabla^{2}\rho(r)>0$ in the convection zone, the consistency condition
becomes

$\frac{\nabla^{2}\rho}{\rho}\left|\lambda_{g}^{2}\left(\frac{\rho}{\rho_{g}}\right)^{\alpha-1}\left[\frac{\alpha}{1-\alpha}-\frac{n_{c}}{n_{c}-1}(\alpha+1)\right]+6q_{1}\right|\ll
1.$ (142)

An expression for $\lambda_{g}$ can be found from the general expression (37)
of $\lambda$ and the integral equation (III.2.1). Let us choose
$|q_{1}|<\varepsilon 10^{-7}R_{\odot}^{2}$, with $\varepsilon<10^{-2}$ (see
Section III.1.3). Then, using inequality (100), for an estimate of the order
of magnitude of $\lambda_{g}$ we can neglect the term involving $q_{1}$ in the
integral equation, obtaining:

$\lambda_{g}^{2}\approx\frac{1-\alpha}{\rho_{g}}\int_{r_{s}}^{R_{\odot}}\rho(r)rdr,$
(143)

where we have neglected the contribution to the integral from density in the
solar atmosphere outside the photosphere. Now, it turns out that both

$\frac{\nabla^{2}\rho}{\rho}\qquad\mbox{and}\qquad\left(\frac{\rho}{\rho_{g}}\right)^{\alpha-1}$
(144)

are increasing functions of $r$ in the convection zone, and the maximum of
$\nabla^{2}\rho/\rho$ is of order of $10^{7}/R_{\odot}^{2}$. Hence, if the
consistency condition is satisfied for $r=r_{s}$, then it is satisfied also
for any $r<r_{s}$ and the verification of such a condition depends mainly on
the factor $(\rho/\rho_{g})^{\alpha-1}$.

Let us now choose for instance the value $r_{s}\approx 0.98R_{\odot}$ for the
screening radius, which saturates the Cassini constraint (see the results in
the previous section). Using the density profile in Section .3.1 of Appendix
B, we have

$\frac{\nabla^{2}\rho(r)}{\rho(r)}<\frac{10^{7}}{R_{\odot}^{2}},\quad\left(\frac{\rho(r)}{\rho_{g}}\right)^{\alpha-1}<10^{-20(1-\alpha)},\quad\mbox{for
}r<r_{s},$ (145)

and using (143) we have
$\lambda_{g}^{2}\approx(1-\alpha)10^{18}R_{\odot}^{2}$. Inserting these
numbers into inequality (142) it turns out that the consistency condition is
verified for the values $\alpha\leq-1$ which satisfy the Cassini constraint.

Analogous results can be found for different values of $r_{s}$ (satisfying the
Cassini constraint) and for $r$ varying in the radiative interior.

##  VI Effect of the extra force in the solar atmosphere

We compute the effect of fifth force and extra force on temperature and
density in the solar atmosphere when the screening radius lies in the
convection zone. We consider the chromosphere-corona transition region where a
steep density gradient takes place Mariska  (see Section .4.3 of Appendix B),
and the extra force is expected to have largest intensity.

Models of the transition region have been derived from the distribution of the
emission measure computed from the intensities of spectral lines in the
ultraviolet region of the solar spectrum Jordan76 . In the approximation of a
plane parallel geometry the average emission measure distribution $EM$, for
the quiet Sun, is parametrized as a function of temperature $T$ by Jordan

$EM=\int_{\Delta T}N_{e}^{2}dh=aT^{b},$ (146)

where $N_{e}$ is the electron number density, $h=h(T)$ is height above the
solar limb, $\Delta T$ denotes the temperature range $\Delta
T=(0.891T,1.122T)$ used in RayDoyle , and $a,b$ are parameters which take
different values in the two zones below and above the temperature
$T_{0}=1.5\times 10^{5}\,{\rm K}$, respectively.

In the sequel we use the following approximations: we assume hydrostatic
equilibrium, complete ionization of hydrogen (which is realized for $T\geq
2\times 10^{4}$ K), single ionization of helium in the zone below $T_{0}$
(lower transition region), and double ionization of helium in the zone above
$T_{0}$ (upper transition region). Moreover, we neglect elements other than
hydrogen and helium both in the equation of hydrostatic equilibrium and in
ionization equilibrium. Assuming the perfect gas law, the electron density is
then given by $N_{e}=p/\xi k_{B}T$, where $p$ is pressure, $k_{B}$ is the
Boltzmann constant, and $\xi$ is related to the degree of ionization of
helium: $\xi=2$ in the lower transition region, and $\xi=1.91$ in the upper
transition region. Equality (146) is then written in the form

$\frac{1}{\xi^{2}k_{B}^{2}}\int_{\Delta
T}\frac{p^{2}}{(T^{\prime})^{2}}\frac{dh}{dT^{\prime}}dT^{\prime}=aT^{b}.$
(147)

Following the method in Jordan76 ; Jordan , both pressure $p$ and the quantity
$(1/T^{2})dh/dT$ are approximated with a constant over the temperature
interval $\Delta T$, so that Eq. (147) yields for a monotone temperature
profile $T(h)$:

$\frac{p^{2}}{T}\approx\frac{\xi^{2}}{0.231}k_{B}^{2}\frac{a}{b+1}\frac{d}{dh}\left(T^{b+1}\right),$
(148)

where, both in this formula and in the following computations, the two zones
in the transition region below and above the temperature $T_{0}$, have to be
considered separately (see Section .4.3 of Appendix B).

Multiplying now by pressure $p$ both sides of Eq. (110) in the case of
hydrostatic equilibrium (i.e. $dv/dt=0$), approximating $dU/dr\approx-
GM_{\odot}/R_{\odot}^{2}$, replacing $d/dr$ with $d/dh$, using $\rho=\mu
m_{p}p/k_{B}T$ with $m_{p}$ the proton mass and $\mu$ the mean molecular
weight, and using expressions (115) with $r^{2}\approx R_{\odot}^{2}$ and
(121) of fifth force and extra force outside the screening radius,
respectively, we obtain

$\displaystyle\frac{d}{dh}\left(\frac{p^{2}}{2}\right)$ $\displaystyle=$
$\displaystyle-\frac{\mu
m_{p}G}{k_{B}R_{\odot}^{2}}\left(M_{\odot}+\frac{1}{3}M_{\rm
eff}(R_{\odot})+8\pi q_{1}\rho^{\prime}_{s}r_{s}^{2}\right)\frac{p^{2}}{T}$
(149) $\displaystyle+$
$\displaystyle\mathcal{F}(\alpha,q_{1},q_{2})\,p\rho^{\frac{\alpha}{1-\alpha}}\,\frac{d\rho}{dh}.$

We use $\mu=0.65$ in the lower transition region and $\mu=0.62$ in the upper
transition region. We denote $p_{0}$ the pressure at the temperature $T_{0}$,
then setting $p=p_{0}+\delta p$, it is known that $\delta p$ is a small
variation in the transition zone Doschek , and we approximate

$\mathcal{F}(\alpha,q_{1},q_{2})\,p\approx\mathcal{F}(\alpha,q_{1},q_{2})\,p_{0},$
(150)

where we have negected the product $\mathcal{F}\delta p$. Then, combining Eqs.
(148) and (149), and expressing again $\rho$ in terms of pressure and
temperature, we find that the following expression has a vanishing derivative
with respect to $h$:

$\frac{p^{2}}{2}+\frac{\xi^{2}}{0.231}\frac{\mu
m_{p}k_{B}a}{b+1}\frac{G}{R_{\odot}^{2}}\left(M_{\odot}+\frac{1}{3}M_{\rm
eff}(R_{\odot})+8\pi
q_{1}\rho^{\prime}_{s}r_{s}^{2}\right)T^{b+1}-(1-\alpha)\mathcal{F}(\alpha,q_{1},q_{2})\left(\frac{\mu
m_{p}}{k_{B}}\right)^{\frac{1}{1-\alpha}}p_{0}^{\frac{2-\alpha}{1-\alpha}}T^{\frac{1}{\alpha-1}}=constant,$
(151)

where we have again approximated
$\mathcal{F}p^{1/1-\alpha}\approx\mathcal{F}p_{0}^{1/1-\alpha}$. The value of
the constant in Eq. (151) is given by the above expression evaluated at
$T_{0},p_{0}$, so that pressure $p$ and electron number density $N_{e}$ can be
computed, obtaining

$N_{e}(\alpha,q_{1},q_{2})=\frac{1}{\xi
T}\left(A_{1}+A_{2}T^{b+1}+A_{3}T^{\frac{1}{\alpha-1}}\right)^{1/2},$ (152)

where we stress the dependence of $N_{e}$ on NMC parameters, and

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\left(\frac{p_{0}}{k_{B}}\right)^{2}-A_{2}T_{0}^{b+1}-A_{3}T_{0}^{\frac{1}{\alpha-1}},$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle-\frac{8.658\xi^{2}}{b+1}a\frac{\mu
m_{p}G}{k_{B}R_{\odot}^{2}}\left(M_{\odot}+\frac{1}{3}M_{\rm
eff}(R_{\odot})+8\pi q_{1}\rho^{\prime}_{s}r_{s}^{2}\right),$ $\displaystyle
A_{3}$ $\displaystyle=$
$\displaystyle(1-\alpha)\frac{2}{k_{B}^{2}}\left(\frac{\mu
m_{p}}{k_{B}}\right)^{\frac{1}{1-\alpha}}p_{0}^{\frac{2-\alpha}{1-\alpha}}\mathcal{F}(\alpha,q_{1},q_{2}).$
(153)

For pressure $p_{0}$ we use $p_{0}/k_{B}=2\times 10^{15}\,{\rm K\,cm^{-3}}$
derived from RayDoyle , while the values of the various constants are reported
in Appendix B.

Accurate determination of electron density in the transition region is
achieved by the observation of the intensity ratio of spectral lines emitted
by the same ion at the temperature of formation of such an ion Doschek .
Particularly, the quiet Sun intensity ratios of lines emitted by the ${\rm
Si}^{+2}$ ion have been used to determine the electron density $N_{e}$ in the
lower transition region Doschek , where $T<T_{0}$, density exhibits the
steepest gradient (see Section .4.3 of Appendix B), and the extra force has
largest intensity.

Uncertainties in the derived electron densities Doschek  are due to
uncertainties in line ratios (typically less than $20\%$ for quiet Sun),
spectral line blending ($25\%$-$30\%$ as maximum possible effect), time
variations of the physical conditions in the solar atmosphere ($15\%$ change
in line intensity for quiet Sun over a time interval of 1.5 hr), and
uncertainties in the atomic physics data (systematic): an overall uncertainty
of a factor of 2 in density is then considered in Doschek .

Since the measurement of electron density from spectral line ratios does not
depend on gravity, while the expression (152) of $N_{e}$ is affected by
gravity, then equating such an expression with density $N_{e,{\rm obs}}$,
observed from line ratios, yields a constraint on $\alpha,q_{1},q_{2}$.

The electron density at the temperature of formation of the ${\rm Si}^{+2}$
ion in ionization equilibrium, determined from intensity ratios of lines
within the Si III multiplet near 1300 Å, is found to be for quiet Sun Doschek

$N_{e,{\rm obs}}=4.3\times 10^{10}\,{\rm cm}^{-3}.$ (154)

Since $\mathcal{F}(\alpha,q_{1},q_{2})<0$, taking into account the overall
density uncertainty considered above, we get the constraint

$N_{e}(\alpha,q_{1},q_{2})\geq 2.15\times 10^{10}\,{\rm cm}^{-3},$ (155)

where formula (152) has to be evaluated using the temperature $T=3.2\times
10^{4}$ K of formation of the ${\rm Si}^{+2}$ ion, and using the values of
$a,b$ reported in Section .4.3 of Appendix B for the lower transition region.

We observe that, using the constraint (130) from Cassini measurement, the
coefficient $A_{2}$ can be approximated by the Newtonian value

$A_{2}\approx-\frac{8.658\xi^{2}}{b+1}a\frac{\mu
m_{p}}{k_{B}}\,\frac{GM_{\odot}}{R_{\odot}^{2}},$ (156)

and, using the expression (115) of the fifth force outside the screening
radius, it follows

$\left|F_{\rm f}(r)\right|<10^{-6}\rho(r)\frac{GM_{\odot}}{R_{\odot}^{2}},$
(157)

so that the fifth force is negligible in comparison with the Newtonian force.

###  VI.1 Constraints on NMC gravity parameters

The constraint from spectroscopic measurements determines an admissible region
in the three-dimensional parameter space with coordinates
$\alpha,q_{1},q_{2}$. We represent the admissible region by means of two-
dimensional exclusion plots obtained using sections with planes
$\alpha=constant$. We use the variables
$\xi,\widetilde{q}_{1},\widetilde{q}_{2}$ introduced in Section V.2 about the
Cassini constraint.

We find convenient to express the electron density $N_{e}$ as a function of
$\alpha$, $q_{1}$ and the screening radius $r_{s}$, moreover, we use
approximation (156). The screening radius appears in coefficient $A_{3}$,
particularly in the term $\eta(R_{\odot})$ inside $\mathcal{F}$, both
explicitly and implicitly through the effective mass $M_{\rm eff}(R_{\odot})$.
In order to obtain the function $N_{e}=N_{e}(\alpha,q_{1},r_{s})$ we eliminate
$q_{2}$ from $\mathcal{F}(\alpha,q_{1},q_{2})$ using the integral equation
(III.2.1):

$\alpha q_{2}\approx\frac{1}{\rho_{g}^{\alpha}}\left(\frac{8\pi
G}{c^{2}}\right)^{1-\alpha}\left[q_{1}\left(\rho_{s}+\rho_{s}^{\prime}r_{s}\right)+\frac{1}{6}\int_{r_{s}}^{R_{\odot}}\rho(r)rdr\right],$
(158)

where we have again neglected the contribution to the integral from density in
the solar atmosphere outside the photosphere. Then we substitute the above
expression of $\alpha q_{2}$ in formula (122) of $\mathcal{F}$ and in formula
(101) of $\eta$, which has to be evaluated in $r=R_{\odot}$ to compute
$\mathcal{F}$. After these substitutions the constraint (155) becomes

$N_{e}(\alpha,q_{1},r_{s})\geq 2.15\times 10^{10}\,{\rm cm}^{-3}.$ (159)

For given $\alpha$ such an inequality determines an admissible region in the
half plane with coordinates $(q_{1},r_{s})$. Such a region is further
restricted by the consistency condition, inequality (95), that has to be
satisfied for any $r<r_{s}$. Hence, using coordinates
$(\widetilde{q}_{1},\xi)$ in the half plane, the admissible region is
restricted by means of the intersection with a region, symmetric with respect
to the axis $\widetilde{q}_{1}=0$, which shrinks as $\xi$ increases, being
$\rho/\nabla^{2}\rho$ a decreasing function of $\xi$.

Figure 8: Constraint from extra force on the half plane
$\widetilde{q}_{1},\xi$ for $\alpha=-1/4$. The solid line yields the lower
bound on the screening radius from spectroscopic measurements. The dotted
colored lines yield the bound from the consistency condition: if
$\varepsilon=10^{-2}$ the zone plotted in light grey between the dotted
colored lines and the bound on $\xi$ is admissible.

Figure 9: Constraint from extra force on the half plane
$\widetilde{q}_{1},\xi$ for $\alpha=-1/3$. The lower bound on the screening
radius increases, but it decreases if $\alpha$ approaches the value
$\alpha=-1$. The dotted colored lines yield the bound from the consistency
condition: if $\varepsilon=10^{-2}$ the zone plotted in light grey between the
dotted colored lines and the bound on $\xi$ is admissible.

Figure 10: Constraint from extra force on the half plane
$\widetilde{q}_{1},|\widetilde{q}_{2}|$ for $\alpha=-1/4$. The solid line
yields the upper bound on $|\widetilde{q}_{2}|$ from spectroscopic
measurements. The dotted colored lines yield the bound from the consistency
condition: if $\varepsilon=10^{-2}$ the zone plotted in light grey between the
dotted colored lines and the bound on $|\widetilde{q}_{2}|$ is admissible.

Figure 11: Constraint from extra force on the half plane
$\widetilde{q}_{1},|\widetilde{q}_{2}|$ for $\alpha=-1/3$. The upper bound on
$|\widetilde{q}_{2}|$ decreases, but it increases if $\alpha$ approaches the
value $\alpha=-1$. The dotted colored lines yield the bound from the
consistency condition: if $\varepsilon=10^{-2}$ the zone plotted in light grey
between the dotted colored lines and the bound on $|\widetilde{q}_{2}|$ is
admissible.

The numerical computation shows that the constraint from spectroscopic
measurements imposes a lower bound on the screening radius which lies in the
convection zone for $-1<\alpha<0$. Conversely, for $\alpha\leq-1$ the
screening radius lies in Sun’s radiative interior, so that for such values of
$\alpha$ the constraint from the extra force is surely much weaker than the
constraint from Cassini measurement.

Our results are graphically reported in Figures 9-11: admissible regions for
parameters are plotted in white, while the excluded regions are plotted in
grey. We denote $r_{\rm conv}$ the radius at the base of the convection zone.
Fig. 9 shows the admissible region in the half plane $(\widetilde{q}_{1},\xi)$
for $\alpha=-1/4$ and for values $\varepsilon=10^{-2}$ and
$\varepsilon=10^{-3}$, while Fig. 9 shows the admissible region for
$\alpha=-1/3$.

Figures 11 and 11 show the admissible region in the half plane with
coordinates $(\widetilde{q}_{1},|\widetilde{q}_{2}|)$ for $\alpha=-1/4$ and
$\alpha=-1/3$, respectively. The value $\alpha=-1/3$ has been used in Ref.
drkmattgal  to model the rotation curves of galaxies. Figure 11 shows some
numerical noise on the left side of the solid curve which marks the upper
bound on $|\widetilde{q}_{2}|$ from the extra force, nevertheless such a noise
does not affect importantly the behavior of the bound, moreover it lies in the
excluded region.

When the screening radius lies deeply in the convection zone the Taylor
approximation of the effective mass and of the integral in the equation
(III.2.1) cannot be used, then the integrals are evaluated analytically by
approximating the polytropic index in the convection zone with $n_{c}\approx
2$ (see Section .3.1 of Appendix B).

Eventually, we argue that the constraint from extra force becomes competitive
with the constraint from Cassini measurement for suitable values of $\alpha$
in the range $-1<\alpha<0$. Further improvements could be expected from new
spectroscopic measurements in solar missions like ESA Solar Orbiter Muller-Cyr
and NASA Parker Solar Probe Fox-Velli .

##  VII Conclusions

In this work we obtain Solar System bounds on an alternative theory of gravity
with nonminimal coupling between curvature and matter, that has been
introduced in Ref. BBHL , by suitably implementing a screening mechanism which
is the NMC version of the well known chameleon mechanism KW ; HS . An analytic
approximate solution of the field equations that exhibits the screening
mechanism has been computed inside and around the Sun. Such a solution shows
the typical features of a chameleon solution, though significantly affected by
the nonminimal coupling: i) a screening radius $r_{s}$ in the Sun’s interior
that has to be determined, and such that GR is approximately satisfied for
$r<r_{s}$ if a suitable consistency condition on the parameters of the model
is also satisfied; ii) a thin shell close to the Sun’s edge where the solution
interpolates between the values in the Sun’s interior and the values in the
outskirts of the Solar System; iii) the existence of the thin shell permits to
satisfy stringent constraints from Solar System gravity experiments.

A specific feature of the nonminimal coupling between geometry and matter is
the presence of a non-Newtonian extra force which, together with a fifth force
of the Yukawa type, appears in the solar interior and atmosphere modeled as a
perfect fluid. Constraints on the parameters of the gravity model have been
computed from the Cassini measurement of PPN parameter $\gamma$, and from the
effect in the solar atmosphere of the extra force by resorting to
spectroscopic measurements.

Our results show that taking into account all known solar physics bounds, the
Cassini measurement of the parameter $\gamma$ constrains the parameters for
the specific model Eq. (13), under conditions Eq. (31) and rescaled according
to Eq. (136), to be order $|\widetilde{q}_{2}|<10^{-12}$ for
$|\widetilde{q}_{1}|<3\times 10^{-6}$ and $-10<\alpha<-1$. These are specific
constraints obtained under the assumption that the galaxy itself is screened
within a distance from its center corresponding to the location of the Solar
System. The assessment of such an assumption requires an extension of the
Solar System to galaxy analysis carried out in the present paper.
Particularly, one should look for a chameleon solution of the field equations
in the Milky Way considering the transition from the outskirts of the galactic
halo to the large-scale structure of the Universe (see Ref. HS  for the
analogous problem in $f(R)$ gravity). Such an analysis will impose additional
constraints on the parameters of the NMC gravity model that have to be
compared with the Solar System bounds obtained in the present paper. Moreover,
it could also help to better understand the role of different powers $\alpha$
of curvature in the specific model Eq. (13) at different scales, as it has
been conjectured in Ref. drkmattgal . Eventually, the extension of the
analysis in the present paper to the gravitational interaction between the
Sun, planets and satellites should also provide further constraints on the
parameters of the gravity model. All these issues will be investigated in
future research.

## Appendix A

We give a proof a posteriori of inequalities Eqs. (59) and (64) that have been
used in order to compute the solution $\eta$ of equation (26). Hence, we prove
that the solution (75) found for function $\eta$ is consistent with such
inequalities.

###  .1 Verification of inequality $|\partial V/\partial\eta|\ll\rho c^{2}/3$

The inequality is equivalent to inequality (60) which we then consider. We
prove that for $r>r_{s}$ and $r$ close to $r_{s}$, the solution (75) for
$\eta$ is such that curvature $R=\omega(\eta,\rho)$ quickly becomes much
smaller than the GR curvature, so that inequality (60) is verified a
posteriori.

We use formula (62) for curvature which we write in the form

$\frac{c^{2}}{8\pi
G}\frac{\omega(\eta,\rho)}{\rho}=\left(\frac{\rho}{\rho_{s}}\right)^{\alpha/(1-\alpha)}\left(\frac{N}{D}\right)^{1/(1-\alpha)}.$
(160)

Using Eqs. (37) and (62) we have

$N=\frac{\lambda_{s}^{2}}{6(1-\alpha)}\rho_{s},\qquad D=\frac{c^{2}}{16\pi
G}-q_{1}\rho-\frac{1}{2c^{2}}\eta.$ (161)

We approximate the denominator $D$ by means of a second order Taylor expansion
around $r_{s}$, so that, using the expression (75) for $\eta$, for
$\rho(r)\gg\rho_{g}$ we have

$\displaystyle-\frac{1}{24\pi}\frac{M_{\rm
eff}(r)}{r}+\frac{1}{6}\int_{r_{s}}^{r}\rho(r^{\prime})r^{\prime}dr^{\prime}$
(162) $\displaystyle=$
$\displaystyle\frac{1}{12}r^{2}\rho_{s}\left(1-\frac{r_{s}}{r}\right)^{2}+O\left(\left(1-\frac{r_{s}}{r}\right)^{3}\right),$

and

$D\approx\frac{\lambda_{s}^{2}}{6(1-\alpha)}\left[\rho_{s}+\alpha\rho_{s}^{\prime}r_{s}\left(1-\frac{r_{s}}{r}\right)\right]+A_{s}\left(1-\frac{r_{s}}{r}\right)^{2},$
(163)

where

$A_{s}=\frac{r_{s}^{2}}{2}\left[\frac{1}{6}\rho_{s}-q_{1}\nabla^{2}\rho(r_{s})\right].$
(164)

Since density $\rho(r)$ is decreasing, we have $\alpha\rho_{s}^{\prime}>0$ for
$\alpha<0$ from which, using the consistency condition (38) for small
$\lambda_{s}$ in the form (95), a lower bound for $D$ follows for
$\varepsilon<1/6$ and $r>r_{s}$:

$D>\frac{1}{12}\rho_{s}r_{s}^{2}(1-6\varepsilon)\left(1-\frac{r_{s}}{r}\right)^{2}.$
(165)

Then, using this lower bound and Eq. (50), we obtain an upper bound on the
ratio $\omega/\rho$:

$\frac{c^{2}}{8\pi
G}\frac{\omega(\eta,\rho)}{\rho}<\frac{\rho_{g}}{\rho_{s}}\left(\frac{\rho}{\rho_{s}}\right)^{\alpha/(1-\alpha)}\left[\frac{2}{(1-\alpha)(1-6\varepsilon)}\left(\frac{\lambda_{g}}{r_{s}}\right)^{2}\right]^{1/(1-\alpha)}\left(1-\frac{r_{s}}{r}\right)^{2/(\alpha-1)}.$
(166)

Now we observe that, for given $\alpha$, $q_{1}$ and $r_{s}$, the Yukawa range
$\lambda_{g}$ is determined by formula (37) and by the integral equation
(III.2.1):

$\displaystyle\frac{1}{6}\int_{r_{s}}^{r_{g}}\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right]rdr$
$\displaystyle=$ $\displaystyle-
q_{1}(\rho_{s}+\rho_{s}^{\prime}r_{s})+\frac{\lambda_{g}^{2}}{6(1-\alpha)}\rho_{g}.$

If we compute $\lambda_{g}$ from Eq. (.1), using approximation (98) and
inequality (100), the effect of the term with $q_{1}$ is a fraction of the
term with $\lambda_{g}$ for $\varepsilon<10^{-2}$ and starts to become
negligible for $\varepsilon<10^{-3}$. Since we are interested in order of
magnitude estimates, we neglect the term with $q_{1}$ in the following
computations. An analogous result holds if $r_{s}$ lies in the solar
atmosphere.

Let us give two examples that show for which values of $r$ the desired
inequality (60) is achieved for given $\alpha$ and $r_{s}$ in the convection
zone. By these examples one can understand that the results obtained by using
inequality (166) are quite general.

Let us consider the value $r_{s}\approx 0.98R_{\odot}$ for the screening
radius which, using the model of density profile in the convection zone
reported in Appendix B, corresponds to $\rho_{s}\approx 3\times 10^{-4}\,{\rm
g\,cm^{-3}}$, and makes the constraint from Cassini measurement satisfied for
$-10\leq\alpha\leq-1$ and $\varepsilon\leq 10^{-2}$.

We fix $\varepsilon=10^{-2}$ since the result does not depend significantly on
$\varepsilon$ if $\varepsilon$ is small enough (see inequality (166)), and we
first set $\alpha=-1$ which implies, using Eq. (.1), $\lambda_{g}\approx
1.87\times 10^{9}R_{\odot}\approx 8.7\times 10^{6}$ AU. In this case
inequality (166) becomes

$\displaystyle\frac{c^{2}}{8\pi G}\frac{\omega(\eta,\rho)}{\rho}$
$\displaystyle<$
$\displaystyle\frac{\rho_{g}}{\sqrt{(1-6\varepsilon)\rho_{s}\rho}}\,\frac{\lambda_{g}}{r_{s}}\left(1-\frac{r_{s}}{r}\right)^{-1}$
(168) $\displaystyle<$ $\displaystyle\frac{2\times 10^{-10}}{1-r_{s}/r}\,,$

where for the last upper bound we have used $\rho(r)>\rho(r_{p})\approx
2.73\times 10^{-7}\,{\rm g\,cm^{-3}}$, where $r_{p}$ is the radius at the top
of the convection zone. Inequality (60) then follows:

$\frac{c^{2}}{8\pi G}\omega(\eta,\rho)<10^{-3}\rho\quad\mbox{for
}r-r_{s}>2\times 10^{-7}r_{s}\approx 140\,{\rm m},$ (169)

with $10^{-7}r_{s}\ll R_{\odot}-r_{s}\approx 13900\,{\rm km}$. Hence, for
$r>r_{s}$ and $r$ very close to $r_{s}$ there is a sharp transition to a low
value of curvature (much smaller than GR curvature) which is a typical
property of a chameleon solution (see also Refs. KW ; HS ).

If we now set $\alpha=-2$ we get $\lambda_{g}\approx 2.29\times
10^{9}R_{\odot}\approx 1.1\times 10^{7}$ AU, and

$\displaystyle\frac{c^{2}}{8\pi G}\frac{\omega(\eta,\rho)}{\rho}$
$\displaystyle<$
$\displaystyle\rho_{g}\left[\frac{2}{3(1-6\varepsilon)\rho_{s}}\right]^{1/3}\left[\frac{\lambda_{g}}{r_{s}\rho\left(1-r_{s}/r\right)}\right]^{2/3}$
(170) $\displaystyle<$
$\displaystyle\frac{10^{-12}}{\left(1-r_{s}/r\right)^{2/3}},$

from which inequality (60) follows:

$\frac{c^{2}}{8\pi G}\omega(\eta,\rho)<10^{-6}\rho\quad\mbox{for
}r-r_{s}>10^{-9}r_{s}\approx 70\,{\rm cm}.$ (171)

Increasing $|\alpha|$ the transition to low curvature becomes sharper.

###  .2 Verification of inequality $|q_{1}|\ll|q_{2}|R^{\alpha-1}$

Using Eq. (62) with $R=\omega(\eta,\rho)$, the inequality is written in the
form

$\frac{16\pi G}{c^{2}}|\alpha q_{1}|\rho\ll 1-\frac{8\pi
G}{c^{4}}\eta-\frac{16\pi G}{c^{2}}q_{1}\rho.$ (172)

Using now expression (75) for $\eta$, neglecting $\lambda_{s}$ and assuming
the screening radius in the solar convection zone, we approximate $M_{\rm
eff}(r)\approx M_{\rm eff}(R_{\odot})$ and we get

$1-\frac{8\pi G}{c^{4}}\eta\approx\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}+\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}r_{s}\left(1-\frac{r_{s}}{r}\right)-\frac{2G}{3c^{2}}\frac{M_{\rm
eff}(R_{\odot})}{r}+\frac{8\pi G}{3c^{2}}\int_{r_{s}}^{R_{\odot}}\rho(r)rdr,$
(173)

from which, using the inequality

$4\pi\int_{r_{s}}^{R_{\odot}}\rho(r)rdr>\frac{4\pi}{R_{\odot}}\int_{r_{s}}^{R_{\odot}}\rho(r)r^{2}dr\approx\frac{M_{\rm
eff}(R_{\odot})}{R_{\odot}},$ (174)

we obtain

$\displaystyle 1-\frac{8\pi G}{c^{4}}\eta$ $\displaystyle>$
$\displaystyle\frac{16\pi G}{c^{2}}q_{1}\rho_{s}+\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}^{\prime}r_{s}\left(1-\frac{r_{s}}{r}\right)$ (175)
$\displaystyle+$ $\displaystyle\frac{2G}{3c^{2}}\,\frac{M_{\rm
eff}(R_{\odot})}{R_{\odot}}\left(1-\frac{R_{\odot}}{r}\right).$

Using now inequality (99), valid for $r_{s}$ in the convection zone, we find
the following lower bounds:

$\displaystyle 1-\frac{8\pi G}{c^{4}}\eta$ $\displaystyle>$
$\displaystyle\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}+\frac{2G}{3c^{2}}\,\frac{M_{\rm
eff}(R_{\odot})}{r_{s}}\left[(1-20\varepsilon)\left(1-\frac{r_{s}}{r}\right)-\left(1-\frac{r_{s}}{R_{\odot}}\right)\right],\qquad\mbox{for
}q_{1}>0,$ $\displaystyle 1-\frac{8\pi G}{c^{4}}\eta$ $\displaystyle>$
$\displaystyle\frac{16\pi
G}{c^{2}}q_{1}\left\\{\rho_{s}+\rho_{s}^{\prime}r_{s}\left[1-\frac{r_{s}}{r}\left(1+\frac{1}{20\varepsilon}\left(1-\frac{r}{R_{\odot}}\right)\right)\right]\right\\},\qquad\mbox{for
}q_{1}<0.$ (176)

Let us first consider the case $q_{1}>0$. In this case the lower bound can be
written in the form

$1-\frac{8\pi G}{c^{4}}\eta>\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}+\phi(r,r_{s},\varepsilon),$ (177)

where the function $\phi(r,r_{s},\varepsilon)$ is positive if the following
condition is satisfied:

$r>\left[1+\delta(\varepsilon,r_{s})\right]R_{\odot},$ (178)

with

$\delta(\varepsilon,r_{s})=\frac{20\varepsilon(1-r_{s}/R_{\odot})}{r_{s}/R_{\odot}-20\varepsilon}.$
(179)

Let now $r$ in the solar atmosphere, $r_{s}$ in the convection zone, and
$\varepsilon>0$ be such that $\phi(r,r_{s},\varepsilon)$ is positive and
$\rho(r)\ll\rho_{s}=\rho(r_{s})$. Then we have

$\frac{16\pi G}{c^{2}}q_{1}\rho(r)\ll\frac{16\pi
G}{c^{2}}q_{1}\rho_{s}<1-\frac{8\pi G}{c^{4}}\eta,$ (180)

so that inequality (172) is satisfied if the condition
$|\alpha|\rho(r)\ll\rho_{s}$ is also satisfied.

For $r_{s}<r_{p}$, where $r_{p}=R_{\odot}-500\,{\rm km}$ is the radius at the
base of the photosphere (see Appendix B), the function
$\delta(\varepsilon,r_{s})$ is decreasing as $r_{s}$ increases and increasing
as $\varepsilon$ increases. Particularly, for $\varepsilon$ small enough, in
order to have inequality (172) satisfied for $r$ varying in the inner regions
of the solar atmosphere, $r_{s}$ has to be close enough to $R_{\odot}$, which
means that the contribution to the effective mass only comes from a thin shell
in the upper part of the solar interior, but this is just the requirement
imposed by the constraint from the Cassini measurement (see Section V).

Let us give some examples that show for which values of
$r,r_{s},\varepsilon,\alpha$, with $r_{s}$ satisfying the Cassini constraint,
the desired inequality (172) is satisfied. Let us consider the value
$r_{s}=0.98R_{\odot}$ which corresponds to $\rho_{s}\approx 3\times
10^{-4}\,{\rm g\,cm^{-3}}$.

For $\varepsilon=10^{-2}$ we have $\delta=0.005$ so that, using the density
model of the solar atmosphere reported in Appendix B,
$r_{0}=(1+\delta)R_{\odot}$ is located in the chromosphere-corona transition
region, and for $r>r_{0}$ we have $\rho(r)<\rho(r_{\rm cr})=9.82\times
10^{-15}\,{\rm g\,cm^{-3}}$ (density at the top of the chromosphere). Then
$|\alpha|\rho(r)\ll\rho_{s}$ for $|\alpha|<10^{8}$ so that for such values of
$\alpha$ inequality (172) is satisfied.

For $\varepsilon=10^{-3}$ we have $\delta=4\times 10^{-4}$, with
$r_{0}=(1+\delta)R_{\odot}$ located in the chromosphere so that, using the
density profile of the chromosphere reported in Appendix B, for $r>r_{0}$ we
have $\rho(r)<\rho(1.0004R_{\odot})\approx 10^{-9}\,{\rm g\,cm^{-3}}$. In this
case inequality (172) is satisfied for $|\alpha|<10^{3}$.

If now $r_{s}$ increases and reaches the top of the convection zone,
$r_{s}=r_{p}$, then $\rho_{s}=2.73\times 10^{-7}\,{\rm g\,cm^{-3}}$ and the
Cassini constraint is satisfied. For $\varepsilon=10^{-2}$ we have
$\delta=1.8\times 10^{-4}$, and now we need to take values of $r$ larger than
the ones implied by the inequality $r>(1+\delta)R_{\odot}$. If we take
$r_{0}=1.001R_{\odot}$ in the chromosphere, then for $r>r_{0}$ we have
$\rho(r)<6.99\times 10^{-11}\,{\rm g\,cm^{-3}}$, and inequality (172) is
satisfied for $|\alpha|<39$.

We consider now the case $q_{1}<0$. We need the expression of
$\rho_{s}^{\prime}r_{s}$ which is obtained by computing the derivative of the
density profile in the convection zone given by Eqs. (185) and (187) in
Appendix B:

$\rho_{s}^{\prime}r_{s}=-n_{c}\frac{GM_{\odot}}{C_{p}}\frac{\rho_{s}}{GM_{\odot}/C_{p}-T_{0}r_{s}},$
(181)

where $n_{c}$ is an effective polytropic index, $C_{p}$ is an averaged value
of the specific heat at constant pressure, and $T_{0}$ is a reference
temperature (see Appendix B for details and numerical values). Then the lower
bound (176) can be written in the form

$1-\frac{8\pi G}{c^{4}}\eta>\frac{16\pi
G}{c^{2}}\left|q_{1}\right|\rho_{s}\left[\psi(r,r_{s},\varepsilon)-1\right],$
(182)

where

$\psi(r,r_{s},\varepsilon)=n_{c}\frac{GM_{\odot}/C_{p}}{GM_{\odot}/C_{p}-T_{0}r_{s}}\left[1-\frac{r_{s}}{r}\left(1+\frac{1}{20\varepsilon}\left(1-\frac{r}{R_{\odot}}\right)\right)\right].$
(183)

The function $\psi(r,r_{s},\varepsilon)$ is positive in the domain $D=\\{r\geq
R_{\odot},0\leq r_{s}\leq r_{p},\varepsilon\geq 0\\}$, and in this domain it
achieves the absolute minimum, independent of $\varepsilon$, at $r=R_{\odot}$,
$r_{s}=r_{p}$ with value ${\rm min}_{D}\psi(r,r_{s},\varepsilon)\approx
1.918$. Let now $r$ in the solar atmosphere and $r_{s}$ in the convection zone
be such that $|\alpha|\rho(r)\ll 0.918\rho_{s}$. Then we have

$\displaystyle\frac{16\pi G}{c^{2}}\left|\alpha q_{1}\right|\rho(r)\ll
0.918\frac{16\pi G}{c^{2}}\left|q_{1}\right|\rho_{s}$ (184) $\displaystyle<$
$\displaystyle\frac{16\pi
G}{c^{2}}\left|q_{1}\right|\rho_{s}\left[\psi(r,r_{s},\varepsilon)-1\right]<1-\frac{8\pi
G}{c^{4}}\eta,$

so that, being $q_{1}<0$, inequality (172) is satisfied for any
$\varepsilon>0$.

We now give two examples. Let $r$ and $r_{s}$ be such that $r>R_{\odot}$ and
$r_{s}=0.98R_{\odot}$. Then, using the density profile in the chromosphere, we
have $\rho(r)<5.9\times 10^{-9}\,{\rm g\,cm^{-3}}$ (density at the bottom of
the chromosphere), and inequality (172) is satisfied for $|\alpha|<10^{2}$. If
now $r$ and $r_{s}$ are such that $r>r_{0}$, with $r_{0}=1.001R_{\odot}$ in
the chromosphere, and $r_{s}=r_{p}$ at the top of the convection zone, then
inequality (172) is satisfied for $|\alpha|<35$.

Analogous results can be obtained if the screening radius is located in the
photosphere or outside of the photosphere.

## Appendix B

We report on a model of solar mass density profile $\rho(r)$, both for the
interior and the atmosphere, which is used in order to find analytical
estimates of the constraints on the NMC gravity parameters at a suitable order
of magnitude. Matter in the Sun is modeled as a perfect gas in hydrostatic
equilibrium, with the exception of the outer corona, where dynamical
equilibrium of a stationary atmosphere is used.

###  .3 Solar interior

The radius of the Sun $R_{\odot}=6.9634\times 10^{5}\,{\rm km}$ is defined to
be the radius of the edge, or limb, of the Sun when observed in white light.
The solar atmosphere begins below the spherical surface of radius $R_{\odot}$
and center in the origin Tan-Han , at a depth of about 500 km, and extends
outward from the Sun. Then $r_{p}=R_{\odot}-500\,{\rm km}$ is the radius at
the base of the photosphere. The density profile is computed for Newtonian
gravity, while a computation based on the solution of the Lane-Emden equation
in the presence of an Yukawa force can be found in Ref. laneemdenyukawa .

####  .3.1 Convection zone

This outer zone of the solar interior is important since the screening radius
that saturates the Cassini bound (125) lies in such a zone. We use a polytrope
model with an effective polytropic index $n_{c}=2.33$ Mullan :

$\rho(r)=K_{c}\left[T(r)\right]^{n_{c}},$ (185)

with $K_{c}=3.44\times 10^{-16}$, and the radius $r$ varying in the range

$r_{\rm conv}\leq r<r_{p},\qquad r_{\rm conv}=5.3185\times 10^{5}\,\mbox{km}.$
(186)

The temperature profile is approximated by Mullan

$T(r)=\frac{GM_{\odot}}{C_{p}r}-T_{o},$ (187)

with $M_{\odot}=1.989\times 10^{33}\,{\rm g}$, $C_{p}=2.95\times 10^{8}\,{\rm
erg\,g^{-1}\,K^{-1}}$ is an averaged value of the specific heat at constant
pressure, and $T_{o}=6.461\times 10^{6}\,{\rm K}$.

The density values computed by means of the above model and by means of solar
models that incorporate helioseismological measurements Bahcall  are of the
same order of magnitude.

In order to verify the consistency condition (39) the expression of the
Laplacian of density is required:

$\nabla^{2}\rho(r)=n_{c}(n_{c}-1)\left(\frac{GM_{\odot}}{C_{p}}\right)^{2}\frac{K_{c}}{r^{4}}\left[T(r\right)]^{n_{c}-2}.$
(188)

We have $\nabla^{2}\rho(r)>0$ in the convection zone and $\nabla^{2}\rho/\rho$
is an increasing function, so that the maximum of the quantity
$|\nabla^{2}\rho|/\rho$ is achieved at $r=r_{p}$ and it is of order of
$10^{7}/R_{\odot}^{2}$.

####  .3.2 Radiative interior

The Cassini bound prevents the screening radius from penetrating inside the
radiative interior, hence the density profile in this region will only be used
for verifying the consistency condition (38). We use the polytrope model with
$n_{r}=3$:

$\rho(r)=K_{r}\left[T_{6}(r)\right]^{n_{r}},\qquad r<r_{\rm conv},$ (189)

with $K_{r}=2.6885\times 10^{-2}$ and $T=10^{6}T_{6}\,{\rm K}$. The
temperature profile is estimated by resorting to an approximation of the
gravitational acceleration profile inside the zone (see Ch. 9 of Ref. Mullan
for details):

$T(r)=T(r_{\rm conv})+T_{\ast}\left[\chi(r)-\frac{R_{\odot}}{r_{\rm
conv}}\right],$ (190)

where

$T_{\ast}=\frac{\mu m_{p}}{(n_{r}+1)k_{B}}\,\frac{GM_{\odot}}{R_{\odot}},$
(191)

and

$\chi(r)=\begin{cases}2\left[3-\left(4r/R_{\odot}\right)^{2}\right]&0\leq
r<R_{\odot}/4,\\\ R_{\odot}/r&R_{\odot}/4\leq r<r_{\rm conv},\end{cases}$
(192)

$T(r_{\rm conv})=2\times 10^{6}\,{\rm K}$, $m_{p}=1.66\times 10^{-24}\,{\rm
g}$ is the proton mass, $\mu=0.62$ is the mean molecular weight, and
$k_{B}=1.3806\times 10^{-16}\,{\rm erg\,K^{-1}}$ is the Boltzmann constant.
This density profile is used up to $r=0$.

Since the temperature profile is continuous across the base $r=r_{\rm conv}$
of the convection zone and the polytropic index changes from the convection
zone to the radiative interior, the density model introduces an artificial
discontinuity at $r=r_{\rm conv}$. Since the density profile in the inner
region of the convection zone and in the radiative interior is only used for
verifying the consistency condition (39), which requires the computation of
the Laplacian of mass density, then the consistency condition is verified in
the separate regions without crossing the discontinuity.

For $0\leq r<R_{\odot}/4$ the Laplacian of density is given by

$\displaystyle\nabla^{2}\rho$ $\displaystyle=$ $\displaystyle 6.4\times
10^{-11}n_{r}\frac{K_{r}T_{\ast}}{R_{\odot}^{2}}\left[T_{6}(r)\right]^{n_{r}-2}$
$\displaystyle\times$
$\displaystyle\left[T_{\ast}\left(224\frac{r^{2}}{R_{\odot}^{2}}-3\frac{R_{\odot}}{r_{\rm
conv}}-18\right)-3T(r_{\rm conv})\right],$

and for $R_{\odot}/4\leq r<r_{\rm conv}$ it is given by

$\nabla^{2}\rho=10^{-12}n_{r}(n_{r}-1)K_{r}T_{\ast}^{2}\left[T_{6}(r)\right]^{n_{r}-2}\frac{R_{\odot}^{2}}{r^{4}}.$
(194)

Using the above expressions we find that the maximum of the quantity
$|\nabla^{2}\rho|/\rho$ in the radiative interior is achieved at $r\approx
0.22R_{\odot}$ and it is of order of $10^{2}/R_{\odot}^{2}$.

###  .4 Solar atmosphere

Large density gradients take place in the solar atmosphere, particularly in
the chromosphere-corona transition region. Then the extra force, which depends
on the density gradient, can become a significant perturbation of hydrostatic
equilibrium in such regions with large gradients.

####  .4.1 Photosphere

For the density in the photosphere we use the following model, adapted from
DeJager , for $r_{p}\leq r<R_{\odot}$:

$\rho(r)=\frac{\mu
m_{p}p_{m}}{k_{B}[T_{m}+A(R_{\odot}-r)^{2}])}\exp\left(\frac{R_{\odot}-r}{H_{p}}\right),$
(195)

where $\mu=1.26$, $H_{p}=117\,{\rm km}$, $A=8.8\times 10^{-3}\,{\rm
km^{-2}K}$, $T_{m}=4.4\times 10^{3}\,{\rm K}$ is the temperature minimum at
the top of the photosphere, $p_{m}$ is pressure corresponding at the
temperature minimum, such that $p_{m}/k_{B}=1.2\times 10^{19}\,{\rm
K\,cm^{-3}}$.

####  .4.2 Chromosphere

The middle chromosphere is characterized by a broad temperature plateau at
$T\approx 6500$ K Vernazza , hence, for simplicity, we adopt an isothermal
model for this layer of the solar atmosphere:

$\rho(r)=\rho_{p}\exp\left(-\frac{r-R_{\odot}}{H_{c}}\right),\quad
R_{\odot}\leq r<1.003R_{\odot},$ (196)

with

$H_{c}=\frac{k_{B}TR_{\odot}^{2}}{\mu m_{p}GM_{\odot}}\approx 157\,{\rm km},$
(197)

$\mu=1.26$, and $\rho_{p}=5.9\times 10^{-9}\,{\rm g\,cm^{-3}}$ is density at
the top of the photosphere (value interpolated from the average quiet Sun
model in Vernazza ). The density profile (196) coincides with the best fitting
formula used in Kontar .

####  .4.3 Chromosphere-corona transition region

Here the steepest density gradient takes place. The density profile in this
region is computed by using hydrostatic equilibrium together with the average
emission measure distribution given in Eq. (146). An analytic estimate of
density is found by approximating the integral in Eq. (146) as in Jordan76 ,
see Section VI. The resulting density profile $\rho\approx\mu
m_{p}p_{0}/(k_{B}T)$ is given by

$\rho(r)=\mu
m_{p}\frac{p_{0}}{k_{B}}\left[T_{0}^{b+2}+0.231\frac{b+2}{\xi^{2}a}\left(\frac{p_{0}}{k_{B}}\right)^{2}(r-r_{0})\right]^{-\frac{1}{b+2}},$
(198)

with $T_{0}=T(r_{0})=1.5\times 10^{5}\,{\rm K}$, $p_{0}/k_{B}=2\times
10^{15}\,{\rm K\,cm^{-3}}$, and we use the values of $a,b$ reported in
RayDoyle  for the average quiet Sun. The transition region is divided into two
zones:

  * (i)

zone with temperature $2\times 10^{4}\,{\rm K}\leq T<1.5\times 10^{5}\,{\rm
K}$; in this zone we use $\mu=0.65$, $\xi=2$, $b=-3.5$, $a=1\times 10^{43}$,
and the zone corresponds to

$0.003R_{\odot}\approx 2089\,{\rm km}\leq r-R_{\odot}<2186\,{\rm km};$ (199)

  * (ii)

zone with temperature $1.5\times 10^{5}\,{\rm K}\leq T<1\times 10^{6}\,{\rm
K}$; in this zone we use $\mu=0.62$, $\xi=1.91$, $b=1.65$, $a=2\times
10^{16}$, and the zone corresponds to

$2186\,{\rm km}\leq r-R_{\odot}<3906\,{\rm km}.$ (200)

Density exhibits the steepest gradient in zone (i) which is denoted the lower
transition region. The thickness of the transition zone is in good agreement
with the empirical atmospheric model computed in Mariska-book  (Ch. 6, Fig.
6.1) using the emission measures tabulated in RayDoyle .

####  .4.4 Inner corona

We consider the inner corona approximately isothermal, so that the density
profile is

$\rho(r)=\rho_{\rm ic}\exp\left[\frac{\mu
m_{p}GM_{\odot}}{k_{B}T}\left(\frac{1}{r}-\frac{1}{r_{\rm ic}}\right)\right],$
(201)

with $\mu=0.62$,

$r_{\rm ic}\leq r<r_{a},\quad r_{\rm ic}=R_{\odot}+3906\,{\rm km},\quad
r_{a}=1.06R_{\odot},$ (202)

$r_{a}$ is the base of the outer corona, $\rho_{\rm ic}=2\times 10^{-15}\,{\rm
g\,cm^{-3}}$ is density at the base of the inner corona, and $T\approx
10^{6}\,{\rm K}$.

####  .4.5 Outer corona

The model is based on the stationary expansion of the outer corona in which
temperature $T$ is assumed constant between the coronal base at $r=r_{a}$ and
an outer boundary at $r=r_{b}$, with adiabatic expansion beyond $r_{b}$ Stix ;
Parker . The flow of the expanding atmosphere turns from subsonic to
supersonic, giving rise to the solar wind, as it passes through a critical
radius $r_{c}$ such that $r_{a}<r_{c}<r_{b}$. Such a model yields values of
solar wind speed and density which are in good agreement with measurements
from the HELIOS solar probes Schwenn  and Parker Solar Probe Allen .

Subsonic zone. This zone corresponds to the range of distances $r_{a}\leq
r<r_{c}$, where the critical radius $r_{c}$ is given by Stix ; Parker

$r_{c}=\frac{\mu m_{p}GM_{\odot}}{2k_{B}T_{a}},$ (203)

where we choose $T_{a}=1\times 10^{6}\,{\rm K}$ for the constant value of
temperature in the range $(r_{a},r_{b})$, and $\mu=0.54$ which corresponds to
the average value 0.032 of the helium abundance Schwenn , the ratio of helium
to proton number density. Then we get $r_{c}\approx 6.19R_{\odot}$. If we set
$\psi=v^{2}/c_{s}^{2}<1$, where $v=v(r)$ is the velocity of the expanding gas
and $c_{s}$ is the isothermal speed of sound, in the region where
$\psi\ll|\ln(\psi)|$ the density profile is approximated by

$\rho(r)=\rho_{a}\exp\left[2r_{c}\left(\frac{1}{r}-\frac{1}{r_{a}}\right)\right],$
(204)

where $\rho_{a}=9.7\times 10^{-16}\,{\rm g\,cm^{-3}}$ is density at the base
of the outer corona.

Supersonic zone. This zone corresponds to distances $r\geq r_{c}$ and it is
further divided in two subzones. In the first subzone the expansion is
isothermal and it takes place in the range $r_{c}\leq r<r_{b}$. We use the
value $r_{b}=192R_{\odot}\approx 0.89$ AU which corresponds to a velocity
$v(r_{b})\approx 460\,{\rm km\,s^{-1}}$ of the solar wind Parker , in good
agreement with observations Schwenn ; Allen .

In the supersonic zone $\psi>1$, and at distances $r$ such that
$\psi\gg\ln(\psi)$ the density profile in the first subzone is approximated by

$\rho(r)=\rho_{a}\,\frac{\exp\left[3/2-2(r_{c}/r_{a})\right]}{2\left[\ln(r/r_{c})\right]^{1/2}}\left(\frac{r_{c}}{r}\right)^{2},$
(205)

and we use this approximation in the range $(r_{c},r_{b})$. In the second
subzone $r\geq r_{b}$, the expansion is adiabatic, and density is approximated
by

$\rho(r)=\rho_{a}\,\frac{\exp\left[3/2-2(r_{c}/r_{a})\right]}{2\left[\ln(r_{b}/r_{c})\right]^{1/2}}\left(\frac{r_{c}}{r}\right)^{2},$
(206)

hence density decreases as $1/r^{2}$.

## Acknowledgments

The work of R.M. is partially supported, and the work of M.M. and S.DA is
fully supported, by INFN (Istituto Nazionale di Fisica Nucleare, Italy), as
part of the MoonLIGHT-2 experiment in the framework of the research activities
of the Commissione Scientifica Nazionale n. 2 (CSN2). C.G. is supported by the
Fundo Regional para a Ciência e Tecnologia and Azores Government Grant No.
M3.2DOCPROF/F/008/2020.

## References

  * (1) S. Capozziello, Int. J. Mod. Phys. D 11, 483 (2002).
  * (2) S. M. Carroll, V. Duvvuri, M. Trodden and M.S. Turner, Phys. Rev. D 70, 043528 (2004).
  * (3) S. Capozziello and M. De Laurentis, Phys. Repts. 509, 167 (2011).
  * (4) A. De Felice and S. Tsujikawa, Liv. Rev. Rel. 13, 3 (2010).
  * (5) O. Bertolami, C. G. Böhmer, T. Harko and F. S. N. Lobo, Phys. Rev. D 75, 104016 (2007).
  * (6) O. Bertolami and J. Páramos, JCAP 009, 1003 (2010).
  * (7) O. Bertolami, P. Frazão and J. Páramos, Phys. Rev. D 81, 104046 (2010).
  * (8) C. Gomes, J. G. Rosa and O. Bertolami, JCAP 06, 021 (2017).
  * (9) O. Bertolami, P. Frazão and J. Páramos, JCAP 05, 029 (2013).
  * (10) O. Bertolami, C. Gomes and F. S. N. Lobo, Eur. Phys. J. C 78, 303 (2018).
  * (11) O. Bertolami and C. Gomes, JCAP 09, 010 (2014).
  * (12) C. Gomes, Eur. Phys. J. C 80, 633 (2020).
  * (13) T. D. Ferreira, N. A. Silva, O. Bertolami, C. Gomes and A. Guerreiro, Phys. Rev. E 101, 023301 (2020).
  * (14) T.D. Ferreira, J. Novo, N.A. Silva, A. Guerreiro and O. Bertolami, Phys.Rev. D 103, 124019 (2021).
  * (15) R. March, J. Páramos, O. Bertolami and S. Dell’Agnello, Phys. Rev. D 95, 024017 (2017).
  * (16) O. Bertolami, F. Lobo and J. Páramos, Phys. Rev. D 78, 064036 (2008).
  * (17) R. March, O. Bertolami, M. Muccino, R. Baptista and S. Dell’Agnello, Phys. Rev. D 100, 042002 (2019).
  * (18) M.A. Zumberge et al., Phys. Rev. Lett. 67, 3051 (1991).
  * (19) J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004).
  * (20) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007).
  * (21) B. Bertotti, L. Iess and P. Tortora, Nature 425, 374 (2003).
  * (22) S.B. Fisher and E.D. Carlson, arXiv:2108.13492 [gr-qc].
  * (23) O. Bertolami and J. Páramos, Class. Quant. Grav. 25, 245017 (2008).
  * (24) T. P. Sotiriou and V. Faraoni, Class. Quant. Grav. 25, 5002 (2008).
  * (25) O. Bertolami and J. Páramos, Phys. Rev. D 77, 084018 (2008).
  * (26) C.M. Will, Theory and Experiment in Gravitational Physics, Revised Ed. (Cambridge University Press, Cambridge, 1993).
  * (27) M. Stix, The Sun - An Introduction (Springer, Berlin, 2002).
  * (28) D.A.Gurnett and W.S. Kurth, Nature Astron. 3, 1024 (2019).
  * (29) A-C. Davies, E.A. Lim, J. Sakstein and D.J. Shaw, Phys. Rev. D 85, 123006 (2012).
  * (30) J. Christensen-Dalsgaard, D.O. Gough and M.J. Thompson, Astrophys. J. 378, 413 (1991).
  * (31) J.T. Mariska, Ann. Rev. Astron. Astrophys. 24, 23 (1986).
  * (32) C. Jordan, Phil. Trans. R. Soc. Lond. A 281, 391 (1976).
  * (33) C. Jordan Astron. Astrophys. 86, 355 (1980).
  * (34) J.C. Raymond and J.G. Doyle, Astrophys. J. 247, 686 (1981).
  * (35) G.A. Doschek, U. Feldman, J.M. Laming, H.P. Warren, U. Schühle and K. Wilhelm, Astrophys. J. 507, 991 (1998).
  * (36) D. Müller et al., Astron. Astrophys. 642, A1 (2020).
  * (37) N.J. Fox et al., Space Sci. Rev. 204, 7 (2016).
  * (38) E. Tandberg-Hanssen, Solar Activity (Blaisdell Publishing Company, Waltham, Massachusetts, 1967).
  * (39) O. Bertolami and J. Páramos, Phys. Rev. D 71, 023521 (2005).
  * (40) D.J. Mullan, Physics of the Sun: A First Course (CRC Press, Boca Raton, Florida, 2010).
  * (41) J.N. Bahcall, A.M. Serenelli and S. Basu, Astrophys. J. Suppl. 165, 400 (2006).
  * (42) C. De Jager, Bull. Astron. Institut. Netherlands 17, 209 (1962).
  * (43) J.E. Vernazza, E.H. Avrett and R. Loeser, Astrophys. J. Suppl. 45, 635 (1981).
  * (44) E.P. Kontar, I.G. Hannah, N.L.S. Jeffrey and M. Battaglia, Astrophys. J. 717, 250 (2010).
  * (45) J.T. Mariska, The Solar Transition Region (Cambridge University Press, New York, 1992).
  * (46) E.N. Parker, Astrophys. J. 128, 664 (1958).
  * (47) R. Schwenn, JPL Solar Wind Five (SEE N84-13067 03-92), 489 (1983).
  * (48) R.C. Allen et al., Astrophys. J. Suppl. 246, 36 (2020).