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/00046361/202141466  
Published online  15 September 2021 
Weighing the Galactic disk using phasespace spirals
II. Most stringent constraints on a thin dark disk using Gaia EDR3
^{1}
Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark
email: axel.widmark@nbi.ku.dk
^{2}
Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEECUB), Martí i Franquès 1, 08028 Barcelona, Spain
^{3}
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, Chiba 2778583, Japan
^{4}
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 10691 Stockholm, Sweden
^{5}
Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l’Université, 67000 Strasbourg, France
Received:
4
June
2021
Accepted:
5
July
2021
Using the method that was developed in the first paper of this series, we measured the vertical gravitational potential of the Galactic disk from the timevarying structure of the phasespace spiral, using data from Gaia as well as supplementary radial velocity information from legacy spectroscopic surveys. For eleven independent data samples, we inferred gravitational potentials that were in good agreement, despite the data samples’ varied and substantial selection effects. Using a model for the baryonic matter densities, we inferred a local halo dark matter density of 0.0085 ± 0.0039 M_{⊙} pc^{−3} = 0.32 ± 0.15 GeV cm^{−3}. We were also able to place the most stringent constraint on the surface density of a thin dark disk with a scale height ≤50 pc, corresponding to an upper 95% confidence limit of roughly 5 M_{⊙} pc^{−2} (compared to the previous limit of roughly 10 M_{⊙} pc^{−2}, given the same scale height). For the inferred halo dark matter density and thin dark disk surface density, the statistical uncertainties are dominated by the baryonic model, which potentially could also suffer from a significant systematic error. With this level of precision, our method is highly competitive with traditional methods that rely on the assumption of a steady state. In a general sense, this illustrates that timevarying dynamical structures are not solely obstacles to dynamical mass measurements, but they can also be regarded as assets containing useful information.
Key words: Galaxy: kinematics and dynamics / Galaxy: disk / solar neighborhood / astrometry
© 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, coplanar 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 subspecies with strong dissipative selfinteractions (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 timevarying dynamical effects.
Although disk formation requires a fairly quiescent accretion history (Freeman & BlandHawthorn 2002), nonequilibrium 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 largescale 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; PriceWhelan et al. 2015; Sheffield et al. 2018; Bergemann et al. 2018), possibly sourced by selfexcitation (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), flybys (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 timeevolving nature. Using the second data release from Gaia, Antoja et al. (2018) revealed the presence of a phasespace spiral in the solar neighbourhood. The presence of the phasespace 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 diskwide disequilibrium phenomenon. This reinforces the connection between the perturbations seen in the solar neighbourhood and the Galactic outskirts, as was originally predicted by preGaia 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 phasespace 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; BlandHawthorn 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 phasespace spiral is a timevarying 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 phasespace 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 onedimensional 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
where Z_{⊙} is the Sun’s height with respect to the Galactic midplane. 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
We adopt values of V_{⊙} = {11.1, 12.24, 7.25} km s^{−1} (Schönrich et al. 2010).
The Poisson equation is written
where we have neglected the contribution in the azimuthal direction (which is zero in the case of rotational symmetry). The second term of the lefthand side is known as the circular velocity term, because the circular velocity of the Galactic plane v_{c} is given by
This term can be written as a matter density correction, as
Because the rotational velocity curve is close to flat, this correction is small; the slope of the circular velocity is roughly ∂v_{c}/∂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
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
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.
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 phasespace. 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’ sixdimensional phasespace 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 phasespace 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 Gband 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:
3.2. Phasespace cuts
For the main data samples, we made the following phasespace 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
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 midpoint and multiplying it by the circular velocity, according to
assuming a flat rotation curve with a circular velocity of v_{c} = 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 L_{z} distribution is not perfectly centred on L_{z, circ}, but rather slightly shifted towards smaller values, as expected due to asymmetric drift (Binney & Tremaine 2008).
Fig. 1. Histogram of angular momentum (L_{z}) relative to the angular momentum of a circular orbit (L_{z, circ.}), for different cuts in Galactocentric radius, where R − R_{⊙} ∈ [100s − 50, 100s + 50] pc. The grey band highlights the region where L_{z} is within 10% of L_{z, 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 kiloparsec.
In summary, the phasespace cuts of the main data samples can be written:
where s is an integer in range −5 to 5.
To validate our choice of phasespace 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.
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 twodimensional 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 d_{i, j}, where the indices (i, j) label the respective bins. Similar to the phasespace 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 singlearmed 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 socalled tilt term accounts for the coupling between the radial and vertical directions, and is a derivative of the bulk phasespace density distribution with respect to the radial direction. Because our method disregards the bulk density, no such term enters our method. The radialvertical coupling of the gravitational potential does induce a radial motion, but for the relevant range of vertical energies (E_{z} < Φ(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 phasespace 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 phasespace densities
In our method, we fit a phasespace distribution to data, consisting of a product of bulk and a spiral phasespace 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 phasespace density parameters, written Ψ_{bulk}, and the spiral phasespace 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.
Free parameters in our model of inference.
We model the bulk density as a Gaussians mixture model according to
where Ψ_{bulk} are the bulk density components, which includes the Gaussian weights (a_{k}) and dispersions (σ_{z, k}, σ_{w, k}).
A star’s vertical oscillation has a total time period of
where E_{z} = Φ(z  ρ_{h})+w^{2}/2 is a star’s vertical energy per mass, z_{max} 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
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 spiralresembling 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 E_{z} evolves according to
where is the initial angle of the perturbation.
The total phasespace density of our analytical model is equal to
In this expression,
is the relative number density of the spiral with respect to the bulk, where α is a unitless 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
where
This sets a lower boundary in E_{z} to the spiral of our analytic model; close to the origin of the (z, w)plane, the spiral is washed out due to selfgravity 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 E_{z} ≳ Φ(300 pc) instead of E_{z} ≳ Φ(400 pc). The reason for this change is that the phasespace spiral of the actual Milky Way is more clearly defined in the inner region than for the onedimensional simulations of Paper I, possibly because the effects of selfgravity are not as strong and cohesive for the more complex, threedimensional kinematics of our Galaxy.
4.2. Data likelihood and masks
In our method, the phasespace 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 phasespace 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
where and M(Z_{i}, W_{j}) are mask functions described below, and f(Z_{i} + Z_{⊙}, W_{j} + W_{⊙}  Ψ) is the model phasespace density as defined in Eq. (16).
The first mask function in Eq. (20) is a Heaviside step function equal to
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
where sigm is the sigmoid function of Eq. (19). The numerical values of Z_{lim.} and W_{lim.} 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 Z_{lim.} = 800 pc and W_{lim.} = 44 km s^{−1}. In the second step, when the spiral was fitted, we set Z_{lim.} = 700 pc and W_{lim.} = 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 midplane 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 online^{1}.
4.3. Jackknifing
In order to estimate the statistical uncertainty of our respective data samples, we employed the technique of “deleted 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
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 preGaia 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 midplane matter density (ρ_{t, 0}) and its vertical velocity dispersion (σ_{w, t}). It is assumed that the respective components are isothermal, such that their matter density profiles decay with height according to
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).
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 isothermality 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 COtoH_{2} conversion factor, while measuring atomic hydrogen depends on corrections for optical depth (e.g., discussed in Hessman 2015). Furthermore, the cold gas is nonuniformly 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 midplane (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 midplane where the radial velocity measurements are obfuscated by stellar crowding.
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 × d_{i, j}/max(d_{i, j})]. 
In Fig. 4, we show the phasespace spirals for the same four data samples as in Fig. 3 (s = { − 4, −1, 0, 3}). The visualised quantity is equal to
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 midplane 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 midplane, 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 phasespace 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.
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 anticentre (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 onedimensional 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 phasespace 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.
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 km^{2} 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 onedimensional 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 uncertainty^{2}. A thin dark disk with a surface density of Σ_{DD} = 1 M_{⊙} pc^{−2} and a scale height of h_{DD} = {20, 50, 100} pc, assuming a matter density profile
contributes with {9.9, 8.9, 7.2} km^{2} 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
with a total statistical variance of
For the respective values of Z_{⊙}, the summary results of data samples s≤3 was
These three different values of Z_{⊙} have a total mean value of 301.8 km^{2} s^{−2}, and the dispersion between them amounts to 4.6 km^{2} 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
Hence, for the joint analysis of samples s≤3, we inferred Φ(400 pc) = 301.8 ± 6.0 km^{2} 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
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 timevarying structure of the Milky Way phasespace 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 onedimensional 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 midplane; conversely, for the distant data samples, most notably s≥4, the midplane 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 midplane value and quickly decreasing with height) and argued that their result must be biased by timevarying 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 (E_{z} ≲ Φ(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 midplane (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 kiloparsec 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 midplane 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 threedimensional galaxy simulations. To reproduce wellresolved phasespace 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 timevarying structure of the Milky Way phasespace spiral. Our method extracts information from the shape of the phasespace 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 timevarying 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.
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 (CF150384). 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. 63820138993 and the Oskar Klein Centre for Cosmoparticle Physics. GM acknowledges funding from the Agence Nationale de la Recherche (ANR project ANR18CE310006 and ANR19CE310017) 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 opensource Python packages: MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (McKinney 2010), TENSORFLOW (Abadi et al. 2015).
References
 Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: LargeScale Machine Learning on Heterogeneous Systems, software available from tensorflow.org [Google Scholar]
 Ablimit, I., Zhao, G., Flynn, C., & Bird, S. A. 2020, ApJ, 895, L12 [Google Scholar]
 Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Antoja, T., Helmi, A., RomeroGómez, M., et al. 2018, Nature, 561, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Asano, T., Fujii, M. S., Baba, J., et al. 2020, MNRAS, 499, 2416 [CrossRef] [Google Scholar]
 Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 Bergemann, M., Sesar, B., Cohen, J. G., et al. 2018, Nature, 555, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Schönrich, R. 2018, MNRAS, 481, 1501 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics (Second Edition (Princeton University Press)) [Google Scholar]
 BlandHawthorn, J., Sharma, S., TepperGarcia, T., et al. 2019, MNRAS, 486, 1167 [NASA ADS] [CrossRef] [Google Scholar]
 Bonaca, A., & Hogg, D. W. 2018, ApJ, 867, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Buch, J., Leung, J. S. C., & Fan, J. 2019, JCAP, 2019, 026 [CrossRef] [Google Scholar]
 Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513 [NASA ADS] [CrossRef] [Google Scholar]
 Buschmann, M., Safdi, B. R., & Schutz, K. 2021, ArXiv eprints [arXiv:2103.05000] [Google Scholar]
 Caputo, A., Zavala, J., & Blas, D. 2018, Phys. Dark Univ., 19, 1 [CrossRef] [Google Scholar]
 Carlin, J. L., DeLaunay, J., Newberg, H. J., et al. 2013, ApJ, 777, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Cautun, M., BenítezLlambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291 [Google Scholar]
 Chakrabarti, S., Chang, P., Lam, M. T., Vigeland, S. J., & Quillen, A. C. 2021, ApJ, 907, L26 [CrossRef] [Google Scholar]
 Chequers, M. H., & Widrow, L. M. 2017, MNRAS, 472, 2751 [NASA ADS] [CrossRef] [Google Scholar]
 Chequers, M. H., Widrow, L. M., & Darling, K. 2018, MNRAS, 480, 4244 [CrossRef] [Google Scholar]
 Cole, D. R., & Binney, J. 2017, MNRAS, 465, 798 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Deng, L.C., Newberg, H. J., Liu, C., et al. 2012, Res. Astron. Astrophys., 12, 735 [NASA ADS] [CrossRef] [Google Scholar]
 de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, JCAP, 2019, 037 [Google Scholar]
 de Salas, P. F., & Widmark, A. 2020, Rep. Prog. Phys., submitted, [arXiv:2012.11477] [Google Scholar]
 Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Nature, 454, 735 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Efron, B. 1982, The Jackknife, the Bootstrap and Other Resampling Plans (Society for industrial and applied mathematics) [Google Scholar]
 Eilers, A.C., Hogg, D. W., Rix, H.W., & Ness, M. K. 2019, ApJ, 871, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Facchinetti, G., Lavalle, J., & Stref, M. 2020, ArXiv eprints [arXiv:2007.10392] [Google Scholar]
 Fan, J., Katz, A., Randall, L., & Reece, M. 2013a, Phys. Rev. Lett., 110, 211302 [Google Scholar]
 Fan, J., Katz, A., Randall, L., & Reece, M. 2013b, Phys. Dark Univ., 2, 139 [CrossRef] [Google Scholar]
 Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiss, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Freeman, K., & BlandHawthorn, J. 2002, ARA&A, 40, 487 [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021a, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Smart, R. L., et al.) 2021b, A&A, 649, A6 [EDP Sciences] [Google Scholar]
 Gilmore, G., Randich, S., Asplund, M., et al. 2012, Messenger, 147, 25 [Google Scholar]
 Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [Google Scholar]
 Gómez, F. A., White, S. D. M., Marinacci, F., et al. 2016, MNRAS, 456, 2779 [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guo, R., Liu, C., Mao, S., et al. 2020, MNRAS, 495, 4828 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [CrossRef] [PubMed] [Google Scholar]
 Hattori, K., Valluri, M., & Vasiliev, E. 2020, ArXiv eprints [arXiv:2012.03908] [Google Scholar]
 Hessman, F. V. 2015, A&A, 579, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Johnston, K. V., Law, D. R., & Majewski, S. R. 2005, ApJ, 619, 800 [Google Scholar]
 Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kafle, P. R., Sharma, S., Lewis, G. F., & BlandHawthorn, J. 2014, ApJ, 794, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kazantzidis, S., Bullock, J. S., Zentner, A. R., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254 [Google Scholar]
 Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kingma, D. P., & Ba, J. 2014, ArXiv eprints [arXiv:1412.6980] [Google Scholar]
 Klasen, M., Pohl, M., & Sigl, G. 2015, Prog. Part. Nucl. Phys., 85, 1 [Google Scholar]
 Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597 [Google Scholar]
 Koposov, S. E., Rix, H.W., & Hogg, D. W. 2010, ApJ, 712, 260 [Google Scholar]
 Kramer, E. D., & Randall, L. 2016a, ApJ, 824, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, E. D., & Randall, L. 2016b, ApJ, 829, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75 [Google Scholar]
 Lallement, R., Welsh, B. Y., Vergely, J. L., Crifo, F., & Sfeir, D. 2003, A&A, 411, 447 [EDP Sciences] [Google Scholar]
 Laporte, C. F. P., Johnston, K. V., Gómez, F. A., GaravitoCamargo, N., & Besla, G. 2018, MNRAS, 481, 286 [Google Scholar]
 Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134 [Google Scholar]
 Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586 [CrossRef] [Google Scholar]
 Li, Z.Z., Qian, Y.Z., Han, J., et al. 2020, ApJ, 894, 10 [Google Scholar]
 Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56 [Google Scholar]
 McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
 Nesti, F., & Salucci, P. 2013, JCAP, 2013, 016 [Google Scholar]
 Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245 [Google Scholar]
 Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001 [CrossRef] [Google Scholar]
 Nobile, E. D. 2021, ArXiv eprints [arXiv:2104.12785] [Google Scholar]
 Petač, M. 2020, Phys. Rev. D, 102, 123028 [Google Scholar]
 PriceWhelan, A. M., Johnston, K. V., Sheffield, A. A., Laporte, C. F. P., & Sesar, B. 2015, MNRAS, 452, 676 [NASA ADS] [CrossRef] [Google Scholar]
 Purcell, C. W., Bullock, J. S., & Kaplinghat, M. 2009, ApJ, 703, 2275 [NASA ADS] [CrossRef] [Google Scholar]
 Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 063101 [NASA ADS] [CrossRef] [Google Scholar]
 Read, J. I., Lake, G., Agertz, O., & Debattista, V. P. 2008, MNRAS, 389, 1041 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [Google Scholar]
 Ruchti, G. R., Read, J. I., Feltzing, S., Pipino, A., & Bensby, T. 2014, MNRAS, 444, 515 [Google Scholar]
 Salomon, J.B., Bienaymé, O., Reylé, C., Robin, A. C., & Famaey, B. 2020, A&A, 643, A75 [EDP Sciences] [Google Scholar]
 Sanders, J. L., & Das, P. 2018, MNRAS, 481, 4093 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., & Dehnen, W. 2018, MNRAS, 478, 3809 [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
 Schutz, K., Lin, T., Safdi, B. R., & Wu, C.L. 2018, Phys. Rev. Lett., 121, 081101 [NASA ADS] [CrossRef] [Google Scholar]
 Sheffield, A. A., PriceWhelan, A. M., Tzanidakis, A., et al. 2018, ApJ, 854, 47 [Google Scholar]
 Sivertsson, S., Silverwood, H., Read, J. I., Bertone, G., & Steger, P. 2018, MNRAS, 478, 1677 [Google Scholar]
 Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685 [NASA ADS] [CrossRef] [Google Scholar]
 Stoehr, F., White, S. D. M., Springel, V., Tormen, G., & Yoshida, N. 2003, MNRAS, 345, 1313 [NASA ADS] [CrossRef] [Google Scholar]
 Stref, M., & Lavalle, J. 2017, Phys. Rev. D, 95, 063003a [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Vogelsberger, M., Helmi, A., Springel, V., et al. 2009, MNRAS, 395, 797 [CrossRef] [Google Scholar]
 Widmark, A. 2019, A&A, 623, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Widmark, A., & Monari, G. 2018, MNRAS, 482, 262 [Google Scholar]
 Widmark, A., Malhan, K., de Salas, P. F., & Sivertsson, S. 2020, MNRAS, 496, 3112 [Google Scholar]
 Widmark, A., de Salas, P. F., & Monari, G. 2021a, A&A, 646, A67 [CrossRef] [EDP Sciences] [Google Scholar]
 Widmark, A., Laporte, C., & de Salas, P. F. 2021b, A&A, 650, A124 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239 [Google Scholar]
 Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.Y. 2012, ApJ, 750, L41 [Google Scholar]
 Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971 [Google Scholar]
 Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, Y., Liu, C., Tian, H., et al. 2020, ApJ, 905, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377 [Google Scholar]
 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.
Fig. A.1. Same as Fig. 3, but for data samples s = { − 5, −3, −2, 1, 2, 4, 5} and . 
Fig. A.1. continued. 
Fig. A.2. Same as Fig. 4, but for s = { − 5, −3, −2, 1, 2, 4, 5} and . 
Fig. A.2. continued. 
Fig. A.3. Same as Fig. 5, but for s = { − 5, −4, −3, −1, 1, 3, 4, 5} and . 
Fig. A.3. continued. 
Fig. A.3. continued. 
All Tables
Stellar number counts of radial velocity information for our eleven main data samples.
Total number count of our main eleven data samples summed together, with the additional constraints that Z< 800 pc.
All Figures
Fig. 1. Histogram of angular momentum (L_{z}) relative to the angular momentum of a circular orbit (L_{z, circ.}), for different cuts in Galactocentric radius, where R − R_{⊙} ∈ [100s − 50, 100s + 50] pc. The grey band highlights the region where L_{z} is within 10% of L_{z, circ.}. 

In the text 
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 
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 × d_{i, j}/max(d_{i, j})]. 

In the text 
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 
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 
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 
Fig. A.1. Same as Fig. 3, but for data samples s = { − 5, −3, −2, 1, 2, 4, 5} and . 

In the text 
Fig. A.1. continued. 

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

In the text 
Fig. A.2. continued. 

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

In the text 
Fig. A.3. continued. 

In the text 
Fig. A.3. continued. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.