Amir Ghari

Benoit Famaey

Chervin Laporte

Hosein Haghi

Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Keywords: galaxies, kinematics and dynamics -galaxies, spiral -dark matter

come   L'archive ouverte pluridisciplinaire

Introduction

The nature of the dark sector of the Universe arguably represents one of the deepest mystery of modern physics. Despite earlier hints (e.g., Zwicky 1933), it is the flattening of rotation curves (RCs) of disk galaxies which provided the first clean observational evidence for mass discrepancies in galactic systems (Bosma 1978;Rubin et al. 1978). Many di erent observational probes -especially on cosmological scales (e.g., Planck Collaboration et al. 2016)then led to the development of the current standard cosmological framework ( CDM), in which one needs to add non-baryonic dark matter (DM) particles to the baryonic content of the Universe, which can thus in principle explain the RC shapes.

However, while arguably very successful on large scales, the current CDM picture is nevertheless, at face value, facing some challenges (e.g., Bullock & Boylan-Kolchin 2017), notably to explain the shapes of RCs in some galaxies.

From a theoretical perspective, the standard paradigm proposes that galaxies form and evolve in virialised haloes (e.g., White & Rees 1978) of cold DM (CDM). Quantum fluctuations from the vacuum state create the seed perturbations in the early universe which are stretched by infla- † CITA National Fellow tion, and because DM forms a collisionless and pressureless fluid, it is then able to collapse under gravity, providing the first potential wells in which the baryons later cool into. The evolution of the DM fluid is modelled by numerical simulations (Davis et al. 1985). DM-only (DMO) cosmological simulations treat 100 percent of the matter in the Universe as collisionless DM. CDM haloes simulated with DMO numerical models have density profiles that are well described to first order by the NFW profile (Navarro et al. 1996b) with a central cusp (d lnfl/d lnr = ≠1 in the center, where fl is the DM density), while observations rather point to large constant density cores of DM in the central parts of a large number (but not all) of rotationally-supported galaxies (de Blok et al. 2001;Gentile et al. 2004). This problem of DMO simulations becomes even more stringent when realizing that galaxies residing in halos of similar maximum circular velocity display a wide range of RC shapes in the central parts (Oman et al. 2015). This diversity of rotation curve shapes is to be contrasted with the uniformity found in the relation between the baryonic gravitational acceleration and the total gravitational acceleration (Mc-Gaugh et al. 2016;Lelli et al. 2017;Li et al. 2018a;Desmond 2017a,b).

Despite the NFW profile being widely used, one has to note that Navarro et al. (2004) and Macciò et al. (2008) Article number, page 1 of 29 arXiv:1811.06554v2 [astro-ph.GA] 29 Jan 2019 have proposed another model that fits the density profiles of halos in DMO simulations better in the inner regions, a profile which had been previously introduced for the distribution of stellar light and mass in galaxies by Einasto (1965), and which does not harbour a central cusp but a slope d lnfl/d lnr à ≠r 1/n in the center, meaning that the slope is indeed zero very close to the nucleus. For halos of the order and below the typical mass of the Milky Way halo in DMO simulations, the Einasto index is n ≥ 6. Nevertheless, as first shown by Chemin et al. (2011), fitting such profiles to observed rotation curves often leads to values of n in disagreement with DMO simulations, implying much larger cores in observations than in simulations, as expected from the previously known core-cusp issue. Here we reproduce such a study by using the latest and most up-to-date database on galaxy rotation curves, the Spitzer Photometry and Accurate Rotation Curves (SPARC) database compiled by Lelli et al. (2016a).

Upon the completion of this manuscript, Li et al. (2018b) presented a study on scaling relations between DM halo parameters and disc galaxy luminosities from the SPARC database. The authors chose two DM density profile parametrisations for their exploration: the Einasto profile, like us, and general (-, -, ") models (Zhao 1996) with parameters motivated by empirical relations as functions of stellar-to-halo mass described in Di Cintio et al. (2014) who analysed results from hydrodynamical cosmological simulations of galaxies drawn from the MAGICC project (Brook et al. 2012;Stinson et al. 2013). Using these models they fit their rotations curves through a Markov Chain Monte Carlo method, taking lognormal priors on near-infrared stellar mass-to-light ratios, Gaussian priors on the inclinations and distances as well as CDM priors for the DM halos. The authors in particular showed that the density at the radius where the slope of the Einasto profiles becomes ≠2 displays no correlation with luminosity.

Our work, while sharing many similarities in data and methodology with Li et al. (2018b), takes a di erent turn instead in examining the correlations of inferred DM halo parameters with the extent of the baryonic components. Such scaling relations are highly relevant to constraining the cusp-core transformation process. We also choose to fix the baryonic galactic parameters in two extreme cases which we propose to compare: (i) the original SPARC parameters with constant mass-to-light ratios for bulges and disks (Lelli et al. 2016a;McGaugh et al. 2016), and (ii) the parameters for which galaxies follow a very tight radial acceleration relation (RAR, Li et al. 2018a). In Sec. 2 we briefly present the SPARC sample, while we present the form of the Einasto profile and our fitting procedure in Sec. 3. We then present our results on scaling relations in Sec. 4 and conclude in Sec. 5.

Data

The SPARC database (Lelli et al. 2016a) comprises ≥ 200 extended HI rotation curves compiled from the literature, for which near-infrared photometry at 3.6 microns is available from Spitzer Sheth et al. (2010) and Schombert & McGaugh (2014). It forms a representative sample of disk galaxies in the local Universe in terms of various properties such as luminosity, size and surface brightness. Here, we exclude from the sample the galaxies that have been given a bad quality label (Q = 3) in the original database, leaving us with a sample of 160 galaxies.

The availability of near-infrared photometry allows one to fix in principle the stellar mass-to-light ratio as it is known that this quantity does not vary much in that band. Indeed, Schombert & McGaugh (2014) constructed stellar population synthesis models of disk galaxies to show that the stellar mass-to-light ratio does not depend strongly on age, metallicity, or color, for a large range of models with di erent star formation histories. As a first step, and following McGaugh et al. (2016), we adopt here a mass-to-light ratio d = 0.5 M § /L § for disks and b = 0.7 M § /L § for bulges. In a subsequent paper, Li et al. (2018a, hereafter L18a) adjusted, with well-chosen priors, the mass-to-light ratio d for disks and b for bulges, distances and inclinations for galaxies in the SPARC sample, showing that the radial acceleration relation (RAR) between the gravitational acceleration generated by baryons and the total gravitational acceleration could have a scatter as small as 0.057 dex, compatible with observational uncertainties on measured velocities themselves. As a second step, we adopt these parameters, such that, in both cases, only the DM halo parameters are left as free parameters.

The Einasto profile fits

The Einasto halo model has been proposed by Navarro et al. (2004) as a density profile for CDM halos in DMO CDM simulations, and had been previously introduced for the distribution of stellar light and mass in early-type galaxies (Einasto 1965(Einasto , 1968(Einasto , 1969)). Its density profile reads:

fl E (r) = fl ≠2 exp C ≠2n A 3 r r ≠2 4 1/n ≠ 1 BD , ( 1
)
where, r ≠2 is the radius at which the density profile has a slope of ≠2, and fl ≠2 is the density at that radius. The third parameter, n, is the Einasto index and sets the general shape of the density profile. As shown in, e.g., Chemin et al. (2011), increasing the Einasto index at fixed r ≠2 and fl ≠2 increases and steepens the density profile in the central part of the halo, making it more cuspy, the central density being given by fl 0 = fl ≠2 exp(2n). Hence n ≥ 6 DMO halos are indeed cuspy, even though they display a tiny core very close to the center. Conversely, at fixed n ≥ 1, the profile is more cored, and the characteristic radius r ≠2 is a good proxy to give an idea of the core-size. Chemin et al. (2011) have used for the first time the Einasto model for mass decompositions of RCs on a sample of spiral galaxies from the THINGS survey. They found that the RCs are significantly better fit with the Einasto halo than with either the isothermal or NFW halo models. Here, in the same spirit, we use the Einasto DM halo profile for the decomposition of rotation curves of the SPARC database. This will allow us to display up-to-date scaling relations for DM and baryons, complementary to the recent study of Li et al. (2018b).

For finding the structural parameters of the fitted DM halo model, we decompose the observed rotation curves into four components: the stellar disk and bulge, the gaseous disk and the DM halo. Only the (spherical) DM halo parameters are free in our fits. The observed rotation curves can be written as:

V 2 rot = d V 2 d + b V 2 b + V 2 gas + V 2 halo (2)
where V d is the rotation velocity of the stellar disk, and V b of the stellar bulge, both for a stellar mass-to-light ratio of one. Finally, V halo is the rotation velocity of the DM halo which we are fitting to the data with three free parameters. We fix all the other galaxy parameters (including distance and inclination) (i) to the original SPARC values (Lelli et al. 2016a) with fixed mass-to-light ratios d = 0.5 M § /L § and b = 0.7 M § /L § , which we call the "MLconst" case hereafter, and (ii) to the L18a galaxy parameters.

DMO-inspired fits with virial mass as only parameter

For illustrative purposes, at first we fixed the Einasto index to n = 6 and use the c ≠ M 200 relation of DMO simulations in Dutton & Macciò (2014),

log(c) = 0.997 ≠ 0.13 ◊ log( M 200 10 12 h ≠1 ). (3)
Then, the virial mass M 200 is the only free parameter, which we can use to fit the last observed point of the rotation curve. Note that there is in fact non-negligible scatter in this DMO relation, as well as small variations of n, which we do not account for here. These fits are mainly made here for illustratory purpose, showing what typical rotation curve shape is expected at a given virial mass in the absence of any feedback from the baryons. As expected, this leads to halos that are generally too cuspy. We show this curve for the L18a case in the Appendix (Fig. A.1) as the blue dashed curve. It is clear that many of these fits are not satisfactory, especially for low-mass galaxies.

Three-parameters fits

We then use the Einasto DM halo profile from Eq. 1 with all three parameters being free in the fit. We use the a ne Invariant Markov chain Monte Carlo Ensemble sampler, from the open source Python package, emcee (Foreman-Mackey et al. 2013) to fit the observational velocity curve to the theoretical model by finding the peak of the likelihood function. Here, we do not use any CDM inspired prior, and thus simply assume a flat linear prior for the three free parameters. The MCMC sampler then estimates the posterior probability distribution for all these parameters.

The best fit values of the parameters fl ≠2 , r ≠2 and n are chosen at the maximum likelihood for each galaxy. We can then compare the quality of fits in the "MLconst" and L18a cases, by comparing the global likelihoods. Interestingly, the "MLconst" case, which is most often used in the literature to compare simulations to data, gives less good fits than the L18a case, despite the L18a marginalization making absolutely no assumption on the dark matter halo parameters themselves. The global likelihood ratio is of the order of 3 -and there are 15 galaxies with a likelihood ratio higher than 10 -in favour of the L18a case, which is not highly significant, but still indicates a slight global preference for the L18a parameters.

The results of the best fit parameters in the L18a case, together with their associated uncertainties, for all 160 galaxies are listed in Table A.1 in the Appendix. The rotation curves of these galaxies are shown in

Scaling relations

We now explore the scaling relations between the free parameters of the DM halo itself, as well as DM-baryons scaling relations.

Dark matter halo scaling relations

We start by plotting the histogram of Einasto indices on Fig. 1, for both cases. The first interesting thing to note here is that the majority of Einasto indices is comprised between 0 and 2, even though the number of very small n is slighlty reduced in the L18a case. This is in line with the findings of Chemin et al. (2011) and indicates that the profiles tend to be more cored than in DMO simulations.

Second, like Chemin et al. (2011), we find a tight and strong anti-correlation between the characteristic halo radius, r ≠2 and the characteristic halo density, fl ≠2 (Fig. 2). Interestingly, the slope of this relation becomes closer to the Chemin et al. (2011) one for the L18a parameters. The scatter of the relation also becomes smaller in the L18a case, decreasing from 0.35 dex in the "MLconst" case to 0.28 dex in the L18a case.

The best fit relation for the L18a parameters is

log(fl ≠2 ) = (≠1.32 ± 0.15) ◊ log(r ≠2 ) ≠ (1.27 ± 0.18). (4)
This anti-correlation is reminiscent of the scaling relation found with Burkert halos between the central density and core radius (Donato et al. 2009), although, as shown also by Li et al. (2018b), it does not imply a constant surface density scale fl ≠2 ◊ r ≠2 independent of luminosity. Indeed, we find a linear relation between the value of log(fl ≠2 ◊r ≠2 ) and the log of the 3.6 micron luminosity, with a slope 0.13 in the "MLconst" case, and a slope of 0.1 in the L18a case. This relation is related to the absence of any significant correlation between the fl ≠2 density scale and the luminosity, already put forward by Li et al. (2018b) (see Fig. 3).

We nevertheless also note that, contrary to fl ≠2 , the central volume density of DM haloes, fl 0 , does display some correlation with luminosity, due to the higher prevalence of high Einasto indices n (and hence more cuspy halos) among galaxies of higher luminosities, as shown on the lower panels of Fig. 3.

The latter result brings us to the correlation of the Einasto index n with fl ≠2 and r ≠2 , displayed on Fig. 4. The correlation between r ≠2 and n is rather weak in the "MLconst" case, with a slope consistent with zero and a scatter of 0.4 dex, but it becomes steeper with a scatter of 0.3 dex for the L18a parameters:

log(n) = (+0.5 ± 0.15) ◊ log(r ≠2 ) ≠ (0.16 ± 0.18).
(5)

Similarly, the anticorrelation between fl ≠2 and n has a small slope, with a large scatter of 0.4 dex in the "MLconst" case, but becomes steeper and with a smaller scatter of 0.23 dex for the L18a parameters:

log(n) = (≠0.4 ± 0.11) ◊ log(fl ≠2 ) ≠ (0.71 ± 0.29). (6)

Dark matter-baryon scaling relations: diversity and uniformity

We now concentrate on scaling relations between DM and baryons. The first well-known such scaling relation, not reproduced here, is of course the baryonic Tully-Fisher relation (BTFR) between the baryonic mass of disk galaxies Fig. 2: The anti-correlation between the characteristic halo radius, r ≠2 , and the characteristic halo density, fl ≠2 for the 160 galaxies of the SPARC sample (left: "MLconst"; right: L18a). In this plot the black dots with error bars are the best fit parameters derived from fitting the rotation curves, the red doted line is the line of Chemin et al. (2011), and the green solid line is the best fit line to the data. The green shaded area represents the 1 sigma uncertainty on this best fit line (confidence interval, which does not represent the intrinsic scatter around the relation).

and the asymptotic circular velocity, which has an intrinsic scatter smaller than 0.1 dex (McGaugh et al. 2000;Lelli et al. 2016b), a priori posing a fine-tuning problem in the CDM context (Desmond 2017b). Semi-empirical models with cored DM haloes can reproduce the observed slope and normalization, but not this small scatter (Di Cintio & Lelli 2016). But what is more, galaxies with the same asymptotic circular velocity -hence twins of identical baryonic mass on the BTFR -can display a very broad range of rotation curve shapes, which we could call the "BTFR twins paradox". This unexpected diversity of shapes of rotation curves has been nicely put forward in, e.g., Oman et al. (2015), where the circular velocity at 2 kpc from the center was shown to be extremely diverse for a given maximum circular velocity of the halo. We reproduce this diversity plot in Fig. 5, for both sets of galaxy parameters. It should be noted that the "MLconst" parameters are often used in the literature to compare directly the data to simulations. However, as we showed hereabove, the L18a parameters yield better fits to the SPARC galaxy rotation curves. It is therefore interesting to compare the diversity of rotation curve shapes in the two cases: whilst the global picture remains the same, it is however interesting to note that some of the most extreme low values of V rot (2 kpc) in the "MLconst" case become more in line with the general population when switching to the L18a parameters. Hence some of the most extreme cored DM profile cases which are challenging to reproduce in simulations can be easier to reproduce with the L18a parameters. Nevertheless, one important information missing from the usually presented diversity plot is the close relation between the value of the observed velocity at 2 kpc and the baryon-induced velocity at this radius, produced in the absence of DM. This is illustrated on the bottom row of Fig. 5, where the points are color-coded by log[V b (2 kpc)]. One can clearly see how smaller baryonic velocities at 2 kpc are associated with smaller total velocities V rot (2 kpc) at this radius. Of course, while this figure clearly shows that the baryons are an important part of the story, the fact remains that the central mass deficit is enormous relative to a cuspy DM halo, sometimes higher than the total baryonic mass of the galaxy itself (Oman et al. 2015).

In order to display a few concrete examples of this relation between the baryon-induced velocities and the diversity of measured velocities at 2 kpc, we show for the L18a parameters (and on the same R-axis scale) on Fig. 6 the rotation curves of three galaxy pairs with similar asymptotic velocities, hence "BTFR twins", but very di erent measured velocities at 2 kpc from the center. In the case of IC 2574, note that we added a prior on the size of the core not being larger than 5 times the last observed point of the rotation curve. We alternatively added to the SPARC data one outer data point, coming from the THINGS survey (Oh et al. 2008) at large radius, recalibrated for the exact same distance and inclination as quoted in L18a. Those two approaches gave almost exactly the same fit. What is immediately striking on Fig. 6 is that the more compact galaxies in terms of baryons are also more compact in terms of their DM distribution. This is of course expected from adiabatic contraction of DM halos, but the problem is that this must happen after an e cient core formation mechanism, and that going back to very steep DM distributions for the highest baryonic surface densities, through adiabatic contraction, can then be complicated, as we dis-cuss below in the context of hydrodynamical simulations of galaxy formation. Here, the gravitational acceleration generated by baryons at 2 kpc is much smaller in F568-3, UGC 5750 and IC 2574 than in NGC 024, UGC 5721 and UGC 7690, and so is the measured circular velocity, and hence the total gravitational acceleration at 2 kpc. This is highly suggestive of the fact that the diversity of rotation curve shapes is driven by the diversity of baryonic distri-Fig. 6: The BTFR twins "paradox". Rotation curves of three pairs of galaxies with L18a parameters, each pair having very similar asymptotic circular velocities, but very di erent velocities at 2 kpc from the center. In these plots the black dots with error bars are the observed rotation curves of each galaxy. The green doted curves are the baryonic contribution to the rotation curves, the violet dot-dashed curves are the rotation curves of the DM halo and the red solid curves are the resulting best-fit rotation curves to the observed data. butions at a given maximum velocity scale, together with a tight baryon-DM scaling relation. This relation is the RAR. What this means is that, while a tight RAR as obeyed in the L18a case slightly decreases the diversity compared to the original "MLconst" case, the diversity of baryon-induced accelerations at 2 kpc of the SPARC galaxies is su cient to induce a large diversity of rotation curve shapes with a tight RAR.

g obs = F(g bar ) = g bar 1 ≠ e ≠ Ô g bar /g † , ( 7
)
with the only free parameter g † , McGaugh et al. 2016 andLelli et al. (2017) found a good fit to the RAR of the SPARC sample of galaxies with g † = (1.20 ± 0.02) ◊ 10 ≠10 m s ≠2 . With the L18a parameters, the RAR scatter is as small as 0.057 dex. In the outer flat part of rotation curves, this relation is equivalent to the BTFR (e.g. Lelli et al. 2016b), but it is important to remember that the BTFR does not imply Fig. 7: Top panels: correlation between the stellar disk scale length R d and the dark matter halo parameter r ≠2 (left: "MLconst"; right: L18a). Bottom panels: for cored profiles with n < 2, correlation between the stellar disk scale-length and the core size r 0 defined as the radius at which the DM density becomes half the central one (left: "MLconst"; right: L18a).

the RAR in the central parts of galaxies, where the rotation curves are rising, and where the RAR is more closely connected to the size of the DM cores.

Indeed, the RAR implies that the rotation curve shape is actually correlated with the surface density of the baryons (see Figure 15 of Famaey & McGaugh 2012), as for instance illustrated in Lelli et al. (2013Lelli et al. ( , 2016c)). Also, the rotation curve shapes of late-type spiral galaxies in the lowacceleration regime should all be similar when expressed in units of disk scale-length (Swaters et al. 2009;McGaugh 2014), at least for perfect exponential disks. Yet in other words, one also expects that the size of the DM core will be directly correlated with the disk scale-length, a finding already put forward by Donato et al. (2004). However, gaseous distributions are not exponential, so that any plotted correlation using the stellar disk scale-length as a proxy for the extent and "flu ness" of a disk with a certain baryonic mass will be much less tight than the RAR itself. Fig. 10: Left: The RAR at 1 kpc with L18a parameters. Observational errors from the rotation curve measurements are also plotted. Note that those large errors seem to be associated with a downward scatter in g obs w.r.t. the RAR. Right: The RAR at 5 kpc. The observational errors are small and the RAR is particularly tight at this radius. The 6 representative galaxies of Fig. 6 are displayed as in Fig. 9. Fig. 11: Diversity plot from simulations: the circular velocity at 2 kpc from the center V rot (2kpc) versus the circular velocity at the last observed point V Rlast from the NI-HAO (cyan triangles) and APOSTLE/EAGLE (red points) projects. The grey shaded area indicates the region spanned by observed galaxies in SPARC (i.e., the black points of Fig. 9). The six representative galaxies are plotted as in Fig. 9.

If we restrict ourselves to cored profiles with n < 2, and define the core size r 0 as the radius at which the density becomes half the central one fl 0 , we end up with a similar correlation, which we display on the lower panels of Fig. 7. In that case, the scatter is reduced from 0.5 dex in the "MLconst" case to 0.4 dex in the L18a case. The L18a correlation of the lower-right panel of Fig. 7 for n < 2 reads log(r 0 ) = (0.59

± 0.3) ◊ log(R d ) + (0.04 ± 0.14). (9)
Finally, a way to illustrate the influence of baryons on the diversity of central DM profiles in rotationallysupported galaxies is to look at the averaged DM density within 2 kpc, and compare it to the baryon-induced velocity at that radius. We indeed found on Fig. 8, in both the "MLconst" and L18a cases, a very well-defined relation for galaxies with low-enough baryonic velocities at 2 kpc, V b (2 kpc) < 80 km s ≠1 , meaning that the baryons are not heavily dominating the measured circular velocity at 2 kpc (which would imply a very bad constraint on the enclosed DM mass). The log-log slope of the relation is 0.93 in the "MLconst" case, with a scatter of 0.3 dex, and 0.87 in the L18a case with a scatter of only 0.08 dex. In the L18a case, the relation for

V b (2 kpc) < 80 km s ≠1 reads log(f(2 kpc)) = (0.87 ± 0.4) ◊ log(V b (2 kpc)) ≠ (2.84 ± 0.58). ( 10
)
Let us note that such a relation can actually directly be inferred from the RAR, a relation which we overplot as a dotted red curve on Fig. 8. Let us note that this Figure is in fact closely related to the color gradient visible on Fig. 5, as the average DM density within 2 kpc also reads as f = (V 2 rot ≠ V 2 b )/GR 2 , where R = 2 kpc.

Eq. 9 and Eq. 10 can be considered to represent the strong version of the old core-cusp problem. Not only do late-type low-mass galaxies often require DM cores, the averaged density and size of these cores are closely correlated to the baryonic distribution. Together with the diversity of baryonic surface density distributions at a given DM halo mass scale, this is at the root of the diversity of rotation curve shapes. This diversity of rotation curve shapes is actually an important constraint to take into account in hydrodynamical simulations of galaxy formation, closely related to the ability of simulations at transforming initially cuspy DM profiles into cores (e.g., Navarro et al. 1996a;Pontzen & Governato 2012). As shown in Oman et al. (2015), cosmological hydrodynamical simulations from the EAGLE and APOSTLE projects are unable to produce large cores as required by the low values of V rot (2 kpc) in a large number of galaxies. This issue was explored further by Benitez-Llambay et al. ( 2018) who showed that the core formation crucially depends on the star formation gas density treshold. They showed that high tresholds are needed for baryons to collapse at the center at densities comparable to DM, before being ejected by feedback, in order for cores to form under the influence of a small number of distinctive gas blowouts. A simulation with high star formation gas density threshold resulted in galaxies with sizeable cores over a limited halo mass range, but still insu cient variety in mass profiles to explain the observed diversity of galaxy rotation curves. Other simulations do create cores more easily than those based on EAGLE. For instance, Di Cintio et al. ( 2014) used a suite of simulated galaxies drawn from the MaGICC project to investigate the e ects of baryonic feedback on the density profiles of DM haloes, producing cored profiles. Read et al. (2016Read et al. ( , 2018) also used high resolution simulations of isolated dwarf galaxies to study the physics of DM cusp-core transformation, which was argued to be linked to short potential fluctuations associated with bursty star formation over a long period (in contradiction to Benitez-Llambay et al. 2018). Most of these simulations also showed how the general shape of the RAR was naturally obtained when including baryonic feedback (e.g., Ludlow et al. 2017;Keller & Wadsley 2017). Perhaps the most impressive simulations in that sense are the recent NIHAO simulations (e.g., Santos-Santos et al. 2018), which are argued (Dutton et al. in prep. 1 ) to produce a tight RAR with a scatter of order 0.08 dex over the entire distance range. Nevertheless, we argue hereafter that the real challenge is to produce both a tight RAR and the observed diversity of rotation curve shapes. To illustrate the relation between a tight RAR and the diversity problem, we plot, on Fig. 9, 124 galaxies with L18a parameters (rejecting galaxies with uncertainty> 10% in V rot (2 kpc) and those with g obs < g bar at 2 kpc), compared to the RAR, at 2 kpc only. This means that, contrary to the traditional RAR plot, each point on this plot represents a di erent galaxy. Eq. 7 is overplotted together with those points. The scatter around the relation is ≥ 0.06 dex. The six representative galaxies plotted on Fig. 6 are displayed as squares and stars on this plot. These selected galaxies can then illustrate the e ect of the RAR on the diversity plot, which we reproduce for the same set of 124 galaxies on the right panel of Fig. 9. Note that, in this case, the scatter of V rot around the median at a given V Rlast is 1 see darkmatter2018.weebly.com/program.html ≥ 0.2 dex. This can then be compared directly to the points from Figure 2 of Santos-Santos et al. (2018). Interestingly, it appears that the NIHAO simulations are very e cient at producing cores reproducing the right V rot (2 kpc) for galaxies sitting at the bottom of the RAR, with a low baryonic acceleration at 2 kpc. But the three high-baryonic acceleration galaxies of Fig. 6, NGC 024, UGC 5721 and UGC 7690 seem to all have too high observed V rot (2 kpc) to be reproduced by the NIHAO simulations. On the left panel of Fig. 9, we also overplot the NIHAO points on the RAR at 2 kpc from Dutton et al. (in prep.), spanning the same baryonic acceleration range as our sample. Interestingly, the NIHAO simulations seem to slightly underpredict the observed RAR at 2 kpc in the regime of baryonic accelerations ranging from ≥ 10 ≠11 to 10 ≠10 ms ≠2 , corresponding to the regime of our three representative high baryonic acceleration galaxies. One interesting thing to note is that ? showed that 'observing' and modelling the same simulated galaxies as in Oman et al. (2015) to measure their rotation curves, rather than simply measuring the rotation curve from the mass profile, introduced an additional downward scatter in the diversity plot. Interestingly, when considering the observed RAR for the same galaxies as in Fig. 9 at di erent radii, one clearly sees that the scatter increases downwards when the quoted observational errors increase (Fig. 10). At 1 kpc, the observational errors become quite important (whilst non-circular motions might also become more important than in the outskirts), and many galaxies scatter down in the RAR. At 5 kpc, on the other hand, the quoted observational errors are smaller, and the tightness of the RAR becomes impressive. In terms of diversity, that could mean that results such as those from NIHAO are actually worse than they appear, since if the galaxies were realistically observed and modelled they could scatter down even more in terms of V rot (2 kpc). Additionally, it would be important to check whether the relation between V max and V Rlast in NIHAO are in line with observational fits. Interestingly, as originally shown by Oman et al. (2015), galaxies such as NGC 024, UGC 5721 and UGC 7690 can in principle be more closely reproduced in the diversity plot by simulations from the EAGLE and APOSTLE projects, which in turn cannot reproduce the smaller V rot (2 kpc) galaxies. This is illustrated in Fig. 11 where the diversity plot of the NIHAO and APOSTLE/EAGLE simulations are compared. Even taking into account observational errors, it is unlikely that the APOSTLE/EAGLE points will scatter down enough. On the other hand, NIHAO points will probably scatter down too much once such errors are taken into account. This might indicate that the answer lies somewhere between the two types of simulations, combining cores slightly more modest than in NIHAO with some contribution from rotation curve errors (to add some scatter and downward shift). But in any case, reproducing both a tight RAR and the observed diversity of rotation curve shapes thus raises an interesting challenge for simulations of galaxy formation.

Discussion and conclusion

In this paper, we used the Einasto DM halo profile for the decomposition of rotation curves, in the same spirit as Chemin et al. (2011) and Li et al. (2018b). Our study is complementary to the latter, as we compared here the quality of the fits and the scaling relations obtained for dif-ferent sets of galaxy parameters, and focused in particular on scaling relations between the characteristics of the DM core of low-mass galaxies and their baryonic distributions.

The two sets of galaxy parameters used in this study were the original SPARC parameters with constant massto-light ratios for bulges and disks (Lelli et al. 2016a;Mc-Gaugh et al. 2016), and the parameters for which galaxies follow the tightest radial acceleration relation (RAR, L18a). We found that the fits were globally slightly better in the latter case, and that the dark halo parameter scaling relations also became a bit tighter. In particular, we found for both sets of parameters a tight and strong anti-correlation between the characteristic halo radius, r ≠2 and the characteristic halo density, fl ≠2 . The latter is also anti-correlated with the Einasto index n, but in a tighter and steeper way for the L18a parameters than for the "MLconst" ones.

In terms of dark matter-baryons scaling relations, we focused on relations between the core properties and the extent of the baryonic component, which are highly relevant to the cusp-core transformation process. For galaxies with Einasto index n < 2, we found a positive correlation between the core size, defined as the radius where the density reaches half of its central value, and the stellar disk scale-length (Eq. 9), albeit with some significant scatter. A tighter relation has been found between the averaged dark matter density within 2 kpc and the baryon-induced rotational velocity at that radius (Eq. 10), for galaxies not too strongly dominated by baryons in the center, namely where this velocity does not exceed 80 km s ≠1 . These scaling relations could be referred to as the strong version of the old core-cusp problem.

These scaling relations are directly related to the consequence of the RAR on the diversity of rotation curve shapes (Oman et al. 2015), quantified by the observed rotational velocity at 2 kpc, at a given maximum or asymptotic velocity scale. This diversity of rotation curve shapes could be referred to as the BTFR twins "paradox", since the a priori expectation, without e cient feedback, would have been a tight relation between the rotational velocity at 2 kpc and the maximum halo velocity for dark matter dominated galaxies. We confirmed this observed diversity of rotation curve shapes, but pointed out how it was directly related to the tightness of the radial acceleration relation together with the diversity of baryon distributions in galactic disks. Indeed, while we showed that a tight RAR with the L18a parameters slightly decreases the observed diversity compared to the "MLconst" case, we also showed that, in the L18a case, the diversity of baryon-induced accelerations at 2 kpc toghether with the RAR are su cient to induce a large diversity of rotation curve shapes.

Reproducing all the scaling relations displayed in this paper together, and in particular the RAR and the diversity of rotation curve shapes, side by side, is still challenging for current hydrodynamical simulations of galaxy formation in a cosmological context, as we illustrated in Figs. 9-11. This challenge could perhaps imply something fundamental on the nature of dark matter (Kamada et al. 2017;Famaey et al. 2018) or even gravity (Milgrom 1983;Famaey & McGaugh 2012;Smolin 2017;Verlinde 2017). In the standard context, it might seem to call for an apparently fine-tuned feedback from the baryons, where galaxies would self-regulate their star formation, implying an interplay between gas inflow, star formation, feedback, and the final DM distribution, which would have to conspire to make galaxies find an attractor acceleration relation with little scatter over the entire acceleration range, while maintaining a high diversity of rotation curve shapes.

Fig.

A.1: Rotation curve fits for SPARC galaxies with L18a parameters. In each panel the black dots represent the measured rotation velocities and their observational uncertainties. The red curves are the best-fit curves to the data with our three-parameters fits. The green dashed curves represent the contribution of the baryonic matter to the rotation curves and the dark matter halo rotation velocities are plotted by violet dot-dashed curves. The blue dashed curves are there for visual comparison only, using n = 6 and the virial mass as the only free parameter to give an idea of the typical shape expected in the DMO case.

Table A.1: Maximum posterior DM halo parameters of individual rotation curve fits to 160 SPARC galaxies with Q < 3 and fixed mass-to-light ratios, inclinations and distances from L18a. The first column is the name of the galaxies. The columns 2, 3 and 4 are the the parameters of the Einasto halo profile, together with the 68% confidence interval. Column 5 is the maximum circular velocity of the DM halo V max (halo) and column 6 is the measured rotational velocity at 2 kpc from the center V rot (2kpc).

Galaxy fl ≠2 r ≠2 n V max (halo) V rot (2kpc) 10 ≠3 M § pc ≠3 kpc

  Fig. A.1.

Fig. 1 :

 1 Fig.1: Histogram of Einasto indices n for the fits to the 160 galaxies with Q < 3 in the SPARC sample. Left panel: original SPARC parameters with constant mass-to-light ratios for bulges and disks ("MLconst",Lelli et al. 2016a;McGaugh et al. 2016). Right panel: parameters fromLi et al. (2018a, L18a).

Fig. 3 :

 3 Fig. 3: Top panels: the characteristic density fl ≠2 at the radius of density slope ≠2 versus the 3.6 micron luminosity (left: "MLconst"; right: L18a), compared to the density scale of Li et al. (2018b). Bottom panels: the central density fl 0 versus the 3.6 micron luminosity (left: "MLconst"; right: L18a), colour-coded by the value of the Einasto index n.

Fig. 4 :

 4 Fig.4: Correlation between Einasto Index n and dark matter halo parameters r ≠2 (top) and fl ≠2 (bottom), in the "MLconst" (left) and L18a (right) cases.

Fig. 5 :

 5 Fig. 5: Top panels: the measured circular velocity at 2 kpc from the center V rot (2kpc) versus the maximum circular velocity of the DM halo V max (halo) (left) and the last observed point of the rotation curve V Rlast (right). The "MLconst" values are denoted with red diamonds and the L18a values with blue circles. The grey dotted line is the 1:1 line. Bottom panels: The same plots for the L18a points color-coded by the log of the baryon-induced velocity at 2 kpc, log[V b (2 kpc)].

  By fitting the following function proposed byMc- Gaugh (2008),Famaey &McGaugh (2012) and McGaugh & Schombert (2014)

Fig. 7

 7 Fig.8: The "strong" version of the old core-cusp problem. Correlation between the averaged DM density within 2 kpc f(2 kpc) and the baryon-induced velocity at 2 kpc, V b (2 kpc) (left: "MLconst"; right: L18a). The red-dotted curve is the correlation induced by the RAR in Eq. 7. Note how this plot is closely related to the color gradient visible on Fig.5, as the average DM density within 2 kpc simply reads as f = (V 2 rot ≠ V 2 b )/GR 2 , where R = 2 kpc.

Fig. 9 :

 9 Fig.9: Left: The radial acceleration relation (RAR) at 2 kpc (see text) with L18a parameters. The 6 representative galaxies plotted on Fig.6are displayed as green, magenta and orange squares and stars on this plot. The grey dotted line is the 1:1 line. The red dashed curve is the RAR from Eq. 7. UGC 7690 and NGC 024 fall particularly well on the relation. Also overplotted as cyan triangles are the points of the RAR at 2 kpc from the NIHAO simulations ofDutton et al. (in prep.). Right: Diversity plot: the circular velocity at 2 kpc from the center V rot (2kpc) versus the circular velocity at the last observed point V Rlast for the same galaxies. The cyan triangles are the predictions of the NIHAO simulations from figure 2 ofSantos-Santos et al. (2018). UGC 7690 and NGC 024 appear to have too high velocities V rot (2kpc) and hence to high central DM concentrations (n = 5.5 and 4.7, respectively) compared to simulations with very e cient core formation.

  Article number, page 28 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves Table A.1: Continued.

			Table A.1: Continued. Table A.1: Continued.
	IC 2574 UGC 04483 D 564-8 F 565-V2 UGC 05750 NGC 3741 UGC 07577 KK 98-251 UGC 00891 NGC 3109 UGC 05005 UGC 07559 UGC 06818 F 571-V1 UGC 08837 UGC 05999 DDO 154 D 512-2 UGC 05918 UGC 07866 UGCA 444 UGC 01281 F 568-3 UGCA 442 D 631-7 DDO 161 UGC 06667 UGC 04278 UGC 05829 F 583-1 UGC 07608 F 563-1 DDO 170 NGC 0055 ESO 444-G084 UGC 07089 NGC 0100 UGC 09992 F 583-4 UGC 02023 Galaxy DDO 064 F 571-8 UGC 05986 NGC 5585 UGC 07524 UGC 06399 UGC 08550 NGC 4068 F 568-1 UGC 00731 UGC 05414 NGC 0300 DDO 168 UGC 07125 UGC 07323 UGC 07603 UGC 00634 UGC 05764 ESO 116-G012 ESO 079-G014 UGC 05716 NGC 2915 UGC 04499 F 574-1 UGC 11820 UGC 12632 UGC 09037 ESO 563-G021 F 568-V1 UGC 12732 UGC 01230 NGC 4010 UGC 10310 UGC 06923 UGC 11557 NGC 1003 NGC 3917 F 563-V2 UGC 07690 UGC 06917 Galaxy NGC 0247 UGC 06446 NGC 4214 UGC 07399 UGC 07151 UGC 07261 UGC 00128 UGC 06930 UGC 06628 UGC 06983 UGC 08286 UGC 02259 UGC 08490 UGC 00191 UGC 03580 NGC 3972 UGC 05721 NGC 4183 UGC 04325 UGC 11455 NGC 2976 NGC 4559 NGC 7793 NGC 3198 F 579-V1 NGC 3726 NGC 3769 NGC 4085 NGC 2403 NGC 0024 NGC 3877 NGC 6503 NGC 1090 NGC 4217 NGC 4088 NGC 4100 NGC 4157 NGC 4051 NGC 3949 NGC 6015 Galaxy NGC 6946 NGC 2903 NGC 3893 IC 4202 UGC 12506 NGC 3953 NGC 4138 UGC 06786 NGC 2683 NGC 2998 UGC 03205 NGC 7331 NGC 0289 UGC 08699 NGC 5985 NGC 3521 NGC 5055 NGC 0891 NGC 5907 NGC 4013 NGC 6674 NGC 5033 UGC 06614 NGC 2841 UGC 03546 NGC 3992 NGC 5371 NGC 2955 UGC 02916 NGC 5005 NGC 0801 UGC 11914 UGC 02953 NGC 6195 UGC 05253 NGC 7814 UGC 06787 UGC 02885 UGC 09133 UGC 02487	0.222 +0.057 ≠0.030 28.642 +21.184 ≠13.217 2.459 +2.059 ≠1.229 1.798 +1.764 ≠0.919 2.793 +1.567 ≠1.205 0.137 +0.060 ≠0.038 1.198 +1.089 ≠0.744 29.707 +8.561 ≠8.084 2.616 +1.012 ≠0.862 2.224 +1.033 ≠0.811 0.856 +0.672 ≠0.392 6.520 +5.000 ≠2.938 3.604 +2.664 ≠1.652 2.037 +1.399 ≠0.882 3.560 +0.967 ≠1.031 2.897 +1.134 ≠1.063 4.948 +0.405 ≠0.447 5.405 +5.798 ≠3.132 4.378 +4.800 ≠2.319 9.026 +3.668 ≠2.920 0.318 +0.247 ≠0.135 8.332 +2.342 ≠2.210 5.381 +1.622 ≠1.437 3.418 +1.452 ≠1.219 5.383 +1.153 ≠0.913 2.039 +0.618 ≠0.561 3.772 +1.454 ≠1.353 0.969 +0.621 ≠0.334 1.065 +0.725 ≠0.416 4.474 +0.904 ≠0.923 4.481 +3.465 ≠1.940 2.774 +1.481 ≠1.161 3.438 +1.080 ≠1.061 4.178 +0.672 ≠0.734 1.776 +1.318 ≠0.692 1.039 +0.783 ≠0.442 2.924 +1.206 ≠1.154 52.805 +29.457 ≠21.937 0.379 +0.405 ≠0.192 9.438 +8.152 ≠6.172 fl ≠2 10 ≠3 M § pc ≠3 14.964 +8.384 ≠5.633 2.004 +0.295 ≠0.323 5.732 +0.583 ≠0.612 0.055 +0.023 ≠0.014 4.619 +1.287 ≠1.324 4.367 +2.589 ≠1.907 5.043 +3.138 ≠2.138 8.744 +4.662 ≠4.121 6.469 +3.656 ≠2.767 3.671 +2.217 ≠1.604 5.383 +3.403 ≠1.979 2.046 +1.057 ≠0.786 21.386 +2.109 ≠2.035 5.525 +2.322 ≠2.076 0.873 +0.582 ≠0.344 8.408 +2.621 ≠2.604 2.062 +0.552 ≠0.652 24.718 +2.077 ≠2.454 4.156 +1.028 ≠1.105 4.281 +0.593 ≠0.646 2.918 +0.626 ≠0.621 6.776 +1.775 ≠1.746 6.743 +2.981 ≠2.732 5.886 +1.508 ≠1.590 0.421 +0.251 ≠0.131 5.153 +2.719 ≠2.193 2.915 +0.708 ≠0.738 4.822 +0.398 ≠0.395 4.169 +2.880 ≠2.058 0.494 +0.317 ≠0.172 5.780 +2.134 ≠1.870 9.404 +1.780 ≠1.956 7.338 +5.073 ≠3.890 9.323 +6.262 ≠4.507 8.403 +4.851 ≠3.906 0.139 +0.040 ≠0.020 7.790 +1.181 ≠1.248 11.762 +7.509 ≠5.497 39.453 +34.945 ≠22.461 5.000 +2.276 ≠2.033 fl ≠2 10 ≠3 M § pc ≠3 0.733 +0.383 ≠0.230 3.942 +2.794 ≠1.686 0.029 +0.010 ≠0.009 3.999 +1.776 ≠1.444 16.330 +3.992 ≠4.088 4.012 +3.755 ≠1.870 0.600 +0.147 ≠0.129 3.510 +2.565 ≠1.777 23.971 +20.973 ≠13.054 4.681 +2.002 ≠1.928 5.913 +1.811 ≠1.708 0.272 +0.242 ≠0.134 6.569 +2.020 ≠1.931 4.650 +1.933 ≠1.504 0.113 +0.047 ≠0.041 2.931 +2.789 ≠1.529 14.790 +2.841 ≠2.883 3.938 +2.536 ≠1.821 25.617 +5.925 ≠6.185 3.126 +0.402 ≠0.413 14.343 +9.634 ≠6.094 3.221 +1.486 ≠1.272 10.009 +5.189 ≠4.278 2.904 +0.264 ≠0.273 12.317 +9.315 ≠6.449 1.169 +1.140 ≠0.637 2.629 +1.182 ≠1.044 16.594 +7.896 ≠6.395 2.150 +0.155 ≠0.158 2.531 +2.309 ≠1.311 30.472 +5.744 ≠5.087 2.789 +0.359 ≠0.365 4.481 +0.849 ≠0.852 6.426 +1.956 ≠1.990 1.045 +1.101 ≠0.595 8.527 +1.980 ≠2.060 0.813 +0.879 ≠0.448 14.806 +11.304 ≠8.174 0.186 +0.191 ≠0.105 0.104 +0.013 ≠0.010 fl ≠2 10 ≠3 M § pc ≠3 6.090 +1.626 ≠1.649 3.443 +0.222 ≠0.224 6.532 +2.323 ≠2.412 103.974 +6.572 ≠6.181 3.022 +1.865 ≠0.962 12.727 +9.177 ≠5.933 20.322 +9.162 ≠6.553 0.795 +0.169 ≠0.154 8.579 +4.530 ≠3.319 0.054 +0.024 ≠0.014 0.162 +0.063 ≠0.039 0.206 +0.066 ≠0.045 1.534 +0.570 ≠0.520 1.019 +0.590 ≠0.419 1.839 +0.572 ≠0.514 0.604 +0.313 ≠0.175 3.016 +0.095 ≠0.097 7.908 +0.487 ≠0.507 0.032 +0.014 ≠0.008 0.092 +0.045 ≠0.024 0.038 +0.014 ≠0.007 4.868 +0.207 ≠0.206 0.771 +0.464 ≠0.249 0.264 +0.104 ≠0.081 4.989 +0.644 ≠0.823 0.209 +0.177 ≠0.093 0.033 +0.019 ≠0.011 6.778 +1.326 ≠1.404 9.038 +1.304 ≠1.138 6.439 +5.837 ≠3.332 0.056 +0.015 ≠0.012 0.513 +0.160 ≠0.084 1.549 +0.088 ≠0.089 0.104 +0.068 ≠0.035 3.355 +0.154 ≠0.160 1.497 +1.029 ≠0.641 0.036 +0.012 ≠0.007 0.357 +0.220 ≠0.119 1.369 +0.084 ≠0.088 1.220 +0.395 ≠0.351	41.010 +5.842 ≠7.758 0.681 +0.292 ≠0.165 3.771 +2.329 ≠1.252 9.336 +5.686 ≠3.021 8.333 +2.726 ≠1.567 26.772 +6.057 ≠5.510 6.286 +3.755 ≠3.295 1.104 +0.186 ≠0.086 6.406 +1.834 ≠1.117 7.553 +3.093 ≠1.778 18.468 +7.796 ≠4.781 2.247 +1.060 ≠0.682 10.396 +6.260 ≠4.449 9.144 +3.389 ≠2.080 8.745 +2.216 ≠2.426 9.709 +2.519 ≠1.223 3.350 +0.174 ≠0.135 2.486 +1.806 ≠0.839 3.016 +1.653 ≠0.987 1.773 +0.400 ≠0.353 17.767 +7.771 ≠6.037 3.160 +0.523 ≠0.272 7.479 +1.215 ≠0.786 4.616 +1.253 ≠0.719 5.021 +0.603 ≠0.503 7.235 +1.464 ≠0.972 6.233 +1.685 ≠0.896 25.370 +9.683 ≠9.014 12.583 +4.696 ≠3.795 6.151 +0.750 ≠0.529 4.565 +2.064 ≠1.333 8.638 +2.634 ≠1.483 4.815 +1.024 ≠0.610 6.462 +0.730 ≠0.491 6.567 +2.447 ≠1.925 14.021 +6.637 ≠4.510 7.715 +2.715 ≠1.325 0.798 +0.240 ≠0.171 20.575 +11.494 ≠7.672 24.411 +17.290 ≠15.996 r ≠2 kpc 1.935 +0.649 ≠0.410 16.404 +1.522 ≠1.033 7.076 +0.328 ≠0.257 77.844 +14.845 ≠15.275 5.887 +1.161 ≠0.671 5.958 +2.387 ≠1.257 3.522 +1.260 ≠0.800 6.027 +2.685 ≠2.587 5.920 +1.933 ≠0.931 5.168 +1.928 ≠1.145 4.404 +1.629 ≠1.257 8.694 +2.755 ≠1.749 2.488 +0.088 ≠0.080 3.859 +1.104 ≠0.624 21.897 +8.929 ≠7.112 3.170 +0.667 ≠0.372 10.893 +2.293 ≠1.127 1.520 +0.070 ≠0.055 7.552 +1.317 ≠0.750 11.476 +0.956 ≠0.658 6.552 +0.882 ≠0.617 4.214 +0.584 ≠0.410 3.985 +1.294 ≠0.655 5.764 +0.942 ≠0.561 22.838 +5.658 ≠5.504 4.389 +1.508 ≠0.824 12.604 +2.098 ≠1.292 15.615 +0.558 ≠0.481 6.219 +2.634 ≠1.451 20.767 +5.953 ≠5.270 8.057 +1.358 ≠1.037 6.450 +0.455 ≠0.366 3.819 +1.924 ≠0.837 3.838 +2.057 ≠0.896 6.400 +2.823 ≠1.172 44.511 +3.921 ≠5.825 6.530 +0.569 ≠0.404 4.080 +1.245 ≠0.626 1.176 +0.585 ≠0.311 6.466 +2.166 ≠1.110 r ≠2 kpc 17.451 +4.670 ≠4.015 5.499 +1.942 ≠1.336 67.969 +14.253 ≠12.906 6.636 +1.894 ≠1.181 2.377 +0.377 ≠0.239 4.825 +1.954 ≠1.446 29.520 +4.067 ≠3.228 8.157 +3.655 ≠1.922 2.072 +0.971 ≠0.600 6.698 +2.094 ≠1.064 4.555 +0.902 ≠0.584 23.210 +10.335 ≠6.857 4.092 +0.762 ≠0.489 5.308 +1.305 ≠0.915 56.846 +16.508 ≠10.210 10.872 +6.563 ≠3.673 2.934 +0.267 ≠0.205 7.063 +2.688 ≠1.544 2.293 +0.288 ≠0.190 19.362 +1.331 ≠1.024 5.013 +3.692 ≠1.988 9.041 +2.790 ≠1.590 4.166 +1.487 ≠0.791 11.875 +0.559 ≠0.464 3.928 +1.733 ≠0.914 19.930 +11.005 ≠6.252 10.083 +2.865 ≠1.741 4.552 +1.833 ≠0.853 11.489 +0.494 ≠0.430 8.745 +4.343 ≠2.546 3.772 +0.272 ≠0.240 9.236 +0.658 ≠0.527 9.712 +1.064 ≠0.781 9.091 +1.951 ≠1.081 21.268 +13.420 ≠7.093 7.043 +0.927 ≠0.649 26.817 +15.265 ≠8.967 4.474 +2.048 ≠1.042 68.550 +32.620 ≠26.605 73.416 +4.529 ≠4.950 r ≠2 kpc 8.014 +1.559 ≠0.966 12.808 +0.432 ≠0.392 8.903 +2.245 ≠1.187 2.963 +0.068 ≠0.067 15.019 +3.225 ≠3.228 5.656 +2.090 ≠1.350 4.875 +0.969 ≠0.960 32.248 +3.875 ≠3.069 7.056 +1.769 ≠1.493 117.195 +21.541 ≠20.608 72.112 +11.996 ≠12.249 73.686 +11.111 ≠10.858 19.070 +4.292 ≠2.655 23.136 +8.173 ≠5.227 23.983 +4.589 ≠3.167 37.370 +8.398 ≠8.344 16.610 +0.223 ≠0.206 10.205 +0.274 ≠0.263 160.201 +26.972 ≠28.459 79.864 +13.980 ≠15.559 170.883 +20.335 ≠25.775 12.658 +0.271 ≠0.275 35.184 +8.277 ≠7.621 70.651 +15.846 ≠11.760 10.845 +0.926 ≠0.556 62.548 +22.943 ≠17.175 157.174 +30.321 ≠37.190 13.523 +0.988 ≠0.719 8.555 +0.612 ≠0.553 13.183 +7.463 ≠5.353 130.564 +19.318 ≠16.852 60.221 +6.956 ≠10.498 26.697 +0.845 ≠0.760 125.123 +36.431 ≠32.578 16.351 +0.393 ≠0.357 23.241 +9.061 ≠6.048 183.521 +24.166 ≠27.329 64.595 +16.304 ≠14.781 27.905 +0.919 ≠0.795 38.987 +7.301 ≠4.915	2.248 +0.120 ≠0.161 1.637 +0.994 ≠0.736 1.523 +0.776 ≠0.555 1.323 +1.012 ≠0.685 0.891 +0.647 ≠0.435 5.106 +0.593 ≠0.532 1.765 +1.734 ≠0.801 0.297 +0.226 ≠0.168 1.116 +0.431 ≠0.347 1.113 +0.320 ≠0.269 1.866 +0.933 ≠0.712 1.160 +0.899 ≠0.584 0.751 +0.498 ≠0.363 1.802 +1.389 ≠0.881 0.757 +0.263 ≠0.174 0.831 +0.563 ≠0.354 1.019 +0.128 ≠0.109 2.417 +1.828 ≠1.252 2.535 +1.582 ≠1.115 2.220 +1.744 ≠1.026 4.680 +0.971 ≠0.758 0.448 +0.248 ≠0.187 0.642 +0.234 ≠0.183 1.403 +0.553 ≠0.432 0.576 +0.168 ≠0.162 1.872 +0.454 ≠0.380 1.311 +0.574 ≠0.412 1.634 +0.404 ≠0.385 2.576 +0.896 ≠0.669 0.872 +0.249 ≠0.191 1.125 +0.766 ≠0.529 1.777 +0.787 ≠0.529 2.017 +0.727 ≠0.515 1.005 +0.220 ≠0.170 3.111 +0.677 ≠0.676 2.534 +0.743 ≠0.590 1.286 +0.527 15.132 +9.907 4.710 +1.428 0.429 +0.771 ≠0.319 ≠1.135 ≠9.657 ≠0.370 n 0.715 +0.557 ≠0.384 0.721 +0.185 ≠0.143 0.671 +0.134 ≠0.108 7.395 +0.400 ≠0.449 1.376 +0.432 ≠0.316 1.348 +0.760 ≠0.570 2.951 +0.978 ≠0.802 0.750 +0.540 ≠0.308 0.996 +0.641 ≠0.445 2.942 +0.962 ≠0.761 1.478 +0.533 ≠0.439 2.339 +0.676 ≠0.573 0.149 +0.069 ≠0.058 2.144 +0.879 ≠0.624 2.734 +0.676 ≠0.535 1.147 +0.474 ≠0.332 1.335 +0.716 ≠0.444 0.903 +0.218 ≠0.161 1.281 +0.402 ≠0.287 0.775 +0.188 ≠0.145 2.043 +0.431 ≠0.348 1.619 +0.656 ≠0.447 1.517 +0.750 ≠0.509 1.230 +0.430 ≠0.305 4.156 +0.648 ≠0.739 2.316 +1.077 ≠0.762 1.339 +0.388 ≠0.280 0.990 +0.090 ≠0.081 2.677 +2.016 ≠1.241 4.671 +0.864 ≠0.816 1.819 +0.930 ≠0.589 0.313 +0.264 ≠0.166 2.144 +1.532 ≠0.957 1.147 +0.891 ≠0.618 0.776 +0.632 ≠0.412 7.182 +0.544 ≠0.610 0.966 +0.212 ≠0.167 0.815 +0.665 ≠0.408 5.507 +5.786 ≠3.357 1.485 +0.716 ≠0.500 n 4.096 +0.475 ≠0.525 3.380 +1.335 ≠1.106 18.336 +7.093 ≠4.497 2.081 +0.679 ≠0.536 1.264 +0.398 ≠0.288 6.241 +5.349 ≠2.666 4.528 +0.625 ≠0.572 3.230 +2.077 ≠1.329 4.718 +3.506 ≠3.020 2.193 +1.158 ≠0.738 2.612 +0.637 ≠0.489 13.833 +3.582 ≠2.767 3.350 +0.944 ≠0.720 2.979 +0.550 ≠0.490 7.941 +0.820 ≠0.648 2.124 +0.855 ≠0.724 1.380 +0.450 ≠0.319 4.321 +1.741 ≠1.300 1.352 +0.513 ≠0.338 1.430 +0.195 ≠0.163 1.050 +0.369 ≠0.312 2.659 +1.020 ≠0.781 2.472 +0.863 ≠0.668 2.591 +0.231 ≠0.199 4.137 +2.348 ≠1.491 4.304 +1.693 ≠1.321 3.471 +1.958 ≠1.156 0.569 +0.418 ≠0.283 3.466 +0.146 ≠0.134 4.738 +1.294 ≠1.142 0.640 +0.239 ≠0.182 3.565 +0.498 ≠0.419 2.420 +0.467 ≠0.369 1.224 +0.487 ≠0.354 5.359 +2.181 ≠1.606 2.225 +1.047 ≠0.668 5.163 +1.885 ≠1.613 3.376 +4.397 ≠2.321 7.252 +4.963 ≠2.458 9.786 +0.359 ≠0.336 n 2.043 +0.485 ≠0.366 2.606 +0.246 ≠0.217 1.412 +0.850 ≠0.512 0.219 +0.044 ≠0.041 10.765 +4.173 ≠3.240 5.522 +5.609 ≠3.088 1.002 +0.662 ≠0.616 3.989 +0.509 ≠0.450 5.722 +6.955 ≠3.014 20.180 +4.199 ≠3.286 9.323 +0.854 ≠0.846 6.427 +0.643 ≠0.597 4.318 +1.746 ≠1.255 4.782 +0.969 ≠0.854 3.354 +0.713 ≠0.566 4.474 +0.882 ≠0.745 1.971 +0.076 ≠0.071 0.417 +0.085 ≠0.070 19.084 +2.311 ≠2.048 13.121 +2.029 ≠1.734 18.839 +2.861 ≠2.376 0.916 +0.120 ≠0.103 2.379 +0.858 ≠0.751 6.311 +0.874 ≠0.785 1.302 +0.463 ≠0.314 23.132 +9.756 ≠7.386 12.911 +3.843 ≠2.734 0.633 +0.304 ≠0.202 0.844 +0.173 ≠0.160 1.032 +1.058 ≠0.684 11.468 +1.474 ≠1.327 3.211 +0.321 ≠0.324 2.829 +0.126 ≠0.118 6.483 +0.803 ≠0.809 1.495 +0.096 ≠0.089 2.633 +0.901 ≠0.795 12.142 +0.606 ≠0.736 4.253 +1.099 ≠1.020 2.921 +0.195 ≠0.181 3.255 +1.113 ≠0.873	V V V	km/s 127.285 23.412 38.074 79.761 85.878 68.390 44.779 33.406 65.114 70.772 111.659 36.165 119.880 85.056 100.299 101.288 46.467 38.402 42.059 35.204 68.960 52.713 103.863 54.597 69.013 67.509 77.061 161.626 86.622 80.075 60.767 93.678 58.635 82.290 59.030 95.240 151.760 87.268 41.518 83.853 max (halo) V km/s 45.246 140.528 101.836 128.583 80.845 79.429 53.195 108.220 93.739 66.580 65.647 82.471 59.457 59.816 137.026 57.898 99.725 46.655 97.831 144.640 73.574 70.941 66.598 88.576 101.448 66.033 137.220 213.375 84.913 100.502 126.330 110.456 68.214 73.806 113.019 116.503 113.179 85.596 51.233 92.922 max (halo) V km/s 102.272 73.976 83.016 87.346 60.944 67.377 157.180 103.301 69.889 95.707 73.960 86.464 71.032 77.015 134.297 122.668 72.081 96.166 74.054 219.389 118.721 108.467 87.708 135.058 94.381 147.807 110.951 109.782 114.309 95.871 124.633 104.811 136.640 145.920 150.608 135.957 167.288 116.645 207.004 167.406 max (halo) V km/s 130.006 158.666 145.681 162.802 185.369 139.974 136.848 196.540 143.579 196.666 204.997 233.249 162.071 161.001 220.272 199.606 189.173 164.879 206.233 173.176 240.395 172.632 205.130 253.045 154.126 206.545 204.544 210.523 157.854 208.888 220.493 291.419 222.908 281.342 192.593 189.943 246.840 264.791 219.488 291.172	km/s 23.021 24.762 28.791 28.991 30.516 30.792 31.869 33.684 33.900 33.976 34.637 35.034 35.042 35.053 35.060 35.160 36.123 37.318 37.660 37.914 38.360 38.733 38.734 39.362 39.407 39.498 40.338 40.932 41.118 41.435 41.580 41.799 42.054 44.019 44.207 44.332 46.397 46.052 45.904 44.392 rot (2kpc) km/s 46.531 46.540 48.169 48.280 48.600 48.721 49.617 49.867 49.959 50.096 50.227 50.683 51.025 51.257 51.708 51.750 53.507 53.689 53.897 54.376 54.384 54.523 54.894 54.968 55.359 55.400 55.637 56.165 56.484 56.508 56.937 57.582 58.734 59.011 59.881 60.004 60.384 60.701 60.807 62.113 rot (2kpc) km/s 62.293 62.666 62.702 63.879 64.602 65.574 65.970 66.710 68.725 68.852 70.313 71.815 73.839 74.831 74.983 75.545 76.546 77.856 78.782 81.435 84.677 85.704 90.253 90.673 91.310 91.435 92.687 93.255 93.323 95.824 99.039 105.152 111.164 111.773 112.843 113.577 115.527 122.867 131.944 138.659 rot (2kpc) km/s 139.752 141.299 145.805 155.362 155.912 165.349 175.357 175.980 184.474 184.833 188.861 197.100 198.795 201.786 203.722 204.271 206.882 209.442 209.942 213.446 213.831 215.180 215.427 221.563 224.497 226.200 230.266 231.898 233.305 235.651 236.328 237.935 240.853 241.926 245.137 248.367 255.110 296.261 311.334 427.767

Article number, page 4 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Article number, page 6 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Article number, page 8 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Article number, page 10 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Article number, page 12 of 29 Amir Ghari et al.: Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves

Acknowledgements

The authors thank the referee, Kyle Oman, for excellent suggestions, and for providing us with the APOS-TLE/EAGLE V Rlast and V rot (2 kpc) simulated datapoints. We also thank Aaron Dutton for providing us with the NI-HAO RAR at 2 kpc. We acknowledge useful discussions with Federico Lelli and Stacy McGaugh.


Tim Haines
email: thaines@astro.wisc.edu

Elena D'onghia

Benoit Famaey

Chervin Laporte

Lars Hernquist

IMPLICATIONS OF A TIME-VARYING GALACTIC POTENTIAL FOR DETERMINATIONS OF THE DYNAMICAL SURFACE DENSITY

Keywords: Galaxy, kinematics and dynamics -stars, kinematics and dynamics

We analyze a high-resolution N-body simulation of a live stellar disk perturbed by the recent passage of a massive dwarf galaxy that induces a wobble, in-plane rings and phase spirals in the disk. The implications of the phase-space structures and the estimate of the matter density through traditional Jeans modeling are investigated. The dwarf satellite excites rapid time-variations in the potential, leading to a significant bias of the local matter surface density determined through such a method. In particular, while the Jeans modelling gives reasonable estimates for the surface density in the most overdense regions of the disk, we show that it fails elsewhere. Our results show that the spiralshape feature in the z -v z plane is visible preferentially in the under-dense regions and vanishes more quickly in the inner parts of the stellar disk or in high-density regions of the disk. While this prediction can be verified with Gaia DR3 our finding is highly relevant for future attempts in determining the dynamical surface density of the outer Milky Way disk as a function of radius. The outer disk of the Milky Way is indeed heavily perturbed, and Gaia DR2 data have clearly shown that such phase-space perturbations are even present locally. While our results show that traditional Jeans modelling should give reliable results in overdense regions of the disk, the important biases in underdense regions call for the development of non-equilibrium methods to estimate the dynamical matter density locally and in the outer disk.

INTRODUCTION

Thanks to the second Gaia data release (Gaia DR2) (Gaia Collaboration et al. 2018), a map of the kinematics of disk stars in the Galaxy can now be provided not only for the solar vicinity but also over a larger spatial extent out to a few kpc from the Sun, allowing us to explore the phase-space of the stellar disk. This will become even more topical when the Gaia DR2 will be combined with future large spectroscopic surveys such as 4MOST (de Jong et al. 2012) and WEAVE [START_REF] Dalton | Proc. SPIE[END_REF][START_REF] Famaey | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]. In principle, this should allow us to map the dynamical surface density as a function of height at various radii in the disk, and compare this with the baryonic density to infer the distribution of dark matter around the disk (e.g. Bovy & Rix 2013). Locally the dynamical surface mass density estimates inferred from stars oscillating vertically above and below the plane (Kuijken & Gilmore 1991;Garbari et al. 2012;Bovy & Tremaine 2012;Bienaymé et al. 2014;Hagen & Helmi 2018;Kuijken & Gilmore 1989a,b;Siebert et al. 2003;Holmberg & Flynn 2004;Bovy & Tremaine 2012;Widmark & Monari 2019;Widmark 2019) have been used to probe the local dark matter density, the dark matter vertical distribution and its possible deviations from a spherical halo profile, as well as the presence of a dark disk, thereby providing potential constraints on the nature of dark matter, and on alternatives (Bienaymé et al. 2009). In particular, if the dark matter density is very low at the solar vicinity, around the midplane, the relative contributions from visible and dark matter can be disentangled by trying to constrain the gravitational potential out to larger heights (Bovy & Tremaine 2012;Holmberg & Flynn 2004).

However, the various values reported for the dark matter density estimates in the solar neighborhood are affected by uncertainties that are associated with the large errors in the estimate of the gas mass density and the local stellar densities (see Silverwood et al. 2016). Furthermore, all these analyses assume that the disk potential is static and in dynamical equilibrium. In particular potential collective modes in the Galactic disk are not accounted for. For example Banik et al. (2017) recently studied models of a Galactic disk out of dynamical equilibrium to show that when a single tracer population is adopted to measure the vertical force at the solar vicinity, the uncertainty in the inferred value may be of the order of 20%.

Gaia DR2 data indicate that the Milky Way stellar disk in the solar neighborhood is far from being smooth but shows evidence of substructure in the kinematic space that manifests in the form of ridges and arches. While such structures had long been known for planar motions (e.g. Dehnen 1998;[START_REF] Famaey | [END_REF]Monari et al. 2017), the finding of a phase-space spiral structure reported in the z -v z plane (Antoja et al. 2018) could potentially have important implications for determinations of the dynamical surface density. This feature appears to extend well beyond the solar neighborhood cylinder in which it was originally discovered and traced out to R ∼ 11 kpc (Laporte et al. 2018b) using Gaia DR2 (see also Bland-Hawthorn et al. 2019, for a GALAH R ± 0.5 kpc exploration).

Previous observational facts already indicated that the stellar disk is out of dynamical equilibrium, namely the vertical bulk motion of stars discovered in the solar neighborhood. These motions manisfest as vertical oscillations of the stellar disk as measured by SEGUE, RAVE and LAMOST (e.g., Widrow et al. 2012;Williams et al. 2013;Xu et al. 2015;Carrillo et al. 2018;Wang et al. 2018b,a) and may be caused by the passage of a massive satellite such as the Sagittarius dwarf galaxy (Gómez et al. 2012(Gómez et al. , 2013;;Widrow et al. 2014;de la Vega et al. 2015;D'Onghia et al. 2016), although Faure et al. (2014), [START_REF] Debattista | [END_REF]Monari et al. (2016) showed that vertical oscillations in the form of breathing modes are also obtained by the passage of the stars through spiral arms. Several independent studies showed that rings, wobbles and in-plane velocity anisotropies may be the dynamical response of the stellar disk to the gravitational disturbance of a satellite galaxy (Minchev et al. 2009;Purcell et al. 2011;D'Onghia et al. 2016) or their collective effects (Kazantzidis et al. 2009;Chequers et al. 2018). In the case of a perturbation by the Sagittarius dwarf, not only do such disturbances affect the inner part of the galaxy, but they can impact the outer disk in regions where the stellar surface density decreases and the disk is more sensitive to external perturbations. In particular, the recent pre-Gaia DR2 Sgr impact models of Laporte et al. (2018a) demonstrated that it could simultaneously account for the morphology of the outer disk and amplitude of density and streaming motion fluctuations in the solar neighbourhood.

Recently, the phase-space spiral was also attributed to the buckling of the Galactic bar (Khoperskov et al. 2019) although Laporte et al. (2018b) showed that it affects stars of all ages, which is a priori more in line with a recent perturbation. It is however possible that perturbations linked to the last pericentric passage of the Sagittarius dwarf and to the buckling of the bar could co-exist in the Milky Way disk.

The extension of the spatial scale of kinematic studies to larger regions of the Galactic disk by Gaia DR2 indicated the evidence of streaming motions in the radial, azimuthal and vertical velocity as well as small-amplitude fluctuations in the velocity dispersions in the region between 3 and 8 kpc far from the Solar location. In particular, the data suggest an increase in median vertical velocity, from the inner to the outer disk, that shows a complex dependence on the azimuthal component (Katz et al. 2018).

In this Letter we use an N-body simulation of a Milky Way-sized galaxy to investigate the implications of timevarying potentials for the phase-space structure and the surface density of the stellar disk induced by the impact of a massive satellite galaxy. We show that the assumptions of dynamical equilibrium and a static potential used to infer the disk dynamical density for different azimuthal angles, and particularly at large distances from galaxy center, based on stellar kinematics may not be robust unless the time-varying effects are taken into account. This paper starts with a brief description of the numerical experiment used in this study and its peculiar phasespace structure (Sec. 2). We then describe the classical Jeans approach that one would be allowed to follow in dynamical equilibrium (Sec. 3), and show how it could mislead an observer in deducing the incorrect dynamical surface density of the disk when the potential is actually time-varying.

THE NUMERICAL SIMULATION

Set-up

We analyse a high-resolution N-body simulation of a Milky Way sized galaxy performed with the parallel TreePM code GADGET-3. A detailed description of the simulation is available in the literature (Laporte et al. 2018a). In this study the galaxy model consists of a live dark matter halo, a rotationally supported disk of stars, and a spherical stellar bulge. The dark matter mass distribution of the galaxy is initially modeled with a Hernquist profile (Hernquist 1990). The total mass of the halo is 1 × 10 12 M and is sampled with 40 million particles. The radial scale length of the halo is 52kpc prior to applying adiabatic contraction according to (Blumenthal et al. 1986). Each dark particle has a mass of 2.6×10 4 M and a gravitational softening length of 60 pc. For the stellar disk we adopt an exponential radial profile with an isothermal vertical distribution:

ρ * (R, z) = Σ(R)ζ(z) = M 4πz 0 h 2 R exp (-R/h R ) sech 2 z z 0 . ( 1
)
where M is the stellar disk mass, h R is the exponential disk scale length and z 0 is the scale-height of the sech 2 vertical distribution. The disk mass is set to M = 6.0 × 10 10 M , h R = 3.5 kpc, and z 0 = 0.53 kpc.

The stellar disk is discretized with 5 million particles.

The radial and vertical stellar velocity dispersionsquared both follow an exponential profile with the same scale-length as the stellar density h σ = 3.5 kpc. The stellar particles have masses of 1.2×10 4 M and gravitational softening lengths of 30 pc. The stellar bulge is described by the Hernquist model, with a total mass 1 × 10 10 M and scale radius 0.7 kpc.

Interaction with a massive satellite

The satellite adopted in the simulation (Model L2 of Laporte et al. 2018a) has a virial mass of 6×10 10 M and sinks through dynamical friction for 6 Gyr into the host galaxy halo. It excites a wake in the dark halo, creating torques acting on the disk, before transitioning to acting directly on the disk through tides at late times. As shown in Laporte et al. (2018b), each pericentric passage excites a new phase-space spiral in the z -v z plane, similar to the one observed in the Solar neighbourhood with Gaia. Each impact induces a prompt response of the disk after the collision, with the stars around the impact region moving inward first and then outwards. A wave of enhanced density is the outcome of this contraction and expansion. Such a wave that propagates through the disk tends to form a two-armed ring-like pattern with over-dense and under-dense regions together with strong bending vertical velocities (see e.g. Gómez et al. 2013;D'Onghia et al. 2016).

Hereafter, we are studying the problem of performing Jeans modelling on such a galaxy which has generically reacted to a satellite perturbation, in a general fashion, and have chosen a random time of t=6 Gyr for the analysis. The mass model of the reacting galaxy is close to that of the Milky Way in its gross properties, but our Jeans modelling should not be quantitatively compared to Gaia data, since the model is not meant to be a perfect realization of the present-day Milky Way. However, as already shown in Laporte et al. (2018a,b), a number of qualitative and quantitative agreements with the data exist. Here, we divide the disk into six sextants, and display in Fig. 1 the breathing and bending velocities in layers of 400 pc height along the vertical component of the disk. Following Katz et al. (2018), these velocities are computed as the half-sum (mean) and half-difference of the median vertical velocities in symmetric layers with respect to the Galactic mid-plane. The breathing (top row) and bending (bottom row) mode velocities within 8 ± 4kpc at t= 6.0 Gyrs.

With the above caveat in mind, if we locate the Sun at the sextant with the azimuth angle between 60 < θ < 0 (which corresponds also to an under-dense region in the disk) the bending velocities are negative (i.e. they are oriented towards the south Galactic pole) for |z| between [0,400] pc and [400, 800] pc but become positive , in qualitative agreement with the Gaia DR2 data. Furthermore, the breathing velocity changes from positive to negative values at larger heights from the midplane, again in qualitative agreement with the current observations.

At the Solar vicinity, the unprecedented precision and the high-sampling offered by the Gaia DR2 data showed that the stars are distributed in the plane defined by the vertical coordinate and the vertical velocity, z -v z with a curled spiral-shape whose density increases towards the leading edge of the feature (Antoja et al. 2018). Such a spiral shape is also clearly visible when mapping the average azimuthal velocity as a function of position in the z -v z surface of section.

This feature is interpreted as the incomplete phasemixing phenomenon occurring in the stellar disk and a manifestation of the disk being locally out of dynamical equilibrium. Recent studies using analytic models (Antoja et al. 2018;Binney & Schönrich 2018) or numerical simulations (D'Onghia et al. 2016;Laporte et al. 2018a;[START_REF] Dalton | Proc. SPIE[END_REF]) interpreted this feature as likely caused by the coupling of vertical and in-plane motions induced by the passage of a satellite galaxy through the Galactic disk. In fact in the presence of rings induced by the impact of a satellite, the intra-ring underdense regions coincide with a local increase of the characteristic scale height (see Fig. 1 in D'Onghia et al. 2016). This is consistent also with the fact that the feature is not traced by the OB stars, which are younger than 100 Myrs and did not complete one orbital period around the Galactic center (Cheng et al. 2019). In normal conditions for disk stars, the vertical motion of the stars is almost fully decoupled from the in-plane motion with the vertical frequency of stars being usually much greater than the radial epicyclic frequency, with the stars oscillating vertically much faster than their radial motions. However, in under-dense regions, stars feel a weaker gravitational restoring force, and have their vertical frequency reduced. Thus, we expect that stars located in the underdense regions of the disk will have the horizontal and vertical motions coupled. Hence the spiral-shape feature in the z -v z plane should be visible preferentially in the under-dense regions and should vanish more quickly in the inner parts of the stellar disk or in locations with enhanced density like density waves.

Fig. 2 displays the azimuthal velocity in the z -v z plane in one of the most underdense regions of the simulated disk at t=6 Gyr, defined by the sextan with the azimuthal angle ranging -120 < θ < -60 (right panel) as compared to the same analysis reported in the enhanced density region in the sextant characterized by -180 < θ < -120 (left panel). Indeed, as expected, the spiral-shape feature is visible in the under-dense region, but has a much less well-defined shape in the over-dense region. Note that this finding suggests that the disk self-gravity for time-varying potentials plays a key role in the interpretation of the spiral feature. Therefore, approaches based merely on the impulse approximation and the phase-wrapping of the stellar orbits to describe the feature are limited.

THE JEANS ANALYSIS

We now proceed with analyzing the simulation at two different times, namely (i) t = 0.75 Gyr, after the disk has settled from the initial conditions, but before the massive satellite has had a significant impact on it, and (ii) t = 6 Gyr, an arbitrary time-step 100 Myr before the last pericentric passage of the massive satellite.

In both cases, we make a vertical Jeans analysis at R = 8 kpc by constructing an annulus 200 pc in width at this radius. This annulus is then broken into equallysized azimuthal patches that approximately correspond to the view used by Katz et al. (2018) to analyze the vertical velocity field of the Gaia DR2 data (see Fig. 1).

As is traditionally done, we then combine the vertical Jeans equation with the Poisson equation to obtain an estimate of the dynamical surface density as a function of height at R = 8 kpc. This Jeans estimate can then directly be compared with the actual surface density of stars and dark matter in the simulation.

The equation we use for the Jeans estimate of the dynamical surface density at R = 8 kpc is (see, e.g., Hagen & Helmi 2018)

2πGΣ dyn (R, z) σ 2 z h z - ∂σ 2 z ∂z (2) -cov (v R , v z ) 1 R - 2 h R +2v c |z| R ∂v c
∂R where the three first terms correspond to the estimate of the vertical force from the vertical Jeans equation and the last term accounts for a non-flat rotation curve at the radius of interest, assuming that it does not vary significantly within the height of interest. Here, h z = -(dln(ρ)/dz) -1 is just a notation corresponding to the local "exponential scale-height" at the height of interest, not to be confused with the (constant) scale-height z 0 of Eq.( 1). We note that h z is calculated independently in each azimuthal bin, as well. The vertical density profile is not exponential. Note that this equation is based on a few important assumptions, notably that the radial and vertical velocity dispersion-squared vary with radius with the same exponential scale length as the density, and that the velocity ellipsoid tilt angle is locally constant with radius. The assumption that the dispersion-squared follows the same radial exponential profile is justified by the set-up of the initial conditions. As we shall see below, this set of assumptions indeed yields reasonable results when the disk is unperturbed (Fig. 3). However, an obvious other assumption is that the disc is axisymmetric indeed. In this case, while the results at t = 0.75 Gyr should logically not be affected, it is not completely clear how much this assumption will affect the results at t = 6 Gyr (see Sects. 4 and 5).

Within each azimuthal patch, we recenter the stars by subtracting the position of the peak vertical density. The stars are then split into vertical bins of equal counts to provide an approximately equal signal-to-noise in each bin when measuring the velocity dispersion. The vertical density is fit by a sech 2 distribution, which allows us to compute h z (z). At t = 0.75 Gyr, the simulated disk is nicely plane-symmetric, which allows us to fit a single sech 2 to the whole distribution and test our framework by computing the total density inferred from Eq.( 2), Σ dyn , and the actual total density, and compare it to Σ tot = Σ + Σ DM . The results are shown in Fig. 3, and demonstrate that the method works reasonably well when the simulated disk is close to equilibrium.

At t = 6 Gyr, even after correcting for the position of the peak vertical density, the important North-South asymmetries in the vertical distribution of stars force us (as an observer confronted to such asymmetries would be) to consider a separate Jeans analysis in the Northern and Southern Galactic hemispheres.

RESULTS

Fig. 4 shows the cumulative stellar surface mass density (SMD) (blue), the cumulative total matter SMD (stars + DM; orange), the SMD derived from Eq. 2 using the stellar kinematics (green at z > 0 and red at z < 0) as described in Section 2.1 as a function of height above the plane.

We note that, in the most over-dense regions of the stellar disk (the upper left and bottom right panels), the matter density dynamically inferred by the Jeans analysis tends to be consistent in both the Northern and Southern hemispheres, and is consistent with the total matter density present in the simulations (albeit with a slight underestimation at large heights in the sextant -180 < θ < 120). On the other hand, in the most under-dense regions of the stellar disk (the left middle and bottom panels) the matter surface density inferred from the Jeans analysis is also quite consistent in both hemispheres, but systematically overestimates the total surface density. This systematic overestimation is in line with the previous analysis of a simpler model of a disk affected by breathing modes (see Figure 6 in Banik et al. 2017). Finally, in the regions of average density, the Northern and Southern estimates tend to disagree with each other, but their average is close to the true value. For patches where the North and South Jeans analyses give consistent answers, one might wonder whether the assumptions of Eq. 2 might be the culprits. The stellar density scale-length is not modified at 6 Gyr compared to the initial conditions, once averaging over the wiggles due to azimuthal over-and under-densities. The velocity dispersion-squared scale length is largely unmodified too, although it varies slightly from one azimuth to the other, between one and twice the density scale-length. This, together with the assumption on the tilt angle (see, e.g., Hagen et al. 2019;Everall et al. 2019, for recent discussions), could modify the results slightly. Moreover, due to the clear non-axisymmetry at 6 Gyr, an axial term should also be included (Silverwood et al. 2016). However, we expect all these terms to be largely subdominant with respect to the first term of the r.h.s. of Eq. ( 2), which largely dominates the result we obtained, as the last two terms only account for modifications at the level of ∼ 1M /pc 2 . Indeed, we can interpret the overestimate of the dynamical surface density in underdense regions mostly in terms of h z = (dln(ρ)/dz) -1 , which is too small compared to its equilibrium value in those underdense regions. Moreover, although unlikely, even if a thorough non-axisymmetric analysis would be able to recover the correct surface density, it could never explain the large North-South differences found in 2 out of 6 sextans, which, if found in a data analysis, would be a smoking gun of disequi- librium. Interestingly, the actual stellar densities are actually rather symmetric in the North and South (with differences only of the order of ∼10 M /pc 2 in those patches where the derived surface densities are strongly asymmetric.

CONCLUSIONS

In this work, we extended the classical Jeans analysis to a simulated stellar disk impacted by a massive satellite. This collision induces a vertical wobble in the disk and in-plane rings, features that are not accounted for in the estimate of the local matter density inferred from the observed kinematics of the stars. Our analysis is a simple zeroth-order one, as is often performed in the solar neighbourhood, without using nonparametric approaches or taking into account the possible effects of non-axisymmetry (Widmark & Monari 2019;Widmark 2019;Silverwood et al. 2016).

Our main finding is that, while this simple Jeans modelling gives reasonable estimates of the total surface density in the most over-dense regions of the disk, it tends to systematically overestimate the surface density, by a factor as large as 1.5, in the most under-dense regions, where the scale-height is too small compared to its equilibrium value. This reinforces our interpretation that the vertical motion of the stars is coupled to their in-plane motion in these under-dense regions, thereby rendering the equilibrium analysis obsolete. We also point out that, in the regions of average density, the Jeans analysis tends to give very discrepant results in the Northern and Southern Galactic hemispheres, whose average value is nevertheless close to the true one. Such North-South asymmetries are a smoking gun of disequilibrium, as a thorough non-axisymmetric analysis, by construction, could never account for such asymmetries.

Our findings are highly relevant for future attempts at determining the dynamical surface density of the outer Milky Way disk as a function of radius, by combining Gaia data with future spectroscopic surveys. The important biases of simple Jeans modelling, especially in under-dense regions call for the development of nonequilibrium methods to estimate the dynamical matter density in the outer disk. For instance, the method recently proposed by Binney & Schönrich (2018) could make it possible to take advantage of the presence of a phase-space spiral to constrain the potential once the timing of the main perturber (for instance the last pericentric passage of the Sgr dwarf) is known. However, we note that such methods rely on the impulse approximation, whereas yet again indicated by our results the selfgravity of the disk is relevant, meaning that the method should be improved to take it into account.

  Figure 1.The breathing (top row) and bending (bottom row) mode velocities within 8 ± 4kpc at t= 6.0 Gyrs.

Figure 2 .

 2 Figure 2. Two-dimensional surface of section in bins of ∆z = 0.05kpc and ∆vz = 1km/s colored as a function of the average v φ for stars within 500pc of R = 8kpc and -1.5 ≤ z ≤ 1.5kpc at t=6Gyrs in the -180 ≤ θ < -120 azimuthal bin (left) and the -120 ≤ θ < -60 bin (right) smoothed using a Gaussian with zero mean and σ = 1.

Figure 3 .

 3 Figure 3. Cumulative stellar surface mass density(blue), cumulative total density (stars and dark matter; orange), and the density derived from Eqn. 2 (grey) as a function of height above the midplane where the lower half-plane has been folded onto the upper half-plane in each azimuthal bin at t=0.75 Gyrs.

Figure 4 .

 4 Figure 4. Cumulative stellar surface mass density(blue), cumulative total density (stars and dark matter; orange), and the density derived from Eq. 2 for z< 0kpc (brown-red) and z> 0kpc (green) as a function of height above the midplane in each azimuthal bin at t=6 Gyrs.

ACKNOWLEDGEMENTS

We thank the anonymous referee for a thoughtful and constructive report. We are grateful to Justin Read and George Lake for their valuable advice. E.D.O acknowledges support from the Vilas Associate Research Fellowship and thanks the Institute for Theory and Computation (ITC) at Harvard and the Center for Computational Astrophysics at the Flatiron Institute for the hospitality during the completion of this work. CL is supported by a CITA National Fellow award. This work is supported by ATP NASA Grant No NNX144AP53G. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575. E.D.O, B.F and C.L acknowledge hospitality at the KITP, supported by the National Science Foundation under Grant No. NSF PHY-1748958, in the final stages of this work.