EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 653 (September 2021) A&A, 653 (2021) A86 Full HTML
Free Access
Issue
A&A
Volume 653, September 2021
Article Number A86
Number of page(s) 19
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202141466
Published online 15 September 2021

© ESO 2021

1. Introduction

Our knowledge about the dynamics, composition, and history of the Milky Way is intimately connected to the determination of its gravitational potential (Dehnen & Binney 1998; Klypin et al. 2002; Widrow et al. 2008; Kafle et al. 2014; Cole & Binney 2017; McMillan 2017; Nitschai et al. 2020; Cautun et al. 2020; Li et al. 2020). Furthermore, direct and indirect dark matter detection experiments rely on precise knowledge on how dark matter is distributed in our Galaxy (Stoehr et al. 2003; Vogelsberger et al. 2009; Klasen et al. 2015; Petač 2020; Nobile 2021). The gravitational potential of the Galaxy is often inferred by fitting a stellar number density distribution to data under the assumption of a steady state, either in the solar neighbourhood (McKee et al. 2015; Widmark & Monari 2018; Sivertsson et al. 2018; Schutz et al. 2018; Buch et al. 2019; Guo et al. 2020; Salomon et al. 2020; Li & Widrow 2021) or a more global spatial volume (Nesti & Salucci 2013; Cole & Binney 2017; McMillan 2017; Cautun et al. 2020; Hattori et al. 2020; Petač 2020). Other studies have explored a measurement of the gravitational potential directly from stellar accelerations (Chakrabarti et al. 2021; Buschmann et al. 2021) or from stellar streams (Koposov et al. 2010; Bonaca & Hogg 2018; Malhan & Ibata 2019; Widmark et al. 2020).

Gravitational probes of the Milky Way can help constrain the particle nature of dark matter through the detection or exclusion of dark substructures, such as the subhalos predicted in the cold dark matter scenario (Diemand et al. 2008; Springel et al. 2008; Stref & Lavalle 2017; Facchinetti et al. 2020). There has also been speculation that a dark disk, co-planar with the stellar disk, could exist in the Milky Way, formed either from the accretion of satellites (Read et al. 2008; Purcell et al. 2009; Ruchti et al. 2014) or, in a less standard scenario, from a dark matter particle sub-species with strong dissipative self-interactions (Fan et al. 2013a,b). This latter type of dark disk could potentially be very thin, with a scale height as small as a few tens of parsecs. Such a disk has been constrained in dynamical mass measurements of the solar neighbourhood (Kramer & Randall 2016b; Caputo et al. 2018; Schutz et al. 2018; Buch et al. 2019), where the strongest constraints use Gaia observations of the very local population of stars. However, in a similar study by Widmark et al. (2021a) (also drawing from the results of Widmark 2019), they argue that such measurements are strongly biased by time-varying dynamical effects.

Although disk formation requires a fairly quiescent accretion history (Freeman & Bland-Hawthorn 2002), non-equilibrium effects are certainly not lacking in the Milky Way. Over the last decades, a growing body of observations, partly enabled by the advent of large automated Milky Way surveys, have revealed signs of large-scale asymmetries in the vertical structure in the solar neighbourhood (Widrow et al. 2012) and the outskirts of the Galactic disk (Newberg et al. 2002), as well as in its velocity distribution (Widrow et al. 2012; Williams et al. 2013; Carlin et al. 2013). This suggests the presence of bending waves in the disk (Xu et al. 2015; Price-Whelan et al. 2015; Sheffield et al. 2018; Bergemann et al. 2018), possibly sourced by self-excitation (Chequers & Widrow 2017), satellite and dark subhalo interactions (Kazantzidis et al. 2008; Gómez et al. 2013; Widrow et al. 2014; Chequers et al. 2018), fly-bys (Gómez et al. 2016), or even bar buckling (Khoperskov et al. 2019). The Gaia satellite reaffirmed the presence of these asymmetries and helped characterise the wavelength of these perturbations (e.g., Schönrich & Dehnen 2018), as well as their time-evolving nature. Using the second data release from Gaia, Antoja et al. (2018) revealed the presence of a phase-space spiral in the solar neighbourhood. The presence of the phase-space spiral at all stellar ages, as shown by Laporte et al. (2019), indicates a recent collective perturbation of the disk. Furthermore, the spiral’s presence over more than two disk scale lengths (Laporte et al. 2019) and as far as 14 kpc from the Galactic centre (Xu et al. 2020) confirms it being a manifestation of a disk-wide disequilibrium phenomenon. This reinforces the connection between the perturbations seen in the solar neighbourhood and the Galactic outskirts, as was originally predicted by pre-Gaia models of satellite interactions (Laporte et al. 2018), such as the Sagittarius dwarf galaxy. Indeed, using toy models, Antoja et al. (2018) and Binney & Schönrich (2018) noted that the perturbation may have been set off as far back as 900 Myr ago, which is similar to the orbital period of the Sagittarius dwarf galaxy derived from stream fitting models (e.g., Johnston et al. 2005; Vasiliev et al. 2021). Although the phase-space spiral was originally identified and mostly studied in its azimuthal and radial velocity moments (e.g., Antoja et al. 2018; Binney & Schönrich 2018; Khoperskov et al. 2019; Bland-Hawthorn et al. 2019), we exploit the spiral’s shape as seen in terms of its relative number density with respect to the stellar bulk background, as first revealed in Laporte et al. (2019). This reduces the exercise of potential fitting to only the vertical dimension, similar to what is typically done when determining the gravitational potential and dark matter density in the solar neighbourhood (Read 2014; de Salas & Widmark 2020).

The phase-space spiral is a time-varying dynamical structure and as such it constitutes an obstacle and a systematic bias to traditional dynamical mass measurement methods that assume a steady state (e.g., Jeans modelling). The method employed in this work is complementary to such traditional methods, in the sense that it extracts information from the shape of the phase-space spiral itself. It does so under the assumption that the winding angle of the spiral is a smooth function with respect to vertical energy; given the shape of spiral in the (z, w)-plane, this sets strong constraints on the vertical gravitational potential. The general principles of our method are discussed at length in the first paper of this series—Widmark et al. (2021b), henceforth referred to as Paper I—where we also tested our method on one-dimensional simulations. In those tests, we were able to retrieve the true gravitational potentials of our simulations with high accuracy.

In this work, we applied our method to the Milky Way, using the early instalment of Gaia’s third data release (EDR3), supplemented with radial velocity measurements from legacy spectroscopic surveys through the catalogue compiled by Sanders & Das (2018). We constructed eleven main data samples using different cuts in Galactocentric radius and angular momentum. We were able to measure the gravitational potential to high precision and, using a model for the baryonic densities, we inferred the local halo dark matter density and placed the most stringent constraints on the surface density of a thin dark disk.

This article is structured as follows. We begin with some definitions in Sect. 2, such as the coordinate system, continuing with a description of the data in Sect. 3. In Sects. 4 and 5, we outline our model of inference and our model for the baryonic matter densities in the solar neighbourhood. In Sect. 6, we present our results. Finally, we discuss and conclude in Sects. 7 and 8.

2. Coordinate system and other definitions

In this paper, we use the following system of coordinates. The spatial coordinates X ≡ {X, Y, Z} denote the position with respect to the Sun, where positive X corresponds to the direction of the Galactic centre, positive Y corresponds to the direction of Galactic rotation, and positive Z corresponds to the direction of Galactic north. Their respective time derivatives correspond to velocities V ≡ {U, V, W} in the solar rest frame.

The spatial coordinate Z, denoting height with respect to the Sun, is related to the height with respect to the Galactic disk, written z, according to

(1)

where Z⊙ is the Sun’s height with respect to the Galactic mid-plane. Similarly, the velocities in the Local Standard of Rest (see Binney & Tremaine 2008), written v ≡ {u, v, w}, are found via the Sun’s peculiar motion, written V⊙ ≡ {U⊙, V⊙, W⊙}, according to

(2)

We adopt values of V⊙ = {11.1, 12.24, 7.25} km s−1 (Schönrich et al. 2010).

The Poisson equation is written

(3)

where we have neglected the contribution in the azimuthal direction (which is zero in the case of rotational symmetry). The second term of the left-hand side is known as the circular velocity term, because the circular velocity of the Galactic plane vc is given by

(4)

This term can be written as a matter density correction, as

(5)

Because the rotational velocity curve is close to flat, this correction is small; the slope of the circular velocity is roughly ∂vc/∂R = −1.5 ± 0.2 km s−1 kpc−1 (−1.7 ± 0.1 km s−1 kpc−1 in Eilers et al. 2019; −1.33 ± 0.1 km s−1 kpc−1 in Ablimit et al. 2020), giving Δρ ≃ −0.0016 ± 0.0006 M⊙ pc−3.

In our model of inference, the gravitational potential is equal to

(6)

where the free parameters ρh = {1, 2, 3, 4} are constrained to lie in the range [0, 0.2] M⊙ pc−3. Via the Poisson equation of Eq. (3) and also accounting for the matter density correction of Eq. (5), the total matter density is equal to

(7)

Using this functional form, the gravitational potential and matter density distribution are very free to vary in shape (although assumed to be symmetric, smooth, and strictly decreasing with |z|). They are not constrained by strong prior information about the baryonic matter density distribution in the Galactic disk, but are flexible enough to emulate such models (see e.g., McKee et al. 2015 and Schutz et al. 2018).

3. Data

In this work we have used data from Gaia EDR3, supplemented with radial velocity information compiled in Sanders & Das (2018) using LAMOST DR3 (Deng et al. 2012), GALAH DR2 (Buder et al. 2018), RAVE DR5 (Kunder et al. 2017), APOGEE DR14 (Abolfathi et al. 2018), SEGUE (Yanny et al. 2009), and GES DR3 (Gilmore et al. 2012). A supplementary radial velocity was used if a star in the Gaia catalogue had missing radial velocity information or had a radial velocity uncertainty larger than 3 km s−1. If there was radial velocity information from several supplementary surveys, we used the value with the smallest associated uncertainty. Additionally, we excluded any star with discrepant radial velocity measurements, requiring less than a 2.5σ tension between all supplementary surveys that had a radial velocity uncertainty smaller than 5 km s−1; this amounted to removing roughly one per cent of the stars with supplementary radial velocity measurements. The stellar number counts of the radial velocity information taken from the respective surveys are listed in Table 1.

Table 1.

Stellar number counts of radial velocity information for our eleven main data samples.

We constructed eleven main data samples on which we applied our method. To begin with, we made cuts in data quality, which all samples were subjected to. After that, we constructed our respective samples by making cuts in phase-space. These cuts are discussed in detail below.

3.1. Data quality cuts

In this work, we applied quite strong constraints with respect to data quality. This was done in order to be able to neglect observational uncertainties and utilise the stars’ six-dimensional phase-space information. At the same time, this gave rise to strong selection effects, most importantly with regards to radial velocity completeness, which has a strong spatial dependence. However, strong selection effects are not detrimental to our method, as long as the shape of the phase-space spiral is robustly extracted (this is discussed further in the beginning of Sect. 4.1).

We required that the stars in our data sets would have a Gaia G-band magnitude smaller than 15 mag (in Gaia EDR3 there are only spurious stars with radial velocity information above 15 mag, see Gaia Collaboration 2021a), and a radial velocity uncertainty (σRV) smaller than 3 km s−1. For the radial velocity measurements coming from the Gaia spectrograph, the uncertainties are typically around 0.3 km s−1 for bright stars and 1.8 km s−1 for dim stars (Katz et al. 2019), with only a thin tail of stars with uncertainties larger than 3 km s−1; as such, our cut in radial velocity uncertainty only removed a small fraction of stars from our data samples.

In terms of the astrometric measurements, we required that the renormalised unit weight error (RUWE) would be smaller than 1.4 and that the parallax uncertainty (σϖ) would be smaller than 0.05 mas. The first of these criteria is already very constraining, such that applying the subsequent second criteria removed only one in a thousand stars.

In summary, the data quality cuts can be written as follows:

(8)

3.2. Phase-space cuts

For the main data samples, we made the following phase-space cuts. We constructed bins in Galactocentric radius that were 100 pc wide, labelled by the index s, and also restricted the spatial extent in the azimuthal direction to |Y|≤400 pc. In this spatial volume, we made cuts in angular momentum, requiring it to be close (±10%) to that of a circular orbit. The angular momentum of a star is equal to

(9)

where vϕ is the velocity in the azimuthal direction in the Galactic rest frame. The angular momentum of a perfectly circular orbit was calculated by taking the data sample’s Galactocentric mid-point and multiplying it by the circular velocity, according to

(10)

assuming a flat rotation curve with a circular velocity of vc = 240 km s−1 (Reid et al. 2014) and a solar Galactocentric radius of R⊙ = 8178 pc (Gravity Collaboration 2019). In Fig. 1, we show the histogram of the angular momentum relative to the angular momentum of a circular orbit. The Lz distribution is not perfectly centred on Lz, circ, but rather slightly shifted towards smaller values, as expected due to asymmetric drift (Binney & Tremaine 2008).

thumbnail Fig. 1.

Histogram of angular momentum (Lz) relative to the angular momentum of a circular orbit (Lz, circ.), for different cuts in Galactocentric radius, where R − R⊙ ∈ [100s − 50, 100s + 50] pc. The grey band highlights the region where Lz is within 10% of Lz, circ..

When making these cuts in data, we neglected any observational uncertainties associated with the measurements. For example, a star’s coordinate Y is taken directly from its angular position on the sky and its parallax value, according to Y = cos(l)cos(b) (mas/ϖ) kpc. This is motivated by our restrictive cuts in data quality; for example, with this parallax precision, the distance is known to a relative uncertainty of only a few percent at the furthest relevant distance of one kilo-parsec.

In summary, the phase-space cuts of the main data samples can be written:

(11)

where s is an integer in range −5 to 5.

To validate our choice of phase-space cuts, we also applied our method to a variety of data samples using different cuts. In particular, we include the results from an additional data sample where we only make a spatial cut according to , with no cut in angular momentum. The total number of stars in our respective data samples is listed in Table 2, counting only the stars for which |Z|< 800 pc.

Table 2.

Total number count of our main eleven data samples summed together, with the additional constraints that |Z|< 800 pc.

3.3. Data reduction

In our model of inference, the data was in the reduced form of a two-dimensional histogram in the (Z, W)-plane. This histogram is the number count of observed stars in bins of size (20 pc) × (1 km s−1), written di, j, where the indices (i, j) label the respective bins. Similar to the phase-space cuts defined in Sect. 3.2, the values for Z and W are taken directly from the astrometric measurements, neglecting any observational uncertainties.

4. Model of inference

The method used in this work is the same as the one presented in Paper I, although with some minor modifications, which are as follows: (i) we constrain ourselves to an asymmetric single-armed spiral; (ii) the spiral is fitted to slightly lower vertical energies; (iii) we mask the data for latitudes lower than |b|≲20°; (iv) the Sun’s vertical position and vertical velocity are fixed rather than free parameters. Our method is described below, with emphasis on its modifications with respect to the original version.

Just as for the tests run in Paper I, we assume separability of the gravitational potential and reduce the dynamics to only the vertical dimension. In Jeans analysis, the so-called tilt term accounts for the coupling between the radial and vertical directions, and is a derivative of the bulk phase-space density distribution with respect to the radial direction. Because our method disregards the bulk density, no such term enters our method. The radial-vertical coupling of the gravitational potential does induce a radial motion, but for the relevant range of vertical energies (Ez < Φ(700 pc)) this motion is small compared to the radial extent of our data samples. In principle, this can still become relevant to the extent that the shape of the spiral has a radial dependence. Because the inferred gravitational potentials and phase-space spiral shapes are so similar between neighbouring data samples (see Sect. 6 for more details), we deem such effects to be negligible, which motivates the assumption of vertical separability when using this method.

4.1. Bulk and spiral phase-space densities

In our method, we fit a phase-space distribution to data, consisting of a product of bulk and a spiral phase-space density distribution. In traditional methods that are based on the assumption of a steady state, the bulk density is the quantity that is used to infer the gravitational potential. Conversely, in our method the bulk is only fitted in order to extract the spiral shape and does not influence, nor is influenced by, the inferred gravitational potential. Because the bulk is fitted as a mere background, to a large extent it absorbs any selection effects pertinent to the data. Such effects are critical to account for in steady state modelling, but can be disregarded in our method as long as the shape of the spiral is robustly extracted.

The free parameters of our model are shown in Table 3. The parameters are split into two groups: the bulk phase-space density parameters, written Ψbulk, and the spiral phase-space density parameters, written Ψspiral. The total number of free parameters is equal to 7 + 3K, where K is the number of Gaussian components in the bulk stellar density, for which we set K = 6.

Table 3.

Free parameters in our model of inference.

We model the bulk density as a Gaussians mixture model according to

(12)

where Ψbulk are the bulk density components, which includes the Gaussian weights (ak) and dispersions (σz, k, σw, k).

A star’s vertical oscillation has a total time period of

(13)

where Ez = Φ(z | ρh)+w2/2 is a star’s vertical energy per mass, zmax is the maximum height that a star reaches, and ρh parametrises the gravitational potential according to Eq. (6). The position of a star in the (z, w)-plane is associated with a temporal angle, which is given by

(14)

In our analytic spiral model, the initial perturbation is assumed to have no initial winding. Because the gravitational potential is anharmonic, this perturbation winds into a spiral with time. Even though the initial perturbation might have a more complicated form, the winding behaviour can be expected to dominate the shape of the eventual spiral, as long as the perturbation is not a spiral-resembling shape to begin with (a more thorough discussion on the underlying assumptions of our spiral model can be found in Paper I). The spiral angle as a function of vertical energy Ez evolves according to

(15)

where is the initial angle of the perturbation.

The total phase-space density of our analytical model is equal to

(16)

In this expression,

(17)

is the relative number density of the spiral with respect to the bulk, where α is a unit-less amplitude in range [0, 1]. This is the first difference with respect to how the method was formulated in Paper I, which included a symmetric spiral component with an amplitude β. This component was not included here, as no second arm is seen for the Milky Way spiral. Furthermore, Eq. (16) contains the quantity

(18)

where

(19)

This sets a lower boundary in Ez to the spiral of our analytic model; close to the origin of the (z, w)-plane, the spiral is washed out due to self-gravity effects. The numerical values of the m(z, w | ρh) function constitute our second modification with respect to the method as formulated in Paper I, in that the inner boundary is somewhat less restrictive, using a limit of Ez ≳ Φ(300 pc) instead of Ez ≳ Φ(400 pc). The reason for this change is that the phase-space spiral of the actual Milky Way is more clearly defined in the inner region than for the one-dimensional simulations of Paper I, possibly because the effects of self-gravity are not as strong and cohesive for the more complex, three-dimensional kinematics of our Galaxy.

4.2. Data likelihood and masks

In our method, the phase-space density model was fitted to data in two separate steps. In the first step, we fitted the bulk density distribution without the spiral; in other words, we minimised the data likelihood with respect to Ψbulk while α = 0. In the second step, we fitted the relative phase-space density spiral; in other words, we minimised the likelihood with respect to Ψspiral while Ψbulk remained fixed.

The data likelihood is given by the Poisson count comparison of the model and data in the (Z, W)-plane, in bins labelled by the indices (i, j). The logarithm of the likelihood is equal to

(20)

where and M(Zi, Wj) are mask functions described below, and f(Zi + Z⊙, Wj + W⊙ | Ψ) is the model phase-space density as defined in Eq. (16).

The first mask function in Eq. (20) is a Heaviside step function equal to

(21)

where . This mask is applied to stars with low heights, roughly corresponding to a cut in Galactic latitude of |b|< 20°. This constitutes the third modification with respect to how the method was formulated in Paper I. It is included in this work due to severe selection effects, especially in terms of the availability of radial velocity information close to the Galactic plane (see Sect. 6 and Fig. 4).

The second mask function in Eq. (20) is M(Z, W), defining a circular outer boundary in the (Z, W)-plane. It is applied in order to mask the high vertical energies where the stellar number density is low and the spiral is less pronounced. It is defined

(22)

where sigm is the sigmoid function of Eq. (19). The numerical values of Zlim. and Wlim. differed for the two minimisation steps of our method, in order to avoid fitting artefacts. In the first step, when the bulk was fitted, we set Zlim. = 800 pc and Wlim. = 44 km s−1. In the second step, when the spiral was fitted, we set Zlim. = 700 pc and Wlim. = 40 km s−1.

The fourth and final modification of our method with respect to Paper I is that we fixed the Sun’s vertical position and velocity (Z⊙ and W⊙), rather than letting them be free parameters in our fitting procedure. This is motivated by significant and spatially dependent selection effects, coming from radial velocity completeness and other data quality cuts; completeness is poor close to the Galactic mid-plane and also has a strong asymmetry with respect to the Galactic north and south for many data samples. These issues made our inference of Z⊙ and W⊙ highly unstable and they are better determined in other studies. For the Sun’s vertical velocity, we set W⊙ = 7.25 km s−1, which is a well established value, known to within an uncertainty smaller than 0.05 km s−1 (Schönrich et al. 2010). Within this uncertainty, varying W⊙ has a negligible effect on our results. For the Sun’s vertical position, the uncertainty is more significant. Different studies have produced somewhat discrepant results, ranging from roughly 0–20 pc (Jurić et al. 2008; Yao et al. 2017; Bennett & Bovy 2019; Buch et al. 2019; Widmark et al. 2021a; Gaia Collaboration 2021b). For this reason, we adopt three different values of Z⊙ = {0, 10, 20} pc when applying our method, which allows us to estimate the uncertainty of our result with respect to the Sun’s position.

Our method was implemented in TENSORFLOW, allowing for efficient minimisation using the Adam optimiser (Kingma & Ba 2014). Still, minimising the spiral likelihood function is computationally expensive, requiring several hundred CPU hours. For a more detailed explanation of our method, as well as illustrative examples highlighting its general principles, we refer back to Paper I. The code used in this work is open source and available online1.

4.3. Jackknifing

In order to estimate the statistical uncertainty of our respective data samples, we employed the technique of “delete-d jackknifing” (Efron 1982). With this technique, the statistical uncertainty of some inferred parameter (e.g., ψ) is evaluated by inferring this parameter for a number of subsamples constructed from the full data sample. The statistical variance, written Var(ψ), is proportional to the variance of the inferred parameter between the different subsamples, written Var*(ψ), according to

(23)

where n is the total number of data objects and d is the number of deleted data objects.

Jackknifing should ideally be performed by running inference on all possible subsample constructions. This was not possible for us due to the computational cost of our algorithm. For each data sample and fixed value of Z⊙, we independently constructed ten subsamples by randomly selecting half of the stars (i.e. d = n/2).

5. Baryonic model

The gravitational potential of our model of inference, as stated in Eq. (6), did not rely on any baryonic model. However, after applying our method, we compared our results with a model for the baryonic matter densities, in order to constrain the halo dark matter density and the surface density of a thin dark disk.

We used the baryonic model from Schutz et al. (2018), based on the pre-Gaia studies by Flynn et al. (2006), McKee et al. (2015), and Kramer & Randall (2016a). This baryonic model is also used in other local dynamical mass measurements, such as Sivertsson et al. (2018), Buch et al. (2019), Widmark (2019), and Widmark et al. (2021a). In this model, the total baryonic density is a sum of twelve different components: four gas components (molecular, cold atomic, warm atomic, and hot ionised gas); six stellar components divided by absolute magnitude; white dwarfs; and brown dwarfs. Each of the twelve components is described in terms of its mid-plane matter density (ρt, 0) and its vertical velocity dispersion (σw, t). It is assumed that the respective components are iso-thermal, such that their matter density profiles decay with height according to

(24)

In Fig. 2, we show the matter density distribution and vertical gravitational potential of this baryonic model. In this figure, the baryons are divided into four different groups: stars and dwarfs, cold atomic gas, molecular gas, warm and hot gas. We also include a halo dark matter density of 0.011 ± 0.03 M⊙ pc−3, which encapsulates the approximate range of recent local dark matter density measurements (de Salas & Widmark 2020).

thumbnail Fig. 2.

Baryonic model from Schutz et al. (2018), including a component of halo dark matter, shown in terms of its matter density distribution (left panel) and vertical gravitational potential (right panel). Details are found in Sect. 5.

It is important to note that this model of baryonic matter densities in the solar neighbourhood could potentially suffer from significant systematic errors. For the stellar components, the assumption of iso-thermality are not ideal descriptions when comparing with Gaia data: the stellar components’ vertical velocity distributions are not strictly Gaussian but have heavier tails and the shape of the stellar number density profiles differ from those predicted by the baryonic model (a comparison can be found in the appendix of Widmark et al. 2021a). This is, at least in part, connected to the fact that the stellar components are categorised in terms of absolute magnitude, which does not differentiate them in terms of age and kinematic properties. Perhaps even more worrisome are the matter density distributions of the gas components, which have larger statistical uncertainties and are arguably more prone to significant systematic bias. For example, measuring molecular hydrogen depends on the CO-to-H2 conversion factor, while measuring atomic hydrogen depends on corrections for optical depth (e.g., discussed in Hessman 2015). Furthermore, the cold gas is non-uniformly distributed, both in terms of density and ionisation (Lallement et al. 2003). As such, it is not at all impossible that the baryonic model suffers from systematic errors significantly larger than the reported statistical uncertainties.

6. Results

In this section, we present the results from having applied our method of inference on our respective data samples. Using a summary of those results, we also place constraints on the local dark matter density and the surface density of a thin dark disk. Most figures in this section only show a smaller number of data samples. The corresponding figures for the remaining data samples are found in Appendix A.

In Fig. 3, we show the data histograms of four representative data samples (s = { − 4, −1, 0, 3}). The stellar number density differs quite dramatically between data samples (their number counts are listed in Table 2), mainly depending on the distance. They vary in terms of their dispersion in Z, where the more distant data samples are wider while nearby data samples seem pinched (seen most clearly when comparing data samples s = −4 and s = 0). Furthermore, the more distant data samples have a number density depression close to the Galactic mid-plane (i.e. for low |Z|, seen most clearly for s = −4 and to a lesser extent also for s = 3). All these features are due to selection effects, most significantly because of the completeness of the radial velocity measurements. For the nearby spatial volume, the radial velocity data is deeper in terms of absolute magnitude. The completeness gets progressively worse with greater distances, especially close to the Galactic mid-plane where the radial velocity measurements are obfuscated by stellar crowding.

thumbnail Fig. 3.

Histograms in the (Z, W)-plane, for data samples s = { − 4, −1, 0, 3}. The scale of the colour bar is not linear, but follows sinh[10 × di, j/max(di, j)].

In Fig. 4, we show the phase-space spirals for the same four data samples as in Fig. 3 (s = { − 4, −1, 0, 3}). The visualised quantity is equal to

(25)

thumbnail Fig. 4.

Spirals in the (Z, W)-plane, for data samples s = { − 4, −1, 0, 3}. They are plotted in terms of the relative difference between the data histogram and the fitted bulk density distribution, as defined in Eq. (25). The vertical dashed lines correspond to the boundaries of the θ mask function defined in Eq. (21).

where we have assumed a value of Z⊙ = 10 pc. Furthermore, this quantity is smoothed to an effective bin size of (40 pc) × (2 km s−1) for better visibility. The number density depression close to the Galactic mid-plane is seen clearly in this figure, and is confined to the boundaries of the θ mask function of Eq. (21). There is also some asymmetry with respect to the Galactic mid-plane, at least for the more distant data samples (seen also in the figures in Appendix A). It is evident from both Figs. 3 and 4 that selection effects vary dramatically between the respective data samples. Despite this, the phase-space spirals that emerge when subtracting the bulk density are qualitatively similar, most importantly in terms of their shapes.

In Fig. 5, we show the inferred matter density distribution and inferred gravitational potential for three representative data samples (s = { − 2, 0, 2}). The coloured band corresponds to the 1σ statistical uncertainty given by jackknifing (as explained in Sect. 4.3), including the variance coming from the three different solar height values. The black lines correspond to the mean result for the three different solar heights (Z⊙ = {0, 10, 20} pc). The grey band shows the baryonic model to within a 1σ uncertainty, including a local dark matter density of 0.011 ± 0.003 M⊙ pc−3. The inferred gravitational potential agrees well with the baryonic model, and is very well consistent between samples.

thumbnail Fig. 5.

Inferred matter density distribution (left panels) and inferred gravitational potential (right panels), for data samples s = { − 2, 0, 2}. The black lines correspond to the mean values for different fixed values of Z⊙. We also show 1σ bands of the inferred results and the baryonic density distribution model.

The inferred matter density distribution has a higher relative statistical uncertainty than the inferred gravitational potential, which is expected in any dynamical mass measurement due to it being derived from the gravitational potential’s second order derivative. It also varies more dramatically between samples, and also for the different values of Z⊙. A general trend seen in Fig. 5, as well as for the other data samples shown in Appendix A, is that the inferred matter density distribution is flatter and agrees better with the baryonic model for the data samples closer to the Galactic centre (s ≤ 1). Conversely, the data samples in the direction of the anti-centre (s ≥ 2) are more pinched, and somewhat discrepant with the baryonic model. There is a similar but less pronounced trend for the different values of Z⊙, where especially Z⊙ = 0 pc gives rise to flatter matter density distributions.

In Fig. 6, we show the inferred value of Φ(400 pc) for all twelve data samples. As suggested by the test on one-dimensional simulations in Paper I, the gravitational potential at this height is the quantity that is most robustly inferred. Indeed, the inferred values for Φ(400 pc) are largely consistent between the different data samples. Interestingly, this is the case even for the samples at greater distances (i.e. large |s|), for which the spiral was much less clearly defined (e.g., s = −4 in Fig. 4). This demonstrates that our method works reasonably well even when the data is noisy and suffers from severe selection effects, which are absorbed by the bulk density distribution. It is expected that the mass of the Galactic disk varies with Galactocentric radius, giving rise to a negative slope with respect to the data sample index s. Such a slope is allowed but not preferred given our results, and varies depending on what data samples are considered reliable enough to include in the fit of such a slope. For some data samples, the size of the statistical uncertainties differ quite a lot depending on Z⊙. This is most notable for data sample s = −4, which has a very large uncertainty for Z⊙ = 0 pc, but rather small uncertainties for Z⊙ = {10, 20} pc. The reason is that the inner mask function, m(z, w | ρh) as defined in Eq. (18), is a function of the phase-space coordinates in the rest frame of the Galactic disk. This means that this masked region shifts somewhat in the (Z, W) histogram, depending on the fixed value of Z⊙. Especially for the noisy histogram of data sample s = −4, this small shift has a significant effect on the likelihood function.

thumbnail Fig. 6.

Inferred gravitational potential at height z = 400 pc for all data samples (s from −5 to 5 and ). The markers and vertical lines show the mean and standard deviations as derived from jackknifing, where the colours represent the three fixed values for the height of the Sun (Z⊙ = {0, 10, 20} pc). The grey band represents the baryonic model, including a dark matter density of 0.011 ± 0.003 M⊙ pc−3.

In order to infer the local dark matter density, we compared the inferred value of Φ(400 pc) with that of the baryonic model. The baryons and the matter density correction of Eq. (5) amounted to a total contribution of 265.0 ± 16.0 km2 s−2 to Φ(400 pc). Any surplus with respect to this value was interpreted as halo dark matter. The reason for only using the value of Φ(400 pc) is that this quantity was shown to be the most accurately inferred in tests on one-dimensional simulations (Paper I), and also agrees well between the different data samples and fixed values of Z⊙ used in this work.

When placing constraints on the surface density of a thin dark disk, we used a similar line of reasoning: a thin dark disk was allowed to the extent that there is room for its contribution to Φ(400 pc). We use the same model for the baryonic matter densities, and a dark matter halo density coming from independent studies of the Milky Way circular velocity curve (to use the halo density derived in this or any other work that analyses the vertical motion of stars would be circular reasoning). In order to be conservative when placing thin dark disk constraints, we assumed a local halo dark matter density of ρDM, ⊙ = 0.009 ± 0.003 M⊙ pc−3, which is a somewhat low value with a large uncertainty2. A thin dark disk with a surface density of ΣDD = 1 M⊙ pc−2 and a scale height of hDD = {20,  50,  100} pc, assuming a matter density profile

(26)

contributes with {9.9,  8.9, 7.2} km2 s−2 to Φ(400 pc).

When inferring the local dark matter density, as well as setting constraints on the thin dark disk surface density, we only use the information coming from the nearby data samples, namely |s|≤3. The reason for this is that the more distant data samples have less statistics and suffer from more extreme selection effects. This is seen in for example Fig. 4, where the spiral of data sample s = −4 is poorly defined. Our method does seem to have worked well even for those data samples, but we still found it more cautious to exclude them. When summarising the results of several data samples and a fixed value of Z⊙, the total weighted mean was calculated according to

(27)

with a total statistical variance of

(28)

For the respective values of Z⊙, the summary results of data samples |s|≤3 was

(29)

These three different values of Z⊙ have a total mean value of 301.8 km2 s−2, and the dispersion between them amounts to 4.6 km2 s−2. In order to account for the statistical variance given a fixed value of Z⊙, as well as the variance that arises from changing Z⊙, we added these measurement uncertainties together according to

(30)

Hence, for the joint analysis of samples |s|≤3, we inferred Φ(400 pc) = 301.8 ± 6.0 km2 s−2 and ρDM, ⊙ = 0.0085 ± 0.0039 M⊙ pc−3. The dominant component of this uncertainty comes from the baryonic model (0.0037 M⊙ pc−3), which was added in quadrature to the statistical measurement uncertainty (0.0014 M⊙ pc−3).

Using the assumptions discussed above and a dark disk scale height of 50 pc, the likelihood of its surface density is a Gaussian with a mean and standard deviation of −0.24 ± 2.40 M⊙ pc−2, which is consistent with zero. Negative values for this quantity are unphysical and therefore excluded; the upper limits are given by the likelihood ratio relative to the null hypothesis (i.e. ΣDD = 0 M⊙ pc−2). We obtain an upper 68% (95%) confidence limit of 2.17 M⊙ pc−2 (4.56 M⊙ pc−2). For scale heights of 20 pc and 100 pc, the limits are 1.96 M⊙ pc−2 (4.13 M⊙ pc−2) and 2.71 M⊙ pc−2 (5.59 M⊙ pc−2), respectively. Even if we would have been extremely conservative and assumed no halo dark matter, we would still have placed the most stringent constraints on a thin dark disk surface density, corresponding to an upper 68% (95%) confidence limit of 6.02 M⊙ pc−2 (7.93 M⊙ pc−2), assuming a scale height of 50 pc.

As a comparison, for the data sample with a spatial cut (), the inferred values were

(31)

Using the same analysis as above, applied to this one data sample, we obtained Φ(400 pc) = 312.9 ± 11.2 and ρDM, ⊙ = 0.0111 ± 0.0045 M⊙ pc−3. This result is statistically consistent with the one presented above, although slightly higher in both value and uncertainty. For a thin dark disk with a scale height of 50 pc, the upper 68% (95%) confidence limit is 3.62 M⊙ pc−2 (6.25 M⊙ pc−2).

7. Discussion

Using a new method that extracts information from the time-varying structure of the Milky Way phase-space spiral, we have been able to infer the vertical gravitational potential of the solar neighbourhood for eleven statistically independent data samples. Using these results and a model for the total matter density of baryons, we have inferred the local halo dark matter density, as well as placed the most stringent constraints on the surface density of a thin dark disk.

Overall, the inferred gravitational potential of the respective data samples agree well, at least in terms of Φ(400 pc), which was deemed to be the most robustly inferred quantity in tests on one-dimensional simulations in Paper I. The results are consistent, despite the fact that the data samples are subject to very different selection effects and differ significantly in their number density profiles as a function of height. For data sample s = 0, we see a high stellar number density, especially close to the Galactic mid-plane; conversely, for the distant data samples, most notably |s|≥4, the mid-plane number density is low. Despite this, the different data samples gave rise to very similar results; this illustrates that our method is robust with respect to such systematics, and that selection effects are absorbed by the bulk density distribution.

The inferred matter density distributions agree fairly well with the baryonic model, especially for Z⊙ = 0 pc and s ≤ 1. This can be compared with Widmark (2019) and Widmark et al. (2021a), where they used a similarly free model for the total density, but weighed the Galactic disk using a distribution function fitting method based on the assumption of a steady state. For several statistically independent stellar samples, they inferred a matter density distribution that was very pinched (i.e. with a high mid-plane value and quickly decreasing with height) and argued that their result must be biased by time-varying dynamics. The reason that the results of this work agree better with the baryonic model is probably explained by the fact that the different methods are subject to different systematic biases. For example, the method used in this work is not as sensitive to the distribution of stars with low vertical energies (Ez ≲ Φ(300 pc)).

Using summary statistics of seven data samples (fulfilling |s|≤3), we infer a local halo dark matter density of ρDM, ⊙ = 0.0085 ± 0.0039 M⊙ pc−3. This is a somewhat low value compared to most other analyses of the vertical motion of stars in the solar neighbourhood (although strictly speaking not discrepant as our uncertainty is fairly large), and more in line with more global estimates, such as those coming from the Milky Way rotational velocity curve (de Salas & Widmark 2020). We also set an upper 68% (95%) confidence limit of 2.17 M⊙ pc−2 (4.56 M⊙ pc−2) to the matter density of a thin dark disk with a scale height ≤50 pc. This is significantly stronger than any previous constraints on a thin dark disk (Schutz et al. 2018 set an upper 95% confidence limit of roughly 10 M⊙ pc−2, for a scale height of 50 pc). The statistical uncertainties on the inferred halo dark matter density and thin dark disk surface density are dominated by the uncertainty associated with the baryonic model. Therefore, in order to make further progress, it is crucial to revisit and improve on the baryonic matter density distributions of the solar neighbourhood. As discussed in Sect. 5, it is not implausible that the baryonic model used in this work suffers from significant systematic errors, potentially larger than the reported statistical uncertainties.

For the data sample that only had a cut in spatial coordinates, according to , the results agree with the main analysis to within 1σ, both in terms of Φ(400 pc) and ρDM, ⊙ (where only the latter includes uncertainties from the baryonic model). The results agree particularly well for Z⊙ = 0 pc. This data sample is not statistically independent from the main analysis and was only included as a test of consistency. We consider this result less reliable, due to not having any cuts in angular momentum and thus including stars with strongly elliptical orbits.

In this work, we did not attempt to infer the height of the Sun with respect to the Galactic mid-plane (Z⊙), mainly because of strong selection effects. Rather, we accounted for the uncertainty of this parameter by producing results for three different fixed values of Z⊙ = {0, 10, 20} pc. Interestingly, studies based on the more local spatial volume, within a few hundred parsec, tend to prefer lower values (0–10 pc, e.g., Buch et al. 2019; Widmark et al. 2021a; Gaia Collaboration 2021b); conversely, higher values typically come from studies that reach several kilo-parsec in height (e.g., Jurić et al. 2008; Yao et al. 2017; Bennett & Bovy 2019). This discrepancy indicates that the Galaxy is not perfectly mirror symmetric and that an estimate of Z⊙ can depend on how this quantity is defined. With this caveat in mind, a comparison between the baryonic model and the inferred matter density distributions seems to suggest that our results are vaguely in favour of Z⊙ = 0 pc, which is in agreement with other local studies.

The spatial reach of the method employed in this work is mainly limited by the completeness of the radial velocity sample. This will improve with Gaia’s future data releases, most immanently with the full third data release (EDR3 only contains a cleaned version of the second data release’s radial velocity sample). Apart from such improvements to the data, it should be possible to use the astrometric information of stars even when the radial velocity is missing; especially for more distant parts of the Galactic disk, stars close to the Galactic mid-plane have a vertical velocity that is well approximated by their latitudinal proper motion. With careful data treatment, we might be able to circumvent the issue of missing radial velocity measurements and weigh different parts of the Galactic disk with high precision.

In order to perform more careful tests and further refine this new method, we plan to apply it to high resolution three-dimensional galaxy simulations. To reproduce well-resolved phase-space spirals requires simulations with roughly a billion particles, which only very recently have become feasible (Asano et al. 2020, Hunt et al., in prep.).

8. Conclusion

For the first time, we have employed a method for weighing the Galactic disk using the time-varying structure of the Milky Way phase-space spiral. Our method extracts information from the shape of the phase-space spiral and is complementary to traditional methods that are based on the assumption of a steady state. Using a baryonic model, we have inferred a local halo dark matter density of ρDM, ⊙ = 0.0085 ± 0.0039 M⊙ pc−3 = 0.32 ± 0.15 GeV cm−3, which is consistent with other recent measurements. Using conservative assumptions, we have been able to place the most stringent constraints on the surface density of a thin dark disk: a 95% confidence limit of roughly 5 M⊙ pc−2, assuming a dark disk scale height ≤50 pc.

For both the halo dark matter density and the surface density of a thin dark disk, the statistical uncertainty is dominated by the baryonic model. As discussed in Sect. 5, this model is somewhat outdated and could potentially suffer from significant systematic errors, both in terms of its stellar and gaseous components. We plan to update this model using new data in the near future.

Our method places strong constraints on the weight of the Galactic disk and its dark sector components. In terms of its precision, it is highly competitive with respect to methods based on the assumption of a steady state. In a general sense, this illustrates that time-varying structures, that break the assumption of a steady state, are not solely obstacles to dynamical mass measurements, but can in fact be assets containing useful information.

2

In a recent review on the local dark matter density by de Salas & Widmark (2020), the measurements based on the Milky Way circular velocity curve lie in the range 0.008–0.013 M⊙ pc−3, which is lower than the range of 0.011–0.016 M⊙ pc−3 for more local analyses. For example, de Salas et al. (2019) report 0.008 ± 0.0008 M⊙ pc−3 or 0.001 ± 0.0001 M⊙ pc−3 depending on the baryonic model, Cautun et al. (2020) report 0.009 ± 0.0005 M⊙ pc−3, and Petač (2020) reports 0.0096 ± 0.0005 M⊙ pc−3 or 0.0010 ± 0.0005 M⊙ pc−3 depending on the shape of the dark matter halo. With this in mind, ρDM, ⊙ = 0.009 ± 0.003 M⊙ pc−3 is a conservative choice when constraining the thin dark disk surface density, because the assumed ρDM, ⊙ has a low mean and a large uncertainty.

Acknowledgments

We thank the anonymous referee for a constructive review. AW acknowledges support from the Carlsberg Foundation via a Semper Ardens grant (CF15-0384). CL acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 852839). PFdS acknowledges support by the Vetenskapsrådet (Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centre for Cosmoparticle Physics. GM acknowledges funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 834148). This work made use of an HPC facility funded by a grant from VILLUM FONDEN (projectnumber 16599). This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research utilised the following open-source Python packages: MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (McKinney 2010), TENSORFLOW (Abadi et al. 2015).

References

  1. Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, software available from tensorflow.org [Google Scholar]
  2. Ablimit, I., Zhao, G., Flynn, C., & Bird, S. A. 2020, ApJ, 895, L12 [Google Scholar]
  3. Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42 [NASA ADS] [CrossRef] [Google Scholar]
  4. Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360 [NASA ADS] [CrossRef] [Google Scholar]
  5. Asano, T., Fujii, M. S., Baba, J., et al. 2020, MNRAS, 499, 2416 [CrossRef] [Google Scholar]
  6. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bergemann, M., Sesar, B., Cohen, J. G., et al. 2018, Nature, 555, 334 [NASA ADS] [CrossRef] [Google Scholar]
  8. Binney, J., & Schönrich, R. 2018, MNRAS, 481, 1501 [NASA ADS] [CrossRef] [Google Scholar]
  9. Binney, J., & Tremaine, S. 2008, Galactic Dynamics (Second Edition (Princeton University Press)) [Google Scholar]
  10. Bland-Hawthorn, J., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 486, 1167 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bonaca, A., & Hogg, D. W. 2018, ApJ, 867, 101 [NASA ADS] [CrossRef] [Google Scholar]
  12. Buch, J., Leung, J. S. C., & Fan, J. 2019, JCAP, 2019, 026 [CrossRef] [Google Scholar]
  13. Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513 [NASA ADS] [CrossRef] [Google Scholar]
  14. Buschmann, M., Safdi, B. R., & Schutz, K. 2021, ArXiv e-prints [arXiv:2103.05000] [Google Scholar]
  15. Caputo, A., Zavala, J., & Blas, D. 2018, Phys. Dark Univ., 19, 1 [CrossRef] [Google Scholar]
  16. Carlin, J. L., DeLaunay, J., Newberg, H. J., et al. 2013, ApJ, 777, L5 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291 [Google Scholar]
  18. Chakrabarti, S., Chang, P., Lam, M. T., Vigeland, S. J., & Quillen, A. C. 2021, ApJ, 907, L26 [CrossRef] [Google Scholar]
  19. Chequers, M. H., & Widrow, L. M. 2017, MNRAS, 472, 2751 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chequers, M. H., Widrow, L. M., & Darling, K. 2018, MNRAS, 480, 4244 [CrossRef] [Google Scholar]
  21. Cole, D. R., & Binney, J. 2017, MNRAS, 465, 798 [NASA ADS] [CrossRef] [Google Scholar]
  22. Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429 [NASA ADS] [CrossRef] [Google Scholar]
  23. Deng, L.-C., Newberg, H. J., Liu, C., et al. 2012, Res. Astron. Astrophys., 12, 735 [NASA ADS] [CrossRef] [Google Scholar]
  24. de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, JCAP, 2019, 037 [Google Scholar]
  25. de Salas, P. F., & Widmark, A. 2020, Rep. Prog. Phys., submitted, [arXiv:2012.11477] [Google Scholar]
  26. Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Nature, 454, 735 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  27. Efron, B. 1982, The Jackknife, the Bootstrap and Other Resampling Plans (Society for industrial and applied mathematics) [Google Scholar]
  28. Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120 [NASA ADS] [CrossRef] [Google Scholar]
  29. Facchinetti, G., Lavalle, J., & Stref, M. 2020, ArXiv e-prints [arXiv:2007.10392] [Google Scholar]
  30. Fan, J., Katz, A., Randall, L., & Reece, M. 2013a, Phys. Rev. Lett., 110, 211302 [Google Scholar]
  31. Fan, J., Katz, A., Randall, L., & Reece, M. 2013b, Phys. Dark Univ., 2, 139 [CrossRef] [Google Scholar]
  32. Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiss, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  33. Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487 [Google Scholar]
  34. Gaia Collaboration (Brown, A. G. A., et al.) 2021a, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Gaia Collaboration (Smart, R. L., et al.) 2021b, A&A, 649, A6 [EDP Sciences] [Google Scholar]
  36. Gilmore, G., Randich, S., Asplund, M., et al. 2012, Messenger, 147, 25 [Google Scholar]
  37. Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [Google Scholar]
  38. Gómez, F. A., White, S. D. M., Marinacci, F., et al. 2016, MNRAS, 456, 2779 [Google Scholar]
  39. Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Guo, R., Liu, C., Mao, S., et al. 2020, MNRAS, 495, 4828 [Google Scholar]
  41. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [CrossRef] [PubMed] [Google Scholar]
  42. Hattori, K., Valluri, M., & Vasiliev, E. 2020, ArXiv e-prints [arXiv:2012.03908] [Google Scholar]
  43. Hessman, F. V. 2015, A&A, 579, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  45. Johnston, K. V., Law, D. R., & Majewski, S. R. 2005, ApJ, 619, 800 [Google Scholar]
  46. Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  47. Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2014, ApJ, 794, 59 [NASA ADS] [CrossRef] [Google Scholar]
  48. Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Kazantzidis, S., Bullock, J. S., Zentner, A. R., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254 [Google Scholar]
  50. Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Kingma, D. P., & Ba, J. 2014, ArXiv e-prints [arXiv:1412.6980] [Google Scholar]
  52. Klasen, M., Pohl, M., & Sigl, G. 2015, Prog. Part. Nucl. Phys., 85, 1 [Google Scholar]
  53. Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597 [Google Scholar]
  54. Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260 [Google Scholar]
  55. Kramer, E. D., & Randall, L. 2016a, ApJ, 824, 116 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kramer, E. D., & Randall, L. 2016b, ApJ, 829, 126 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75 [Google Scholar]
  58. Lallement, R., Welsh, B. Y., Vergely, J. L., Crifo, F., & Sfeir, D. 2003, A&A, 411, 447 [EDP Sciences] [Google Scholar]
  59. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286 [Google Scholar]
  60. Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134 [Google Scholar]
  61. Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586 [CrossRef] [Google Scholar]
  62. Li, Z.-Z., Qian, Y.-Z., Han, J., et al. 2020, ApJ, 894, 10 [Google Scholar]
  63. Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995 [NASA ADS] [CrossRef] [Google Scholar]
  64. McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13 [NASA ADS] [CrossRef] [Google Scholar]
  65. McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56 [Google Scholar]
  66. McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
  67. Nesti, F., & Salucci, P. 2013, JCAP, 2013, 016 [Google Scholar]
  68. Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245 [Google Scholar]
  69. Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001 [CrossRef] [Google Scholar]
  70. Nobile, E. D. 2021, ArXiv e-prints [arXiv:2104.12785] [Google Scholar]
  71. Petač, M. 2020, Phys. Rev. D, 102, 123028 [Google Scholar]
  72. Price-Whelan, A. M., Johnston, K. V., Sheffield, A. A., Laporte, C. F. P., & Sesar, B. 2015, MNRAS, 452, 676 [NASA ADS] [CrossRef] [Google Scholar]
  73. Purcell, C. W., Bullock, J. S., & Kaplinghat, M. 2009, ApJ, 703, 2275 [NASA ADS] [CrossRef] [Google Scholar]
  74. Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 063101 [NASA ADS] [CrossRef] [Google Scholar]
  75. Read, J. I., Lake, G., Agertz, O., & Debattista, V. P. 2008, MNRAS, 389, 1041 [NASA ADS] [CrossRef] [Google Scholar]
  76. Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [Google Scholar]
  77. Ruchti, G. R., Read, J. I., Feltzing, S., Pipino, A., & Bensby, T. 2014, MNRAS, 444, 515 [Google Scholar]
  78. Salomon, J.-B., Bienaymé, O., Reylé, C., Robin, A. C., & Famaey, B. 2020, A&A, 643, A75 [EDP Sciences] [Google Scholar]
  79. Sanders, J. L., & Das, P. 2018, MNRAS, 481, 4093 [NASA ADS] [CrossRef] [Google Scholar]
  80. Schönrich, R., & Dehnen, W. 2018, MNRAS, 478, 3809 [Google Scholar]
  81. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
  82. Schutz, K., Lin, T., Safdi, B. R., & Wu, C.-L. 2018, Phys. Rev. Lett., 121, 081101 [NASA ADS] [CrossRef] [Google Scholar]
  83. Sheffield, A. A., Price-Whelan, A. M., Tzanidakis, A., et al. 2018, ApJ, 854, 47 [Google Scholar]
  84. Sivertsson, S., Silverwood, H., Read, J. I., Bertone, G., & Steger, P. 2018, MNRAS, 478, 1677 [Google Scholar]
  85. Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685 [NASA ADS] [CrossRef] [Google Scholar]
  86. Stoehr, F., White, S. D. M., Springel, V., Tormen, G., & Yoshida, N. 2003, MNRAS, 345, 1313 [NASA ADS] [CrossRef] [Google Scholar]
  87. Stref, M., & Lavalle, J. 2017, Phys. Rev. D, 95, 063003a [NASA ADS] [CrossRef] [Google Scholar]
  88. Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279 [NASA ADS] [CrossRef] [Google Scholar]
  89. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  90. Vogelsberger, M., Helmi, A., Springel, V., et al. 2009, MNRAS, 395, 797 [CrossRef] [Google Scholar]
  91. Widmark, A. 2019, A&A, 623, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Widmark, A., & Monari, G. 2018, MNRAS, 482, 262 [Google Scholar]
  93. Widmark, A., Malhan, K., de Salas, P. F., & Sivertsson, S. 2020, MNRAS, 496, 3112 [Google Scholar]
  94. Widmark, A., de Salas, P. F., & Monari, G. 2021a, A&A, 646, A67 [CrossRef] [EDP Sciences] [Google Scholar]
  95. Widmark, A., Laporte, C., & de Salas, P. F. 2021b, A&A, 650, A124 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239 [Google Scholar]
  97. Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41 [Google Scholar]
  98. Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971 [Google Scholar]
  99. Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [NASA ADS] [CrossRef] [Google Scholar]
  100. Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105 [NASA ADS] [CrossRef] [Google Scholar]
  101. Xu, Y., Liu, C., Tian, H., et al. 2020, ApJ, 905, 6 [NASA ADS] [CrossRef] [Google Scholar]
  102. Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377 [Google Scholar]
  103. Yao, J. M., Manchester, R. N., & Wang, N. 2017, MNRAS, 468, 3289 [Google Scholar]

Appendix A: Plots of the remaining data samples

In this appendix, we present the plots corresponding to Figs. 3, 4, and 5, for the data samples that were not already shown in Sect. 6. We show the data histograms in Fig. A.1, the extracted spirals in Fig. A.2, and the inferred matter density distribution and gravitational potential in Fig. A.3.

thumbnail Fig. A.1.

Same as Fig. 3, but for data samples s = { − 5, −3, −2, 1, 2, 4, 5} and .

thumbnail Fig. A.1.

continued.

thumbnail Fig. A.2.

Same as Fig. 4, but for s = { − 5, −3, −2, 1, 2, 4, 5} and .

thumbnail Fig. A.2.

continued.

thumbnail Fig. A.3.

Same as Fig. 5, but for s = { − 5, −4, −3, −1, 1, 3, 4, 5} and .

thumbnail Fig. A.3.

continued.

thumbnail Fig. A.3.

continued.

All Tables

Table 1.

Stellar number counts of radial velocity information for our eleven main data samples.

In the text
Table 2.

Total number count of our main eleven data samples summed together, with the additional constraints that |Z|< 800 pc.

In the text
Table 3.

Free parameters in our model of inference.

In the text

All Figures

thumbnail Fig. 1.

Histogram of angular momentum (Lz) relative to the angular momentum of a circular orbit (Lz, circ.), for different cuts in Galactocentric radius, where R − R⊙ ∈ [100s − 50, 100s + 50] pc. The grey band highlights the region where Lz is within 10% of Lz, circ..
In the text
thumbnail Fig. 2.

Baryonic model from Schutz et al. (2018), including a component of halo dark matter, shown in terms of its matter density distribution (left panel) and vertical gravitational potential (right panel). Details are found in Sect. 5.
In the text
thumbnail Fig. 3.

Histograms in the (Z, W)-plane, for data samples s = { − 4, −1, 0, 3}. The scale of the colour bar is not linear, but follows sinh[10 × di, j/max(di, j)].
In the text
thumbnail Fig. 4.

Spirals in the (Z, W)-plane, for data samples s = { − 4, −1, 0, 3}. They are plotted in terms of the relative difference between the data histogram and the fitted bulk density distribution, as defined in Eq. (25). The vertical dashed lines correspond to the boundaries of the θ mask function defined in Eq. (21).
In the text
thumbnail Fig. 5.

Inferred matter density distribution (left panels) and inferred gravitational potential (right panels), for data samples s = { − 2, 0, 2}. The black lines correspond to the mean values for different fixed values of Z⊙. We also show 1σ bands of the inferred results and the baryonic density distribution model.
In the text
thumbnail Fig. 6.

Inferred gravitational potential at height z = 400 pc for all data samples (s from −5 to 5 and ). The markers and vertical lines show the mean and standard deviations as derived from jackknifing, where the colours represent the three fixed values for the height of the Sun (Z⊙ = {0, 10, 20} pc). The grey band represents the baryonic model, including a dark matter density of 0.011 ± 0.003 M⊙ pc−3.
In the text
thumbnail Fig. A.1.

Same as Fig. 3, but for data samples s = { − 5, −3, −2, 1, 2, 4, 5} and .
In the text
thumbnail Fig. A.1.

continued.
In the text
thumbnail Fig. A.2.

Same as Fig. 4, but for s = { − 5, −3, −2, 1, 2, 4, 5} and .
In the text
thumbnail Fig. A.2.

continued.
In the text
thumbnail Fig. A.3.

Same as Fig. 5, but for s = { − 5, −4, −3, −1, 1, 3, 4, 5} and .
In the text
thumbnail Fig. A.3.

continued.
In the text
thumbnail Fig. A.3.

continued.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.