EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 663 (July 2022) A&A, 663 (2022) A15 Full HTML
Free Access
Issue
A&A
Volume 663, July 2022
Article Number A15
Number of page(s) 18
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202142819
Published online 01 July 2022

© ESO 2022

1. Introduction

The dynamics of stars can be related to the gravitational potential that they inhabit via the collisionless Boltzmann equation. For systems in a steady state with certain symmetry properties (typically spherical or axisymmetric) it is possible to find solutions to the phase-space distribution of a stellar tracer population, either through distribution function modelling or via the moments of the Boltzmann equation (Kapteyn 1922; Oort 1932; Bahcall 1984a,b; Kuijken & Gilmore 1989; Crézé et al. 1998; Holmberg & Flynn 2000; Bienayme et al. 2006; Binney & Tremaine 2008; Garbari et al. 2012; Bovy & Rix 2013; Salomon et al. 2020; Guo et al. 2020; Widmark et al. 2021a). Given the relatively quiet conditions necessary to form disk galaxies, the assumption of equilibrium for near-equilibrium systems has been widely and favourably applied to the Milky Way and other galaxies (McMillan 2011; Binney & McMillan 2011). Our place in the Milky Way makes it ideal to accurately measure its gravitational potential and mass distribution, since it is the only system where we have access to full six-dimensional phase space information, from its inner regions all the way to its outermost edge (e.g. Deason et al. 2021). A precise and robust measurement of the gravitational potential is crucial for our general understanding of the Milky Way (Dehnen & Binney 1998; Klypin et al. 2002; Widrow et al. 2008; Weber & de Boer 2010; McMillan 2011, 2017; Kafle et al. 2014; Cole & Binney 2017), and also for probing its distribution of dark matter (Read 2014; Nitschai et al. 2020; Cautun et al. 2020; Li et al. 2020; de Salas & Widmark 2021). The local dark matter density is of fundamental importance for direct and indirect dark matter detection experiments (Jungman et al. 1996; Klasen et al. 2015). In a broader sense, dark matter’s gravitational signatures, studied via stellar dynamics and gravitational lensing, is one of the most competitive avenues for constraining its thus far elusive particle nature (Bertone & Tait 2018; Ferreira 2021). The Gaia satellite has been instrumental to this field, pushing the size of the astrometric sample from a few hundred thousand stars (Perryman et al. 1997) to roughly two billion (Gaia Collaboration 2018a).

With Gaia, it has become evident that the Milky Way disk is not in a steady state; rather, it is a dynamical system with clear time-varying features, for example in the form of radial and vertical asymmetries (Widrow et al. 2012; Williams et al. 2013; Gaia Collaboration 2018b). This is expected from a theoretical perspective, considering that Milky Way-like galaxies do undergo interactions with satellite galaxies which can warp and heat the Galactic disk, induce spiral arms and corrugations (Velazquez & White 1999; Villalobos & Helmi 2008; Purcell et al. 2011; Gómez et al. 2013; Quillen et al. 2020), which are also in interplay with giant molecular clouds (D’Onghia et al. 2013), as well as induce bar formation (Hohl 1971). In fact, such types of dynamical features were seen tentatively also in the pre-Gaia era (e.g. Minchev et al. 2009; now confirmed with Ramos et al. 2018). The broken steady state is perhaps most clearly exemplified by the recently discovered phase-space spiral (Antoja et al. 2018), seen in the plane of vertical position and vertical velocity for stars in the solar neighbourhood. It is probably the remnant of a perturbation event a few hundred million years ago, although the precise time is highly uncertain (see e.g. Darling & Widrow 2019). Pre-Gaia self-consistent models of the Milky Way’s interaction with a Sagittarius-like satellite (Laporte et al. 2018) have been used to interpret the formation of the phase-space spiral as a global phenomenon (Laporte et al. 2019), demonstrating that the spiral is present many kilo-parsec beyond the solar neighbourhood and that it has a similar shape for stars of essentially all ages (at least in the range of roughly 1–9 billion years). Its global nature was further demonstrated in subsequent works (e.g. Xu et al. 2020; Gaia Collaboration 2021a). This indicates that the perturbation should probably have a recent origin, ruling out bar buckling (Khoperskov et al. 2019) which would also violate constraints on the Galactic bar’s chemical structure (Ness et al. 2013; Debattista et al. 2019).

In order to deepen our understanding of the Galaxy, it seems all the more fruitful and necessary to go beyond the common assumptions of symmetry and a steady state (what is referred to as an Ideal Galaxy in de Salas & Widmark 2021). There has been significant effort in testing the traditional steady state based methods against simulations, in order to control for the systematic biases that might arise due to the breaking of symmetry or a steady state. Generally, such methods perform well (Haines et al. 2019; Salomon et al. 2020; Sivertsson et al. 2022), but in principle we can go further than that. As an example, Li & Widrow (2021) apply steady state modelling to the solar neighbourhood but use the phase-space spiral as a consistency check, by comparing a spiral model with the residual that emerges in their inferred equilibrium distribution. There is also the possibility to extract information directly from time-varying dynamical structures, which this work is an example of.

This article is the third part of a longer series about a new method for weighing the Galactic disk, in which the vertical gravitational potential is inferred from the shape of the phase-space spiral in the plane of vertical position and vertical velocity. A way to think about how this method gains its power of inference is to consider the following. Given that the spiral has winded into its current shape in a fairly stable disk, its spiral angle, defined in Eq. (14), is a smooth and monotonic function of vertical energy. In other words, if we trace a curve along the arm of the spiral, from the inside out in the vertical phase-space plane, the vertical energy should be smoothly and strictly increasing. Given this property, the shape of this curve places a very strong constraint on the vertical gravitational potential. This property is expected to hold even though the stars that collectively make up the spiral perturbation have varying dynamical histories and Galactocentric guiding radii. Given that a well-defined single-armed spiral is visible, it is difficult and contrived to imagine a scenario where the property of a strictly increasing vertical energy is broken.

In the two previous papers in this series – Widmark et al. (2021b,c), henceforth referred to as Paper I and Paper II – we have tested our method on one-dimensional simulations and applied our method to the immediate solar neighbourhood using the radial velocity sample of Gaia’s early instalment of its third data release (EDR3). In a subsequent fourth article, we have also tested our method on a high-resolution three-dimensional N-body simulation – Widmark et al. (2022), henceforth referred to as Paper IV – which further motivates and supports this work and its results. Our method produced accurate results for the rich and complex dynamics of an externally perturbed three-dimensional disk galaxy. We also demonstrated our method’s robustness with respect to severe and unknown spatially dependent selection effects, as well as a biased height of the disk mid-plane.

In this third article, we applied our method to distant regions of the Galactic disk, using the Gaia EDR3 proper motion sample. While radial velocity information is essential for seeing the phase-space spiral in the immediate solar neighbourhood, it is less important when observing the spiral at a distance of a few kilo-parsec. The reason is that disk stars at such distances have a small Galactic latitude (b ≃ 0 deg), such that the proper motion in the latitudinal direction has a close projection to vertical velocity (proportional to cos b), while the contribution coming from the radial velocity component is small (proportional to sin b). We constructed 360 separate data samples, by dividing the Galactic disk into different area cells in the directions parallel to the Galactic plane, with a distance bin length of 200 pc in the range 1.4–3.4 kpc and a Galactic longitude bin length of ten degrees. This analysis is the first of its kind, in the sense that we are weighing distant regions of the Galactic disk with a high spatial resolution in the directions parallel to the Galactic plane. Because all data samples are fitted individually, our inference is not subject to otherwise commonly made assumption about the large-scale spatial structure of the Galactic potential (mainly the assumptions of axisymmetry and a disk matter density that decays exponentially with Galactocentric radius). While we constructed 360 data samples, we could only apply our method to roughly half of them, mainly due to severe selection effects in the general direction of the Galactic centre. In order to extract the shape of the spiral in the remaining spatial volume, we modified the method used in this work as compared to Paper I and Paper II by adding a simple extinction model as a function of spatial position. Unlike traditional methods that are based on the assumption of a steady state, our modelling of selection effects is not required to be very precise, but only good enough in order to robustly extract the shape of the phase-space spiral.

This article is structured as follows. In Sect. 2, we define a coordinate system and a few other key quantities. We describe the data sample construction in Sect. 3 and our model of inference in Sect. 4. In Sect. 5, we present our results. In Sects. 6 and 7, we discuss and conclude.

2. Coordinate system and gravitational potential

In this article, we use a Cartesian system of coordinates centred on the Sun’s location and rest frame, whose spatial coordinates X ≡ {X, Y, Z} correspond to the direction of the Galactic centre, the direction of Galactic rotation and the direction of Galactic north, respectively. The time derivatives in the three directions give the velocities V ≡ {U, V, W}. How these phase-space coordinates are related to the Gaia observables is described in Appendix A.

The height, also referred to as vertical position, with respect to the Galactic plane is written as

(1)

where Z⊙ is the height of the Sun relative to the Galactic mid-plane. The velocity in vertical direction in the rest frame of the Galactic disk is written as

(2)

where W⊙ corresponds to the Sun’s velocity relative to the disk’s bulk motion. In the immediate viscinity of the Sun, we have that W⊙, local = 7.25 km s−1 (from e.g. Schönrich et al. 2010). However, because we are studying distant parts of the disk over which the mode of the vertical velocity distribution can vary, we let W⊙ be a free parameter, fitted individually for the respective regions of Galactic disk where we applied our method. Ideally, Z⊙ should have been made a free parameter as well, but this quantity cannot be robustly inferred for distant disk regions due to severe selection effects. Instead, we assume Z⊙ to be given directly by the height of the Sun in the immediate solar neighbourhood; in other words, we approximated the Galactic disk mid-plane to be perfectly flat within the considered distances (< 3.4 kpc). The Sun’s height is typically evaluated to lie in the range 0–20 pc. A broken Galactic plane symmetry is indicated by the fact that lower values for Z⊙ are preferred in more local studies (e.g. Buch et al. 2019; Widmark et al. 2021a; Gaia Collaboration 2021b), as opposed to studies that reach several kilo-parsec in distance (e.g. Juric 2008; Yao et al. 2017; Bennett & Bovy 2019). For this work, we chose a fixed value of Z⊙ = 10 pc.

In our model, the vertical gravitational potential is written as

(3)

which corresponds to a total matter density consisting of a sum of four matter density components with different scale heights of {200, 400, 800, 1600} pc (this is found via the Poisson equation, see Paper II for further details). Using this functional form, the gravitational potential is free to vary in shape, and is flexible enough to emulate models for the total matter density of the solar neighbourhood (e.g. McKee et al. 2015; Schutz et al. 2018).

The vertical energy per mass is given by the vertical gravitation potential and vertical kinetic energy per mass, according to

(4)

However, throughout this paper, we refer to any “energy per mass” as simply “energy”, for shorthand.

Although it does not directly enter our model of inference, we used the Galactocentric cylindrical radius when interpreting our results. This quantity is given by

(5)

where we assumed a solar position of R⊙ = 8178 pc (GRAVITY Collaboration 2019).

3. Data

The main difference between this work and previous articles in this series is that we have used the parallax and proper motion information of Gaia EDR3 without requiring the availability of radial velocities. When extracting the shape of the spiral in the vertical phase-space plane, we need the vertical velocity (W), while other velocity components (U and V) are largely irrelevant. For distant regions of the Galactic disk, the Galactic latitude is close to zero, such that the proper motion in the latitudinal direction has a close projection to vertical velocity. While this approximation would not be feasible in the immediate solar neighbourhood, it is reasonable for disk stars at kilo-parsec distances. This is discussed more carefully in Sects. 3.2 and 3.3 below.

The radial velocity information of Gaia EDR3 was included when available, and was also supplemented with legacy spectroscopic surveys (compiled in Sanders & Das 2018, including 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). Similar to in Paper II, a supplementary radial velocity was used if the Gaia radial velocity information was missing or had an uncertainty larger than 3 km s−1. When there were radial velocity measurements from several supplementary surveys, we used the measurement with the smallest uncertainty. Additionally, we disregarded the supplementary radial velocity information altogether for stars with discrepant measurement values in the supplementary surveys; this was the case if two separate measurement, each with a precision of at least σRV < 5 km s−1, had a statistical tension of greater than 2.5σ between them.

The data used to construct the data samples of this work consisted of parallax and its associated uncertainty (ϖ and σϖ), Galactic longitude and latitude (l and b), proper motions (μb, μl), the astrometric renormalised unit weight error (RUWE), radial velocity and its associated uncertainty (vRV and σRV, only for a subset of stars). We corrected for the parallax zero-point offset in Gaia EDR3 according to the bias function described in Lindegren et al. (2021a). In order to compute this, we also included the data columns phot_g_mean_mag, nu_eff_used_in_astrometry, pseudocolour, ecl_lat, astrometric_params_solved.

3.1. Data cuts

When constructing our data samples, we made cuts in data quality by requiring RUWE < 1.4 and σϖ < 0.05 mas. We included the radial velocity information only for stars with σRV < 3 km s−1. These cuts in data quality set an effective limit in apparent magnitude, due to dimmer sources having worse astrometric precision; the number of included sources drop quickly for Gaia G-band apparent magnitudes in range 16–17. These quality cuts also induce strong spatially dependent selection effects, because of how dust extinction and stellar crowding affect Gaia’s completeness and astrometric precision.

We divided the disk plane into area cells which are 10 degrees wide in Galactic longitude (l) and 200 pc in terms of distance parallel to the Galactic plane (), ranging from 1400 to 3400 pc. This amounts to a total number of 36 × 10 = 360 separate area cells.

We also performed a cut in the longitudinal proper motion divided by parallax (μl/ϖ), which was performed separately for each individual data sample. We excluded any star for which μl/ϖ was more than one standard deviation away from the mode of that data sample’s distribution. This was done in order to exclude stars that are evidently not moving in unison with the bulk of the stellar disk, such as stars with highly eccentric orbits. This phase-space cut is not ideal, in the sense we are cutting in only one out of two velocities parallel to the Galactic plane; regardless, making this crude cut is helpful in isolating the shape of the phase-space spiral.

The area cells in the (X, Y)-plane, and their respective stellar number counts, are visible in Fig. 1. However, due to severe selection effects, we could only analyse about half of those data samples using our phase-space spiral method. This is discussed further in Sect. 5.

thumbnail Fig. 1.

Stellar number counts in our 360 data samples, shown in their respective locations in the (X, Y)-plane. The Galactic centre is located to the right, while the direction of Galactic rotation is upwards.

3.2. Assigning radial velocities

A star’s vertical velocity in the solar rest frame is equal to

(6)

It is the sum of the projected velocities in the latitudinal and line-of-sight directions, which are proportional to sin b and cos b, respectively. For distant regions of the Galactic disk, the latitude (b) is close to zero, such that the proper motion term dominates in the above equation. For this reason, we can approximate the vertical velocity to decent precision even if we only have proper motion information.

For the stars in our data samples that were missing radial velocity information, we assigned these missing values. We did so by estimating the group velocity in the line-of-sight direction, using the sub-sample of stars with available radial velocities. For each respective data sample, we fitted a second degree polynomial of vRV as a function of b, and assigned the missing radial velocity values according to that function. The second degree polynomial was completely free to vary in this fit and not required to be symmetric with respect to for example b = 0 deg; as such, this fit is independent of the assumed value for Z⊙. Because the data samples were fairly small in the X and Y directions, fitting this simple function was enough to capture the most important dependency of vRV. This interpolation is shown in Fig. 2 for a representative data sample, which contained a subset of 3528 stars with radial velocity information. A discussion about possible biases that can arise from this procedure is found below, in the second to last paragraph of Sect. 3.3.

thumbnail Fig. 2.

Radial velocity (vRV) as a function of Galactic latitude (b) for the data sample with l ∈ [90, 100] deg and . The dashed line corresponds to the fitted function and the scatter point correspond to stars with available radial velocity information. The range in b on the horizontal axis corresponds to a range in height spanning roughly Z ∈ [ − 1, 1] kpc at the relevant distance.

3.3. Observational and systematic uncertainties

A useful way to understand how our method can infer the vertical gravitational potential is that the shape of the phase-space spiral in the (z, w)-plane informs us about what stars have similar vertical energies (see e.g. Fig. 1 in Paper I). Hence, the gravitational potential difference between stars can be known, to the extent that the difference in vertical kinetic energy (w2/2) is known as well. If the height or vertical velocity information is biased, this will propagate into a bias in terms of the inferred gravitational potential. The most significant sources of potential bias are discussed below.

In this work, we made the approximation that the Galactic plane is perfectly flat in the studied volume, such that Eq. (1) holds for all data samples. However, the height of the Galactic mid-plane most likely varies with X and Y. In order to ameliorate this systematic uncertainty we would need a deeper and more detailed understanding of the Galactic disk’s spatial structure, which is unfortunately beyond the scope of this work. Due to incompleteness effects and other data systematics, inferring this spatial dependence is difficult, especially given the data quality cuts that were used in this work. As we saw in our analysis of the immediate solar neighbourhood in Paper II, varying Z⊙ between the values {0, 10, 20} pc affected the inferred value of Φ(400 pc) by only a few per cent. In this work, the disk mid-plane is probably confined to vary by a few tens of parsecs within the studied volume; looking at the larger scale warp of the Galactic disk, as studied using Cepheids in for example Skowron et al. (2019) and Chen et al. (2019), vertical variations of the Galactic mid-plane seem to reach scales of 100 pc and above only at distances around 5 kpc from the Sun. To the extent that the disk mid-plane does vary, it will bias our result by distorting the shape of the inferred phase-space spiral. Roughly speaking, a mid-plane bias of 20 pc could lead our method to mistake Φ(400 pc) for Φ(380 pc) or Φ(420 pc). It seems possible that our assumption about a perfectly flat Galactic plane might bias our inferred potential of the order of ten per cent, especially for the more distant data samples and for potential values close to the Galactic mid-plane. This was further tested on a three-dimensional simulation in Paper IV, were we found that the induced bias was contained within the numerical estimate given above.

We have restricted ourselves to parallax uncertainties smaller than 0.05 mas. The parallax bias that we correct for is typically of the order of 0.02 mas, with potential uncontrolled systematic errors of similar magnitude (Lindegren et al. 2021b,a). For the most distant data samples at , the relative statistical and systematic uncertainties in parallax are of the order of ten per cent. For our method of inference, a biased parallax propagates into a systematic bias in both height and vertical velocity, which are affected with roughly the same numerical constant. In terms of the inferred gravitational potential, the systematic errors in height and vertical velocity counteract each other. As such, a relative bias in parallax of ten per cent should translate to a smaller relative bias in the inferred gravitational potential of at most a few per cent.

For stars with missing radial velocities, we assign those values according to the procedure presented in Sect. 3.2. This assignment is of course not perfect, which unavoidably introduces an error in the resulting vertical velocities. To the extent that these vertical velocity errors are distributed in a symmetric and unbiased manner, they only serve to soften the resolution of the observed phase-space spiral. However, the distribution of radial velocities is actually somewhat skewed and our estimate of its mean is possibly slightly biased, which propagates into a bias in the stars’ vertical velocities (W), according to

(7)

This bias in the radial velocity assignment should be sub-dominant with respect to the total dispersion in radial velocity, which is σRV ≃ 30 km s−1 (varies somewhat between data samples). We make the conservative (i.e. large) estimate that bias (vRV)≲10 km s−1. How this affects the shape of the observed phase-space spiral is illustrated in Fig. 3, where we have assumed a distance of and Z⊙ = 0 pc for simplicity. This spiral inhabits a gravitational potential that follows Eq. (3) with parameter values ρh = {0.06, 0.03, 0, 0} M⊙pc−3, which is representative of the actual data samples and observed spirals that are analysed in this work. We show the shape of the spiral as seen when the radial velocity assignment is positively and negatively biased by 10 km s−1. This bias can of course have a more complicated behaviour, for example vary as a function of z, but should in either case be roughly constrained to lie within the dashed lines of in Fig. 3. In the mid-plane, no bias is propagated from the radial velocity assignment, so the three spirals all cross through the same points along the vertical axis defined by z = 0 pc. Along the horizontal axis, defined by w = 0 km s−1, the three spirals cross through roughly the same points as well. The main difference between the three cases is along the diagonals of the (z, w)-plane, where the distance to the origin of the (z, w)-plane is biased by around five per cent. While a lot of information about the gravitational potential can be gathered from looking at where the spiral crosses the horizontal and vertical axes, our method uses the full shape of the phase-space spiral in the (z, w)-plane and can therefore be biased (although probably most severely in terms of the precise shape of the inferred gravitational potential and matter density distribution, rather than e.g. Φ(400 pc)). The relative bias in the inferred gravitational potential, as propagated from a biased assignment of radial velocities, should be well contained within

(8)

thumbnail Fig. 3.

Shape of the phase-space spiral, represented as a one-dimensional line in the (z, w)-plane, when biased by faulty radial velocity assignments. The dashed blue and orange lines show the spiral when radial velocity assignment is biased by ±10 km s−1, assuming that the spiral is observed at a distance of . The spiral is plotted for vertical energies in range Ez ∈ [Φ(250 pc),Φ(800 pc)]. The plus sign marks the origin of the (z, w)-plane.

In summary, there are potential significant biases coming from the approximation of a flat Galactic plane, from parallax measurements, and from the assignment of radial velocities. There is also a trade-off, in the sense that the former two types of bias are likely most severe at greater distances, while the latter type of bias is most severe at smaller distances. All in all, systematic biases in terms of Φ(400 pc) should be contained to a relative error smaller than about ten per cent, where the more severe biases apply to the most nearby and most distant data samples.

4. Model of inference

The method of inference used in this work builds upon Paper I and Paper II; it more closely resembles the latter, which was also on analysis on Gaia EDR3 data. The most significant modification with respect to previous papers is that we include a mask function that depends on z, which models the severe selection effects that are present in distant regions of the Galactic disk. Secondly, the gravitational potential is modelled as a sum of four components with scale heights of {200, 400, 800, 1600} pc (see Eq. (3)), which is twice that of previous papers; due to strong selection effects and a potentially warping disk mid-plane, we did not expect to robustly infer the precise shape of the gravitational potential, especially at low heights, and thus opted for a less flexible functional form. A third important modification is that Z⊙ is a fixed parameter, while W⊙ is free. These parameters were free to vary in Paper I, but fixed parameters in Paper II. The method of inference used in this work is described below, with an emphasis on its differences with respect to the previous articles in this series, especially the extinction mask function.

4.1. Phase-space densities

The full phase-space density of our model of inference is equal to

(9)

It consists of three distributions that are free to vary: a bulk density distribution, B(z, w | Ψbulk); an extinction mask in height, Ξ(z | Ψz-mask); and a relative spiral density perturbation, S(z, w | Ψspiral). These three distributions depend on the model’s free parameters Ψ = {Ψbulk,  Ψz-mask,  Ψspiral}, which are listed in Table 1. The quantity m(z, w) is a fixed inner mask function. All these distributions are defined below.

Table 1.

Free parameters in our model of inference.

The bulk density distribution is equal to

(10)

which is a Gaussian mixture model, with the constraints that all Gaussians are centred on the same point in the (z, w)-plane and have zero correlations between the z and w directions. The respective Gaussians are labelled by an index k = {1, 2, …, K}. In traditional methods that are based on the assumption of a steady state, the bulk density is used to infer the gravitational potential, for example by fitting it to data under the requirement that it fulfils the stationary collisionless Boltzmann equation. In our modelling, the bulk density does not inform, nor is informed by, the gravitational potential; as such, it is solely a background distribution fitted in order to extract the shape of the spiral.

The data samples that are studied in this work suffer from severe selection effects, mainly due to stellar crowding and dust extinction. Especially because of dust, these selection effects are difficult to estimate. While incompleteness is to a significant extent induced by our cuts in parallax uncertainty (see Sect. 3.1), it is not enough to model the incompleteness as a function of apparent magnitude and position on the sky. Dust extinction for disk stars at these distances can be several magnitudes in the Gaia G-band, and in the relevant spatial volumes current three-dimensional dust maps are not very precise (e.g. Fig. 9 in Lallement et al. 2019). Thankfully, due to the nature of this selection, we can assume that it depends mainly on spatial position, rather than velocity, and that it is most severe close to the Galactic mid-plane. Because of this, we modelled the selection effects by using an extinction mask function written as

(11)

This mask function is fitted to data and its free parameters are all written with hats and with an index l = {1, 2, …, L}. The inclusion of this mask function constitutes the most significant modification with respect to the method used in Papers I and II. The purpose of this mask function is not to model selection effects very accurately; in fact, we expect there to be degeneracies between the fitted extinction mask and bulk density distribution. Rather, the purpose of this mask is to facilitate the extraction of the phase-space spiral’s shape in the (z, w)-plane. This aspect of our method has been thoroughly tested in Paper IV, where we subjected the data to selection effects similar to what is seen in this work. In these tests, could accurately extract the spiral and vertical gravitational potential, despite degeneracies between the fitted bulk and extinction mask. This point is further supported by the results of Paper II, where the respective data samples were inconsistent in terms of the scale heights of the fitted bulk distributions (due to the strong distance dependence of the Gaia radial velocity sample) but consistent in terms of their extracted spirals and vertical gravitational potentials.

The spiral relative density perturbation is identical to how it was defined in Paper II. It is written as

(12)

This depends on the angles of a given phase-space point, which is defined as

(13)

as well as that of the spiral,

(14)

In these expressions, zmax is the maximum height that a star reaches, and P is the period of vertical oscillation, given by

(15)

The quantity m(z, w) is a inner mask function. It defines the inner boundary around the origin of the (z, w)-plane, within which the phase-space spiral is not seen. The mask function is equal to

(16)

where Std(w) is the data sample’s standard deviation in vertical velocity, and

(17)

is a sigmoid function. This differs slightly from how our method was formulated in Papers I and II, where the inner mask function was defined in terms of a boundary in vertical energy equal to Φ(300 pc); in those earlier paper, the inner mask’s boundary in velocity varied with the free parameters ρh. We opted to change this in order to make our algorithm somewhat more stable; in this updated version, the region of the (z, w)-plane that the fitted spiral is sensitive to does not change with the free parameters of our model.

4.2. Data likelihood and masks

The fitting procedure is essentially the same is in the previous papers in this series. We fitted the bulk density and mask function (Ψbulk and Ψz-mask) in a first step, without any spiral density perturbation present. In a second step, we fitted the relative phase-space density of the spiral (Ψspiral) while the bulk and mass function remain fixed. The data is reduced to a histogram in the (z, w)-plane, written di, j, where the respective bins are labelled by (i, j) and have widths of 20 pc and roughly 1 km s−1. The logarithm of the likelihood is written as

(18)

where

(19)

is an outer mask function, which defines the circular region in the (z, w)-plane in which we perform our fit. The two steps of our fitting procedure use different outer boundaries: for the first step, when fitting the bulk density, we set zlim. = 7/2 × 300 pc and wlim. = 7/2 × Std(w); for the second step, when fitting the spiral, we set zlim. = 7/3 × 300 pc and wlim. = 7/3 × Std(w). This is similar to the outer bounds that we applied in previous papers in this series, although the boundary in vertical velocity varies between data samples, mainly because Std(w) becomes smaller with greater Galactocentric radius.

In order to make the minimisation algorithm computationally tractable, it is implemented in TENSORFLOW. Readers can refer to Paper I for further details.

5. Results

When running our minimisation algorithm, we used a total number of K = 6 Gaussians for the bulk density distribution and L = 6 Gaussians for the mask function. We found that modelling of the selection effects did not improve when increasing L.

The mode of the vertical velocity distribution, which was inferred in the first step of our minimisation algorithm, is shown in Fig. 4. It is expressed in terms of W⊙, local − W⊙, where W⊙, local = 7.25 km s−1 is the vertical velocity of the Sun with respect to the Galactic disk in the immediate solar neighbourhood, while W⊙ is the corresponding free parameter for the Galactic disk at the position of the respective data samples. There is significant dependence on X and Y with regards to this parameter, with especially low values in the direction of l ≃ 320 deg and high values in the direction of l ≃ 190 deg. These result agree well with those found by Gaia Collaboration (2018b, see their Fig. 10), Poggio et al. (2018, see their Fig. 3), and Martinez-Medina et al. (2022, see their Fig. 6).

thumbnail Fig. 4.

Bulk vertical velocity relative to that of the immediate solar neighbourhood. This is equivalent to W⊙, local − W⊙, where W⊙, local = 7.25 km s−1 and W⊙ is a free parameter fitted in the first step of our minimisation algorithm. The Galactic centre is located to the right, while the direction of Galactic rotation is upwards.

In Fig. 5, we show the data histogram, fitted bulk density, and spiral as seen in the data and best fit, for the data sample with l ∈ [120, 130] deg and . This specific data sample was chosen as a representative and illustrative example of how our algorithm can extract and fit the phase-space spiral despite severe selection effects. Looking at panels a and b, which show the data histogram and the fitted bulk times z-mask, there are clear extinction features in the form of vertical stripes, seen most clearly at the height of z ≃ 200 pc. In panel a, it is very difficult to see the shape of the spiral by eye. However, a spiral clearly emerges in the data after removing the bulk times z-mask in panel c, and its shape is well reproduced by the fitted spiral in panel d. The asymmetry of these extinction features with respect to the Galactic mid-plane illustrate why Z⊙ is not a free parameter in our model, as this quantity would be strongly biases by selection effects.

thumbnail Fig. 5.

Data and fitted phase-space density of the data sample with l ∈ [120, 130] deg and . All four panels span the same range of the (z, w)-plane and show: (a) the data histogram; (b) the fitted bulk and z-mask; (c) the spiral as seen in the data after removing the bulk and z-mask, where the dashed black line shows the boundary of the inner mask function; (d) the relative phase-space density perturbation of the best fit spiral.

5.1. Disqualified and dubious data samples

After running our minimisation algorithm, we inspected the results of the 360 data samples by eye. Many data samples had to be removed altogether or be marked as dubious, due to severe selection effects or some other systematic issue.

The most significant reason for disqualifying data samples is that of selection effects. This is most severe in the approximate direction of the Galactic centre, where we had to disqualify the region where |l|≲50 deg altogether. For the majority of these disqualified data samples, no convincing spiral could be seen by eye after fitting the first step of our minimisation procedure (i.e. in plots corresponding to panel c in Fig. 5). We also disqualified a few data samples were the phase-space spiral could be seen by eye, but our algorithm did not manage to fit that spiral correctly, despite being run several times with different initial conditions. For some cases the fitted spiral had some qualitative differences with respect to the spiral seen by eye, in which case we marked them as dubious. Additionally, in some cases our fit agreed reasonably well with a spiral that was somewhat vaguely discernible by eye, in which case we marked those data samples as dubious.

The direction of the Galactic anti-centre is less plagued by selection, but suffers from another type of systematic issue. At a distance of about 2600 pc, the phase-space spiral undergoes a fairly dramatic change in shape over quite small distance scales. This is illustrated in Fig. 6, showing the phase-space spiral seen in the data after removing the bulk and z-mask (corresponding to panel c in Fig. 5), for data samples within the spatial volume defined by l ∈ [170, 190] deg and . If we compare the closest and most distant data samples in Fig. 6, the arm of the phase-space spiral protrudes from either the left or right. It seems possible that this is due to poor data accuracy, as an effect of the kinematics at low heights being blurred. If this is indeed something physical, we expect to find out with a more dedicated analysis using future Gaia data releases. Either way, we did not account for the spiral to have any dependence on the spatial direction parallel to the Galactic plane; this can only be a reasonable approximation if the shape of the phase-space spiral is fairly constant between neighbouring data samples. For this reason, we either disqualified or marked the data sample as dubious in this spatial region.

thumbnail Fig. 6.

Phase-space spiral seen in the direction of the Galactic anti-centre, for data samples with l ∈ [170, 180] degrees (left column) and l ∈ [180, 190] degrees (right column), and in range 2200–3200 pc (increasing from top to bottom). The respective panels show the relative difference between the data histogram and fitted bulk density distribution: M × (d − B)/B. The horizontal axis is identical for all panels.

In summary, we disqualified a total number of 181 data samples, and marked 85 data samples as dubious. That left us with a total number of 94 well behaved data samples, referred to as good data samples below. We show a few examples of dubious data samples in Appendix C.

5.2. The inferred gravitational potential

We present our main results for the vertical gravitational potential in terms of Φ(400 pc) and Φ(500 pc), because the potential in this height range was found to be the most robustly inferred quantity in our tests on simulations in Paper I. In Figs. 7 and 8, we show the inferred gravitational potential values for our good and dubious data samples, at their respective positions in the (X, Y)-plane. This map is more or less split in half in terms of what data samples produced useful results, where the direction of the Galactic centre is left completely blank. Overall, the dubious data samples agree quite well with the general spatial dependence of the good data samples. As expected, there is a clear trend of lower values with greater Galactocentric radius (i.e. in the direction of negative X). We also see some variations as a function of azimuth, where the direction of l ≃ 225 deg (bottom left quadrant) has somewhat lower values compared to l ≃ 135 deg (top left quadrant); this is discussed further below and in Sect. 6. We also show results for other heights in Appendix B, as well as the inferred vertical acceleration and the inferred time since the perturbation was produced (the model parameter t).

thumbnail Fig. 7.

Inferred gravitational potential at a height of |z| = 400 pc, for the respective data samples. Disqualified data samples are left blank and dubious data samples are marked with a white cross.

thumbnail Fig. 8.

Same as Fig. 7, but for |z| = 500 pc.

The maps shown in Figs. 7 and 8 are fairly smooth. If we compare neighbouring good data samples (in both spatial directions), the relative difference in Φ(500 pc) between them has a median and mean value of 10% and 13%, respectively (where the latter is significantly higher due to a few strong outliers). This could very well be explained, at least in part, by intrinsic variability in the vertical gravitational potential between the spatial locations of the respective data samples. Perhaps most importantly, the cold gas present in the Galactic disk, constituting roughly one fourth of the thin disk surface density (McKee et al. 2015), is highly structured on smaller spatial scales. The dust distribution, which traces the most significant gas component of cold atomic gas, has significant structure on scales of around 100 pc (Lallement et al. 2003, 2019), as does the second-most significant component of cold molecular gas (Dame et al. 2001). Considering that our method has a statistical accuracy of a few per cent at best, which is what we had when we applied it to simulations in Paper I, the difference in gravitational potential values between neighbouring data samples is reasonable, at the very least when discounting the few strong outliers.

The dependence of Φ(400 pc) and Φ(500 pc) with respect to Galactocentric radius is visible in Figs. 9 and 10. The circular markers, representing good data samples, are colour coded according to their spatial Y coordinate, which highlights the broken axisymmetry seen at higher Galactocentric radii. The difference between these regions is of the order of 10–20%. Furthermore, the results for the immediate solar neighbourhood in Paper II, shown as a diamond marker in Figs. 9 and 10, have a smaller value with respect to the data samples of this work that are at a similar Galactocentric radius, with a similar relative difference. These variations are roughly consistent with our 10% estimate of possible systematic biases that could be present in our study (see in Sect. 3.3). When comparing the results presented in Figs. 9 and 10, as well as the supplementary plots in Appendix B, its clear that the variations as a function of azimuth are larger (in a relative sense) for gravitational potential values at lower heights. For this reason, it seems plausible that our results could be biased by a warping of the Galactic disk; this is further discussed in Sect. 6.

thumbnail Fig. 9.

Inferred gravitational potential at a height of |z| = 400 pc, as a function of Galactocentric radius. The results of the good data samples are marked with circles, which are colour coded according to the data samples’ respective Y coordinate. The dashed line shows a best fit exponential curve with respect to the good data samples. We also show the results coming from the dubious data samples, and the result of Paper II for the immediate solar neighbourhood.

thumbnail Fig. 10.

Same as Fig. 9, but for |z| = 500 pc.

To the good data samples, we fitted an analytic function proportional to exp(−R/hL), where R is the Galactocentric radius and hL is a disk scale length. We also fitted separate curves for the groups of good data samples with either only positive or only negative Y-coordinates. These best fitted curves are seen in Figs. 9 and 10 as grey lines. Their respective inferred scale lengths–labelled hL, hL, Y > 0, and hL, Y < 0–are also written out in the figure legends, where the uncertainty comes from assuming a measurement uncertainty for all data samples that makes the scaled χ2 value equal to unity. Clearly, when choosing different spatial cuts, in this case splitting the spatial volume in half along the solar azimuth, we obtain very discrepant results, regardless of what gravitational potential or vertical acceleration value we consider.

6. Discussion

In this work, we have weighed the Galactic disk for distances in the range 1.4–3.4 kpc, with a high spatial resolution in the directions parallel to the Galactic plane, using the Gaia EDR3 proper motion sample. We were able to relax otherwise common assumptions about the Milky Way’s large scale gravitational potential and matter density distribution, most importantly the assumptions of axisymmetry and the exponential decay of the disk mass as a function of Galactocentric radius.

Despite using the Gaia EDR3 proper motion sample, the data samples we constructed were still subject to severe selection effects, due to the combination of their large spatial distance and our cuts in data quality. For this reason, we expanded our model of inference by adding a simple data-driven extinction function that models incompleteness as a function of height (written Ξ, see Sect. 4.1). While our modelling of incompleteness effects was far from perfect, it was largely sufficient for extracting the shape of the phase-space spiral, which is what we need in order to infer the gravitational potential. In this manner, our method differs qualitatively from traditional methods that are based on the assumption of a steady state; such methods are highly sensitive to selection and the gravitational potential that they infer can only be as accurate as the underlying completeness model. This fundamental difference and robustness of our method is exemplified in Paper II, where nearby data samples (100 pc wide in Galactocentric radius) produced stable results despite very varied selection effects (for example, the scale height of the included stars differed by up to 50% between data samples). We could only apply our complete method of inference to 179 out of 360 data samples, where the majority of disqualified data samples were in the direction of the Galactic centre, where stellar crowding and dust extinction is most severe. Out of the 179 data samples, an additional 85 were marked as dubious when comparing the observed and fitted spiral by eye (see Sect. 5.1 for details).

The first step of our minimisation algorithm, where we fit the bulk density distribution, was applied to all 360 data samples. Thus we could infer the vertical velocity of the Sun with respect to the bulk, written W⊙, for the full region of the Galactic disk, which agreed well with the results of more dedicated analyses (e.g. Gaia Collaboration 2018b; Martinez-Medina et al. 2022). In the studied region, the bulk vertical velocity relative to that of the immediate solar neighbourhood (the quantity shown in Fig. 4) lies in the range [ − 2.4, 2.8] km s−1. However, the data samples with the most significant outlier values are either disqualified or marked as dubious; for our good data samples, the same quantity lies in the range [ − 1.0, 1.6] km s−1. As such, it seems like we are constraining ourselves to a spatial region of the Galactic disk with less significant warping.

In terms of the inferred gravitational potential, we present our results in terms of Φ(400 pc) and Φ(500 pc), which were found to be the most robustly inferred quantities in our tests on simulations in Paper I. Because it seems like our results are somewhat biased especially at lower heights, we consider Φ(500 pc) to be more robust for this specific work. The gravitational potential value is close to proportional to the total surface density of the thin disk; this linear relationship holds to the extent that the general shape of the matter density distribution as a function of height is the same for all relevant data samples (i.e. they differ only in terms of their respective amplitudes). In terms of our inferred results, Φ(500 pc) is proportional to the inferred total surface density within |z|< 300 pc, to a relative precision of 0.5%. Due to this strong linear relation, we use the quantities Φ(500 pc) and thin disk surface density more or less interchangeably in the remainder of this paper, as well as in the abstract.

Our results are not perfectly axisymmetric, but have relative variations of the order of 10–20% with respect to the azimuth. This is seen most clearly for data samples at the Galactocentric radius of around 10 kpc, but also in comparison with the immediate solar neighbourhood as analysed in Paper II. This variation is roughly the same order of magnitude as the considered systematic biases, evaluated to roughly 10% in Sect. 3.3. Another interesting facet of our results, seen even more clearly in the supplementary plots of Appendix B, is that these discrepancies are more severe (in a relative sense) when comparing gravitational potential values at smaller heights. This could be explained by a systematic bias coming from an erroneous assumption of a perfectly flat Galactic plane, which would affect the inferred potential especially at low heights. In summary, the broken axisymmetry that we observe could well be explained by a systematic bias, most likely due to variations in height of the Galactic disk mid-plane. For future analysis similar to this work, we would need to include, and possibly produce ourselves, detailed maps of vertical variations to the Galactic disk mid-plane.

Using the inferred value of Φ(500 pc) and assuming an exponential decay as a function of Galactocentric radius, we infer a thin disk scale length of 2.2 ± 0.1 kpc. This is in decent agreement with previous studies by for example Bovy & Rix (2013, 2.15 ± 0.14 kpc), Bland-Hawthorn & Gerhard (2016, 2.6 ± 0.2 kpc), and Mackereth et al. (2017, 1.9 ± 0.1 kpc). However, if we fit the scale length to the spatial volumes with positive or negative Y-coordinates separately, we get very discrepant values of hL, Y > 0 = 2.7 ± 0.2 kpc and hL, Y < 0 = 1.9 ± 0.1 kpc, showing how simplistic assumptions can produce biased results, which in this case was highly sensitive to the chosen spatial volume. This highlights the importance of focusing on accuracy and not only precision, given the enormous statistical power that is granted in the current Gaia era. In order to not be mislead or misleading when it comes to measured properties of the Milky Way, it is crucial to not simply take such numbers at face value, but to carefully reflect on what those numbers mean in the context that they were produced and how they fit into the broader picture of ongoing and future analysis. In a similar vein, symmetry and equilibrium assumptions can bias measurements of for example stars’ angle-action coordinates (e.g. Beane et al. 2019) and the phase-space position of the Sun with respect to a Galactic rest-frame (e.g. Bland-Hawthorn & Gerhard 2016 discuss measurements and potential biases to the height of the Galactic mid-plane and the local standard of rest).

We plan to revisit and extend this analysis with future Gaia data releases. The full third data release is planned for the first half year of 2022 and will contain 33 million radial velocity measurements1 (compared to the current 7.6 million). Additionally, there are complementary distance information such as the photo-astrometric distance measurements produced with StarHorse, where Anders et al. (2022) claim a distance precision as good as 3%. With these improvements to the data, we will be able to apply this method with significantly greater depth, precision, and accuracy. This would give us a better understanding of the large-scale structure of the Galactic disk and the Milky Way’s dynamics and history. This information will also propagate into other more global dynamical mass measurements (e.g. circular velocity curve analysis) and assist in constraining both the large-scale distribution of baryons as well as dark matter. In the longer term, we might be able to resolve variations in the Galactic potential also on smaller spatial scales, for example sourced by the Milky Way’s spiral arms or variations intrinsic to the nature of dark matter.

7. Conclusion

This article is the third part of a longer series about a new method for weighing the Galactic disk, by using the time-varying dynamical structure of the phase-space spiral. In this work, we have applied our method to the Gaia EDR3 proper motion sample in distant regions (1.4–3.4 kpc) of the Galactic disk. We can observe the spiral at this great depth without requiring radial velocity measurements, due to the fact that distant disk stars have a small Galactic latitude (to absolute value), such that their latitudinal proper motion has a close projection with vertical velocity.

This is the first analysis of its kind, in the sense that we are able to weigh the Galactic disk at large distances with a high spatial resolution. In our inference, we do not make strong assumptions about the spatial dependence of the Galactic disk surface density, for example in terms of axisymmetry or an exponential decay with Galactocentric radius. In our post-inference results, we do observe a decay of the disk surface density and fit a disk scale length of 2.2 ± 0.1 kpc. We also observe variations in the surface density on smaller spatial scales, of the order of 10–20%, which are possibly explained by systematic biases. Given these smaller scale variations, it’s clear that a different spatial cut would give a significantly different result for the fitted scale length.

We plan to revisit the analysis made in this paper with future Gaia data releases. This work is the first of its kind and although our analysis seems to suffer from some uncontrolled systematics of the order of roughly 10%, we have demonstrated that we can produce useful results for distant regions of the Galactic disk. We are confident that we will be able to improve and expand our analysis with Gaia’s full third data release, also using supplementary photo-parallax information from the recently updated StarHorse catalogue (Anders et al. 2022), in order to go even further in distance and achieve greater accuracy.

Acknowledgments

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). 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 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 et al. 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. Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360 [Google Scholar]
  5. Bahcall, J. N. 1984a, ApJ, 287, 926 [Google Scholar]
  6. Bahcall, J. N. 1984b, ApJ, 276, 169 [Google Scholar]
  7. Beane, A., Sanderson, R. E., Ness, M. K., et al. 2019, ApJ, 883, 103 [Google Scholar]
  8. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [Google Scholar]
  9. Bertone, G., & Tait, T. M. P. 2018, Nature, 562, 51 [Google Scholar]
  10. Bienayme, O., Soubiran, C., Mishenina, T. V., Kovtyukh, V. V., & Siebert, A. 2006, A&A, 446, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Binney, J., & McMillan, P. 2011, MNRAS, 413, 1889 [NASA ADS] [CrossRef] [Google Scholar]
  12. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, Second Edition (Princeton University Press) [Google Scholar]
  13. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  14. Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115 [Google Scholar]
  15. Buch, J., Leung, J. S. C., & Fan, J. 2019, JCAP, 2019, 026 [Google Scholar]
  16. Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513 [Google Scholar]
  17. Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291 [Google Scholar]
  18. Chen, X., Wang, S., Deng, L., et al. 2019, Nat. Astron., 3, 320 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cole, D. R., & Binney, J. 2017, MNRAS, 465, 798 [Google Scholar]
  20. Crézé, M., Chereul, E., Bienayme, O., & Pichon, C. 1998, A&A, 329, 920 [NASA ADS] [Google Scholar]
  21. Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792 [Google Scholar]
  22. Darling, K., & Widrow, L. M. 2019, MNRAS, 484, 1050 [CrossRef] [Google Scholar]
  23. de Salas, P. F., & Widmark, A. 2021, Rep. Prog. Phys., 84 [Google Scholar]
  24. Deason, A. J., Erkal, D., Belokurov, V., et al. 2021, MNRAS, 501, 5964 [Google Scholar]
  25. Debattista, V. P., Gonzalez, O. A., Sanderson, R. E., et al. 2019, MNRAS, 485, 5073 [Google Scholar]
  26. Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429 [Google Scholar]
  27. Deng, L.-C., Newberg, H. J., Liu, C., et al. 2012, Res. Astron. Astrophys., 12, 735 [Google Scholar]
  28. D’Onghia, E., Vogelsberger, M., & Hernquist, L. 2013, ApJ, 766, 34 [Google Scholar]
  29. Ferreira, E. G. M. 2021, A&ARv., 29, 7 [NASA ADS] [CrossRef] [Google Scholar]
  30. Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Gaia Collaboration (Katz, D., et al.) 2018b, A&A, 616, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Gaia Collaboration (Antoja, T., et al.) 2021a, A&A, 649, A8 [EDP Sciences] [Google Scholar]
  33. Gaia Collaboration (Smart, R. L., et al.) 2021b, A&A, 649, A6 [EDP Sciences] [Google Scholar]
  34. Garbari, S., Liu, C., Read, J. I., & Lake, G. 2012, MNRAS, 425, 1445 [Google Scholar]
  35. Gilmore, G., Randich, S., Asplund, M., et al. 2012, The Messenger, 147, 25 [NASA ADS] [Google Scholar]
  36. Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [Google Scholar]
  37. GRAVITY Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Guo, R., Liu, C., Mao, S., et al. 2020, MNRAS, 495, 4828 [Google Scholar]
  39. Haines, T., D’Onghia, E., Famaey, B., Laporte, C., & Hernquist, L. 2019, ApJ, 879, L15 [Google Scholar]
  40. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  41. Hohl, F. 1971, ApJ, 168, 343 [NASA ADS] [CrossRef] [Google Scholar]
  42. Holmberg, J., & Flynn, C. 2000, MNRAS, 313, 209 [Google Scholar]
  43. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  44. Jungman, G., Kamionkowski, M., & Griest, K. 1996, Phys. Rept., 267, 195 [Google Scholar]
  45. Juric, M., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [Google Scholar]
  46. Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2014, ApJ, 794, 59 [Google Scholar]
  47. Kapteyn, J. C. 1922, ApJ, 55, 302 [NASA ADS] [CrossRef] [Google Scholar]
  48. Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Klasen, M., Pohl, M., & Sigl, G. 2015, Prog. Part. Nucl. Phys., 85, 1 [Google Scholar]
  50. Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597 [Google Scholar]
  51. Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571 [NASA ADS] [CrossRef] [Google Scholar]
  52. Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75 [Google Scholar]
  53. Lallement, R., Welsh, B. Y., Vergely, J. L., Crifo, F., & Sfeir, D. 2003, A&A, 411, 447 [EDP Sciences] [Google Scholar]
  54. Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286 [Google Scholar]
  56. Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134 [Google Scholar]
  57. Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586 [CrossRef] [Google Scholar]
  58. Li, Z.-Z., Qian, Y.-Z., Han, J., et al. 2020, ApJ, 894, 10 [Google Scholar]
  59. Lindegren, L., Bastian, U., Biermann, M., et al. 2021a, A&A, 649, A4 [EDP Sciences] [Google Scholar]
  60. Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021b, A&A, 649, A2 [EDP Sciences] [Google Scholar]
  61. Mackereth, J. T., Bovy, J., Schiavon, R. P., et al. 2017, MNRAS, 471, 3057 [Google Scholar]
  62. Martinez-Medina, L., Pérez-Villegas, A., & Peimbert, A. 2022, MNRAS, 512, 1574 [Google Scholar]
  63. McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13 [Google Scholar]
  64. McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56 [Google Scholar]
  65. McMillan, P. J. 2011, MNRAS, 414, 2446 [Google Scholar]
  66. McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
  67. Minchev, I., Quillen, A. C., Williams, M., et al. 2009, MNRAS, 396, L56 [NASA ADS] [CrossRef] [Google Scholar]
  68. Ness, M., Freeman, K., Athanassoula, E., et al. 2013, MNRAS, 430, 836 [NASA ADS] [CrossRef] [Google Scholar]
  69. Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001 [Google Scholar]
  70. Oort, J. H. 1932, Bull. Astron. Inst. Netherlands, 6, 249 [Google Scholar]
  71. Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 500, 501 [NASA ADS] [Google Scholar]
  72. Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21 [Google Scholar]
  73. Purcell, C. W., Bullock, J. S., Tollerud, E. J., Rocha, M., & Chakrabarti, S. 2011, Nature, 477, 301 [Google Scholar]
  74. Quillen, A. C., Pettitt, A. R., Chakrabarti, S., et al. 2020, MNRAS, 499, 5623 [CrossRef] [Google Scholar]
  75. Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 063101 [Google Scholar]
  77. Salomon, J.-B., Bienaymé, O., Reylé, C., Robin, A. C., & Famaey, B. 2020, A&A, 643, A75 [EDP Sciences] [Google Scholar]
  78. Sanders, J. L., & Das, P. 2018, MNRAS, 481, 4093 [CrossRef] [Google Scholar]
  79. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [Google Scholar]
  80. Schutz, K., Lin, T., Safdi, B. R., & Wu, C.-L. 2018, Phys. Rev. Lett., 121 [CrossRef] [Google Scholar]
  81. Sivertsson, S., Read, J. I., Silverwood, H., et al. 2022, MNRAS, 511, 1977 [NASA ADS] [CrossRef] [Google Scholar]
  82. Skowron, D. M., Skowron, J., Mróz, P., et al. 2019, Acta Astron., 69, 305 [NASA ADS] [Google Scholar]
  83. Velazquez, H., & White, S. D. M. 1999, MNRAS, 304, 254 [NASA ADS] [CrossRef] [Google Scholar]
  84. Villalobos, Á., & Helmi, A. 2008, MNRAS, 391, 1806 [Google Scholar]
  85. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  86. Weber, M., & de Boer, W. 2010, A&A, 509, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Widmark, A., de Salas, P. F., & Monari, G. 2021a, A&A, 646, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Widmark, A., Laporte, C., & de Salas, P. F. 2021b, A&A, 650, A124 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Widmark, A., Laporte, C., de Salas, P. F., & Monari, G. 2021c, A&A, 653, A86 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Widmark, A., Hunt, J. A. S., Laporte, C. F. P., & Monari, G. 2022, A&A, 663, A16 (Paper IV) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239 [Google Scholar]
  92. Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41 [Google Scholar]
  93. Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [Google Scholar]
  94. Xu, Y., Liu, C., Tian, H., et al. 2020, ApJ, 905, 6 [NASA ADS] [CrossRef] [Google Scholar]
  95. Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377 [Google Scholar]
  96. Yao, J. M., Manchester, R. N., & Wang, N. 2017, MNRAS, 468, 3289 [Google Scholar]

Appendix A: Coordinate transformations

In this appendix, we describe how we transform the Gaia observables to the solar rest-frame Cartesian coordinates X and V, as defined in the beginning of Sect. 2. The Gaia observables include the parallax ϖ, the Galactic longitude and latitude (l and b), their corresponding proper motions (μl and μb), and the radial velocity (vRV).

The spatial position is given by

(A.1)

The height with respect to the Galactic mid-plane is then found via a translation according to Eq. (1).

The velocities in the solar rest frame are given by

(A.2)

where kμ = 4.74057 yr × km s−1 and

(A.3)

is a rotational matrix. The vertical velocity in the rest frame of the Galactic disk is then found via a translation according to Eq. (2).

There is further important information on how the data is processed, such as accounting for the Gaia EDR3 zero-point offset and assigning missing radial velocities. This is described in detail in Sect. 3.

Appendix B: Supplementary plots

In this appendix, we show some additional plots of our results. In Figs. B.1–B.5, we show how the inferred vertical gravitational potential values Φ(300 pc) and Φ(600 pc), as well as the inferred vertical acceleration values Kz(200 pc), Kz(300 pc), and Kz(400 pc), depend on Galactocentric radius. Their azimuthal dependence is also visible via the marker colour coding. In all these figures, we see very similar disk scale lengths and a similar general behaviour when splitting the data samples into two separate groups along the solar azimuth.

thumbnail Fig. B.1.

Same as Fig. 9, but for |z| = 300 pc.

thumbnail Fig. B.2.

Same as Fig. 9, but for |z| = 600 pc.

thumbnail Fig. B.3.

Same as Fig. 9, but for the vertical acceleration at |z| = 200 pc.

thumbnail Fig. B.4.

Same as Fig. 9, but for the vertical acceleration at |z| = 300 pc.

thumbnail Fig. B.5.

Same as Fig. 9, but for the vertical acceleration at |z| = 400 pc.

In Fig. B.6 we show the inferred time since the perturbation (the model parameter t). As we saw in our tests on simulations in Paper I, this parameter is not very robustly inferred, due to its strong degeneracy with respect to the precise shape of the gravitational potential. We do in fact see a manifestation of this degeneracy, in that the inferred time has a strong correlation with the broken axisymmetry of the inferred gravitational potential at low heights (seen most clearly in Φ(300 pc)). This is further indication that the observed broken axisymmetry is likely due to a systematic bias. We also see two very strong outliers in the direction of l ∈ [100, 110] deg. The inferred gravitational potentials of these two data samples are close to harmonic, which pushes the inferred time to high values. Due to these two outliers, the distribution of inferred times is best captured by its 16th, 50th, and 84th percentiles; for our good data samples, they are 350, 622, and 767 Myr, respectively. This is in decent agreement with the results of for example Antoja et al. (2018) and Darling & Widrow (2019).

thumbnail Fig. B.6.

Inferred time since the perturbation was produced (model parameter t) for the respective data samples. Disqualified data samples are left blank and dubious data samples are marked with a white cross.

Appendix C: Examples of dubious data samples

In this appendix, we show a few examples of dubious data samples, in terms of their data histograms and inferred phase-space densities. This is done in order to illustrate our rationale for disqualifying and marking data samples as dubious.

In Fig. C.1, we show the data histogram and inferred phase-space density of the data sample with l ∈ [60, 70] deg and . This data sample is representative of the larger spatial region of l ∈ [60, 90] deg, which was all disqualified or marked as dubious. The data samples in this region suffer from fairly strong extinction but also a skewed spiral, as seen in panel (c) of this figure. The skewed shape of the spiral is likely produced by some systematic bias, plausible related to the assignment of radial velocities, as that could indeed skew the spiral in this manner (see Fig. 3). Even though this region (l ∈ [60, 90] deg) had no good data samples, its dubious data samples still produced reasonable results that agreed well with our inferred disk scale length (see especially Figs. 9 and 10).

thumbnail Fig. C.1.

Same as Fig. 5, but for the data sample with l ∈ [60, 70] deg and . This data sample was marked as dubious.

In Fig. C.2, we show the data histogram and inferred phase-space density of the data sample with l ∈ [270, 280] deg and , which is representative of the region with the same Galactic longitude. Looking at panel (c), the phase-space spiral is clearly present in the data as a perturbation with respect to the fitted bulk density. However, the velocity distribution does seem rather asymmetric and skewed (in addition to the spiral perturbation), with a very red (blue) region at the bottom (top) of the panel. If this is due to an actual skewed velocity distribution or some systematic bias is difficult to say, but either way it lead us to mark the data samples in this region as dubious. Again, most data samples in this region agree fairly well with the general trends of our results.

thumbnail Fig. C.2.

Same as Fig. 5, but for the data sample with l ∈ [270, 280] deg and . This data sample was marked as dubious.

In Fig. C.3, we show the data histogram and inferred phase-space density of the data sample with l ∈ [130, 140] deg and . This is an example of a data sample which overall seemed well behaved but suffered from fairly strong selection effects. The inferred spiral seen in panel (d) seems like a good fit, but the spiral in panel (c) is not seen as a clearly continuous structure. This lead us to mark this data sample as dubious.

thumbnail Fig. C.3.

Same as Fig. 5, but for the data sample with l ∈ [130, 140] deg and . This data sample was marked as dubious.

All Tables

Table 1.

Free parameters in our model of inference.

In the text

All Figures

thumbnail Fig. 1.

Stellar number counts in our 360 data samples, shown in their respective locations in the (X, Y)-plane. The Galactic centre is located to the right, while the direction of Galactic rotation is upwards.
In the text
thumbnail Fig. 2.

Radial velocity (vRV) as a function of Galactic latitude (b) for the data sample with l ∈ [90, 100] deg and . The dashed line corresponds to the fitted function and the scatter point correspond to stars with available radial velocity information. The range in b on the horizontal axis corresponds to a range in height spanning roughly Z ∈ [ − 1, 1] kpc at the relevant distance.
In the text
thumbnail Fig. 3.

Shape of the phase-space spiral, represented as a one-dimensional line in the (z, w)-plane, when biased by faulty radial velocity assignments. The dashed blue and orange lines show the spiral when radial velocity assignment is biased by ±10 km s−1, assuming that the spiral is observed at a distance of . The spiral is plotted for vertical energies in range Ez ∈ [Φ(250 pc),Φ(800 pc)]. The plus sign marks the origin of the (z, w)-plane.
In the text
thumbnail Fig. 4.

Bulk vertical velocity relative to that of the immediate solar neighbourhood. This is equivalent to W⊙, local − W⊙, where W⊙, local = 7.25 km s−1 and W⊙ is a free parameter fitted in the first step of our minimisation algorithm. The Galactic centre is located to the right, while the direction of Galactic rotation is upwards.
In the text
thumbnail Fig. 5.

Data and fitted phase-space density of the data sample with l ∈ [120, 130] deg and . All four panels span the same range of the (z, w)-plane and show: (a) the data histogram; (b) the fitted bulk and z-mask; (c) the spiral as seen in the data after removing the bulk and z-mask, where the dashed black line shows the boundary of the inner mask function; (d) the relative phase-space density perturbation of the best fit spiral.
In the text
thumbnail Fig. 6.

Phase-space spiral seen in the direction of the Galactic anti-centre, for data samples with l ∈ [170, 180] degrees (left column) and l ∈ [180, 190] degrees (right column), and in range 2200–3200 pc (increasing from top to bottom). The respective panels show the relative difference between the data histogram and fitted bulk density distribution: M × (d − B)/B. The horizontal axis is identical for all panels.
In the text
thumbnail Fig. 7.

Inferred gravitational potential at a height of |z| = 400 pc, for the respective data samples. Disqualified data samples are left blank and dubious data samples are marked with a white cross.
In the text
thumbnail Fig. 8.

Same as Fig. 7, but for |z| = 500 pc.
In the text
thumbnail Fig. 9.

Inferred gravitational potential at a height of |z| = 400 pc, as a function of Galactocentric radius. The results of the good data samples are marked with circles, which are colour coded according to the data samples’ respective Y coordinate. The dashed line shows a best fit exponential curve with respect to the good data samples. We also show the results coming from the dubious data samples, and the result of Paper II for the immediate solar neighbourhood.
In the text
thumbnail Fig. 10.

Same as Fig. 9, but for |z| = 500 pc.
In the text
thumbnail Fig. B.1.

Same as Fig. 9, but for |z| = 300 pc.
In the text
thumbnail Fig. B.2.

Same as Fig. 9, but for |z| = 600 pc.
In the text
thumbnail Fig. B.3.

Same as Fig. 9, but for the vertical acceleration at |z| = 200 pc.
In the text
thumbnail Fig. B.4.

Same as Fig. 9, but for the vertical acceleration at |z| = 300 pc.
In the text
thumbnail Fig. B.5.

Same as Fig. 9, but for the vertical acceleration at |z| = 400 pc.
In the text
thumbnail Fig. B.6.

Inferred time since the perturbation was produced (model parameter t) for the respective data samples. Disqualified data samples are left blank and dubious data samples are marked with a white cross.
In the text
thumbnail Fig. C.1.

Same as Fig. 5, but for the data sample with l ∈ [60, 70] deg and . This data sample was marked as dubious.
In the text
thumbnail Fig. C.2.

Same as Fig. 5, but for the data sample with l ∈ [270, 280] deg and . This data sample was marked as dubious.
In the text
thumbnail Fig. C.3.

Same as Fig. 5, but for the data sample with l ∈ [130, 140] deg and . This data sample was marked as dubious.
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.