The Radial Acceleration Relation in Disk Galaxies in the MassiveBlack-II Simulation

TL;DR

This study examines the radial acceleration relation in disk galaxies within the MassiveBlack-II simulation, finding a strong correlation between baryonic and total accelerations, with some differences from observed data.

Contribution

It provides the first detailed comparison of the radial acceleration relation in a hydrodynamic galaxy formation simulation with observational data, highlighting similarities and differences.

Findings

Strong correlation between baryonic and total accelerations.

Radial acceleration relation can be described by a simple power law.

Differences observed between simulation results and real galaxy data.

Abstract

A strong correlation has been measured between the observed centripetal accelerations in galaxies and the accelerations implied by the baryonic components of galaxies. This empirical radial acceleration relation must be accounted for in any viable model of galaxy formation. We measure and compare the radial accelerations contributed by baryons and by dark matter in disk galaxies in the MassiveBlack-II hydrodynamic galaxy formation simulation. The sample of 1594 galaxies spans three orders of magnitude in luminosity and four in surface brightness, comparable to the observed sample from the Spitzer Photometry & Accurate Rotation Curves (SPARC) dataset used by McGaugh et al. (2016). We find that radial accelerations contributed by baryonic matter only and by total matter are highly correlated, with only small scatter around their mean or median relation, despite the wide ranges of galaxy luminosity and surface brightness. We further find that the radial acceleration relation in this simulation differs from that of the SPARC sample, and can be described by a simple power law in the acceleration range we are probing.

Tables

Table 1. Table 1: Best-fit values for ( A , B ) 𝐴 𝐵 (A,B) in Equation (missing)  3 and for g † subscript 𝑔 † g_{\dagger} in Equation (missing)  1 , and the corresponding residuals, for the full disk galaxy sample, disk samples split into bright and faint galaxies, and also central galaxies with subhalo mass, M halo > 10 12 ​ h − 1 ​ M ⊙ subscript 𝑀 halo superscript 10 12 superscript ℎ 1 subscript M direct-product M_{\text{halo}}>10^{12}\,h^{-1}\text{M}_{\odot}{} . The best-fit value of g † subscript 𝑔 † g_{\dagger} and residual of SPARC sample from (McGaugh et al., 2016 ) is also listed for comparison.

	Power-law fit (Eq. 3)			McGaugh et al. fit (Eq. 1)
Galaxy sample	$A$ (slope)	$B$ (intercept)	Residual	$g_{†}$ [${10}^{-10}\text{m}{\text{s}}^{-2}$]	Residual	# of data points
Full Disk Sample	$0.6849\pm 0.0004$	$-2.759\pm 0.004$	${0.013}_{-0.091}^{+0.061}$	$2.01\pm 0.004$	${0.053}_{-0.17}^{+0.068}$	159400
Bright Sample	$0.6811\pm 0.0006$	$-2.784\pm 0.006$	${0.006}_{-0.056}^{+0.047}$	$3.05\pm 0.007$	${0.023}_{-0.084}^{+0.052}$	71600
Faint Sample	$0.5523\pm 0.0012$	$-4.311\pm 0.014$	${0.007}_{-0.12}^{+0.11}$	$1.151\pm 0.003$	${0.007}_{-0.12}^{+0.11}$	87800
Centrals ($M_{\text{halo}}>{10}^{12}{h}^{-1}{\text{M}}_{\odot }$)	$0.6374\pm 0.0005$	$-3.215\pm 0.005$	${0.0008}_{-0.049}^{+0.050}$	$3.491\pm 0.005$	${0.0008}_{-0.041}^{+0.047}$	32648
SPARC sample	-	-	-	$1.2\pm 0.02$	$0{.}_{-0.055}^{+0.055}$	2693

Equations

$$g_{\text{obs}}=\frac{g_{\text{bar}}}{1-\exp \left(-\sqrt{g_{\text{bar}}/g_{†}}\right)},$$

$$j_{\text{circ}}(r)=r{V}_{\text{circ}}(r)=\sqrt{GM(<r)r}.$$

$${\log }_{10}\left(\frac{g_{\text{obs}}}{\text{m}{\text{s}}^{-2}}\right)=A{\log }_{10}\left(\frac{g_{\text{bar}}}{\text{m}{\text{s}}^{-2}}\right)+B,$$

Full text

††pubyear: 2017††pagerange: The Radial Acceleration Relation in Disk Galaxies in the MassiveBlack-II Simulation–References

1 Introduction

In the current standard cosmological model ($\Lambda$CDM), nearly 85% of the Universe’s matter content is composed of dark matter. This picture has been consistently corroborated by a variety of cosmological probes, such as the cosmic microwave background (Planck Collaboration et al., 2014), the growth of large-scale structure (Blumenthal et al., 1984), and Big Bang nucleosynthesis yields for light elements (Walker et al., 1991). On galaxy and cluster scales, gravitational lensing provides further evidence for the existence of dark matter with roughly the same relative abundance as that inferred from cosmology (Kilbinger et al., 2013; Abbott et al., 2016; Joudaki et al., 2017).

Early evidence of dark matter was found in the context of galactic dynamics. Zwicky (1933) suggested its existence as an explanation for the discrepancy between the luminous matter and the dynamical mass in the Coma cluster. A case for dark matter via measurements of galaxy rotation curves was made later by Rubin & Ford (1970). Rotation curves with the same characteristic shape have now been observed for thousands of spiral galaxies (e.g., Persic et al., 1996).

Recently, McGaugh et al. (2016) employed the new Spitzer Photometry & Accurate Rotation Curves (SPARC) database (Lelli et al., 2016) to report a correlation between the radial acceleration traced by rotation curves ($g_{\text{obs}}$) and the radial acceleration inferred from the observed distribution of baryons ($g_{\text{bar}}$). This correlation has long been known (see e.g., Sanders, 1990; McGaugh et al., 2000); however, the new analysis is particularly revealing because previous uncertainties in baryonic mass-to-light ratios are mitigated using infrared imaging. The SPARC analysis exploited a sample of 153 galaxies with a wide range of luminosities, effective surface brightnesses, and gas fractions. The authors found that the Radial Acceleration Relation (RAR) displayed by the data is well described by

[TABLE]

where $g_{\dagger}=1.20\pm 0.26\times 10^{-10}\,\text{m}\,\text{s}^{-2}$ corresponds to the acceleration scale below which $g_{\text{obs}}$ deviates from $g_{\text{bar}}$, with little scatter ($0.11$ dex). Notably, the correlation persists when dark matter dominates $g_{\text{obs}}$. This implies that the distribution of dark matter is completely determined by that of baryons or vice versa.

The RAR over a wide range of galaxy masses and surface brightnesses is not an obvious result of standard galaxy formation scenarios in dark matter haloes. Though it may not be difficult to imagine that a relation of this kind is obtained in standard $\Lambda$CDM cosmological models (e.g., Kaplinghat & Turner, 2002; Navarro et al., 2016), the tightness of the observed RAR may be more difficult to understand (e.g., Desmond, 2017). In any case, the observed RAR has led to speculation that its origin lies in an underlying departure from Newtonian gravity (e.g., Milgrom, 1983a, b).

Of course, it is certainly possible that the RAR is a natural outcome of galaxy formation in the standard cosmological model and any claims to the contrary should be based upon specific predictions of the standard cosmological model. Only recently have simulations including dark matter, hydrodynamics, star formation, and various forms of baryonic feedback attained sufficient mass resolution to form realistic disk galaxies (see e.g., Vogelsberger et al., 2014; Schaye et al., 2015; Khandai et al., 2015) such that the RAR can be studied within the context of the $\Lambda$CDM model for structure formation.

Indeed, Keller & Wadsley (2016) and Ludlow et al. (2016) studied the RAR in the MUGS2 (Stinson et al., 2010; Keller et al., 2014; Keller et al., 2015) and EAGLE (Schaye et al., 2015) simulations, respectively. While Keller & Wadsley (2016) were able to recover the RAR found in McGaugh et al. (2016), their sample of 18 galaxies with masses ranging from $1.8\times 10^{10}$ to $2.7\times 10^{11}\,\text{M}_{\odot}$ does not include representatives of the fainter galaxies present in SPARC. Moreover, the MUGS2 sample is sufficiently small that deviations from and/or scatter about the RAR could be missed simply due to finite sampling.

Ludlow et al. (2016) obtained a tight RAR in the EAGLE simulation, with an optimal fitting parameter $g_{\dagger}=3.0\times 10^{-10}\,\text{m}\,\text{s}^{-2}$, which is 2.5 times larger than that of the observed SPARC sample. However, Ludlow et al. (2016) do not quote the range of galaxy surface brightness of their simulated galaxy sample, and do not select disk galaxies, therefore it is difficult to compare the EAGLE sample with the SPARC sample in a robust manner.

The authors of the MUGS2 and EAGLE studies reach different conclusions regarding the evolution of the RAR with redshift. While Keller & Wadsley (2016) find the RAR to evolve with redshift in the MUGS2 simulation, Ludlow et al. (2016) report that the RAR remains unchanged at earlier times in the EAGLE simulation.

In this study, we use the MassiveBlack-II (MB-II) hydrodynamic simulation to analyze the RAR for disk galaxies spanning a wide range of masses and surface brightnesses, roughly comparable to those present in the SPARC data. We calculate the radial accelerations contributed by baryonic matter only ($g_{\text{bar}}$) and by all matter ($g_{\text{obs}}$), and find them tightly correlated over the full range of disk galaxies, with only small scatter of around $0.1$ dex, despite the wide range of luminosities and surface brightnesses spanned by the simulated galaxy sample. However, the RAR in this simulation differs from that of the SPARC sample. In particular, we find that the RAR in MB-II can be described by a single power law for accelerations $g_{\text{bar}}\lesssim 10^{-9}\,\text{m}\,\text{s}^{-2}$. This power law is steeper than the low-acceleration asymptote $g_{\text{obs}}\rightarrow\sqrt{g_{\dagger}\,g_{\text{bar}}}$ of Equation (missing) 1. Moreover, the SPARC data prefer a clear break from this asymptotic behavior at an acceleration of $g_{\dagger}\simeq 10^{-10}\,\text{m}\,\text{s}^{-2}$.

This paper is organized as follows. We introduce the MB-II simulation, the galaxy sample used in this work, and the measurement of the RAR in Section 2. We present our results in Section 3, and further discuss and summarize these results in Section 4.

2 Methods

2.1 MassiveBlack-II Simulation

In this study, we inspect the RAR in the MassiveBlack-II (MB-II) hydrodynamic simulation. MB-II is a state-of-the-art high-resolution, large-volume, cosmological hydrodynamic simulation of structure formation. This simulation has been performed with p-gadget, which is a hybrid version of the parallel code, gadget2 (Springel et al., 2005), upgraded to run on Peta-flop scale supercomputers. In addition to gravity and smoothed-particle hydrodynamics, the p-gadget code also includes the physics of multiphase interstellar medium model with star formation (Springel & Hernquist, 2003), black hole accretion, and AGN feedback (Springel et al., 2005; Di Matteo et al., 2012). Radiative cooling and heating processes are also included (as in Katz et al., 1996), as is photoheating due to an imposed ionizing ultraviolet background. The details of this simulation can be found in Khandai et al. (2015).

MB-II contains $N_{\text{part}}=2\times 1792^{3}$ dark matter and gas particles in a cubic periodic box with a side length of $100\,h^{-1}\text{Mpc}$, and a (Plummer-equivalent) gravitational softening length $\epsilon=1.85\,h^{-1}\text{kpc}$ in comoving units. A single dark matter particle has a mass $m_{\text{DM}}=1.1\times 10^{7}\,h^{-1}\text{M}_{\odot}{}$ and the initial mass of a gas particle is $m_{\text{gas}}=2.2\times 10^{6}\,h^{-1}\text{M}_{\odot}{}$, with the mass of each star particle being $m_{\text{star}}=1.1\times 10^{6}\,h^{-1}\text{M}_{\odot}{}$. The cosmological parameters used in the simulation are as follows: amplitude of matter fluctuations $\sigma_{8}=0.816$, spectral index $n_{s}=0.96$, mass density parameter $\Omega_{m}=0.275$, cosmological constant density parameter $\Omega_{\Lambda}=0.725$, baryon density parameter $\Omega_{b}=0.046$, and Hubble parameter $h=0.702$ as per WMAP7 (Komatsu et al., 2011).

2.2 Galaxy samples

To identify the morphological type of a galaxy in the MB-II simulation, we follow the procedure from Scannapieco et al. (2009) and define a circularity parameter, $\epsilon\equiv j_{z}/j_{\text{circ}}(r)$, for each star within 10 times the stellar half-mass radius. Here $j_{z}$ is the component of the specific angular momentum of the star in the direction of the total angular momentum of the galaxy, and $j_{\text{circ}}(r)$ is the specific angular momentum of a circular orbit at the same radius as the star,

[TABLE]

The total angular momentum of the galaxy is calculated using all star particles within 10 times the stellar half-mass radius.

All stars with $\epsilon>0.7$ are then identified as disk stars. The bulge-to-total ratio is hence defined as $\text{BTR}=1-f_{\epsilon>0.7}$, where $f_{\epsilon>0.7}$ is the fraction of disk stars within 10 times the stellar half-mass radius. We then classify galaxies with $\text{BTR}<0.7$ as disk galaxies.

In the left panel of Figure 1, we show the distribution of luminosity and surface brightness of all disk galaxies in this simulation. The luminosity is in the Spitzer $3.6\,\mu\text{m}$ band. To compare our sample with the observed SPARC sample and other simulated samples, we also outline the range of luminosities ($10^{7}<L_{[3.6]}<5\times 10^{11}L_{\odot}$) and surface brightnesses ($5<\Sigma_{*}<3\times 10^{3}\,\text{pc}^{-2}\,\text{L}_{\odot}$) of the SPARC sample with the black parallelogram, and show the range of luminosities in the MUGS2 (Keller & Wadsley, 2016) and the EAGLE (Ludlow et al., 2016) samples in red and green bands respectively. Note that the stellar masses are converted into luminosities based on the mass-to-light ratio of $0.5$ in the SPARC sample.

Among the samples shown in Figure 1, the SPARC sample has the largest dynamical range. Our galaxy sample in MB-II spans a much larger range compared to the other two simulated samples, and is broadly comparable with the SPARC sample. Compared with the SPARC sample, the MB-II sample has fewer high luminosity and high surface brightness galaxies. This difference is likely due to both selection effects and the limited volume of the MB-II simulation. We also note that in the MB-II sample, there seems to be a sharp decrease in population above $2\times 10^{9}\,\text{L}_{\odot}$. This is caused by our disk selection criteria. We also divide the disk sample into bright and faint galaxies along the first principal component (axis with greatest variance) in the $\log L_{*}-\log\Sigma_{*}$ plane, as shown in the right panel of Figure 1, to study the dependence of the RAR on luminosity.

2.3 Measuring and fitting of the RAR

For each galaxy, we calculate the radial acceleration at locations on the disk plane that are logarithmically spaced in distance from the center of the galaxy. The range of distances spans $0.25-4.0$ times the stellar half-light radius. At each of these locations, we estimate $g_{\text{obs}}$ (in the notation of McGaugh et al., 2016), which denotes the “observed” radial component of the gravitational acceleration contributed by all particles in the galaxy. The quantity $g_{\text{obs}}$ includes accelerations induced by dark matter as well as those induced by gas, stars, and black holes. At each of these locations, we also compute $g_{\text{bar}}$, which denotes the gravitational acceleration contributed by all particles except for dark matter particles. We calculate both $g_{\text{obs}}$ and $g_{\text{bar}}$ via direct summation of the gravitational accelerations projected onto the radial direction in the plane of the galaxy disk. We have also confirmed that, for all the disk galaxies in this sample, more than 85% of the particles have a radial acceleration that is larger than 90% of the magnitude of its total acceleration.

While $g_{\text{obs}}$ and $g_{\text{bar}}$ as defined here share the same notations and physical motivations as in those defined in McGaugh et al. (2016), the operational definitions are distinct. McGaugh et al. (2016) determine $g_{\text{obs}}$ by examining rotation speeds, and $g_{\text{bar}}$ by estimating the gravitational potential induced by the observed baryons. By contrast, the velocities of star particles never enter our calculation of $g_{\text{obs}}$ since accelerations can be determined by direct summation and do not need to be inferred from measured velocities. It would be interesting to examine, in a high-resolution zoom-in simulation, how well one can recover the actual mass distribution from the observed velocity dispersion. In the present paper, we focus our study on our idealized measurements of $g_{\text{obs}}$ and $g_{\text{bar}}$.

For both of our high- and low-luminosity samples, we fit the measured RAR with two different functional forms and inspect the residuals. The first function to which we fit our data is given in Equation (missing) 1 and is the relation used in McGaugh et al. (2016). This form has only one free parameter, $g_{\dagger}$. The other functional form we use here is a power-law relation,

[TABLE]

with two free parameters: the log–log slope $A$, and the intercept $B$. While there is no physical motivation to use a power-law relation, we find that the power-law relation gives good fits to the RAR (over a restricted range of accelerations) as we will demonstrate in the next section. However, we note that the power law cannot describe the relation at large radial acceleration as we expect $g_{\text{obs}}\geq g_{\text{bar}}$ to always hold, and hence the power law must break down at large acceleration.

Before proceeding to our primary results, it is worth noting that the fits can be affected significantly by the distribution of galaxy properties. The galaxies in the SPARC sample are (approximately) uniformly distributed within the parallelogram in Figure 1, whereas the galaxies in our simulated samples are volume limited and thus distributed quite differently in the $L_{*}-\Sigma_{*}$ plane. We weight the galaxies in our samples so as to mimic the weighting in the SPARC sample. To be specific, we identify the principal components of the galaxy distribution in the $\log L_{*}-\log\Sigma_{*}$ plane, based on a linear transformation that defines the axis with greatest variance as the first principal axis. The principal components are shown as red lines in the right panel of Figure 1. We then bin the galaxies in our sample according to the first principal component, and then weight the data from each galaxy by $1/N_{i}$, where $N_{i}$ is the population in the corresponding bin. While this improves the comparison between the MB-II and SPARC samples, it certainly does not account for the fact that the high-luminosity and high-surface brightness regions of the SPARC sample are not present in the simulated MB-II sample.

3 Results

Figure 2 shows the RAR for all the disk galaxies in the MB-II simulation. We fit the $g_{\text{obs}}-g_{\text{bar}}$ relation with the two functional forms introduced in Section 2.3, and compare them with the relation for the SPARC sample. The RAR of the simulated galaxies clearly deviates from that of the SPARC sample. In particular, the RAR of the simulated galaxies is a significantly steeper at low acceleration than that of the SPARC sample. As a consequence, Equation (missing) 1 is not a good description of the simulated data because it asymptotes to a power law $g_{\text{obs}}\propto\sqrt{g_{\text{bar}}}$ for $g_{\text{bar}}\ll g_{\dagger}$. The MB-II disk galaxies obey a relation closer to $g_{\text{bar}}\propto g_{\text{bar}}^{0.7}$. In fact, over the range of accelerations represented in the simulations, the RAR of the simulated disk galaxies is much better described by a single power-law relation than by Equation (missing) 1. At most, there is only a slight hint of a break from the power-law near an acceleration of $g_{\text{bar}}=10^{-9}\,\text{m}\,\text{s}^{-2}$, an order of magnitude larger than the transition acceleration $g_{\dagger}\simeq 10^{-10}\,\text{m}\,\text{s}^{-2}$ of the SPARC sample. The scatter in the simulated sample is relatively small around the best-fit power-law relation, with a width in the residual distribution of only around 0.1 dex. In detail, this scatter is a function of acceleration and is larger for lower accelerations. The residual distribution is also slightly skewed, as is evident in the inset in Figure 2.

To further study how the RAR depends on luminosity and surface brightness, we also divide the disk sample into bright and faint galaxies along the first principal component in the $L_{*}-\Sigma_{*}$ plane, as shown in the right panel of Figure 1. Figure 3 shows the RARs for both samples. The scatter in the bright sample is much smaller than that of the full sample, and the median relation (fitted by a power law) of the bright sample is similar to that of the full simulated disk sample and distinct from the SPARC sample. The RAR of the bright sample is what drives the fits to be formally distinct from the fit to the SPARC sample, because the RAR has a much steeper slope than the asymptotic low-acceleration slope of the SPARC sample.

The faint sample, shown in the right panel of Figure 1, spans a significantly narrower range of accelerations. However, over the range of accelerations where the faint simulated galaxy sample overlaps with the SPARC data, the two samples exhibit similar RARs. In particular, the slope of the RAR of faint galaxies in MB-II is much closer to that of the SPARC sample, with approximately $g_{\text{obs}}\propto\sqrt{g_{\text{bar}}}$. This is evident through the proximity of the green and purple curves in the right panel of Figure 1. The faint sample also exhibits much larger scatter in the RAR, which is likely due to the effect of finite resolution.

For a direct comparison with the EAGLE result presented in Ludlow et al. (2016), we analyze a sample of central galaxies in MB-II with a halo mass $M_{\text{halo}}>10^{12}\,h^{-1}\text{M}_{\odot}{}$, shown in Figure 4. The method to calculate radial acceleration also matches the calculation in Ludlow et al. (2016), where the $g_{\text{obs}}$ and $g_{\text{bar}}$ accelerations at a given radial distance from the center are obtained from the total mass and baryonic mass enclosed at that radius respectively. Nevertheless, the difference between the “enclosed mass” method and the “direct summation” method is negligible. The RAR of this sample is in good agreement with the result obtained using the EAGLE simulation, and also has a small scatter around the median relation. Yet the RAR of this sample still differs from that of the SPARC sample.

In Table 1 we list all the values of the best-fit parameters for each sample with both the functional form introduced in Section 2.3, and also the corresponding residuals. In many cases, the residual values in simulation do not closely follow a normal distribution; hence we report the median and the $16^{\text{th}}$ and the $84^{\text{th}}$ percentiles (in the notation of $m^{+1\sigma}_{-1\sigma}$). The error values of each of the fitting parameters reported in the table are the random (statistical) error, not the scatter in the distribution. Because of the large number of data points (roughly 100 locations for each galaxy in our sample), we obtain much smaller random error values.

4 Discussion and Summary

In this paper, we investigated the RAR in the disk galaxies of the MassiveBlack-II simulation at redshift $z=0.06$. A recent study by McGaugh et al. (2016), found a strong correlation between the observed radial acceleration and that due to baryons in the SPARC sample of galaxies for a wide range of luminosities, surface brightnesses, galaxy sizes and morphologies. This relation has been investigated in galaxy formation simulations by Keller & Wadsley (2016) in the MUGS2 cosmological zoom-in simulations and Ludlow et al. (2016) in the EAGLE simulation. However, the samples of simulated galaxies in these studies span a relatively limited range of luminosities and surface densities. In contrast, the galaxies in our sample have been morphologically identified based on a dynamical bulge-disc decomposition and span a significantly wider range of luminosities and surface brightnesses than previous theoretical studies.

In accord with the results of both observations (the SPARC sample) and other simulations, we find that the radial acceleration due to baryons is strongly correlated with the total radial acceleration at a variety of points throughout galactic disks and that this correlation exhibits fairly limited scatter. This result seems to imply that galaxy formation in a $\Lambda$CDM universe can, plausibly, induce a remarkably tight correlation between $g_{\text{obs}}$ and $g_{\text{bar}}$. Previous work has argued that the dark matter haloes, which dominate the local gravitational potential, may lead to a tight relation of this kind simply because galaxy sizes do not vary significantly and because haloes are nearly self-similar (e.g., Kaplinghat & Turner, 2002; Navarro et al., 2016). Indeed, a large body of circumstantial evidence suggests that the properties of dark matter haloes play the critical role in determining the properties of galaxies that lie within them.

However, the RARs of our simulated galaxy samples are also clearly distinct from the observed RAR of the SPARC sample. The slope of the RAR in the MB-II simulation is steeper than the asymptotic slope of the SPARC sample, and we do not see clear evidence for a transition from this asymptotic behavior in the acceleration range that we can probe. When being fitted to the SPARC functional form, Equation (missing) 1, the MB-II galaxies exhibit a significantly larger $g_{\dagger}$ value than the SPARC result. The slight deviation from the power-law fit at the high-acceleration end in Figure 2 does hint at a transition to the regime where $g_{\text{obs}}=g_{\text{bar}}$, yet in the regime where we have most data points there is not a clear transition.

The lack of this transition in our simulation is in part due to a very limited number of data points that have large acceleration, which reside in the inner region of massive galaxies. Although we do the best to mimic the SPARC population, the distribution of disk luminosities and surface brightnesses available in the MB-II simulation is still significantly different from that of the SPARC sample, particularly with a deficit of high surface brightness galaxies. Also, the finite gravitational softening length prevents accurate simulation of the high acceleration inner regions of massive disks. A statistical sample of individual disks which span the transition region can provide one of the most incisive tests of the observed RAR and pin down the transition acceleration. This test, however, would require simulations that have both better resolution and larger volumes, which are still beyond our current reach.

The effect of the finite resolution is twofold. First, resolution limits the number of particles in smaller galaxies and results in larger shot noise. Second, the finite spatial resolution limits the baryonic physics that can be resolved, which creates larger uncertainty (and possibly systematic biases). For example, the resolution limit manifests in our investigation of the bright and faint samples, as we find that the scatter (or residual with respect to the fitted relation) is significantly smaller in the bright sample. These two samples also have different median RARs. The best-fit RAR of the bright sample is very close to that of the full disk sample. When we remove data points that are within four times of the softening length ($4\epsilon$), the resulting RAR stays approximately unchanged, but with much smaller scatter at the low-acceleration range, since the best-fit RAR of the full disk sample is dominated by the bright sample.

In this work we also made a direct comparison with the result from the EAGLE simulation (Ludlow et al., 2016), by selecting a similar sample of galaxies. The results from the MB-II and EAGLE simulations are in good agreement, and both have smaller scatter. It appears that for this sample of central galaxies, the difference in the detail baryonic feedback between the MB-II and EAGLE simulations does not affect the RAR significantly. This is an interesting result in and of itself, as it suggests that the RAR is a quantity that is robust to these choices. However, the RARs in both simulations differ from that of the observational measurements using the SPARC sample. The EAGLE and MB-II simulations also have similar mass resolution (while EAGLE adopts a smaller softening length), and hence both are subject to the impact of finite resolution as discussed above.

The degree to which the RAR is sensitive to baryonic physics implementations in different simulations remains somewhat unclear. For example, Keller & Wadsley (2016) found that the RAR evolves with redshift in the MUGS2 sample, and Ludlow et al. (2016) find that the relation remains unchanged in the EAGLE sample. At face value this discrepancy seems to suggest that the evolution of the RAR could be sensitive to the implementation of baryonic physics. However, we should note that the MUGS2 and EAGLE samples span very different galaxy populations, as shown in Figure 1. Hence, it is possible that the aforementioned discrepancy is caused by the mass (or luminosity) dependence of the evolution of the RAR. The MUGS2 simulations, though with a small sample size and a limited luminosity range, have better resolution compared to the EAGLE simulation, which may also contribution to the different evolution of the RAR. In this work we limit our study to the present-day population, and leave the redshift evolution of the RAR to follow-up studies.

In conclusion, the RAR does pose a significant challenge for the understanding of galaxy formation, whether within the context of the standard $\Lambda$CDM scenario or in some other models. Any viable model of galaxy formation must yield the observed RAR. Therefore, the RAR represents a new and interesting test with which to confront simulations of galaxy formation in the hope of improving our understanding of the mechanisms that lead to the illumination of the universe. Though current simulations do lead to an RAR with little scatter, in detail the simulated RAR differs significantly from that in the observed SPARC sample. Questions such as whether the RAR has a clear transition to $g_{\text{obs}}=g_{\text{bar}}$ and whether the RAR depends on galaxy luminosity and surface brightness also need to be answered. Further investigation of these issues will require larger-volume galaxy formation simulations and higher-resolution zoom-in simulations of individual galaxies.

Acknowledgements

We thank Stacy McGaugh and Jeffery Newman for helpful discussions. The MassiveBlack-II simulation was run on the Cray XT5 supercomputer Kraken at the National Institute for Computational Sciences. YYM is supported by the Samuel P. Langley PITT PACC Postdoctoral Fellowship. RACC and TDM acknowledge support from NSF AST-1616168, NSF ACI-1614853, NSF ACI-1036211, NSF AST-1517593, NSF AST-1009781, NSF OCI-0749212, the NSF PAID program and RACC from NSF AST-1412966. AK and FZ have been supported by NSF AST-1312380. FZ acknowledges support from the Andrew Mellon Predoctoral Fellowship. ARZ has been supported by NSF AST-1517563 and NSF AST-1516266.

This research made use of many community-developed or community-maintained software packages, including (in alphabetical order): Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), and SciPy (Jones et al., 01). We also use the NASA Astrophysical Data System for bibliographic information.