© ESO 2020

1. Introduction

One of the most outstanding open questions in astrophysics is the mass discrepancy problem: the apparent amount of mass in the Universe is roughly ten times larger than the mass that is visible through its electromagnetic emission (Ostriker & Peebles 1973; Ostriker et al. 1974; Sanders 2010). This discrepancy occurs all the way from the largest scale of cosmic microwave background radiation to the galactic scale whenever we model the observed dynamics of astrophysical systems or their gravitational lensing features, using general relativity or its Newtonian weak field limit (Rubin & Ford 1970; Sanders 1990; Paraficz et al. 2016; Planck Collaboration VI 2020). Such observations are usually reconciled with expectations from the theory of gravity by assuming the existence of non-baryonic dark matter, the specific properties of which are still under debate (van Albada et al. 1985; Vikhlinin et al. 2006; Kunz et al. 2016; Bode et al. 2001; Spergel & Steinhardt 2000; Roszkowski et al. 2018).

The most widely investigated interpretation is the cold dark matter (CDM) paradigm (e.g. Dodelson et al. 1996), where the cosmic structure forms as a result of the aggregation of smaller structures. In this context of substantially stochastic merging processes, some regularities in the observed properties of disk galaxies do not appear to occur naturally; they might, rather, require a substantial fine-tuning between the properties of the baryonic matter and the expected properties of the dark matter halo embedding the galaxy (McGaugh 2005; Famaey et al. 2018). For example, we might not naively expect a tight relation between the flat rotation velocity v_f and the baryonic mass, the so-called baryonic Tully–Fisher relation (Tully & Fisher 1977; McGaugh et al. 2000; Lelli et al. 2016a), because v_f is set by the depth of the gravitational potential of the dark matter halo that contains ∼90% of the total mass of the galaxy and should hardly be affected by the 10% baryonic mass: in the CDM framework, a very careful balance between star formation efficiency and stellar feedback might be necessary (McGaugh 2012; Lelli et al. 2016a). Similarly, the observed centripetal acceleration implied by the rotation curve, $g_{\text{obs}}=v_{\text{obs}}^2(R)/R$, is tightly correlated with the Newtonian acceleration due to the baryonic matter distribution, g_bar, and the two accelerations perfectly coincide only above a single acceleration scale which is common to all galaxies (McGaugh et al. 2016), whereas, at decreasing accelerations, the discrepancy between g_obs and g_bar monotonically increases. This radial acceleration relation (RAR), although recently disputed (Rodrigues et al. 2018), could be particularly relevant because both g_obs and g_bar and their uncertainties are completely independent of each other (Li et al. 2018).

These observed regularities, although some of them appear reproducible in the CDM model (e.g. Ludlow et al. 2017), might suggest that an alternative solution to the dark matter paradigm is a modification of the theory of gravity. MOdified Newtonian Dynamics (MOND, Milgrom 1983a) is capable of modeling, even sometime predicting (Sanders & McGaugh 2002), these observations by assuming a breakdown of the Newtonian gravity in low-acceleration environments, where the acceleration threshold is set to a₀ = 1.2 × 10⁻¹⁰ m s⁻² by observations (e.g. McGaugh 2004).

More recently, refracted gravity (RG, Matsakos & Diaferio 2016) has proposed a different modification of the theory of gravity that is based on a completely different idea from MOND but it is expected to share most of its successes. RG is a classic theory of gravity, whose modified Poisson equation includes the gravitational permittivity, ϵ(ρ), a monotonic function of the local mass density, ρ, that boosts the gravitational field in low-density environments. RG can be reformulated as a scalar-tensor theory (Sanna et al., in prep.) and would thus share most of their general properties (e.g. Quiros 2019; Kobayashi 2019).

Specifically, the scalar field, which is non-minimally coupled to the gravitational field, is responsible both for the gravitational permittivity and could, thus, remove the need for dark matter and for the accelerated expansion of the universe. This feature is particularly attractive because both dark matter and dark energy can be mimicked by a single scalar field, similarly to other models that attempt to unify dark matter and dark energy, for example, f(R) theories (e.g. Sotiriou & Faraoni 2010), quartessence theories (e.g. Brandenberger et al. 2019), mimetic gravity (Sebastiani et al. 2017), or generalised Chaplygin gas and beyond (Hernández-Almada et al. 2019). Since the scalar field in RG is a dynamical quantity, RG predicts a time evolution of the equation of state of the effective dark energy that can, in principle, be measured by upcoming space missions like Euclid¹.

Here, we test the viability of RG by modelling the observed dynamics of disk galaxies. To provide the most stringent tests of the full dynamics of a disk galaxy, rather than modelling rotation curves alone, we consider a sample of galaxies where both the rotation curves and the velocity dispersion profiles, in the direction perpendicular to the disk, are available.

The DiskMass Survey (DMS, Bershady et al. 2010a) provides a sample that is fitting for our purpose. It contains 46 galaxies from the Uppsala General Catalogue (UGC) whose disks appear close to face-on; for 30 galaxies the measures of both the rotation curves and the vertical velocity dispersion profiles are publicly available. The DMS collaboration modelled the galaxy dynamics with Newtonian gravity by adopting the disk-scale heights derived from the observations of edge-on galaxies; they obtained sub-maximal disks, with mass-to-light ratios that are systematically smaller than what is expected from stellar population synthesis (SPS) models (Bell & de Jong 2001).

The DMS sample was also used by Angus et al. (2015) to test MOND; they found mass-to-light ratios consistent with the SPS values but with disk-scale heights systematically smaller than those inferred from the observations of edge-on galaxies. Milgrom (2015), however, pointed out that this inconsistent disk thickness might originate from an observational bias: the measured velocity dispersion is inferred from the absorption lines near the V-band of the integrated spectra, which are dominated by the younger stellar population (Aniyan et al. 2016); this measured velocity dispersion is thus smaller than the velocity dispersion of the older stellar population which sets the estimate of the disk-scale height from near infrared photometry of edge-on galaxies (Kregel et al. 2002; Pohlen et al. 2000; Schwarzkopf & Dettmar 2000; Xilouris et al. 1997, 1999; Bershady et al. 2010b).

This bias would also explain the low mass-to-light ratios estimated by the DMS collaboration because the disk surface mass density is proportional to the ratio between the velocity dispersion and the disk-scale height and it is thus underestimated (Aniyan et al. 2016). The analysis of the dynamics of the DMS galaxies that we present here could also be affected by this bias. In Sect. 4.1 below, we estimate the amount of this bias and find that it is indeed consistent with the estimates of Milgrom (2015) and Aniyan et al. (2016).

The outline of this paper is as follows. Section 2 summarises the main features of the RG theory. In Sect. 3, we test RG by modelling the rotation curves of the DMS galaxies alone, whereas in Sect. 4, we model both their rotation curves and their vertical velocity dispersion profiles. We describe our model of the galaxy mass distribution and our Poisson solver in Appendixes A and B, respectively.

Modelling the dynamics of each galaxy in RG requires two parameters for the galaxy, namely, its mass-to-light ratio and its disk-scale height, along with three RG free parameters. In Sect. 5, we show that the DMS sample could, in principle, be modelled by a single set of these three free RG parameters. In Sect. 6, we show that RG can also model the RAR of the DMS sample, although some tensions indeed exist. We present our conclusions in Sect. 7. We adopt the Hubble constant H₀ = 73 km s⁻¹ Mpc⁻¹, as per Martinsson et al. (2013a), throughout this paper.

2. Refracted gravity

RG is a classical theory of gravity inspired by the behaviour of electric fields in matter (Matsakos & Diaferio 2016): when an electric field line crosses a dielectric medium with a non-uniform permittivity, it suffers a change both in direction, namely, a refraction, and in magnitude. To mimic this behaviour in a gravitational field, RG adopts the modified Poisson equation:

$$\begin{array}{c}\hfill \mathrm{\nabla }·[ϵ(\rho )\mathrm{\nabla }\varphi ]=4\pi G\rho ,\end{array}$$(1)

where ϕ is the gravitational potential, and ϵ(ρ) is the gravitational permittivity, an arbitrary monotonically increasing function of the mass density ρ. By adopting for ϵ(ρ) the asymptotic limits,

$$\begin{array}{c}\hfill ϵ(\rho )=\{\begin{array}{cc}1,\hfill & \rho \gg {\rho }_{\mathrm{c}}\hfill \\ {ϵ}_0,\hfill & \rho \ll {\rho }_{\mathrm{c}},\hfill \end{array}\end{array}$$(2)

the RG Poisson equation is reduced to the Newtonian form,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\varphi =4\pi G\rho ,\\ \hfill \end{array}$$(3)

in environments where the local density is much larger than the critical density ρ_c. On the contrary, for a constant ϵ₀ <  1, the gravitational field is boosted in low-density environments.

For a spherically symmetric system, with mass M(<r) within the radius r, the integration of Eq. (1) yields ∂ϕ/∂r = [G/ϵ(ρ)]M(<r)/r²; therefore, the gravitational field has the same direction and the same dependence on r as the Newtonian field, but it is enhanced by the factor 1/ϵ(ρ). For systems that are not spherically symmetric, expanding the divergence in Eq. (1),

$$\begin{array}{c}\hfill \frac{\partial ϵ}{\partial \rho }\mathrm{\nabla }\rho ·\mathrm{\nabla }\varphi +ϵ(\rho ){\mathrm{\nabla }}^2\varphi =4\pi G\rho ,\end{array}$$(4)

shows that ϕ depends both on the density field ρ, according to the second term in the left-hand side of the equation, as in the spherically symmetric case, and on its variation, according to the first term in the left-hand side of the equation.

We thus see that the analogy between the gravitational field and the electric field in a dielectric medium occurs for non-spherical systems. For example, in disk galaxies, the boost of the gravitational field can be visualised as a focussing of the gravitational field lines towards the disk plane, as they are refracted by the low-density regions above and below the disk. According to this feature, RG predicts that increasingly flatter systems should require increasing dark matter content when interpreted in Newtonian gravity, as suggested by preliminary studies of elliptical galaxies (Deur 2014).

As Matsakos & Diaferio (2016) show, this focussing yields, for the gravitational acceleration g in low-density regions at large distances R from the disk centre, the asymptotic behaviour g ∼ (|g_N|a₀)^1/2 ∝ R⁻¹, where g_N is the Newtonian acceleration and a₀ coincides with the MOND critical acceleration that is set by the observed normalisation of the Tully–Fisher relation. This asymptotic behaviour is identical to the MOND limit in low-acceleration environments and suggests that the successes of MOND on the scale of galaxies should be shared by RG.

Adopting RG, rather than MOND, as the modified theory of gravity is advantageous mostly from a theoretical point of view. In MOND, the transition from a Newtonian regime to a regime of modified gravity is driven by the gravitational acceleration generated by the ordinary matter. The acceleration scale for this transition is indeed supported by extended observational evidence (e.g. Milgrom 1983b; Sanders & McGaugh 2002; McGaugh 2004; McGaugh et al. 2016). From a theoretical perspective, however, this feature has made the construction of a covariant formulation of MOND considerably challenging (e.g. Skordis 2009; Bekenstein 2011). On the contrary, the adoption of a scalar quantity, like the density ρ, appears to simplify this task for RG. In addition, as suggested in Matsakos & Diaferio (2016), RG might in principle reproduce the phenomenology properly described by MOND without necessarily explicitly inserting an acceleration scale in the theory.

For simplicity, RG assumes that the permittivity ϵ only depends on the density ρ of ordinary matter. In principle, the gravitational sources are characterised by other scalar quantities, including their total mechanical and thermodynamical energy or their entropy. However, these quantities partly depend on the mass density and we might expect that adopting a more complex dependence of the permittivity ϵ might return a phenomenology that is comparable to the one we investigate here by assuming a simple dependence on ρ.

Clearly, all these issues remain unsettled at this stage of the investigation of RG: suggestions on how they could be properly tackled might originate from a better understanding of the connection between the permittivity ϵ and the scalar field φ appearing in the covariant formulation of RG. Before exploring this connection, however, we need to investigate whether RG is indeed comparable to MOND in the description of the kinematic properties of disk galaxies. This is the task we intend to accomplish here.

As a test case for the quantitative analysis we describe in this work, following Matsakos & Diaferio (2016), we adopt a smooth step function for the gravitational permittivity,

$$\begin{array}{c}\hfill ϵ(\rho )={ϵ}_0+(1-{ϵ}_0)\frac{1}{2}\{\tanh \left[\text{ln}{\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)}^Q\right]+1\},\end{array}$$(5)

which depends on three parameters that we expect to be universal: the permittivity of the vacuum ϵ₀, the power index Q, and the critical density ρ_c. The parameter ϵ₀ is limited in the range of [0, 1] by the definition of ϵ(ρ) and its asymptotic limits in Eq. (2). The parameter Q sets the steepness of the transition between the Newtonian and the RG regimes. The critical density, ρ_c, sets the local density where this transition occurs. Figure 1 shows an example of the gravitational permittivity for different values of Q and for ϵ₀ = 0.25. We emphasise that Eq. (5) is just an arbitrary expression for ϵ(ρ) that we choose here to test the viability of RG. Other expressions for ϵ(ρ), ones that can still increase monotonically with ρ and have the asymptotic limits of Eq. (2), are clearly possible.

We conclude this section with a brief comment on RG, MOND, and electrostatics. RG was inspired by the behaviour of electric fields in matter, but the connection between RG and electrostatics does not go beyond the phenomenological formulation of the modified Poisson Eq. (1). A completely different idea, still based on electrostatics, has instead been developed in a number of papers (Blanchet 2007; Blanchet & Le Tiec 2008a,b, 2009; Blanchet & Heisenberg 2017): the phenomenology described by MOND is interpreted by introducing, in addition to the standard CDM particles, a dark fluid subject to a polarisation in a gravitational field, similarly to the electrostatic polarisation of a dieletric medium. In this dipole dark matter model, the mechanism of gravitational polarisation is guaranteed by the presence of a vector field (Blanchet & Heisenberg 2017). This dipole dark matter model has no connection nor any similarity with RG; moreover, RG, at least at this stage of its development, clearly benefits from a much simpler framework in both its phenomenological and covariant formulations. A different issue is whether RG can indeed describe the observed properties of real systems, as we intend to investigate in the present work.

3. Modelling the rotation curves alone of the DMS galaxies

To test whether RG can describe the dynamics of disk galaxies, we first consider the rotation curves of the 30 published galaxies of the DMS catalogue on their own (Bershady et al. 2010a,b; Westfall et al. 2011a,b; Martinsson et al. 2013a,b).

For an axisymmetric mass density distribution ρ(R, z), the Poisson equation (1) returns the gravitational potential ϕ(R, z) that, in turn, yields the rotation curve,

$$\begin{array}{c}\hfill v(R,z=0)={\left[R\frac{\partial \varphi (R,z)}{\partial R}\right]}^{1/2},\end{array}$$(6)

on the disk plane, z = 0.

We model each disk galaxy with a stellar disk, a stellar bulge, and an interstellar gas disk separated into an atomic and a molecular component. The stellar disk is described by a linear interpolation of the measured radial surface brightness and by an exponentially decreasing density profile along the vertical axis. The stellar bulge is described by a Sérsic profile. The details of our model of the mass distribution are given in Appendix A.

To model the galaxy rotation curves, we estimate the RG potential by numerically solving the Poisson Eq. (1) with a successive over relaxation algorithm described in Appendix B. We then perform the numerical derivative of the potential to obtain the rotation curve from Eq. (6). In our model, the rotation curve depends on two parameters describing the galaxy, namely the disk mass-to-light ratio, Υ, and the disk-scale height, h_z, and on the three parameters of the RG gravitational permittivity: ϵ₀, Q, and ρ_c.

We adopt the same mass-to-light ratio for the bulge and the stellar disk because, in our sample, the galaxy luminosity is dominated by the disk (see Appendix A); assuming a different mass-to-light ratio for the bulge only introduces an additional free parameter without substantially improving the galaxy model. We explore this five-dimensional parameter space with a Monte Carlo Markov chain (MCMC) algorithm.

We assume a Gaussian prior for the two galaxy parameters Υ and h_z and a flat prior for the three RG parameters. Specifically:

1. For the mass-to-light ratio Υ in the K-band, adopted for the measurement of the surface brightness of the DMS galaxies, we use a Gaussian prior; the first moment of the Gaussian is the value derived from the SPS models of Bell & de Jong (2001) applied to the DMS galaxies; for these galaxies the B−K colour ranges from 2.7 (UGC 7244) to 4.2 (UGC 4458) (see Martinsson et al. 2013a, Table 1). We set the second moment of the Gaussian to three times the maximum errors derived from the SPS models for each galaxy. This choice yields a second moment in the range of 0.21–0.36 dex². We set the Gaussian tail to zero where Υ <  0.

2. For the disk-scale height h_z, we adopt a Gaussian prior with mean h_z, SR, where h_z, SR is the disk-scale height derived from the relation, reported in Eq. (A.12) in Appendix A, between the observed disk-scale heights and the disk-scale lengths, inferred from the observations of edge-on galaxies (Bershady et al. 2010b). We set the standard deviation of the Gaussian to three times the errors on h_z, SR, which basically coincide with the intrinsic scatter of the relation (A.12). On average, the error on h_z, SR is 0.11 kpc for the DMS galaxies. We set the Gaussian tail to zero where h_z <  0.

3. For the vacuum permittivity ϵ₀, we adopt a flat prior in the range of [0.10, 1]. In principle, the full allowed range is [0, 1]; however, we do not explore values smaller than 0.10 because the boost of the gravitational field would yield unphysically large rotation velocities, of the order of ∼600−1000 km s⁻¹.

4. For Q, we adopt a flat prior in the range of [0.01, 2]. This parameter regulates the steepness of the transition between the Newtonian and the RG regimes. Our range explores from very smooth (Q = 0.01) to very steep transitions (Q = 2).

5. For log₁₀ρ_c, we adopt a flat prior in the range of [ − 27, −23], with the critical density ρ_c in units of g cm⁻³. This range includes the two extreme values −27 and −24 considered by Matsakos & Diaferio (2016).

We adopt the Metropolis–Hastings acceptance criterion in our MCMC algorithm. The random variate x at step i + 1 is drawn from the probability density G(x|x_i), which depends on the random variate x_i at the previous step. For the probability density G(x|x_i), we adopt a multi-variate Gaussian density distribution with mean value x_i; its multiple standard deviations are 1/3 the standard deviation of the Gaussian priors, for Υ and h_z, and 10% of the prior ranges, for the three RG parameters. In our case, x is the five-dimensional vector x = (Υ,h_z,ϵ₀,Q, log₁₀ ρ_c). We adopt the likelihood,

$$\begin{array}{c}\hfill \mathcal{L}(\mathbf{x})=exp[-\frac{{\chi }_{\mathrm{red},\mathrm{RC}}^2(\mathbf{x})}{2}],\end{array}$$(7)

where,

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red},\mathrm{RC}}^2(\mathbf{x})=\frac{1}{n_{\mathrm{dof},\mathrm{RC}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{RC}}}}\frac{{[v_{\mathrm{mod}}(R_i;\mathbf{x})-v_{\mathrm{data}}(R_i)]}^2}{v_{\mathrm{data},\mathrm{err}}^2(R_i)},\end{array}$$(8)

N_RC is the number of data points of the rotation curve, v_data are the velocity measures at the projected distance R_i with their uncertainty v_data, err, v_mod is estimated with Eq. (6) and n_dof, RC = N_RC − 5 is the number of degrees of freedom. If p(x) is the product of the priors of the components of x, the Metropolis–Hastings ratio is:

$$\begin{array}{c}\hfill A=\frac{p(\mathbf{x})\times \mathcal{L}(\mathbf{x})}{p({\mathbf{x}}_i)\times \mathcal{L}({\mathbf{x}}_i)}\frac{G(\mathbf{x}|{\mathbf{x}}_i)}{G({\mathbf{x}}_i|\mathbf{x})}·\end{array}$$(9)

If A ≥ 1, we set x_i+1 = x, otherwise we set either x_i+1 = x, with probability A, or x_i+1 = x_i, with probability 1 − A.

For the chain starting points, we adopt the values found by the DMS collaboration for Υ and h_z (see Angus et al. 2015, Table 1); for the three RG parameters ϵ₀, Q and log₁₀ρ_c, we set 0.30, 1.00, and −24.0, respectively. We run the MCMC algorithm for 19 000 steps, after a burn-in chain of 1000 steps. This number of steps guarantees the achievement of the chain convergence. We check the chain convergence using the Geweke diagnostic (Geweke 1992): we compare the means of the first 10% of the chain steps, after the burn-in, with the last 50% of the chain steps; we compare the means with a Gaussian test, adopting the standard deviations of the two portions of the chains as the errors on the two means. In every case, the Gaussian test shows that the two means coincide at a significant level larger than 5%, suggesting that the chains converge.

As an example, Fig. 2 shows the posterior distributions for four galaxies. The posterior distributions for the remaining galaxies are qualitatively very similar. The posterior distributions show a single peak and we can thus adopt the medians of the posterior distributions as our parameter estimates; the range between the 15.9 and the 84.1 percentiles, thus including 68% of the posterior cumulative distribution centred on the median, defines our 1σ uncertainty range on the parameter estimates. Table 1 lists the medians of the parameters and their associated uncertainties.

We use these parameters to compute our rotation curve models. We collect all the figures showing our results in Appendix D. We arrange the figures by galaxy, so that the outcomes of the various analyses we perform here can be compared more easily.

The rotation curves estimated in this section are shown in sub-panels (d) of Figs. D.1 and D.2 as blue solid lines. The red dots with error bars are the DMS data. The vertical lines show the size of the bulge we adopt. In some galaxies, the presence of the bulge produces, at small radii, a relevant spike in the modelled rotation curve that is not present in the data. Other than these cases, the observed rotation curves are modelled relatively well. Because we model the surface brightness of the disk with a linear interpolation of the data (Appendix A), the model rotation curves capture some of the features appearing in the measured rotation curves that have a correspondence in the surface brightness profile of the disk; the galaxies UGC 1635, UGC 4555, UGC 6903, or UGC 7244 are some examples of this correspondence. Nevertheless, the rotation curves of a few galaxies still have some features that are not described by the model, for example UGC 4036, UGC 4622 or UGC 8196.

The ${\chi }_{\text{red},\text{RC}}^2$ of Eq. (8) are listed in Table 1; they quantify the quality of the RG model. In most cases, the large values of ${\chi }_{\text{red},\text{RC}}^2$ originate from a possible underestimation of the error bars of the data, as suggested by the visual inspection of the sub-panels (d) of Figs. D.1 and D.2, rather than from an inappropriate modelling.

Table 1 also lists the mass-to-light ratio Υ_SPS that we estimate for each galaxy with the SPS models of Bell & de Jong (2001) from the B−K colours listed in Table 1 of Martinsson et al. (2013a)³. Specifically, we adopt the SPS model with a mass-dependent formation epoch with bursts and a scaled Salpeter initial mass function (IMF) (see Bell & de Jong 2001, Table 1). According to Bell & de Jong (2001), this model better reproduces (1) the trends in colour-based stellar ages and metallicities (Bell & Bower 2000), (2) the decrease in the colour-stellar mass-to-light ratio slope caused by modest bursts of star formation, and (3) it has an IMF consistent with maximum disk constraints.

To somehow quantify the uncertainty on Υ_SPS, for the error bar we adopt the range covered by all the different models investigated by Bell & de Jong (2001) based on a scaled Salpeter IMF (see Bell & de Jong 2001, Table A3). The resulting uncertainties are asymmetric and, except for three galaxies, the lower limit of the error bar is 0 because, in these cases, our preferred SPS model yields the smallest mass-to-light ratio among all the other models with the same IMF. We neglect the models that adopt different IMFs, under the assumption that the scaled Salpeter IMF can reasonably be considered universal (Bell & de Jong 2001). We do not estimate the SPS mass-to-light ratio for the galaxy UGC 3997, since Martinsson et al. (2013a) does not provide its B−K colour.

The two mass-to-light ratios for each galaxy, namely, the Υ_SPS and Υ derived from our RG model, are shown in the left panel of Fig. 3. In Table 1 we list the difference between these ratios in units of the uncertainties on the mass-to-light ratios. Our Υ’s nearly cover the same range as the Υ_SPS’s and the two mass-to-light ratios for the same galaxy tend to agree with each other: for 23 out of 29 galaxies, they agree within 2σ; for an additional five galaxies, they agree within 2.5–3.5σ; and for only one galaxy, they are discrepant at more than 6.5σ.

The disk-scale height h_z derived in RG tends to be larger than the scale height h_z, SR calculated with Eq. (A.12), as shown in the left panel of Fig. 4. Table 1 also lists their ratios. Specifically, the RG h_z is larger than h_z, SR for 20 galaxies out of 30, and for nine galaxies h_z/h_z, SR ≥ 2; among these nine galaxies, two have h_z/h_z, SR ≥ 5. The ten galaxies with h_z/h_z, SR <  1 have thinner disks than expected because their gravitational potential wells are shallow: in fact, their central disk surface brightness I_d0 and their disk-scale length h_R are among the smallest in the sample, similarly to their rotation velocity and vertical velocity dispersion profiles.

4. Modelling the rotation curves and the vertical velocity dispersion profiles

The DMS sample enables a more comprehensive investigation of the dynamics of disk galaxies because, in addition to the rotation curves, we have a measurement of their stellar vertical velocity dispersion profiles. In fact, the DMS galaxies are close to face-on: as illustrated in Appendix A, based on the Tully–Fisher relation (Eqs. (A.5) and (A.6)), the estimated inclinations of the DMS galaxies are in the range of 5.8−45.3 deg (see Martinsson et al. 2013a, Table 5). Therefore, combining the two kinematic pieces of information can provide unique constraints on the dynamical properties of the galaxies and can represent a stringent test for theories of modified gravity.

To model the stellar vertical velocity dispersion profile, we use the fact that our axisymmetric model of the galaxy, illustrated in Appendix A, is described by a two-integral distribution function f(E, L_z) and thus the velocity dispersions, σ(R, z), along the vertical and radial axes, z and R, coincide (Nagai & Miyamoto 1976; Nipoti et al. 2007). The system also satisfies the Jeans equation,

$$\begin{array}{c}\hfill \frac{\partial [\rho (R,z){\sigma }^2(R,z)]}{\partial z}+\rho (R,z)\frac{\partial \varphi (R,z)}{\partial z}=0,\end{array}$$(10)

where ρ(R, z) is the stellar density. The observed vertical velocity dispersion profile, weighted by the local stellar surface density Σ_*(R), is

$$\begin{array}{c}\hfill {\sigma }_z^2(R)=\frac{1}{{\mathrm{\Sigma }}_{\ast }(R)}{\displaystyle {\int }_{-\infty }^{+\infty }}\rho (R,z){\sigma }^2(R,z)\mathrm{d}z.\end{array}$$(11)

By considering the contribution of the disk alone to Σ_*(R) and ρ(R, z), thus neglecting the luminosity contribution of the bulge (see Appendix A), the stellar surface mass density profile is:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\ast }(R)={\displaystyle {\int }_{-\infty }^{+\infty }}{\rho }_{\ast }(R,z)\mathrm{d}z={\mathrm{{\rm Y}}}_{\mathrm{d}}{I}_{\mathrm{d}}(R),\end{array}$$(12)

where ρ(R, z) = ρ_*(R, z) is our disk model, namely the surface brightness profile I_d(R) multiplied by exp(−|z|/h_z)/(2h_z) (Appendix A.1). Thus, combining Eqs. (11), (12), and (10) yields:

$$\begin{array}{c}\hfill {\sigma }_z^2(R)=\frac{1}{h_z}{\displaystyle {\int }_0^{+\infty }}[{\displaystyle {\int }_z^{+\infty }}exp(-\frac{|z^{\prime }|}{h_z})\frac{\partial \varphi (R,z^{\prime })}{\partial z^{\prime }}\mathrm{d}{z}^{\prime }]\mathrm{d}z.\end{array}$$(13)

We model the two kinematic profiles, the rotation curve (Eq. (6)), and the vertical velocity dispersion profile (Eq. (13)), with the same MCMC algorithm illustrated in Sect. 3. We adopt the same priors of Sect. 3.

The MCMC analysis considers the two kinematic profiles at the same time. We thus define the single figure of merit,

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red},\mathrm{tot}}^2(\mathbf{x})=\frac{{\chi }_{\mathrm{RC}}^2(\mathbf{x})+{\chi }_{\mathrm{VVD}}^2(\mathbf{x})}{n_{\mathrm{dof},\mathrm{tot}}},\end{array}$$(14)

where ${\chi }_{\text{RC}}^2$ is Eq. (8) multiplied by n_dof, RC, and

$$\begin{array}{c}\hfill {\chi }_{\mathrm{VVD}}^2(\mathbf{x})={\displaystyle \sum _{j=1}^{N_{\mathrm{VVD}}}}\frac{{[{\sigma }_{z,\text{mod}}(R_j,\mathbf{x})-{\sigma }_{z,\text{data}}(R_j)]}^2}{{\sigma }_{z,\text{data,err}}^2(R_j)}\end{array}$$(15)

is derived from the vertical velocity dispersion profile, with N_VVD the number of data points of the vertical velocity dispersion profile, σ_z, data the velocity dispersion measured at the projected distance R_j with their uncertainty σ_{z, data,err}, and σ_z, mod the velocity dispersion model estimated with Eq. (13); finally n_dof, tot = N_RC + N_VVD − 5 is the total number of degrees of freedom.

We estimate the vertical velocity dispersion error bars, σ_{z, data,err}, by summing their random (∼0.1−10 km s⁻¹) and systematic (∼1−10 km s⁻¹) uncertainties in quadrature (Martinsson et al. 2013b).

Figure 5 shows the posterior distributions for the same four galaxies shown in Fig. 2 for comparison. For the remaining galaxies, the posterior distributions are qualitatively similar. Because the posterior distributions show a single peak, we adopt the medians of the posterior distributions as our parameter estimates as in Sect. 3.

Table 2 lists the medians and the 15.9 and the 84.1 percentiles of the posterior distributions. We adopt these percentiles as the 1σ uncertainty range. The blue solid lines in the sub-panels (e) and (f) in Figs. D.1 and D.2 show the model rotation curves and the vertical velocity dispersion profiles with the median parameters listed in Table 2. The red dots with error bars show the DMS data. As in Sect. 3, the features of the observed rotation curves are captured by the models in most cases although some discrepancies remain. Remarkably, the modelling of the rotation curves generally improves, like in UGC 1087, UGC 3997, and UGC 9837.

The right panel of Fig. 4 shows that the disk-scale heights h_z required to model the rotation curves and the vertical velocity dispersion profiles are substantially smaller than the disk-scale heights h_z, SR derived from Eq. (A.12), suggested by the observations of edge-on galaxies. In sub-panels (f) of Figs. D.1 and D.2, the cyan solid lines are the vertical velocity dispersion profiles when we adopt the parameters of Table 2 except for h_z replaced by h_z, SR. Comparing the left and the right panels of Fig. 4 confirms that including the vertical velocity dispersion profiles in the modelling is responsible for requiring substantially thinner disks than h_z, SR; in fact, modelling the rotation curve alone may actually require disks thicker than expected from Eq. (A.12).

This comparison supports the conclusion that the discrepancies between h_{z, RC + VVD} and h_z, SR are derived from the fact that the two scale heights are inferred from two different stellar populations: a younger stellar population dominating the observed vertical velocity dispersion, hence h_{z, RC + VVD}, and an older stellar population dominating the surface brightness in the edge-on galaxy sample, hence h_z, SR (Aniyan et al. 2016).

The general agreement between our results and those found by Angus et al. (2015) for MOND also confirms that, as anticipated in the introduction, if this disagreement between h_{z, RC + VVD} and h_z, SR is derived from an overlooked observational bias, it does not suggest a possible tension between the data and the theory of gravity, either MOND or RG.

This observational bias in the vertical velocity dispersion measure brings on consequences for the estimate of the mass-to-light ratio. In fact, the velocity dispersion increases with the disk thickness and with the intensity of the gravitational field originated by the baryonic mass. By reducing the disk thickness, a larger mass must be attributed to the baryonic component to reproduce the observed velocity field. Therefore, our mass-to-light ratios are larger than the DMS values, as also occurs in MOND (Angus et al. 2015). In fact, our mass-to-light ratios agree with the MOND ratios within 2σ for 29 out of 30 galaxies and within 2.13σ for UGC 11318.

Our estimated mass-to-light ratios are corroborated by the comparison with the SPS values (Sect. 3). The right panel of Fig. 3 compares the mass-to-light ratios required to model the rotation curves and the vertical velocity dispersion profiles with the mass-to-light ratios derived from the SPS model. The RG mass-to-light ratios span the same range (∼0.3−1.2 M_⊙/L_⊙) as the RG ratio derived by modelling the rotation curve alone (left panel of Fig. 3), which nearly coincides with the range covered by the mass-to-light ratios derived from the SPS models. For each individual galaxy, our estimated mass-to-light ratio agrees with the mass-to-light ratio derived from the SPS model most of the time: for 20 out of 29 galaxies, the two values agree within 2σ; for six galaxies, they agree within 2–3σ; and for the remaining three galaxies, they are consistent within 3–5σ. The largest discrepancy, 4.6σ, occurs for UGC 1908; this discrepancy is smaller, however, than the 6.7σ discrepancy found in Sect. 3 with the modelling of the rotation curve alone.

On the contrary, for 70% of the galaxy sample, the mass-to-light ratios estimated by the DMS collaboration are at least a factor of 2 smaller than the SPS values and their disks turn out to be sub-maximal (Angus et al. 2015; Aniyan et al. 2016). Despite these discrepancies, the large error bars associated with the DMS mass-to-light ratios make them agree within 3σ with the SPS values (see Angus et al. 2015, Table 1).

4.1. The observational bias in the vertical velocity dispersion profile

As mentioned above, the estimate of the vertical velocity dispersion profile is likely to be affected by an observational bias. Therefore, the disk thickness obtained by modelling this profile might give a severe underestimation of the real value. Here, we want to quantify how correcting the vertical velocity dispersions by an appropriate factor can give disk thicknesses that are consistent with the observations of edge-on galaxies and with mass-to-light ratios still consistent with the SPS models.

Aniyan et al. (2016) model what we would observe in the disk of a typical external face-on galaxy by selecting a sample of giant stars simulated in the Milky Way from the on-line Besançon model (Robin et al. 2003). This model describes the main structural components and stellar populations of the Milky Way and includes recipes for galactic reddening, star formation history, and dynamical evolution (Robin et al. 2003; Aniyan et al. 2016). For their simulated star sample, Aniyan et al. (2016) choose giant stars of spectral type G8III – K5III, with luminosity in the range of −3 ≤ M_V ≤ +3, and colour in the range of 0.8 ≤ B−V ≤ 1.8, within a cylindric volume, centred on the Sun and perpendicular to the Galactic disk, of radius 2 kpc and height 10 kpc. They assume a star formation and a dynamical history similar to the solar neighbourhood.

From the simulated sample, they extract the vertical velocity dispersion integrated vertically through the disk by considering either the entire stellar population (σ₁) or the dynamically hotter and older giant component alone (σ₂). They estimate that η = σ₂/σ₁ = 1.55 ± 0.02. Now, observationally, the disk-scale heights are inferred from measures of the surface brightness that is dominated by the old hot giants, which have a larger vertical velocity dispersion, whereas the younger stellar population mostly contributes to the spectral signal used to derive the observed velocity dispersion. Therefore, according to the result of Aniyan et al. (2016), assuming a single stellar population to interpret the vertical velocity dispersions and the disk-scale height at the same time introduces a bias that originates an underestimate of the disk thickness.

Here, we estimate this observational bias with the ratio between the expected vertical velocity dispersions computed with the parameters inferred from the rotation curves alone (Table 1) and the measured vertical velocity dispersion profiles. We adopt the following procedure: (1) for each galaxy, we use Eq. (13) to compute the expected vertical velocity dispersion, according to the parameters of Table 1; (2) at each radial coordinate, we compute the ratio between the expected and the measured vertical velocity dispersions and we estimate the mean of these ratios along the radial profile of each of the 30 galaxies; (3) we then estimate the mean, η′, of these 30 means and their standard deviation. Thus, we obtain η′ = 1.63 ± 0.65, which agrees within 0.12σ with the value 1.55 ± 0.02 found by Aniyan et al. (2016).

To quantify the effect of this bias on the estimate of the disk thickness, we now repeat the analysis presented in Sect. 4 for five galaxies where we artificially increase the vertical velocity dispersion profile by a factor of 1.63. We choose five galaxies with a h_z-to-h_z, SR ratio smaller than 0.5: UGC 1635, UGC 3091, UGC 4107, UGC 4555, and UGC 9965. We also model the kinematics of these galaxies in MOND where, as free parameters, we only have the galactic parameters Υ and h_z. For the MCMC analysis, we use the same priors as above.

To model the galaxies in MOND, we derive the MOND potential by solving the Poisson equation in the QUMOND formulation of MOND (Milgrom 2010):

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\varphi =\mathrm{\nabla }·\left[\nu \left(\frac{|\mathrm{\nabla }{\varphi }_{\mathrm{N}}|}{a_0}\right)\mathrm{\nabla }{\varphi }_{\mathrm{N}}\right],\end{array}$$(16)

where ϕ_N is the Newtonian potential, a₀ = 1.2 × 10⁻¹⁰ m s⁻² = 3600 kpc⁻¹ (km  s⁻¹)² is the MOND critical acceleration, and ν is the interpolating function regulating the transition between the Newtonian and the MOND regimes. We use the “simple ν-function”,

$$\begin{array}{c}\hfill \nu (y)=\frac{1}{2}(1+\sqrt{1+\frac{4}{y}})·\end{array}$$(17)

Tables 3 and 4 list the medians and percentile ranges of the posterior distributions in RG and QUMOND, respectively. As expected, h_z increases by a factor between 2.07 (found for UGC 3091 in QUMOND) and 2.89 (found for UGC 4555 in QUMOND) with respect to the values obtained in Sect. 4 (see Tables 2–4). For UGC 1635 and UGC 3091, the disk-scale heights are still smaller than the values expected from the observations of edge-on galaxies, but their h_z-to-h_z, SR ratios still increase to values greater than 0.5.

The agreement between the RG mass-to-light ratios estimated in the current section and the SPS values worsens compared to Sect. 4, but for four out of five galaxies, it is within 4σ. QUMOND mass-to-light ratios tend to be closer to the SPS values than the RG mass-to-light ratios: their agreement with their SPS expectations is within 3σ for four out of five galaxies. However, since the QUMOND mass-to-light ratios uncertainties are smaller than in RG, the agreement between the QUMOND mass-to-light ratios and the SPS values is formally worse than in RG for two out of five galaxies (UGC 1635 and UGC 3091).

The sub-panels (i)–(l) of Fig. D.2 show the RG and QUMOND models of the rotation curves and vertical velocity dispersion profiles (blue solid lines). The vertical velocity dispersion data are increased by the factor 1.63 (red dots with error bars). RG describes the rotation curves slightly better than QUMOND whereas both theories describe the vertical velocity dispersion profiles equally well. Clearly, this better performance of RG follows from the version of RG we adopt at this stage that, for each galaxy, has three more free parameters than QUMOND.

Our QUMOND results of these five galaxies are comparable to the results of Angus et al. (2015). Specifically, we also find that the QUMOND rotation curves of UGC 3091 and UGC 9965 properly describe the data, whereas the QUMOND rotation curves of UGC 1635, UGC 4107 and UGC 4555 tend to underestimate the inner profile and to overestimate the outer profile. Our vertical velocity dispersions are also comparable to Angus et al. (2015), although their models slightly uderestimate the most inner data point of UGC 1635 and UGC 3091.

Figure 6 compares the parameters estimated with RG with the original values of the vertical velocity dispersion profile, σ_z(R) (Sect. 4), with the parameters obtained with σ_z(R) increased by the factor of 1.63. These figures show that the three RG parameters are not systematically affected by the increase of σ_z(R), whereas both Υ and h_z tend to be larger. Figure 7 compares Υ and h_z obtained in RG and QUMOND with the σ_z(R) values increased by the factor 1.63. The disk-scale heights are comparable in the two theories of gravity whereas the mass-to-light ratios in RG are systematically larger than in QUMOND. This systematic difference will be relevant in describing the RAR we illustrate below.

4.2. On the error bars of the rotation curves

We conclude this section with a digression on the rotation curve error bars. In their modelling of the DMS data with MOND, rather than maintaining the original DMS error bars as we do here, Angus et al. (2015) arbitrarily increase the error bars from ∼1−5 km s⁻¹ to 10 km s⁻¹. Angus et al. (2015) were concerned by the possibility that disk warping could introduce systematic errors which can indeed affect the rotation curves; by this approach, Angus et al. (2015) clearly intend to increase the statistical weight of the vertical velocity dispersion profiles which are more dependent on the disk thickness.

We verified that this modification is in fact unnecessary. By setting the rotation curves error bars to 10 km s⁻¹ and performing once again the MCMC analysis described above, we find that the new disk-scale heights, h_z, substantially agree with our original results: the distributions of the ratios between the new and the original h_z is peaked around ∼1, with median $1.{09}_{-0.20}^{+0.11}$, where the uncertainties are the 15.9 and 84.1 percentiles. Similarly, the mass-to-light ratios remain unbiased, with their ratios between the new and the original values having median $0.{95}_{-0.11}^{+0.09}$.

5. A universal combination of the RG parameters

In the formulation of RG, ϵ₀, Q and ρ_c are universal parameters. Here we show that, in principle, a single set of these parameters could indeed be able to model the dynamics of our entire sample of disk galaxies.

For this task, the ideal approach would be to use the MCMC algorithm to explore the 63-dimensional space of these three RG parameters and the 2 × 30 parameters describing the galaxies, namely the mass-to-light ratios Υ and the disk-scale heights h_z. However, this approach requires an extraordinary computational effort which, at this stage of our investigation of the RG viability, appears unreasonably large.

We prefer a simpler strategy. Our analyses above show that modelling each individual galaxy by keeping the three RG parameters free returns values of these parameters that are roughly consistent from galaxy to galaxy. In fact, the ϵ₀, Q and log₁₀ρ_c values listed in Table 2, whose distributions are shown in Fig. 8, have mean and standard deviations ϵ₀ = 0.56 ± 0.19, Q = 0.92 ± 0.24, and log₁₀(ρ_c/g cm⁻³) = − 25.30 ± 0.70. These standard deviations are either smaller than or comparable to the mean uncertainties of the values listed in Table 2, which are 0.16, 0.71, and 1.22, for ϵ₀, Q and log₁₀ρ_c, respectively. We can thus reasonably conclude that the different values that we find for different galaxies can in principle be ascribed to statistical fluctuations.

Therefore, to check whether a single set of RG parameters could be able to describe the entire galaxy sample, rather than performing the computationally demanding exploration of the 63-dimensional parameter space, we assume that the values of the mass-to-light ratios Υ and the disk-scale heights h_z for each galaxy, estimated in our previous analysis, are appropriate, and we only explore the 3-dimensional space of the parameters ϵ₀, Q and log₁₀ρ_c for the entire galaxy sample at the same time.

We adopt the priors for the three RG parameters as in Sects. 3 and 4. The figure of merit is now:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red}}^2(\mathbf{x})={\displaystyle \sum _{i=1}^{N_{\text{gal}}}}{\chi }_{\text{red,tot},i}^2(\mathbf{x}),\end{array}$$(18)

where x = (ϵ₀,Q,log₁₀ρ_c), N_gal is the number of DMS galaxies and ${\chi }_{\text{red,tot},i}^2(\mathbf{x})$ is Eq. (14) where the number of free parameters is now equal to three instead of five.

As in Sects. 3 and 4, we run the MCMC for 19 000 steps in addition to the 1000 burn-in steps⁴. To better assess the convergence of the chains we run three chains with three different starting points and we test their convergence with the variance ratio method of Gelman & Rubin (1992), described in Appendix C. We find that 13 000 steps are already sufficient to have the chains converging. We also test that each chain converges according to the Geweke diagnostic (Geweke 1992).

Figure 9 shows the posterior distributions of the three parameters. The green dots show the median values, and the yellow, red, and black curves show the 1, 2, and 3σ contour levels. The medians and 68% confidence intervals are ${ϵ}_0=0.{661}_{-0.007}^{+0.007}$, $Q=1.{79}_{-0.26}^{+0.14}$ and ${log}_{10}{\rho }_{\mathrm{c}}=-24.{54}_{-0.07}^{+0.08}$. The purple points show the means of the three distributions shown in Fig. 8; the error bars show the mean errors listed in Table 2.

As expected, the posterior distributions in Fig. 9 show that considering all the DMS galaxies at the same time provides much tighter constraints on the RG parameters than in our previous analyses. Because of the large errors found in Sect. 4, the universal parameters estimated here are consistent with our previous analyses: the means of the distributions of Fig. 8 agree within 0.63, 1.19 and 0.62 σ with the medians, found here, of ϵ₀, Q and log₁₀ρ_c, respectively.

However, the agreement between the data and the models of the rotation curves and the vertical velocity dispersion profiles derived with the mass-to-light ratios and disk-scale heights listed in Table 2 and the universal combination of RG parameters found here worsens with respect to the agreement obtained with the individual RG parameters found in Sect. 4 as shown in sub-panels (g) and (h) of Figs. D.1 and D.2. Figure 10 compares the reduced χ² found in Sect. 4 with the reduced χ² found here. Despite the presence of a number of galaxies whose χ² substantially increases, the bulk of the sample maintains its χ² close to, albeit still larger than, the original χ².

The increase of the χ² compared to Sect. 4 is mainly due to the rotation curves. In fact, all the vertical velocity dispersion profiles are rather well interpolated with this universal combination of RG parameters: their χ² are close to those of Sect. 4, except for a few outliers, like UGC 3701 and UGC 3997. On the contrary, the rotation curves of only about half of the sample are still well described with this unique combination of RG parameters (e.g. UGC 1081, UGC 1635, UGC 4036, UGC 4380 and UGC 9965), whereas the models worsen for the remaining sample.

The good agreement of this universal combination of RG parameters with the parameters of Sect. 4 and the somewhat limited worsening of the χ² shown in Fig. 10 suggest that finding a unique set of RG parameters that accurately describes the kinematics of the DMS galaxies might be feasible. This simple exercise in fact appears to indicate that properly exploring the full 63-dimensional parameter space might return substantially smaller χ² with still reasonable, albeit different from the results of Sect. 4, mass-to-light ratios and disk-scale heights of the galaxies. We might also expect that the unique set of RG parameters will be statistically equivalent to the set we find here.

6. The radial acceleration relation

To test the viability of RG as a gravity theory describing the dynamics of disk galaxies, we need to consider an additional relevant observational piece of evidence that very clearly quantifies the mass discrepancy on galaxy scales: the RAR. McGaugh et al. (2016) and Lelli et al. (2017) pointed out that the observed centripetal acceleration traced by the rotation curves,

$$\begin{array}{c}\hfill g_{\text{obs}}(R)=v_{}^2\text{obs}(R)/R,\end{array}$$(19)

is tightly correlated with the Newtonian acceleration g_bar(R) due to the baryonic matter distribution alone.

McGaugh et al. (2016) found that the function:

$$\begin{array}{c}\hfill g_{\text{obs}}(R)=\frac{g_{\text{bar}}(R)}{1-exp(-\sqrt{\frac{g_{\text{bar}}(R)}{g_{†}}})}\end{array}$$(20)

provides a good fit for the entire SPARC sample, made of 153 galaxies (Lelli et al. 2016b). The fit has only one free parameter g_† whose single value g_† = 1.20 ± 0.02 (random) ±0.24 (systematic) ×10⁻¹⁰ m s⁻² is appropriate for all the galaxies in the sample. This value also is consistent with the MOND acceleration scale a₀. The asymptotic limit of Eq. (20) for small g_bar returns the acceleration in the MOND regime $g_{\text{obs}}\sim \sqrt{g_{\text{bar}}(R)a_0}$.

For the galaxies of the SPARC sample, this correlation has a relatively small root-mean-square scatter of 0.13 dex, mostly due to possible variations of the stellar mass-to-light ratio from galaxy to galaxy and to observational uncertainties (McGaugh et al. 2016; Lelli et al. 2017). Indeed, for this result, McGaugh et al. (2016) adopt a single 3.6 μm mass-to-light ratio of 0.50 M_⊙/L_⊙ for all the galaxy disks, under the assumption that the stellar mass-to-light ratio does not vary much in this band (McGaugh & Schombert 2014, 2015; Meidt et al. 2014; Schombert & McGaugh 2014). Similarly, for the galaxy bulges, which are present in 31 out of 153 galaxies, McGaugh et al. (2016) adopt a mass-to-light ratio equal to 0.70 M_⊙/L_⊙.

To explore the intrinsic scatter of the RAR due to the mass-to-light ratio variations, Li et al. (2018) fit galaxy mass-to-light ratios to individual galaxies with Eq. (20) and g_† fixed to 1.2 × 10⁻¹⁰ m s⁻². They find a RAR tighter than McGaugh et al. (2016), with an intrinsic root-mean-square scatter of 0.057 dex and mass-to-light ratios generally consistent with the SPS model predictions.

Here, we estimate the RAR for our DMS sample. The observed acceleration, g_obs, is directly derived from the measured rotation curve according to Eq. (19). The Newtonian acceleration attributed to the baryonic matter alone, g_bar = | − ∂ϕ/∂R|, is derived from the numerical solution of the Newtonian Poisson Eq. (3) where the density, ρ(R, z), is Eq. (A.10). For the galaxy disk-scale height, h_z, which appears in ρ(R, z), we use the values derived from Eq. (A.12) inferred from the observations of edge-on galaxies (see Appendix A.3). For the numerical solution of the Poisson equation, we adopt the successive over relaxation algorithm described in Appendix B.

To estimate the mass density ρ(R, z) from the galaxy surface brightness, we adopt the mass-to-light ratios derived with our MCMC analysis (Sect. 4) and listed in Table 2. As discussed above, for 26 galaxies out of 29, these values agree within 3σ with the values of the SPS models (see Bell & de Jong 2001, Table 1) and for all the galaxies, the difference is within 5σ.

Red symbols and error bars in both panels of Fig. 11 show the RAR of all the DMS galaxies in our sample. The black curve is Eq. (20). The blue curves in the left panel of Fig. 11 are the RAR of each DMS galaxy obtained from the RG parameters, the mass-to-light ratios and the disk-scale heights derived from our MCMC analysis of the rotation curves and vertical velocity dispersion profiles (Sect. 4). The blue curves are not fits to the observed RAR, but just the relations between the Newtonian g_bar and the expected RG centripetal acceleration g_obs based on the galaxy parameters estimated with our previous analysis. The dashed line shows the relation g_obs = g_bar for comparison.

We also estimate the RAR expected in MOND, adopting the QUMOND formulation described in Sect. 4.1. For QUMOND, we adopt the mass-to-light ratios Υ and the disk-scale heights h_z that we derive for RG in Sect. 4. These mass-to-light ratios agree within 2σ with the values estimated by Angus et al. (2015), who model the rotation curves and vertical velocity dispersions in QUMOND with the simple interpolating function. Similarly, their values of h_z agree with our estimates within 1σ. The QUMOND RAR curves are shown as green solid lines in the right panel of Fig. 11.

RG correctly reproduces the asymptotic limits of the observed RAR and, on average, it interpolates the data, although it tends to underestimate the relation (20) at low g_bar; on the contrary, QUMOND properly reproduces the shape of the RAR relation (20) along the entire range of g_bar. The fact that RG slightly underestimates the observed RAR while at the same time providing a good fit to the kinematics of the individual galaxies, as shown in the previous sections, suggests that RG might attribute more mass to the luminous matter than QUMOND. In fact, Fig. 12 shows that the RG mass-to-light ratios are systematically larger than in QUMOND, although the mass-to-light ratios in the two models agree with each other within 2σ. This result is consistent with the left panel of Fig. 7 of Sect. 4.1.

A more serious issue for RG is the scatter of the curves along the g_obs axis. Figure 13 shows the distributions of the deviations of each curve from Eq. (20). We only consider the deviations of each curve from Eq. (20) within the horizontal axis range of log₁₀[g_bar(m s⁻²)] = [−11.28,−8.81] covered by the data. This approach makes the comparison with the data more sensible. In passing, we note that we use the disk-scale heights h_z, SR for the data and our estimated h_z for the models. These values can be different by a factor of two, as shown in the right panel of Fig. 4 and in Table 2. However, adopting, for the models, h_z, SR rather than h_z, leaves the distributions of the deviations basically unaffected.

The width of the residual distribution in Fig. 13 for the data, which quantifies the observed scatter of the RAR, is 0.12 dex. We estimate this scatter by removing four outlying points, clearly visible in the right part of both panels of Fig. 11. These points belong to the galaxies UGC 1081, UGC 1862, UGC 3997, and UGC 6903, and they correspond to the innermost point of their rotation curves; these values are 17.00, 0.05, 1.19, and 17.69 km s⁻¹, which are unusually small. If we include these four points, the root-mean-square scatter increases from 0.12 dex to 0.32 dex. The observed scatter of the DMS sample is thus comparable to the value 0.13 dex found by McGaugh et al. (2016) and Lelli et al. (2017).

The widths of the distributions of the residuals for the RG and QUMOND models shown in Fig. 13 are 0.11 and 0.017 dex, respectively. We can identify these widths with the intrinsic scatter of the RAR predicted by the two models. The small intrinsic scatter predicted by QUMOND, consistent with the expectations (Lelli et al. 2017; Brada & Milgrom 1995), is almost an order of magnitude smaller than the RG scatter. MOND actually appears with different formulations: in the version of modified inertia, which modifies the Newtonian second law of dynamics, the intrinsic scatter is predicted to be zero if the orbits are circular (Milgrom 1994); similarly, in the version of modified gravity, which modifies the Poisson equation, like QUMOND does, the intrinsic scatter is predicted to be zero only for spherically symmetric systems (Bekenstein & Milgrom 1984), whereas in flat systems, like disk galaxies, a small not-null intrinsic scatter should appear.

To investigate the nature of the intrinsic scatter predicted by RG, we explore the possible correlation of the RAR residuals with the global and the radially-dependent properties of the galaxies. We plot these residuals in Figs. 14 and 15. The first column shows, in cyan, the residuals of the RG models. The second and the third columns show the residuals for QUMOND, in green, and the DMS data, in pink⁵.

Table 5 lists the Kendall statistic τ (Kendall 1938) and the Spearman statistic ρ (Spearman 1904) with their corresponding p-values, namely, the significance levels of the lack of correlation. For large size N of the sample, the density distributions of the random variates τ and ρ are excellently approximated by the Gaussian distribution with zero mean and variance (4N + 10)/(9N² − 9N) and 1/(N − 1) for τ and ρ, respectively (Best 1973; Kendall 1975). In our analysis of RG, N = 4560 and, assuming Gaussian distributions, the corresponding standard deviations σ’s are ⟨τ²⟩^1/2 = 0.0099 and ⟨ρ²⟩^1/2 = 0.0148.

We can interpret the results of Table 5 by arbitrarily choosing the threshold log₁₀p = −3: p-values smaller than the threshold of p = 10⁻³ indicate that the listed values of τ or ρ have a probability smaller than 0.1% of occurring by chance for an uncorrelated sample. For this probability, the values |τ|> 0.049 and |ρ|> 0.033 are thus more than 3.3σ away from the expected means ⟨τ⟩ = 0 and ⟨ρ⟩ = 0. For RG, the only two parameters that have p-values larger than 10⁻³, namely, −log₁₀p <  3, and therefore their uncorrelation with the residuals is statistically significant, are the central surface brightness I_d0 of the disk and its scale length h_R. The remaining parameters show significant correlations. This result might appear at odds with the observed RAR because the observed residuals do not seem to correlate with the galaxy properties in the SPARC sample (Lelli et al. 2017).

This issue requires additional clarification. In fact, unlike the SPARC sample, the DMS sample also shows some correlations: the p-values listed in Table 5 indicate that the RAR residuals are not significantly correlated (−log₁₀p <  3) only with h_R, I_d0, R_e, and the stellar surface brightness profile, Σ_*(R) (see also Figs. 14 and 15). Moreover, similarly to RG, correlations are found between the residuals of the QUMOND models and all the galaxy properties but the bulge central surface brightness I_e⁶. For QUMOND, this result might not be surprising, however, because QUMOND is a modified-gravity version of MOND, where a non-null intrinsic scatter for the RAR, and thus correlated residuals, are expected for non-spherical systems, unlike modified-inertia versions of MOND, that predict a null intrinsic scatter, and thus uncorrelated residuals (Bekenstein & Milgrom 1984).

The correlations we find for the residuals of the DMS data might partly generate the correlations we find for the RG residuals. In addition, the correlations for the DMS data, at odds with the claimed uncorrelations for the SPARC sample, might suggest a difference between the DMS and the SPARC samples. Unfortunately, we cannot quantify this difference here, if there even is any, because Lelli et al. (2017) only mention in their analysis that the Kendall and Spearman coefficients are in the range of [ − 0.2, 0.1] but they do not report their corresponding significance levels, namely their p-values. Therefore, we cannot assess the statistical significance of the lack of correlation. In fact, many coefficients in Table 5 from our analysis are in the range of [ − 0.2, 0.1], but their p-values clearly indicate that they are at many σ’s away from the null expected means and, thus, they demonstrate the presence of a statistically significant correlation.

The significance levels listed in Table 5 suggest that the correlations for the RG models are much stronger than for the DMS data for all the galaxy properties but I_d0. In addition, RG shows a very strong correlation, at largely more than 5σ, namely −log₁₀p >  6.24, with the radially-dependent properties of the galaxies, whereas the data show a significant correlation, between 4 and 5σ, namely 4.20 <  −log₁₀ p <  6.24, for R and f_gas(R), but no correlation for Σ_*(R).

Therefore, a possible serious tension between RG and the data might indeed be present. However, given the possible tension between the DMS and the SPARC samples, which is yet to be confirmed, we conclude that this issue remains open at this stage of our testing of RG. Further investigations with multiple data samples are necessary to clarify whether reproducing the observed properties of the RAR is indeed a challenge for RG.

7. Conclusions

We test the viability of RG, a theory of modified gravity that does not require the existence of dark matter to describe the dynamics of cosmic structures. We test RG on the scale of galaxies with 30 disk galaxies from the DMS (Bershady et al. 2010a). These galaxies appear almost face-on, having an inclination between 5 deg and 46 deg with respect to the line of sight. We can thus model, using a MCMC approach, both the rotation curves and the vertical velocity dispersion profiles.

The Poisson equation in RG contains a universal monotonic function of the local mass density, the gravitational permittivity, for which we adopt a smooth step function containing three free parameters: the permittivity of the vacuum ϵ₀, the critical density ρ_c, which sets the transition between the RG and the Newtonian regimes, and the transition power index Q.

By modelling each galaxy with two free parameters, the mass-to-light ratio, Υ, and the disk-scale height, h_z, and the additional three RG parameters, our MCMC analysis shows that RG is indeed able to describe the dynamics of these DMS galaxies with sensible values of Υ and h_z. Specifically, the mass-to-light ratios are consistent with those expected by SPS models (Bell & de Jong 2001); the agreement improves when we model both the rotation curves and the vertical velocity dispersion profiles rather than the rotation curves alone.

Similarly, the disk-scale heights, h_z, are consistent with the disk thicknesses, h_z, SR, inferred from the observation of edge-on galaxies. When modelling both the rotation curves and the vertical velocity dispersion profiles, h_z appears to prefer values that are a factor of ∼2 smaller than h_z, SR. However, this discrepancy might originate from an observational bias pointed out by Milgrom (2015) and quantified by Aniyan et al. (2016): the estimate of h_z, SR is based on near infrared photometry coming from old giant stars, with a larger vertical velocity dispersion and disk-scale height, whereas the vertical velocity dispersion profiles are estimated from integrated spectra, where the signal is dominated by young giants that have a smaller velocity dispersion. The observed velocity dispersion thus underestimates the velocity dispersion that would correspond to the observed h_z, SR. We estimate a bias factor of 1.63 for the DMS galaxies, a value that is consistent within 1σ with the value suggested by Aniyan et al. (2016).

In the formulation of RG, ϵ₀, Q, and ρ_c should be universal parameters, rather than free parameters for each individual galaxy as we assume above. To verify that a single set of these three parameters might indeed model the entire galaxy sample, we adopt a simple strategy. We assume that the mass-to-light ratio, Υ, and the disk-scale height, h_z, set by the previous analysis are appropriate and we only perform a MCMC exploration of the three-dimensional RG parameter space for the entire galaxy sample. The most likely set of the RG parameters is, within 1.2σ, consistent with the distributions of the RG parameters from the previous analysis. Although, as expected, the models of the rotation curves and of the vertical velocity dispersion profiles worsen, they still appear to be reasonable for most galaxies. This result suggests that a unique set of RG parameters describing the kinematics of the DMS galaxies with reasonable Υ and h_z might indeed exist.

Finally, we show that RG can in principle describe the RAR, namely the relation between the observed acceleration and the Newtonian acceleration due to the baryonic matter alone (McGaugh et al. 2016). The RAR expected in RG interpolates the data reasonably well and has the two asymptotic limits of the observed RAR. However, the models slightly underestimate the observed accelerations at low Newtonian accelerations. Moreover, the intrinsic scatter originating from the RG models appears not to be consistent with zero, as galaxy samples different from DMS seem to suggest. In fact, McGaugh et al. (2016) use the SPARC sample to show that the observed scatter of 0.13 dex in the RAR, obtained by adopting the same mass-to-light ratio for all the galaxies, is comparable to the scatter of 0.12 dex due to rotation curves, disk inclinations and galaxy distance uncertainties and to possible variations in mass-to-light ratios; this agreement leaves negligible room for intrinsic scatter.

An additional tension might be the correlation between the galaxy properties and the residuals of the RG models from the RAR described by Eq. (20), which might appear at odds with the uncorrelations claimed by Lelli et al. (2017) for the SPARC sample. In particular, for the three radially-dependent properties of the galaxies, radius, surface brightness and gas fraction, the correlations are statistically significant at largely more than 5σ. However, we also find significant correlations with radius and gas fraction, at more than 4σ, for the residuals of the DMS data. Further investigations are thus required to settle the issue. Modelling the SPARC data with RG, as we plan to do next, might shed light on whether the correlations we find here are, in fact, a feature of RG and, therefore represent a serious failure for RG or whether they are partly derived from the observed galaxy sample.

---

Acknowledgments

We sincerely thank the referee, whose detailed comments largely improved and clarified the presentation of our results, and especially our discussion of the RAR. We also thank Alistair Hodson, Michal Bílek, Luisa Ostorero, and Stefano Camera for inspiring discussions. We additionally thank Andrea Mignone for his crucial help with the Poisson solver, Emille Ishida for her help with the MCMC analysis with the software JAGS, and Marco Aldinucci for his help with the parallelisation of the code with OpenMP. We thank Compagnia di San Paolo (CSP) for funding the graduate-student fellowship of VC. We also acknowledge partial support from the INFN grant InDark and the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016.

References

Appendix A: Modelling the mass distribution in the disk galaxies

Here, we describe how we model (1) the surface brightness of the bulge and the stellar disk (Appendix A.1); (2) the mass distribution of the gas (Appendix A.2); and (3) the total three-dimensional mass distribution of the galaxy (Appendix A.3). We adopt the Hubble constant H₀ = 73 km s⁻¹ Mpc⁻¹, as per Martinsson et al. (2013a).

A.1. Surface brightness of the bulge and the stellar disk

The surface brightness of disk galaxies is the sum of the contribution of the bulge, which dominates the inner regions, and the contribution of the disk, which dominates the external regions. The surface brightness profile usually shows a distinct change of slope at the radius that separates the central region dominated by the bulge and the outer region dominated by the disk (see sub-panels (a) of Figs. D.1 and D.2).

To preserve the specific features of the distribution of the luminous matter, which in general correspond to features of the rotation curve, according to “Renzo’s rule” (Sancisi 2004), we model the disk surface brightness with a linear interpolation of the surface brightness data points only beyond the central region. To estimate the disk contribution within the central region, we fit the surface brightness data point beyond the central region with an exponential profile (de Vaucouleurs 1959; Freeman 1970):

$$\begin{array}{c}\hfill I_{\mathrm{d}}(R)=I_{\mathrm{d}0}exp(-\frac{R}{h_R}),\end{array}$$(A.1)

where h_R is the disk-scale length.

Therefore, to describe the disk surface brightness on the full radial range of the galaxy, we adopt the exponential model for the disk contribution within the central region, dominated by the bulge, and the linear interpolation of the data points in the outer region. We also adopt the exponential model in the very outer region if the grid of the Poisson solver we use in Sects. 3–6 goes beyond the measured profile.

We model the bulge surface brightness with a Sérsic spherical profile:

$$\begin{array}{c}\hfill I_{\mathrm{b}}(R)=I_{\mathrm{e}}exp\{-7.67[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{\frac{1}{n_{\mathrm{s}}}}-1]\},\end{array}$$(A.2)

with R_e the effective radius, I_e the surface brightness at radius R_e, and n_s the Sérsic index. The average bulge-to-total luminosity ratio in the K-band in the DMS galaxies is 0.09 (see Table 1 in Martinsson et al. 2013a); therefore, ignoring the triaxiality of the bulge should introduce negligible systematic errors on the modelling of the galaxy luminosity.

For the bulge, we adopt the Sérsic model rather than interpolating the surface brightness data points as we do for the disk, for two reasons: (1) the seeing affects the central region more than the outer region, as we detail below, and (2) unlike a two-dimensional disk, the bulge cannot be trivially deprojected without an analytical approximation. The deprojection is a required step to model the galaxy mass distribution we illustrate in Appendix A.3.

To model the surface brightness profile, we also consider the seeing that affects the measurements made with the 3.5 m diameter ground telescope at the Calar Alto Observatory. We model the seeing with a Gaussian point spread function,

$$\begin{array}{c}\hfill p_{\mathrm{PSF}}(R)=\frac{1}{2\pi {\sigma }^2}exp[-\frac{{(R-R_{\text{peak}})}^2}{2{\sigma }^2}],\end{array}$$(A.3)

where σ is the average effective seeing, as listed in Tables 3 and 4 of Martinsson et al. (2013a), and R_peak is the location of the surface brightness peak.

According to Martinsson et al. (2013a), we only convolve the bulge profile (Eq. (A.2)) with Eq. (A.3):

$$\begin{array}{c}\hfill I_{\mathrm{obs}}(R)={\displaystyle {\int }_{-\infty }^{+\infty }}{I}_{\mathrm{true}}(R^{\prime })p_{\mathrm{PSF}}(R^{\prime }-R)\mathrm{d}{R}^{\prime }.\end{array}$$(A.4)

Ignoring the effect of the seeing on the disk should introduce negligible systematic errors because the disk is a factor of ten more luminous than the bulge, on average, and spatially more extended.

The galaxy disks in the DMS sample appear almost face-on, with small inclination angles i. If we neglect the dust extinction within the galaxy, the observed surface brightness of the disk is brighter than the intrinsic surface brightness because the observed flux appears to come from an area cos i times smaller than the actual disk area. Therefore, to derive the light distribution of the disk from the measured surface brightness, μ_K, we need to correct for the inclination of the disk with respect to the line of sight.

By assuming the Tully–Fisher relation (Verheijen 2001),

$$\begin{array}{c}\hfill {log}_{10}{v}_{\mathrm{TF}}=-0.30103+(5.12-M_K)/11.3,\end{array}$$(A.5)

Martinsson et al. (2013a) estimate the inclination angle,

$$\begin{array}{c}\hfill sin{i}_{\mathrm{TF}}=\frac{v_{\mathrm{obs}}}{v_{\mathrm{TF}}},\end{array}$$(A.6)

from the observed luminosity, M_K, and the observed rotation velocity, v_obs. The face-on surface brightness of the disk is thus:

$$\begin{array}{c}\hfill {\mu }_K^i={\mu }_K-2.5C_{k}log(cos{i}_{\mathrm{TF}}),\end{array}$$(A.7)

where C_k = 1, for a transparent disk, as assumed by Martinsson et al. (2013a).

For the conversion of the surface brightness from the astronomical units mag arcsec⁻² to units L_⊙ pc⁻², we use the equation (Binney & Merrifield 1998):

$$\begin{array}{c}\hfill \mu \left(\frac{\text{mag}}{{\text{arcsec}}^2}\right)=M_{\odot ,K}+21.572-2.5{log}_{10}I\left(\frac{L_{\odot }}{{\text{pc}}^2}\right),\end{array}$$(A.8)

where M_⊙, K = 3.28 is the absolute magnitude of the Sun in the K-band. This equation neglects the expansion of the Universe because all the DMS galaxies are nearby, with the farthest galaxy (UGC 4622) at z = 0.043, namely at the distance 178.2 Mpc (see Martinsson et al. 2013a, Table 1).

We now illustrate the steps we adopt to model the entire surface brightness profile. We apply this procedure to the surface brightness profiles that were already corrected for inclination with Eq. (A.7) by Martinsson et al. (2013a).

As anticipated, most observed surface brightness profiles in the DMS sample show a distinct change of slope that suggests how extended the bulge is. We remove the central data points at radii smaller than the location of the slope change. The removed data points for each galaxy are shown in sub-panels (a) of Figs. D.1 and D.2 as green dots. We assume that the remaining profile, given by the red dots in sub-panels (a) of Figs. D.1 and D.2, is only due to the surface brightness of the disk. We linearly interpolate this remaining profile. The disk contribution in the most inner region, where the bulge contribution is dominant, is estimated with the exponential model of Eq. (A.1), with the parameters of the model set by the least-square best fit to the data points of the outer region alone.

We now subtract the extrapolated surface brightness profile of the disk from the observed total surface brightness profile of the central region, namely the green data points in sub-panels (a) of Figs. D.1 and D.2. The remaining profile is the surface brightness profile of the bulge that we now model with Eq. (A.4). Because of the presence of the convolution integral, to estimate the three parameters I_e, R_e and n_s of the bulge, it is more convenient to apply a MCMC approach. We run the MCMC algorithm with JAGS⁷, a free software that adopts the Gibbs sampling algorithm to generate the Markov chains.

We adopt Gaussian priors on the three free parameters; we set the Gaussian tails to zero for the unphysical negative values of the parameters. The Gaussian dispersions are set to $1/\sqrt{(200)}{L}_{\odot }$ pc⁻², $1/\sqrt{(200)}$ arcsec, and $1/\sqrt{(200)}$ for I_e, R_e and n_s, respectively. Assuming a Gaussian prior avoids the choice of an upper limit required in a uniform prior. To determine the means of the Gaussian priors, we compare the surface brightness profiles with the model profiles with a number of different choices of the parameters set by hand. We pick up the set of parameters that qualitatively best reproduce the data. This simple approach is sufficient to set the means of the Gaussian priors to reasonable values. The burn-in chain has 10 000 steps and we then run the chain for 100 000 steps.

For the best parameter estimates and their uncertainties, we adopt the medians and the standard deviations of the posterior distributions obtained from the MCMC runs. This choice is justified by the fact that the posterior distributions show a single peak and are basically symmetric.

By adopting this procedure, we find that two galaxies, UGC 1862 and UGC 9965, have no surface brightness in excess to the disk in the central region and we thus consider them bulgeless, a finding that is in agreement with Martinsson et al. (2013a). Martinsson et al. (2013a) find two additional bulgeless galaxies: UGC 3091 and UGC 7244; however, for these galaxies, we do find a non-negligible light excess in the central region.

Table A.1 lists the parameters estimated from the MCMC analysis for the bulge profiles and the best-fit parameters of the exponential profile used to estimate the disk surface brightness in the central region. The uncertainties on the parameters I_d0 and h_R are derived from the covariance matrix obtained from the least-square fit. Sub-panels (a) of Figs. D.1 and D.2 show the 30 measured profiles, corrected for inclination, in green the region where the bulge contribution is dominant, and in red the region where the disk contribution is dominant; the models are in blue and are the sum of Eq. (A.4) and the surface brightness of the disk.

To convert the fit parameters from angular (arcsec) to physical (kpc) radial units we use the relation:

$$\begin{array}{c}\hfill R(\text{kpc})=4.84814\times {10}^{-6}D(\text{kpc})R(\text{arcsec}),\end{array}$$(A.9)

where 4.84814 × 10⁻⁶ is the conversion factor from arcsec to radians and D(kpc) is the galaxy distance reported in Table 1 of Martinsson et al. (2013a). This relation is strictly valid in a non-expanding Euclidean geometry, but it can be applied to our DMS galaxies because they all are at redshift smaller than 0.043.

A.2. Gas surface mass density

The gas component of each galaxy is distributed in a disk-like structure that is thinner and more extended than the stellar disk. In addition, both the atomic and the molecular gas contribute to the gas component.

The observed atomic profile is set to Σ_atom = 1.4Σ_HI (Martinsson et al. 2013b), where Σ_HI is the measured HI gas surface mass density estimated from 21 cm radio synthesis observations (see Martinsson 2011, Sect. 2.5). Similarly to the surface brightness of the disk, we linearly interpolate the data points of the surface mass density of the atomic gas.

Martinsson (2011) measured Σ_HI only for 24 galaxies out of the original sample of 30 galaxies. For the remaining six (UGC 1081, UGC 1529, UGC 1862, UGC 1908, UGC 3091, UGC 12391), Martinsson (2011) modelled the Σ_HI profiles with a Gaussian with mean and dispersion taken from the average of the other galaxies with measured Σ_HI (see Sect. 3.1 of Martinsson et al. 2013b, for details), obtaining synthetic data. For these six galaxies, we clearly adopt their synthetic Gaussian profiles.

The molecular gas surface mass density, Σ_mol, is indirectly derived from 24 μm Spitzer observations, based on the CO line detection (see Martinsson et al. 2013b, Sect. 3.2). Again, we linearly interpolate the measures of the surface mass density of the molecular gas.

In sub-panels (b) and (c) of Figs. D.1 and D.2, the red dots with error bars are the measures of the surface mass density of the atomic and molecular gas; the blue lines show the linearly interpolated profiles.

A.3. Three-dimensional mass density model

We model the total baryonic mass density as the sum of the disk mass density, Υ_dj_d(R, z), the bulge mass density, Υ_bj_b(r), with r = (R² + z²)^1/2, and the atomic and molecular gas mass densities, ρ_atom(R, z) and ρ_mol(R, z):

$$\begin{array}{c}\hfill \rho (R,z)={\mathrm{{\rm Y}}}_{\mathrm{d}}{j}_{\mathrm{d}}(R,z)+{\mathrm{{\rm Y}}}_{\mathrm{b}}{j}_{\mathrm{b}}(r)+{\rho }_{\text{atom}}(R,z)+{\rho }_{\text{mol}}(R,z),\end{array}$$(A.10)

where Υ_d and Υ_b are the mass-to-light ratios of the stellar populations of the disk and the bulge, respectively. We assume that the mass-to-light ratios are independent of R and z. In addition, we set Υ_b = Υ_d = Υ because the bulge is, on average, an order of magnitude less luminous than the disk, as mentioned in Appendix A.1. Equation (A.10) is the total mass distribution of the galaxy in modified gravity models where the dark matter component is assumed to be absent.

We model the three-dimensional luminosity density of the disk by multiplying the disk surface brightness profile, I_d(R), which is the sum of the exponential model and the interpolated profile, by an exponentially decreasing density profile along the vertical axis z,

$$\begin{array}{c}\hfill j_{\mathrm{d}}(R,z)=\frac{I_{\mathrm{d}}(R)}{2h_z}exp(-\frac{|z|}{h_z}),\end{array}$$(A.11)

where the disk-scale height h_z is a free parameter. The factor 1/(2h_z) provides the correct normalisation. In the main body of the paper, we compare our estimate of h_z with the scale height h_z, SR obtained from the relation derived in Bershady et al. (2010b) from a combined sample of 60 edge-on late-type galaxies,

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{h_R}{h_{z,\text{SR}}}\right)=0.367{log}_{10}\left(\frac{h_R}{\text{kpc}}\right)+0.708\pm 0.095,\end{array}$$(A.12)

where the term ±0.095 quantifies the ∼25% intrinsic scatter. We estimate the uncertainty on h_z, SR as the sum in quadrature of the intrisic scatter of the relation (A.12) and the uncertainty on h_R. The latter contribution is negligible with respect to the former, so the error on h_z, SR nearly coincides with its intrinsic scatter.

We model the three-dimensional luminosity density of the bulge with the Abel integral:

$$\begin{array}{c}\hfill j_{\mathrm{b}}(r)=-\frac{1}{\pi }{\displaystyle {\int }_r^{\infty }}\frac{\mathrm{d}{I}_{\mathrm{b}}}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}},\end{array}$$(A.13)

where I_b(R) is the surface brightness profile modelled as described in Appendix A.1, R is the radius projected on the sky and r is the three-dimensional radius. The equation above assumes a spherically symmetric bulge. As mentioned in Appendix A.1, neglecting the triaxial structure of the bulge should introduce negligible systematic errors in the mass estimates because the galaxy luminosity is dominated by the disk (Angus et al. 2015). As anticipated, adopting an analytical model for the two-dimensional surface brightness of the bulge, rather than interpolating the data points as we do for the disk, facilitates its deprojection into three dimensions.

We consider both the atomic and molecular gas distributions as razor-thin disks (Martinsson et al. 2013b),

$$\begin{array}{c}\hfill {\rho }_{\mathrm{atom},\mathrm{mol}}(R,z)={\mathrm{\Sigma }}_{\mathrm{atom},\mathrm{mol}}(R)\delta (z),\end{array}$$(A.14)

where Σ(R) is the linearly interpolated mass surface density of the gas disk and δ(z) is the Dirac δ function.

Appendix B: Successive over relaxation Poisson solver

B.1. Numerical solution of the Poisson equation

For the theories of gravity we consider here, deriving the gravitational potential ϕ that originates from the mass density distribution ρ requires solving the Poisson equation,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\varphi (R,z)+S(\rho ;R,z)=0,\end{array}$$(B.1)

where the source term S is a generic function of the density ρ. For axisymmetric disk galaxies, we limit the equation to cylindrical coordinates R and z. We only consider static models and Eq. (B.1) is thus an elliptic partial differential equation independent of time.

We solve the Poisson equation with a successive over relaxation (SOR) algorithm, an iterative procedure based on the Jacobi and the Gauss–Seidel algorithms (Young 1954). We find the solution on a rectangular grid of size L_R × L_z, with (N_R + 1)×(N_z + 1) grid points and steps Δ_R = L_R/N_R and Δ_z = L_z/N_z in the two dimensions (Table B.1).

Given the solution of the gravitational potential ${\varphi }_{i,k}^n$ at the nth iteration on the grid point of indexes (i, k), the solution at the following (n + 1)th iteration is:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\varphi }_{i,k}^{n+1}=& {\varphi }_{i,k}^n(1-{\omega }_{\text{SOR}})\hfill \\ \hfill & +\frac{{\omega }_{\text{SOR}}}{2R_i({\mathrm{\Delta }}_R^2+{\mathrm{\Delta }}_z^2)}[{\varphi }_{i+1,k}^{n+1}(R_i+\frac{1}{2}{\mathrm{\Delta }}_R){\mathrm{\Delta }}_z^2\hfill \\ \hfill & +{\varphi }_{i-1,k}^{n+1}(R_i-\frac{1}{2}{\mathrm{\Delta }}_R){\mathrm{\Delta }}_z^2+({\varphi }_{i,k+1}^{n+1}+{\varphi }_{i,k-1}^{n+1})R_i{\mathrm{\Delta }}_R^2\hfill \\ \hfill & +R_i{S}_{i,k}{\mathrm{\Delta }}_R^2{\mathrm{\Delta }}_z^2],\hfill \end{array}\end{array}$$(B.2)

where S_i, k is the source term on the grid and ω_SOR is a parameter in the range of (0, 2) to guarantee the convergence. The value of ω_SOR which guarantees the fastest convergence is ω_SOR = 2/(1 + π/N) for a square grid and ω_SOR = 2/(1 + π/N_min) for a rectangular grid, where N_min is the smallest number between N_R and N_z. In general, it is computationally convenient to choose the number of grid points N such that ω_SOR is in the range of (1, 2); with this choice, the number of iterations necessary to reach convergence is linearly proportional to N, whereas for ω_SOR in the range of (0, 1) the number of iterations is proportional to N² (Young 1954).

We adopt L_R = 2 × 12h_R for almost all galaxies and L_z = 2 × 100h_z, SR for all galaxies, where h_z, SR is derived from Eq. (A.12). With this choice, the grid domain is substantially larger than the galaxy size, and we can adopt the asymptotic behaviour of the gravitational potential to set the proper boundary conditions, as we describe in Appendix B.3 below. For UGC 4368 and UGC 6918, we adopt L_R = 2 × 20h_R and L_R = 2 × 18h_R, respectively, because L_R = 2 × 12h_R is smaller than the extension of the measured rotation curve. UGC 4458 has h_R = 9.27 kpc and, with this large scale length, we adopt L_R = 2 × 10h_R, which is already sufficient to reach the asymptotic behaviour of the gravitational potential.

The radial resolution of the rotation curves and of the vertical velocity dispersion profiles measured for each galaxy in the DMS sample differs from galaxy to galaxy. For each galaxy, we thus choose N_R and N_z that yield both ω_SOR in the range of (1, 2) and the numerical resolution Δ_R = L_R/N_R and Δ_z = L_z/N_z comparable to the observed resolution.

For most galaxies, the grid where we compute the mass distribution and the galactic potential is more extended than the measured surface brightness of the disk and the measured gas surface mass density. To estimate the mass distribution in these regions, we extrapolate the disk surface brightness with Eq. (A.1) with the parameters listed in Table A.1; we instead set to zero the gas mass density because its contribution to the galaxy mass is negligible in these outer regions.

We centre the computational domain on the origin R = z = 0. The coordinate R appears at the denominator in Eq. (B.2). Therefore, to avoid divergences, we choose the grid so that the R = 0 axis is not a grid strand, unlike the z = 0 axis. We can thus compute both the rotation curve in the plane z = 0 and the vertical velocity dispersion at z = 0 for any R ≠ 0.

We set the initial values of the potential to ${\varphi }_{i,k}^0=0$ over the entire domain except the boundaries (see Appendix B.3) and stop the iteration when:

$$\begin{array}{c}\hfill {\epsilon }^{n+1}=\frac{1}{(N_R-1)(N_z-1)}{\displaystyle \sum _{i,k}}|{\varphi }_{i,k}^{n+1}-{\varphi }_{i,k}^n|<{10}^{-9}.\end{array}$$(B.3)

We test our algorithm with mass density distributions where the Poisson equation in Newtonian gravity can be solved analytically: the Miyamoto–Nagai disk (see Eqs. (B.12) and (B.13)), Satoh disk, logarithmic potential, Plummer sphere, isochrone potential, Hernquist sphere, and Navarro–Frenk–White potential (Binney & Tremaine 2008). We compare the numerical and the analytical potentials as a function of R and z.

Within the half-scale length from the centre, the numerical solution is within 1% of the analytic solution for the Miyamoto–Nagai and the Satoh disks, 0.05% for the logarithmic potential, 0.25% for the Plummer sphere, 0.13% for the isochrone potential, 4% for the Hernquist sphere, and 2% for the Navarro–Frenk–White potential. The agreement improves outwards: beyond two scale lengths, it is smaller than 0.5% for the Miyamoto–Nagai and Satoh disks, Hernquist and Navarro–Frenk–White spheres, 0.1% for the Plummer sphere, 0.05% for the isochrone potential, and 0.015% for the logarithmic potential).

B.2. The source term S(R,z)

In this work, we consider two gravity theories: MOND and RG.

In the QUMOND formulation of MOND (Milgrom 2010), we have:

$$\begin{array}{c}\hfill S_{\text{MOND}}(R,z)=-\mathrm{\nabla }·\left[\nu \left(\frac{|\mathrm{\nabla }{\varphi }_{\mathrm{N}}|}{a_0}\right)\mathrm{\nabla }{\varphi }_{\mathrm{N}}\right],\end{array}$$(B.4)

where ϕ_N is the Newtonian gravitational potential due to the mass density distribution of the baryonic matter, and ν(y) is the MOND interpolating function, with y = |∇ϕ_N|/a₀. Here, we adopt the simple ν-function, given by Eq. (17) (see Famaey & McGaugh 2012, Eq. (50) with n = 1).

In cylindrical coordinates, Eq. (B.4) becomes:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S_{\text{MOND}}(R,z)=& -(\frac{\nu }{R}\frac{\partial {\varphi }_{\mathrm{N}}}{\partial R}+\frac{\partial \nu }{\partial R}\frac{\partial {\varphi }_{\mathrm{N}}}{\partial R}+\nu \frac{{\partial }^2{\varphi }_{\mathrm{N}}}{\partial R^2}\hfill \\ \hfill & +\frac{\partial \nu }{\partial z}\frac{\partial {\varphi }_{\mathrm{N}}}{\partial z}+\nu \frac{{\partial }^2{\varphi }_{\mathrm{N}}}{\partial z^2})·\hfill \end{array}\end{array}$$(B.5)

Deriving the MOND gravitational potential clearly requires to solve the Poisson equation twice: first, to estimate the Newtonian potential ϕ_N, where, in Eq. (B.2), we use the standard source term:

$$\begin{array}{c}\hfill S_{\mathrm{N}}(R,z)=-4\pi G\rho (R,z),\end{array}$$(B.6)

and, subsequently, to compute ϕ_MOND with the source term S_MOND(R, z) of Eq. (B.5).

To derive the source term in the RG case, we need to recast the RG Poisson equation (1),

$$\begin{array}{c}\hfill \mathrm{\nabla }·[ϵ(\rho )\mathrm{\nabla }{\varphi }_{\text{RG}}]=4\pi G\rho ,\end{array}$$(B.7)

as:

$$\begin{array}{c}\hfill \mathrm{\nabla }ϵ(\rho )·\mathrm{\nabla }{\varphi }_{\text{RG}}+ϵ(\rho ){\mathrm{\nabla }}^2{\varphi }_{\text{RG}}=4\pi G\rho ,\end{array}$$(B.8)

so that the source term is:

$$\begin{array}{c}\hfill S_{\mathrm{RG}}(R,z)=-\frac{4\pi G\rho (R,z)-\mathrm{\nabla }ϵ(\rho )·\mathrm{\nabla }{\varphi }_{\text{RG}}(R,z)}{ϵ(\rho )}·\end{array}$$(B.9)

Here, the source term contains the unknown ϕ_RG. At each iteration step, in the source term, we insert the potential ϕ_RG estimated at the previous step.

The form of Eq. (B.9) requires some additional care in the numerical algorithm: the source term increases when the vacuum permittivity ϵ₀ decreases, and the Poisson solver does not necessarily converge with the optimal value ω_SOR = 2/(1 + π/N) ∼ 1.97−1.98, for our typical N_R and N_z. To make the Poisson solver converge for ϵ₀ in the flat prior range of [0.10−1], we need to set ω_SOR to values smaller than ∼1.97−1.98. For example ω_SOR = 2/(1 + π/25) guarantees convergence for every ϵ₀ in our flat prior range. Yet, setting ω_SOR to this unique value would increase the total computational time by at least a factor of 4. We thus vary ω_SOR from 2/(1 + π/25), for ϵ₀ close to 0.10, to 2/(1 + π/200), for ϵ₀ close to 1, according to the explored value of ϵ₀ within the flat prior range.

In MOND, this issue is not present and we can use the same value of ω_SOR for all the parameter combinations, namely 2/(1 + π/N) for a square grid or 2/(1 + π/N_min) for a rectangular grid.

Numerically, we compute all the derivatives present in all the above equations with the leap-frog method. The first and second derivatives of a function f at a point x are

$$\begin{array}{c}\hfill \frac{\partial f}{\partial x}=\frac{f(x+\mathrm{\Delta }x)-f(x-\mathrm{\Delta }x)}{2\mathrm{\Delta }x}\end{array}$$(B.10)

and

$$\begin{array}{c}\hfill \frac{{\partial }^{2}f}{\partial x^2}=\frac{f(x+\mathrm{\Delta }x)+f(x-\mathrm{\Delta }x)-2f(x)}{\mathrm{\Delta }{x}^2},\end{array}$$(B.11)

where Δx is the grid step. With the grid steps we adopt, the leap-frog method guarantees relative errors smaller than 5% on the estimated derivatives.

B.3. Boundary conditions

Elliptic equations, like our Poisson equations, require boundary conditions to be set. To facilitate the choice of the boundary conditions we put the disk galaxy at the centre of the grid domain and choose the size of the domain substantially larger than the galaxy size. We could use the axial symmetry of the problem and only consider a fourth of that domain, with two boundary conditions along the R and z axes. This choice would clearly reduce the computation time by a factor of four, but it would pose the non-trivial problem of setting the boundary conditions on the two axes cutting through the galaxy centre. We thus prefer to compute the potential in all the four quadrants and to set the boundary conditions far from the galaxy centre.

Since the bulge and gas components are not dynamically dominant in the DMS galaxies (Bershady et al. 2010a), we use only the disk component to set the boundary conditions. In the external regions of the domain we model this component with a double exponential disk (Eq. (A.11)), since we use the exponential profile (A.1) to extrapolate the disk surface brightness in these regions. In Newtonian gravity, this mass distribution generates a gravitational potential that does not have a close analytic expression, but it is well approximated by the sum of the gravitational potentials of three Miyamoto–Nagai disks (Smith et al. 2015). The Miyamoto–Nagai disk has mass density distribution (Miyamoto & Nagai 1975):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{MN}}(R,z)=\left(\frac{b^2{M}_{\text{MN}}}{4\pi }\right)\frac{a{R}^2+(a+3\sqrt{z^2+b^2}){(a+\sqrt{z^2+b^2})}^2}{{[R^2+{(a+\sqrt{z^2+b^2})}^2]}^{5/2}{(z^2+b^2)}^{3/2}},\end{array}$$(B.12)

and generates the Newtonian gravitational potential (e.g. Binney & Tremaine 2008),

$$\begin{array}{c}\hfill {\varphi }_{\text{MN}}(R,z)=-\frac{G{M}_{\text{MN}}}{\sqrt{R^2+{(a+\sqrt{z^2+b^2})}^2}},\end{array}$$(B.13)

where M_MN is the disk mass and a and b are its scale length and scale height.

The gravitational potential of a double exponential disk can thus be approximated by the equation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\varphi }_{\text{exp}}(R,z)=& -\frac{G{M}_1}{\sqrt{R^2+{(a_1+\sqrt{z^2+b^2})}^2}}\hfill \\ \hfill & -\frac{G{M}_2}{\sqrt{R^2+{(a_2+\sqrt{z^2+b^2})}^2}}\hfill \\ \hfill & -\frac{G{M}_3}{\sqrt{R^2+{(a_3+\sqrt{z^2+b^2})}^2}},\hfill \\ \hfill \end{array}\end{array}$$(B.14)

where M₁, M₂ and M₃ are the three Miyamoto–Nagai disk masses, a₁, a₂, a₃ are their scale lengths, and b is the scale height which is equal in all the three disks. The three disks are only mathematical artefacts to generate the physical gravitational potential, but they do not correspond to three physical disks. Indeed, the mass M₂ of the second disk is always negative.

The parameters of the Miyamoto–Nagai disks are related to the parameters of the double exponential disk (Eq. (A.11)) and its total mass M_d according to the equations (Smith et al. 2015)

$$\begin{array}{cc}\hfill \frac{b}{h_R}=& -0.269{\left(\frac{h_z}{h_R}\right)}^3+1.080{\left(\frac{h_z}{h_R}\right)}^2+1.092\frac{h_z}{h_R},\hfill \end{array}$$(B.15)

$$\begin{array}{cc}\hfill \frac{M_1}{M_{\mathrm{d}}}=& -0.0090{\left(\frac{b}{h_R}\right)}^4+0.0640{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & -0.1653{\left(\frac{b}{h_R}\right)}^2+0.1164\frac{b}{h_R}+1.9487,\hfill \end{array}$$(B.16)

$$\begin{array}{cc}\hfill \frac{M_2}{M_{\mathrm{d}}}=& 0.0173{\left(\frac{b}{h_R}\right)}^4-0.0903{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & +0.0877{\left(\frac{b}{h_R}\right)}^2+0.2029\frac{b}{h_R}-1.3077,\hfill \end{array}$$(B.17)

$$\begin{array}{cc}\hfill \frac{M_3}{M_{\mathrm{d}}}=& -0.0051{\left(\frac{b}{h_R}\right)}^4+0.0287{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & -0.0361{\left(\frac{b}{h_R}\right)}^2-0.0544\frac{b}{h_R}+0.2242,\hfill \end{array}$$(B.18)

$$\begin{array}{cc}\hfill \frac{a_1}{h_R}=& -0.0358{\left(\frac{b}{h_R}\right)}^4+0.2610{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & -0.6987{\left(\frac{b}{h_R}\right)}^2-0.1193\frac{b}{h_R}+2.0074,\hfill \end{array}$$(B.19)

$$\begin{array}{cc}\hfill \frac{a_2}{h_R}=& -0.0830{\left(\frac{b}{h_R}\right)}^4+0.4992{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & -0.7967{\left(\frac{b}{h_R}\right)}^2-1.2966\frac{b}{h_R}+4.4441,\hfill \end{array}$$(B.20)

$$\begin{array}{cc}\hfill \frac{a_3}{h_R}=& -0.0247{\left(\frac{b}{h_R}\right)}^4+0.1718{\left(\frac{b}{h_R}\right)}^3\hfill \\ \hfill & -0.4124{\left(\frac{b}{h_R}\right)}^2-0.5944\frac{b}{h_R}+0.7333.\hfill \end{array}$$(B.21)

In Newtonian gravity, we impose that the gravitational potential equals Eq. (B.14) on the borders of the rectangular domain. In MOND, when the acceleration drops below the critical acceleration, a₀ = 3600 kpc⁻¹(km s⁻¹)², the gravitational field has the asymptotic behaviour (Milgrom 1983a):

$$\begin{array}{c}\hfill g=\sqrt{a_0|g_{\mathrm{N}}|},\end{array}$$(B.22)

which, in cylindrical components, becomes:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial R}=\sqrt{a_0\left|\frac{\partial {\varphi }_{\mathrm{N}}}{\partial R}\right|},\end{array}$$(B.23)

and

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial z}=\sqrt{a_0\left|\frac{\partial {\varphi }_{\mathrm{N}}}{\partial z}\right|}·\end{array}$$(B.24)

To solve the MOND Poisson equation, we thus impose the Neumann boundary conditions (B.23) and (B.24), with ∂ϕ_N/∂R and ∂ϕ_N/∂z derived analytically from Eq. (B.14).

As illustrated in Matsakos & Diaferio (2016), at sufficiently large distances from the galaxy centre, RG recovers the MOND asymptotic behaviour of the radial gravitational field. We thus impose the Neumann condition (B.23) on the four borders of the domain for the radial component of the gravitational field. For the vertical component, we resort to the fact that the field lines are refracted and, at large distances from the source, they are almost parallel to the disk plane. We thus set:

$$\begin{array}{c}\hfill \frac{\partial \varphi }{\partial z}=0\end{array}$$(B.25)

on the domain borders.

In the analysis of the DMS galaxies, both the disk-scale height b and the disk mass M_d, which appear in Eqs. (B.14)–(B.21), are initially unknown. These parameters depend on h_z and Υ, which are two free parameters of the fit. Therefore, we update the boundary conditions derived from Eqs. (B.23) and (B.24) at each step of the MCMC chain.

Appendix C: Convergence tests of the MCMC chains

C.1. The variance ratio method

The variance ratio method of Gelman & Rubin (1992) is a convergence diagnostic test for monitoring the convergence of MCMC chains: in other words, it estimates how close to convergence a chain is and if the convergence can be improved with additional steps of the chain (Brooks & Roberts 1998).

The test compares the output of m ≥ 2 independent chains that have different starting points. Each chain has n₀ + n elements, where the first n₀ are the discarded burn-in chain. Each chain i, at each step t, returns an estimate of a given parameter of interest ${\theta }_i({\mathbf{x}}^t)\equiv {\theta }_i^t$, where x are the variables of the problem updated at each chain step t. For each chain, we compute the mean of the parameter estimates,

$$\begin{array}{c}\hfill {\overline{\theta }}_i=\frac{1}{n}{\displaystyle \sum _{t=n_0+1}^{n_0+n}}{\theta }_i^t,\end{array}$$(C.1)

the mean of the means of the m chains,

$$\begin{array}{c}\hfill \overline{\theta }=\frac{1}{m}{\displaystyle \sum _{i=1}^m}{\overline{\theta }}_i,\end{array}$$(C.2)

and their variance,

$$\begin{array}{c}\hfill \frac{B}{n}=\frac{1}{m-1}{\displaystyle \sum _{i=1}^m}{({\overline{\theta }}_i-\overline{\theta })}^2.\end{array}$$(C.3)

In addition, we compute the variance of ${\theta }_i^t$ within each chain,

$$\begin{array}{c}\hfill s_i^2=\frac{1}{n-1}{\displaystyle \sum _{t=n_0+1}^{n_0+n}}{({\theta }_i^t-{\overline{\theta }}_i)}^2,\end{array}$$(C.4)

and their mean,

$$\begin{array}{c}\hfill W=\frac{1}{m}{\displaystyle \sum _{i=1}^m}{s}_i^2.\end{array}$$(C.5)

Gelman & Rubin (1992) assume that the mean of the posterior distribution of the parameter θ is $\widehat{\mu }=\overline{\theta }$ and that its variance is:

$$\begin{array}{c}\hfill {\widehat{\sigma }}^2=\frac{n-1}{n}W+\frac{B}{n}·\end{array}$$(C.6)

Gelman & Rubin (1992) model the variability of $\widehat{\mu }$ and ${\widehat{\sigma }}^2$ due to sampling with an approximate Student’s t distribution with mean $\widehat{\mu }=\overline{\theta }$, variance

$$\begin{array}{c}\hfill \widehat{V}={\widehat{\sigma }}^2+\frac{1}{m}\frac{B}{n},\end{array}$$(C.7)

and number of degrees of freedom $\nu =2{\widehat{V}}^2/\widehat{\text{var}}(\widehat{V})$, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{\text{var}}(\widehat{V})=& {\left(\frac{n-1}{n}\right)}^2\frac{1}{m}\widehat{\text{var}}(s_i^2)+{\left(\frac{m+1}{\mathit{mn}}\right)}^2\frac{2}{m-1}{B}^2\hfill \\ \hfill & +2\frac{(m+1)(n-1)}{m{n}^2}\frac{n}{m}[\widehat{\text{cov}}(s_i^2,{\overline{\theta }}_i^2)-2\overline{\theta }\widehat{\text{cov}}(s_i^2,{\overline{\theta }}_i)],\hfill \end{array}\end{array}$$(C.8)

and where the variance $\widehat{\text{var}}(s_i^2)$ and the covariances $\widehat{\text{cov}}(s_i^2,{\overline{\theta }}_i^2)$ and $\widehat{\text{cov}}(s_i^2,{\overline{\theta }}_i)$ are estimated with the m values ${\overline{\theta }}_i$ and $s_i^2$.

To monitor the convergence of the chains, Gelman & Rubin (1992) compute the potential scale reduction factor

$$\begin{array}{c}\hfill \widehat{R}=\frac{\widehat{V}}{W}\frac{\nu }{\nu -2}·\end{array}$$(C.9)

If $\widehat{R}\to 1$ as n → ∞, the estimated variance $\widehat{V}$ of the expected posterior distribution of θ is increasingly closer to the variance W estimated from the chains for an increasingly large number of steps. In other words, values of $\widehat{R}\gg 1$ suggest that the estimate of the target distribution can be improved with additional steps, whereas values of $\widehat{R}\sim 1$ suggest that the chains are close to the target distribution.

C.2. The MCMC convergence for the universal combination of the RG parameters

In our MCMC analysis we adopt a burn-in chain with n₀ = 1000 and a chain with n = 19 000 iterations. To assess the convergence of the chains, in Sect. 5 we adopt the variance ratio method.

We run the test for the first n = 13 000 iterations. If the test is positive for this n, we are confident that running the chains for n = 19 000 iterations will provide a posterior distribution close to the target distribution. We consider m = 3 chains.

The values we obtain for the potential scale reduction factor, $\widehat{R}$, for each RG parameter, ϵ₀, Q and ρ_c, are 1.01, 1.04 and 1.004, respectively. According to the interpretation of Gelman & Rubin (1992), these values suggest that our distributions estimated with n = 13000 iterations already are reasonably close to the target distributions. The results we show in the text are for n = 19 000 iterations and we are thus confident that our values for ϵ₀, Q and ρ_c are robust estimates of the target values.

Appendix D: Figures of the individual DMS galaxies

In this appendix, we collect all the figures showing the relevant quantities of each DMS galaxy. Figures D.1 and D.2 show the profiles of the surface brightness, the surface mass density of the atomic and molecular gas and the kinematic profiles. Each galaxy appears in an individual panel. Figure D.2 shows the five galaxies that we analysed both in RG and in QUMOND. Each panel contains 8 sub-panels in Fig. D.1 and 12 sub-panels in Fig. D.2.

The sub-panels (a) show the surface brightness. The filled circles with error bars are the data in the K-band, corrected for inclination (Martinsson et al. 2013a) and the blue solid curves are our models (Appendix A.1). The green data points are removed before performing the fit with the exponential profile (A.1). We only use this exponential model to estimate the disk surface brightness both in the central region, to separate the disk and the bulge contributions, and in the outer regions not covered by the data but still within the numerical domain of our Poisson solver. The grey vertical lines show the radius for our bulge-disk decomposition.

The sub-panels (b) show the surface mass density profile of the atomic gas. The filled circles with error bars are the estimates according to Martinsson et al. (2013b) and the blue solid curves are our linear interpolations (Appendix A.2). The sub-panels (c) show the surface mass density profile of the molecular gas (Appendix A.2). Symbols and lines are as in the sub-panels (b).

The sub-panels (d) show the rotation curves. The filled circles with error bars are the data; the blue solid lines show the RG models whose parameters are the medians of their posterior distributions estimated from the rotation curve alone (Sect. 3). The dashed magenta vertical lines show the bulge effective radius, whereas the dashed green vertical lines show the bulge radius we adopt in our disk-bulge decomposition for the surface brightness fit (Appendix A.1). For the galaxy UGC 1087, these two radii coincide. The dashed green vertical lines coincide with the grey lines of sub-panels (a). These two vertical lines are repeated in all sub-panels (d)–(h).

The sub-panels (e) again show the rotation curves and the sub-panels (f) show the vertical velocity dispersion profiles. In the sub-panels (e), the data are the same as in the sub-panels (d). In the sub-panels (f), the data are the filled circles with error bars. In the sub-panels (e) and (f), the blue solid lines show the RG models whose parameters are the medians of their posterior distributions estimated from the rotation curve and the vertical velocity dispersion profile at the same time (Sect. 4). In the sub-panels (f), the cyan solid lines show the vertical velocity dispersion profile when we adopt the same parameters as for the blue lines except for the disk-scale height h_z; for the cyan solid lines, h_z is the value derived from Eq. (A.12). For the galaxy UGC 6918 (Fig. D.1), the cyan line overlaps the blue line because the estimated h_z coincides with the value derived from Eq. (A.12).

The sub-panels (g) and (h) again show the rotation curves and the vertical velocity dispersion profiles. In these sub-panels, the blue solid lines show the RG models where the mass-to-light ratios and the disk-scale heights are set to the values derived in Sect. 4 and listed in Table 2, whereas the values of the three RG parameters ϵ₀, Q, and ρ_c are set to those of the unique combination found in Sect. 5.

Figure D.2 shows four additional sub-panels. The sub-panels (i) and (k) show the rotation curves and the sub-panels (j) and (l) show the velocity dispersion profiles. The measured velocity dispersion profiles in these sub-panels are artificially increased by the factor 1.63, which is our estimate of the observational bias suggested by Aniyan et al. (2016) (Sect. 4.1). The blue curves are the models whose parameters are the medians of their posterior distributions estimated from the rotation curve and the vertical velocity dispersion profile at the same time (Sect. 4.1). The sub-panels (i) and (j) show the RG models and the sub-panels (k) and (l) show the QUMOND models. The cyan lines in the sub-panels (j) and (l) show the models where the scale height h_z is set by Eq. (A.12). For the galaxy UGC 4107 (Fig. D.2), the cyan lines in sub-panels (j) and (l) overlap the blue lines because the estimated h_z in both RG and QUMOND are almost identical to the values obtained with Eq. (A.12). The vertical magenta and green lines are the same as in the sub-panels (d)–(h).

All Tables

All Figures