NIHAO XVIII: Origin of the MOND phenomenology of galactic rotation curves in a LCDM universe
aa r X i v : . [ a s t r o - ph . GA ] F e b MNRAS , 1–15 (2019) Preprint 20 February 2019 Compiled using MNRAS L A TEX style file v3.0
NIHAO XVIII: Origin of the MOND phenomenology ofgalactic rotation curves in a Λ CDM universe
Aaron A. Dutton, ⋆ Andrea V. Macci`o, , Aura Obreja, , Tobias Buck New York University Abu Dhabi, PO Box 129188, Abu Dhabi, United Arab Emirates Max-Planck-Institut f¨ur Astronomie, K¨onigstuhl 17, 69117 Heidelberg, Germany University Observatory Munich, Scheinerstraße 1, D-81679 Munich, Germany
20 February 2019
ABSTRACT
The phenomenological basis for Modified Newtonian Dynamics (MOND) is the radial-acceleration-relation (RAR) between the observed acceleration, a = V ( r )/ r , and theacceleration accounted for by the observed baryons (stars and cold gas), a bar = V ( r )/ r .We show that the RAR arises naturally in the NIHAO sample of 89 high-resolution Λ CDM cosmological galaxy formation simulations. The overall scatter from NIHAO isjust 0.079 dex, consistent with observational constraints. However, we show that thescatter depends on stellar mass. At high masses ( ∼ < M star ∼ < M ⊙ ) the simulatedscatter is just ≃ . dex, increasing to ≃ . dex at low masses ( ∼ < M star ∼ < M ⊙ ).Observations show a similar dependence for the intrinsic scatter. At high masses theintrinsic scatter is consistent with the zero scatter assumed by MOND, but at lowmasses the intrinsic scatter is non-zero, strongly disfavoring MOND. Applying MONDto our simulations yields remarkably good fits to most of the circular velocity profiles.In cases of mild disagreement the stellar mass-to-light ratio and/or “distance” can betuned to yield acceptable fits, as is often done in observational mass models. In dwarfgalaxies with M star ∼ M ⊙ MOND breaks down, predicting lower accelerations thanobserved and in our Λ CDM simulations. The assumptions that MOND is based on(e.g., asymptotically flat rotation curves, zero intrinsic scatter in the RAR) are onlyapproximately true in Λ CDM . Thus if one wishes to go beyond Newtonian dynamicsthere is more freedom in the observed RAR than assumed by MOND.
Key words: methods: numerical – galaxies: fundamental parameters – galaxies:haloes – galaxies: kinematics and dynamics – dark matter
The observational evidence for dark matter goes back manydecades. Zwicky (1933) showed that if the Virgo cluster isbound, its total mass must greatly exceed the sum of themasses of the individual members. In the 1970’s galaxy ro-tation curves were found to level off at large radii, implyingthe presence of substantial mass outside the optically visibledimensions of galaxies (Faber & Gallagher 1979).Further evidence for missing mass comes from thegrowth of cosmic structure. The Cosmic Microwave Back-ground (CMB) is remarkably uniform with temperature fluc-tuations ∆ T / T = ∆ ρ / ρ ∼ − (Smoot et al. 1992). In or-der for these density fluctuations to grow into the observedgalaxies ∆ ρ / ρ ∼ , they need to increase by a factor of !In the baryon-dominated universe, in order to have nonlinear ⋆ [email protected] structure now, the induced fluctuations in the CMB wouldhave to be much larger than is observed. In a dark matterdominated universe, dark matter perturbations grow priorto recombination, and this induces a rapid growth in baryonperturbations shortly after recombination, allowing nonlin-ear structure to form now from a smaller baryon fluctuationat recombination (Bond et al. 1980).An alternative to dark matter was proposed by Milgrom(1983) who considered the possibility that rather than therebeing missing mass in galaxies, there is modification to New-tonian dynamics (MOND) at low accelerations a < a , where a ∼ − m s − . The observed centripetal acceleration (inthe plane of the disk) is defined as a ( r ) = V ( r )/ r = − ∂ Φ ( r )/ ∂ r , (1)where Φ is the total gravitational potential. Likewise, theacceleration predicted due to the observed baryons, a bar = © Dutton et al.
Figure 1.
Commonly used MOND interpolation functions of the form a µ ( a / a ) = a bar , and a = a bar / ν ( a bar / a ) . Left: a vs a bar on alog-scale (radial acceleration relation), right: − a bar / a vs a bar ( mass discrepancy - acceleration relation) with a linear vertical scale tohighlight the differences between fitting functions. The dotted line shows the 1:1 relation. V / r . The following assumptions restrict the form of themodification to Newtonian dynamics. Assumption 1: there is no missing mass at high accel-erations. Thus a ( r ) = a bar ( r ) , ( a ≫ a ) . (2) Assumption 2: rotation curves are asymptotically flatat low accelerations. At large radii the total acceleration isgiven by a ( r ) = V / r , and the enclosed baryonic mass is aconstant: a bar ( r ) = V ( r )/ r = GM bar / r . Requiring that a isa function of a bar and not radius implies V = GM bar a / a bar . (3)Since V flat , G , and M bar are constants, this implies a ( r )/ a bar ( r ) ≡ a , ( a ≪ a ) (4)is also a constant. Note that Eq. 3 is a Baryonic Tully-Fisherrelation with slope of 4. Assumption 3: there is a unique interpolation func-tion, µ ( x ) , between the two acceleration extremes: a ( r ) µ ( a / a ) = a bar ( r ) . (5)However, this interpolation function is not specified by thetheory. A popular choice is µ ( x ) = x / p + x . (6)This can be inverted to express a in terms of a bar (van den Bosch & Dalcanton 2000): a ( r ) = a bar ( r ) (cid:18) + q + a / a (cid:19) / . (7)Other choices and their inversions include µ ( x ) = x /( + x ) ; (8) a ( r ) = a bar ( r ) (cid:18) + p + a / a bar (cid:19) , (9)and µ ( x ) = − exp (− x ) ; (10)Alternatively McGaugh et al. (2016) adopt a ( r ) = a bar ( r )/ ν ( a bar / a ) (11)with ν ( y ) = ( − exp (−√ y )) (12)and find a = . ± . (ran.) ± . (sys.) × − m s − .Fig. 1 shows these interpolation functions. On a log-arithmic scale over 4 orders of magnitude the differencesappear small (left panel). Upon closer inspection the differ-ences are significant (right panel). At a bar = a the massdiscrepancy expressed as the dark matter fraction variesfrom . ∼ < f DM ∼ < . . At a bar = a , the mass discrep-ancy varies from . ∼ < f DM ∼ < . . These differences high-light the non-uniqueness of predictions from MOND. Therelation between total acceleration and baryonic accelera-tion (left panel) is often referred to as the radial accelera-tion relation (RAR), while the relation between the massdiscrepancy and baryonic acceleration (right panel) is of-ten referred to as the mass discrepancy acceleration relation(MDAR, Sanders 1990; McGaugh 2004).There have been numerous theoretical studies of theRAR in a Λ CDM context, using both analytic mod-els and cosmological hydrodynamical simulations (e.g.,van den Bosch & Dalcanton 2000; Di Cintio & Lelli 2016;Desmond 2017; Keller & Wadsley 2017; Ludlow et al. 2017;Navarro et al. 2017). These authors have shown that galaxyformation in a Λ CDM universe results in a RAR close to, butnot identically equal to that commonly ascribed to MOND.
MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM In this paper we use the NIHAO suite of cosmologi-cal galaxy formation simulations (Wang et al. 2015) to gaininsight into the origin of the RAR in a Λ CDM universe.Compared to previous cosmological hydrodynamical simula-tions, NIHAO has more galaxies with higher resolution overa wider range of galaxy stellar masses ∼ < M star ∼ < M ⊙ .In § § § § § As our observational sample we use the Spitzer Photom-etry and Accurate Rotation Curves (SPARC) database(Lelli et al. 2016). SPARC is the largest sample of galaxyrotation curves with spatially resolved data on the distribu-tion of both stars and gas. The full sample of 175 galaxiesspans the full range of stellar masses of rotationally sup-ported galaxies ( ∼ < M star ∼ < M ⊙ ). It includes near-infrared ( . µ m) observations to trace the distribution ofstellar mass and 21 cm observations that trace the atomichydrogen gas ( Hi ). The 21 cm data also provide velocityfields from which rotation curves are derived. In some casesthe rotation curves are supplemented with higher spatial res-olution rotation curves from ionized gas. All rotation curvesare nominally corrected for inclination.Following McGaugh et al. (2016) we apply some qualitycriteria. Ten galaxies with inclination i < are removed tominimize sin ( i ) corrections to the observed velocities. Twelvegalaxies with asymmetric rotation curves are rejected (theseare flagged with Q = ). This leaves a sample of 153 galaxies.Additionally, only data points with relative velocity errorsless than 10% are kept.Fig. 2 shows histograms of the global stellar and atomichydrogen gas properties of the SPARC galaxies (blue lines),and NIHAO simulations (red shaded, see § M star ; half-light radius at 3.6 micron, R star ; atomic hydro-gen mass, M HI ; and radius at which Hi density equals 1[ M ⊙ pc − ], R HI . The upper left panel of Fig. 3 shows therotation curves of this sample, where the color code corre-sponds to the stellar mass. The boundaries of the mass binsare log ( M star / M ⊙ ) = . , . , . , . , . . The lower panelshows average rotation curves in bins of stellar mass. Foreach mass bin we interpolate each rotation curve onto a ra-dial grid, then we average in log ( V ) . We also calculate theaverage minimum and maximum radius, which defines therange over which the average rotation curves are plotted.The mean mass is given to the right of each average rota-tion curve. We see the shape and extent of the rotation curvesystematically changes from low to high mass.Here we briefly discuss uncertainties in the observations.(i) Conversion from stellar light to stellar mass. The nom-inal value is assumed to be M / L . = . for the disk and . for the bulge (McGaugh et al. 2016), but could plausi-bly vary from 0.2 to 1.0 depending on initial mass function(IMF), age and metallicity (McGaugh & Schombert 2015). Figure 2.
Histograms of masses and sizes for atomic hydrogenand stars from NIHAO simulations (red, shaded), and SPARCobservations (blue). For SPARC the stellar masses and projectedhalf-light radii are obtained from . µ m photometry, while forNIHAO they are calculated from the stellar particles. (ii) Distance. Since these galaxies are nearby, distance er-rors can be significant (up to 30%). Only a subset of theSPARC sample has accurate (5%) distances from e.g., usingthe Tip of the Red Giant Branch (TRGB). Since physicalsize scales as D , and luminosity as D , the acceleration ofthe baryons ( ∼ M bar / r ) is independent of distance. Theobserved rotation velocity is also independent of distance(ignoring beam smearing which is worse for more distantgalaxies), so that the total acceleration depends on distanceas D − .(iii) Rotation curves. The conversion from observed ro-tation to circular velocity depends on inclination, pressuresupport and warps (especially at large radii).(iv) Gas budget. The majority of gas in the galaxy is inatomic hydrogen, which is observed for SPARC galaxies, butthere is also a small amount of molecular and ionized gas,which is not easily observed, and generally missing from theSPARC data. For theoretical predictions for the RAR in a Λ CDM uni-verse we use the NIHAO galaxy formation simulations(Wang et al. 2015). These are a sample of ∼ cosmolog-ical hydrodynamical galaxy formation simulations run withthe SPH code gasoline2 (Wadsley et al. 2017).Haloes are selected at redshift z = from parent dissipa-tionless simulations of size 60, 20, & 15 Mpc / h , presented inDutton & Macci`o (2014) which adopt a flat Λ CDM cosmol-ogy with parameters from the Planck Collaboration et al.(2014): Hubble parameter H = 67.1 km s − Mpc − , matterdensity Ω m = . , dark energy density Ω Λ = − Ω m = MNRAS , 1–15 (2019)
Dutton et al.
Figure 3.
Rotation curves from SPARC observations (left) and circular velocity profiles in the plane of the disk from NIHAO simulations(right). Upper panels show all the individual profiles color coded by the stellar mass. The boundaries of the mass bins are log ( M star / M ⊙ ) = . , . , . , . , . , . , and are chosen to have equal numbers of simulated galaxies in each bin. For NIHAO these are shown from twicethe force softening of the dark matter particles to the Hi radius (which encloses 90% of the Hi mass). Lower panels show mean profilesin bins of stellar mass, and are plotted between the average minimum and maximum radii of the profiles. The mean stellar masses areshown to the right of each line (this can be slightly different between NIHAO and SPARC due to the different distribution of galaxieswithin each mass bin). . , baryon density Ω b = . , power spectrum nor-malization σ = . , power spectrum slope n = . .Haloes are selected uniformly in log halo mass from ∼ to ∼ without reference to the halo merger history, con-centration or spin parameter. Star formation and feedbackis implemented as described in Stinson et al. (2006, 2013).Mass and force softening are chosen to resolve the mass pro-file of the target halo at ∼ < the virial radius, which resultsin ∼ dark matter particles inside the virial radius of allhaloes at z = . The motivation of this choice is to ensurethat the simulations resolve the galaxy dynamics on the scale of the half-light radii, which are typically ∼ . of the virialradius (Kravtsov 2013).Each hydro simulation has a corresponding dark mat-ter only (DMO) simulation of the same resolution. Thesesimulations have been started using the identical initial con-ditions, replacing baryonic particles with dark matter par-ticles. The full sample of hydro simulations we use is 89.When comparing hydro to DMO we remove 5 simulationsfor which either the hydro or DMO simulation is undergoinga major merger and is thus out of equilibrium.Haloes in NIHAO zoom-in simulations were identified MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM Figure 4.
Total acceleration, a , vs acceleration due to baryons, a bar , in NIHAO simulations (left red/orange) and SPARC observations(right, blue/cyan) compared to MOND (black solid line) with an acceleration scale a = . × − m s − . Small points show the individualdata points. Large circles with error bars show the mean and σ scatter in bins of a bar . The lines show the spline interpolated meanrelations. The inset panel shows the residuals (orange/cyan histogram) relative to the interpolated mean relations, which has a standarddeviation of just 0.079 dex in NIHAO and 0.132 dex in SPARC. The corresponding Gaussian is shown in red (for NIHAO) and blue (forSPARC). The range of galaxy stellar masses, M star , is as indicated, while N is the number of data points. using the MPI+OpenMP hybrid halo finder AHF (Gill et al.2004; Knollmann & Knebe 2009). AHF locates local over-densities in an adaptively smoothed density field as prospec-tive halo centers. In this study we use one halo from eachzoom-in simulation (the one with the most particles). Thevirial masses of the haloes are defined as the masses within asphere whose average density is 200 times the cosmic criticalmatter density, ρ crit = H / π G . The virial mass, size andcircular velocity of the hydro simulations are denoted: M , R , V . The corresponding properties for the dark matteronly simulations are denoted with a superscript, DMO . Forthe baryons we calculate masses enclosed within spheres ofradius r gal = . R , which corresponds to ∼ to ∼ kpc.The stellar mass inside r gal is M star , the atomic hydrogen, Hi , inside r gal is computed following Rahmati et al. (2013)as described in Gutcke et al. (2017).The NIHAO simulations are the largest set of cosmo-logical zoom-ins covering the halo mass range ∼ < M ∼ < M ⊙ . Their uniqueness is in the combination of high spa-tial resolution coupled to a statistical sample of haloes. Inthe context of Λ CDM they form the “right” amount of starsboth today and at earlier times (Wang et al. 2015). Theircold gas masses and sizes are consistent with observations(Stinson et al. 2015; Macci`o et al. 2016; Dutton et al. 2019),they follow the gas, stellar, and baryonic Tully-Fisher rela-tions (Dutton et al. 2017). On the scale of dwarf galaxies thedark matter haloes expand yielding cored dark matter den-sity profiles consistent with observations (Tollet et al. 2016),and resolve the too-big-to-fail problem of field galaxies(Dutton et al. 2016). They reproduce the diversity of dwarf http://popia.ft.uam.es/AMIGA galaxy rotation curve shapes Santos-Santos et al. (2018),and the Hi linewidth velocity function (Macci`o et al. 2016;Dutton et al. 2019). They match the observed clumpy mor-phology of galaxies seen in CANDELS (Buck et al. 2017).As such they provide a good template with which to studythe RAR in a Λ CDM context.Histograms of the stellar and atomic hydrogen massesand sizes are shown with red shading in Fig. 2. Broadlyspeaking NIHAO and SPARC have similar distributions, butNIHAO extends to lower masses. The right panels of Fig. 3shows the circular velocity profiles from the NIHAO galaxyformation simulations. These have been calculated using thepotential in the plane of the disk using Eq.1. The NIHAOsimulations have been plotted from two dark matter forcesoftening lengths ( ≈ the convergence radius), to the radiusenclosing 90% of the HI ( ≈ the edge of the observable galac-tic HI). The radial range of the velocity profiles is similar,but slightly narrower. The change in velocity profile shapesfrom low to high mass is similar between simulations andobservations. In detail there are differences, caused at leastin part by differences in the distribution of baryons. Fig. 4 shows the relation between total acceleration, a , andthe acceleration due to the baryons, a bar for NIHAO galaxies(left) and SPARC observations (right). The dots show allof the individual measurements for observed and simulatedgalaxies. The total number of data points, N ∼ forNIHAO and N ∼ for SPARC. For NIHAO we samplethe circular velocity profiles linearly in units of the virialradius, but we only include points which are observable (i.e., MNRAS , 1–15 (2019)
Dutton et al.
Figure 5.
Scatter in the RAR. The total observed scatter is . (grey dashed vertical line). The observational errors contribute . ± . (blue dotted histogram). Subtracting the observa-tional errors gives the intrinsic scatter (red histogram). For caseswhere the observed error is greater than the total scatter we set σ = , which is the spike. The intrinsic scatter has a median of0.06 and a 90% confidence interval of 0.00 to 0.09. within the HI radius). This is similar to the observations,which tend to sample the rotation curves linearly in radius,and with similar numbers of data points for different massgalaxies. The large points with error bars show the mean of log ( a ) in bins of log ( a bar ) , while the error bar shows thescatter about the spline interpolated mean relation.The solid black line shows the relation fitted to SPARCdata by McGaugh et al. (2016). The functional form waschosen to follow the asymptotic relations for high and lowaccelerations in MOND (dotted lines). Overall the NIHAOsimulations provide a very good match to the observed RAR.The agreement is especially good when one considers thatthere are observational systematics, such as the stellar mass-to-light ratio and distance scale that can shift the observedrelation (McGaugh et al. 2016; Lelli et al. 2017). Further-more, the NIHAO simulations are necessarily an approxi-mation for how galaxies form in a Λ CDM universe.The inset panels show the scatter about the interpolatedmean relation (histogram). The standard deviation and cor-responding Gaussian is given in red (for NIHAO) and blue(for SPARC). The scatter from NIHAO is just 0.079 dex.Similar small scatters were reported by other Λ CDM simu-lation studies (Keller & Wadsley 2017; Ludlow et al. 2017).This is less than the total observed scatter of 0.132 dex.However, there are observational errors which means the in-trinsic scatter is lower. Note one should be aware that inprinciple observational errors can correlate with offsets fromthe RAR leading to lower observed than intrinsic scatter.For example if galaxies with higher stellar M / L have lower a , this would result in the RAR measured with constant M / L having less scatter than the true RAR. Figure 6.
Mass dependence of scatter in the RAR. Red circlesshow results from NIHAO simulations. Blue squares show the to-tal and intrinsic scatter from SPARC galaxies. The vertical errorbars show the intrinsic scatter from SPARC adopting measure-ment errors of 0.09 and 0.14 dex, for the upper and lower limit,respectively. The horizontal lines show the range of M star in eachbin. As discussed above, observational sources of error includedistance, stellar mass-to-light ratio, disk inclination, ro-tation velocity, and atomic hydrogen flux. According toMcGaugh et al. (2016) these contribute 0.08, 0.06, 0.05,0.03, 0.01 dex to the scatter in a at a bar , yielding a totalobservational error of 0.116 dex. Subtracting this (in quadra-ture) from the total observed scatter of 0.132 dex yields anintrinsic scatter of 0.063 dex.However, it should be noted that the observational er-rors are just estimates. The true observational errors couldbe larger or smaller. For example, McGaugh et al. (2016)adopt an error on the stellar mass-to-light ratio of 0.11 dex.The true error could plausibly be anywhere in the range 0.05to 0.20 dex, and it could depend systematically on the galaxymass, or other galaxy property. To get an idea of the plausi-ble range of the intrinsic scatter we assume the observationalerrors are drawn from a Gaussian with the mean specifiedabove, and a standard deviation of 20% of the mean. Forexample, if the distance errors are reported to be are 0.08dex, we adopt a standard deviation of 0.016 dex.Assuming the uncertainties in the errors of the fivesources are uncorrelated, using a Monte Carlo simulation thetotal observational error is . ± . (Fig. 5). Subtract-ing this from the total observed scatter of 0.132 dex yieldsan intrinsic scatter with a wide variation (red histogram).The 90% confidence interval ranges from zero to . dex.Thus the scatter in the RAR from our simulations of 0.079dex is consistent with the observations. The prediction fromMOND for zero intrinsic scatter is also consistent with theobservations. In order to use the scatter in the RAR to dis-tinguish between Λ CDM and MOND, either a more accurate
MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM Figure 7.
As Fig. 4 but for galaxies in narrow stellar mass ranges as indicated at the top right of each panel. At these (high) massesthe simulated RAR is very similar to the observed one, and has very small scatter making the individual data points are hard to see. measurement is required, or as we shall see below we needto look at the RAR scatter for different mass galaxies.The largest source of observational uncertainty are thedistances, which account for half of the total variance. Soif these can be reduced a significant improvement in the in-trinsic scatter is possible. The galaxies with the largest dis-tance errors from the SPARC survey use the Hubble flow, D = v sys / H , so the unknown peculiar velocities provide arandom source of error, while the value of the Hubble con-stant provides a systematic error. The peculiar velocity errorcan be reduced by simply observing galaxies at larger dis-tances. This is easier said than done, because galaxies thatare further away have smaller angular sizes. It is possible toobtain rotation curves using optical emission lines for suchgalaxies (e.g., Courteau et al. 2007), but atomic hydrogengas maps will require future radio telescopes such as the Square Kilometer Array (SKA). These will have the resolu-tion and sensitivity to measure HI density maps and rota-tion curves of galaxies at large enough distances such thatpeculiar velocity errors are negligible. The dependence of the scatter in the RAR on stellar mass isshown in Fig. 6. In both simulations and observations highermass galaxies have smaller scatter than lower mass galaxies.In the highest mass bin centered on M star ≃ . M ⊙ , thescatter in NIHAO is just 0.036 dex (0.045 dex relative toMOND), compared with 0.110 dex for SPARC. In the massbin centered on M star ≃ . M ⊙ the scatter in NIHAO in-creases to 0.10 dex (0.14 dex relative to MOND), comparedwith 0.172 dex for SPARC. MNRAS , 1–15 (2019)
Dutton et al.
Figure 8.
As Fig. 4 but for galaxies in narrow stellar mass ranges as indicated at the top right of each panel. At these (low) massesthe simulated RAR shows small but significant deviations from the observed one. Both simulations and observations have larger scatterthan at higher masses.
The quadratic differences between the simulated andobserved scatter from high to low mass are 0.104, 0.087,0.138, and 0.140 dex. In other words, these are the size ofthe measurement errors needed to reconcile the observedscatter with the simulated scatter. Based on Fig. 5 and thereported measurement errors for SPARC galaxies, whichdo not vary significantly with galaxy mass, these are rea-sonable numbers. So we conclude that the mass depen-dent scatter in the NIHAO RAR is consistent with the in-trinsic scatter from observations. The zero scatter assumedby MOND is also consistent with the high mass galaxies . ∼ < M star ∼ < . M ⊙ , but inconsistent with the lowmass galaxies . ∼ < M star ∼ < . M ⊙ . Thus using the samedata that has been used as evidence in favor of MOND byMcGaugh et al. (2016), we have shown that the intrinsic scatter in the RAR actually disfavors MOND, and is in ex-cellent agreement with predictions from Λ CDM . Our resultechos that of Rodrigues et al. (2018) who concluded that theobserved rotation curves from SPARC could not be fittedwith MOND assuming a universal value of the accelerationscale, a .Figs. 7 & 8 shows the RAR divided into four bins ofstellar mass centered on M star ≃ . , . , . , and . M ⊙ . These bins are the same as used in Fig. 3, andhave been chosen to have roughly the same number of NI-HAO galaxies per bin.A feature that stands out in both simulations and ob-servations is that higher mass galaxies span a wider rangeof accelerations than lower mass galaxies. This is due to thesystematic change in the shape of the circular velocity pro- MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM Figure 9.
Total acceleration, a , vs acceleration due to baryons, a bar , for the lowest stellar mass galaxies in NIHAO. Black cir-cles show the acceleration at the projected stellar half-mass radii.These galaxies deviate significantly from the MOND RAR, butare consistent with the observed relation for dwarf spheroidalgalaxies (green line, Lelli et al. 2017). files with galaxy mass (see Fig. 3). The highest mass galaxieshave velocity profiles close to constant, and thus accelerationvaries inversely with radius, whereas low mass galaxies haverising velocity profiles with slopes close to 0.5, which resultsin acceleration independent of radius.For the two highest mass bins . ∼ < M star ∼ < . M ⊙ (Fig. 7) the mean NIHAO RAR follows the MOND RARextremely closely. For the two lowest mass bins . ∼ < M star ∼ < . M ⊙ (Fig. 8) the mean NIHAO RAR showsdepartures from the MOND RAR. At baryon accelerationsof a bar ∼ − ms − MOND over-predicts NIHAO, while at a bar ∼ < − ms − MOND under-predicts NIHAO. This trendcontinues when we look at even lower mass galaxies.Fig. 9 shows the RAR for the lowest mass galaxies inNIHAO, with stellar masses between . M ⊙ and . M ⊙ .These galaxies deviate significantly from the MOND relationbelow an acceleration scale of a bar = − m s − . The greenline shows a fit to the observed RAR for dwarf spheroidalgalaxies from Lelli et al. (2017). While the observed and sim-ulated dwarfs have similar ranges of stellar masses, thereare some differences that should be considered. The observa-tions are primarily satellite galaxies, whereas the simulationsare all centrals. The observed total accelerations are basedon stellar kinematics, whereas for the simulations we showpoints sampled from the radii observable with atomic hydro-gen (using the same procedure as for more massive galaxies).The black circles in Fig. 9 show the RAR for NIHAO galax-ies at the location of the projected half-mass radii of thestars. These points tend to sample higher accelerations, butfollow the same trend as the full velocity curve points. Eventhough there are some caveats in the comparison betweenour Λ CDM based simulations and observations, the agree-ment is an encouraging sign.
Figure 10.
Dark matter acceleration, a dm , vs acceleration dueto baryons, a bar , for NIHAO. The points with lowest and highestaccelerations deviate the most from the MOND RAR (black line).The dotted lines show the 1:1 relation and the asymptotic MONDrelation at low a bar . The deviation of the observed RAR for dwarf Spheroidalgalaxies from the MOND prediction is another apparent fail-ure of MOND. The particular assumption that breaks is thatthe RAR has a slope of 1/2 at low baryon accelerations. Re-call, that this assumption was made to ensure that rotationcurves are asymptotically flat at large radii.
Observationally the two axes of the RAR are independent,however in simulations the two axes are correlated when thedark matter fraction is low, because a ≡ a bar + a dm . Fig. 10shows an alternative way to look at the RAR: the dark mat-ter acceleration vs the baryon acceleration. This is exactlythe same data as in Fig. 4, but with independent axes weclearly see that NIHAO and MOND diverge at high baryonaccelerations. The scatter in a dm | a bar is also larger than thatof a | a bar : 0.122 vs 0.079. Observationally it is not advisedto calculate the dark matter acceleration, as measurementerrors can result in negative values. However, it would bepossible to determine the scatter in a dm by forward mod-elling the observed RAR. Having established that the NIHAO simulations reproducethe slopes, normalization, and small scatter in the RAR,we now ask: where does the small scatter come from in oursimulations?Since the RAR is a prescription for how the total ac-celeration is related to the acceleration due to baryons, wecan apply the RAR to our simulated baryon circular ve-locity profiles (blue dashed lines in Fig. 11). The baryon
MNRAS , 1–15 (2019) Dutton et al.
Figure 11.
Examples of NIHAO circular velocity curves (red circles) where the MOND prediction (black lines) works well. The verticalgrey lines show twice the dark matter softening length (dotted) and the HI radius (dashed). The standard deviation of the velocityresiduals, σ V , and acceleration residuals, σ loga , is computed for data points between these two lines and is given in the top right cornerof each panel. These galaxies span over three orders of magnitude in stellar masses. circular velocity profiles are the (quadratic) sum of the cir-cular velocity due to the stars and gas, calculated from thegravitational potential in the plane of the disk.Fig. 11 shows examples of four simulations for whichthe MOND prescription (fitted to the SPARC data) accu-rately predicts the circular velocity profile of the simulation.The standard deviation of V and log a is given in the topright corner. The agreement, to within − , is remark-able! We stress that there are no free parameters. The stellarmasses of these examples range over more than three ordersof magnitude from . to . M ⊙ .However, we have chosen these 4 examples out of a sam-ple of 89 examples. Fig. 12 shows examples of two simula-tions for which the MOND prescription clearly fails, one where MOND under-predicts the velocity (g3.61e11, up-per left) and another where MOND over-predicts (g6.77e10,lower left). For g3.61e11 an excellent fit can be obtainedsimply by rescaling the stellar mass profile by a factor of 1.3(i.e., 0.11 dex), which is a typical uncertainty on an observedstellar mass-to-light ratio. For g6.77e10 changing the stellarmass normalization has little impact because the baryons aredominated by gas. In this case a good fit can be obtainedby reducing the distance by 40%. This would be a σ un-certainty in distances. One can also change the accelerationscale, a , to obtain a similar effect.Fig. 13 shows the standard deviation of MOND fits toNIHAO circular velocity profiles. The fiducial (no free pa-rameter) fit is shown with red open circles. The median error MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM Figure 12.
Examples of NIHAO circular velocity curves where MOND fails. Top left: MOND under predicts the velocity. An excellentfit can be obtained by increasing the stellar mass by 30% (top right), mimicking the observational uncertainty in the stellar mass-to-lightratio. Bottom left: MOND over predicts the velocity. A good fit can be obtained by reducing the distance by 40% (bottom right). is just − . If we allow the stellar mass and distance tovary, we can obtain even better fits (black circles), with amedian error of just . − . This analysis is of relevanceto the study of Li et al. (2018) who are able to fit the SPARCsample with a single RAR, allowing for variation in M / L ofthe disk and bulge, distance and inclination. At first thisseems an impressive result for MOND. However, our similaranalysis of NIHAO galaxies shows that MOND fits have toomuch freedom, and can fit galaxy velocity profiles that theyshould not be able to.These figures also highlight some of the issues with mea-suring the scatter in the RAR. Because the scatter is mea-sured in logarithmic space, galaxies with flat rotation curvestend to have smaller scatter than those with rising rotationcurves where errors at small velocities get magnified. The other issue is spatial sampling. Here we sample galaxies uni-formly in units of the virial radius, but we only include pointsthat are observable (i.e., within the HI radius). Thus for twogalaxies that live in similar mass haloes but have differentHI sizes, the larger galaxy will have more points in the RAR.This implies that galaxies effectively have different weightsin the RAR. Given the current large uncertainties in the in-trinsic scatter this is a detail, but in the future it will benecessary to sample the rotation curves in observations andtheory in a consistent way (see e.g., Desmond 2017). Since the RAR is a prescription of predicting the total ac-celeration given the baryon acceleration, it is thus a pre-
MNRAS , 1–15 (2019) Dutton et al.
Figure 13.
Residuals of MOND fits to NIHAO circular velocitycurves. Fiducial fits (red open circles) which have no free param-eters, result in a median scatter of just 6 km/s. When the stellarmass and “distance” are fitted for (black filled circles) residualsare reduced even further to 2.5 km/s. scription for predicting the mass discrepancy or in otherwords the dark matter acceleration. When the RAR workswell in NIHAO this means that RAR is accurately predict-ing the dark matter profile. By construction, for the galax-ies where the RAR under-predicts the circular velocity, theRAR under-predicts the dark matter halo, and for the galax-ies where the RAR over-predicts the circular velocity, RARover-predicts the dark matter halo. It is then interesting tocompare the dark matter profiles from the RAR with thosefrom the NIHAO and DMO simulations. This comparison ismade in Fig. 14.Each panel shows a different stellar mass range (corre-sponding to those used in Figs. 7 - 9), except we have splitthe highest mass bin into two to make clearer the differentdark matter profiles at the highest masses we simulate. Themean log ( M star ) is indicated in the top left of each panel.The thick lines show mean (of log V ) dark matter circularvelocity profiles for NIHAO (red solid), DMO (cyan longdashed) and NIHAO MOND (black short-dashed). The thinlines show the standard deviation (of log V ).One of the striking differences between Λ CDM (bothDMO and NIHAO) and MOND is that at large radii (i.e.,close to the virial radius) MOND always over-predicts thecircular velocity. This is because of the assumption that a ∝ a bar at low accelerations. If all the baryons were near thecenter of the halo, the velocity profile would indeed becomeflat (which would not be as bad, but still not in agreementwith the declining velocity profiles of the simulations), but itrises at large radii due to the gas in the halo. Recall that inour comparison to SPARC data we only considered radii thatcan be traced with atomic hydrogen gas, which typically onlyextends up to 5-10% of the virial radius (see the red circles inFig. 14 which show the average Hi radius). Restricting to the observable radii, the MOND dark matter profiles typicallyfall between the DMO and NIHAO lines.At stellar masses between ∼ < M star ∼ < M ⊙ theNIHAO simulations have resulted in halo expansion to aconstant density (compare with the reference lines for ρ ∝ r − α with α = and α = ), see also Tollet et al. (2016).MOND also tends to predicts expanded haloes, but not asmuch as in NIHAO.For the highest masses ( M star = . M ⊙ , top left panel)the NIHAO simulations have contracted dark haloes, whileMOND predicts expansion. These are examples similar tothat shown in the top left of Fig. 12, where this deficit indark matter can be compensated for with increased stellarmass. In spite of these differences, because the centers ofthese galaxies are dominated by baryons, they are only off-set by a small amount from the SPARC RAR. The averagerms in a of these 8 galaxies is 0.056 dex, making them typicalgalaxies. Determining the dark matter profiles observation-ally, and thus distinguishing between the predictions fromMOND and NIHAO requires an accurate determination ofthe stellar mass-to-light ratio, which is a challenge.A summary of the density profiles is shown in Fig. 15.This shows the dark matter circular velocity at 1% of thevirial radius vs stellar mass (left panel) and the slope of thedark matter circular velocity between 1 and 2% of the virialradius [ γ V = log ( V . / V . )/ log ( ) ] vs stellar mass (rightpanel). Both plots show the same qualitative trends. DMOsimulations (cyan) have structure that is almost scale free,with velocity log [ V dark ( . R )/ V ] ≃ − . and velocityslope γ V ≃ . . Note that we choose to show velocity sloperather than (local) density slope because the latter involvesan extra derivative of the velocity profile, and so is noisierin both observations and simulations. The circular velocityprofile, or equivalently the cumulative mass profile, or cu-mulative enclosed mass density profile, also has advantagesin terms of analytic models (Dekel et al. 2017). The velocityslope can be easily converted into an enclosed mass densityslope using γ ρ = γ V − . So for example a constant den-sity core, γ ρ = corresponds to a velocity slope of γ V = ,a NFW density slope of γ ρ = − corresponds to a velocityslope of γ V = . , and an isothermal density slope of γ ρ = − corresponds to a velocity slope of γ V = .In NIHAO hydro simulations (red points and lines)the halo response is strongly mass dependent (See alsoTollet et al. 2016; Dutton et al. 2016). At low stellar masses M star ∼ < M ⊙ there is very little halo response. As the stel-lar mass increases the halo expands (i.e., lower velocity andhigher velocity slope), reaching a maximum expansion atstellar masses of M star ∼ M ⊙ , where the velocities are ≃ . dex lower and the velocity slopes are 0.4 higher thanDMO and close to a constant density core ( γ V = ). Bya mass of M star ∼ . M ⊙ the haloes are on average un-changed at small radii, and at the highest stellar masses M star ∼ M ⊙ the haloes contract at small radii. TheMOND dark matter haloes (black points and lines) typicallyshow behavior intermediate between NIHAO and DMO (i.e.,the haloes are not as cuspy and DMO, and not as expandedas NIHAO). The exception is at the highest masses, whereMOND haloes have lower velocities and higher slopes, theopposite of NIHAO, which has higher velocities and lowerslopes. MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM Figure 14.
Dark matter circular velocity profiles from simulations. Dark matter only simulation (long-dashed cyan), hydrodynamical(red solid), and MOND prediction using the baryon circular velocity from the hydrodynamical simulation (short-dashed black). Lines areplotted between twice the dark matter softening to the virial radius. The red dot shows the average Hi radius (enclosing 90% of the Hi mass). Each panel shows galaxies with a different range of stellar mass, the average stellar mass h log ( M star / M ⊙ )i is indicated in the topleft of each panel. The standard deviation between the NIHAO and MOND dark matter velocity profiles is given by σ . For reference weshow lines with density slopes α = and α = , where ρ ∝ r − α . One can see systematic deviations between MOND and hydro simulations.MNRAS , 1–15 (2019) Dutton et al.
Figure 15.
Dark matter circular velocity at 1% of the virial radius (left) and dark matter circular velocity slope between 1-2% of thevirial radius (right) for NIHAO hydro (red solid), DMO (cyan long-dashed), and MOND (black short-dashed). The right axis of the rightpanel shows the enclosed DM density slope. Stellar mass bins are the six panels shown in Fig. 14. DMO simulations have a velocity andvelocity slope independent of the stellar mass. NIHAO hydro simulations have contraction at the highest stellar masses, and maximumexpansion for stellar masses of M star ∼ M ⊙ . MOND haloes are typically intermediate between DMO and Hydro, except for the higheststellar mass, where MOND predicts expansion to a core while NIHAO predicts a contracted cusp. We use 89 cosmological galaxy formation simulations fromthe NIHAO project (Wang et al. 2015) to study the origin ofthe radial acceleration relation (RAR) in a Λ CDM universe.The unique combination of halo mass range, resolution, andnumber of haloes allows us to compare with observations ofthe RAR from the SPARC survey (McGaugh et al. 2016).The NIHAO galaxies have stellar masses spanning . ∼ < M star ∼ < . M ⊙ .We summarize our results as follows: • The RAR exists in the NIHAO simulations with a sim-ilar slope and normalization as observationally determinedby SPARC (Fig. 4). • The RAR in NIHAO has a scatter of 0.079 dex. This isconsistent with estimates of the intrinsic scatter from obser-vations ( . ∼ < σ int ∼ < . ) (Fig. 5). • The scatter in the RAR depends on the stellar massof the galaxy, with lower scatter in higher mass galaxies inboth observations and simulations (Fig. 6). For high massgalaxies ( M star ∼ > . M ⊙ ) the intrinsic scatter ( σ int ∼ < . )is consistent with both MOND and NIHAO. However, forlow mass galaxies ( M star ∼ < . M ⊙ ) the intrinsic scatter( . ∼ < σ int ∼ < . ) is consistent with NIHAO, but inconsis-tent with MOND. • The mass dependence of the RAR scatter in simulationsis partially explained by the correlation between a ≡ a bar + a dm and a bar at low dark matter fractions. • In the lowest mass galaxies we simulate ( M star ∼ M ⊙ )The RAR deviates significantly from the MOND prediction,but is in agreement with the observed RAR from dwarfspheroidal galaxies (Fig. 9). • We can use the RAR to accurately predict (to within a few km s − ) the circular velocity profiles of individual NI-HAO simulations by just using the baryon circular velocityprofile as input (Fig. 11). • In cases where the fiducial values provide a bad fit, thestellar masses can be re-scaled, and/or the distances changedto provide excellent fits (Fig. 12). Given realistic uncertain-ties in baryon profiles and galaxy distances, MOND has toomuch flexibility in fitting individual galaxy rotation curves.Thus the success of MOND at fitting individual galaxy ro-tation curves is not as impressive as it is often claimed to be(e.g., Sanders & McGaugh 2002; Li et al. 2018). • The RAR predicts dark matter circular velocity profilesthat are on average similar to that found in Λ CDM simula-tions (Fig. 14). • In detail, the RAR generally predicts mild halo expan-sion at small radii (relative to dissipationless CDM), butnot as much expansion as found in the NIHAO simulations(Fig. 15). This contradicts the claim by Navarro et al. (2017)that explaining the RAR in CDM does not require modifi-cations to the cuspy inner mass profiles of dark haloes. • The largest differences in the dark matter profiles are inthe highest mass galaxies we simulate. For ( M star ∼ M ⊙ )NIHAO predicts halo contraction, while the RAR predictsexpansion to a dark matter core (Figs. 14 & 15).It may seem paradoxical that the MOND RAR appearsto be in the NIHAO Λ CDM galaxy formation simulations.The resolution is both MOND and NIHAO are approximat-ing the same thing (the observable Universe).Since the phenomenological basis for MOND is natu-rally reproduced by Λ CDM galaxy formation simulations,it seems unnecessary to invoke the simplifying assumptionsand radical new force law required by MOND. A useful anal-ogy is with solar system. A phenomenological theory might
MNRAS , 1–15 (2019) rigin of MOND phenomenology in LCDM assume that the orbits of planets are circles, since they areobserved to be roughly circles, and this is the simplest ge-ometry. But we know that with sufficiently accurate obser-vations the orbits are actually ellipses, and the circle modelis an over-simplification.The assumptions MOND is based on are only approx-imately true in Λ CDM : 1) At high accelerations the darkmatter fraction is low, but non-zero; 2) Circular velocityprofiles are only approximately flat at the HI radius of agalaxy. At larger radii the profiles decline; 3) The scatterin the RAR is small, but non-zero, primarily reflecting non-universality of the dark matter density profiles. Finally, wenote that if one wishes to go beyond Newtonian dynamicsthere is more freedom in the observed RAR than imposedby MOND.
ACKNOWLEDGMENTS theo cluster of the Max-Planck-Institut f¨ur As-tronomie; and on the hydra cluster at the Rechenzentrumin Garching. We greatly appreciate the contributions of allthese computing allocations. AO is funded by the DeutscheForschungsgemeinschaft (DFG, German Research Founda-tion) – MO 2979/1-1. TB was supported by the Sonder-forschungsbereich SFB 881 “The Milky Way System” (sub-projects A1 and A2) of the DFG. The analysis made use ofthe pynbody package (Pontzen et al. 2013).
REFERENCES
Bond, J. R., Efstathiou, G., & Silk, J. 1980, Physical ReviewLetters, 45, 1980Buck, T., Macci`o, A. V., Obreja, A., et al. 2017, MNRAS, 468,3628Courteau, S., Dutton, A. A., van den Bosch, F. C., MacArthur,L. A., Dekel, A., McIntosh, D. H., & Dale, D. A. 2007, ApJ,671, 203Dekel, A., Ishai, G., Dutton, A. A., & Maccio, A. V. 2017, MN-RAS, 468, 1005Desmond, H. 2017, MNRAS, 464, 4160Di Cintio, A., & Lelli, F. 2016, MNRAS, 456, L127Dutton, A. A., & Macci`o, A. V. 2014, MNRAS, 441, 3359Dutton, A. A., Macci`o, A. V., Frings, J., et al. 2016a, MNRAS,457, L74Dutton, A. A., Obreja, A., Wang, L., et al. 2017, MNRAS, 467,4937Dutton, A. A., Obreja, A., & Macci`o, A. V. 2019, MNRAS, 482,5606Faber, S. M., & Gallagher, J. S. 1979, ARA&A, 17, 135Gill, S. P. D., Knebe, A., & Gibson, B. K. 2004, MNRAS, 351,399Gutcke, T. A., Stinson, G. S., Macci`o, A. V., Wang, L., & Dutton,A. A. 2017, MNRAS, 464, 2796Keller, B. W., & Wadsley, J. W. 2017, ApJ, 835, L17Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608 Kravtsov, A. V. 2013, ApJ, 764, L31Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S.2017, ApJ, 836, 152Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615,A3Ludlow, A. D., Ben´ıtez-Llambay, A., Schaller, M., et al. 2017,Physical Review Letters, 118, 161103Macci`o, A. V., Udrescu, S. M., Dutton, A. A., et al. 2016, MN-RAS, 463, L69McGaugh, S. S. 2004, ApJ, 609, 652McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, PhysicalReview Letters, 117, 201101Milgrom, M. 1983, ApJ, 270, 365Navarro, J. F., Ben´ıtez-Llambay, A., Fattahi, A., et al. 2017, MN-RAS, 471, 1841Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014,A&A, 571, A16Pontzen, A., Roˇskar, R., Stinson, G., & Woods, R. 2013, Astro-physics Source Code Library, 1305.002Rahmati, A., Schaye, J., Pawlik, A. H., & Raiˇcevi ` c , M. 2013,MNRAS, 431, 2261Rodrigues, D. C., Marra, V., del Popolo, A., & Davari, Z. 2018,Nature Astronomy, 2, 668Sanders, R. H. 1990, A&ARv, 2, 1Sanders, R. H., & McGaugh, S. S. 2002, ARA&A, 40, 263Santos-Santos, I. M., Di Cintio, A., Brook, C. B., et al. 2018,MNRAS, 473, 4392Smoot, G. F., Bennett, C. L., Kogut, A., et al. 1992, ApJ, 396,L1Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074Stinson, G. S., Brook, C., Macci`o, A. V., et al. 2013, MNRAS,428, 129Stinson, G. S., Dutton, A. A., Wang, L., et al. 2015, MNRAS,454, 1105Tollet, E., Macci`o, A. V., Dutton, A. A., et al. 2016, MNRAS,456, 3542van den Bosch, F. C., & Dalcanton, J. J. 2000, ApJ, 534, 146Wadsley, J. W., Keller, B. W., & Quinn, T. R. 2017, MNRAS,471, 2357Wang, L., Dutton, A. A., Stinson, G. S., Macci`o, A. V., Penzo,C., Kang, X., Keller, B. W., Wadsley, J. 2015, MNRAS, 454,83Zwicky, F. 1933, Helvetica Physica Acta, 6, 110This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000