# Accurate Computation of the Screening of Scalar Fifth Forces in Galaxies

Clare Burrage  Bradley March  and Aneesh P. Naik

###### Abstract

Screening mechanisms allow light scalar fields to dynamically avoid the
constraints that come from our lack of observation of a long-range fifth
force. Galactic scale tests are of particular interest when the light scalar
is introduced to explain the dark matter or dark energy that dominates our
cosmology. To date, much of the literature that has studied screening in
galaxies has described screening using simplifying approximations. In this
work, we calculate numerical solutions for scalar fields with screening
mechanisms in galactic contexts, and use these to derive new, precise
conditions governing where fifth forces are screened. We show that the
commonly used binary screened/unscreened threshold can predict a fifth force
signal in situations where a fuller treatment does not, leading us to conclude
that existing constraints might be significantly overestimated. We show that
various other approximations of the screening radius provide a more accurate
proxy to screening, although they fail to exactly reproduce the true screening
surface in certain regions of parameter space. As a demonstration of our
scheme, we apply it to an idealised Milky Way and thus identify the region of
parameter space in which the solar system is screened.

##  1 Introduction

Attempts to resolve long-standing complications within the dark sector often
introduce a new scalar field which naturally couples to the gravitational
sector, causing modifications to standard gravity, General Relativity (GR), in
the form of scalar fifth forces. However, such deviations from GR are tightly
constrained by lab-based experiments testing the inverse-square law and tests
of Solar System dynamics [1]. To avoid such constraints, scalar-tensor
theories with screening propose that the scalar field is environment-
dependent, allowing for modifications to GR on certain scales whilst
suppressing fifth forces on other scales [2, 3, 4, 5, 6, 7, 8].

In this paper, we study chameleon [9, 10] and symmetron [11, 12] screening
models (for related work on other scalar theories see refs. [13, 14, 15, 16,
17, 18, 19, 20]). Chameleon and symmetron theories introduce a density-
dependent scalar potential that screens the fifth force in areas of high
density or deep gravitational potential, such as in stellar interiors or the
centres of galaxies.

Much of the recent literature has studied the behaviour of these theories on
galactic scales, where the screening phenomenon leads to a rich phenomenology,
such as equivalence principle violations and other telltale kinematic
signatures [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37]. The recent interest in galactic scales is a result of two circumstances:
the fact that galaxies inhabit an unconstrained “desert” in modified gravity
parameter space [38], and the recent advent of vast datasets from various
galaxy surveys.

However, investigating these scalar fifth forces on galactic scales is a
challenging problem: it is difficult to isolate any potential fifth force
signal from the numerous astrophysical effects that can obscure it, and
forward-modelling the signal involves the solution of a non-linear equation of
motion (EoM). As a consequence, much of the literature has relied on
approximations to solve the equations of motion, assuming a simplistic setup
and dropping terms to simplify the non-linear dynamics. The most common
approximation is to simply assume that a galaxy, depending on its mass and
environment, is either fully screened or fully unscreened. This disregards the
possibility of partial screening, under which the central region of the galaxy
is screened while its outskirts remain unscreened. An example of a work using
this approximation is that of ref. [27], which uses it to derive the current
strongest constraints on chameleon $f(R)$ gravity, altogether ruling out the
astrophysically-relevant part of parameter space111In ref. [27] a Bayesian
analysis is used so that a galaxy is assigned a probability of being fully
screened, or fully unscreened. It is possible that this captures some of the
variation due to galaxies being partially screened, but we are unable to say
how much.. The purpose of the present work is to test these approximations,
and the possible impact that relaxing them might have on such constraints.

In this work, we investigate the non-linear nature of screening theories in
the galactic regime by utilising a numerical solver, based on refs. [39, 31,
40], to determine the true scalar field profile inside our static model
galaxies. Guided by these field solutions, we obtain empirical rules
describing the locations of screening surfaces and compare these results with
the results of the common approximations. We discuss the implications of the
proper non-linear treatment of these theories and how our findings can impact
observables and constraints established using simplified methods. These
results will form the foundation for a reevaluation of galaxy-based
constraints in the near future, using methods that we have thoroughly explored
and tested in this paper. This reanalysis will enable us to provide a
quantitative estimate of the impact of using standard approximations to
screening, compared to our precise numerical methods that consider partial
screening.

The outline for this paper is as follows. In section 2 we present an overview
of the general scalar-tensor theory, as well as our specific chameleon and
symmetron models. In section 3 we describe our methods to solve the theories
numerically, including the galactic model we use, the mechanics behind the
numerical solver, and the criteria used to define the screened and unscreened
regions. Section 4 outlines the results of our investigation. Our conclusions
are laid out in section 5.

##  2 Screened Modified Gravity Theories

###  2.1 General Scalar-Tensor Model

In the present study, our investigation into modified theories of gravity is
exclusively focused on scalar-tensor theories. These theories are described by
the Einstein frame action, given by

$S=\frac{c^{3}}{8\pi G}\int
d^{4}x\sqrt{-g}\left[\frac{R}{2}-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi-V(\phi)\right]+S_{m}\left[\tilde{g}_{\mu\nu},\psi^{SM}_{i}\right]\,,$
(2.1)

where $g$ is the determinant of the Einstein frame metric, $g_{\mu\nu}$, $R$
is the Ricci scalar, $G$ is the gravitational constant, $\phi$ is a
canonically normalised (dimensionless) scalar field with a potential
$V(\phi)$, and $S_{m}$ represents the action for the standard model fields,
$\psi_{i}^{SM}$. The theories also include a conformal coupling to matter,
meaning that matter fields move on geodesics of the Jordan frame metric

$\tilde{g}_{\mu\nu}=A^{2}(\phi)g_{\mu\nu}\,.$ (2.2)

Deviations from GR arise from the non-minimal coupling, $A(\phi)$. Although
both frames have identical classical observables, for ease of calculation,
this analysis focuses solely on the Einstein frame, where there is a direct
coupling between matter and the scalar field, and the pure gravitational
action does not depend on the scalar field.

Scalar-tensor theories of gravity modify gravitational forces by introducing
_fifth forces_ , which can be obtained by taking the Newtonian limit of the
geodesic equations (see for example ref. [5]) to arrive at the following
expression for the acceleration

$\frac{\bm{a}_{5}}{c^{2}}=-\beta(\phi)\bm{\nabla}\phi\,,$ (2.3)

where $\beta(\phi)\equiv\mathrm{d}\ln A/\mathrm{d\phi}$. Meanwhile, the scalar
field’s EoM is derived by extremising the general scalar-tensor action, eq.
(2.1), with an energy-momentum tensor for the matter fields that assumes
matter is a non-relativistic, perfect fluid with mass density $\rho$

$\Box\phi=\frac{dV_{\mathrm{eff}}}{d\phi}\equiv\frac{dV(\phi)}{d\phi}+\frac{8\pi
G\rho}{c^{2}}\beta(\phi)\,,$ (2.4)

where the effective potential, $V_{\mathrm{eff}}(\phi)$, is given by

$V_{\mathrm{eff}}(\phi)\equiv V(\phi)+\frac{8\pi G\rho}{c^{2}}\ln A(\phi)\,.$
(2.5)

For certain choices of $V(\phi)$ and $A(\phi)$, the effective potential has a
stable minimum and this density-dependent EoM admits screening solutions, i.e.
solutions in which the fifth force, eq. (2.3), is rendered negligible in
certain environments but unleashed elsewhere. Screening can arise via several
mechanisms: a small matter coupling, $\beta(\phi)$, a large scalar mass
resulting in short-ranged forces, or not all the mass sourcing (or coupling
to) the scalar field [5].

###  2.2 Astrophysical Screening Approximations

The crux of this paper is to determine the location of the screened region
inside galaxies, defined by what we term the ‘screening surface’: the
interface between the screened central region and the unscreened outskirts of
a partially screened object. To determine this, we use numerical solutions of
scalar EoM. These results will then be compared with the approximate methods
that are commonly used in the literature.

One prevalent method to approximate screening surfaces considers a static,
spherically-symmetric source, $\delta\rho(r)$, in a constant density
background, $\bar{\rho}$, e.g. a star embedded in the galactic medium or a
dark matter halo embedded in the cosmological background, as in refs. [41, 42,
43, 44]. For such an object, the screening surface will necessarily be a
sphere with radius $r_{s}$, the ‘screening radius’. The general EoM, eq.
(2.4), is then solved in a piecewise manner, described in full in appendix A,
to obtain the expression

$\displaystyle\chi\equiv\frac{\bar{\phi}}{2\beta(\bar{\phi})}$
$\displaystyle=\frac{4\pi G}{c^{2}}\int^{\infty}_{r_{s}}r\delta\rho(r)dr\,,$
(2.6)
$\displaystyle=-\frac{1}{c^{2}}\left(\Phi_{N}(r_{s})+r_{s}\Phi_{N}^{\prime}(r_{s})\right)\,,$
(2.7)

where $\bar{\phi}=\phi(\bar{\rho})$ is the background field solution and
$\Phi_{N}(r)$ is the Newtonian potential. To obtain this $\chi$ parameter,
hereby referred to as the critical potential, several assumptions are
required, specifically that the field reaches its minimum value inside the
object, that this minimum value can be neglected, that the theory can be
linearised in the region containing the screening radius, that the scalar mass
is negligible in this region, and that the coupling function $\beta(\phi)$ is
constant throughout.

If eq. (2.6), or equally eq. (2.7), has no solution for $r_{s}$ this
corresponds to a fully unscreened object, i.e. none of the object is screened.
Since $\Phi_{N}<0$ and $\Phi^{\prime}_{N}>0$, an object must be fully
unscreened when $\chi>|\Phi_{N}|/c^{2}$, where the Newtonian potential is
measured at the surface of the object. In reverse, this indicates only objects
with $\chi<|\Phi_{N}|/c^{2}$ are partially screened (or fully screened, i.e.
$r_{s}$ is larger than the extent of the object, in the case
$\chi\ll|\Phi_{N}|/c^{2}$).

The critical potential is a starting point for two approaches commonly seen in
the literature for approximating the screening of astrophysical objects. In
the first approach, the implicit equation (2.6) is solved for the screening
radius, which appears as the lower limit of the integral. As we shall see,
this can readily be done numerically, but analytic solutions can also be found
for certain choices of the density distribution (e.g., for idealised dark
matter haloes; as in ref. [45]).

The second approach, which we refer to as the ‘binary screening condition’,
disregards partial screening altogether and classifies objects as either fully
screened or fully unscreened by comparing $\chi$ to an estimate of $\Phi_{N}$
at the surface of the object. Neglecting the contribution of its environment,
an object is classed as screened if $|\Phi_{N}|/c^{2}>\chi$, and unscreened
otherwise. For example, at the level of current constraints in $f(R)$ gravity,
refs. [32, 24, 25, 27], objects with deep potentials, e.g. main sequence stars
with $|\Phi_{N}|/c^{2}\sim 10^{-6}$ are likely screened, while objects with
shallower potentials e.g. dwarf galaxies with $|\Phi_{N}|/c^{2}\sim 10^{-8}$
can still be unscreened. In the context of galaxies, this approach was
suggested by ref. [46], who proposed using
$\Phi_{\mathrm{vir}}=GM_{\mathrm{vir}}/R_{\mathrm{vir}}$, $M_{\mathrm{vir}}$
and $R_{\mathrm{vir}}$ will be defined in sec. 3.1. The same approach was then
subsequently used in a string of recent papers deriving modified gravity
constraints on galactic scales (refs. [47, 24, 25, 27]). However, it is worth
noting that these works also included environmental contributions, classifying
objects as screened when $(|\Phi_{N}|+|\Phi_{\mathrm{env.}}|)/c^{2}>\chi$,
where $|\Phi_{\mathrm{env.}}|$ is the environmental contribution to the local
gravitational potential.

One of the main aims of the present work is to compare the results obtained
from these two approximate approaches with more accurate results obtained from
exact numerical solutions of the scalar equations of motion.

###  2.3 Chameleon $f(R)$

Chameleon models are a class of scalar-tensor theories that screen primarily
via a density-dependent scalar mass, leading to the presence of short-range
fifth forces. Screened regions in chameleon models are characterised by flat
scalar field profiles, which result in the suppression of the fifth force, eq.
(2.3). We choose to work with $f(R)$ chameleon screening because of its
simplicity as a model and its ubiquity in the literature.

$f(R)$ theories modify GR by replacing the Ricci scalar, $R$, in the Einstein-
Hilbert action with the generalised $R+f(R)$, allowing for higher-order
curvature terms in the action (see ref. [48] for a general review),

$S=\frac{c^{3}}{16\pi G}\int
d^{4}x\sqrt{-\tilde{g}}\left[R+f(R)\right]+S_{m}\left[\tilde{g}_{\mu\nu},\psi^{SM}_{i}\right]\,,$
(2.8)

where $R=R(\tilde{g})$. Here, $f_{R}\equiv\frac{\mathrm{d}f(R)}{\mathrm{d}R}$
acts as the scalar field of the theory. This action can be recast to the
general scalar-tensor action eq. (2.1) with the field redefinition [49]

$\phi=-\frac{\ln(1+f_{R})}{2\beta(\phi)}\,,\qquad
V(\phi)=\frac{Rf_{R}-f(R)}{2(1+f_{R})^{2}}\qquad\mathrm{and}\qquad\beta(\phi)=\frac{1}{\sqrt{6}}\,.$
(2.9)

For specific choices of the functional form of $f(R)$, the theory can exhibit
the chameleon screening mechanism. Several forms of $f(R)$ can achieve this,
refs. [50, 51, 52], we adopt the commonly used Hu-Sawicki model, [53]. In this
model

$f(R)=-\frac{am^{2}}{1+(R/m^{2})^{-b}}\,,$ (2.10)

where $a$ and $b$ are positive dimensionless parameters, and $m$ is a
dimensionful parameter. For simplicity, we choose $b=1$, in line with the bulk
of the current literature. The two remaining free parameters $a$ and $m$ can
be related by requiring that gravity returns to GR + $\Lambda$CDM in the high
curvature limit (i.e. $f(R)\approx-2\Lambda$ when $R\gg m^{2}$), resulting in

$am^{2}=2\Lambda\,,$ (2.11)

where $\Lambda$ is the cosmological constant. This allows us to describe the
theory in terms of one free model parameter, $a$ or $m$, which we instead
express as the present-day background field value $f_{R0}$:

$a=-\frac{4\Omega_{\Lambda}^{2}}{(\Omega_{m}+4\Omega_{\Lambda})^{2}}\frac{1}{f_{R0}}\quad\mathrm{and}\quad
m^{2}=-\frac{3H_{0}^{2}(\Omega_{m}+4\Omega_{\Lambda})^{2}}{2\Omega_{\Lambda}c^{2}}f_{R0},$
(2.12)

where $\Omega_{m}$ and $\Omega_{\Lambda}$ are the cosmological density
parameters for matter and dark energy respectively. Throughout this article we
assume values of $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$.

Extremising the action eq. (2.8) with respect to the metric results in a set
of modified Einstein field equations. Taking the trace and then the Newtonian
limit, with the assumption $|f_{R}|\ll 1$ and the quasi-static approximation
$|\nabla f_{R}|\gg\partial_{t}f_{R}$, produces the $f(R)$ EoM,

$\nabla^{2}f_{R}=\frac{1}{3}\left(\delta R-\frac{8\pi
G}{c^{2}}\delta\rho\right)\,.$ (2.13)

Here, $\delta\rho$ is the density perturbation and $\delta R$ is the curvature
perturbation, given by

$\delta R=R_{0}\left[\sqrt{\frac{f_{R0}}{f_{R}}}-1\right]\,,$ (2.14)

under the Hu-Sawicki model, with zero subscripts denoting background values
today.

The fifth force acceleration, eq. (2.3), caused by our scalar field is given
by

$\frac{\bm{a}_{5}}{c^{2}}=\frac{1}{2}\bm{\nabla}f_{R}\,,$ (2.15)

where we have assumed $f_{R}\ll 1$, and the background Compton wavelength is
given by

$\lambda_{C}\approx 32\sqrt{\frac{f_{R0}}{10^{-4}}}\,\mathrm{Mpc}\,,$ (2.16)

as in ref. [46], based on ref. [53]. Lastly, using our field redefinition eq.
(2.9) we can define the critical potential eq. (2.6) in terms of our
background scalar field value

$\chi=-\frac{3}{2}f_{R0}\,.$ (2.17)

###  2.4 Symmetron

In contrast with chameleon models, where screening arises due to a variable
scalar mass, symmetron models screen by suppressing the coupling to matter in
dense regions. This is achieved via spontaneous symmetry breaking (SSB), given
the correct form of potential and coupling in the scalar-tensor action, eq.
(2.1). In our case, we choose to work with a quartic potential and quadratic
coupling,

$\displaystyle V(\phi)$
$\displaystyle=\frac{1}{L_{\mathrm{pl}}^{2}}\left[-\frac{1}{2}\left(\frac{\mu}{M_{\mathrm{pl}}}\right)^{2}\phi^{2}+\frac{1}{4}\lambda\phi^{4}\right]\,,$
(2.18) $\displaystyle A(\phi)$
$\displaystyle=1+\frac{1}{2}\left(\frac{M_{\mathrm{pl}}}{M_{\mathrm{sym}}}\right)^{2}\phi^{2}\,,$
(2.19)

where $M_{\mathrm{pl}}^{2}\equiv\hbar c/8\pi G$ is the reduced Planck mass and
$L_{\mathrm{pl}}^{2}\equiv 8\pi G\hbar/c^{3}$ is the reduced Planck length.
This choice of potential and coupling results in an effective potential eq.
(2.5) with the necessary SSB properties. This model has three free parameters:
the dimensionless $\lambda$ governing the scalar self-interactions, the mass
scale $M_{\mathrm{sym}}$ controlling the strength of the coupling to matter,
and $\mu$, representing the bare tachyonic mass of the scalar field.

Extremising the action, eq. (2.1), with respect to the scalar field gives us
the static symmetron EoM, for non-relativistic matter,

$\nabla^{2}\phi=\frac{1}{2L_{c}^{2}}\left[\left(\frac{\rho}{\rho_{\mathrm{SSB}}}-1\right)\phi+\lambda\left(\frac{M_{\mathrm{pl}}}{\mu}\right)^{2}\phi^{3}\right]\,,$
(2.20)

where $L_{c}$ and $\rho_{\mathrm{SSB}}$ are the background Compton wavelength
and SSB density scale, expressed as

$\displaystyle L_{c}^{2}$
$\displaystyle=\frac{1}{2}\left(\frac{M_{\mathrm{pl}}}{\mu}\right)^{2}L_{\mathrm{pl}}^{2}\,,$
(2.21) $\displaystyle\rho_{\mathrm{SSB}}$
$\displaystyle=\left(\frac{M_{\mathrm{sym}}\mu}{M_{\mathrm{pl}}^{2}}\right)^{2}\frac{M_{\mathrm{pl}}}{L_{\mathrm{pl}}^{3}}\,.$
(2.22)

It is possible to reduce this three-parameter EoM to a two-parameter one by
rescaling the field by its vacuum expectation value. We find the vacuum field
value by considering solutions to eq. (2.20) in an infinite box of density
$\rho$. When $\rho>\rho_{\mathrm{SSB}}$ the field is screened and $\phi=0$. If
$\rho<\rho_{\mathrm{SSB}}$, then our scalar field falls into one of two
solutions

$\phi=\pm\phi_{\infty}\equiv\pm\phi_{\mathrm{vev}}\sqrt{1-\frac{\rho}{\rho_{\mathrm{SSB}}}}\,,$
(2.23)

where $\phi_{\mathrm{vev}}$ is the vacuum expectation value given by

$\phi_{\mathrm{vev}}=\frac{1}{\lambda}\frac{\mu}{M_{\mathrm{pl}}}\,.$ (2.24)

Now, rescaling the scalar field, $\varphi=\phi/\phi_{\mathrm{vev}}$, in eq.
(2.20) returns the simpler two-parameter EoM,

$\nabla^{2}\varphi=\frac{1}{2L_{c}^{2}}\left[\left(\frac{\rho}{\rho_{\mathrm{SSB}}}-1\right)+\varphi^{2}\right]\varphi\,.$
(2.25)

The acceleration a test particle experiences due to the scalar fifth force is
given by

$\frac{\bm{a}_{5}}{c^{2}}=-\phi_{\mathrm{vev}}^{2}\left(\frac{M_{\mathrm{pl}}}{M_{\mathrm{sym}}}\right)^{2}\varphi\bm{\nabla}\varphi\,.$
(2.26)

Note that despite reducing the EoM to two parameters, the theory remains a
three-parameter one as the scaling of the fifth force with
$\phi_{\mathrm{vev}}^{2}$ reintroduces a dependence on $\lambda$. The critical
potential eq. (2.7) for this symmetron model is

$\chi=\frac{1}{2}\left(\frac{M_{\mathrm{sym}}}{M_{\mathrm{pl}}}\right)^{2}\,.$
(2.27)

Note that $\chi$ depends only on $M_{\mathrm{sym}}$. However, as we shall see
later, screening is a function of both $M_{\mathrm{sym}}$ and
$L_{\mathrm{c}}$. This discrepancy is a result of the approximations used to
derive the general form of $\chi$, eq. (2.6), primarily that $\beta(\phi)$ is
assumed to be constant.

The symmetry-breaking nature of the symmetron provides a third possible
approximation for the screening surface — in addition to the two described in
section 2.2 for scalar-tensor theories in general. One can assume that the
screening surface universally corresponds to the surface where
$\rho(r=r_{s})=\rho_{\mathrm{SSB}}$, neglecting other relevant factors such as
the coupling strength and Compton wavelength of the theory. We will also
examine the efficacy of this approximation in section 4.

##  3 Numerical Screening Solutions

This section describes the numerical methods we use to obtain isolated galaxy
solutions for the scalar field equations of motion under the chameleon $f(R)$,
eq. (2.13), and symmetron theories, eq. (2.25). Section 3.1 describes our
choice of galaxy model, section 3.2 then describes our numerical solver for
the equations of motion, and section 3.3 outlines the conditions we use to
identify the screening surface within our scalar field solutions.

Our methods are summarised in figure 1, the individual panels of which show a
typical density profile under our galaxy model, the coordinate grid of our
scalar field solver, the scalar field solutions for the same galaxy, and the
screening surfaces identified from the scalar field solutions. These various
panels of the figure will be examined in more detail in the corresponding
parts of this section.

Figure 1:  Various panels illustrating the different stages of our numerical
procedure, as described throughout section 3. All plots are zoomed in to
$5R_{\mathrm{disc}}$, showing only the disc-dominated central region. Top
left: input density profile, for an axisymmetric galaxy with
$M_{\mathrm{vir}}=10^{11.5}\,M_{\odot}$, rescaled by the cosmic mean density.
Top right: coordinate grid, showing every fifth radial and azimuthal angle
gridline. Bottom left, left half: $f_{R}$ field profile rescaled by the
background field value, $|f_{R0}|=10^{-6.4}$. Bottom left, right half: the
curvature-density ratio, eq. (3.10a). Bottom right, left half: symmetron field
profile rescaled by its background field value, for the input parameters
$M_{\mathrm{sym}}=10^{-4.5}\,M_{\mathrm{pl}}$ and
$L_{\mathrm{c}}=10^{-1}\,\mathrm{kpc}$. Bottom right, right half: truncated
scalar field profile to illustrate the symmetron screening condition, eq.
(3.10b).

###  3.1 Galactic Model

Given the $f(R)$ and symmetron EoMs, equations (2.13) and (2.25) respectively,
we can solve the scalar field profile given a suitable galactic density model.
For this, we use a two-component axisymmetric model comprising a dark matter
halo and a galactic disc (note that we are limiting our study to late-type
disc galaxies, not ellipticals). We model the dark matter halo with a Navarro-
Frenk-White (NFW) profile [54]

$\rho(r)=\frac{\rho_{\mathrm{NFW}}}{\dfrac{r}{r_{\mathrm{NFW}}}\left(1+\dfrac{r}{r_{\mathrm{NFW}}}\right)^{2}}\,,$
(3.1)

where $r$ is the spherical radial coordinate, and $\rho_{\mathrm{NFW}}$ and
$r_{\mathrm{NFW}}$ are the characteristic density and length scale. The
galactic disc is represented by a simple double exponential profile222Using a
continuous density profile for a discrete distribution of stars is technically
only correct if the Compton wavelength of our theory is larger than the
typical separation of stellar objects. This is true for all $f(R)$ parameters
we consider, but in the symmetron case our adopted parameter range approaches
the parsec scale at its extreme end. [55]

$\rho(R,z)=\frac{\Sigma_{\mathrm{disc}}}{2z_{\mathrm{disc}}}e^{-R/R_{\mathrm{disc}}}e^{-|z|/z_{\mathrm{disc}}}\,,$
(3.2)

where $R$ represents the cylindrical radial coordinate and $z$ represents the
coordinate perpendicular to the galactic disc, and $\Sigma_{\mathrm{disc}}$,
$R_{\mathrm{disc}}$ and $z_{\mathrm{disc}}$ are the characteristic density and
length scales. In actuality, a more sophisticated galaxy model (e.g. [56])
could also include various other components, such as central bulges, stellar
haloes, and multiple discs; however, we have chosen to exclude these for the
sake of keeping our model as simple as possible while still capturing the key
characteristics.

Including these density components will introduce five extra input parameters
used to scale the density profiles: $\rho_{\mathrm{NFW}}$, $r_{\mathrm{NFW}}$,
$\Sigma_{\mathrm{disc}}$, $R_{\mathrm{disc}}$ and $z_{\mathrm{disc}}$.
Nonetheless, for the remainder of this article, we shall quote our galactic
model in terms of one input parameter: the virial mass, $M_{\mathrm{vir}}$.
This is defined as the mass enclosed within the virial radius,
$R_{\mathrm{vir}}$, the radius within which the mean density inside is 200
times the cosmological critical density. We reduce the five-dimensional
parameter space into a one-dimensional space by using empirical relations
between the various parameters. In doing so, we typically adopt median values
on the relations, ignoring any scatter. Thus, for a given $M_{\mathrm{vir}}$,
the galaxy we model is in some sense a maximally typical one. We describe this
process in further detail in appendix B. The top left panel of figure 1 shows
an example density profile for a galaxy with
$M_{\mathrm{vir}}=10^{11.5}\,M_{\odot}$, corresponding to the following values
for the parameters appearing in our density model:
$\rho_{\mathrm{NFW}}=5.2\times 10^{6}\,M_{\odot}\,\mathrm{kpc}^{-3}$,
$r_{\mathrm{NFW}}=15\,\mathrm{kpc}$, $\Sigma_{\mathrm{disc}}=6.4\times
10^{8}\,M_{\odot}\,\mathrm{kpc}^{-2}$, $R_{\mathrm{disc}}=1.6\,\mathrm{kpc}$
and $z_{\mathrm{disc}}=0.26\,\mathrm{kpc}$.

The adoption of the NFW profile for the dark matter halo presents two
numerical challenges: the central density singularity and the logarithmically
divergent mass. To resolve these issues, we set a constant density within the
central 50 pc, and introduce an outer cutoff at 2.2$R_{\mathrm{vir}}$. A sharp
cutoff at the halo’s edge is well motivated astrophysically, e.g. ref. [57].

In this work, we assume that the galaxy is embedded in a constant-density
background equal to the cosmological critical density. Therefore, we do not
consider environmental effects, such as inhomogeneities in the matter
distribution throughout the Universe, which can manifest as regions of lower
density (voids) and higher density (clusters, walls, and filaments), as well
as neighbouring galaxies and satellites. This simplification enables us to
better investigate other variables that impact screening in galaxies. Overall,
higher-density environments tend to enhance screening effects, resulting in
more fully screened galaxies and larger screened regions in partially screened
galaxies, while the opposite is true for lower-density environments. We
comment on where this will impact our results throughout section 4, and plan
to account for these effects in follow-up work.

###  3.2 Numerical Solver

Using the axisymmetric galaxy density profiles from the previous section, we
solve the $f(R)$ and symmetron equations of motion, given by eqs. (2.13) and
(2.25) respectively. We do this using a numerical solver based on that
employed by ref. [31], updated to incorporate some optimisations suggested by
[40] and to include the symmetron. The solver is ultimately a two-dimensional
version of the $f(R)$ solver used in the cosmological simulation code mg-
gadget [39]: it constructs a coordinate grid, discretises the EoM on this
grid, then solves this discretised EoM using an iterative Gauss-Seidel method.

Although a cylindrical polar coordinate system might seem the natural choice
for an axisymmetric system, we use a spherical polar system (neglecting the
azimuthal coordinate). This is because the radial range of our grid should
stretch from the galactic centre to several Compton wavelengths (eqs. (2.16)
and (2.21)) beyond the galaxy’s edge; the dynamic range is such that a
spherical coordinate system better enables us to use a variable grid
resolution with smaller cells in the central region and larger cells further
out. Specifically, we use a logarithmic radial discretisation, i.e. cells are
equally spaced in the coordinate $u\equiv\ln r$, producing higher resolution
at smaller radii, as required. Similarly, in the polar direction, $\theta$,
because the galaxy has a disc, we need a higher resolution near the disc-plane
$\theta\approx\pi/2$ than at the poles. We can achieve this by discretising in
the coordinate $v\equiv\cos\theta$. In principle, one could assume mirror
symmetry about the disc plane and restrict $v$ to the range $[0,1]$, but we
are often interested in the behaviour of the scalar field in the disc plane,
and so wish to avoid pushing the disc plane to the boundary of the grid. We
thus discretise $v$ in the range $[-1,1]$. A visualisation of the coordinate
grid, out to a maximum radius of $5R_{\mathrm{disc}}$, is given in the top-
right panel of figure 1.

In $u-v$ coordinates, the Laplace operator acting on a function $f$ (dropping
the azimuthal term) is

$D\equiv\nabla^{2}f=\frac{1}{r^{3}}\frac{\partial}{\partial
u}\left(r\frac{\partial f}{\partial
u}\right)+\frac{1}{r^{2}}\frac{\partial}{\partial
v}\left(\sin^{2}\theta\frac{\partial f}{\partial v}\right)\,.$ (3.3)

For later convenience, we have only expressed the derivatives in terms of $u$
and $v$, while keeping other factors that contain $r$ and $\theta$, as their
conversion to $u$ and $v$ is straightforward

We construct a regular grid in $u-v$ space. We use 512 radial ($u$) cells in
the range $[0.05,10^{4}]$ kpc, and 201 angular ($v$) cells. The maximum radial
extent of $10^{4}$ kpc is several times larger than any Compton wavelength we
consider in this work. By using an odd number of angular cells, we ensure that
one row of the grid is centred on the disc plane. Each grid cell is labelled
with indices $(i,j)$, where $i$ represents the $u$-coordinate and $j$
represents the $v$-coordinate, i.e. $u=u_{i}$ and $v=v_{j}$ at cell $(i,j)$.
All other quantities (density, scalar field etc.) can be discretised on this
grid. For example, we can write $\rho_{ij}$ to represent the density at the
centre of cell $(i,j)$.

The Laplace operator acting on a function $f$, eq. (3.3), can then be
discretised as

$\begin{split}D_{ij}&=\frac{1}{h_{u}^{2}r_{i}^{3}}\left[f_{i+1,j}r_{i+1/2}+f_{i-1,j}r_{i-1/2}-f_{ij}\left(r_{i+1/2}+r_{i-1/2}\right)\right]\\\
&+\frac{1}{h_{v}^{2}r_{i}^{2}}\left[f_{i,j+1}s_{j+1/2}^{2}+f_{i,j-1}s_{j-1/2}^{2}-f_{ij}\left(s_{j+1/2}^{2}+s_{j-1/2}^{2}\right)\right]\,,\end{split}$
(3.4)

where $h_{u}$ and $h_{v}$ are the (constant) grid spacings in $u$ and $v$
respectively, and $s\equiv\sin\theta$, so that $s_{j}^{2}=1-v_{j}^{2}$. Note
that half-integer subscripts (e.g. $r_{i+1/2}$) indicate quantities evaluated
at the cell edges rather than the cell centres.

We can now proceed to discretising the EoMs, eqs. (2.13) and (2.25). In the
$f(R)$ case, we follow reference [40] in redefining the scalar field as
$\varphi\equiv(f_{R}/f_{R0})^{1/2}$. Consequently, the discretised forms of
the symmetron and $f(R)$ EoM both turn out to be cubic equations of the form

$\mathcal{L}_{ij}\equiv\varphi_{ij}^{3}+p_{ij}\varphi_{ij}+q_{ij}=0\,,$ (3.5)

where the coefficients $p_{ij}$ and $q_{ij}$ are given in the $f(R)$ case by

$\begin{split}p_{ij}&\equiv-\frac{1}{A_{ij}}\left(B_{ij}(\varphi^{2})+\frac{8\pi
G\delta\rho_{ij}}{3c^{2}f_{R0}}+\frac{R_{0}}{3f_{R0}}\right)\,;\\\
q_{ij}&\equiv\frac{1}{A_{ij}}\frac{R_{0}}{3f_{R0}}\,,\end{split}$ (3.6)

and in the symmetron case by

$\begin{split}p_{ij}&\equiv\frac{\rho_{ij}}{\rho_{\mathrm{SSB}}}+2L_{c}^{2}A_{ij}-1\,;\\\
q_{ij}&\equiv-2L_{c}^{2}B_{ij}(\varphi)\,,\end{split}$ (3.7)

and the objects $A_{ij}$ and $B_{ij}$ are given in turn by

$\begin{split}A_{ij}&=\frac{r_{i+1/2}+r_{i-1/2}}{r_{i}^{3}h_{u}^{2}}+\frac{s_{j+1/2}^{2}+s_{j-1/2}^{2}}{r_{i}^{2}h_{v}^{2}}\,;\\\
B_{ij}(\psi)&=\frac{r_{i+1/2}\psi_{i+1,j}+r_{i-1/2}\psi_{i-1,j}}{r_{i}^{3}h_{u}^{2}}+\frac{s_{j+1/2}^{2}\psi_{i,j+1}+s_{j-1/2}^{2}\psi_{i,j-1}}{r_{i}^{2}h_{v}^{2}}\,,\end{split}$
(3.8)

where $\psi$ represents the argument of $B_{i,j}$, i.e. $\varphi^{2}$ in the
$f(R)$ case and $\varphi$ in the symmetron case.

Before proceeding to solve eq. (3.5) we must specify our boundary conditions
and initial field profile guesses, which we will iteratively converge to the
true solution. Our grid has four boundaries: $v=-1,1$ (i.e. $\theta=\pi,0$)
and $u=\ln r_{\mathrm{min}},\ln r_{\mathrm{max}}$. At the first three, we
impose $\partial\varphi/\partial v,\partial\varphi/\partial u=0$. At the outer
radial boundary, we instead impose $\varphi=\varphi_{\infty}$, i.e. the cosmic
background value of the scalar field. In the $f(R)$ case, this is simply
$f_{R0}$; in the symmetron case, this is the field value that minimises the
effective potential eq. (2.23) evaluated at the cosmic matter density. For
initial guesses, in the $f(R)$ case we simply set $\varphi=1$ everywhere, and
allow the solver to converge on the correct profile. For the symmetron case,
we speed up the convergence by setting the field to 1 everywhere, except in
the region where $\rho>\rho_{\mathrm{SSB}}$ which is assumed to be screened
and thus set to 0.

We solve eq. (3.5) using a Gauss-Seidel iteration method. Here, our 2D grid is
decomposed into a chequerboard pattern, i.e. we alternately assign cells with
‘red’ and ‘black’ labels. Each iteration $n$ is then comprised of two halves:
a red sweep and a black sweep. During the ‘red sweep’, the scalar field value
is changed in all of the red cells using a Newton-Raphson update,

$\varphi_{i,j}^{n}=\varphi_{i,j}^{n-1}-\frac{\mathcal{L}_{i,j}}{\mathcal{L}^{\prime}_{i,j}}\,,$
(3.9)

where
$\mathcal{L}^{\prime}\equiv\partial\mathcal{L}/\partial\varphi=3\varphi^{2}+p$.
Note that $\mathcal{L}_{i,j}$ here is evaluated using $\varphi^{n-1}$, i.e.
the values of $\varphi$ calculated at the previous Gauss-Seidel iteration (or
the initial guess in the case $n=1$).

Next, we perform a ‘black sweep’ in which the scalar field is updated in all
of the black cells using the same expression. There is however one key
difference in the black sweep: rather than using the leftover values from the
previous Gauss-Seidel iteration ($n-1$) as was done during the red sweep, one
uses the neighbouring red cell values that have just been updated during the
red sweep, i.e. in the first half of iteration $n$.

This iterative process repeats, updating the field values at each grid point
until a converged solution is obtained. The iteration stops when the
difference between the field in the previous step and the current one is less
than a given threshold. In the $f(R)$ case we use a threshold of $10^{-7}$,
and in the symmetron case we use $10^{-14}$. We find these values give robust,
converged results.

The left semicircles in the bottom row of figure 1 show outputs of the solver
for the same example galaxy shown in the upper left panel of the figure, with
a mass $M_{\mathrm{vir}}=10^{11.5}\,M_{\odot}$, theory parameters of
$f_{R0}=-10^{-6.4}$ in the $f(R)$ case,
$M_{\mathrm{sym}}=10^{-4.5}\,M_{\mathrm{pl}}$ and
$L_{\mathrm{c}}=10^{-1}\,\mathrm{kpc}$ in the symmetron case. In both cases,
the field solutions are plotted to a maximum radius of $5R_{\mathrm{disc}}$.

###  3.3 Screening Conditions

This section aims to establish a universal and reliable screening condition to
locate the (generally non-spherical) screening surface, $r_{s}(\theta)$, in a
galaxy’s scalar field profile. We will find that our conditions to calculate
$r_{s}(\theta)$ are given by

$\displaystyle\left.\frac{c^{2}\delta R(r,\theta)}{8\pi
G\delta\rho(r,\theta)}\right|_{r=r_{s}(\theta)}=0.9\,,$ (3.10a)
$\displaystyle\left.\frac{\varphi(r,\theta)}{\varphi_{\infty}}\right|_{r=r_{s}(\theta)}=0.1\,,$
(3.10b)

for the $f(R)$ and symmetron models respectively. In other words, for a given
angle, the screening surface is the surface at which the $f(R)$ curvature-
density ratio or the symmetron field value matches a given threshold. The
right sides of the bottom row panels in figure 1 display the $f(R)$ curvature-
density ratio and the symmetron field with a colour map truncated at the
threshold. The sharp transition in both of these cases vindicates our chosen
threshold.

In the $f(R)$ case, the screening condition, eq. (3.10a), is a threshold on
the curvature-density ratio, motivated by the two terms on the right-hand side
of the $f(R)$ EoM eq. (2.13). Specifically, inside the screened region the
field is flat and the Laplacian is zero, thus the two terms on the right-hand
side of the EoM are approximately equal in magnitude, and the curvature-
density ratio is unity. Conversely, in the unscreened regime, the field adopts
its background value, causing $\delta R$, eq. (2.14), and thus the curvature-
density ratio to become zero. We find that the curvature-density ratio
provides a clearer demarcation between the screened and unscreened region than
the scalar field, as can be seen in the bottom left plot of figure 1.
Nevertheless, there are still instances where the transition is more gradual,
particularly for smaller values of $|f_{R0}|$. We find that a threshold on the
curvature-density of 0.9 correctly separates the screened/unscreened regions
and produces universal behaviour across the whole parameter space.

For the symmetron condition, we use a threshold value on the field itself.
This is motivated by the coupling function, eq. (2.19), having a first
derivative proportional to the scalar field; therefore the density-dependent
term in the EoM is then also proportional to the scalar field. Setting this
threshold as a ratio of the background scalar field allows for a screening
condition that can be universally applied over the whole parameter space. We
determine that a threshold value of 10% of the background scalar field value
provides consistent results, with the scalar field profile being visually
flat, and thus suppressing fifth forces, within a substantial region of
$r_{s}$.

Figure 2 shows the value of the $f(R)$ curvature-density parameter and the
symmetron field along the plane of the galactic disc, as well as the magnitude
of the fifth forces of both theories. We rescale the radial coordinates by the
derived $r_{s}$ along the disc plane to better visualise the behaviour over
the whole parameter space. Taking figures 1 and 2 together, we are able to
justify our choice of thresholds, namely 0.9 for the curvature-density
parameter and 0.1 for the symmetron field, finding that this choice allows for
consistent results, that suppress the fifth forces inside of the screened
region, across our full range of parameters.

Figure 2:  Values of the $f(R)$ curvature-density ratio (first row), the
$f(R)$ fifth force (second row, the symmetron field (third row) and the
symmetron fifth force (fourth row), in the plane of the galactic disc, plotted
against radial coordinates rescaled by $r_{s}$ (the radius at which the
screening surface meets the disc plane), across a range of theory parameters
and galaxy masses. In the $f(R)$ case, $\log_{10}(|f_{R0}|)\in[-8,-5]$ and
$\log_{10}(M_{\mathrm{vir}}/M_{\odot})\in[10,13.5]$, both in increments of
0.2. In the symmetron case,
$\log_{10}(M_{\mathrm{sym}}/M_{\mathrm{pl}})\in[-6.5,-3]$,
$\log_{10}(L_{\mathrm{c}}/\mathrm{kpc})\in[-3,3]$ and
$\log_{10}(M_{\mathrm{vir}}/M_{\odot})\in[10,13.5]$, all in increments of 0.5.
Profiles are cut off at the extent of the galaxy, $2.2R_{\mathrm{vir}}$.
Dashed vertical lines represent the location of the screening radius. Dashed
horizontal lines correspond to the screening condition threshold. Purple lines
on the symmetron plot indicate solutions that evade the fully unscreened
central scalar field definition but are properly considered fully unscreened
by the central Laplacian definition.

We classify a galaxy as fully screened if the screening surface determined
according to eqs. (3.10) is greater than the extent of the galaxy, or if SSB
has not yet occurred in the symmetron case, i.e.
$\bar{\rho}>\rho_{\mathrm{SSB}}$. Conversely, a galaxy is considered fully
unscreened if there is no screened region within the galaxy, i.e., if the
conditions given by eq. (3.10) are not satisfied anywhere. However, we also
find a small subset of galaxies for which the screening conditions are
satisfied somewhere but nonetheless have non-negligible fifth forces in their
central regions. To catch these cases, we consider the innermost field value.
In the $f(R)$ case, solutions satisfying $f_{R}(r\approx 0)>10^{-3}f_{R0}$ are
fully unscreened. We find that an analogous threshold on the field central
field value can not be used in the symmetron case, as some field solutions
yield small central field values but quickly grow away from the centre, giving
an unscreened central region. So, for the symmetron, we instead use the
Laplacian of the field: solutions satisfying the condition
$\nabla^{2}\varphi(r\approx
0)>10^{-1}\mathrm{max}\left(\nabla^{2}\varphi\right)$ are fully unscreened.

The symmetron solutions that avoid the central field threshold condition (i.e.
satisfying $\phi(r\approx 0)<10^{-3}\phi_{\infty}$) but are correctly
classified as fully unscreened by the central Laplacian threshold condition
are highlighted in figure 2. Here we see that these edge case field profiles
grow to be large well within the screening radius, causing a non-zero fifth
force in this region, as shown in the bottom panel.

##  4 The Screening of Galaxies

###  4.1 Position and Shape of Screening Surfaces

Figure 3:  Screening surfaces across a range of theory parameters, in a galaxy
with mass $M_{\mathrm{vir}}=10^{11.5}\,M_{\odot}$, for $f(R)$ (top) and
symmetron (bottom) models. In the symmetron case, we vary $M_{\mathrm{sym}}$
holding $L_{\mathrm{c}}$ fixed (bottom left) and vary $L_{\mathrm{c}}$ holding
$M_{\mathrm{sym}}$ fixed (bottom right). Sky blue dashed lines are the
calculated screening surfaces with theory parameters labelled. Grey dashed-
dotted lines are approximate solutions: numerically integrating the critical
potential (2.6) in the $f(R)$ case; isodensity contours equal to
$\rho_{\mathrm{SSB}}$ in the symmetron case. In each case, the right-hand
semi-circle plots the radial range of the galaxy out to $5r_{\mathrm{NFW}}$
(approximately 75 kpc), while the left-hand semi-circle zooms into the central
region, out to $5R_{\mathrm{disc}}$ (approximately 8 kpc). The coloured
background regions represent the adopted matter density profile, plotted as a
fraction of the cosmic mean density.

Figure 3 shows an array of solutions for the screening surface (dashed lines),
over a range of $f(R)$ and symmetron parameters, for a galaxy with
$M_{\mathrm{vir}}=10^{11.5}\,M_{\odot}$. Generally, decreasing the adopted
value of $|f_{R0}|$ or $M_{\mathrm{sym}}$, or increasing $L_{\mathrm{c}}$,
pushes the screening surface further out. We separate these solutions into
those for which the screening surfaces lie within the disc-dominated region
(left-hand semi-circles) and the halo-dominated region (right-hand semi-
circles). The figure shows that screening surfaces that occur in the disc-
dominated region are elongated along the galactic disc, whereas larger
screening surfaces are spherical and trace the dark matter component.

We also include approximate screening surfaces (dashed-dotted lines) for
comparison. In the case of $f(R)$, the approximate screening surface is
computed by numerically solving eq. (2.6), where $\chi$ is given by the $f(R)$
critical potential, eq. (2.17), and $\delta\rho(r)$ is the galactic density
averaged over spherical shells. The approximate symmetron screening surfaces
are isodensity contours with a value equal to the symmetry breaking scale,
i.e. $\rho(r=r_{s}(\theta))=\rho_{\mathrm{SSB}}$.

In the $f(R)$ case (top left panel), the true screening surfaces (dashed
lines) in the disc-dominated region follow the isodensity contours away from
the disc plane, but deviate nearer the disc, where the isodensity contours are
flattened to a larger extent. We find that the screening surfaces lie
somewhere between the highly flattened isodensity contours of our galactic
model and the isopotential contours (not shown), which are significantly more
spherical than the screening surfaces.

The non-spherical nature of the screening surface in the disc-dominated region
of the galaxy (left-hand side of the $f(R)$ plot) indicates that using
approximations (dashed-dotted lines) that result in spherical screening radii,
such as eq. (2.6), can fail to capture the true nature of the screened region.
This is especially relevant when searching for observables along the plane of
the disc, where $r_{s}$ will be underestimated as the spherical screening
surface meets the disc plane at a smaller radius than the true, flattened
screening surface. As a result, for a given value of the theory parameter, one
may expect the presence of fifth force observables at a certain location in
the disc plane because it lies outside the assumed spherical screening
surface, whereas, in reality, the location lies within the true, flattened,
screening surface, and so the fifth force is suppressed. This in turn could
lead to exaggerated constraints on the theory parameters.

In the halo-dominated region of the galaxy (right-hand side of the $f(R)$
plot) we see the approximation overestimates the screening radius for smaller
values of $|f_{R0}|$ — corresponding to larger values of $r_{s}$. This leads
to a significant overestimation of the amount of dark matter that is screened
by this approximation. This behaviour is likely due to a breakdown of the
approximations utilised to derive the critical potential at large $r_{s}$
values, namely that the field can be linearised around the screening surface,
and that the scalar mass is negligible in this region. The effect that this
overestimation of $r_{s}$ has on observables is discussed further in section
4.2.

By contrast with the $f(R)$ screening surfaces, the screening surfaces under
the symmetron theory (bottom plots) follow the isodensity contours very
closely. However, this behaviour breaks down in the large $L_{\mathrm{c}}$,
large $M_{\mathrm{sym}}$ regime, as seen with the larger screening radii in
the left-sided symmetron plots, where the screening surface fails to trace the
finer details of the isodensity contour, especially near the disc plane. The
reasoning for this breakdown is due to the physics associated with each
parameter. The Compton wavelength, $L_{\mathrm{c}}$, is the characteristic
scale on which the field can vary, i.e. the field will not respond to changes
on scales much smaller than the Compton wavelength, but in the large
$L_{\mathrm{c}}$ regime the scale height of the disc profile is much smaller
than $L_{\mathrm{c}}$. In other words, a larger Compton wavelength responds to
mass averaged over a larger distance, causing the field profile to smooth out
the finer details of the disc density profile. Meanwhile, for the coupling
constant, $M_{\mathrm{sym}}$, the coupling function, and thus the density-
dependent term in the EoM, eq. (2.25), have an inverse proportionality to
$M_{\mathrm{sym}}^{2}$, meaning a large change in density has less impact on
the field profile.

The approximation of symmetron screening surfaces as isodensity contours
matching $\rho_{\mathrm{SSB}}$ is a natural consequence of the highly
sensitive relationship between the screened region and the symmetry-breaking
scale. These approximate solutions are shown as grey dashed-dotted lines on
the symmetron plots of figure 3. For the most part, this approximation does a
good job of calculating the true value of $r_{s}$, especially in the small
$M_{\mathrm{sym}}$, small $L_{\mathrm{c}}$ regime. Outside of this region of
parameter space, the approximation always overestimates the value of $r_{s}$
when compared to the true value. This overestimation is especially prevalent
in the large $L_{\mathrm{c}}$, large $M_{\mathrm{sym}}$ regime.

###  4.2 Amounts of Galactic Mass Screened

One of the main results of our analysis is contained within figure 4, where we
show the percentage of mass contained within the screening surface over a
range of galactic and theory parameters. The left column represents the disc
mass and the right column is the halo mass. In the symmetron plots, we vary
one parameter with the other parameter fixed. This figure allows us to study
how our system transitions from screened to unscreened across the parameter
space, which we can compare with the transition from screened to unscreened
given by the binary screening condition, eq. (2.17) in the $f(R)$ case and eq.
(2.27) in the symmetron case. As a proxy for the ‘true’ answer, we include a
dashed contour line to show the parameter boundary where 50% of the mass is
screened in each case333There is ultimately a degree of arbitrariness in our
choice of 50% as a proxy. It intuitively feels the fairest way to divide the
spectrum of solutions in two, but one could instead argue for a lower (higher)
threshold such that the partially screened solutions are preferentially
grouped with the fully screened (unscreened) solutions.. This line can then be
directly compared with the binary screening condition boundary, which we
represent with a dotted line in the $f(R)$ and symmetron $M_{\mathrm{sym}}$
plots respectively (not included in the $L_{\mathrm{c}}$ plot since the
symmetron critical potential is a function of $M_{\mathrm{sym}}$ only).

Figure 4:  Fraction of disc mass (left) and dark matter (right) screened, for
a range of $f(R)$ (top) and symmetron parameters, with fixed
$L_{\mathrm{c}}=10^{-1}\,\mathrm{kpc}$ (middle) and with fixed
$M_{\mathrm{sym}}=10^{-4.5}\,M_{\mathrm{pl}}$ (bottom). Sky blue dashed
contours show where 50% of mass is screened. Grey dashed-dotted contours
indicate where 50% of mass is screened when $r_{s}$ is approximated: by
numerically integrating eq. (2.6) for the $f(R)$ critical potential, or by
finding isodensity contours equal to $\rho_{\mathrm{SSB}}$ for the symmetron
model. Grey dotted lines correspond to the binary screening condition;
solutions below this line would be considered ‘fully screened’. Solid white
lines separate the fully unscreened and partially screened solutions.

All the sub-figures contained in figure 4 exhibit some similar behaviours. The
galactic disc component transitions sharply from being fully screened to fully
unscreened, with the two regimes separated in the parameter space of galaxy
and model parameters by only a narrow band in which the disc is partially
screened. In contrast, the fraction of screened dark matter varies gradually,
due to the significant dynamical range of the halo, which spans several orders
of magnitude in density. Another notable feature is the discrepancy between
the regions considered fully screened and unscreened by the binary condition
and the region where half of the mass is screened. This is explored further
throughout this section. In all cases, the $r_{s}$ approximations (dashed-
dotted lines) are closer to the true 50% screened solution than the binary
screening condition.

Focusing first on the $f(R)$ case, we see that the $f(R)$ binary screening
condition, eq. (2.17), does not correctly characterise the screening of the
galactic disc. Specifically, the region of parameter space considered ‘fully
unscreened’ by the binary condition (above the dotted line in the top left
panel) includes roughly an order of magnitude in which the galactic disc is
actually entirely screened. This ultimately means that searching for
signatures of screening theories in ‘fully unscreened’ galaxies may lead to
false negative results since the visible stellar and gas components are fully
screened. For example, when looking at stellar-gas offsets, ref. [24, 26], in
this region of parameter space the disc is fully screened and there is
therefore no displacement between the stellar and gas components. However,
using eq. (2.17) as a binary screening threshold would lead one to assume that
the disc is unscreened and thus incorrectly predict the presence of an offset.
In practice, this could entail relaxing current constraints by roughly an
order of magnitude.

By contrast with the disc, for the halo the binary screening condition kicks
in at a point where the halo is only partially screened, so there will be
haloes (below the dotted line in the top right panel) which will be classed as
screened but actually carry a significant unscreened component in their
outskirts, as much as 60-70% of the total halo mass. Similarly, there will be
haloes (above the dotted line) which will be classed as unscreened but
actually have screened central regions, comprising as much as 30-40% of the
total halo mass. This is the fundamental flaw in the binary screening
condition: it assumes a step change in a quantity which actually changes
gradually over the course of several orders of magnitude in parameter space.
In summary, the $f(R)$ binary screening condition underestimates unscreened
mass when it classifies an object as screened, and it underestimates screened
mass when it classifies an object as unscreened.

The various panels in figure 4 also include dashed-dotted lines to indicate
when 50% of the mass is screened when $r_{s}$ is derived by approximate means.
In the $f(R)$ case (top row), this is when $r_{s}$ is derived by numerically
integrating eq. (2.6) for a spherically averaged version of our galaxy density
profile. This approximation solution yields results significantly closer to
the exact solution than the binary screening condition. For the galactic disc
(top left panel), this approximation contour is almost identical to the
equivalent numerical solution contour, and the difference in the proportion of
mass screened between the two calculations is no more than 5% across the range
of parameters we consider. However, as we saw in figure 3, the true and
approximate solutions for $r_{s}$ can have vastly different shapes and the
spherical nature of the approximate screening surface results in an
underestimation of $r_{s}$ along the galactic disc and an overestimation along
the perpendicular axis.

For the dark matter component (top right panel), the approximate value of
$r_{s}$, and therefore the percentage of mass screened, is very close to the
actual value for low values of $r_{s}$. However, as $r_{s}$ becomes larger,
the approximate solution begins to overestimate $r_{s}$ compared to the actual
solution, as previously seen in figure 3. In the extreme case (i.e. the lower
right of the $f(R)$ parameter space in figure 4) the estimated proportion of
screened dark matter is as much as 20% larger than the actual solution. In
other words, if the 50% threshold were increased to, e.g., 80%, then the two
lines would be further separated. Ultimately this will result in the predicted
fifth forces sourced by the galaxy being larger than in the true case, but not
to the same extent as the binary screening condition would predict.

For the symmetron model (middle and bottom rows), the binary screening
condition is much less accurate. As discussed in section 2.4, the condition
depends only on $M_{\mathrm{sym}}$, but even in the middle row, with varying
$M_{\mathrm{sym}}$ and fixed $L_{\mathrm{c}}$, this estimate of screening
disagrees with our solution by several orders of magnitude. Looking at the row
with varying $L_{\mathrm{c}}$ and fixed $M_{\mathrm{sym}}$, it is clear that
the location of the screening surface is a strong function of
$L_{\mathrm{c}}$. Therefore, this simplified version of the symmetron binary
screening condition should be avoided.

Instead, a more useful approximation for screening is obtained by assuming the
screening surface is coincident with the surface where the total density
equals $\rho_{\mathrm{SSB}}$. The parameter space that results in 50% of each
mass component being enclosed by this approximate screening surface is shown
in the symmetron panels of figure 4 as a dashed-dotted contour. As discussed
in section 4.1, this approximation will never underestimate the value of
$r_{s}$, and will overestimate $r_{s}$ for large values of $M_{\mathrm{sym}}$
or $L_{\mathrm{c}}$. Looking first at the disc mass, in both cases (varying
$M_{\mathrm{sym}}$ and varying $L_{\mathrm{c}}$) the approximation contour
follows a similar behaviour to the true contour, but somewhat overestimates
the amount of disc mass screened. Notably, the distance in parameter space
between the true and approximate 50% screened contours is reduced in the case
where $L_{\mathrm{c}}$ is smaller (towards the right of the bottom-left panel
with varying $L_{\mathrm{c}}$ and fixed $M_{\mathrm{sym}}$). For the dark
matter component, we see that the approximate and true 50% screening lines are
indistinguishable in the middle right panel with varying $M_{\mathrm{sym}}$
and fixed $L_{\mathrm{c}}$. This is due to the fact that $L_{\mathrm{c}}$ is
fixed at a low value of $10^{-1}\,\mathrm{kpc}$, and so to get a large
screening radius that screens 50% of the halo we need a small value of
$M_{\mathrm{sym}}$, thus falling into the regime where this approximation
works best (small $M_{\mathrm{sym}}$, small $L_{\mathrm{c}}$). In the bottom
right panel with varying $L_{\mathrm{c}}$ and fixed $M_{\mathrm{sym}}$, the
50% line is a good match for larger galaxies, but for smaller galaxies the
approximation overestimates the amount of mass screened. As length scales get
smaller they become comparable with $L_{\mathrm{c}}$ and the field can not
respond to the density changes quickly enough to screen exactly along the
isodensity contour. In summary, the approximation does not account for the
full effects of varying $M_{\mathrm{sym}}$ and $L_{\mathrm{c}}$, so it falters
in the regimes in which these effects are strongest.

One significant caveat to this discussion is that we do not consider any
environmental effects, e.g., the presence of over-/underdensities within a few
Compton wavelengths of our target galaxies. In general, these would serve to
vary the location of the screening surface of a given galaxy from that
calculated in isolation. However, because our exact solutions and the
approximate solutions we compare them against should all shift under such
environmental effects, the conclusions we draw from figure 4 remain robust.

###  4.3 Milky Way Observables

Figure 5:  Positions where non-spherical screening surface intersects the disc
plane of a Milky Way-like galaxy over $f(R)$ (top) and symmetron (bottom)
parameter space. In the $f(R)$ plot, the dashed line shows the true $r_{s}$
value, whereas the dashed-dotted line shows the approximate $r_{s}$ values,
derived by numerically integrating eq. (2.6) for the $f(R)$ critical
potential. In the symmetron plots, the colourmap shows the true $r_{s}$ value,
the dashed line indicates where $r_{s}$ is equal to the Milky Way centre-Solar
System separation of 8 kpc, and the dashed-dotted line shows where the
approximate $r_{s}$, derived by finding isodensity contours equal to
$\rho_{\mathrm{SSB}}$, is equal to the Milky Way centre-Solar System
separation. Dotted lines indicate the binary screening condition. In the
symmetron plot, a solid white line separates the fully unscreened and
partially screened solutions, and the white region shows where SSB does not
occur, i.e. the cosmic background density, $\bar{\rho}$, is greater than
$\rho_{\mathrm{SSB}}$.

Figure 5 shows the derived screening radius along the plane of the galactic
disc, over the entire model parameter spaces, for a Milky Way-like galaxy with
$M_{\mathrm{vir}}=1.5\times 10^{12}M_{\odot}$ [56]. The approximate distance
from the Milky Way centre to the Solar System (8 kpc) is marked by the dashed
lines. The dashed-dotted lines show the screening radius derived via
approximate methods, namely, numerically integrating eq. (2.6) given the
$f(R)$ critical potential, eq. (2.17), or defining the screening surface as
the isodensity contour equal to $\rho_{\mathrm{SSB}}$ in the symmetron case.
The dotted lines represent the binary screening boundary, separating the
parameters that cause solutions to be considered ‘fully screened’ and ‘fully
unscreened’. Specifically, all $|f_{R0}|$ values less than $10^{-6.7}$ and all
$M_{\mathrm{sym}}$ values less than $10^{-3.1}\,M_{\mathrm{pl}}$ are
considered ‘fully screened’ by the binary screening conditions.

Stringent constraints based on tests of celestial gravitational dynamics
require that the Milky Way screening radius extends to at least 8 kpc,
encompassing the entire Solar System [58]. Figure 5 allows us to determine the
parameter space permitted by such constraints, namely $|f_{R0}|\lesssim
10^{-6.0}$ and the portion of the symmetron parameter space satisfying
$M_{\mathrm{sym}}\lesssim 10^{-3.6}M_{\mathrm{pl}}$ and to the right of the
dashed line on the lower panel. There are a few caveats to these constraints.
Firstly, they only apply to this particular galaxy model and choice of mass.
We do not expect the Milky Way to be a maximally-typical galaxy but for
consistency with our previous analysis, we have modelled it as such. In
addition, there remains significant empirical uncertainty regarding the total
mass of the Milky Way halo [59], which would have a non-negligible effect on
the bounds we obtain. We have also neglected any environmental effects, such
as the influence of the Local Group and the galaxy’s wider cosmographic
context.

Given the strict constraint on the minimum value of the Milky Way’s screening
radius, the question arises as to whether there can still be observable
effects in the periphery of the galactic disc. For this galaxy, the disc
component becomes subdominant along the disc plane after roughly 14 kpc.
Looking at either panel of figure 5, between the $r_{s}=8$ kpc contour and the
parameter space where $r_{s}=14$ kpc, there is only a small slice of parameter
space in which the Solar System is screened but the bulk of the galactic disc
is not. This means it would require a great deal of fortune to be in a
universe adopting these parameter values and allowing us to explain observable
effects in the disc via fifth forces.

Nevertheless, there remains several orders of magnitude of unconstrained
parameter space in which the screening radius lies between the edge of the
galactic disc and the virial radius of $\sim$ 300 kpc. As such, observable
consequences of scalar fields with screening mechanisms located in the Milky
Way’s dark matter halo, are entirely plausible. For example, the kinematics of
halo tracers or the detection of equivalence-principle-violating asymmetries
in the leading and trailing tails of stellar streams, from tidally disrupted
dwarf galaxies in the dark matter halo, as in refs. [28, 60, 61, 62].

An additional feature in figure 5 worth mentioning is the turn-off in the top
right of the symmetron plot. Away from this region, lines of equal screening
radius are diagonals across parameter space. These diagonals correspond with
lines of equal $\rho_{\mathrm{SSB}}$ ($\rho_{\mathrm{SSB}}\propto
M_{\mathrm{sym}}/L_{\mathrm{c}}^{2}$; cf. eq. (2.22)), and as we saw in
figure. 3, screening surfaces closely match the $\rho=\rho_{\mathrm{SSB}}$
surfaces in the low $M_{\mathrm{sym}}$, low $L_{\mathrm{c}}$ regime. The turn-
off is another symptom of the breakdown of this behaviour in the large
$M_{\mathrm{sym}}$, large $L_{\mathrm{c}}$ limit.

##  5 Conclusions

In this work, we have thoroughly investigated the effects of scalar
modifications of gravity on galactic scales, for both chameleon $f(R)$ and
symmetron screening theories, with a particular focus on how screening differs
between full numerical solutions and common approximations used to derive
constraints. To determine the screened region, we solved the quasi-static
scalar field EoM numerically, given an input density described by a maximally-
typical galaxy consisting of a double-exponential disc profile and a dark
matter NFW halo. We developed new and robust screening conditions to compute
the position of the screening surface based on the scalar field values, eqs.
(3.10).

We examined the shape of our calculated screening surfaces. We find that
screening surfaces in the halo-dominated region follow the spherical
isodensity contours. In the disc-dominated region, $f(R)$ screening surfaces
almost follow the isodensity contours, but are overall more spherical than the
disc profile. Symmetron screening surfaces closely follow isodensity contours
in the disc-dominated region, except in the large $M_{\mathrm{sym}}$, large
$L_{\mathrm{c}}$ limit where the individual microphysics of each parameter
becomes important.

The simplest approximation we explore is the binary screening condition, ref.
[46], where if a galaxy has a virial potential,
$\Phi_{\mathrm{vir}}=GM_{\mathrm{vir}}/R_{\mathrm{vir}}$, greater/smaller than
the critical potential then it is classified as fully screened/unscreened,
ignoring the effects of partial screening. We note that the symmetron binary
screening condition, eq. (2.27), is not suitable as a proxy for screening.
This is because the approximations required to derive the critical potential,
eq. (2.6), are incompatible with the symmetron theory.

We explored the effectiveness of the binary screening condition as a proxy to
$f(R)$ screening, for observables in the visible disc and dark matter halo. We
found that the binary screening condition drastically underestimates the
amount of parameter space in which the galactic disc is fully screened. In
this fully screened region of parameter space, using the binary screening
condition would erroneously define the galaxy as ‘fully unscreened’, i.e. it
would lead one to predict an observable signature that would not occur. It
follows that constraints set using the binary screening condition in this
contested region of parameter space would be tighter than if one numerically
computed the true screening surface.

We also compare our true solutions against approximate screening surfaces. In
the $f(R)$ case, $r_{s}$ can be approximated by numerically solving the
critical potential, eq. (2.6), for a spherically averaged galactic density. In
this case, the resulting approximate screening surfaces are spherical. For the
symmetron model, it is possible to approximate $r_{s}$ as the location where
the density is equal to the symmetron breaking scale, i.e.
$\rho(r=r_{s})=\rho_{\mathrm{SSB}}$.

For the halo, we found the binary screening condition sits in the region of
parameter space where the haloes are partially screened. This means objects
classified as “screened” actually have 60-70% unscreened mass, while objects
classified as “unscreened” actually have 30-40% screened mass. This is the
fundamental flaw of the binary screening condition, the fact that the binary
screening condition assumes a discrete step change between fully screened and
unscreened, but actually the degree of screening varies gradually over many
orders of magnitude in parameter space.

Direct approximation of $r_{s}$ allows for partial screening solutions and is
overall a closer approximation to the true screening solutions than the binary
screening condition. In regions where the disc component is not negligible,
the spherical approximate screening surface will underestimate $r_{s}$ along
the disc and overestimate $r_{s}$ perpendicular to the disc. As a result of
this, we might expect to see an observable feature at a location along the
plane of the disc, but in fact, the location is within the true screening
surface and therefore there is no observable feature. In addition, this
approximation overestimates $r_{s}$ for small values of $|f_{R0}|$,
corresponding to large screening surfaces, and ultimately overestimates the
amount of screened dark matter, however not to the same extent as the binary
screening condition.

We have shown that defining the symmetron $r_{s}$ as the radius at which the
total density equals $\rho_{\mathrm{SSB}}$ is a useful approximation for a
large portion of the parameter space. However, in the large
$M_{\mathrm{sym}}$, large $L_{\mathrm{c}}$ limit, varying the two parameters
separately has consequences for the scalar field profile that are not captured
purely through a dependence on $\rho_{\mathrm{SSB}}$ and the approximation
breaks down. In general, the approximate screening surface will never
underestimate $r_{s}$ and the amount of overestimation depends on how large
$M_{\mathrm{sym}}$ and $L_{\mathrm{c}}$ are. In smaller galaxies, the
approximation may break down if the galactic scale length becomes comparable
to $L_{\mathrm{c}}$.

Finally, we explored the possibility of finding observables in a Milky Way-
like galaxy. Assuming that the Solar System is screened by the Milky Way, we
found that (for our isolated test Milky-Way galaxy) $|f_{R0}|$ must be less
than $10^{-6.0}$ and $M_{\mathrm{sym}}$ must be less than
$10^{-3.6}\,M_{\mathrm{pl}}$ with $L_{\mathrm{c}}$ being to the right of the
dashed line in the bottom panel of figure 5. We saw that there is only a small
slice of model parameter space in which there could be observable effects in
the galactic disc. This suggests we would have to be rather fortunate to exist
in a universe where we could detect the modifications to GR considered in this
work within our own galactic disc — however, searches for observables in the
halo remain promising, thanks in part to the large dynamical range of the halo
density.

As cosmological surveys continue to improve the precision of their
measurements it is important that theoretical predictions continue to improve
in accuracy alongside them. It is in this spirit that we have presented, in
this paper, new definitions of screening for both the $f(R)$ chameleon model
and the symmetron model of screening. We have shown that it is important to
differentiate between situations where the galactic disc is screened, but much
of the dark matter halo is unscreened, and situations where the majority of
the total mass of the galaxy is screened. We have shown that approximations
currently used in the literature could lead to order-of-magnitude errors in
the constraints placed on model parameters from galactic observations. We
intend to present a reanalysis of current constraints in the near future.

##  Appendix A $\chi$ Parameter Derivation

This appendix provides a complete derivation of the critical potential, eq.
(2.6), introduced in section 2.2. This prevalent approximation to screening
comes from solving the general EoM eq. (2.4) for a static spherically-
symmetric source, $\delta\rho(r)$, of radius $R$, in a constant density
background, $\bar{\rho}$. This method can be applied to scenarios such as a
star embedded in the galactic medium or a dark matter halo embedded in the
cosmological background.

Far from the object, the field settles into the minimum of the effective
potential, i.e.
$\phi\rightarrow\bar{\phi}\equiv\phi_{\mathrm{min}}(\bar{\rho})$. Near the
object, if the field is able to reach the minimum of the potential inside the
object then $\nabla^{2}\phi=0$, meaning $\phi=\phi_{\mathrm{min}}(\rho)$, and
the fifth force is screened due to the vanishing field gradients. Otherwise,
the EoM needs to be solved, which can be done if we assume the object is
static and spherically symmetric, and therefore
$\Box\rightarrow\nabla^{2}(r)$. If the field remains close to the background
value, $\bar{\phi}$, we linearise as $\phi=\bar{\phi}+\delta\phi$, resulting
in the EoM

$\nabla^{2}\phi=m_{0}^{2}\delta\phi+\frac{8\pi
G}{c^{2}}\beta(\bar{\phi})\delta\rho(r)\,,$ (A.1)

where $m_{0}^{2}=V^{\prime\prime}(\bar{\phi})$. If we assume that $\phi$ is
sufficient light in the background environment, such that $m_{0}^{2}R\ll 1$,
we find that we can consistently neglect the first term on the RHS of our
linearised EoM. In summary, well inside the object, the field is at the
minimum of the effective potential, then, near $r_{s}$, it begins to roll
towards the background value.

Now we can solve the EoM outside the screened region by integrating eq. (A.1)
(with the assumptions laid out in the previous paragraph) to get

$\frac{d\phi}{dr}=\frac{2G\beta(\bar{\phi})}{c^{2}}\frac{\mathcal{M}(r)-\mathcal{M}(r_{s})}{r^{2}}\,,$
(A.2)

where $\mathcal{M}(r)=\int^{r}_{0}4\pi r^{\prime
2}\delta\rho(r^{\prime})dr^{\prime}$. Then we determine the screening radius
by integrating again, this time from $r_{s}$, where
$\phi=\phi_{s}\approx\phi_{\mathrm{min}}(\rho)$, to infinity, where the field
assumes the background value $\phi=\bar{\phi}$,

$\bar{\phi}-\phi_{s}=\frac{2G\beta(\bar{\phi})}{c^{2}}\left(\frac{\mathcal{M}(r_{s})}{r_{s}}+\int^{\infty}_{r_{s}}\frac{\mathcal{M}(r)}{r^{2}}dr\right)\,.$
(A.3)

Finally, integrating by parts and using our expression for $\mathcal{M}(r)$,
we find the expression

$\chi\equiv\frac{\bar{\phi}}{2\beta(\bar{\phi})}=\frac{4\pi
G}{c^{2}}\int^{\infty}_{r_{s}}r\delta\rho(r)dr\,,$ (A.4)

where the $\phi_{s}$ term has been dropped since it is assumed the screening
mechanism acts to force the field to smaller values inside high-density
objects, thus $\phi_{s}\ll\bar{\phi}$.

##  Appendix B Galaxy Pipeline

In this appendix we outline our procedure for deriving all five of our galaxy
parameters ($\rho_{\mathrm{NFW}},\ r_{\mathrm{NFW}},\ \Sigma_{\mathrm{disc}},\
R_{\mathrm{disc}}$ and $z_{\mathrm{disc}}$), introduced in section 3.1 via
eqs. (3.1) and (3.2), in terms of one input parameter, $M_{\mathrm{vir}}$. To
achieve this we take advantage of several empirical scaling relations. Note
that all non-analytic relations come with a related dispersion. As discussed
in the main text, section 3.1, we ignore the scatter about the scaling
relations, and simply adopt median values in each case.

###  B.1 NFW Parameters

On the dark matter side, an NFW profile eq. (3.1) is defined in terms of one
of two sets of input parameters, $\rho_{\mathrm{NFW}}$ and $r_{\mathrm{NFW}}$,
or $M_{\mathrm{vir}}$ and $c_{\mathrm{vir}}$. We choose the latter since we
already have one of the required input parameters, $M_{\mathrm{vir}}$,
equivalent to the commonly used $M_{200}$. As such, we only need to employ a
single empirical relation between the virial mass and concentration parameter.
For this, we use equation 8 of ref. [63],

$\log_{10}c_{\mathrm{vir}}=0.905-0.101\log_{10}\left(\frac{M_{\mathrm{vir}}}{10^{12}h^{-1}M_{\odot}}\right)\,,$
(B.1)

where $h$ is the Hubble parameter, taken to be $0.7$. This relation reflects
the well-known result from $\Lambda$CDM simulations that more massive dark
matter haloes are less centrally concentrated. The intrinsic scatter is given
as log-normal with $\sigma_{\log_{10}c_{\mathrm{vir}}}=0.11$, (see end of §4.1
in ref. [63]).

At this point, $\rho_{\mathrm{NFW}}$ and $r_{\mathrm{NFW}}$ can be obtained
from $M_{\mathrm{vir}}$ and $c_{\mathrm{vir}}$ from simultaneous solution of

$M_{\mathrm{vir}}=\frac{4}{3}\pi\Delta\rho_{\mathrm{crit}}R_{\mathrm{vir}}^{3}\,,$
(B.2) $R_{\mathrm{vir}}=c_{\mathrm{vir}}r_{\mathrm{NFW}}\,,$ (B.3)
$M_{\mathrm{vir}}=4\pi\rho_{\mathrm{NFW}}r_{\mathrm{NFW}}^{3}\left[\ln(1+c_{\mathrm{vir}})-\frac{c_{\mathrm{vir}}}{1+c_{\mathrm{vir}}}\right]\,,$
(B.4)

which derive respectively from the definition of the virial mass, the
definition of the concentration, and the relation for the enclosed mass at a
given radius under the NFW profile.

###  B.2 Galactic Disc Parameters

To derive the three parameters ($\Sigma_{\mathrm{disc}}$, $R_{\mathrm{disc}}$
and $z_{\mathrm{disc}}$) needed to construct the galactic disc profile, eq.
(3.2), we must pass from the virial mass to the stellar mass (contributing to
the bulk of the galactic disc). To go from halo to disc mass, we employ the
empirical stellar-halo mass relation (SHMR), namely eq. 2 in ref. [64]

$M_{\mathrm{disc}}=2NM_{\mathrm{halo}}\left[\left(\frac{M_{\mathrm{halo}}}{M_{1}}\right)^{-\beta}+\left(\frac{M_{\mathrm{halo}}}{M_{1}}\right)^{\gamma}\right]^{-1}\,,$
(B.5)

where the four fitting parameters and related errors are in table 3 of ref.
[64].

Now that we have the disc mass, we pass via the half-light radius to the disc
scale length, $R_{\mathrm{disc}}$. To achieve this, we use the observational
relation defined by eq. 18 in ref. [65]. Specifically,

$\displaystyle\frac{\overline{R}_{hl}}{\mathrm{kpc}}=\gamma\left(\frac{M_{\mathrm{disc}}}{M_{\odot}}\right)^{\alpha}\left(1+\frac{M_{\mathrm{disc}}}{M_{0}}\right)^{\beta-\alpha}\,,$
(B.6)

where $\overline{R}_{hl}$ denotes the median half-light radius and fitting
parameters ($\alpha$, $\beta$, $\gamma$, $M_{0}$, $\sigma_{1}$ and
$\sigma_{2}$) take the values from the bottom row of table 1 in ref. [65].
Then, by integrating the density, eq. (3.2), to calculate the mass, and
assuming a constant mass-to-light ratio so that we can utilise
$M_{\mathrm{disc}}(R=\overline{R}_{hl})=\frac{1}{2}M_{\mathrm{disc}}(R\rightarrow\infty)$,
we arrive at the relation

$R_{\mathrm{disc}}=-\overline{R}_{hl}\left(1+W_{-1}\left(-\frac{1}{2e}\right)\right)\approx
0.595824\overline{R}_{hl}\,,$ (B.7)

where $W_{-1}(x)$ is the $n=-1$ solution to the Lambert $W_{n}$ function.

The disc mass can be determined by integrating the density profile, eq. (3.2),

$M_{\mathrm{disc}}(R\rightarrow\infty)=2\pi\Sigma_{\mathrm{disc}}R_{\mathrm{disc}}^{2}\,,$
(B.8)

at which point it is easily reversed and combined with our value of
$R_{\mathrm{disc}}$ to calculate $\Sigma_{\mathrm{disc}}$.

Finally, we employ one last observational relation to determine the disc scale
height, $z_{\mathrm{disc}}$, in terms of the disc scale radius
$R_{\mathrm{disc}}$. For this we utilise eq. 1 of ref. [66],

$\log_{10}\left(\frac{R_{\mathrm{disc}}}{z_{\mathrm{disc}}}\right)=0.367\log_{10}\left(\frac{R_{\mathrm{disc}}}{\mathrm{kpc}}\right)+0.708\pm
0.095\,.$ (B.9)

With this, we have defined all five density and length scales
($\rho_{\mathrm{NFW}},\ r_{\mathrm{NFW}},\ \Sigma_{\mathrm{disc}},\
R_{\mathrm{disc}}$ and $z_{\mathrm{disc}}$) in terms of a single input,
$M_{\mathrm{vir}}$.

## Acknowledgments

C.B. is supported by a Research Leadership Award from the Leverhulme Trust and
by the STFC under grant ST/T000732/1. B.M. is supported by an STFC
studentship. A.P.N. is supported by an Early Career Fellowship from the
Leverhulme Trust.

## References

  * [1] E.G. Adelberger, B.R. Heckel and A.E. Nelson, _Tests of the Gravitational Inverse-Square Law_ , _Annual Review of Nuclear and Particle Science_ 53 (2003) 77 [hep-ph/0307284].
  * [2] A. Joyce, B. Jain, J. Khoury and M. Trodden, _Beyond the Cosmological Standard Model_ , _Phys. Rept._ 568 (2015) 1 [1407.0059].
  * [3] K. Koyama, _Cosmological Tests of Modified Gravity_ , _Rept. Prog. Phys._ 79 (2016) 046902 [1504.04623].
  * [4] M. Ishak, _Testing General Relativity in Cosmology_ , _Living Rev. Rel._ 22 (2019) 1 [1806.10122].
  * [5] C. Burrage and J. Sakstein, _Tests of chameleon gravity_ , _Living Reviews in Relativity_ 21 (2018) 1 [1709.09071].
  * [6] P. Brax, S. Casas, H. Desmond and B. Elder, _Testing Screened Modified Gravity_ , _Universe_ 8 (2021) 11 [2201.10817].
  * [7] E.N. Saridakis, R. Lazkoz, V. Salzano, P. Vargas Moniz, S. Capozziello, J. Beltrán Jiménez et al., _Modified Gravity and Cosmology: An Update by the CANTATA Network_ (May, 2021), 10.48550/arXiv.2105.12582, [2105.12582].
  * [8] V. Vardanyan and D.J. Bartlett, _Modeling and Testing Screening Mechanisms in the Laboratory and in Space_ , _Universe_ 9 (2023) 340 [2305.18899].
  * [9] J. Khoury and A. Weltman, _Chameleon fields: Awaiting surprises for tests of gravity in space_ , _Phys. Rev. Lett._ 93 (2004) 171104 [astro-ph/0309300].
  * [10] J. Khoury and A. Weltman, _Chameleon cosmology_ , _Phys. Rev. D_ 69 (2004) 044026 [astro-ph/0309411].
  * [11] K. Hinterbichler and J. Khoury, _Screening Long-Range Forces through Local Symmetry Restoration_ , _Phys. Rev. Lett._ 104 (2010) 231301 [1001.4525].
  * [12] K. Hinterbichler, J. Khoury, A. Levy and A. Matas, _Symmetron Cosmology_ , _Phys. Rev. D_ 84 (2011) 103521 [1107.2112].
  * [13] D.F. Mota and D.J. Shaw, _Evading Equivalence Principle Violations, Cosmological and other Experimental Constraints in Scalar Field Theories with a Strong Coupling to Matter_ , _Phys. Rev. D_ 75 (2007) 063501 [hep-ph/0608078].
  * [14] H. Dehnen, H. Frommert and F. Ghaboussi, _Higgs field and a new scalar - tensor theory of gravity_ , _Int. J. Theor. Phys._ 31 (1992) 109.
  * [15] E. Gessner, _A new scalar tensor theory for gravity and the flat rotation curves of spiral galaxies_ , _Astrophys. Space Sci._ 196 (1992) 29.
  * [16] T. Damour and A.M. Polyakov, _The String dilaton and a least coupling principle_ , _Nucl. Phys. B_ 423 (1994) 532 [hep-th/9401069].
  * [17] M. Pietroni, _Dark energy condensation_ , _Phys. Rev. D_ 72 (2005) 043535 [astro-ph/0505615].
  * [18] K.A. Olive and M. Pospelov, _Environmental dependence of masses and coupling constants_ , _Phys. Rev. D_ 77 (2008) 043524 [0709.3825].
  * [19] G. Dvali, G. Gabadadze and M. Porrati, _4D gravity on a brane in 5D Minkowski space_ , _Physics Letters B_ 485 (2000) 208 [hep-th/0005016].
  * [20] A.I. Vainshtein, _To the problem of nonvanishing gravitation mass_ , _Phys. Lett. B_ 39 (1972) 393.
  * [21] L. Hui, A. Nicolis and C.W. Stubbs, _Equivalence principle implications of modified gravity models_ , _Phys. Rev. D_ 80 (2009) 104002 [0905.2966].
  * [22] B. Jain and J. VanderPlas, _Tests of modified gravity with dwarf galaxies_ , _JCAP_ 2011 (2011) 032 [1106.0065].
  * [23] D.J. Bartlett, H. Desmond and P.G. Ferreira, _Calibrating galaxy formation effects in galactic tests of fundamental physics_ , _Phys. Rev. D_ 103 (2021) 123502 [2103.10356].
  * [24] H. Desmond, P.G. Ferreira, G. Lavaux and J. Jasche, _Fifth force constraints from the separation of galaxy mass components_ , _Phys. Rev. D_ 98 (2018) 064015 [1807.01482].
  * [25] H. Desmond, P.G. Ferreira, G. Lavaux and J. Jasche, _Fifth force constraints from galaxy warps_ , _Phys. Rev. D_ 98 (2018) 083010 [1807.11742].
  * [26] H. Desmond, P.G. Ferreira, G. Lavaux and J. Jasche, _The fifth force in the local cosmic web_ , _Mon. Not. Royal Astron. Soc._ 483 (2019) L64 [1802.07206].
  * [27] H. Desmond and P.G. Ferreira, _Galaxy morphology rules out astrophysically relevant Hu-Sawicki f (R ) gravity_ , _Phys. Rev. D_ 102 (2020) 104060 [2009.08743].
  * [28] A.P. Naik, N.W. Evans, E. Puchwein, H. Zhao and A.C. Davis, _Stellar streams in chameleon gravity_ , _Phys. Rev. D_ 102 (2020) 084066 [2002.05738].
  * [29] V. Vikram, A. Cabré, B. Jain and J.T. VanderPlas, _Astrophysical tests of modified gravity: the morphology and kinematics of dwarf galaxies_ , _JCAP_ 2013 (2013) 020 [1303.0295].
  * [30] E.M. Pedersen and C.W. Stubbs, _Using Elliptical Galaxy Kinematics to Compare of the Strength of Gravity in Cosmological Regions of Differing Gravitational Potential – A First Look_ , 2304.02123.
  * [31] A.P. Naik, E. Puchwein, A.-C. Davis and C. Arnold, _Imprints of Chameleon f(R) gravity on Galaxy rotation curves_ , _Mon. Not. Royal Astron. Soc._ 480 (2018) 5211 [1805.12221].
  * [32] A.P. Naik, E. Puchwein, A.-C. Davis, D. Sijacki and H. Desmond, _Constraints on chameleon f(R)-gravity from galaxy rotation curves of the SPARC sample_ , _Mon. Not. Royal Astron. Soc._ 489 (2019) 771 [1905.13330].
  * [33] V. Vikram, J. Sakstein, C. Davis and A. Neil, _Astrophysical Tests of Modified Gravity: Stellar and Gaseous Rotation Curves in Dwarf Galaxies_ , _arXiv e-prints_ (2014) arXiv:1407.6044 [1407.6044].
  * [34] L. Lombriser and J. Peñarrubia, _How chameleons core dwarfs with cusps_ , _Phys. Rev. D_ 91 (2015) 084022 [1407.7862].
  * [35] C. Burrage, E.J. Copeland and P. Millington, _Radial acceleration relation from symmetron fifth forces_ , _Phys. Rev. D_ 95 (2017) 064050 [1610.07529].
  * [36] P. Cataldi, S. Pedrosa, N. Padilla, S. Landau, C. Arnold and B. Li, _Fingerprints of modified gravity on galaxies in voids_ , _Mon. Not. Royal Astron. Soc._ 515 (2022) 5358 [2207.12917].
  * [37] C.A.J. O’Hare and C. Burrage, _Stellar kinematics from the symmetron fifth force in the Milky Way disk_ , _Phys. Rev. D_ 98 (2018) 064019 [1805.05226].
  * [38] T. Baker, A. Barreira, H. Desmond, P. Ferreira, B. Jain, K. Koyama et al., _The Novel Probes Project – Tests of Gravity on Astrophysical Scales_ , _arXiv e-prints_ (2019) arXiv:1908.03430 [1908.03430].
  * [39] E. Puchwein, M. Baldi and V. Springel, _Modified-Gravity-GADGET: a new code for cosmological hydrodynamical simulations of modified gravity models_ , _Mon. Not. Royal Astron. Soc._ 436 (2013) 348 [1305.2418].
  * [40] S. Bose, B. Li, A. Barreira, J.-h. He, W.A. Hellwing, K. Koyama et al., _Speeding up N-body simulations of modified gravity: chameleon screening models_ , _JCAP_ 2017 (2017) 050 [1611.09375].
  * [41] L. Hui, A. Nicolis and C.W. Stubbs, _Equivalence principle implications of modified gravity models_ , _Phys. Rev. D_ 80 (2009) 104002.
  * [42] A.-C. Davis, E.A. Lim, J. Sakstein and D.J. Shaw, _Modified gravity makes galaxies brighter_ , _Phys. Rev. D_ 85 (2012) 123006.
  * [43] J. Sakstein, _Stellar oscillations in modified gravity_ , _Phys. Rev. D_ 88 (2013) 124013.
  * [44] C. Burrage and J. Sakstein, _A compendium of chameleon constraints_ , _JCAP_ 2016 (2016) 045 [1609.01192].
  * [45] C. Arnold, V. Springel and E. Puchwein, _Zoomed cosmological simulations of Milky Way-sized haloes in f(R) gravity_ , _Mon. Not. Royal Astron. Soc._ 462 (2016) 1530 [1604.06095].
  * [46] A. Cabré, V. Vikram, G.-B. Zhao, B. Jain and K. Koyama, _Astrophysical tests of gravity: a screening map of the nearby universe_ , _JCAP_ 2012 (2012) 034 [1204.6046].
  * [47] H. Desmond, P.G. Ferreira, G. Lavaux and J. Jasche, _Reconstructing the gravitational field of the local Universe_ , _Mon. Not. Royal Astron. Soc._ 474 (2018) 3152 [1705.02420].
  * [48] T.P. Sotiriou and V. Faraoni, _f(R) theories of gravity_ , _Reviews of Modern Physics_ 82 (2010) 451 [0805.1726].
  * [49] P. Brax, C. van de Bruck, A.-C. Davis and D.J. Shaw, _f(R) gravity and chameleon theories_ , _Phys. Rev. D_ 78 (2008) 104021 [0806.3415].
  * [50] T. Faulkner, M. Tegmark, E.F. Bunn and Y. Mao, _Constraining f(R) gravity as a scalar-tensor theory_ , _Phys. Rev. D_ 76 (2007) 063505 [astro-ph/0612569].
  * [51] A.A. Starobinsky, _Disappearing cosmological constant in f( R) gravity_ , _Soviet Journal of Experimental and Theoretical Physics Letters_ 86 (2007) 157 [0706.2041].
  * [52] I. Navarro and K. Van Acoleyen, _f(R) actions, cosmic acceleration and local tests of gravity_ , _JCAP_ 2007 (2007) 022 [gr-qc/0611127].
  * [53] W. Hu and I. Sawicki, _Models of f(R) cosmic acceleration that evade solar system tests_ , _Phys. Rev. D_ 76 (2007) 064004 [0705.1158].
  * [54] J.F. Navarro, C.S. Frenk and S.D.M. White, _The Structure of Cold Dark Matter Halos_ , _Astrophys. J._ 462 (1996) 563 [astro-ph/9508025].
  * [55] J. Binney and S. Tremaine, _Galactic Dynamics: Second Edition_ , Princeton University Press, Princeton, NJ USA (2008).
  * [56] P.J. McMillan, _The mass distribution and gravitational potential of the Milky Way_ , _Mon. Not. Royal Astron. Soc._ 465 (2017) 76 [1608.00971].
  * [57] A.J. Deason, A. Fattahi, C.S. Frenk, R.J.J. Grand, K.A. Oman, S. Garrison-Kimmel et al., _The edge of the Galaxy_ , _Mon. Not. Royal Astron. Soc._ 496 (2020) 3929 [2002.09497].
  * [58] E.G. Adelberger, J.H. Gundlach, B.R. Heckel, S. Hoedl and S. Schlamminger, _Torsion balance experiments: A low-energy frontier of particle physics_ , _Prog. Part. Nucl. Phys._ 62 (2009) 102.
  * [59] T.M. Callingham, M. Cautun, A.J. Deason, C.S. Frenk, W. Wang, F.A. Gómez et al., _The mass of the Milky Way from satellite dynamics_ , _Mon. Not. Royal Astron. Soc._ 484 (2019) 5453 [1808.10456].
  * [60] M. Kesden and M. Kamionkowski, _Galilean Equivalence for Galactic Dark Matter_ , _Phys. Rev. Lett._ 97 (2006) 131303 [astro-ph/0606566].
  * [61] M. Kesden and M. Kamionkowski, _Tidal tails test the equivalence principle in the dark-matter sector_ , _Phys. Rev. D_ 74 (2006) 083007 [astro-ph/0608095].
  * [62] J.A. Keselman, A. Nusser and P.J.E. Peebles, _Galaxy satellites and the weak equivalence principle_ , _Phys. Rev. D_ 80 (2009) 063517 [0902.3452].
  * [63] A.A. Dutton and A.V. Macciò, _Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles_ , _Mon. Not. Royal Astron. Soc._ 441 (2014) 3359 [1402.7073].
  * [64] B.P. Moster, T. Naab and S.D.M. White, _Galactic star formation and accretion histories from matching galaxies to dark matter haloes_ , _Mon. Not. Royal Astron. Soc._ 428 (2013) 3121 [1205.5807].
  * [65] S. Shen, H.J. Mo, S.D.M. White, M.R. Blanton, G. Kauffmann, W. Voges et al., _The size distribution of galaxies in the Sloan Digital Sky Survey_ , _Mon. Not. Royal Astron. Soc._ 343 (2003) 978 [astro-ph/0301527].
  * [66] M.A. Bershady, M.A.W. Verheijen, K.B. Westfall, D.R. Andersen, R.A. Swaters and T. Martinsson, _The DiskMass Survey. II. Error Budget_ , _Astrophys. J._ 716 (2010) 234 [1004.5043].