Constraints on Chameleon f(R)-Gravity from Galaxy Rotation Curves of the SPARC Sample

TL;DR

This study tests chameleon $f(R)$-gravity using galaxy rotation curves from the SPARC sample, finding no convincing evidence for $f(R)$ effects at current sensitivity levels and constraining the scalar field parameter $f_{R0}$.

Contribution

It provides new constraints on $f(R)$-gravity from galaxy rotation curves, highlighting the degeneracy with core/cusp issues and comparing different halo profiles.

Findings

Models with $oxed{ ext{log}_{10}|f_{R0}| > -6.1}$ are excluded.

Favored $f(R)$ range is $-7.5< ext{log}_{10}|f_{R0}|<-6.5$, peaking at -7.

No convincing evidence of $f(R)$ gravity at $|f_{R0}| o 6 imes 10^{-8}$.

Abstract

In chameleon $f(R)$-gravity, the fifth force will lead to `upturns' in galaxy rotation curves near the screening radius. The location of the upturn depends on the cosmic background value of the scalar field $f_{R0}$, as well as the mass, size and environment of the galaxy. We search for this signature of modified gravity in the SPARC sample of measured rotation curves, using an MCMC technique to derive constraints on $f_{R0}$. Assuming NFW dark matter haloes and with $f_{R0}$ freely varying for each galaxy, most galaxies prefer $f(R)$ gravity to $\Lambda$CDM, but there is a large spread of inferred $f_{R0}$ values, inconsistent with a single global value. Requiring instead a consistent $f_{R0}$ value for the whole sample, models with $\log_{10}|f_{R0}| > -6.1$ are excluded. On the other hand, models in the range $-7.5<\log_{10}|f_{R0}|<-6.5$ seem to be favoured with respect to $\Lambda$CDM, with a significant peak at -7. However, this signal is largely a result of galaxies for which the $f(R)$ signal is degenerate with the core/cusp problem, and when the NFW profile is replaced with a cored halo profile, $\Lambda$CDM gives better fits than any given $f(R)$ model. Thus, we find no convincing evidence of $f(R)$ gravity down to the level of $|f_{R0}| \sim 6 \times 10^{-8}$, with the caveat that if cored halo density profiles cannot ultimately be explained within $\Lambda$CDM, a screened modified gravity theory could possibly provide an alternative solution for the core/cusp problem. However, the $f(R)$ models studied here fall short of achieving this.

Figures

Tables

Table 1. Table 1: Summary of the models investigated in this work. Columns are (1) Model: the alphabetical/alphanumerical ID given in this work to a given model, (2) Theory: Λ Λ \Lambda CDM or f ​ ( R ) 𝑓 𝑅 f(R) , i.e. whether or not a fifth force contribution is included in the rotation curve model, (3) f ¯ R0 subscript ¯ 𝑓 R0 \bar{f}_{\mathrm{R0}} : for an f ​ ( R ) 𝑓 𝑅 f(R) model, f ¯ R0 subscript ¯ 𝑓 R0 \bar{f}_{\mathrm{R0}} is either left as a freely ‘varying’ parameter, or in the case of Models H0-19 and I0-19, imposed as a series of 20 fixed values given in Table 2 (see discussion in § 4.1 ), (4) DM Halo:, NFW or DC14 dark matter halo profile (§ 4.1 ), (5) Large-scale Env.: whether a large-scale environment is added to the scalar field solver (§ 4.3 ). Note that in all models, we take into account environmental screening effects with the final sample cut discussed in § 3 , and this additional test is merely to ensure that the remaining sample is largely unaffected by environmental screening. (6) Stellar Screening: whether the stellar component is excluded from sourcing a fifth force in the scalar field solver (§ 4.4 ), (7) Υ Υ \Upsilon : Whether the mass-to-light ratio is treated with one free parameter (using the same ratio for disc and bulge),‘1’, two free parameters (allowing for different values for disc and bulge),‘2’, or as ‘fixed’ empirical values,‘0’ (§ 4.1 ), and (8) Parameters: a list of the free parameters for a given model (see main text for details).

Model	Theory	${\overline{f}}_{\mathrm{R0}}$	DM Halo	Large-scale Env.	Stellar Screening	${\rm Y}$	Parameters
A	$\Lambda$CDM	-	NFW	-	-	1	$v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
B	$f\left(R\right)$	Varying	NFW	$\times$	$\times$	1	${\overline{f}}_{\mathrm{R0}}$, $v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
C	$f\left(R\right)$	Varying	NFW	✓	$\times$	1	${\overline{f}}_{\mathrm{R0}}$, $v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
D	$f\left(R\right)$	Varying	NFW	$\times$	✓	1	${\overline{f}}_{\mathrm{R0}}$, $v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
E	$f\left(R\right)$	Varying	NFW	$\times$	$\times$	0	${\overline{f}}_{\mathrm{R0}}$, $v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\sigma }_g$
F	$f\left(R\right)$	Varying	NFW	$\times$	$\times$	2	${\overline{f}}_{\mathrm{R0}}$, $v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${{\rm Y}}_{\mathrm{disc}}$, ${{\rm Y}}_{\mathrm{bulge}}$, ${\sigma }_g$
G	$\Lambda$CDM	-	DC14	-	-	1	$v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
H0-19	$f\left(R\right)$	20 imposed values	NFW	$\times$	$\times$	1	$v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$
I0-19	$f\left(R\right)$	20 imposed values	DC14	$\times$	$\times$	1	$v_{\mathrm{vir}}$, $c_{\mathrm{vir}}$, ${\rm Y}$, ${\sigma }_g$

Table 2. Table 2: Values of log 10 ⁡ | f ¯ R0 | subscript 10 subscript ¯ 𝑓 R0 \log_{10}|\bar{f}_{\mathrm{R0}}| used in Models H0-19 and I0-19. They form a sequence of values that are uniformly spaced in the logarithm from log 10 ⁡ ( 1.563 × 10 − 8 ) subscript 10 1.563 superscript 10 8 \log_{10}\left(1.563\times 10^{-8}\right) to log 10 ⁡ ( 2 × 10 − 6 ) subscript 10 2 superscript 10 6 \log_{10}\left(2\times 10^{-6}\right) .

Model	${\log }_{10}\left|{\overline{f}}_{\mathrm{R0}}\right|$	Model	${\log }_{10}\left|{\overline{f}}_{\mathrm{R0}}\right|$
H0 (or I0)	-7.806	H10	-6.697
H1	-7.695	H11	-6.586
H2	-7.584	H12	-6.475
H3	-7.473	H13	-6.364
H4	-7.362	H14	-6.253
H5	-7.252	H15	-6.143
H6	-7.141	H16	-6.032
H7	-7.030	H17	-5.921
H8	-6.919	H18	-5.810
H9	-6.808	H19	-5.699

Equations

$$S=\int d^{4}x\sqrt{-g}\frac{1}{16\pi G}\left[R+f\left(R\right)\right]+S_m[g_{\mu \nu },{\psi }_i],$$

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

$$\delta R=\bar{R}(a)\left(\sqrt{\frac{{\bar{f}}_{\mathrm{R}}(a)}{f_{\mathrm{R}}}}-1\right),$$

$${\bar{f}}_{\mathrm{R}}(a)=a^6{\bar{f}}_{R0}{\left(\frac{1+4\frac{{\Omega }_{\Lambda }}{{\Omega }_m}}{1+4\frac{{\Omega }_{\Lambda }{a}^3}{{\Omega }_m}}\right)}^2,$$

$$\bar{R}(a)=3\frac{H_0^2{\Omega }_m}{a^3}\left(1+4\frac{{\Omega }_{\Lambda }{a}^3}{{\Omega }_m}\right).$$

$${\mathbf{a}}_5=\frac{1}{2}{c}^2\nabla f_{\mathrm{R}}.$$

$$v_{\mathrm{model}}^2(r)=a_{\mathrm{tot}}r=(a_{\mathrm{N}}+a_5)r,$$

$$a_{\mathrm{N}}(r)=\frac{v_{\mathrm{gas}}^2+v_{\mathrm{DM}}^2+{{\rm Y}}_{\mathrm{disc}}{v}_{\mathrm{disc}}^2+{{\rm Y}}_{\mathrm{bulge}}{v}_{\mathrm{bulge}}^2}{r},$$

$$c_{\mathrm{vir}}=\frac{r_{\mathrm{vir}}}{r_{-2}},$$

$$v_{\mathrm{vir}}=\sqrt{\frac{G{M}_{\mathrm{vir}}}{r_{\mathrm{vir}}}},$$

$${\nabla }^2{e}^u+\frac{1}{3c^2{\bar{f}}_R(a)}\left[\bar{R}\left(a\right)\left(1-e^{-\frac{u}{2}}\right)+8\pi G\delta \rho \right]=0.$$

$${\mathcal{L}}_i\equiv {\left({\nabla }^2{e}^u\right)}_i+\frac{1}{3c^2{\bar{f}}_R(a)}\left[\bar{R}\left(a\right)\left(1-e^{-\frac{u_i}{2}}\right)+8\pi G\delta {\rho }_i\right],$$

$${\mathcal{L}}_i=0.$$

$$u_i^{n+1}=u_i^n-\frac{{\mathcal{L}}_i^n}{\frac{\partial {\mathcal{L}}_i^n}{\partial u_i^n}}.$$

$${\nabla }^{2}f=\frac{1}{r^3}\frac{\partial }{\partial x}\left(r\frac{\partial f}{\partial x}\right)$$

$${\left({\nabla }^{2}f\right)}_i=\frac{1}{r_i^3{h}_x^2}\left(r_{i+\frac{1}{2}}(f_{i+1}-f_i)-r_{i-\frac{1}{2}}(f_i-f_{i-1})\right),$$

$$\frac{\partial {\mathcal{L}}_i}{\partial u_i}=\frac{\bar{R}(a)}{6c^2{\bar{f}}_R(a)}{e}^{-\frac{u_i}{2}}-e^{u_i}\left(\frac{r_{i+\frac{1}{2}}+r_{i-\frac{1}{2}}}{r_i^3{h}_x^2}\right).$$

$${\lambda }_C\approx 32\sqrt{\frac{|{\bar{f}}_{R0}|}{10^{-4}}}\mathrm{Mpc}.$$

$$\mathcal{P}(\theta |\mathcal{D})\propto \mathcal{L}(D|\theta )\pi (\theta ),$$

$$\ln \mathcal{L}=-\frac{1}{2}\sum _j\left[{\left(\frac{v_{{\mathrm{data}}_j}-v_{\mathrm{model}}(r_j)}{{\sigma }_j}\right)}^2+\ln \left(2\pi {\sigma }_j^2\right)\right],$$

$${\sigma }_j^2={\sigma }_{\mathrm{obs},j}^2+{\sigma }_g^2.$$

$$v_5(r)=\sqrt{a_5(r)r}$$

$$\mathrm{BIC}=\ln (n)k-2\ln (\mathcal{L}),$$

$${\mathcal{L}}_{ij}\equiv {\left({\nabla }^2{e}^u\right)}_{ij}+\frac{1}{3c^2{\bar{f}}_R(a)}\left[\bar{R}\left(a\right)\left(1-e^{-\frac{u_{ij}}{2}}\right)+8\pi G\delta {\rho }_{ij}\right],$$

$${\mathcal{L}}_{ij}=0.$$

$${\nabla }^{2}f=\frac{1}{r^3}\frac{\partial }{\partial x}\left(r\frac{\partial f}{\partial x}\right)+\frac{1}{r^{2}s}\frac{\partial }{\partial \theta }\left(s\frac{\partial f}{\partial \theta }\right),$$

$$\begin{gathered} {\left({\nabla }^{2}f\right)}_{ij}=\frac{1}{r_i^3{h}_x^2}\left(r_{i+\frac{1}{2}}(f_{i+1,j}-f_{ij})-r_{i-\frac{1}{2}}(f_{ij}-f_{i-1,j})\right) \\ +\frac{1}{r_i^2{s}_j{h}_{\theta }^2}\left(s_{j+\frac{1}{2}}(f_{i,j+1}-f_{ij})-s_{j-\frac{1}{2}}(f_{ij}-f_{i,j-1})\right). \end{gathered}$$

$$\begin{gathered} \frac{\partial {\mathcal{L}}_{ij}}{\partial u_{ij}}=\frac{\bar{R}}{6c^2{\bar{f}}_R}{e}^{-\frac{u_{ij}}{2}} \\ -e^{u_{ij}}\left[\frac{r_{i+\frac{1}{2}}+r_{i-\frac{1}{2}}}{r_i^3{h}_x^2}+\frac{s_{j+\frac{1}{2}}+s_{j-\frac{1}{2}}}{r_i^2{s}_j{h}_{\theta }^2}\right]. \end{gathered}$$

Full text

††pubyear: 2019††pagerange: Constraints on Chameleon $f(R)$-Gravity from Galaxy Rotation Curves of the SPARC Sample–A

1 Introduction

In the previous two decades, the $\Lambda$CDM cosmological model has proven remarkably successful at describing the large-scale structure of the Universe (Planck Collaboration, 2016). Despite this success, there remain some outstanding theoretical questions, such as the precise nature of dark matter, and the large discrepancy between (naively) expected and observed values of the cosmological constant $\Lambda$ (Bull et al., 2016).

In addition to these theoretical issues, there also remain a few unresolved tensions with observations on galactic scales, a review of which can be found in the article by Bullock & Boylan-Kolchin (2017). There are three problems in particular that have been the subjects of widespread attention. Firstly, the ‘missing satellites’ problem (Moore et al., 1999; Klypin et al., 1999), which arises from a mismatch between the observed number of satellite galaxies of the Milky Way and that expected from the number of substructures in $\Lambda$CDM cosmological simulations. Secondly, the ‘too-big-to-fail’ problem, which arises from the unexpectedly low central densities of observed Milky Way satellites (Boylan-Kolchin et al., 2011).

The (likely related) third problem, the so-called ‘core/cusp’ problem, is of particular relevance here. A key prediction of dark matter-only $\Lambda$CDM simulations is that of ‘cuspy’ density profiles for dark matter haloes, i.e. profiles with steep central slopes (Navarro et al., 1997), roughly $\propto r^{-1}$. On the other hand, some observations (Moore, 1994; Flores & Primack, 1994; Walker & Peñarrubia, 2011; Oh et al., 2015), e.g., of galaxy rotation curves, have been found to suggest flatter, ‘cored’ profiles in the inner regions of haloes. As a possible resolution, some simulations including the effects of baryonic feedback have suggested that such mechanisms can induce cored central density distributions in dark matter haloes (Mashchenko et al., 2008; Governato et al., 2010; Pontzen & Governato, 2012; Di Cintio et al., 2014). This is thought to happen when feedback, e.g., from supernovae pushes around large amounts of gas in the central region of a galaxy, the resulting fluctuation in the gravitational potential of the gas could then also affect the dark matter distribution. Other simulations, however, find no such effect (e.g, Bose et al., 2019). In addition, the recent work of Genina et al. (2018) found that the inference of a cored profile from observations might be the result of an incorrect assumption of spherical symmetry, and that observed galaxies are consistent with cusps. The core/cusp debate hence remains far from resolved.

Because of these theoretical and observational tensions, alternatives and extensions to $\Lambda$CDM warrant investigation. There have been many works investigating viable substitutes for cold dark matter, such as ‘self-interacting dark matter’ which has been shown to produce cored haloes in simulations (Vogelsberger et al., 2012). Alternatively, a generic and widely-adopted method of extending $\Lambda$CDM is through the introduction of a scalar field coupled to gravity, giving rise to gravitational-strength couplings to matter: so-called ‘fifth forces’. Overviews of these modified gravity theories can be found in, e.g. the works by Amendola & Tsujikawa (2010); Clifton et al. (2012); Joyce et al. (2015); Koyama (2016).

For the scalar fields to retain cosmological relevance, they need sufficiently low masses (i.e. large Compton wavelengths), but their fifth forces will consequently violate stringent Solar System tests of gravity unless a ‘screening mechanism’ is introduced (Jain & Khoury, 2010; Khoury, 2010). Most relevant to the present work is the chameleon mechanism (Khoury & Weltman, 2004), in which the mass of the scalaron increases in dense environments, where the fifth force is consequently rendered undetectable.

The theory investigated in this work falls within the widely-studied class of $f(R)$ gravity models. First described in Buchdahl (1970), the idea behind $f(R)$ gravity is simple: the Ricci scalar $R$ in the Einstein-Hilbert action is generalised to $R+f(R)$. The Hu-Sawicki form of $f(R)$ (Hu & Sawicki, 2007) has been shown to give a theory that is formally equivalent to a subclass of chameleon theories (Brax et al., 2008). The theory has one key parameter, the cosmic background value of the scalar field, $\bar{f}_{\mathrm{R0}}$, constraints on which can therefore be translated into constraints on generalised chameleon parameters.

It is worth noting that the significant recent ‘multi-messenger’ observation of GW170817 placed severe constraints on many modified gravity theories, which posited differing speeds for gravitational and electromagnetic waves. $f(R)$ gravity, however, is consistent with such constraints, predicting a gravitational wave speed equal to that of light (Lombriser & Lima, 2017; Sakstein & Jain, 2017; Ezquiaga & Zumalacárregui, 2017).

Chameleon gravity has been tested using a multitude of methods, on a range of length scales from the laboratory to the CMB. Review articles summarising these tests and corresponding constraints can be found in the works by Lombriser (2014); Burrage & Sakstein (2016, 2018). In particular, astrophysical tests have been found to give some of the strongest constraints on $f(R)$ gravity to date. For instance, distances inferred using Cepheid variables and stars at the tip of the Red Giant Branch would be different under modified gravity than those in $\Lambda$CDM, and can therefore be used to constrain modified gravity theories. Using such a method, Jain et al. (2013) find an upper bound of $|\bar{f}_{\mathrm{R0}}|=5\times 10^{-7}$. As another example, constraints at around the $10^{-6}$ level were found by Xu (2015), considering the effect of $f(R)$ gravity on redshift-space distortions. A similar upper bound was calculated by Vikram et al. (2018), who searched for evidence of significant differences in the gaseous and stellar rotation curves of isolated dwarf galaxies; the self-screening of stars would mean that in contrast to the gas, there would be no enhancement of the stellar rotation velocity due to the fifth force.

Vikram et al. (2013) look for signatures of displacement between galactic stellar and gaseous components, and signatures of warping as a result of this displacement, but find no significant constraints on modified gravity. Employing similar principles, Desmond et al. (2019, 2018b) use updated methods and data to search for these stellar-gas offsets, calculating an upper bound on $|\bar{f}_{\mathrm{R0}}|$ of a few $\times 10^{-8}$. Desmond et al. (2018a), meanwhile, searches for the warping signature. Interestingly, both of these studies independently find a possible fifth force signal. In the language of these studies, which vary the fractional change in effective Newton’s constant $\Delta G/G_{N}$ and Compton wavelength $\lambda_{C}$, there is a signal at $\Delta G/G_{N}\approx 10^{-2}$ and $\lambda_{C}\approx 2$ Mpc. This is irrelevant to $f(R)$, which predicts $\Delta G/G_{N}=1/3$ in unscreened regions, but interesting nonetheless. However, the authors add the cautionary note that the signal could well be a result of a number of other effects, including unaccounted-for galaxy formation physics.

In Naik et al. (2018)—henceforth Paper I—it was found that the presence of a chameleon fifth force would lead to observable imprints on galactic rotation curves. Specifically, there will be a ‘screening radius’ beyond which the fifth force is active. Within the radius, standard gravity is restored. This has the consequence that an ‘upturn’ appears in the rotation curve at this screening radius, with a fifth force enchancement beyond the screening radius. The location of this screening radius within a given galaxy depends on the global value of $\bar{f}_{\mathrm{R0}}$ and the environment of the galaxy, as well as its mass and size.

In this work, we use the SPARC sample of high-quality rotation curve measurements (Lelli et al., 2016) to search for this signature. After several cuts to the sample, 85 galaxies remain, each of which we model with a gaseous disc, a stellar disc, a dark matter halo, and where appropriate, a stellar bulge. Using these components, we solve the $f(R)$ scalar field equations for the fifth force contribution, thereby constructing a forward model for the galaxy’s rotation curve. With these models, we then employ an MCMC technique to explore the posterior probabilities for the parameters for each galaxy, looking in particular at the inferred values for the background value of the scalar field $\bar{f}_{\mathrm{R0}}$. We also perform model comparisons between $f(R)$ and $\Lambda$CDM using the likelihood ratios calculated for our best fit models.

Several recent studies have worked with the SPARC sample, e.g. McGaugh et al. (2016); Desmond (2017a, b); Lelli et al. (2017); Li et al. (2018). Of particular relevance to the present work, however, is the study of Katz et al. (2017). In that work, the SPARC galaxies were used to investigate the core/cusp problem. It was found in that work that fits to the rotation curves typically improved when cuspy ‘NFW’ haloes (Navarro et al., 1997) were replaced with the cored ‘DC14’ haloes of Di Cintio et al. (2014), derived empirically from $\Lambda$CDM simulations incorporating subgrid stellar feedback mechanisms.

We argued in Paper I that the NFW profile provides not only a good fit to dark matter haloes in $\Lambda$CDM simulations, but also in $f(R)$ simulations. In fact in $f(R)$ gravity, depending on halo mass, a higher concentration of the NFW profile is needed (e.g. Arnold et al., 2019). Thus, $f(R)$ gravity does not reduce the central densities of haloes and consequently does not offer a direct solution to the core/cusp problem. However, Lombriser & Peñarrubia (2015) made an interesting observation: a chameleon fifth force in the outer regions of a cuspy halo can mimic the observational signature of a cored halo profile in measurements of the kinematics of galaxies, we will discuss this possibility further in §5 and §6.

Baryonic feedback mechanisms that could in principle lead to the formation of a halo core, would likely be as able to do so under modified gravity as under $\Lambda$CDM. We, therefore explore the impact of a cored density profile by computing results both for NFW and DC14 haloes (similar to Katz et al., 2017). We utilise a model in which $\bar{f}_{\mathrm{R0}}$ is allowed to vary freely for each galaxy, as well as a grid of models covering a range of $\bar{f}_{\mathrm{R0}}$ values fixed for the whole sample. In the former case, we analyse the marginal posteriors on $\bar{f}_{\mathrm{R0}}$ for each galaxy to find any evidence of significant clustering around a single preferred value. In the latter case, we are able to compare each model with $\Lambda$CDM, and thus search more directly for any evidence of a globally preferred value for $\bar{f}_{\mathrm{R0}}$.

We also explore several different prescriptions for the mass-to-light ratios, and the impacts of stellar and environmental screening.

All of the code, analysis tools, and plotting scripts used in this work have been made into a publicly available python 3 package111https://github.com/aneeshnaik/spam.

This work is structured as follows. Firstly, §2 briefly describes $f(R)$ gravity, providing the theoretical details relevant to this work. Next, §3 describes the SPARC data set and §4 describes the methodology of the present work. Results are then outlined in §5, before discussion and concluding remarks in §6.

2 $f(R)$ Gravity

In this section, we summarise the most salient theoretical details about $f(R)$ gravity. A more in-depth description of the relevant aspects can be found in Paper I and the references therein.

$f(R)$ gravity was first introduced in Buchdahl (1970), more recent reviews can be found in the works of Amendola & Tsujikawa (2010); Clifton et al. (2012); Joyce et al. (2015); Koyama (2016).

The action of this theory is given by:

[TABLE]

where $S_{m}$ is the matter action, governing the matter fields $\psi_{i}$. As the name of the theory would suggest, the Ricci scalar $R$ in the classical Einstein-Hilbert action is replaced with a generalised $R+f\left(R\right)$. The theory reduces to $\Lambda$CDM in the case that $f(R)=-2\Lambda$.

In $f(R)$ gravity, the derivative $f_{\mathrm{R}}\equiv\frac{\mathrm{d}f}{\mathrm{dR}}$ plays the role of a scalar degree of freedom of the theory, similar to a scalar field. $f_{\mathrm{R}}$ is governed by the equation

[TABLE]

which can be derived by extremising the action in Eq. (1) with respect to the metric, and making a quasistatic approximation (Sawicki & Bellini, 2015; Noller et al., 2014). Here $\delta\rho$ is the perturbation of the matter density from its background value, while $\delta R$ denotes the perturbation of the scalar curvature.

The specific $f(R)$ model we consider in this work is that of Hu & Sawicki (2007), which is known to exhibit a ‘chameleon’ screening mechanism, whereby gradients of the scalar field are suppressed in dense regions. Equivalently, this can be framed as the mass of the scalaron being environment dependent, and growing unobservably large in dense environments. In these ‘screened’ regions, the range of the fifth force plummets, rendering it unobservable, and thus allowing such theories to escape stringent laboratory and Solar System constraints. In this theory, $\delta R$ is given by

[TABLE]

where $\bar{f}_{\mathrm{R}}(a)$ and $\bar{R}(a)$ are the background values of the scalar field and scalar curvature respectively, at a given cosmic scale factor $a$. These are given in turn by

[TABLE]

[TABLE]

Equation (2), combined with (3), (4), and (5), can be used to solve for the scalar field $f_{R}$ everywhere. The only inputs needed are the mass distribution (via the density perturbation $\delta\rho$) and the cosmic background value of the scalar field $\bar{f}_{\mathrm{R0}}$.

The fifth force, i.e. the resultant deviation from standard gravity, is then given by gradients of the scalar field $f_{R}$. The corresponding acceleration reads

[TABLE]

In Paper I, it was found by post-processing simulations that for certain values of $\bar{f}_{\mathrm{R0}}$, galaxies exhibit screening radii, beyond which rotation curves would show enhancements due to this fifth force.

3 Data

For this investigation, we have used the rotation curves from the SPARC sample (Lelli et al., 2016). The SPARC sample consists of 175 high-quality HI/H$\alpha$ rotation curves, along with accompanying $3.6\mu m$ photometry from Spitzer. The sample is diverse, spanning 5 orders of magnitude in mass, and encompassing a variety of morphologies. $3.6\mu m$ surface brightness is believed to be a good tracer of stellar mass, so its pairing with rotation curve data enables detailed modelling of the different components of a given galaxy.

For each galaxy, the SPARC sample includes a wealth of information. Most relevant is $v_{\mathrm{circ}}(R)$, the rotation speed as a function of radius, along with its corresponding error bars. In addition to this, the photometry information is converted into contributions to the rotation speed due to the stellar disc and stellar bulge, $v_{\mathrm{disc}}$ and $v_{\mathrm{bulge}}$. The bulge component is only included for galaxies for which the photometry profile departs significantly from an exponential disc profile in the central regions, a minority of the overall sample (13 of our final 85 galaxies). The quantities $v_{\mathrm{disc}}$ and $v_{\mathrm{bulge}}$ are based on a mass-to-light ratio of $1\,M_{\odot}/L_{\odot}$, and therefore need re-scaling for any other assumed ratio. Finally, for each galaxy, gas surface brightness information has been converted into a gas contribution $v_{\mathrm{gas}}$. A more detailed description of the data and the derivation of the velocity components can be found in Lelli et al. (2016).

From the original sample of 175 galaxies, we make a series of cuts, yielding a final sample of 85 galaxies. The cuts are summarised below, along with numbers in parenthesis indicating how many galaxies were removed at each successive cut.

We include galaxies for which:

• •

There are $\geq$ 5 data points. $(175-4=171)$

• •

Inclination $i\geq 30^{\circ}.$ $(171-12=159)$

• •

Quality flag $Q=1,2$. This excludes galaxies with major asymmetries, for which $Q=3$. $(159-10=149)$

• •

Information about the HI gas distribution is available. This excludes galaxies for which the $v_{\mathrm{gas}}$ contribution can only be estimated approximately. $(149-2=147)$

• •

There is not a significant environmental screening effect. See §4.3 for a more detailed description of this cut. $(147-62=85)$

4 Methods

4.1 Rotation Curve Models

In this study, we have experimented with a range of different models for galaxy rotation curves, e.g., models with different dark matter halo profiles, fixed or varying $\bar{f}_{\mathrm{R0}}$, as well as different treatments of mass-to-light ratio, environmental effects and a potential self-screening of stars. A summary of the models investigated can be found in Table 1. The terms used in this table will be explained over the course of the following few subsections.

For each galaxy, we construct a circular velocity model $v_{\mathrm{model}}$, given by

[TABLE]

where $r$ is the distance in the disc plane from the centre of the galaxy, $a_{\mathrm{N}}$ is the standard gravity (Newtonian) acceleration, and $a_{5}$ is the additional acceleration due to the fifth force.

The Newtonian acceleration $a_{\mathrm{N}}$ is sourced by the various components of the galaxy: the gas, stellar disc, stellar bulge, and dark matter. This relationship can be written as

[TABLE]

where $v_{x}$ is the contribution of component $x$ to the overall velocity curve as a function of radius $r$, and $\Upsilon_{\mathrm{disc}}$ and $\Upsilon_{\mathrm{bulge}}$ are the mass-to-light ratios for the stellar disc and bulge. The gaseous and stellar contributions ($v_{\mathrm{gas}}$, $v_{\mathrm{disc}}$, and $v_{\mathrm{bulge}}$) have been calculated by the SPARC team and are included within the SPARC data. Note also that the bulge term only applies to galaxies for which the SPARC team have found evidence for a stellar bulge, and have included a corresponding $v_{\mathrm{bulge}}$ contribution in the data; this corresponds to 13 of the 85 galaxies in our final sample.

All that remains to complete the model is to calculate the dark matter contribution $v_{\mathrm{DM}}$, the mass-to-light ratios $\Upsilon$, and the acceleration due to the fifth force $a_{5}$.

Firstly, we calculate the dark matter contribution $v_{\mathrm{DM}}$ by modelling the dark matter halo as a 2-parameter profile. We consider two different choices for the profile: the widely used Navarro-Frenk-White (NFW) profile (Navarro et al., 1997), as well as the ‘DC14’ profile of Di Cintio et al. (2014). The latter profile is derived empirically from $\Lambda$CDM simulations incorporating baryonic feedback, and provides a more ‘cored’ profile than NFW for galaxies with a significant stellar component. It was found by Katz et al. (2017) to provide a better fit to the SPARC galaxies than NFW. As indicated in Table 1, Models A-F and H0-19 all use NFW haloes, while Models G and I0-19 employ DC14 haloes.

In either case, the profile is determined by 2 parameters, which we adopt as free parameters in our fit: the concentration parameter $c_{\mathrm{vir}}$ and the virial velocity $v_{\mathrm{vir}}$. $c_{\mathrm{vir}}$ is defined by

[TABLE]

where $r_{\mathrm{vir}}$ is the virial radius, the radius enclosing a region of average density 93.6 times the cosmic critical density, and $r_{-2}$ is the radius at which the logarithmic slope of the halo density profile is -2, which in the case of NFW is the same as the scale radius. The virial velocity, $v_{\mathrm{vir}}$ is in turn given by

[TABLE]

where $M_{\mathrm{vir}}$ is the virial mass: the mass contained within the virial radius, approximated here by the total dark matter mass contained within the virial radius. Details of the translation between these parameters and $v_{\mathrm{DM}}$ can be found in, e.g., the appendix of Katz et al. (2018).

Secondly, we have experimented with three approaches for the mass-to-light ratios $\Upsilon_{\mathrm{disc}}$ and $\Upsilon_{\mathrm{bulge}}$. In the most general case, we treat them as two free parameters for each galaxy. Alternatively, we treat them as a single free parameter (taking $\Upsilon_{\mathrm{disc}}=\Upsilon_{\mathrm{bulge}}$), mirroring the approach of Katz et al. (2017). The final option is to treat them as fixed values ($\Upsilon_{\mathrm{disc}}=0.5M_{\odot}/L_{\odot}$ and $\Upsilon_{\mathrm{bulge}}=0.7M_{\odot}/L_{\odot}$). This was the approach taken in McGaugh et al. (2016). As indicated in Table 1, we mostly adopt the middle approach, that of a single free parameter (denoted as ‘1’ in the $\Upsilon$ column of Table 1). However, Model E takes the fixed values for the mass-to-light ratios (denoted as ‘0’), while Model F treats them as two free parameters (denoted as ‘2’). As mentioned previously, only 13 of the galaxies have a bulge component. For the remaining 72 galaxies, Models F and B are identical.

Finally, the fifth force contribution, $a_{5}$, is given by gradients of the $f(R)$ scalar field, $f_{R}$, via Eq. (6). In order to compute this, the scalar field $f_{R}$ needs to be calculated across the galactic disc by solving Equation (2). For this purpose, we use a spherical 1D scalar field solver, which takes as an input the mass distribution and the parameter $\bar{f}_{\mathrm{R0}}$. The details of this solver are discussed in the next subsection.

As indicated in Table 1, we adopt two different, but complementary approaches for the parameter $\bar{f}_{\mathrm{R0}}$:

• •

$\bar{f}_{\mathrm{R0}}$ is a free parameter in the fit of the rotation curve for each galaxy. The resulting spread of marginal $\bar{f}_{\mathrm{R0}}$ posteriors values across the galaxies of the sample can then be inspected for any significant clustering around a single preferred global value. This option is denoted as ‘varying’ under the $\bar{f}_{\mathrm{R0}}$ column in Table 1, and is used for all $f(R)$ models except H0-19 and I0-19.

• •

By contrast, Models H0-19 and I0-19 take a log-space grid of 20 values of $|\bar{f}_{\mathrm{R0}}|$, ranging from $\log_{10}\left(1.563\times 10^{-8}\right)$ to $\log_{10}\left(2\times 10^{-6}\right)$. These values are given in Table 2. Each of these 20 values are in turn globally imposed over the whole sample, for example Model H17 imposes a global $\log_{10}|\bar{f}_{\mathrm{R0}}|$ value of -5.921. These 20 models can then be compared with each other and with $\Lambda$CDM, providing constraints on $\bar{f}_{\mathrm{R0}}$.

4.2 Scalar Field Solver

This subsection describes the calculation of the scalar field $f_{R}$ for a given mass distribution and cosmic background value $\bar{f}_{\mathrm{R0}}$. The technique is essentially a 1D version of the Newton-Gauss-Seidel scalar field solver implemented in the $f(R)$ N-body code MG-GADGET (Puchwein et al., 2013) which was used in Paper I.

Given the finding of Paper I that the scalar field profiles within galaxies typically adopt a discoid shape, the spherically-symmetric 1D approximation is not ideal but is made for reasons of computational cost; the jump to 2D corresponds to roughly two orders of magnitude in computational time, typically making it prohibitively expensive to perform a sufficient number of iterations to achieve MCMC convergence.

On the other hand, an analytic calculation of the screening radius and fifth force for a spherical body, such as that given in Sakstein (2013), would make the computation instantaneous. However, this prescription was found to be too inaccurate and sensitive to the (somewhat arbitrary) choice of outer limit of integration. Thus, the 1D numerical computation described in this subsection represents a reasonably accurate and reasonably fast compromise. The robustness of this 1D approximation is examined in detail in Appendix §A, where a comparison to an axisymmetric 2D solver is presented.

To avoid unphysical positive values of the scalar field $f_{R}$ (e.g., due to finite step sizes in the iterative scheme), the calculation is performed in terms of the quantity $u\equiv\ln\frac{f_{R}}{\bar{f}_{R}\left(a\right)}$. Equation (2) can then be written as

[TABLE]

The role of the scalar field solver is to solve Equation (11) for $f_{R}$ across the galaxy, given an input density profile $\delta\rho$. This scalar field profile can then be used to calculate the fifth force contribution to the rotation curve via Eq. (6).

Equation (11) is discretised assuming spherical symmetry on a 1D radial grid. The outer edge of the grid needs to be at a radius larger than the Compton wavelength of the theory. For all of the $\bar{f}_{\mathrm{R0}}$ values within the bounded prior (see §4.6), 5 Mpc is a sufficiently large value. Finer resolution is required at smaller radii than at these large radii, so it is appropriate to use logarithmically spaced grid cells, i.e. the radial gridlines are equally spaced in the coordinate $x\equiv\ln r$, with constant grid spacing $h_{x}$. It was found that 175 cells between $r_{\mathrm{min}}=0.05$kpc and $r_{\mathrm{max}}=5$Mpc gave sufficiently accurate, converged results.

Defining

[TABLE]

where the index $i$ denotes the radial grid cells, Equation (11) is then

[TABLE]

This is solved, as in MG-GADGET, with an iterative Newton-Gauss-Seidel approach, where at iteration $n$ the scalar field is updated via

[TABLE]

In order to do this, we need discretised expressions for the Laplace operator on our grid, as well as the quantity $\frac{\partial\mathcal{L}_{i}^{n}}{\partial u_{i}^{n}}$. The Laplace operator in the coordinate $x\equiv\ln r$ is given by

[TABLE]

which is discretised as

[TABLE]

where $r_{i}$ indicates the radial position of the cell centre of cell $i$, while $r_{i-\frac{1}{2}}$ and $r_{i+\frac{1}{2}}$ are the positions of the inner and outer cell boundaries. Finally, the quantity $\frac{\partial\mathcal{L}_{i}}{\partial u_{i}}$ is given by

[TABLE]

Iterations of Eq. (14) are performed until $\Delta u_{i}\equiv u_{i}^{n+1}-u_{i}^{n}$ is less than $10^{-7}$ at all grid cells $i$. Through experimentation, this tolerance level has been found to give sufficiently accurate, converged results.

All that remains now is to provide a density profile $\rho_{i}$ for Eq. (12). For a given galaxy, the density profile depends on the parameter choices for the mass-to-light ratio(s) and the two dark matter halo parameters. Density data for the baryonic components are not provided with the SPARC data, so we have instead fitted $v_{\mathrm{gas}}$ and $v_{\mathrm{disc}}$ with exponential disc profiles and $v_{\mathrm{bulge}}$ with a Hernquist profile. In the cases of the 2D exponential discs, we spherically average the profiles for use in the 1D solver. These density profiles are then fed to the scalar field solver (multiplied by the mass-to-light ratio in the case of the stellar disc and bulge), along with the spherical dark matter halo profile, either NFW or DC14.

4.3 Environmental Screening

The environment of an object plays a role in its scalar field profile. As verified in Paper I, a galaxy embedded within a large-scale overdensity experiences an effective $\bar{f}_{\mathrm{R0}}$ that is lower than the cosmic background value, and as a result will have a larger screening radius than that of the case in which the galaxy occupied a region of cosmic mean density.

It has been shown in simulations (Zhao et al., 2011a, b; Cabré et al., 2012) that as a first approximation, the degree of environmental screening of a given galaxy can be quantified by the gravitational potential due to external sources $\Phi_{\mathrm{ext}}$.

For this very purpose of measuring the impact of environmental screening on tests of gravity, two groups—Cabré et al. (2012) and Desmond et al. (2018c)—have constructed ‘screening maps’: 3D maps of $\Phi_{\mathrm{ext}}$ throughout the local universe. In this work, we use the screening map of Desmond et al. (2018c) which uses updated techniques and catalogues compared to that of Cabré et al. (2012). Full details regarding the construction of the map can be found in the original paper, but it is worth noting that the value of $\Phi_{\mathrm{ext}}$ given by the screening map at a given point in space depends on the adopted value for $\bar{f}_{\mathrm{R0}}$, because the map sums contributions from all mass within the Compton wavelength $\lambda_{C}$ of the considered point, which relates to $\bar{f}_{\mathrm{R0}}$ via the approximate equation (Cabré et al., 2012)

[TABLE]

In this work, we use the screening map of Desmond et al. (2018c) to find $\Phi_{\mathrm{ext}}$ for each of the 147 galaxies that remain in the sample after making the first four cuts described in §3. The screening map uses a Bayesian methodology to reconstruct full posteriors on $\Phi_{\mathrm{ext}}$ at any point in space within $\sim 200$Mpc, by propagating uncertainties on the various inputs. For each of these galaxies, the posterior distribution for $\Phi_{\mathrm{ext}}$ is calculated for a log-spaced sequence of 20 values of $\bar{f}_{\mathrm{R0}}$, the same values used in Models H0-19 (see Table 2). This information is then used in two ways.

Firstly, we exclude from the sample any galaxy for which the $1\sigma$ upper bound on $|\Phi_{\mathrm{ext}}|/c^{2}$ is greater than $\frac{3}{2}|\bar{f}_{\mathrm{R0}}|$ (which is an approximate criterion for environmental screening, see Cabré et al. 2012), for any of the 20 considered values of $\bar{f}_{\mathrm{R0}}$. This is tantamount to removing from the sample any galaxies for which we believe environmental screening would be a significant effect in an $f(R)$ Universe. As indicated in §3, this cuts 62 of the 147 remaining galaxies once the other sample cuts have already been performed, leaving our final sample of 85 galaxies.

Secondly, we investigate the effect of the environmental contribution to the scalar field profile of the remaining galaxies in our sample. We consider a spherical top hat model, with a radius equal to the Compton wavelength $\lambda_{C}$, and a gravitational potential equal to the maximum posterior value of $\Phi_{\mathrm{ext}}$. As indicated by the quantity $\rho_{\mathrm{env}}$ in Figure 1, this structure is added to the overall density profile before it is fed to the scalar field solver. As shown in Table 1, this change is only made in the case of Model C, which is otherwise identical to Model B, i.e. $f(R)$+varying $\bar{f}_{\mathrm{R0}}$+NFW. The results of Model C can be compared to those of Model B to check the robustness of our results against the residual effects of environmental screening in our final sample. Note that because $\Phi_{\mathrm{ext}}$ is a function of $\bar{f}_{\mathrm{R0}}$, and $\bar{f}_{\mathrm{R0}}$ is a free parameter in Model C, the mass and size of this large-scale structure needs to be recalculated with each Monte Carlo realisation.

4.4 Stellar Self-screening

An assumption underlying our analysis is that the scalar field $f_{R}$ is sourced by all of the mass components: the dark matter, gas, and stars. However, in an $f(R)$ universe, one might expect a significant population of stars to be ‘self-screened’, and thus neither source nor couple to the fifth force (Davis et al., 2012).

This effect was the underlying principle for the aforementioned works by Vikram et al. (2018) and Desmond et al. (2019, 2018a, 2018b), which all searched for signals predicted to result from this differential coupling of the fifth force to stars and gas.

For the purposes of the present work, the expected effect of the stellar self-screening is that for a given $\bar{f}_{\mathrm{R0}}$, the screening radius will be smaller than if the effect is ignored. This was demonstrated in Paper I.

The contribution of a given star depends on its mass, radius, and environment. Thus, performing a fully self-consistent calculation incorporating the detailed contribution of the entire stellar population is difficult. We instead repeat the approach adopted in Paper I: we implement the extreme scenario in which all stars are assumed to be fully self-screened, and not sourcing the scalar field at all. In practice, this simply amounts to omitting the stellar input to the scalar field solver.

As in the test for environmental screening in the previous subsection, we employ this technique only for one model, Model D, which is otherwise identical to Model B, i.e. $f(R)$+varying $\bar{f}_{\mathrm{R0}}$+NFW. The ‘true’ result should then be bookended by these two extremes, so the results of Model D can be compared to those of Model B in order to gain an undestanding of the error induced in our inference by the effects of stellar self-screening.

4.5 Markov Chain Monte Carlo (MCMC)

We use an affine-invariant, parallel-tempered Markov chain Monte Carlo (MCMC) technique, using the publicly available python package emcee (Foreman-Mackey et al., 2013) to explore the posteriors of the free parameters of each model.

For a given galaxy with rotation curve data $\mathcal{D}$, the posterior probability $\mathcal{P}(\theta|\mathcal{D})$ for a given set of parameters $\theta$ is given by

[TABLE]

where $\mathcal{L}(D|\theta)$ represents the likelihood function, and $\pi(\theta)$ represents the prior probability distributions for the parameters $\theta$. The choices of priors will be discussed further in §4.6.

The pipeline for calculating this posterior probablity, for a given galaxy and a set of parameters $\theta$, is illustrated heuristically in Figure 1. The diagram specifically illustrates the pipeline for an $f(R)$ model with freely varying $\bar{f}_{\mathrm{R0}}$ and an environmental contribution (i.e. Model C); some of the details of the diagram would change for other models, and examples of such changes are given in the accompanying caption.

Assuming that the errors $\sigma$ on the data are Gaussian, the log-likelihood function $\ln\mathcal{L}(\mathcal{D}|\theta)$ for a given SPARC galaxy is given by

[TABLE]

where the sum is over the data points of an individual rotation curve.

In addition to the observational errors provided with the data points, for each galaxy we add in quadrature an additional constant error component $\sigma_{g}$, i.e.,

[TABLE]

$\sigma_{g}$ is then left as a free parameter in the fit for each galaxy. This additional error term is included in order to account for galaxy features not included in the model, e.g. spiral arms and other baryonic features, but a fit is penalised for adopting too large a value of $\sigma_{g}$, via the second term in the log-likelihood Eq. (20).

For each MCMC run, we use 30 walkers, 4 temperatures, and 5000 iterations (after burn-in). The chains have all been checked for convergence, both visually and by ensuring that the Gelman-Rubin statistic (Gelman & Rubin, 1992) $\mathcal{R}$ is sufficiently close to 1. In particular, for the overwhelming majority of the fits, $|\mathcal{R}-1|<0.01$, apart from one or two galaxies in some models, which exhibit some bimodality.

An example of a fit generated with this method is shown in Figure 2. The figure shows the observed rotation curve of NGC3741, along with the best fit (i.e. maximum a posteriori) rotation curve under the $f(R)$ model with freely varying $\bar{f}_{\mathrm{R0}}$ (i.e. Model B). Also shown are the various components: $v_{\mathrm{disc}}$, $v_{\mathrm{gas}}$, $v_{\mathrm{DM}}$, which combine in quadrature via Eq. (7) to give the model.

A fourth component is shown in Figure 2: $v_{5}$. This is the contribution of the fifth force to the rotation curve model, given by

[TABLE]

The only contribution to the model that is not a function of fit parameters is the gas curve $v_{\mathrm{gas}}$, which is provided with the data as a fixed quantity. $v_{\mathrm{disc}}$ is also provided with the data, but its overall normalisation is set by the mass-to-light ratio $\Upsilon$. It can be seen in the $v_{5}$ curve that the fit places a screening radius at around 2 kpc. This results in a mild upturn in the model rotation curve, which visibly appears to give a good fit to the data. It should be noted that this rotation curve was chosen because of this ideal behaviour, and is not necessarily representative of the sample at large.

4.6 Priors

As indicated in Table 1, different models have different free parameters. For the models with a freely varying $\bar{f}_{\mathrm{R0}}$, we adopt a flat prior in $\log_{10}|\bar{f}_{\mathrm{R0}}|$, between $|\bar{f}_{\mathrm{R0}}|=10^{-9}$ and $|\bar{f}_{\mathrm{R0}}|=2\times 10^{-6}$. For even larger values of $|\bar{f}_{\mathrm{R0}}|$ the model would struggle to screen the Milky Way at the solar radius, while even smaller values would be of limited interest due to the minuteness of the deviations from $\Lambda$CDM.

For any given galaxy, we adopt a flat prior for $\sigma_{g}$, between 0 and twice the maximum observed error for that galaxy.

The priors for the remaining parameters—the two dark matter parameters, $c_{\mathrm{vir}}$ and $v_{\mathrm{vir}}$, and where applicable, the one or two mass-to-light ratios $\Upsilon$—are less agnostic. Using flat priors was found to have the result that the fits are able to artificially inflate or deflate their halo masses so as to have a screening radius within the radial range of the rotation curve. As a consequence, unphysical vertical clustering features would appear in the stellar mass-halo mass diagrams and concentration-halo mass relation diagrams of the best fit models.

To avoid this undesirable behaviour, we follow the approach of Katz et al. (2017) and use empirical stellar mass-halo mass and concentration-halo mass relations from Moster et al. (2013) and Dutton & Macciò (2014) respectively as lognormal priors. These relations are depicted in Figure 3. For a fixed mass-to-light ratio, these relations translate directly to priors on $v_{\mathrm{vir}}$ and $c_{\mathrm{vir}}$ respectively. In the case of a freely varying $\Upsilon$, there is an increased freedom, as the stellar mass depends on the value of $\Upsilon$.

It should be noted that these relations are derived from $\Lambda$CDM simulations, so their applicability in an $f(R)$ universe is not entirely clear. $f(R)$ gravity simulations, such as those of Mitchell et al. (2019) suggest that for $|\bar{f}_{\mathrm{R0}}|=10^{-6}$ (denoted F6) and below, halo concentrations are enhanced by 0.1 dex or so. Taking this into account, we increase the width of our adopted concentration-halo mass relation by 0.1 dex; approximately equivalent to a doubling of the width. The priors should then encompass the ‘true’ relations of both a $\Lambda$CDM and an $f(R)$ universe.

For the stellar mass-halo mass relation, galaxy formation simulations are not yet available in the $f(R)$ paradigm, so the $f(R)$ relation is unknown. However, the work of e.g. Cataneo et al. (2016) suggests that the maximum increase to the halo mass function at the F6 level is small; approximately 10% or so. Assuming baryonic feedback mechanisms to have largely similar effects in an $f(R)$ universe, one might then expect the stellar mass-halo mass relation to be changed by a similar degree. Thus, a doubling of the width of the relation should also comfortably encompass the ‘true’ relation in the $f(R)$ paradigm. We approximately achieve this by adding 0.2 dex to the width.

As a final note, we follow Katz et al. (2017) in imposing some additional constraints. Firstly, we require that the baryon mass fraction is always less than the cosmological value. Furthermore, the parameters $\log_{10}c_{\mathrm{vir}}$, $\log_{10}\frac{v_{\mathrm{vir}}}{1\mathrm{m/s}}$, and $\log_{10}\frac{\Upsilon}{1M_{\odot}/L_{\odot}}$ are constrained to lie within the bounds (0, 2), (4, 5.7), and (-0.52, -0.1) respectively. The constraints on $c_{\mathrm{vir}}$ and $v_{\mathrm{vir}}$ are rather loose, while the constraint on $\Upsilon$ is somewhat tighter, but informed by the stellar synthesis models of Meidt et al. (2014); McGaugh & Schombert (2014); Schombert & McGaugh (2014).

5 Results

5.1 Constraints on $\bar{f}_{\mathrm{R0}}$

A widely used method for model comparison is the Bayesian Information Criterion (Schwarz, 1978), given by

[TABLE]

where $n$ is the number of data points, $k$ is the number of parameters of the model, and $\mathcal{L}$ is the maximised likelihood of the model. The first term acts on behalf of Occam’s razor, penalising overcomplicated models with too many parameters. Two models, 1 and 2, can then be compared by calculating the difference in the BIC: $\Delta\mathrm{BIC}\equiv\mathrm{BIC}_{1}-\mathrm{BIC}_{2}$. A positive (negative) value for this quantity indicates a preference for model 2 (1). We define values in the range $|\Delta\mathrm{BIC}|>2$ as ‘mildly significant’, and values in the range $|\Delta\mathrm{BIC}|>6$ as ‘strongly significant’ (Kass & Raftery, 1995).

Figure 4 shows a first comparison of an $f(R)$ model with freely varying $\bar{f}_{\mathrm{R0}}$ and a $\Lambda$CDM model, in both cases using NFW haloes (i.e. Models B and A respectively). The cumulative distribution function of $\Delta\mathrm{BIC}\equiv\mathrm{BIC}_{B}-\mathrm{BIC}_{A}$ is displayed. Model B treats $\bar{f}_{\mathrm{R0}}$ as a free parameter and so has one more free parameter than Model A, for which it is penalised by the first term in Eq. (23). Nonetheless, the majority of galaxies show some preference for $f(R)$, with 64/85 galaxies having $\Delta\mathrm{BIC}<0$. Perhaps more striking are the numbers of galaxies that have mildly significant ($|\Delta\mathrm{BIC}|>2$) preferences either way, with 3 galaxies preferring $\Lambda$CDM and 49 galaxies preferring $f(R)$. No galaxies have a strongly significant ($|\Delta\mathrm{BIC}|>6$) preference for $\Lambda$CDM, while 27 galaxies do for $f(R)$.

It is worth noting that the galaxies with the strongest preferences for $f(R)$ also typically infer much lower values for $\sigma_{g}$ under Model B than under Model A; Model A tries to compensate for a poor fit by increasing the scatter.

These numbers, and the preference for $f(R)$ implied by them, is perhaps unsurprising. Given a freely varying $\bar{f}_{\mathrm{R0}}$, the fits are able to place a screening radius anywhere within the radial range of the rotation curve, and are thus able to choose a favourable position such that the fit is improved. It is only in a minority of cases that no such position can be found. For these galaxies, a low $|\bar{f}_{\mathrm{R0}}|$ is chosen such that the galaxy is fully screened. This results in $\mathcal{L}_{f(R)}/\mathcal{L}_{\Lambda\mathrm{CDM}}\rightarrow 1$ and thus $\Delta\mathrm{BIC}\rightarrow\ln(n)$. For the SPARC galaxies, this will typically fall in the range $2-5$, and so it is impossible for a galaxy to have a ‘strongly significant’ preference for $\Lambda$CDM over an $f(R)$ model with a freely varying $\bar{f}_{\mathrm{R0}}$. This just reflects the fact that $\Lambda\mathrm{CDM}$ is contained within the $f(R)$ model as a limiting case.

Of course, in an $f(R)$ Universe, $\bar{f}_{\mathrm{R0}}$ would not vary freely from galaxy to galaxy as we have allowed it to here. Instead, there would be a single global value for $\bar{f}_{\mathrm{R0}}$. The natural question then is whether the individually inferred values for $\bar{f}_{\mathrm{R0}}$ show any clustering around a single global value. The upper panel of Figure 5 addresses this question, showing the marginal posterior distributions for $\bar{f}_{\mathrm{R0}}$ across all the galaxies, fitting again with Model B (freely varying $\bar{f}_{\mathrm{R0}}$+NFW). The overlaid coloured points show the values of $\bar{f}_{\mathrm{R0}}$ for the best fit models, which typically coincide with the peaks of the marginal posteriors, except in cases where there is some multimodality or $\bar{f}_{\mathrm{R0}}$ is not well constrained.

The key result shown in this panel is that there does not appear to be any immediate sign of clustering around a single value of $\bar{f}_{\mathrm{R0}}$. For most galaxies, $\bar{f}_{\mathrm{R0}}$ values in a reasonably narrow range are preferred, but these inferred values are almost evenly spread across two orders of magnitude, which is inconsistent with a single global value.

As an additional note, the Hubble types of the galaxies are indicated in the upper panel of Figure 5. The classifications used are those that accompany the SPARC data, and match the classification scheme of de Vaucouleurs et al. (1991): Sa-Sd types are spiral galaxies with spiral arms of decreasing tightness, Sm and Im types are Magellanic spirals and irregulars, and BCD types are ‘blue compact dwarfs’. There appears to be a correlation between inferred $\bar{f}_{\mathrm{R0}}$ and Hubble type. This is to be expected, as dwarf galaxies would be sensitive to lower values of $|\bar{f}_{\mathrm{R0}}|$, for which larger galaxies would be entirely screened. Conversely, larger galaxies would be sensitive to higher values, for which dwarf galaxies would be fully unscreened.

This finding demonstrates an important caveat to the above result that there is no obvious global value for $\bar{f}_{\mathrm{R0}}$: different galaxies have different sensitivity ranges for $\bar{f}_{\mathrm{R0}}$, so there would naturally be a spread in fitted $\bar{f}_{\mathrm{R0}}$ values when comparing across the whole sample. Even in an $f(R)$ Universe, only the galaxies with the true $\bar{f}_{\mathrm{R0}}$ within their sensitivity ranges would fit this value, and the other galaxies would instead just choose the incorrect $\bar{f}_{\mathrm{R0}}$ values that most improved their fits.

It is important then to invert the question; rather than asking whether all galaxies are consistent with one $\bar{f}_{\mathrm{R0}}$ value, we can ask whether one $\bar{f}_{\mathrm{R0}}$ value is consistent with all galaxies, i.e. whether a specific $\bar{f}_{\mathrm{R0}}$ value can be imposed globally that gives rotation curve models consistent with the observations of all galaxies, even if it is not the most preferred $\bar{f}_{\mathrm{R0}}$ value for a significant subpopulation of galaxies, possibly lying outside of their sensitivity range.

To address this question, we consider a series of 20 $f(R)$ models with globally imposed $\bar{f}_{\mathrm{R0}}$, i.e. Models H0-19. The $\bar{f}_{\mathrm{R0}}$ values utilised are given in Table 2. For each of the 20 $\bar{f}_{\mathrm{R0}}$ values and each of the 85 galaxies, the lower panel of Figure 5 shows the likelihood ratio $\ln\mathcal{L}_{f(R)}/\mathcal{L}_{\Lambda\mathrm{CDM}}$, taking the best fit parameters under Models H$x$ and A. This panel conveys similar information to the top panel; visually, there does not appear to be any obvious global value for $\bar{f}_{\mathrm{R0}}$. However, there is additional information to be gleaned. In the higher regions, $\log_{10}|\bar{f}_{\mathrm{R0}}|\gtrsim-6.5$, there are several galaxies which prefer such values of $\bar{f}_{\mathrm{R0}}$. However, there are also a great many galaxies for which these values provide a worsening of the fit compared with $\Lambda$CDM. On the other hand, in the region $\log_{10}|\bar{f}_{\mathrm{R0}}|\sim-7$, there are several galaxies which prefer such values, and the galaxies which prefer higher values are indifferent towards it, as a result of being completely screened in this region of $\bar{f}_{\mathrm{R0}}$. This region, therefore, appears to be the most promising within our sensitivity range.

A side-note is that the lower panel of Figure 5 also gives an understanding of the sensitivity of the sample. Beneath the diagonal green region is a white region, where galaxies are screened and $\Lambda$CDM and $f(R)$ are indistinguishable. As we go to progressively lower values of $|\bar{f}_{\mathrm{R0}}|$, more galaxies inhabit the white region. We define $\log_{10}|\bar{f}_{\mathrm{R0}}|\sim-7.2$ or $|\bar{f}_{\mathrm{R0}}|\sim 6\times 10^{-8}$ as the ‘sensitivity level’ of our sample, being the point at which roughly the half of the sample has become insensitive.

Rather than considering each galaxy individually, we can perform a model comparison across the sample as a whole. For each of the Models H0-19, we can perform an overall model comparison with $\Lambda$CDM by calculating the total likelihood ratio $\ln\mathcal{L}_{f(R)}/\mathcal{L}_{\Lambda\mathrm{CDM}}$. The log-likelihoods $\ln\mathcal{L}$ are calculated with Eq. (20) taking the best fit (i.e. maximum a posteriori) parameters for each galaxy, then summing the contributions of all galaxies in our sample. A more fully Bayesian approach would be to marginalise over the full parameter space rather than taking the best fit parameters. The likelihoods then become marginal likelihoods, and their ratios are known as Bayes factors. The calculation of these Bayes factors, however, is notoriously expensive so we instead consider the more easily calculable likelihood ratios, which are nonetheless robust tools for model comparison. The green curve in Figure 6 shows the result of this likelihood ratio calculation for the 20 $\bar{f}_{\mathrm{R0}}$ values listed in Table 2.

Figure 6 shows that $f(R)$ is disfavoured for values $\log_{10}|\bar{f}_{\mathrm{R0}}|\gtrsim-6.1$, reaching a decrease in log-likelihood of $\sim 30$ at its most extreme. Note that at higher values still, the likelihood for $f(R)$ starts to approach that of $\Lambda$CDM again. However, this corresponds to the extreme and unphysical case in which every galaxy is fully unscreened, and so the rotation curve models exactly resemble those of $\Lambda$CDM, but with an overall mass reduced to $3/4$ of the inferred $\Lambda$CDM mass, to account for the gravitational accelerations being $4/3$ times that of standard gravity.

At the opposite end of the spectrum, $\log_{10}|\bar{f}_{\mathrm{R0}}|\lesssim-7.6$, the likelihood ratios again approach identity. This corresponds to the regime in which every galaxy is fully screened, so the fits are now exactly identical to those of $\Lambda$CDM.

At intermediate values of $\bar{f}_{\mathrm{R0}}$, $f(R)$ appears to be favourable compared to $\Lambda$CDM. The peak of the signal is for the model with $\log_{10}|\bar{f}_{\mathrm{R0}}|=-7.03$, where there is a rather significant increase in the log-likelihood of around 66.

The lower panels of Figure 6 show the rotation curves and fits of the three galaxies with the highest individual likelihood ratios for Model H7 ($\log_{10}|\bar{f}_{\mathrm{R0}}|=-7.03$) compared with Model A. In other words, the three galaxies that contribute the most to the peak in the upper panel of Figure 6. These galaxies are DDO161 ($\ln\mathcal{L}_{f(R)}/\mathcal{L}_{\Lambda\mathrm{CDM}}=14.96$), UGC00891 (6.14), and NGC3109 (6.13).

It is striking that the three galaxies are all rather similar; dwarf galaxies (Hubble types all Im or Sm) with similar radial extents and maximum rotation speeds. Furthermore, in each case, the behaviour of the fits is similar. The observed rotation curves are poorly fit with cuspy NFW profiles, which place too much mass in the inner regions in order to fit the velocity data in the outer regions. The reason the $f(R)$ fits provide a significant improvement is that the presence of the fifth force in the outer regions allows a reduced mass in the inner regions.

The DC14 profile for dark matter haloes Di Cintio et al. (2014) is empirically derived from $\Lambda$CDM simulations incorporating stellar feedback. The inner slope of the profile depends on the stellar content of the galaxy, giving a more cored profile to galaxies with intermediate ($\sim 0.5\,\%$) stellar mass fractions. Katz et al. (2017) found that fitting the SPARC galaxies with a DC14 profile gives a consistent improvement compared to NFW, particularly in cases such as the galaxies depicted in the lower panels of Figure 6, i.e. dwarf galaxies requiring a cored profile.

The upper panel of Figure 6 also shows a comparison of Models G and A (the dashed purple line), i.e. a model comparison of $\Lambda$CDM models with DC14 and NFW haloes. The increase in log-likelihood is nearly 300; a very significant improvement, and far more significant than the increase in log-likelihood for any given $f(R)$ model with an NFW halo. The reason for this improvement is readily apparent in the rotation curve fits in the lower panels of Figure 6. Echoing the findings of Katz et al. (2017), the DC14 fits to the three galaxies give consistently better agreement with the data than NFW. We find that this holds for both $\Lambda$CDM and $f(R)$ gravity.

The final model comparison shown in the upper panel of Figure 6 is between Models I0-19 and Model A, i.e. a comparison of DC14+$f(R)$ and NFW+$\Lambda$CDM. Model I does not give a marked improvement to the fits compared to Model G. In other words, when cored halo profiles are assumed, the $f(R)$ signal at $\log_{10}|\bar{f}_{\mathrm{R0}}|\sim-7$ largely disappears, so that models with a cored halo profile and a fifth force are not significantly better than those with a cored halo profile and no fifth force.

All of this is not to say, however, that no galaxy prefers an $f(R)$ fifth force to a cored halo. We repeat the analysis shown in Figure 4, i.e. the distribution of $\Delta\mathrm{BIC}$ for a freely varying $\bar{f}_{\mathrm{R0}}$ model versus $\Lambda$CDM, but replacing the $\Lambda$CDM+NFW model with the $\Lambda$CDM+DC14 model (i.e. $\Delta\mathrm{BIC}\equiv\mathrm{BIC}_{B}-\mathrm{BIC}_{G}$). In this case, we find that it is now only a minority of galaxies that show some preference for $f(R)$, with $29/85$ galaxies having $\Delta\mathrm{BIC}<0$, down from 64 in the NFW case. The numbers at the extremes have also shifted, with 14 galaxies showing a strongly significant ($|\Delta\mathrm{BIC}|>6$) preference for $f(R)$ and 21 for $\Lambda$CDM, compared with 27 and 0 respectively in the NFW case. These numbers reinforce the idea conveyed by Figure 6: $\Lambda$CDM with a cored halo profile appears to be preferable over a cuspy halo and $f(R)$ for most of the galaxies in the sample.

Figure 7 shows rotation curve fits of the 3 galaxies with the most extreme preferences for $f(R)$ over $\Lambda$CDM+DC14 (i.e. the three galaxies with the most negative values of $\Delta\mathrm{BIC}\equiv\mathrm{BIC}_{B}-\mathrm{BIC}_{G}$). These are UGC11820 ($\Delta\mathrm{BIC}=-22.1$), NGC2403 (-31.2), and UGC05253 (-47.1). In the case of UGC11820, it appears that the $f(R)$ model does indeed do a markedly better job of fitting the rotation curve. There is not a visible ‘upturn’ feature in the rotation curve, but the the overall shape of the rotation curve is much better described by an NFW halo enchanced with a fifth force in the outer regions than by a cored dark matter halo. Indeed, considering the goodness of fit in each case, the reduced chi-squared statistic $\chi^{2}_{\nu}$, drops from 13.4 to 0.8.

However, in the case of UGC05253, it would appear that there are large-scale features in the rotation curve, both in the central (magnified) region and in the outer stretches, that are not adequately captured by our baryonic models. This picture is borne out by the values of $\chi^{2}_{\nu}$: 17.7 for Model G, and 9.5 for Model B. The screening radius is being imposed to markedly improve the quality of the fits, but in actuality neither model does particularly well. NGC2403 shows similar behaviour, although perhaps not quite so starkly. A significant upturn feature at $\sim 4$ kpc is clearly visible. Both models, with and without a screening radius, approximate this feature. However, the baryonic models are unable to fully fit its magnitude, so a screening radius is imposed to bridge the gap. Again, neither model performs particularly well: $\chi^{2}_{\nu}$ drops from 5.5 to 2.6.

In addition to these caveats, it should be once again borne in mind that the galaxies preferring Model B to Model G are not necessarily showing any preference for a single global $\bar{f}_{\mathrm{R0}}$ value. In fact, the inferred values of $\log_{10}|\bar{f}_{\mathrm{R0}}|$ for the extreme galaxies shown in Figure 7 are rather discrepant: -7.07 for UGC11820, -6.57 for NGC2403, and -6.19 for UGC05253. As shown in Figure 6, no single, global $\bar{f}_{\mathrm{R0}}$ model with a NFW profile is in better agreement with the data than a $\Lambda$CDM model with a cored DC14 profile.

5.2 Environmental Screening

As described in §4.3, we investigate a model with freely varying $\bar{f}_{\mathrm{R0}}$, NFW haloes, and an additional large-scale overdensity to account for the environmental contribution to the scalar field profile. This is Model C, which is identical to Model B in every respect other than the presence of this large-scale overdensity. For a given $\bar{f}_{\mathrm{R0}}$ as drawn in the Monte Carlo sampling, the structure will have gravitational potential equal to $\Phi_{\mathrm{ext}}$ calculated from the screening maps of Desmond et al. (2018c), and size equal to the Compton wavelength calculated using Eq. (18).

Figure 8 shows a comparison of the $\bar{f}_{\mathrm{R0}}$ values inferred from the maximum of the posterior of models Models B and C. Results for the 85 galaxies in our sample are shown, as well as results for the 62 galaxies that were excluded from the sample because of their significant environmental contribution.

Upon including the environmental contribution, the majority of the galaxies in Figure 8 move towards higher $|\bar{f}_{\mathrm{R0}}|$ and to higher $\Phi_{\mathrm{ext}}$. This makes physical sense: for a given $f_{R0}$, the addition of an overdense environment will move the screening radius of a galaxy outward. For a rotation curve fit that is improved by a screening radius at a specific location, this effect is compensated by a higher value of $|\bar{f}_{\mathrm{R0}}|$, hence the shift towards higher $|\bar{f}_{\mathrm{R0}}|$. The shift towards higher $\Phi_{\mathrm{ext}}$ is due to the implicit relation between $\bar{f}_{\mathrm{R0}}$ and $\Phi_{\mathrm{ext}}$: $\Phi_{\mathrm{ext}}$ is integrated up to the Compton wavelength, which in turn is determined by Eq. (18). A higher $|\bar{f}_{\mathrm{R0}}|$ will therefore lead to a higher $\Phi_{\mathrm{ext}}$ (note that this is also the reason for the correlation between inferred $\bar{f}_{\mathrm{R0}}$ and $\Phi_{\mathrm{ext}}$ evident in Figure 8). For these cases, this change is visually much larger for the excluded (environmentally screened) galaxies than for the included galaxies.

There is however, also a subset of galaxies for which the jumps are rather large, and these appear to be equally likely to be towards higher or lower $|\bar{f}_{\mathrm{R0}}|$. In these cases, $f_{R0}$ is typically very poorly constrained and the marginal posteriors are rather broad, flat distributions (both with and without the environmental contribution). In such scenarios, the peak of the posterior can undergo quite large shifts (in either direction) between models, but the distributions are nonetheless mutually consistent.

The finding that the excluded galaxies exhibit a large change in inferred $\bar{f}_{\mathrm{R0}}$ (median $|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||=0.37$) gives a post hoc justification for their exclusion from the sample. Meanwhile, the finding that the included galaxies typically do not change their inferred $\bar{f}_{\mathrm{R0}}$ values significantly (median $|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||=0.04$) gives a post hoc justification for our neglecting the environmental screening effect in most of our models (see Table 1).

5.3 Stellar Self-Screening

Figure 9 shows a comparison of the $\bar{f}_{\mathrm{R0}}$ values from the best fit models for all galaxies under Models D and B, i.e. $f(R)$ models with and without stellar self-screening, as a function of total stellar luminosity. In Model D, stars do not act as a source of the fifth force, as would be expected if they self-screen completely. Self-screening would also prevent a fifth force from accelerating stars, this is however not relevant here as the measured rotation curves are based on observations of the gas component.

The mean (absolute) change in $|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||$ between Models D and B is 0.11, but this number is dominated by a few outliers. In general, the change is not too large (median $|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||=0.02$), with a larger change taking place for more luminous galaxies.

The galaxies with the two largest changes in inferred $\bar{f}_{\mathrm{R0}}$ are UGC05918 ($|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||=1.1$) and UGC07866 (1.7). In the case of UGC05918, under both Models B and D, the best fit models are those in which the galaxy is fully screened. In such a scenario, $\bar{f}_{\mathrm{R0}}$ is poorly constrained as any $|\bar{f}_{\mathrm{R0}}|$ value below a certain threshold will fully screen the galaxy, and all such values have equal posterior probability.

Meanwhile, in the case of UGC07866, the posterior in both Models B and D is bimodal. In one mode, the galaxy is fully screened, while in the other mode, the galaxy is fully unscreened, and the mass of the halo and stellar components are renormalised by a factor of 3/4 to give a very similar rotation curve model. Assuming stellar self-screening causes a slight shift in the relative weights of the two modes, which caused the best fit model to spontaneously ‘jump’ from one mode to the other.

Neither of these cases, nor indeed the other outliers with large changes in inferred $\bar{f}_{\mathrm{R0}}$, should be causes for concern. In each of these cases, the best fit model either fully screens or fully unscreens the galaxy. So, while the best fit values are different under the two models, the marginal posterior distributions for $\bar{f}_{\mathrm{R0}}$ in both cases are wide, flat distributions that are consistent with each other. The majority of galaxies, however, exhibit the more typical behaviour, where a screening radius in a particular location is favourable, which corresponds to a specific value of $\bar{f}_{\mathrm{R0}}$. For these galaxies, the required value of $\bar{f}_{\mathrm{R0}}$ is very similar for Models B and D, which is reassuring and justifies neglecting the effect of stellar self-screening in most of the models listed in Table 1.

5.4 Mass-to-light Ratios

As described in §4.1, Model B (as well as all other models considered besides E and F), performs fits with a single freely varying $\Upsilon$, which applies to both the bulge and disc components. Model E instead fixes $\Upsilon$ to be equal to $0.5M_{\odot}/L_{\odot}$ for the disc and 0.7 for the bulge, while Model F fits two freely varying ratios: $\Upsilon_{\mathrm{disc}}$ and $\Upsilon_{\mathrm{bulge}}$.

It should be noted that of the 85 galaxies in the sample, only 13 have a detected bulge component. For the remaining 72 galaxies, Models B and F are identical, and should therefore give identical results.

Figure 10 shows a comparison of results from Models B, E, and F. For all 85 galaxies in the sample, the figure shows the $\bar{f}_{\mathrm{R0}}$ value for the best fit model, for the B, E and F models. The mean (absolute) change is $|\Delta\log_{10}|\bar{f}_{\mathrm{R0}}||=0.14$ between Models B and E, but this is dominated by a few outliers, and the median change is 0.03. For Models B and F, considering only the 13 galaxies with bulge components, the mean change is 0.07 and the median is 0.01.

These results suggest that fixing the values of the mass-to-light ratios might be too crude an approximation. Indeed, most of the visibly large changes in Figure 10 are due to Model E. On the other hand, the reassuringly small changes between Models B and F imply that taking the stellar disc and bulge to have the same mass-to-light ratio is a sufficient approximation.

6 Discussion and Conclusions

In this work, we have searched for potential signatures of Chameleon-$f(R)$ modified gravity in the high-quality HI/H$\alpha$ rotation curve measurements of the SPARC sample. After several cuts to the sample, including an exclusion of all galaxies likely to be environmentally screened, 85 of the original 175 galaxies remained.

Each galaxy is modelled as a stellar disc, gaseous disc, dark matter halo, and (where appropriate) a stellar bulge. Then, for a given $\bar{f}_{\mathrm{R0}}$ value, we solve the $f(R)$ equations assuming spherical symmetry, and calculate the implied fifth force contribution to the rotation curve. Assuming a Gaussian likelihood, we are then able to use an MCMC technique to explore suitable models and search for evidence of modified gravity.

As mentioned in the Introduction, all of the code, analysis tools, and plotting scripts used in this work have been made into a publicly available python 3 package222https://github.com/aneeshnaik/spam.

Our main findings are:

• •

In models with an NFW halo and an $\bar{f}_{\mathrm{R0}}$ that is allowed to vary freely for each galaxy, most of the galaxy sample (64/85) indicated some preference for $f(R)$, according to the Bayesian Information Criterion. Furthermore, 27 galaxies showed a strongly significant ($|\Delta\mathrm{BIC}|>6$) preference for $f(R)$.

• •

Looking at the marginal posterior distributions of $\bar{f}_{\mathrm{R0}}$ of all galaxies in our final sample, most galaxies have reasonably tight posteriors and $\bar{f}_{\mathrm{R0}}$ appears well constrained, but the spread of these inferred $\bar{f}_{\mathrm{R0}}$ values across the sample is very broad, spanning roughly two orders of magnitude. This is inconsistent with a single global value for $\bar{f}_{\mathrm{R0}}$.

• •

This finding is confirmed when analysing models with fixed, globally imposed $\bar{f}_{\mathrm{R0}}$ values. Models with $\log_{10}|\bar{f}_{\mathrm{R0}}|\gtrsim-6.1$ are highly unfavourable compared to $\Lambda$CDM. However, models with lower $|\bar{f}_{\mathrm{R0}}|$ appear to be preferred over $\Lambda$CDM. This signal reaches a peak at $\log_{10}|\bar{f}_{\mathrm{R0}}|\sim-7$, where the log-likelihood ratio $\ln\mathcal{L}_{f(R)}/\mathcal{L}_{\Lambda\mathrm{CDM}}\sim 70$.

• •

The galaxies that dominate this signal are dwarf galaxies, which have previously been noted to exhibit the core-cusp problem.

• •

Modelling with cored DC14 haloes without a fifth force, the overall log-likelihood shows an increase of nearly 300, far more significant than any $f(R)$-model with a single, global $\bar{f}_{\mathrm{R0}}$ value.

• •

Modelling with cored DC14 haloes with a fifth force does not provide a significant improvement over the case of DC14+$\Lambda$CDM, and the signal at $|\bar{f}_{\mathrm{R0}}|\sim 10^{-7}$ has largely vanished.

• •

There are nonetheless some individual galaxies which show a preference for NFW+$f(R)$ over DC14+$\Lambda$CDM. However, for some of these galaxies it appears to be the case that neither model performs particularly well, as a result of large-scale baryonic features (e.g., due to prominent spiral arms) in the rotation curves inadequately captured by our baryonic models.

• •

These results appear to be robust to the effects of environmental screening. Having removed a significant fraction of the galaxies for which environmental screening was believed to play a significant role, the results of the $\bar{f}_{\mathrm{R0}}$ inference for the remaining 85 galaxies were not affected significantly when a large-scale environment was added to the modelling pipeline.

• •

Similar analysis was performed to test for the effects of stellar self-screening and a more general treatment of the mass to light ratio. It was also found here that the results of our inference are robust to these effects.

We thus end by reporting an absence of convincing evidence of modified gravity down to the sensitivity level of our sample at $|\bar{f}_{\mathrm{R0}}|\sim 6\times 10^{-8}$, as the improvements to the rotation curve fits due to the $f(R)$ fifth force (which peak at around $|\bar{f}_{\mathrm{R0}}|\sim 10^{-7}$) are more readily explained by galaxies having cored dark matter profiles than by a genuine signal of modified gravity.

The only work, to our knowledge, that has previously addressed this degeneracy between cored dark matter profiles and screened modified gravity theories is that of Lombriser & Peñarrubia (2015), which considered the velocity dispersions of local dwarf spheroidal galaxies Fornax and Sculptor. The radial slopes of their velocity dispersions have previously been shown to be consistent with cored mass profiles (Walker & Peñarrubia, 2011), and Lombriser & Peñarrubia (2015) show that they can alternatively be interpreted as being due to a chameleon fifth force which is unscreened only in the outer parts of these dwarfs. However, the $\bar{f}_{\mathrm{R0}}=-10^{-7}$ models that we prefer may not be suitable for this as they would likely, depending on the exact Milky Way mass and the environmental effects of the Local Group, screen the Milky Way halo to beyond the positions of the Fornax and Sculptor dwarfs.

It would therefore be interesting to extend the remit of the present work by investigating a wider range of screened modified gravity theories to investigate the constraints that can be derived on them, and whether any gravity model can bring cuspy NFW profiles into good agreement with galaxy rotation curve data.

Acknowledgements

We would like to thank Federico Lelli, Harley Katz, Cameron Lemon, and the anonymous reviewer for helpful comments and discussions.

APN thanks the Science and Technology Facilities Council (STFC) for their PhD studentship. EP acknowledges support by the Kavli foundation. ACD acknowledges partial support from STFC under grants ST/L000385 and ST/L000636. DS acknowledges ERC starting grant 638707 and support from the STFC. HD is supported by St John’s College, Oxford and acknowledges financial support from ERC Grant No 693024 and the Beecroft Trust.

This work used the DiRAC (www.dirac.ac.uk) system: Data Analytic at the University of Cambridge [funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1]. DiRAC is part of the National E-Infrastructure.

Appendix A 1D Approximation

As discussed in §4.2, when solving Equation (2) to calculate the scalar field, we assume spherical symmetry to reduce the computational cost. In Paper I, it was found that the scalar field typically takes a discoid shape in galaxies, so assuming spherical symmetry could potentially lead to our calculations differing from the ‘true’ scalar field profiles, introducing error into our inferences for $\bar{f}_{\mathrm{R0}}$. In particular, when a galaxy is partially screened, a given $\bar{f}_{\mathrm{R0}}$ might correspond to an appreciably different location for the screening radius when the assumption of spherical symmetry is introduced.

To test the robustness of this 1D approximation, we compare some of our results with those that would have been obtained if the 1D solver was replaced with a 2D solver, in which the assumption of spherical symmetry has been replaced with the assumption of azimuthal symmetry. This 2D solver has been found to give galactic scalar field profiles in excellent agreement with those calculated by the full 3D solver in MG-GADGET.

In the 2D solver that we use for comparison, Equation (11) is discretised on a 2D grid in the spherical polar coordinates $r$ and $\theta$. As in the 1D case, we use logarithmically spaced grid cells for the coordinate $r$, i.e. the radial gridlines are equally spaced in $x\equiv\ln r$, with constant grid spacing $h_{x}$.

For the polar angle $\theta$, we adopt a linear grid, ranging from $\theta=0$ to $\theta=\pi$, with $n_{\theta}$ grid cells and constant grid spacing $h_{\theta}$. If $n_{\theta}$ is required to be an odd number, then the $\left(\frac{n_{\theta}+1}{2}\right)$th grid cell corresponds to the disc plane. Through experimentation, it was found that $n_{\theta}=101$ gave sufficiently accurate results.

Defining

[TABLE]

where $i$ and $j$ denote grid cells along the radial and polar directions respectively, we can rewrite Equation (11) as

[TABLE]

This is then solved, with the same Newton-Gauss-Seidel technique, Equation (14). In order to do this, we need discretised expressions for the Laplace operator on our grid, as well as the quantity $\frac{\partial\mathcal{L}_{ij}^{n}}{\partial u_{ij}^{n}}$. Writing $\sin\theta$ as $s$, the Laplace operator in our coordinates $x\equiv\ln r$ and $\theta$ is given by

[TABLE]

which is discretised as

[TABLE]

Finally, the quantity $\frac{\partial\mathcal{L}_{ij}}{\partial u_{ij}}$ is given by

[TABLE]

With a solution for $f_{R}$ everywhere, the fifth force is then given by Equation (6), evaluating the gradients of the scalar field in the disc plane, $\theta=\pi/2$.

We take 10 randomly chosen galaxies from our sample, and rerun our MCMC analysis with this 2D solver for Model B, i.e. the $f(R)$ model with freely varying $\bar{f}_{\mathrm{R0}}$, a single mass-to-light ratio, and NFW haloes. As in the 1D case, we use 30 walkers operating at 4 temperatures. With the 1D solver, the walkers each performed 5000 iterations after a burn-in period of 5000 iterations, but here we limit the computation to 500 iterations after a burn-in of 500, for reasons of computational cost. As a consequence, the convergence is more uncertain: whereas the Gelman-Rubin statistic $\mathcal{R}$ was previously confined to the region $|\mathcal{R}-1|<0.01$ for the vast majority of fits, the statistics for several of the galaxies are now in the range $0.01<|\mathcal{R}-1|<0.05$.

Figure 11 shows the 10 marginal posterior distributions for $\bar{f}_{\mathrm{R0}}$, calculated with 2D and 1D solvers. There is remarkably good agreement between the two solvers, demonstrating the validity of our 1D approximation. In two cases, NGC6015 and NGC2841, the posterior is bimodal for the 1D solver, but unimodal for 2D. However, as the missing mode in both cases falls within the fully unscreened regime, this is unlikely to be related to any difference between the scalar field solutions of the two solvers, but to the convergence issues described above, which have likely prevented the 2D MCMC from finding the second mode.