# Structural Implications of the Chameleon Mechanism on White Dwarfs

Joan Bachs-Esteban<EMAIL_ADDRESS>Ilídio Lopes
<EMAIL_ADDRESS>Centro de Astrofísica e Gravitação - CENTRA,
Departamento de Física, Instituto Superior Técnico - IST, Universidade de
Lisboa - UL, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal  Javier Rubio
<EMAIL_ADDRESS>Departamento de Física Teórica and Instituto de Física de
Partículas y del Cosmos (IPARCOS), Facultad de Ciencias Físicas, Universidad
Complutense de Madrid, 28040 Madrid, Spain

###### Abstract

We study the impact of the chameleon mechanism on the structure of white
dwarfs. Using a shooting method of our design, we solve the corresponding
scalar-tensor equilibrium equations for a Chandrasekhar equation of state,
exploring various energy scales and couplings of the chameleon field to
matter. We find the chameleon field to be in a thick-shell configuration,
identifying for the first time in the literature a similarity relation for the
radially normalised scalar field gradient. Our analysis reveals that the
chameleon mechanism significantly alters the internal pressure of white
dwarfs, leading to a notable reduction in the stellar radii and masses and
shifting the mass-radius curves below those predicted by Newtonian gravity.
Finally, we derive parametric expressions from our results to expedite future
analyses of white dwarfs in scalar-tensor theories.

##  I Introduction

General Relativity (GR) stands strong as a gravity theory, verified by
numerous experiments and observations [1]. Nevertheless, it fails to explain
several astrophysical and cosmological phenomena like the rotation curves of
galaxies or the fundamental origin of dark energy. This allows gravity not to
be described by GR, but by a modified gravity (MG) theory [2].

Among the numerous extensions of GR, scalar-tensor (ST) theories of gravity
[3, 4, 5] stand out as one of the most simple yet elegant proposals. On
general grounds, these types of settings introduce one or more scalar fields
potentially mediating a fifth force among matter components, which, if
sufficiently long-range, could potentially contradict the local tests of
gravity [1] or enhance structure formation in the early [6, 7, 8] and late
Universe [9, 10, 11, 12]. Many ST theories include, however, screening
mechanisms that render the scalar field properties environment-dependent. For
instance, in the chameleon mechanism [13], the mass of the scalar field varies
with the environment. In the symmetron [14] and dilaton scenarios [15], it is
the coupling to matter that changes, while in the Vainshtein [16] and
k-mouflage implementations [17], it is the kinetic function that plays a role.
These screening mechanisms ensure that the scalar field becomes ineffective on
astrophysical scales while remaining potentially relevant at cosmological
scales.

The partial breaking of screening mechanisms within massive sources is
expected to influence the equilibrium structure of stars, altering with it
fundamental properties such as their mass-radius relations and cooling times
(see e.g. [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]). Among the various compact
objects suitable for studying the effects of screening mechanisms, white
dwarfs (WDs) are particularly promising yet relatively unexplored targets [28,
29, 30]. This is due to two key reasons. Firstly, the equation of state (EoS)
describing the microscopic behaviour of matter inside WDs is fairly well
understood [31]. Secondly, the extensive observational data now available from
Gaia’s data releases [32, 33, 34] provides a wealth of information on WDs,
facilitating comprehensive studies of their spatial distribution, kinematics,
and fundamental properties such as luminosity, temperature, and radius.

In this study, we focus on exploring the impact of the chameleon mechanism on
the structure of WDs. To this end, we numerically solve the ST equilibrium
equations in the relativistic and Newtonian limits and employ a Chandrasekhar
EoS. By exploring a broad range of energy scales and conformal chameleon
couplings to matter, we determine the corresponding mass-radius relations to
be confronted with observations, deriving also a set of ready-to-use fitting
formulae aiming to streamline future analyses of WDs in ST theories.

This paper is structured as follows: Section II introduces the physical
framework, discussing the ST theory and the equilibrium equations for static
and spherically symmetric WDs together with the EoS and observational data for
these stars. Section III describes the specific screening mechanism under
consideration and outlines the employed numerical methods, including boundary
conditions and our customized shooting method. Section IV presents our
findings regarding chameleon-screened WDs, followed by a discussion of
implications and future directions in Section V. Finally, Appendix A provides
supplementary details on the validity of the Newtonian approximation in our ST
framework. We use the metric signature $(-,+,+,+)$ and consider $c=\hbar=1$
unless otherwise stated.

##  II Framework

###  II.1  Scalar-Tensor Theory

Several ST theories with environmentally dependent screening mechanisms such
as chameleons, symmetrons or dilatons can be described by a general action
[20]

$\begin{split}S=&\int
d^{4}x\sqrt{-g}\left[\frac{M_{P}^{2}}{2}R-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi-V(\phi)\right]\\\
&+S_{m}\left[\Psi_{m};A^{2}(\phi)g_{\mu\nu}\right],\end{split}$ (1)

with $M_{P}=(8\pi G)^{-1/2}=2.43\times 10^{18}$ GeV the reduced Planck mass,
$g$ and $R$ the determinant and Ricci scalar of the Einstein frame metric
$g_{\mu\nu}$ and $\phi$ a scalar field. Each model belonging to this class of
theories is characterised by a self-interacting potential $V(\phi)$ and a
conformal coupling $A(\phi)$ to the matter fields $\Psi_{m}$. In particular,
the scalar field is taken to be gravitationally coupled to the matter fields
through a conformally rescaled Jordan frame metric $\tilde{g}_{\mu\nu}\equiv
A^{2}(\phi)g_{\mu\nu}$. This modifies the Newtonian force in the non-
relativistic limit, making these ST theories MG theories [20].

By varying Eq. (1) with respect to the metric, one obtains the field equations

$G_{\mu\nu}=\kappa^{2}\left[T_{\mu\nu}+\nabla_{\mu}\phi\nabla_{\nu}\phi-
g_{\mu\nu}\left(\frac{1}{2}\nabla_{\sigma}\phi\nabla^{\sigma}\phi+V(\phi)\right)\right]\,,$
(2)

where $G_{\mu\nu}$ is the Einstein tensor, $\kappa\equiv M_{P}^{-1}$, and

$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta g^{\mu\nu}}$
(3)

is the energy-momentum tensor of the matter fields, which we assume to be
described by a perfect fluid, namely

$T^{\mu\nu}\equiv(\epsilon+p)u^{\mu}u^{\nu}+pg^{\mu\nu}\,,$ (4)

with $u^{\mu}$ the four-velocity of fluid elements, and $\epsilon$ and $p$ the
total energy density and pressure in the fluid’s rest frame, respectively.
Analogously, if we vary Eq. (1) with respect to the field, we obtain the
scalar field equation

$\Box\phi=\frac{dV(\phi)}{d\phi}-\frac{d\text{ ln
}A(\phi)}{d\phi}T\equiv\frac{dV_{\rm eff}(\phi)}{d\phi}\,,$ (5)

with $V_{\rm eff}(\phi)$ the potential effectively governing $\phi$. The
matter equation of motion is determined by the divergence of Eq. (1), namely

$\nabla^{\nu}T_{\mu\nu}=\frac{d\text{ ln }A(\phi)}{d\phi}T\nabla_{\mu}\phi\,,$
(6)

with $T\equiv g^{\mu\nu}T_{\mu\nu}$ denoting the trace of the energy-momentum
tensor. The above equation means that particles do not follow geodesics in the
Einstein frame metric $g_{\mu\nu}$, being their trajectories also affected by
the scalar field gradient. Nevertheless, in an alternative, but equally valid
description, we can analogously define the Jordan frame matter energy-momentum
tensor as

$\tilde{T}_{\mu\nu}\equiv-\frac{2}{\sqrt{-\tilde{g}}}\frac{\delta
S_{m}}{\delta\tilde{g}^{\mu\nu}}\,.$ (7)

Comparing the latter expression with Eq. (3), we see that both tensors are
related through $T_{\mu\nu}=A^{2}(\phi)\tilde{T}_{\mu\nu}$. Using Eq. (1), one
can show that $\tilde{T}_{\mu\nu}$ is indeed covariantly conserved – i.e.,
$\tilde{\nabla}^{\nu}\tilde{T}_{\mu\nu}=0$ – and that free particles follow
the $\tilde{g}_{\mu\nu}$ geodesics. In addition, from the four-velocity
normalisation condition $g_{\mu\nu}u^{\mu}u^{\nu}=-1$, we get the relation
$u^{\mu}=A(\phi)\tilde{u}^{\mu}$. This conformal transformation, together with
that for the energy-momentum tensor above, allows us to find the
correspondence between the fluid variables in both frames, namely
$\epsilon=A^{4}(\phi)\tilde{\epsilon}$ and $p=A^{4}(\phi)\tilde{p}$.

###  II.2 Equation of State

The equation of state (EoS) condenses the microphysics of the stellar interior
in a relation between pressure and density. As explained now and discussed in
Appendix A, WDs can be adequately described as non-relativistic objects, both
in GR and the ST scenarios considered here. Therefore, it is useful to
introduce the rest-mass density $\tilde{\rho}$ and the internal energy density
$\tilde{\Pi}$, which are related to the total energy density
$\tilde{\epsilon}$ as

$\tilde{\epsilon}\equiv\tilde{\rho}\left(1+\frac{\tilde{\Pi}}{\tilde{\rho}c^{2}}\right)\,,$
(8)

where we have explicitly written the speed of light $c$ to evince that
$\tilde{\Pi}$ is a first-order relativistic correction.

Moreover, the pressure and energy in WDs are nonthermal, in the sense that
thermal effects can be modelled as small perturbations on top of the fluid
dynamics [31]. Consequently, the EoS reduces to single parameter functions,
namely $\tilde{P}(\tilde{\rho})$ and $\tilde{\epsilon}(\tilde{\rho})$. Notice
that we are relating the EoS to the Jordan-frame variables, as the common
thermodynamic relation for energy conservation,
$d(\tilde{\epsilon}/\tilde{\rho})=-\tilde{P}d(1/\tilde{\rho})$, holds only in
this frame.

In WDs, the electrostatic energy of the matter structure is negligible as
compared to the Fermi energies, hence the Coulomb forces are too. Therefore,
the electron pressure is given by [31]

$\tilde{P}=\frac{2}{(2\pi)^{3}\hbar^{3}}\int_{0}^{p_{F,e}}\frac{p^{2}c^{2}}{\sqrt{p^{2}c^{2}+\left(m_{e}c^{2}\right)^{2}}}4\pi
p^{2}dp=\frac{m_{e}c^{2}}{\lambda_{e}^{3}}\psi(x)\,,$ (9)

with

$\psi(x)=\frac{1}{8\pi^{2}}\left[x\sqrt{1+x^{2}}\left(\frac{2x^{2}}{3}-1\right)+\text{ln}\left(x+\sqrt{1+x^{2}}\right)\right],$
(10)

where, for the sake of clarity, we have made explicit again the different
$\hbar$ and $c$ factors. The factor 2 in Eq. (9) is due to the electron spin
degeneracy, $p_{F,e}$ is the Fermi momentum of the electrons, $m_{e}$ is the
electron mass, $\lambda_{e}\equiv\hbar/(m_{e}c)$ is the electron Compton
wavelength, and $x\equiv p_{F,e}/m_{e}c$ is the dimensionless Fermi momentum.

Even though the main contribution to WD pressure comes from the degenerate
electrons, the energy is dominated by the ions. Since these are non-
relativistic for densities below the neutron drip
$\tilde{\rho}_{\text{n-drip}}\simeq 4\times 10^{11}\text{ g}\,\text{cm}^{-3}$
[31], the energy density can be expressed in terms of the rest-mass density

$\tilde{\epsilon}=\tilde{\rho}=\frac{m_{B}n_{e}}{Y_{e}}\,,$ (11)

where $m_{B}$ is the mean nucleon mass, $n_{e}=8\pi p_{F,e}^{3}/(3h^{3})$ is
the electron number density, and $Y_{e}=Z/A$ (with $Z$ the atomic number and
$A$ the atomic weight) is the mean number of electrons per nucleon. WDs are
usually modelled as a cold, degenerate matter star made of helium, carbon, or
oxygen 111The first EoS for such stars was derived by Chandrasekhar [44].
Hamada and Salpeter added temperature corrections to it [45].. For any of
these elements, $Y_{e}=0.5$ when they are fully ionised [31]. The mean nucleon
mass of carbon is $m_{B,C}=1.66057\times 10^{-24}\,\text{g}$. Thus, we will
express the density $\tilde{\rho}$ and the pressure $\tilde{P}$ as

$\displaystyle\tilde{\rho}(x)=$ $\displaystyle 1.9479\times
10^{6}\,x^{3}\,\text{g}\,\text{cm}^{-3}\,,$ (12) $\displaystyle\tilde{P}(x)=$
$\displaystyle 1.4218\times 10^{25}\,\psi(x)\,\text{dyn}\,\text{cm}^{-2}\,.$
(13)

Although more precise EoS exist to describe the internal structure of a WD
[31], considering the exploratory nature of this work and the advantage of
simplicity, we opted to use this equation for a more efficient calculation
method. Moreover, using a more precise EoS would not significantly alter our
results and conclusions.

###  II.3 Equilibrium Equations

In the Newtonian description, the gravitational field is weak and static and
particles move slowly compared to the speed of light [36], which means that
pressure is negligible with respect to energy density. Thus, the WD line
element can be written as

$ds^{2}=-\left(1+2\Phi(r)\right)dt^{2}+\left(1-2\Phi(r)\right)dr^{2}+r^{2}d\Omega^{2}\,,$
(14)

with $\Phi(r)$ the Newtonian potential and
$d\Omega^{2}=d\theta^{2}+\text{sin}^{2}\theta d\varphi^{2}$ is the two-sphere
line element. Replacing this metric into Eqs. (2), (6), and (5), one has that

$\displaystyle\Phi^{\prime}=$ $\displaystyle\frac{m}{r^{2}}\,,$ (15)
$\displaystyle m^{\prime}=$
$\displaystyle\frac{\kappa^{2}}{2}r^{2}A^{4}\tilde{\rho}\,,$ (16)
$\displaystyle\tilde{p}^{\prime}=$
$\displaystyle-\tilde{\rho}\left(\Phi^{\prime}+\frac{A_{,\phi}}{A}\sigma\right)\,,$
(17) $\displaystyle\phi^{\prime}=$ $\displaystyle\sigma\,,$ (18)
$\displaystyle\sigma^{\prime}=$
$\displaystyle-\frac{2}{r}\sigma+V_{,\phi}+A_{,\phi}A^{3}\tilde{\rho}\,,$ (19)

where we have ignored pressure contributions as to energy contributions and
neglected second-order terms. We also have replaced $\tilde{\epsilon}$ with
$\tilde{\rho}$ since these two quantities are equivalent in the Newtonian
limit (recall Eq. (8)).

Once one chooses the model functions $V(\phi)$ and $A(\phi)$, and a suitable
EoS, the system of differential of ordinary differential equations (ODE) can
be numerically integrated from the origin. For each central density
$\tilde{\rho}_{0}=\tilde{\rho}(0)$, we obtain the stellar mass $M$ and the
stellar radius $R$. Therefore, if we perform the integration for a range of
central densities, we can get a so-called mass-radius (MR) curve. This family
of stars is unique for each EoS and is parametrised by the central density
$\tilde{\rho}_{0}$ [31]. We provide further details of the integration process
and the boundary conditions in Sec. III.2.

###  II.4 Observational Data

Figure 1:  The 43362 WDs candidates in the [32] catalogue for which their
masses and radii (dark blue) were calculated, together with their respective
measurement errors (light blue). We also display the theoretical MR curve for
WDs in GR (black) in the Newtonian approximation, calculated through the
numerical integration of Eqs. (15), (16), and (17) with $\Lambda=\beta=0$,
using the (12)-(13) EoS for $\tilde{\rho}_{0}=7\times 10^{4}-10^{10}\text{
g}\,\text{cm}^{-3}$.

Thanks to recent astrophysical probes, we have at our disposal a vast amount
of data about WD observations. We use the mass-radius data of the catalogue
presented in [32]. It comprises 73,221 WD candidates from the astrometric and
photometric data of the Gaia DR2 (Data Release 2) catalogue. They used the
Virtual Observatory tool VOSA (see [32] and references therein) to obtain the
effective temperatures and luminosities of the WDs, from which they derived
the WD radii and the corresponding errors. They determined the masses of the
WDs and their respective errors interpolating the radii and effective
temperatures in hydrogen-rich WD cooling sequences.

We want to confront these WD data coming from astrophysical observations with
the theoretical mass and radii for WD in chameleon screening. This will help
us to discern which parameter values of our model – defined in (20) – are
compatible with current observational data. In Fig. 1 we show the
observational masses and radii available for the WDs in the cited catalogue
and the measurement errors for both magnitudes. We also display the
theoretical mass-radius curve obtained using the WD EoS given by (12) and (13)
considering Newtonian gravity. As explained in Sec. II.3, the family of stars
corresponding to a specific EoS can be parametrised through the central
density $\tilde{\rho}_{0}$. To cover the radius range in Fig. 1, we ought to
consider densities in the range $\tilde{\rho}_{0}=7\times 10^{4}-10^{10}\text{
g}\,\text{cm}^{-3}$. Notice that the upper bound does not surpass the neutron
drip density $\tilde{\rho}_{\text{n-drip}}\simeq 4\times 10^{11}\text{
g}\,\text{cm}^{-3}$ [31]. From the figure, we see that some observational
points go the above Chandrasekhar limit of $1.46\,M_{\odot}$ [31], even though
the error bars include it.

##  III Model

###  III.1 Chameleon Screening

In this work, we focus on the chameleon field [13], a scalar field equipped
with a screening mechanism through its effective mass. The environment-
dependent mass of the chameleon comes from the synergy between the self-
interacting potential $V(\phi)$ and the conformal coupling $A(\phi)$. The
potential should be monotonically decreasing and of runaway form, such that it
does not have a minimum but that, together with the coupling, they bestow a
minimum to the effective potential $V_{\rm eff}(\phi)$. We then consider a
classic inverse power-law potential and an exponential conformal coupling

$V(\phi)=\Lambda^{4}\left(\frac{\Lambda}{\phi}\right)^{n}\,,\quad\quad
A(\phi)=e^{\beta\phi/M_{P}}\,,$ (20)

where $n$ is a positive constant, $\Lambda$ has mass units, and $\beta$ is a
dimensionless constant. The $\beta$ parameter is the coupling strength between
the scalar and matter fields, and $\Lambda$ controls the scalar field
contribution to the energy density of the universe, hence we shall refer to it
as the chameleon energy scale. For $n$ and $\beta$ of order unity, equivalence
principle tests impose that $\Lambda\lesssim 10^{-30}\,M_{P}\approx
1\,\text{meV}$ [13], which remarkably coincides with the dark energy scale
causing the current accelerated expansion of the universe. Nevertheless, we do
not regard the chameleon field studied in this work as the force driving the
cosmological expansion.

The minima $\bar{\phi}$ of the chameleon field are the roots of Eq. (5), which
are determined by the transcendental equation
$nM_{P}\Lambda^{n+4}+\beta\tilde{T}\bar{\phi}^{n+1}e^{4\beta\bar{\phi}/M_{P}}=0$
since the traces of the energy-momentum tensor in the Einstein and Jordan
frames are related through the expression $T=A^{4}(\phi)\tilde{T}$. In the
$\beta\bar{\phi}/M_{P}\ll 1$ limit, we can approximate the solution by

$\bar{\phi}\approx\sqrt[n+1]{\frac{nM_{P}\Lambda^{n+4}}{-\beta\tilde{T}}}\approx\sqrt[n+1]{\frac{nM_{P}\Lambda^{n+4}}{\beta\tilde{\rho}}}\,,$
(21)

which will be real whenever $\tilde{T}<0$ since $n$, $\Lambda$, and $\beta$
are positive. This is the case in the Newtonian limit, for the trace is
$\tilde{T}=3\tilde{p}-\tilde{\rho}$ and pressure is negligible in front of
energy density, thus one has that $\tilde{T}\approx-\tilde{\rho}$, hence the
second expression.

Let us study how the chameleon fifth force is screened in WDs. For explanatory
purposes, we consider a static, spherically symmetric WD of total mass $M$,
density $\tilde{\rho}_{0}$, and radius $R$, surrounded by a medium whose
density $\tilde{\rho}_{\infty}$ is much smaller than that of the star – i.e.,
$\tilde{\rho}_{\infty}\ll\tilde{\rho}_{0}$ – for instance, the cosmological
background. From Eq. (21), we deduce that the chameleon will set to a minimum
value within the star, $\bar{\phi}_{0}$, that will be lower than the minimum
outside of it, $\bar{\phi}_{\infty}$, since $\bar{\phi}$ is inversely
proportional to the environment density.

One calculates the chameleon’s mass of small fluctuations around a potential
minimum $\bar{\phi}$ by evaluating the second derivative of the effective
potential with respect to the scalar field at $\bar{\phi}$. Therefore,
deriving Eq. (5) and considering the chameleon model in (20), we have that

$m_{{\rm eff}}^{2}\equiv\frac{d^{2}V_{{\rm
eff}}}{d\phi^{2}}=n(n+1)\frac{\Lambda^{n+4}}{\phi^{n+2}}-4\frac{\beta^{2}}{M_{P}^{2}}e^{4\beta\phi/M_{P}}\tilde{T}\,.$
(22)

Then, when we replace the scalar field value for the expression in Eq. (21),
we see that the chameleon’s effective mass for a WD is proportional to the
density since it is a non-relativistic object

$m_{{\rm
eff}}^{2}|_{\bar{\phi}}\approx\frac{n+1}{\sqrt[n+1]{n\Lambda^{n+4}}}\left(\frac{\beta}{M_{P}}\tilde{\rho}\right)^{\frac{n+2}{n+1}}\,.$
(23)

We have considered the approximation $\beta\bar{\phi}/M_{P}\ll 1$ and
neglected the $\beta^{2}$ term since the typical densities for WDs (see Sec.
II.4 and Fig. 1) are much smaller than the Planck scale, so
$\beta^{2}\tilde{\rho}/M_{P}^{2}\ll 1$. We can estimate the appropriate
chameleon energy scale for a screened WD using Eq. (23) and imposing that the
interaction range of the scalar field is of the size of the star. For a
typical WD of radius $R\sim 10^{4}$ km, central density $\tilde{\rho}_{0}\sim
10^{6}\text{ g}\,\text{cm}^{-3}$, and $n=\beta=1$, one has that $\Lambda\sim
10^{-18}\,M_{P}$. We explore intervals around these reference values in the
numerical results of Sec. IV.

Since the interaction range of the scalar field is inversely proportional to
the effective mass, Eq. (23) means that the chameleon fifth force will be
short-range in dense environments like the one under consideration and will be
acting as a long-range force on cosmological scales. For instance, the
effective mass of the chameleon inside the star will be much higher than the
chameleon mass at cosmological scales, that is $m_{{\rm eff},0}\gg m_{{\rm
eff},\infty}$. This is the key to the screening mechanism.

Qualitatively, one can distinguish, two different screening regimes according
to the behaviour of the field inside the star [13]. In the so-called thin-
shell regime, the chameleon field remains approximately constant within the
star, changing only in a very thin region close to the stellar radius. On the
contrary, in the alternative thick-shell regime, the scalar field evolves
right from the very centre of the star. In the latter situation, the solution
for the scalar field can be approximately written as

$\displaystyle\phi(r)\approx$
$\displaystyle\frac{\beta\rho_{c}r^{2}}{6M_{P}}+\bar{\phi}_{0},\hskip
81.09052pt0<r<R,$ (24) $\displaystyle\phi(r)\approx$
$\displaystyle-\frac{\beta}{4\pi M_{P}}\frac{Me^{-m_{{\rm
eff},\infty}(r-R)}}{r}+\bar{\phi}_{\infty},\quad r>R.$ (25)

Note that, while capturing the essence of the thick-shell regime, these
analytical expressions are based on three assumptions, which are not
guaranteed to be satisfied in the problem under consideration. First, the star
of mass $M$ and radius $R$ is taken to have a homogeneous density $\rho_{c}$.
Second, the contribution of the potential $V$ is assumed to be negligible as
compared to that of the coupling $A$ inside the star. Last, the scalar field
gradient is required to be large enough as compared to the curvature of the
potential outside the star.

###  III.2 Boundary Conditions

Note that we do not need to specify boundary conditions for the gravitational
potential $\Phi(r)$ – as long as we are in an equilibrium configuration –
since the system of equations (15)-(19) depends only on its radial
derivatives. For the pressure, we set $\tilde{p}(0)=\tilde{p}_{0}$, where
$\tilde{p}_{0}$ is the pressure at the WD centre. In our numerical
integration, this value will cover a range of pressures belonging to the EoS
validity domain. Actually, we will consider a range of central densities and
calculate the corresponding central pressures through the EoS introduced in
Sec. II.2 (Eqs. (12) and (13)). The central densities we will employ are the
ones we discussed in Sec. II.4, namely $\tilde{\rho}_{0}=7\times
10^{4}-10^{10}\text{ g}\,\text{cm}^{-3}$, since we want to compare our results
with the observations. Regarding the mass, one should set $m(0)=0$. However,
since our code starts from a certain initial radius $r=r_{0}>0$, at which we
consider the density to be $\tilde{\rho}_{0}$, the initial condition of the
mass is $m(r_{0})=(4/3)\pi r_{0}^{3}\tilde{\rho}_{0}$.

We assume that the universe is permeated by dark energy, meaning
$\tilde{p}=-\tilde{\epsilon}$ outside the star. This condition is not only
necessary because we have required spacetime to become Schwarzschild-de Sitter
far from the star, but also because we need a background density outside the
star for the chameleon to achieve an effective potential minimum at infinity.
Therefore, the stellar radius $R$ is determined by the condition
$\tilde{p}(R)=\tilde{p}_{\infty}=-\tilde{\epsilon}_{\infty}$, although
numerically this will translate into $\tilde{p}(R)$ being close to $0$ given a
certain tolerance since $\tilde{\epsilon}_{\infty}$ is very small. In
practice, we integrate up to a certain distance, which is big enough compared
to the typical WDs radii that it can be thought of as infinite, and then we
look for the radial coordinate at which the pressure is smaller than the
tolerance. That coordinate is the stellar radius $R$, and the stellar mass $M$
is defined as the total mass within $R$ following Eq. (16).

Since the solution for the scalar field must be regular at the centre of the
star, we know that the scalar field gradient fulfils $\sigma(0)=0$. We cannot
know the value of the scalar field at the centre before the integration, but
we do know that at infinity it will reach the exterior minimum, that is
$\phi\rightarrow\bar{\phi}_{\infty}$ as $r\rightarrow\infty$. This condition
implies that the solution is also regular at infinity, i.e. $\sigma\rightarrow
0$ as $r\rightarrow\infty$. Thus, to solve this ODE system, we ought to
implement a shooting method, so we can find out the adequate value of $\phi$
at $r=0$ that leads to $\bar{\phi}_{\infty}$ when we are away from the star.

###  III.3 Shooting Method

The scalar field $\phi$ is governed by a second-order differential equation,
Eq. (5), which we have split into two first-order differential equations, Eqs.
(18) and (19). As explained in Sec. III.2, we know both boundary conditions
for the scalar field gradient $\sigma$, but we only know the infinity one for
$\phi$, leaving the one at the origin to find.

We can estimate the value of the scalar field at the centre of the WD – let us
call it $\phi_{0}$ – from Eq. (21). This is possible because WDs are Newtonian
astrophysical objects and, as we explained in Sec. III.1, the chameleon
effective potential has a minimum inside such a star. Thus, we take it as a
sensible guess. We then perform the numerical integration of the ODE system,
either the relativistic (Eqs. (35)-(39)) or the Newtonian (Eqs. (15)-(19))
one. Afterwards, we compute the relative error between the scalar field
minimum at infinity, $\bar{\phi}_{\infty}$, and the value provided by our code
provides, $\phi(r_{\text{max}})$, where $r_{\text{max}}$ is the maximum radial
coordinate. If the relative error is smaller than a given tolerance – that is,
if
$|\bar{\phi}_{\infty}-\phi(r_{\text{max}})|/\bar{\phi}_{\infty}<\phi_{\text{tol}}$
– we have achieved convergence and, consequently, we store the output.

If the tolerance criterion is not met, we increase or decrease $\phi_{0}$ by a
small amount $\delta\phi$ depending on whether the difference between the
theoretical and the computed value – that is
$\bar{\phi}_{\infty}-\phi(r_{\text{max}})$ – is positive or negative. At every
step of the shooting method, we check the sign of the mentioned difference
and, whenever it changes (indicating that we have gone beyond the desired
value), we reduce $\delta\phi$. In this way, we boost convergence and achieve
higher precision in fewer steps.

##  IV Results

In this section, we present the results we have obtained from the numerical
integration of Eqs. (15)-(19) with the model functions $V(\phi)$ and $A(\phi)$
from (20), considering the EoS given by Eqs. (12) and (13), and using the
shooting method detailed in Sec. III.3. We have considered central densities
ranging from $\tilde{\rho}_{0,\text{min}}=7\times 10^{4}\text{
g}\,\text{cm}^{-3}$ to $\tilde{\rho}_{0,\text{max}}=10^{10}\text{
g}\,\text{cm}^{-3}$ and a background density
$\tilde{\rho}_{\infty}=10^{-4}\tilde{\rho}_{0,\text{min}}$. For each
$\tilde{\rho}_{0}$ value, we obtain the radius $R$ and the mass $M$ of the WD,
as we discussed in Sec. II.3.

We consider $n=1$, coupling strengths $\beta=0.1,\,0.05,\,0.01$, and energy
scales between $\Lambda=10^{-19}\,M_{P}$ and $\Lambda=1.5\times
10^{-18}\,M_{P}$. It should be mentioned that we do not study all possible
combinations of the values for the three parameters since preliminary runs
showed that some $\\{n,\beta,\Lambda\\}$ sets were incompatible with
astronomical observations, as we will discuss in Sec. IV.3. Let it be noted
that $\Lambda\lesssim 10^{-30}\,M_{P}$ should be reached to satisfy
equivalence principle constraints [13], but reaching such small values is
numerically very expensive. In our shooting method, we set
$\phi_{\text{tol}}=10^{-10}$, and to achieve such precision we already need to
work with a considerable amount of decimal places, even for $\Lambda\sim
10^{-18}\,M_{P}$. Nonetheless, our main conclusions would also stand for
realistic values of $\Lambda$, as we shall discuss.

###  IV.1 Stellar Structure

Fig. 2 shows the pressure radial profiles $\tilde{P}(r)$ for WDs in two
different realisations of the chameleon model characterised by a different
coupling strength, namely $\beta=0.1$ and $\beta=0.01$, while fixing the other
two parameters to $n=1$ and $\Lambda=10^{-18}\,M_{P}$. As for a star in GR or
Newtonian gravity, the pressure $\tilde{P}$ decreases along the radius. More
precisely, it drops when $r\sim 10^{3}-10^{4}$ km, which is the range of WDs
radii. We observe that the lower the coupling strength $\beta$, the longer it
takes for the pressure to decrease. From a mathematical point of view, this is
easily understood from Eq. (17). For our choice of chameleon functions, the
term $A_{,\phi}/A$ is simply $\beta/M_{P}$. Hence, the rate at which the
pressure diminishes is directly proportional to the product of $\beta$ and
$\sigma$.

Figure 2: Pressure radial profiles of chameleon-screened WDs. We consider the
chameleon model (20) with $n=1$, $\Lambda=10^{-18}\,M_{P}$, and $\beta=0.1$
(solid) and $\beta=0.01$ (dashed). Colours from purple to green indicate
increasing central densities, specifically for $\tilde{\rho}_{0}=7.0\times
10^{4},\,8.5\times 10^{5},\,1.0\times 10^{7},\,1.2\times 10^{8},\,1.5\times
10^{9}\text{ g}\,\text{cm}^{-3}$. The lower the coupling strength $\beta$, the
longer it takes for the pressure to decrease.

From a physical perspective, it is obvious that the pressure decrease will be
less affected by the scalar field if the coupling between the latter and the
matter is weaker, which is defined here as positive. Regarding the scalar
field gradient $\sigma$, we know it will also be positive. The scalar field
has a minimum inside the star and another one outside of it, the latter being
higher than the former since the outside density is lower than the inside one
(recall Eq. (21)). Plus, since there are no other extrema between these two
minima, the scalar field always increases. So, $\beta\sigma$ will always be
positive (see Figs. 3 and 4 for computational evidence). Then, since the
hydrostatic equation (17) has a global negative sign, the scalar contribution
to it will always boost the pressure decrease.

As we can already imagine, this pressure drop will cause the chameleon-
screened WDs to have different masses and radii than the ones in GR or
Newtonian gravity. Since the pressure falls earlier, we achieve the condition
$\tilde{p}(R)=-\tilde{\epsilon}_{\infty}$ sooner, thus we get a smaller
stellar radius $R$ in chameleon screening than in GR. This reduction is
translated also to the stellar mass since it is defined as the mass contained
in $R$. Accordingly, we obtain less massive stars in our chameleon model (see
Sec. IV.3).

###  IV.2 Scalar Profiles

Figure 3: Top panels: Scalar field radial profiles. Notice that the labels
indicate the difference between the scalar field and a reference value since
the variation is small compared to the latter. Bottom panels: Scalar field
gradient radial profiles.
The radial coordinate has been normalised to the respective stellar radius $R$
for each curve. We consider the chameleon model (20) with
$\Lambda=10^{-19}\,M_{P}$, $\beta=0.05$, and $n=1$ (left panels, green tones)
and $n=2$ (right panels, orange tones). Tones from dark to bright indicate
increasing central densities from $\tilde{\rho}_{0,\text{min}}=7\times
10^{4}\text{ g}\,\text{cm}^{-3}$ to $\tilde{\rho}_{0,\text{max}}=10^{10}\text{
g}\,\text{cm}^{-3}$.

Figure 4: Top panels: Scalar field radial profiles. Notice that the labels
indicate the difference between the scalar field and a reference value since
the variation is small compared to the latter. Bottom panels: Scalar field
gradient radial profiles.
The radial coordinate has been normalised to the respective stellar radius $R$
for each curve. We consider the chameleon model (20) with
$\Lambda=10^{-19}\,M_{P}$, $\beta=0.01$, and $n=1$ (left panels, green tones)
and $n=2$ (right panels, orange tones). Tones from dark to bright indicate
increasing central densities from $\tilde{\rho}_{0,\text{min}}=7\times
10^{4}\text{ g}\,\text{cm}^{-3}$ to $\tilde{\rho}_{0,\text{max}}=10^{10}\text{
g}\,\text{cm}^{-3}$.

In Fig. 3, we display the radial profiles for the scalar field $\phi(r)$ and
its gradient $\sigma(r)$ for two different parameter choices for the chameleon
model (20), namely $n=1$ and $n=2$ for fixed $\beta=0.05$ and
$\Lambda=10^{-19}\,M_{P}$. In Fig. 4 we plot the same radial profiles but for
$\beta=0.01$. In all four cases, the scalar field profile is approximately
flat for low central densities (darker colours). As the central density
increases, the scalar field becomes suppressed inside the star, giving place
to the characteristic thin-shell pattern [13]. This is the essence of the
chameleon screening mechanism, as we explained in Sec. III.1.

It is worth mentioning the difference between scalar field values between
$n=1$ and $n=2$ scenarios. In both cases – see panels (a) and (b) in Figs. 3
and 4 –, the scalar field in the $n=2$ model is approximately two orders of
magnitude bigger than in the $n=1$ model. Yet, this means that the chameleon
potential $V(\phi)$ (recall (20)) has higher values in the $n=1$ model. For
instance, for the scalar field values in Fig. 3, we have that $V\sim
10^{-92}\,M_{P}^{4}$ in the $n=1$ model and that $V\sim 10^{-113}\,M_{P}^{4}$
for $n=2$. For Fig. 4, one has that $V\sim 10^{-92}\,M_{P}^{4}$ for $n=1$ and
$V\sim 10^{-113}\,M_{P}^{4}$ for $n=2$. It seems that, for a fixed pair of
$\beta$ and $\Lambda$, the scalar field potential is much more significant for
$n=1$ than for $n=2$. This could lead us to think that the chameleon screening
will affect much more the WD in the $n=1$ case. Still, one must not forget the
gradient contribution, which, as discussed in Sec. IV.1, plays a crucial role
in the stellar structure.

The radially normalised scalar field gradient – normalised in the sense that
we have multiplied each scalar field gradient profile $\sigma(r)$ by the
corresponding stellar radius $R$ – displays similar values between the $n=1$
and the $n=2$ cases for the two considered coupling strengths (cf. the bottom
panels of Figs. 3 and 4). We also encounter the same values in both
realisations if we change the chameleon energy scale $\Lambda$. This
coincidence is qualitatively explained by the approximate solution (24)
presented in Sec. III.1. In particular, from the top panels of Figs. 3 and 4,
we realise that the scalar field changes throughout the whole stellar profile,
being therefore in the thick-shell regime where the referred solution applies.
This contrasts with solutions for more compact objects, such as neutron stars,
which exhibit a thin-shell behaviour [23]. Deriving Eq. (24) with respect to
$r$,

$\sigma(r)=\frac{\beta\rho_{c}r}{3M_{P}}\,,$ (26)

and taking into account that $\sigma R$ reaches its maximum somewhere between
$R/4$ and $R/2$, we have

$\sigma(xR)=\frac{\beta\rho_{c}xR}{3M_{P}}=\frac{\beta xM}{M_{P}4\pi
R^{2}}=\frac{2\beta xM_{P}\Phi_{c}}{R}\,,$ (27)

with $x$ the corresponding fraction of $R$. Note that in deriving this
expression, we have replaced the average density with the stellar mass and
radius, $\rho_{c}\equiv 3M/(4\pi R^{3})$, and introduced the Newtonian
potential at the surface of the star, $\Phi_{c}\equiv M/(8\pi M_{P}^{2}R)$.
Hence, in this approximation, the maximum of $\sigma R$ depends only on the
coupling strength $\beta$ and the Newtonian potential $\Phi_{c}$, that is

$\sigma(xR)R=2xM_{P}\beta\Phi_{c}\,.$ (28)

Due to this result, which as explained in Sec. III.1 is essentially based on
neglecting the contribution of the potential as compared to the chameleon
coupling function, we expected no dependence on the potential parameters $n$
and $\Lambda$. What we did not expect is that the numerical results would also
be independent of such parameters. It must be said though that Eq. (28) gives
values an order of magnitude below the ones we have obtained computationally.
Still, the ratio between both values of the $\sigma R$ maxima is consistent
through all the realisations of the chameleon model that we have studied. We
find this agreement between our results and the analytical solution
noteworthy.

###  IV.3 Mass-Radius Relation

Figure 5: Theoretical (solid) and parametric (dashed) MR curves for chameleon-
screened WDs for $n=1$ and different values of $\beta$ and $\Lambda$, the
latter changing in $10^{-19}M_{P}$ steps. We also include the WD MR curve in
Newtonian gravity (black) and the astrophysical data discussed in Sec. II.4.
Left panel: $\beta=0.1$, $\Lambda=1-10\times 10^{-19}M_{P}$ (green to blue).
Left panel: $\beta=0.05$, $\Lambda=1-12\times 10^{-19}M_{P}$ (yellow to
green).

Fig. 5 contains numerically computed MR curves for chameleon realisations
defined by the parameter values discussed at the beginning of this section,
with $n=1$. As anticipated from the pressure and mass profiles in Fig. 2, all
MR curves of chameleon-screened WDs are below (or practically on top of) the
MR curve predicted by Newtonian gravity. As expected, the smaller the coupling
parameter $\beta$, the closer the MR curves for our ST theory are between them
and to the Newtonian one.

Let us now compare our numerical results with the observational data in Sec.
II.4. To this end, we have included the data points coming from Gaia DR2 in
Fig. 5 so that one can rule out the parameter values excluded by observations.
We see that some of the computed MR curves are excluded from the currently
available data, even taking into account the error bars. Particularly, the
energy scale must be $\Lambda\lesssim 8\times 10^{-19}\,M_{P}$ if $\beta=0.1$,
$\Lambda\lesssim 10^{-18}\,M_{P}$ if $\beta=0.05$, and $\Lambda\lesssim
1.4\times 10^{-18}\,M_{P}$ if $\beta=0.01$ 222Nevertheless, these bounds are
irrelevant against cosmological constraints if we want the chameleon field to
play the role of dark energy.. However, a considerable amount of them are not.
We observe that, for each value of $\beta$ explored, the MR curves tend to
come together as $\Lambda$ decreases, as if they were reaching an asymptotic
curve that does not meet the Newtonian one. This leads us to think that, if we
were to consider much lower values of $\Lambda$, the MR curve for chameleon-
screened WDs would never meet that for Newtonian WDs. Unfortunately, we need
to significantly increase the numerical precision to explore such small energy
scales. Yet, the precision of the currently available data is not high enough
to be useful in that case.

Since the degeneracy between curves is evident – in the sense that we could
get the same masses and radii with various combinations of $n$, $\beta$, and
$\Lambda$ – it is useful to obtain a parametric function for them. A relation
like $R=R(M,n,\beta,\Lambda)$ can bring order into chaos, allowing not only to
compare different chameleon realisations between themselves but also between
other models, such as other ST theories with screening mechanisms. Since all
the solid curves in Fig. 5 appear to have the same shape, we assume this
formula for the mass-radius relation [38]

$R=R_{*}\left(\frac{M}{M_{\odot}}\right)^{-\frac{1}{3}}\sqrt{1-\left(\frac{M}{M_{*}}\right)^{\frac{4}{3}}},$
(29)

which perfectly fits the MR curve for WDs in Newtonian gravity. In that case,
we have that $R_{*}=8.83\times 10^{3}$ km and the $M_{*}$ parameter naturally
coincides with the Chandrasekhar mass, i.e. $M_{*}=M_{Ch}=1.45$ $M_{\odot}$.
To parameterise the MR curve of the chameleon-screened WDs, we turn the
constants $R_{*}$ and $M_{*}$ into functions of the parameters of our
chameleon model. After having explored the dependence of maximum radii and
maximum masses with $\beta$ and $\Lambda$, we found that parabolic functions
are a good ansatz, namely

$\displaystyle R_{*}$
$\displaystyle=R_{0}+A_{1}\beta+A_{2}\beta^{2}+B_{1}\frac{\Lambda}{\Lambda_{0}}+B_{2}\left(\frac{\Lambda}{\Lambda_{0}}\right)^{2}\,,$
(30) $\displaystyle M_{*}$
$\displaystyle=M_{0}+C_{1}\beta+C_{2}\beta^{2}+D_{1}\frac{\Lambda}{\Lambda_{0}}+D_{2}\left(\frac{\Lambda}{\Lambda_{0}}\right)^{2}\,,$
(31)

with $\Lambda_{0}=10^{-19}\,M_{P}$ a normalisation factor and $R_{0}$,
$M_{0}$, $A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$, $C_{1}$, $C_{2}$, $D_{1}$,
$D_{2}$ constants. The best-fit values for these parameters are summarized in
Tab. 1, with a coefficient of determination $R^{2}=0.992$. The corresponding
curves for values compatible with observational data are shown in Fig. 5, with
the same colouring as their numerical counterparts. As apparent in these
plots, the agreement between numerics and fitting formulae is enough for all
practical purposes.

Parameter | Value (km) | Parameter | Value ($M_{\odot}$)
---|---|---|---
$R_{0}$ | $9.638\times 10^{3}$ | $M_{0}$ | $1.532$
$A_{1}$ | $-2.005\times 10^{4}$ | $C_{1}$ | $-1.919$
$A_{2}$ | $1.066\times 10^{5}$ | $C_{2}$ | $7.528$
$B_{1}$ | $-1.434\times 10^{2}$ | $D_{1}$ | $-1.179\times 10^{-2}$
$B_{2}$ | $-5.941$ | $D_{2}$ | $-1.080\times 10^{-3}$
Table 1: Best-fit values for the parabolic parameterization of MR curves for
chameleon-screened WDs with $n=1$, defined by Eqs. (29), (30), and (31), with
a $R^{2}=0.992$ coefficient of determination. Note that $\beta\ll 1$ results
in the larger values for the $A_{i}$ and $C_{i}$ coefficients.

##  V Conclusion

In this work, we have studied the effect that chameleon screening can have on
the structure of WDs. They are auspicious candidates to test alternative
theories of gravity since there is extensive observational data available and
we have a sensible grasp of the EoS describing the matter within them. We have
shown that the Newtonian approximation accurately describes WDs in this
particular kind of ST theory, as it happens in GR. We have considered a
Chandrasekhar EoS and solved the equilibrium equations with a shooting method
of our design.

We have seen that the presence of the chameleon field significantly affects
the WD’s internal pressure, causing it to drop prematurely compared to when
there is no scalar field. This leads to smaller stellar masses and radii
which, in turn, shifts the MR curves below the MR relation predicted by
Newtonian gravity. However, those stars above the theoretical curve – which
are the majority – cannot come from any chameleon realisation, being therefore
necessary to invoke other mechanisms such as strong magnetic fields [39, 40]
or inverse chameleon settings [41] to explain them.

We have examined the radial profiles of the scalar field and its gradient,
considering a wide range of values for the model parameters – the energy scale
$\Lambda$ and the coupling strength $\beta$ – and showing, for the first time
in the literature, $n=2$ realisations of the chameleon. This allowed us to
identify a similarity relation for the radially normalised scalar field
gradient, that is $\sigma R$. In the thick-shell regime of the chameleon
screening, the $\sigma R$ maximum is determined by $\beta$, and we have
encountered the same behaviour in our numerical results, independently of the
$\Lambda$ and $n$ values considered.

After exploring the chameleon energy scales and coupling strengths compatible
with astronomical observations, we have inferred parametric expressions for
the MR relations that depend on the mentioned parameters. This result allows
us to analyse the compatibility of our model with current and future
observations without having to numerically solve the ODE systems again, a
computationally time-consuming task. Moreover, this allows us to check for
degeneracies between other classes of screening mechanisms once the analogous
formulae are derived.

A continuation of this work could be to explore the stability of chameleon-
screened WDs through radial perturbations, even studying the stable regime of
WD oscillations. Furthermore, the effect that such a scalar field has on other
properties of WDs, such as their cooling time, remains unexplored. Naturally,
an interesting prospect would be to apply the presented framework and the
developed computational tools to other ST theories.

###### Acknowledgements.

We are grateful to Raissa F. P. Mendes for detailed discussions on the
numerical methods employed in this work. The numerical part of this work has
been performed with the support of the Infraestrutura Nacional de Computação
Distribuída (INCD), funded by the Fundação para a Ciência e a Tecnologia (FCT)
and FEDER under the project 01/SAICT/2016 nº 022153. JBE (ORCID
0000-0003-4121-3179) and IL (ORCID 0000-0002-5011-9195) acknowledge the FCT,
Portugal, for the financial support to the Center for Astrophysics and
Gravitation - CENTRA, Instituto Superior Técnico, Universidade de Lisboa,
through Project No. UIDB/00099/2020. JBE is grateful for the support of this
agency through grant No. SFRH/BD/150989/2021 in the framework of the IDPASC-
Portugal Doctoral Program. IL also acknowledges FCT for the financial support
through grant No. PTDC/FIS-AST/28920/2017. JR (ORCID 0000-0001-7545-1533) is
supported by a Ramón y Cajal contract of the Spanish Ministry of Science and
Innovation with Ref. RYC2020-028870-I. This work was supported by the project
PID2022-139841NB-I00 of MICIU/AEI/10.13039/501100011033 and FEDER, UE.

##  Appendix A Relativistic and Newtonian Descriptions

In this appendix, we compare the accuracy of the Newtonian description against
the purely relativistic one for a static, spherically symmetric WD described
by a perfect fluid energy-momentum tensor. We have already discussed the
Newtonian approximation in Sec. III. For the relativistic description, we
adopt the following line element to describe the spacetime

$ds^{2}=-e^{2\nu(r)}dt^{2}+e^{2\lambda(r)}dr^{2}+r^{2}d\Omega^{2}\,,$ (32)

where $d\Omega^{2}=d\theta^{2}+\text{sin}^{2}\theta d\varphi^{2}$. Far away
from the star, the spacetime must become Schwarzschild-de Sitter, that is

$ds^{2}=-f(r)dt^{2}+f(r)^{-1}dr^{2}+r^{2}d\Omega^{2},$ (33)

where $f(r)=1-2a/r-br^{2}$. For the present case, we require that $a=m$ and
$b=0$ [23], which means that

$\lambda(r)\rightarrow-\frac{1}{2}\text{ln}\left(1-\frac{2m(r)}{r}\right)$
(34)

for $r\gg R$, with $R$ being the stellar radius. Then, if we insert the metric
of Eq. (32) into Eqs. (2), (6), and (5), we obtain the
Tolman–Oppenheimer–Volkoff (TOV) equation with a scalar field contribution
333Taking Eq. (35) until the second term, replacing it in Eq. (37), and
setting $A(\phi)=1$, one recovers the well-known TOV equation.

$\displaystyle\nu^{\prime}=$
$\displaystyle\frac{1}{r-2m}\bigg{\\{}\frac{m}{r}$ (35)
$\displaystyle+\frac{\kappa^{2}}{2}r^{2}\left[A^{4}\tilde{p}+\left(\frac{1}{2}-\frac{m}{r}\right)\sigma^{2}-V\right]\bigg{\\}},$
$\displaystyle m^{\prime}=$
$\displaystyle\frac{\kappa^{2}}{2}r^{2}\left[A^{4}\tilde{\epsilon}+\left(\frac{1}{2}-\frac{m}{r}\right)\sigma^{2}+V\right],$
(36) $\displaystyle\tilde{p}^{\prime}=$
$\displaystyle-(\tilde{p}+\tilde{\epsilon})\left(\nu^{\prime}+\frac{A_{,\phi}}{A}\sigma\right),$
(37) $\displaystyle\phi^{\prime}=$ $\displaystyle\sigma,$ (38)
$\displaystyle\sigma^{\prime}=$
$\displaystyle\left[\frac{1}{r-2m}\left(m^{\prime}+3\frac{m}{r}-2\right)-\nu^{\prime}\right]\sigma$
(39)
$\displaystyle+\frac{r}{r-2m}\left[V_{,\phi}-A_{,\phi}A^{3}(3\tilde{p}-\tilde{\epsilon})\right],$

where primes denote derivatives with respect to $r$. Once one chooses the
model functions $V(\phi)$ and $A(\phi)$, and a suitable EoS, this ODE system
can be numerically integrated, as it happened for Eqs. (15)-(19). To obtain a
MR curve, we follow the same procedure we explained in the main body of the
work (see Sec. III.2).

In Fig. 6, we display the MR curve for several realisations of the chameleon
model as well as for GR. Specifically, we have integrated Eqs. (35)-(39) and
Eqs. (15)-(19) for $n=1$, $\Lambda=1.2\times 10^{-18}\,M_{P}$, and different
values of $\beta$. We have considered $\Lambda=\beta=0$ too, which effectively
reduces the two ODE systems to the TOV equation and its Newtonian counterpart,
respectively. As discussed in Sec. II.2, WDs are non-relativistic objects, and
it is well-known that a Newtonian approximation is more than sufficient to
describe them. Consequently, we were expecting what we see from the two solid
curves in Fig. 6, which show us that the masses and radii predicted by the TOV
equation (in blue) coincide with those calculated in the Newtonian limit (in
red).

We observe a small discrepancy just in the more massive end of the curve,
close to the Chandrasekhar limit of $1.4\,M_{\odot}$, which is slightly
surpassed in the Newtonian limit. One can say that Newtonian physics
overestimates radii and underestimates surface gravity, thus demonstrating the
significance of general relativistic effects in determining the physical
properties of these compact stars, only for particularly massive WDs. For
instance, in [43], they found that the radius predicted by GR for a WD with a
mass of 1.415 $M_{\odot}$ is approximately 33% smaller than that calculated in
Newtonian physics. Still, it should be noted that they consider a different
EoS from ours and that WD so massive represent a very small portion of the
observed stars (recall Fig. 1).

Regarding the MR curves computed in our ST theory, we find the same behaviour:
the discrepancy between GR and Newtonian gravity manifests itself in the most
massive WDs. As a result, we have shown that the Newtonian description is
completely sufficient for chameleon-screened WDs, thus it was no blunder from
our side to only consider the latter in the main body of the paper.
Fortunately, this boosted our work from a computational perspective since Eqs.
(15)-(19) are much simpler than Eqs. (35)-(39), hence another reason in favour
of having ignored the modified TOV equation.

Figure 6: Theoretical MR curves for WDs in GR (solid) and in the presence of
chameleon field for $n=1$, $\Lambda=1.2\times 10^{-18}\,M_{P}$ and $\beta=0.1$
(dotted), $\beta=0.05$ (dashed), and $\beta=0.01$ (dash-dotted), both in the
relativistic (TOV) approach (blue) and the Newtonian (NWT) approximation
(red).

## References

  * Will [2018] C. M. Will,  _Theory and Experiment in Gravitational Physics_ (Cambridge University Press, 2018).
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified Gravity and Cosmology, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Burrage and Sakstein [2018] C. Burrage and J. Sakstein, Tests of chameleon gravity, Living Reviews in Relativity 21, 10.1007/s41114-018-0011-x (2018).
  * Brax  _et al._ [2021] P. Brax, S. Casas, H. Desmond, and B. Elder, Testing Screened Modified Gravity, Universe 8, 11 (2021), arXiv:2201.10817 [gr-qc] .
  * Fischer  _et al._ [2024] H. Fischer, C. Käding, and M. Pitschmann, Screened Scalar Fields in the Laboratory and the Solar System,  , arXiv:2405.14638 (2024).
  * Amendola  _et al._ [2018] L. Amendola, J. Rubio, and C. Wetterich, Primordial black holes from fifth forces, Phys. Rev. D 97, 081302 (2018), arXiv:1711.09915 [astro-ph.CO] .
  * Savastano  _et al._ [2019] S. Savastano, L. Amendola, J. Rubio, and C. Wetterich, Primordial dark matter halos from fifth forces, Phys. Rev. D 100, 083518 (2019), arXiv:1906.05300 [astro-ph.CO] .
  * Goh  _et al._ [2024] L. W. K. Goh, J. Bachs-Esteban, A. Gómez-Valent, V. Pettorino, and J. Rubio, Observational constraints on early coupled quintessence, Phys. Rev. D 109, 023530 (2024), arXiv:2308.06406 [astro-ph.CO] .
  * Wetterich [1995] C. Wetterich, The Cosmon model for an asymptotically vanishing time dependent cosmological ’constant’, Astron. Astrophys. 301, 321 (1995), arXiv:hep-th/9408025 .
  * Wetterich [2007] C. Wetterich, Growing neutrinos and cosmological selection, Phys. Lett. B 655, 201 (2007), arXiv:0706.4427 [hep-ph] .
  * Amendola  _et al._ [2008] L. Amendola, M. Baldi, and C. Wetterich, Quintessence cosmologies with a growing matter component, Phys. Rev. D 78, 023015 (2008), arXiv:0706.3064 [astro-ph] .
  * Casas  _et al._ [2016] S. Casas, V. Pettorino, and C. Wetterich, Dynamics of neutrino lumps in growing neutrino quintessence, Phys. Rev. D 94, 103518 (2016), arXiv:1608.02358 [astro-ph.CO] .
  * Khoury and Weltman [2004] J. Khoury and A. Weltman, Chameleon cosmology, Physical Review D 69, 10.1103/physrevd.69.044026 (2004).
  * Hinterbichler  _et al._ [2011] K. Hinterbichler, J. Khoury, A. Levy, and A. Matas, Symmetron Cosmology, Phys. Rev. D 84, 103521 (2011), arXiv:1107.2112 [astro-ph.CO] .
  * Brax  _et al._ [2010] P. Brax, C. van de Bruck, A.-C. Davis, and D. Shaw, The Dilaton and Modified Gravity, Phys. Rev. D 82, 063519 (2010), arXiv:1005.3735 [astro-ph.CO] .
  * Vainshtein [1972] A. I. Vainshtein, To the problem of nonvanishing gravitation mass, Phys. Lett. B 39, 393 (1972).
  * Babichev  _et al._ [2009] E. Babichev, C. Deffayet, and R. Ziour, k-Mouflage gravity, Int. J. Mod. Phys. D 18, 2147 (2009), arXiv:0905.2943 [hep-th] .
  * Babichev and Langlois [2010] E. Babichev and D. Langlois, Relativistic stars in f(R) and scalar-tensor theories, Phys. Rev. D 81, 124051 (2010), arXiv:0911.1297 [gr-qc] .
  * Chang and Hui [2011] P. Chang and L. Hui, Stellar Structure and Tests of Modified Gravity, Astrophys. J. 732, 25 (2011), arXiv:1011.4107 [astro-ph.CO] .
  * Sakstein [2013] J. Sakstein, Stellar oscillations in modified gravity, Physical Review D 88, 10.1103/physrevd.88.124013 (2013).
  * Brito  _et al._ [2014] R. Brito, A. Terrana, M. Johnson, and V. Cardoso, Nonlinear dynamical stability of infrared modifications of gravity, Phys. Rev. D 90, 124035 (2014), arXiv:1409.0886 [hep-th] .
  * Babichev  _et al._ [2016] E. Babichev, K. Koyama, D. Langlois, R. Saito, and J. Sakstein, Relativistic Stars in Beyond Horndeski Theories, Class. Quant. Grav. 33, 235014 (2016), arXiv:1606.06627 [gr-qc] .
  * de Aguiar and Mendes [2020] B. F. de Aguiar and R. F. Mendes, Highly compact neutron stars and screening mechanisms: Equilibrium and stability, Physical Review D 102, 10.1103/physrevd.102.024064 (2020).
  * de Aguiar  _et al._ [2021] B. F. de Aguiar, R. F. P. Mendes, and F. T. Falciano, Neutron Stars in the Symmetron Model, Universe 8, 6 (2021), arXiv:2112.03823 [gr-qc] .
  * Panotopoulos  _et al._ [2023] G. Panotopoulos, J. Rubio, and I. Lopes, On the impact of nonlocal gravity on compact stars, Int. J. Mod. Phys. D 32, 2250139 (2023), arXiv:2106.10582 [gr-qc] .
  * ter Haar  _et al._ [2021] L. ter Haar, M. Bezares, M. Crisostomi, E. Barausse, and C. Palenzuela, Dynamics of Screening in Modified Gravity, Phys. Rev. Lett. 126, 091102 (2021), arXiv:2009.03354 [gr-qc] .
  * Dima  _et al._ [2021] A. Dima, M. Bezares, and E. Barausse, Dynamical chameleon neutron stars: Stability, radial oscillations, and scalar radiation in spherical symmetry, Physical Review D 104, 10.1103/physrevd.104.084017 (2021).
  * Saltas  _et al._ [2018] I. D. Saltas, I. Sawicki, and I. Lopes, White dwarfs and revelations, Journal of Cosmology and Astroparticle Physics 2018 (05), 028–028.
  * Alam and Islam [2023] K. Alam and T. Islam, White dwarf mass-radius relation in theories beyond general relativity, JCAP 08, 081, arXiv:2301.08677 [gr-qc] .
  * Kalita and Uniyal [2023] S. Kalita and A. Uniyal, Constraining fundamental parameters in modified gravity using gaia-dr2 massive white dwarf observations, The Astrophysical Journal 949, 62 (2023).
  * Camenzind [2007] M. Camenzind,  _Compact Objects in Astrophysics. White Dwarfs, Neutron Stars and Black Holes_ (Springer, Berlin Heidelberg, 2007).
  * Jiménez-Esteban  _et al._ [2018] F. M. Jiménez-Esteban, S. Torres, A. Rebassa-Mansergas, G. Skorobogatov, E. Solano, C. Cantero, and C. Rodrigo, A white dwarf catalogue from gaia-dr2 and the virtual observatory, Monthly Notices of the Royal Astronomical Society 480, 4505–4518 (2018).
  * Tremblay  _et al._ [2018] P.-E. Tremblay, E. Cukanovaite, N. P. Gentile Fusillo, T. Cunningham, and M. A. Hollands, Fundamental parameter accuracy of DA and DB white dwarfs in Gaia Data Release 2, Monthly Notices of the Royal Astronomical Society 482, 5222 (2018).
  * Kilic  _et al._ [2020] M. Kilic, P. Bergeron, A. Kosakowski, W. R. Brown, M. A. Agüeros, and S. Blouin, The 100 pc white dwarf sample in the sdss footprint, The Astrophysical Journal 898, 84 (2020).
  * Note [1] The first EoS for such stars was derived by Chandrasekhar [44]. Hamada and Salpeter added temperature corrections to it [45].
  * Misner  _et al._ [1973] C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, San Francisco, 1973).
  * Note [2] Nevertheless, these bounds are irrelevant against cosmological constraints if we want the chameleon field to play the role of dark energy.
  * Walter Greiner [1995] H. S. Walter Greiner, Ludwig Neise,  _Thermodynamics and Statistical Mechanics_ (Springer, Berlin Heidelberg, 1995).
  * Das and Mukhopadhyay [2014] U. Das and B. Mukhopadhyay, Maximum mass of stable magnetized highly super-Chandrasekhar white dwarfs: stable solutions with varying magnetic fields, JCAP 06, 050, arXiv:1404.7627 [astro-ph.SR] .
  * Roy  _et al._ [2019] S. K. Roy, S. Mukhopadhyay, J. Lahiri, and D. N. Basu, Relativistic Thomas-Fermi equation of state for magnetized white dwarfs, Phys. Rev. D 100, 063008 (2019), arXiv:1907.13480 [nucl-th] .
  * Wei and Yu [2021] H. Wei and Z.-X. Yu, Inverse chameleon mechanism and mass limits for compact stars, JCAP 08, 011, arXiv:2103.12696 [gr-qc] .
  * Note [3] Taking Eq. (35) until the second term, replacing it in Eq. (37), and setting $A(\phi)=1$, one recovers the well-known TOV equation.
  * Carvalho  _et al._ [2018] G. A. Carvalho, R. M. Marinho, and M. Malheiro, General relativistic effects in the structure of massive white dwarfs, General Relativity and Gravitation 50, 10.1007/s10714-018-2354-8 (2018).
  * Chandrasekhar [1935] S. Chandrasekhar, The highly collapsed configurations of a stellar mass (Second paper), Mon. Not. Roy. Astron. Soc. 95, 207 (1935).
  * Hamada and Salpeter [1961] T. Hamada and E. E. Salpeter, Models for Zero-Temperature Stars., Astrophys. J. 134, 683 (1961).