The Intrinsic Scatter of the Radial Acceleration Relation

TL;DR

This paper develops a Monte Carlo model to accurately estimate the intrinsic scatter in the stellar radial acceleration relation of galaxies, finding it slightly exceeds LCDM predictions, thus informing galaxy formation theories.

Contribution

It introduces a new Monte Carlo method to disentangle observational errors from intrinsic scatter in galaxy scaling relations.

Findings

Intrinsic scatter of the stellar RAR is 0.11±0.02 dex.

Observed scatter of the stellar RAR is 0.17 dex.

Results are consistent with LCDM predictions, slightly exceeding them.

Abstract

We present a detailed Monte Carlo model of observational errors in observed galaxy scaling relations to recover the intrinsic (cosmic) scatter driven by galaxy formation and evolution processes. We apply our method to the stellar radial acceleration relation (RAR) which compares the local observed radial acceleration to the local Newtonian radial acceleration computed from the stellar mass distribution. The stellar and baryonic RAR are known to exhibit similar scatter. Lelli+2017 (L17) studied the baryonic RAR using a sample of 153 spiral galaxies and inferred a negligible intrinsic scatter. If true, a small scatter might challenge the LCDM galaxy formation paradigm, possibly favoring a modified Newtonian dynamics interpretation. The intrinsic scatter of the baryonic RAR is predicted by modern LCDM simulations to be ~0.06-0.08 dex, contrasting with the null value reported by L17. We have assembled a catalog of structural properties with over 2500 spiral galaxies from six deep imaging and spectroscopic surveys (called PROBES for the "Photometry and Rotation curve OBservations from Extragalactic Surveys") to quantify the intrinsic scatter of the stellar RAR and other scaling relations. The stellar RAR for our full sample has a median observed scatter of 0.17 dex. We use our Monte Carlo method, which accounts for all major sources of measurement uncertainty, to infer a contribution of 0.12 dex from the observational errors. The intrinsic scatter of the stellar RAR is thus estimated to be 0.11$\pm$0.02 dex, in agreement with, though slightly greater than, current LCDM predictions.

Figures

Tables

Table 1. Table 1: Parametrization of the RAR Fit to Each PROBES Data set

1	2	3	4
Survey	$g_{†}$	${\sigma }_{\mathrm{RAR}}$	$N$
	(${10}^{-10}m{s}^{-1}$)	(dex)	#
SV	$0.813\pm 0.044$	$0.166\pm 0.007$	1380
SP	$0.839\pm 0.020$	$0.131\pm 0.003$	2561
SF	$0.558\pm 0.008$	$0.139\pm 0.001$	16,648
C97	$5.224\pm 0.042$	$0.182\pm 0.001$	17,626
M92	$1.790\pm 0.016$	$0.173\pm 0.001$	18,355
M96	$1.345\pm 0.012$	$0.170\pm 0.001$	21,460

Table 2. Table 2: First-order Scatter Predictions for all PROBES Surveys, with Corresponding Intrinsic Scatters

1	2	3	4
Survey	${\sigma }_{\mathrm{RAR}}$	${\sigma }_{1^{\mathrm{st}}\mathrm{Order}}$	${\sigma }_{\mathrm{int}}$
	dex	dex	dex
SV	0.166	0.114	0.121
SP	0.130	0.097	0.087
SF	0.139	0.089	0.107
C97	0.182	0.087	0.160
M92	0.173	0.092	0.146
M96	0.170	0.090	0.145

Table 3. Table 3: Parameterization for the RAR Fit to Each Mock Data Set.

1	2	3	4
Survey	$g_{†}$	${\sigma }_{\mathrm{sim}}$	$N$
	${10}^{-10}m{s}^{-1}$	dex	#
SV	$0.768\pm 0.020$	$0.134\pm 0.004$	1685
SP	$1.051\pm 0.021$	$0.125\pm 0.003$	2416
SF	$0.611\pm 0.004$	$0.117\pm 0.001$	17,173
C97	$5.626\pm 0.024$	$0.100\pm 0.001$	18,065
M92	$1.810\pm 0.009$	$0.130\pm 0.001$	18,401
M96	$1.340\pm 0.007$	$0.121\pm 0.001$	21,137

Table 4. Table 4: Scatter measurements for the RAR.

1	2	3	4
Survey	${\sigma }_{\mathrm{RAR}}$	${\sigma }_{\mathrm{sim}}$	${\sigma }_{\mathrm{int}}$
	(dex)	(dex)	(dex)
SV	$0.166\pm 0.007$	$0.134\pm 0.004$	$0.098\pm 0.013$
SP	$0.131\pm 0.003$	$0.125\pm 0.003$	$0.040\pm 0.013$
SF	$0.139\pm 0.001$	$0.117\pm 0.001$	$0.076\pm 0.002$
C97	$0.182\pm 0.001$	$0.100\pm 0.001$	$0.153\pm 0.002$
M92	$0.173\pm 0.001$	$0.130\pm 0.001$	$0.114\pm 0.002$
M96	$0.170\pm 0.001$	$0.121\pm 0.001$	$0.120\pm 0.002$

Table 5. Table 5: Scaling factors for the uncertainty in each observed parameter.

Parameter	$V$	${\rm Y}$	$L$	$q$	$q_0$	$D$	$m_0$
$ϵ$	3.0	1.7	$>10$	7	$>10$	2.1	5

Table 6. Table 6: Inverse RAR scatter measurements.

1	2	3	4
Survey	${\sigma }_{\mathrm{RAR}}^i$	${\sigma }_{\mathrm{sim}}^i$	${\sigma }_{\mathrm{int}}^i$
	(dex)	(dex)	(dex)
SV	$0.276\pm 0.011$	$0.188\pm 0.006$	$0.202\pm 0.016$
SP	$0.229\pm 0.005$	$0.184\pm 0.004$	$0.137\pm 0.010$
SF	$0.175\pm 0.002$	$0.177\pm 0.001$	$\text{-}0.027\pm 0.014$
C97	$0.251\pm 0.002$	$0.176\pm 0.001$	$0.180\pm 0.003$
M92	$0.277\pm 0.002$	$0.207\pm 0.002$	$0.184\pm 0.004$
M96	$0.237\pm 0.002$	$0.192\pm 0.001$	$0.140\pm 0.003$

Table 7. Table 7: Scatter measurements for the RAR as a function of mass bin, column (1).

1	2	3	4	5	6	7
Mass Bin	Survey	${\sigma }_{\mathrm{RAR}}$	${\sigma }_{\mathrm{sim}}$	${\sigma }_{\mathrm{int}}$	$N_{\mathrm{obs}}$	$N_{\mathrm{sim}}$
log${}_{10}{}^{}\left(\frac{M_{\mathrm{dyn}}}{M_{\odot }}\right)$		dex	dex	dex	#	#
8-10	SV	0.178	0.093	0.152	879	822
	SP	0.144	0.091	0.112	376	402
	SF	-	-	-	0	0
	C97	0.179	0.107	0.144	357	594
	M92	0.160	0.153	0.047	788	1056
	M96	0.210	0.124	0.170	512	829
10-11	SV	0.144	0.145	-0.018	501	862
	SP	0.116	0.125	-0.047	796	768
	SF	0.130	0.094	0.090	2222	1329
	C97	0.193	0.100	0.164	9036	8281
	M92	0.173	0.130	0.115	10,963	11,708
	M96	0.172	0.125	0.119	13,170	13,449
11-13	SV	-	-	-	0	0
	SP	0.134	0.142	-0.048	1387	1246
	SF	0.138	0.119	0.070	14,426	15,844
	C97	0.172	0.095	0.143	8233	9190
	M92	0.161	0.124	0.102	6604	5637
	M96	0.147	0.113	0.095	7778	6859

Equations

$${\mathrm{g}}_{\mathrm{obs}}=\frac{{\mathrm{g}}_{\ast }}{1-{\mathrm{e}}^{-\sqrt{{\mathrm{g}}_{\ast }/{\mathrm{g}}_{†}}}},$$

$${\mathrm{g}}_{\mathrm{obs}}=\frac{{\left({\mathrm{V}}_{\mathrm{obs}}/\sin (\mathrm{i})\right)}^2}{\theta \mathrm{D}},$$

$${\cos }^2(i)=\frac{q^2-q_0^2}{1-q_0^2},$$

$${\left[\frac{{\sigma }_{{\mathrm{g}}_{\mathrm{obs}}}}{{\mathrm{g}}_{\mathrm{obs}}}\right]}^2={\left[\frac{2{\sigma }_{V_{\mathrm{obs}}}}{V_{\mathrm{obs}}}\right]}^2+{\left[\frac{2q{\sigma }_q}{1-q^2}\right]}^2+{\left[\frac{2q_0{\sigma }_{q_0}}{1-q_0^2}\right]}^2+{\left[\frac{{\sigma }_{\mathrm{D}}}{\mathrm{D}}\right]}^2,$$

$${\mathrm{g}}_{\ast }=\frac{\mathrm{G}{\mathrm{{\rm Y}}}_{\mathrm{x}}{\mathrm{L}}_{\mathrm{x}}}{(\theta \mathrm{D}{)}^2}=\frac{\mathrm{G}{\mathrm{{\rm Y}}}_{\mathrm{x}}10^{-({\mathrm{m}}_{\mathrm{x}}-{\mathrm{M}}_{\odot ,\mathrm{x}})/2.5}}{{\theta }^2},$$

$${\left[\frac{{\sigma }_{{\mathrm{g}}_{\ast }}}{{\mathrm{g}}_{\ast }}\right]}^2={\left[\frac{\ln (10)}{2.5}{\sigma }_{m_x}\right]}^2+{\left[\frac{\ln (10)}{2.5}{\sigma }_{m_0}\right]}^2+{\left[\frac{{\sigma }_{{\rm Y}}}{{\rm Y}}\right]}^2.$$

$${\rho }_{\mathrm{disk}}(R,Z)=\Sigma (R)\frac{{\mathrm{sech}}^2(Z/{\mathrm{Z}}_0)}{2Z_0},$$

$$V=\sin (i)\sqrt{{\mathrm{g}}_{\mathrm{obs}}\theta \mathrm{D}}$$

$$L_x=\frac{{\mathrm{g}}_{\ast }(\theta \mathrm{D}{)}^2}{G{{\rm Y}}_x},$$

$$\mathcal{L}({\mathrm{g}}_{\ast ,1}^{\prime },{\mathrm{g}}_{\ast ,2}^{\prime },...,{\mathrm{g}}_{\ast ,\mathrm{N}}^{\prime },q^{\prime },q_0^{\prime },{{\rm Y}}_x^{\prime },D^{\prime })=$$

$$P(q^{\prime }|q,\sigma q)\cdot P(q_0^{\prime }|q_0,\sigma q_0)\cdot P({{\rm Y}}_x^{\prime }|{{\rm Y}}_x,\sigma {{\rm Y}}_x)\cdot$$

$$P(D^{\prime }|D,\sigma D)\cdot {\Pi }_{i=1}^{N}P(V_i^{\prime }|V,\sigma V)\cdot P(L_i^{\prime }|L_i,\sigma L_i),$$

$${\sigma }_{1^{\mathrm{st}}\mathrm{Order}}=\frac{1}{N}\sum ^N\frac{1}{{\mathrm{g}}_{\mathrm{obs}}}\sqrt{{\sigma }_{{\mathrm{g}}_{\mathrm{obs}}}^2+{\left(\frac{\partial \mathcal{F}}{\partial {\mathrm{g}}_{\ast }}{\sigma }_{{\mathrm{g}}_{\ast }}\right)}^2},$$

Full text

galaxies: general – galaxies: spiral – galaxies: kinematics and dynamics – galaxies: statistics – methods: statistical – techniques: miscellaneous

††journal: ApJ

1 Introduction

The correlation between two observables resulting from a physical causation is often referred to as a scaling relation. In the context of galaxies, these relations may be driven by formation and evolution processes, often revealing a subtle interplay between the baryonic and dark matter components. The slope, normalization, and scatter of these scaling relations are sensitive to galaxy formation parameters and mechanisms such as star formation efficiencies, merger histories, and the coupling between baryons and dark matter. For instance, the velocity-luminosity relation and fundamental planes of late- and early-type galaxies can be reproduced through a fine balance between gravitational forces and radiative processes (Courteau et al. 2007; Dutton et al. 2007, 2011; Trujillo-Gomez et al. 2011).

The scatter of an observed relation is the sum of observational and intrinsic (cosmic) components. Observational uncertainties arise from a variety of measurement limitations, including low signals, foreground/background interlopers, imprecise distance indicators, fuzzy parameter definitions, ill-defined data quality control, selection biases, and more. Conversely, the intrinsic scatter results from physical processes that have shaped galaxies since formation (Dutton et al. 2007; Somerville & Davé 2015). Intrinsic scatter can be readily compared with predictions from numerical simulations of galaxies in order to discriminate viable formation models. The accurate assessment of the intrinsic scatter of scaling relations is thus of utmost value for galaxy formation studies.

In this paper, we present a technique to extract the intrinsic scatter of a scaling relation by constructing a detailed model of the observational uncertainties. Doing so requires following a few operational steps. Firstly, a scaling relation is constructed with observational data. Two ideally uncorrelated variables are paired, and a global relation is inferred via an unbiased fitting method (see Sec. 3.3). Second, all observed data points are projected onto the fitted scaling relation. This provides new values for the core measurements under an assumption of zero intrinsic scatter. Third, the measurements are then resampled from their observational uncertainty distributions, centered on the new zero-scatter values. Fourth, the scaling relation including only observational errors is constructed using the resampled data points. It is now possible to directly compare the observed scaling relation with the one containing only observational errors, allowing one to infer the intrinsic variations.

We demonstrate the efficacy of our method by applying it to the radial acceleration relation (RAR) (McGaugh et al. 2016, hereafter M16). The latter consists of recasting the velocity–mass relation, often referred to as the baryonic Tully–Fisher relation (BTFR), of galaxies into a form that involves centripetal accelerations. While the BTFR uses global quantities, the RAR is constructed from spatially resolved parameters. Specifically, the observed centripetal acceleration from spectroscopic measurements is compared with that inferred from the baryonic matter distribution inferred via photometry. As our database does not include gas masses, we focus on the stellar RAR which exhibits the same scatter as the baryonic RAR. The scatter of the stellar and baryonic RARs has been shown to be comparable for both Spitzer Photometry and Accurate Rotation Curves (SPARC) observations (L17) and NIHAO simulations (A. Dutton 2019, private communication). Discrepancies between the observed and expected radial accelerations of galaxies have typically been ascribed to a dark matter particle component (Bertone et al. 2005; Courteau et al. 2014; Bertone & Hooper 2018).

Lelli et al. (2017, hereafter L17) examined the scatter of the residuals in the baryonic RAR ($\sigma_{\rm{RAR}}$) for a heterogeneous compilation of 153 spiral galaxies, labeled “SPARC,” and found it to be $\sigma_{\rm{RAR}}=0.13\,\rm{dex}$. They also found that a Gaussian fit to the residuals returns $\sigma_{\rm{RAR}}=0.11\,\rm{dex}$. This indicates that the data are not normally distributed, otherwise the two scatter measurements would be equal. A first-order analysis, assuming Gaussianity, of the observational uncertainties yields $\sigma_{\rm{obs}}\approx 0.12\,\rm{dex}$. As this estimate is close to their total observed scatter value, L17 claimed that the baryonic RAR is consistent with zero intrinsic scatter. Moreover, since $\Lambda$CDM theory predicts a nonzero intrinsic scatter about the relation Ludlow et al. (2017); Keller & Wadsley (2017); Dutton et al. (2019); L17 suggested that their inferred null scatter result favors Modified Newtonian Dynamics (MOND; Milgrom 1983).

Given the widespread success of $\Lambda$CDM, the claim by L17 remains controversial and calls for further scrutiny. Navarro et al. (2017) showed that the baryonic RAR arises naturally from $\Lambda$CDM using a standard galaxy formation model. The EAGLE, MUGS2, and NIHAO $\Lambda$CDM simulations also yield estimates for the intrinsic baryonic RAR scatter of $\sigma_{\rm{int},\Lambda\rm{CDM}}\approx 0.06-0.08\,\rm{dex}$ (Ludlow et al. 2017; Keller & Wadsley 2017; Dutton et al. 2019).

It was pointed out by L17 and Wheeler et al. (2018) that a galaxy that obeys the BTFR  (Walker 1999; McGaugh et al. 2000; Brook et al. 2016) typically obeys the baryonic RAR, except at the extremes of large and small radii not easily accessible with observations. When comparing the BTFR and RAR, the choice of axes for each relation is also relevant. The BTFR is normally depicted as the stellar mass (larger relative error) on the Y-axis and the circular velocity (smaller relative error) on the X-axis (McGaugh et al. 2000; McGaugh 2012; Ponomareva et al. 2018; Lelli et al. 2019). The corresponding representation for the RAR therefore calls for the inferred acceleration, $\rm{g}_{\rm{bar}}$ (larger relative error), on the Y-axis and the measured acceleration, $\rm{g}_{\rm{obs}}$ (small relative error), on the X-axis. The opposite approach was adopted by L17. Their choice of axes for the RAR thus naturally favors a small forward scatter, $\Delta Y(X)$. Indeed, the inverse scatter of the baryonic RAR, $\Delta X(Y)$, for the SPARC sample is nearly twice as large, on average, as the forward scatter (see Appendix A). The ratio $\Delta Y(X)/\Delta X(Y)$ grows with radius for the SPARC sample, as the relation flattens to a slope of 0.5 for small $\rm{g}_{\rm{bar}}$. As we also show in Sec. 5.2, the largest contribution to the errors is the uncertainty of the stellar mass-to-light ratio, $M_{*}/L$. A forward scatter measurement thus calls for the variables $M_{*}$ and $\rm{g}_{\rm{bar}}$ on the Y-axis for the BTFR and RAR, respectively; this is discussed further in Appendix A.

On a related note, Rodrigues et al. (2018) found that the individual galaxies in the SPARC sample each favor slightly different acceleration scales for $\rm{g}_{\dagger}$ indicating that its value is not universal, which MOND requires. Clearly, this is a case where the intrinsic scatter of a scaling relation bears immediate consequences for the interpretation of its physical underpinnings.

In this paper, we examine the scatter of the stellar RAR in detail, using a statistically compelling collection of resolved photometric and kinematic profiles for spiral galaxies coupled with a comprehensive error analysis. Whereas simulations seem to converge on $\sigma_{\rm{int},\Lambda\rm{CDM}}\approx 0.06-0.08\,\rm{dex}$, our approach consists of inferring $\sigma_{\rm{int},\rm{universe}}$ through careful elimination of all observational errors from the observed scatter.

In Sec. 2, we briefly review the stellar RAR and discuss the various sources of error that contribute to its overall scatter. This is followed in Sec. 3 by a description of six large galaxy surveys, collectively referred to as the Photometry and Rotation curve OBservations from Extragalactic Surveys (PROBES) catalog, providing structural parameters for 2500 galaxies to be used for our extensive RAR analysis. Our description of each galaxy survey includes salient features that are relevant for the stellar RAR. Sec. 4 addresses the dominant sources of uncertainty in the stellar RAR and how these are modeled for our analysis, whereas Sec. 5 describes the technique used to determine the observational error contribution to the stellar RAR scatter. In Sec. 6, the stellar RAR scatter is computed for each survey using a full Monte Carlo uncertainty model. Here the intrinsic scatter, $\sigma_{\rm{int}}$, is assessed by comparing the observed total scatter and the observational uncertainties in quadrature. The value of $\sigma_{\rm{int}}$ is then compared with estimates from numerical galaxy formation models, as well as values obtained by L17. Our measurement of $\sigma_{\rm{int}}$ serves as an empirical validation of $\Lambda$CDM models. For simplicity, unless otherwise stated, “RAR” means “stellar RAR” throughout. The “baryonic RAR” will be explicitly stated when needed.

2 The RAR

First reported in M16, the RAR is a tight relationship between the observed spatially resolved radial acceleration profile of a galaxy and that expected from baryonic matter alone. The RAR is essentially a translation of the mass discrepancy acceleration relation (MDAR; McGaugh 2004). The RAR is also closely related to the BTFR, as the two can be shown to match if $V(r)$ is constant and $M_{\rm bar}(r)\approx M_{\rm bar,total}$, which is typically where the BTFR is constructed (Lelli et al. 2017; Wheeler et al. 2018). As stated in Sec. 1, the BTFR is typically represented with baryonic mass on the Y-axis, while L17’s RAR has the baryon-dependent quantity on the X-axis, hence contributing to the remarkably tight forward scatter, $Y(X)$, observed by L17 as discussed in Sec. 1.

The baryonic mass distribution (and resulting acceleration) of a galaxy is primarily composed of two parts, stellar and gaseous, such that $\rm{g}_{\rm{bar}}=\rm{g}_{\rm{*}}+\rm{g}_{\rm{gas}}$, where $\rm{g}_{\rm{bar}}$ is the observed baryonic radial acceleration, and $\rm{g}_{\rm{*}}$ and $\rm{g}_{\rm{gas}}$ are the contributions to $\rm{g}_{\rm{bar}}$ from the stars and gas. Save for the SPARC data set, the other PROBES surveys (Sec. 3) all lack neutral gas measurements. For those, a stellar RAR ($\rm{g}_{\rm{obs}}$ versus $\rm{g}_{\rm{*}}$) can be readily constructed. The scatters in the stellar and baryonic RAR are comparable, as seen in L17 observations and NIHAO simulations (A. Dutton 2019, private communication). The stellar RAR is thus also a strong constraint on galaxy formation and dark matter models, but with the benefit of being more easily accessible observationally than the baryonic RAR.

The baryonic RAR was parametrized by M16 with a MOND-inspired single-parameter function. We present the equation from M16 in Eq. 1, now using $\rm{g}_{\rm{*}}$ instead of $\rm{g}_{\rm{bar}}$, as we will be examining the stellar RAR,

[TABLE]

where $\rm{g}_{\rm{*}}$ is the stellar radial acceleration, $\rm{g}_{\rm{obs}}$ is the observed radial acceleration, and $\rm{g}_{\dagger}$ is the expected universal MONDian acceleration scale equal to approximately 10^{-10}\,\text{m\,s^{-2}}. This MOND-motivated fitting function has a slope of unity (1:1) for large $\rm{g}_{\rm{*}}$ commonly found in the central regions of a galaxy. In the limit that $\rm{g}_{\rm{*}}$ is small, where standard galaxy formation models typically ascribe high dark matter fractions, that function becomes $\rm{g}_{\rm{obs}}=\sqrt{\rm{g}_{\dagger}\rm{g}_{\rm{*}}}$ . Thus, the representation of the RAR in Eq. 1 is an alternative representation of MOND (Milgrom 1983).

The scatter of the RAR, or any other scaling relation, can be decomposed as the sum in quadrature of its intrinsic and observed components: $\sigma_{\rm{RAR}}^{2}=\sigma_{\rm{int}}^{2}+\sigma_{\rm{obs}}^{2}$, where $\sigma_{\rm{RAR}}$ is the scatter of the scaling relation residuals, $\sigma_{\rm{int}}$ is the intrinsic scatter of the relation, and $\sigma_{\rm{obs}}$ is the scatter due to observational uncertainties/errors. The observational uncertainties of the RAR can be approximated to first order as: $\sigma_{\rm{obs}}^{2}\approx\sigma_{\rm 1^{st}Order}^{2}=\sigma_{\rm{g}_{\rm{obs}}}^{2}+\left(\frac{\partial\mathcal{F}}{\partial\rm{g}_{\rm{*}}}\sigma_{\rm{g}_{\rm{*}}}\right)^{2}$, where $\sigma_{\rm{g}_{\rm{obs}}}$ is the observational uncertainty in $\rm{g}_{\rm{obs}}$, and $\left(\frac{\partial\mathcal{F}}{\partial\rm{g}_{\rm{*}}}\sigma_{\rm{g}_{\rm{*}}}\right)^{2}$ is the observational uncertainty in $\rm{g}_{\rm{*}}$ modulated by the functional form of the RAR relation (Eq. 1). The values $\sigma_{\rm{g}_{\rm{obs}}}$ and $\sigma_{\rm{g}_{\rm{*}}}$ can be further broken down based on the variables used in the calculation of $\rm{g}_{\rm{*}}$ and $\rm{g}_{\rm{obs}}$. As shown below, numerous variables enter this calculation. A first-order uncertainty propagation is also presented.

2.1 Observed Radial Acceleration

The local observed radial acceleration $\rm{g}_{\rm{obs}}$ is computed from a galaxy rotation curve (RC) with the formula

[TABLE]

where $V_{\rm{obs}}$ is the measured line-of-sight velocity, $\theta$ is the angular radius of the velocity measurement, $D$ is the distance to the galaxy, and the inclination $i$ is inferred from the axial ratios of the galaxy. The latter, which represents the galaxy’s photometric tilt relative to the line of sight, is computed according to

[TABLE]

where $q=b/a$ is the measured axis ratio of the semi-major, $a$, and semi-minor, $b$, axes of the isophote, and $q_{0}=h_{Z}/h_{R}$ is a parameter representing the intrinsic flattening of a galaxy (the ratio of disk scale height, $h_{Z}$, to scale length, $h_{R}$). If $q<q_{0}$, an inclination of $90^{\circ}$ is used. Typical values for $q_{0}$ are in the range 0.1 to 0.25 (Kregel et al. 2002; Dutton et al. 2005; Hall et al. 2012). While the value of $q_{0}$ for a given galaxy cannot be measured directly, some correlations with other galaxy parameters do exist. For instance, $h_{z}$ scales with $V_{max}$ in edge-on galaxies (Kregel et al. 2005); however $V_{max}$ is an inclination- (and thus $q_{0}$-) dependent quantity and ill-suited for this analysis. Here $q_{0}$ also correlates with morphological type (T), and we use three flattening formulations based on T-Type, $q_{0}^{[1-3]}=0.20\pm 0.03$, $q_{0}^{[4]}=0.17\pm 0.03$, and $q_{0}^{[5-10]}=0.12\pm 0.02$ where the superscript represents the T-Type (Haynes & Giovanelli 1984; de Grijs 1998).

The uncertainty for each introduced variable plays a significant role in determining $\sigma_{\rm{g}_{\rm{obs}}}$, except $\sigma_{\theta}$ which is assumed to be negligibly small. Distance uncertainties enter into the calculation of $\rm{g}_{\rm{obs}}$ only via $D$, as the velocity and inclination measurements are distance-independent. Using Eqs. 2 and 3, one may compute the first-order uncertainty on $\rm{g}_{\rm{obs}}$,

[TABLE]

where $\sigma_{x}$ is the observational uncertainty on variable $x$. While this is formally the uncertainty on a single $\rm{g}_{\rm{obs}}$ measurement, it includes $D$ and $i$ which are shared variables between all $\rm{g}_{\rm{obs}}$ values for a given galaxy. As the first order calculation does not account for shared variables, it cannot yield accurate predictions. However, this simplified calculation is still useful for comparison with a full Monte Carlo model (Sec. 5), and the analysis performed in L17. Depending on the relative values of $\sigma_{V_{\rm{obs}}}$ versus $\sigma_{\rm{D}}$ and $\sigma_{\rm{i}}$, $\sigma_{\rm{g}_{\rm{obs}}}$ could be dominated by the local velocity measurement or shared galaxy variables ($D$ and $i$). Once again, the RAR data are highly correlated, and a first-order error analysis will incorrectly predict $\sigma_{\rm{obs}}$.

2.2 Stellar Radial Acceleration

The stellar radial acceleration, $\rm{g}_{\rm{*}}$, requires the specification of a three-dimensional spatial model for the stars, as the observed photometric values only yield a projected mass distribution. In the simplest case, the stellar mass can be assumed to lie in spherical shells, and, for a given curve of growth, $\rm{g}_{\rm{*}}$ can be represented as

[TABLE]

where $G$ is the gravitational constant, $\Upsilon_{x}$ is the stellar mass-to-light ratio in a photometric band $x$ in solar units, and $L_{x}$ and $m_{x}$ are the luminosity and apparent magnitude in band $x$ respectively, enclosed by an isophote at angular radius $\theta$. Here, $M_{\odot,x}$ is the absolute magnitude of the Sun in band $x$. The quantity $\rm{g}_{\rm{*}}$ is fortunately independent of distance errors, as the distance dependence of the luminosity $L_{x}$ is canceled by the acceleration formula. The mass-to-light ratio, $\Upsilon_{x}$, is either constant (for the 3.6\,\text{\mum} band) or determined by a colour mass-to-light ratio formula and thus independent of distance. The gravitational constant $G$ and photometric normalization $M_{\odot,x}$ are considered to have negligible error. However, $\Upsilon_{x}$ can have considerable uncertainty (Conroy 2013; Courteau et al. 2014; Roediger & Courteau 2015), and the uncertainty on $m_{x}$ can be significant, especially for fainter galaxies. Using Eq. 5, the first-order uncertainty on $\rm{g}_{\rm{*}}$ is computed to be

[TABLE]

As in the $\rm{g}_{\rm{obs}}$ equation (Eq. 2), some variables have shared uncertainty at all points in a single galaxy. For instance, all isophotal magnitudes share the same photometric zero-point $m_{0}$ and its uncertainty. Also, $m_{x}$ is correlated with all magnitude estimates interior to it, and assuming a constant $\Upsilon_{x}$ means that its uncertainty is shared across all points in the galaxy. This means, like $\rm{g}_{\rm{obs}}$, that the photometric data points $\rm{g}_{\rm{*}}$ are not independent, and a first-order analysis will incorrectly predict $\sigma_{\rm{g}_{\rm{*}}}$.

Eq. 6 is generated under the assumption that a galaxy’s mass distribution is spherical. In reality, disk galaxies are flattened structures whose potential is calculated by solving Poisson’s equation, $\nabla^{2}\phi=4\pi G\rho$, where $\rho$ is the three-dimensional mass density. The acceleration at each location on the disk is then computed as $\rm{g}_{\rm{*}}=-\nabla\phi$. The value of $\rm{g}_{\rm{*}}$ at a given radius depends strongly on the full mass distribution. Therefore, individual $\rm{g}_{\rm{*}}$ values for disk galaxies will be even more correlated than in the spherical case. In our analysis, we use the three-dimensional density distribution from van der Kruit & Searle (1981),

[TABLE]

where $\Sigma(R)$ is the radial projected surface mass density, $Z$ is the height above the disk, and $Z_{0}$ is the sech scale height of the disk. Here $h_{Z}$ can be determined from $q_{0}$ and $h_{R}$ by definition ($h_{Z}=q_{0}h_{R}$; see Sec. 2.1) and we use $Z_{0}=2h_{Z}$ to find the sech scale height (Kregel et al. 2002). With a density distribution specified, Poisson’s equation is solved by galpy (Bovy 2015) using a self-consistent field method (Hernquist & Ostriker 1992). This treatment is only applied to observed galaxies and thus the observed RAR. For the Monte Carlo model described in Sec. 5, the mock galaxies are treated as spherically symmetric, given the difficulty in inverting Poisson’s equation. In the next section, we discuss the galaxy samples used to construct and analyze the RAR.

3 Data sets

Our study relies on the compilation of six distinct surveys of 2634 disk galaxies, including the SPARC dataset of 163 galaxies assembled by Lelli et al. (2016) from various heterogeneous small datasets (see Sec. 3.1.5 below). The full sample, collectively referred to as PROBES will be presented elsewhere (C. Stone et al. 2019, in preparation). The six surveys presented below were collated in order to overcome small sample limitations, such as SPARC used on its own, and to mitigate selection biases. For instance, with only 163 galaxies (albeit with extended HI RCs), SPARC suffers from small number statistics. Numerous SPARC galaxies were also handpicked, leaving open the potential for unintentional bias. Our description of the six samples below highlights various features, such as selection criteria, data quality cuts, and distance estimates, that may either benefit or hinder the determination of an unbiased RAR. The maximal extent of available surface brightness (SB) profiles and RCs are reported below in terms of $R_{23.5}$, or the isophotal radius corresponding to 23.5\,\text{mag\,arcsec^{-2}}, in the relevant photometric band.

Each selected survey includes spatially resolved light profiles, in at least one photometric band, as well as spatially resolved RCs usually extracted from H$\alpha$ long-slit spectra or HI synthesis maps. If only one photometric band is available for the spatially resolved light profiles, stellar mass-to-light ratios are recovered via global colors estimated by the authors or retrieved from the NASA Extragalactic Database (NED). Spatially resolved HI fluxes to compute gas masses at all galactocentric radii and, ultimately, baryonic masses are only available for the SPARC data set. Therefore, mass computations include only stellar masses inferred via suitably chosen $M_{*}/L$ transformations (Roediger & Courteau 2015; Zhang et al. 2017). Missing gas masses in the RAR analysis affect the relation normalization but not its scatter, which is central to this paper, as discussed briefly in Sec. 2.2.

3.1 Our Survey Collection

3.1.1 M92

Mathewson et al. (1992, hereafter M92) acquired H$\alpha$ RCs and $I$-band photometry for 744 galaxies, mostly in the southern hemisphere. Each galaxy has a measured global B-I color, except for 51 galaxies where color information was retrieved from NED. Uncertainties for spatially resolved quantities are not reported for the SB profiles or RCs. Instead, we use simple models, described in Sec. 4, to estimate reasonable uncertainties. Repeat measurements by M92 yielded typical SB errors less than \sim 0.05\,\text{mag\,arcsec^{-2}} and velocity errors less than 10\,\text{km\,s^{-1}}. The median SB profile extends to 1.2$R_{23.5}$, and the median RC profile extends to 0.9$R_{23.5}$. Galaxies were selected for this survey primarily from the ESO-Uppsala catalog (Lauberts 1982, 1998) with morphological types of Sb to Sd, diameters greater than $1.7^{\prime}$, radial velocities typically below 7000\,\text{km,$s$^{-1}}, inclination above $40^{\circ}$, latitude $|b|$ greater than $11^{\circ}$, and a small number of galaxies from other surveys (M92).

3.1.2 M96

Mathewson & Ford (1996, hereafter M96) extended the M92 sample with 1216 additional galaxies with H$\alpha$ RCs and $I$-band photometry, primarily in the southern hemisphere. Many of these galaxies do not have B-I global colors and so global colors were once again retrieved from NED for a total of 399 galaxies. Similar to M92, uncertainties for spatially resolved quantities are not provided, and the models described in Sec. 4 are used. The sampling criteria were similar to those of M92, except with radial velocities in the range 4000 to 14,000\,\text{km\,s^{-1}}, and apparent diameters between $1.0^{\prime}$ and $1.6^{\prime}$, again with a small number taken from other surveys, such as the Uppsala General Catalog. The median SB profile extends to 1.2$R_{23.5}$, and the median RC profile extends to 0.8$R_{23.5}$.

3.1.3 C97

The Courteau (1997, hereafter C97) sample is a collection of 296 Sb- and Sc- type galaxies with H$\alpha$ RCs and $r$-band photometry. These were collected largely for cosmic flow studies for which systematic and random uncertainties were of great interest (Courteau et al. 1993, hereafter C93). Thus, many galaxies have repeat measurements, with some having as many as four remeasured RCs. More than half of the galaxies have multiple SB profiles. For this analysis, wherever multiple integrations exist, the deepest RC or SB profile is used. Repeat measurements are still valuable for the purpose of assessing uncertainties. The median SB profile uncertainty is 0.03\,\text{mag\,arcsec^{-2}}, and the median RC profile uncertainty is 6\,\text{km\,s^{-1}}. The median SB profile extends to 1.3$R_{23.5}$, and the median RC profile extends to 1.0$R_{23.5}$. The sample was selected from the Uppsala Catalog of Galaxies (Nilson 1995; Lauberts 1998), and the catalog of cluster galaxies from Bothun et al. (1985). The galaxies were selected to have Hubble types Sb-Sc, Zwicky magnitude $m_{B}\leq 15.5$, blue galactic extinction less than $0.5\,\text{mag}$ (based on Burstein & Heiles 1984), inclinations between $55^{\circ}$ and $75^{\circ}$, blue major axes less than $4^{\prime}$, and to be noninteracting/merging and have no overlapping bright stars (Courteau 1996).

3.1.4 Courteau et al. (2000)

Shellflow, from Courteau et al. (2000), is a sample of 171 galaxies with H$\alpha$ RCs and both $V$- and $I$-band photometry. The survey was designed to study an all-sky shell in redshift space to measure a cosmological bulk flow of galaxies with high precision. The Shellflow sample geometry meant that a large fraction of the galaxies could be observed from both the northern and southern hemisphere observatories, namely, KPNO and CTIO, thus mitigating calibration errors from using different instrumentation. The multiband photometry enables radially resolved mass-to-light ratios. The median SB profile uncertainty is 0.04\,\text{mag\,arcsec^{-2}} and the median RC profile uncertainty is 6\,\text{km\,s^{-1}}. The full distributions can be found in Secs. 4.4 and 4.3. The median SB profile extends to 1.4$R_{23.5}$, and the median RC profile extends to 0.9$R_{23.5}$. The Shellflow sample was selected from the Optical Redshift Survey by Santiago et al. (1995). Galaxies were chosen to be noninteracting and of morphological types of Sb and Sc, with radial velocities between 4500 and 7000\,\text{km\,s^{-1}}, inclinations between $45^{\circ}$ and $78^{\circ}$, $A_{B}$ extinctions less than $0.3\,\text{mag}$ (as determined by Burstein & Heiles 1982), and having no bright overlapping foreground stars or tidal disturbances.

3.1.5 L16

The Spitzer Photometry and Accurate Rotation Curves (SPARC) sample compiled by Lelli et al. (2016, hereafter L16) is an amalgamation of over 50 smaller samples totaling 163 galaxies222an additional 12 “low quality” galaxies are not used in this paper with Spitzer 3.6\,\text{\mum} photometry and HI RCs. Approximately one-third of the SPARC galaxies have hybrid HI and H$\alpha$ RCs, to combine the higher spatial resolution provided by H$\alpha$, with the extensive radial extent of synthesis HI radio maps, where available. The distances to SPARC galaxies rely on a number of methods, including Hubble flow, tip of the red giant branch, Cepheids, Ursa Major cluster distance, and supernovae. All distances and their uncertainties are included in the survey; we use them directly for our analysis. The median SB profile uncertainty is 0.01\,\text{mag\,arcsec^{-2}} and the median RC profile uncertainty is 4.0\,\text{km\,s^{-1}}. The full velocity uncertainty distribution can be found in Sec. 4.3; resolved magnitude uncertainties are not reported for SPARC and are therefore lacking in Sec. 4.4. The median SB profile extends to 0.7$R_{23.5}$, and the median RC profile extends to 1.2$R_{23.5}$. The sample was carefully chosen to have a wide range of morphology, luminosity, and SB from the limited selection of available galaxies with HI measurements (L16). The photometry is homogeneously collected in the Spitzer 3.6\,\text{\mum} band, RCs are compiled from 56 separate studies (L16) excluding THINGS (de Blok et al. 2008) and LITTLE-THINGS (Oh et al. 2015). Each RC in SPARC is assigned a quality flag, with $Q=1$ and $2$ considered acceptable and $Q=3$ (12 objects) having strong non-circular motions and/or asymmetric features, making them unsuitable for the analysis in L17.

3.1.6 Ouellette et al. (2017)

The Spectroscopy and $H$-band Imaging of Virgo Cluster Galaxies (SHIVir) survey presented in Ouellette et al. (2017) is a dedicated survey of galaxies in the Virgo Cluster with dynamical and multiband information. The sample was selected to examine the impact of the cluster environment on galaxy properties using the nearest cluster in the sky. For our RAR analysis, we focus on the subset of 44 SHIVir spiral galaxies with H$\alpha$ RCs from long-slit spectra and ugriz-H photometry. While SHIVir is the smallest survey in our RAR analysis, its spatially resolved multiband imaging makes it valuable for studying the importance of resolved versus global mass-to-light ratios, as mass-to-light ratios are a substantial source of uncertainty in the RAR. Distances to the SHIVir galaxies are measured using SB fluctuations where available (Jerjen et al. 2004; Blakeslee et al. 2009); otherwise a standard value of $16.5\,\text{Mpc}$ is assumed (Mei et al. 2007). Uncertainties are reported where distance measurements are available; for the rest, a median value of $3\,\text{Mpc}$ is used. The median SB profile uncertainty is 0.03\,\text{mag\,arcsec^{-2}} and the median RC profile uncertainty is 4\,\text{km\,s^{-1}}, the full distributions can be seen in Secs. 4.4 and 4.3. The median SB profile extends to 1.4$R_{23.5}$, and the median RC profile extends to 0.6$R_{23.5}$. The SHIVir galaxies are a subset of 286 galaxies drawn from the Virgo Cluster Catalog (VCC), based on the intersection of the VCC catalog and SDSS 6th Data Release, which is volume complete in a spatial subset of the Virgo Cluster to an absolute magnitude of $M_{B}\leq-15.15$, along with several fainter galaxies, to ensure a broad morphology coverage (McDonald et al. 2011).

3.2 Data Quality Selections

A number of data quality criteria were implemented in our compilation. Only galaxies with inclinations greater than $30^{\circ}$ were considered; rotational velocities are corrected by $\sin(i)^{-1}$ and would therefore diverge at lower inclination. This study focuses solely on late-type galaxies, and only galaxies with morphological types from Sa to Im are considered (T-Type of 1 - 10). Only a few other types were available in each survey and were not beneficial to this analysis. Individual RC data points were discarded if $\sigma_{V}/V>0.1$, as in L17. Individual SB data points were discarded if \sigma_{\mu}>0.1\,\text{mag\,arcsec^{-2}} or $\mu>25$. Data points within $5^{\prime\prime}$ of the center of a galaxy were also discarded to avoid seeing effects.

3.3 The RAR for Each Survey

We now present the RAR for each survey in Fig. 1. Each panel shows a fit to the data using Eq. 1; the fit uses an orthogonal distance regression  (Jones et al. 2001). The axes are calculated using the analysis described in Secs. 2.1 and 2.2 plus an extra $0.33\,\text{dex}$ added to the $\rm{g}_{\rm{*}}$ axis to account for the missing gas mass to the measured stellar mass values. This is done largely for aesthetic and fitting reasons (properly positioned and reasonable fitted $\rm{g}_{\dagger}$ values) and plays no role in our study of the RAR scatter. The RAR scatter, $\sigma_{\rm RAR}$, is measured relative to a running median in a window chosen to include $\sim 100$ datapoints. A parametric assessment of the scatter would be biased by the arbitrariness of the RAR functional form and suffer from an artificial increase to the scatter due to any misalignment with the data. It is thus avoided here.

The values for $\rm{g}_{\dagger}$ in Table 1 differ by many standard deviations. It is, however, important to note that that table only includes bootstrap random errors. The large systematic uncertainties that affect variables such as the mass-to-light ratios, which could exceed a factor of 2, are not included in this treatment. Accounting for them would easily reconcile all of the values of $\rm{g}_{\dagger}$ with each other (Courteau et al. 2014; Roediger & Courteau 2015). Therefore, strong statements regarding the universality of $\rm{g}_{\dagger}$ based on our six surveys cannot be made. As in Fig. 1 and Table 1, each survey is analyzed separately throughout the paper. Systematic differences between the surveys that could artificially inflate the observed scatter are thus minimized.

4 Uncertainty Model

Many parameters enter the computation of the stellar and observed radial accelerations, and they must each be characterized when adopting an appropriate uncertainty model. Once an error function is determined for each parameter, it can be sampled randomly to form the basis of the Monte Carlo model described in Sec. 5. We examine each major parameter and its related uncertainty below.

4.1 Distances

Depending on the method, distances may carry significant uncertainties and represent an important source of scatter in the RAR. Most galaxies in the SHIVir and SPARC surveys are close enough to get distance measurements with variable stars, SB fluctuations, cluster distance, tip of the red giant branch, or supernovae. Their different uncertainties are reported in each survey. For the more distant galaxies in M92, M96, C97, and Shellflow, Hubble flow distances are used  (Riess et al. 2016). A primary source of distance uncertainty when using this method is due to peculiar motions. Here we assume the peculiar motion to be of order \sigma_{pec}=300\,\text{km\,s^{-1}}, where $\sigma_{pec}$ is the peculiar velocity dispersion. Many galaxies in the surveys are far enough that this $\sigma_{pec}$ value is relatively small, so a lower relative uncertainty bound of $\sigma_{D}=0.15D$ is applied. Thus, for distance uncertainties, we have $\sigma_{D}=max(0.15D,~{}300)$ for Hubble flow distances.

4.2 Inclinations

Inclination uncertainties are modeled in multiple steps. First, a per-galaxy axis ratio uncertainty is computed using the scatter in isophotal ellipticity of the outer regions of the SB profile (ellipticities beyond $R_{e}$). These cover a large range of isophotal fitting uncertainties, as shown in Fig. 2. Second, the intrinsic flattening parameter $q_{0}$ is selected with a Gaussian uncertainty, truncated to the range $[0.05-0.25]$, representing typical acceptable values for the intrinsic flattening of a disk (Kregel et al. 2002).

These are combined through Eq. 3 to determine the adopted inclinations. Fig. 3 demonstrates the results of this model. For a given axis ratio $q$, a vertical slice in Fig. 3 gives the corresponding uncertainty distribution in inclination $i$.

4.3 Rotational Velocities

Rotational velocity uncertainties are reported for each survey except M92 and M96. For these, a uniform velocity uncertainty of \sigma_{V}=6\,\text{km\,s^{-1}} is used, representing the median uncertainty for the other surveys. Fig. 4 demonstrates the distribution of relative rotational velocity uncertainties for each survey. The black dashed line indicates the adopted cutoff threshold of $10\,\text{\%}$ (L17), where velocity measurements are not considered in the analysis. Velocity uncertainties account for measurement errors but ignore noncircular motions (see Sec. 6.2).

4.4 Curves of Growth

As with rotational velocities, magnitude uncertainties are reported for each survey except M92 and M96. The reported uncertainties are depicted in Fig. 5 as a function of the SB values. We find that magnitude uncertainties correlate more tightly with SB than magnitude measurements, hence this choice of variable. The dashed red line in the Shellflow panel represents the SB uncertainties used for M92 and M96. This is done by fitting $\sigma_{\rm SB}=a+e^{b(\mu-c)}$ to the values that meet the data quality threshold (black dashed lines). The Shellflow distribution is adopted, as it relies on the same photometric band as M92 and M96; the resulting parameterization for $(a,b,c)$ is (0.00075, 0.52, 31.97). Magnitude uncertainties are small compared to other sources of uncertainty; however, they are retained for consistency.

4.5 Stellar Mass-to-light Ratios

Stellar mass-to-light ratios also carry a significant model uncertainty. The exact normalization of the mass-to-light ratio is inconsequential to our analysis; only the intrinsic variability about the relation matters. Roediger & Courteau (2015) found the random uncertainty of the stellar mass-to-light ratio to be of order $0.09-0.13\,\text{dex}$ for optical mass-to-light ratios. A representative conservative uncertainty of $0.13\,\text{dex}$ is therefore used. Systematic uncertainties (e.g. model-to-model differences) for the stellar mass-to-light ratio of the order of $0.3\,\text{dex}$ (Conroy 2013; Courteau et al. 2014; Roediger & Courteau 2015) are of course larger; however, this uncertainty is not used for our analysis, as a systematic shift (constant factor in log space) does not impact scatter measurements. For the uncertainty in the 3.6\,\text{\mum} stellar mass-to-light ratio, an uncertainty of $0.11\,\text{dex}$ is used, as suggested by Meidt et al. (2014). The effect of choosing even larger stellar mass-to-light ratio uncertainties is partially explored in Sec. 6.2.

5 Monte Carlo Scatter Model

5.1 Projecting onto the Stellar RAR

To assess the scatter induced by observational uncertainties, one must eliminate the effects of intrinsic scatter (if any) in the RAR. As the observational and intrinsic scatter are mixed, let us first eliminate all scatter. To this end, all data points in the $\rm{g}_{\rm{*}}$ $\rm{g}_{\rm{obs}}$ space are projected onto the fitted RAR to impose a zero-scatter assumption. To project a data point onto the RAR, its final scatter-free location must first be identified. This is equivalent to asking what new values of ${\rm g}_{\rm*,zs}$ and ${\rm g}_{\rm obs,zs}$, where “zs” means zero scatter (constrained to the RAR), should be chosen. One option is to project along the $\rm{g}_{\rm{*}}$ or $\rm{g}_{\rm{obs}}$ axis; however, this means artificially selecting an axis. A second option is to adopt an orthogonal projection, where each ${\rm g}_{\rm*,zs}$ and ${\rm g}_{\rm obs,zs}$ is the minimum distance (in $\rm{g}_{\rm{*}}$, $\rm{g}_{\rm{obs}}$ space) from the measured $\rm{g}_{\rm{*}}$, $\rm{g}_{\rm{obs}}$. However, this assumes that the uncertainties in each axis are identical, which may not be true either. Instead, a third and favored option is to consider the full uncertainty distributions from Sec. 4 for each parameter, construct a likelihood function by multiplying the probabilities, and choose the location that maximizes the likelihood whilst imposing a scatter-free RAR (see Eq. 1).

As every parameter is given a probability distribution in Sec. 4, it is possible to construct a likelihood that is a function of each variable but constrained to lay on the RAR. This removes one degree of freedom for each point on the RAR. The equations from Secs. 2.1 and 2.2 can be used to connect $\rm{g}_{\rm{obs}}$ and $\rm{g}_{\rm{*}}$ to their respective measurements,

[TABLE]

where $\rm{g}_{\rm{obs}}$ is determined using Eq. 1, thus constraining it to lay on the RAR. Some of the variables are shared between multiple observations in the same galaxy (distance, axis ratio, intrinsic thickness, and mass-to-light ratio), and so the optimization must be performed for all observations simultaneously. For a galaxy with $N$ observations on the RAR, the likelihood function can be represented as

[TABLE]

where the primed quantities are the new zero-scatter values and the unprimed quantities are the measurements and their uncertainties. Here $V_{i}^{\prime}$ and $L_{i}^{\prime}$ are computed from Eqs. 8 and 9 respectively, using the primed quantities; $\rm{g}_{\rm{obs}}$ is computed from Eq. 1. Ultimately, this process finds the values for all quantities that are most consistent with the original measurements while remaining constrained on the RAR. The values for each variable are now taken as the “true” values, upon which the observational uncertainties can scatter measurements off the RAR. This process is described in the next section.

5.2 Reintroducing Observational Uncertainties

Using the uncertainty models from Sec. 4, new “observations” can be sampled about the zero-scatter values as described in Sec. 5.1. To visualize the effect of this Monte Carlo model, Fig. 6 shows how the different uncertainties affect the final location of a single point off of the RAR. New data points are sampled for each galaxy and compiled into the same data structure as the original measurements. As the RAR is a local scaling relation, each galaxy places multiple data points on the relation; however, some variables are global for the whole galaxy. For example, a single distance value is resampled for each galaxy. Thus, it is clear that the data points in the RAR are not independent, as multiple measured values with large uncertainties are shared.

To match the statistics of the original observations, each point is resampled only once. Thus, the simulated data have exactly the same number of initial entries as the measured data. To include the effects of data quality cuts, all measured values are resampled before the data quality cuts are applied. This means that some points in the RAR that were initially cut may be newly introduced, for example, if a galaxy had an inclination below $30^{\circ}$ but was resampled above that value.

The mock data are then processed with the same code as the original data, thus duplicating any data quality cuts, interpolations, and conditional statements. This critical element is missing from a first-order analysis (e.g. L17) and can potentially impact the final result, as will be seen in Sec. 6.2. An element of L17’s analysis that is not reproduced here is the flat-disk model for the baryonic acceleration. Inverting Poisson’s equation adds unnecessary complexity without impacting the scatter measurements that are the primary focus of this paper. In Sec. 5.1, the projection onto the RAR is made under the assumption that the galaxies are spherical; thus, the calculations remain internally consistent. In both the flat-disk analysis and spherical assumption, the baryonic acceleration depends linearly on luminosity and mass-to-light ratio, which are the two dominant sources of uncertainty in $\rm{g}_{\rm{bar}}$. The uncertainty model for these quantities should behave identically for the two baryon distributions.

6 Results

The PROBES catalog provides a consistent set of data quality cuts and a detailed uncertainty model, enabling a robust characterization of the intrinsic scatter in the RAR. In order to do so, we perform a first-order analysis of all uncertainties using a simplified error model as in L17. This model incorporates the same uncertainties as our Monte Carlo model (Sec. 5), but only treated to first order. Next, we consider the full uncertainty model with Monte Carlo sampling to get an accurate measure of the observational uncertainties. The first-order analysis and the full uncertainty model can then be used to extract the intrinsic scatter from the data; the results are finally compared for each method.

6.1 First-Order Scatter Prediction

A first-order calculation of the average observational uncertainty for the SPARC data set was presented by L17. Here we reproduce this calculation for the full PROBES sample. Eqs. 4 and 6 present the propagation of uncertainties through the calculation of $\rm{g}_{\rm{*}}$ and $\rm{g}_{\rm{obs}}$. In L17 the assumed luminosity uncertainty is $\sigma_{L}=0.04\,\text{dex}$, and the assumed mass-to-light ratio uncertainty is $\sigma_{\Upsilon}=0.1\,\text{dex}$, meaning that the baryonic acceleration uncertainty is always $\sigma_{\rm{g}_{\rm{bar}}}/\rm{g}_{\rm{bar}}=0.11\,\text{dex}$. In this paper, we consider the same luminosity uncertainty and more conservative values for the mass-to-light ratio uncertainty of $\sigma_{\Upsilon}^{[3.6\mu m]}=0.11\,\text{dex}$ (Meidt et al. 2014). For optical wavelengths, we use $\sigma_{\Upsilon}^{\rm optical}=0.13\,\text{dex}$ (Roediger & Courteau 2015; Zhang et al. 2017). While both $\rm{g}_{\rm{*}}$ and $\rm{g}_{\rm{obs}}$ contribute to the scatter in the forward direction, the stellar acceleration uncertainty is modified by the RAR as the slope changes across the relation. The average scatter in the forward residuals takes the form

[TABLE]

where $\sigma_{\rm 1^{st}Order}\approx\sigma_{\rm obs}$ is an estimate of the scatter in the residuals due to observational uncertainties, $N$ is the number of observations, and $\mathcal{F}$ is the functional form of the RAR (Eq. 1). Table 2 presents the results of this calculation for each survey, where the first column represents the survey, the second column reports the observed scatter of the RAR, the third column gives the first-order scatter result from Eq. 11, and the fourth column reports the intrinsic scatter resulting from the difference in quadrature ($\sigma_{\rm int}^{2}=\sigma_{\rm RAR}^{2}-\sigma_{\rm 1^{st}Order}^{2}$). Each survey does indeed have different observational uncertainties, as expected for surveys with such a range of properties. The median value for the intrinsic scatter is then $0.13\,\text{dex}$; this can be compared to the results found in $\Lambda$CDM simulations, which range from 0.06 to $0.08\,\text{dex}$ (Keller & Wadsley 2017; Ludlow et al. 2017). The value determined from first-order calculations is substantially larger than expected from $\Lambda$CDM simulations and certainly not consistent with zero.

The values in Table 2 are close to those from L17. Our respective values of $\sigma_{\rm RAR}$ are essentially identical. For $\sigma_{\rm 1^{st}Order}$, we find $0.1\,\text{dex}$, while L17 quoted $0.12\,\text{dex}$. This discrepancy originates in slightly different choices in uncertainty values for some parameters. We used larger uncertainty values than L17 for the stellar mass-to-light ratios; for the inclination, we typically found smaller uncertainties with our model from Sec. 4.2. These differences account for most of the $0.02\,\text{dex}$ discrepancy.

6.2 Monte Carlo Simulation Results

Here we present the Monte Carlo simulation of each survey with only observational uncertainties about the RAR. With the full dataset resampled, any scatter metric (such as the orthogonal scatter) could be chosen to compare with the original data. For consistency with other studies and with the first-order analysis in Sec. 6.1, the forward residuals will be used as the scatter metric. Fig. 7 presents the simulated data with the same analysis as in Sec. 3.3. Many qualitative elements from Fig. 1 are retained; however, the distributions in Fig. 7 are clearly smoother and tighter. This visual impression is confirmed in Table 3 which lists the same parameters as Table 1, where the scatter measured for each simulated survey is below that of the observed one. Table 4 presents the scatter measurements for the mock data in the same format as Table 2.

As with our first-order scatter calculations in Sec. 6.1, one may compute the intrinsic scatter by taking the difference in quadrature with the observed scatter. Using the mock data (Table 4), the median intrinsic scatter is $0.11\pm 0.02$, where the uncertainty is the weighted (by points in the RAR) standard deviation of the intrinsic scatter values. These results are in closer agreement with $\Lambda$CDM expectations (Ludlow et al. 2017; Keller & Wadsley 2017; Dutton et al. 2019) than the first-order analysis. Dutton et al. (2019) noted that the baryonic RAR scatter depends on galaxy mass. We observe the same effect with PROBES galaxies but not our mock data. Appendix B explores the RAR as a function of mass bins. This sheds light on the range of scatter measurements amongst surveys that sample different mass regimes.

Several factors may potentially cause the Monte Carlo simulated scatter to underrepresent the full observational uncertainties. For instance, noncircular motions, which might increase the observed scatter, especially in the inner disk, are not modeled here. We find that the median intrinsic scatter for RAR points beyond one $R_{e}$, i.e. for regions relatively free of noncircular motions, is $0.10\pm 0.01\,\text{dex}$, still in excellent agreement with $\Lambda$CDM simulations (Ludlow et al. 2017; Keller & Wadsley 2017; Dutton et al. 2019).

It is also possible that some of the observational uncertainties have been underestimated. To assess that possibility, we consider a scaling factor that is applied to all uncertainty values for a given parameter. We consider one parameter $x$ at a time, with its uncertainty being scaled as $\sigma_{x}^{\prime}=\epsilon\sigma_{x}$ and $\epsilon$ as the scale factor, whilst all the other uncertainties are kept at their fiducial value. The scaling factor is tuned until the Monte Carlo simulated scatter is in agreement with the observed scatter (indicating zero intrinsic scatter). The final scaling factor value is reported in Table 5.

This exercise suggests that, in order to explain away all intrinsic scatter in the RAR, the parameter errors would need to be inflated by factors of $2$ (stellar mass-to-light ratio) to $10$ (total luminosity and disk flattening ratio). Considering the already conservative values adopted for the uncertainty in each variable, it seems unlikely that the scaled uncertainties represent a reasonable path to explain the nonzero scatter. We reiterate that the nonzero scatter is most likely a true feature of the data, in agreement with $\Lambda$CDM simulations.

7 Conclusions

We have presented a method for determining the intrinsic scatter about any scaling relation using the stellar RAR as a test case. The database used in this paper represent the largest collection of extended spatially resolved RCs and SB profiles to date. With over 2500 galaxies sampled over a representative range of galaxy properties, the detailed study of numerous scaling relations such as the stellar RAR becomes possible.

As demonstrated in M16, the baryonic RAR presents an interesting local scaling relation with the capacity to test alternative dark matter models. We have examined the intrinsic scatter in the stellar RAR with great care in order to accurately model observational uncertainties. All six PROBES galaxy surveys, each with different selection criteria, wavelength coverage, and instrumental setups, were analyzed simultaneously to minimize systematic biases. Our Monte Carlo modeling of the observational uncertainties enables a simultaneous treatment of all nonlinear effects. Ultimately, we find that the intrinsic scatter in the stellar RAR is of order $0.11\pm 0.02\,\text{dex}$, which broadly agrees with $\Lambda$CDM simulations that predict $0.06-0.08\,\text{dex}$ (Keller & Wadsley 2017; Ludlow et al. 2017). In order to explain away any intrinsic scatter that we ascribe to $\Lambda$CDM effects, the observational uncertainties would have to be considerably inflated from their already conservative estimates (Table 5). Ultimately, our conservative assessment of errors in the stellar RAR yields full consistency with $\Lambda$CDM expectations.

We are grateful to the Natural Sciences and Engineering Research Council of Canada, Ontario Government, and Queen’s University for support through various scholarships and grants. The referee is thanked for a thoughtful and constructive report. Nikhil Arora and Larry Widrow are also thanked for illuminating discussions. The software packages galpy and astroquery were used extensively in this analysis. We also acknowledge the NASA/IPAC Extragalactic Database for the wealth of information that it provides.

Appendix A Inverse RAR Scatter

We have presented the RAR with $\rm{g}_{\rm{*}}$ (the stellar mass-dependent quantity) on the x-axis and $\rm{g}_{\rm{obs}}$ (the velocity-dependent quantity) on the y-axis. This choice impacts scatter measurements, which have so far been measured from forward (vertical) residuals. Inverse scatter is measured horizontally in the RAR diagram. Rather than measuring the scatter of $\rm{g}_{\rm{*}}$ for a given $\rm{g}_{\rm{obs}}$ (forward method), we now compute the inverse scatter of $\rm{g}_{\rm{obs}}$ for a given $\rm{g}_{\rm{*}}$ (inverse method). At high acceleration, the slope of $\rm{g}_{\rm{obs}}$ versus $\rm{g}_{\rm{*}}$ is roughly 1, meaning the forward and inverse scatters are the same. However, at low acceleration, the slope of $\rm{g}_{\rm{obs}}$ versus $\rm{g}_{\rm{*}}$ is $\sim$0.5; thus, the forward scatter is about half that of the inverse relation.

Table 6 shows the inverse scatter calculated using the same techniques as in Table 4. In some cases (e.g. Shellflow), the simulated scatter can be slightly larger than the observed scatter; this unphysical result simply results from the statistical nature of our analysis, and the intrinsic scatter is consistent with zero.

The median of all of the intrinsic scatter measurements is $\sigma_{\rm int}^{i}=0.16\pm 0.04\,\text{dex}$, where the uncertainty $\pm 0.04$ is the weighted standard deviation of the $\sigma_{\rm int}^{i}$ values. Thus, we find that the intrinsic scatter value is statistically equal to the forward scatter, albeit with a much larger variance indicating that the scatter is poorly constrained.

Overall, we adopt a median observed forward and inverse scatter of $\sigma_{\rm RAR}=0.17$ and $\sigma_{\rm RAR}^{i}=0.24\,\text{dex}$, respectively. The choice of forward or inverse residuals greatly impacts the resulting scatter. Needless to say, this consideration affects all galaxy scaling relations.

A comparison between RAR and BTFR scatters is warranted once the different units are accounted for. Hall et al. (2012) examined the BTFR scatter in both axes, finding $\sigma_{VM_{\rm bar}}=0.078\,\text{dex}$ and $\sigma_{M_{\rm bar}V}=0.274\,\text{dex}$. The forward and inverse scatters are again quite different, with the forward scatter being numerically smaller than the baryonic RAR scatter from L17. However, in order to properly compare with the RAR, the forward scatter of the BTFR must be multiplied by  $\sim 1.5$, since on one axis, the BTFR scales with $V$ and the RAR with $V^{2}$. The scatter of the pseudo-BTFR is thus $0.078*1.5=0.12\,\text{dex}$, in very close agreement with the baryonic RAR scatter found by L17 ($0.11\,\text{dex}$).

Appendix B Mass-dependent Stellar RAR Scatter

This appendix presents a study of the RAR scatter as a function of mass bins. Dutton et al. (2019) found that the baryonic RAR scatter decreases with increasing mass, finding a scatter of $\sigma_{\rm RAR,\Lambda CDM}=0.11\,\text{dex}$ in the $M_{*}$ range [$10^{7},10^{9.3}$] $M_{\odot}$, and $\sigma_{\rm RAR,\Lambda CDM}=0.04\,\text{dex}$ in the range [$10^{9.3},10^{11}$] $M_{\odot}$. A number of factors may explain this trend. Firstly, as noted in Dutton et al. (2019), the low- and high-mass galaxies exhibit different RC shapes, and the extra scatter at low mass is likely a reflection of the diversity of dwarf galaxy RCs (Oman et al. 2015, 2019). Secondly, covariant errors between the RAR axes (such as with $q_{0}$ seen in Fig. 6) can result in data points moving along the relation in the high-acceleration regime; however, in the low-acceleration regime, where the RAR flattens, these errors boost the scatter. Thirdly, the higher observed scatter at lower masses also reflects the higher relative velocity errors on the corresponding galaxy RCs. Most velocity measurements have errors of order 3-10\,\text{km\,s^{-1}} (Fig. 4).

Since high- and low-mass galaxies have peak velocities around 200 and 50 km s*-1*, respectively, the relative velocity (and thus acceleration) errors are higher for low-mass galaxies.

We can quantify the mass dependence of the RAR scatter and use the mock data to determine if the dependence is intrinsic. Table 7 addresses this issue by examining three dynamical mass bins, where the total mass is computed within the maximal extent of each RC (typically near $R_{23.5}$). Some of the measurements produce negative intrinsic scatter values (e.g. for SHIVir and SPARC). These are cases where our conservatively large uncertainty estimates break down, and the negative intrinsic scatter can be viewed as statistical fluctuations.

A clear trend of decreasing scatter with increasing mass is detected for the observed data; for these, the median values of the RAR scatter for each bin are 0.178, 0.158, and 0.147, respectively. No clear trend is detected for the simulated data where, the median scatter values per mass bin are 0.107, 0.125, and 0.119, respectively. Thus, the intrinsic scatter, calculated by subtracting the observed scatter from the simulated scatter in quadrature, shows a distinct dependence on mass with median values of 0.144, 0.103, and 0.095, respectively. We conclude that there is a clear trend of decreasing intrinsic scatter for the observed RAR as a function of mass, independent of signal-to-noise variations, likely bolstering the case for enhanced RC diversity at low masses.