© The Authors 2024

1. Introduction

In recent years, observations from the Gaia satellite (Gaia Collaboration 2016) have revolutionised our understanding of the structure, history, and composition of our Galaxy. Gaia has measured the positions and motions of more than a billion stars, enabling unprecedented analyses of Milky Way stellar dynamics and stellar population properties.

As an example of this progress, many recent studies have sought a better understanding of the spiral structure of the Milky Way (Xu et al. 2018; Shen & Zheng 2020), which can yield insights into the Galaxy’s formation history. The spiral arms have been probed using various tracers, such as hydrogen gas (Levine et al. 2006), masers (Xu et al. 2013; Reid et al. 2014, 2019), and young stars (Poggio et al. 2021; Castro-Ginard et al. 2021; Lin et al. 2022; Gaia Collaboration 2023a,b). Different studies report slightly different arm morphologies in terms of pitch angle and spatial extent (e.g. comparing Reid et al. 2019 and Poggio et al. 2021; sometimes such differences can, at least in part, be explained by intrinsic differences between tracers, as discussed by Miyachi et al. 2019). The spiral arms have also been associated with various kinematic signatures (Eilers et al. 2020; Widmark et al. 2022b; Khoperskov & Gerhard 2022; Martinez-Medina et al. 2022; Palicio et al. 2023).

Despite the aforementioned advances in data precision and depth, the Milky Way’s spiral structure has not been detected in a dynamical mass measurement¹. To detect such variations in the disk surface density requires precise measurements of distant disk regions. This is challenging, mainly because of stellar crowding and dust extinction, which cause strong selection effects and other data biases that are difficult to model. Data incompleteness is especially severe for the radial (i.e. line-of-sight) velocity measurements, which are only available for a bright subsample of stars.

In the present work we used state-of-the-art data processing techniques to overcome these obstacles and performed a vertical Jeans analysis of the thin disk out to a distance of 3 kpc from the Sun. The vertical Jeans analysis was formulated a century ago by Jeans (1922) and remains a widely used probe of the dynamical properties of the solar neighbourhood, such as the local density of dark matter (Sivertsson et al. 2018; Salomon et al. 2020; Guo et al. 2020; de Salas & Widmark 2021; Widmark et al. 2021a). Under this technique, one assumes dynamical equilibrium and that stellar motions in the directions parallel and perpendicular to the disk plane can be approximately decoupled. Then, a stellar tracer population’s number density distribution and vertical velocity distribution are interrelated via the vertical gravitational potential.

Using the Gaia Data Release 3 (DR3) catalogue, we constructed four separate data samples of bright stars, defined by cuts in absolute magnitude. We divided the disk plane into a square grid with a 100 pc spacing and applied the vertical Jeans analysis to each area cell and data sample separately. For the two key ingredients of vertical Jeans analysis, we did the following: For the stellar number density field, we used the 3D maps produced by Widmark et al. (2022b) using Gaussian processes (GPs), based on Gaia and StarHorse spectro-astrometric distances (Anders et al. 2022). For the vertical velocity fields, we used Gaia measurements supplemented by Bayesian neural network (BNN) predictions of the radial velocity (i.e. line-of-sight velocity) from Naik & Widmark (2024)². These radial velocity predictions are especially useful for our application because radial velocity is a sub-dominant component of the vertical velocity in distant regions of the disk, allowing us to make full use of the highly informative Gaia proper motion sample. These data products, in combination, allowed us to make a detailed map of how the surface density varies in the directions parallel to the Galactic plane.

For the first time, we find evidence for the imprint of spiral structure in the dynamically measured vertical gravitational potential of the Milky Way disk. In all four of our data samples, we see the same non-axisymmetric over-density, which is very similar to that of the Local Arm as traced by stellar age and chemistry (e.g. Poggio et al. 2021).

2. Method

2.1. Coordinate system and data samples

We used solar rest-frame Cartesian coordinates X = (X, Y, Z) pointing in the directions of the Galactic centre, Galactic rotation, and Galactic north. The vertical coordinate in the disk frame is equal to

$$\begin{array}{c}\hfill z=Z+Z_{\odot },\end{array}$$(1)

where Z_⊙ is the height of the Sun with respect to the Galactic disk mid-plane. The vertical velocity in the solar rest frame is W ≡ dZ/dt. The galactocentric radius is given by

$$\begin{array}{c}\hfill R=\sqrt{{(R_{\odot }-X)}^2+Y^2},\end{array}$$(2)

where R_⊙ = 8.2 kpc (McMillan 2016). The galactocentric angle, ϕ, fulfils

$$\begin{array}{c}\hfill sin\varphi =\frac{Y}{R},cos\varphi =\frac{R_{\odot }-X}{R}.\end{array}$$(3)

We constructed four separate stellar samples by making cuts in the Gaia G-band absolute magnitude, according to M_G ∈ (0, 1], (1, 2], (2, 3], and (3, 4] (the same as in Widmark et al. 2022b). We divided the (X, Y) plane into a square grid with a cell spacing of 100 pc, for the spatial volume within $\sqrt{X^2+Y^2}<3\mathrm{kpc}$.

2.2. Vertical Jeans analysis

The vertical Jeans equation gives

$$\begin{array}{c}\hfill \frac{1}{n}\frac{\partial (n{\sigma }_W^2)}{\partial Z}+\frac{\partial \mathrm{\Phi }}{\partial Z}=0,\end{array}$$(4)

where Φ is the gravitational potential, n and ${\sigma }_W^2$ are the number density and vertical velocity variance of a stellar tracer population, and Z is the direction perpendicular to the Galactic plane.

In the equation above, we made a few simplifying assumptions. In particular, we assumed dynamical equilibrium and dropped a term representing a coupling between radial and vertical motions (the ‘tilt’ term), which gives a sizable contribution only at heights closer to 1 kpc (Read 2014).

2.3. Stellar number density field

For the stellar number count density, we used the results from Widmark et al. (2022b), where the stellar number density fields of the four data samples were modelled with GP regression. This was based on StarHorse spectro-astrometric distances (Anders et al. 2022), using data from Gaia DR3 (Gaia Collaboration 2023c), Pan-STARRS1 (Scolnic et al. 2015), SkyMapper (Casagrande et al. 2019), Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006), and Wide-field Infrared Survey Explorer mission (AllWISE, Marrese et al. 2019).

As discussed in Sect. 1, selection effects are a significant obstacle to dynamical mass measurements in distant regions of the disk. The Gaia selection function is highly complex, and depends on brightness limits, stellar crowding and the Gaia scanning law (Boubert & Everall 2020; Everall & Boubert 2022; Rybizki et al. 2021; Cantat-Gaudin et al. 2023). Furthermore, the stellar number density field can be biased by open clusters and dust extinction. All these issues were carefully considered and mitigated in Widmark et al. (2022b). To avoid issues with the bright and dim end, they limited the G-band apparent magnitudes to the range of roughly 6–17 mag, where the Gaia catalogue is very close to complete, with deviations of a few percent at most (Everall & Boubert 2022; Cantat-Gaudin et al. 2023). Furthermore, spatial volumes were carefully masked to avoid open clusters and incompleteness issues associated with dust extinction. By using GPs, the stellar number density field was modelled as a smooth and differentiable function. This approach was also fairly model-independent, and did not impose constraints such as Galactic axisymmetry.

The GP regression used correlation lengths of (l_X, l_Y, l_Z) = (300, 300, 100) pc, which imposes a degree of smoothness in all three spatial directions. In this manner, even though the area cells have a 100 pc spacing in (X, Y) plane, the number density field is informed by a larger area. This is especially useful when a larger fraction of an area cell is masked (e.g. due to an open cluster); in that case, its number density field is still informed by the data in its surrounding spatial volume. However, in order to avoid area cells that are very poorly sampled, we masked area cells that have an average effective volume fraction smaller than 0.5 for |Z|< 250 pc.

We wanted to retain the useful properties of GP regression discussed above. However, there are time-varying perturbations to the stellar number population (which the GP regression encapsulates); such perturbations can create biases in dynamical mass measurements since they are not in a steady state (see Appendix A and Sect. 4 for further details). Therefore, we fitted an analytic function to the stellar number density that was obtained via GP regression, separately for each data sample and area cell, assuming a functional form which is mirror symmetry and monotonically decreasing with height. This analytic function is a mixture model of three disk components, according to

$$\begin{array}{c}\hfill \widehat{n}(Z|a_i,h_i,Z_{\odot })={\displaystyle \sum _{i=1}^3{a}_i\times {\text{cosh}}^{-2}\left(\frac{Z+Z_{\odot }}{h_i}\right),}\end{array}$$(5)

where a_i are the disk component normalisations, h_i are the scale heights, and Z_⊙ is the height of the Sun with respect to the disk mid-plane. This functional form is flexible enough to faithfully model the vertical number density profile of our stellar tracer populations. We used wide flat box priors and the GP regression uncertainties in our fit.

In Fig. 1 we show an example of the stellar number density distribution, in terms of the underlying data counts, GP regression, and symmetric analytic function, for a small area of nine neighbouring area cells of the 2 < M_G ≤ 3 data sample. We chose this particular spatial volume in order to illustrate a few important features. In the top-right panel, the raw data number count has a very large uncertainty, due to open cluster masking. The GP curve is still inferred there, informed by its surrounding spatial volume. In all panels, the stellar over-densities at roughly Z = 200 pc and Z = −400 pc (and their mirrored under-densities) are projections of the phase space spiral (see Fig. 7 in Widmark et al. 2022b). Because the phase-space spiral is a time-varying structure, it constitutes a bias to Jeans analysis, which we averaged away by fitting the symmetric and analytic function $\widehat{n}$.

Fig. 1. Stellar number density distribution for a few area cells of the data sample with 2 < MG ≤ 3. Each panel is labelled by the area cell’s mid-point (X, Y) coordinates. The function obtained by GP regression is shown in solid black; the fitted symmetric function ($\widehat{n}$) of Eq. (5) is shown in dashed red; the raw number count (N) with its associated Poisson noise uncertainty is shown as a 1σ band.

2.4. Vertical velocity field

For the vertical velocity field, we used all Gaia DR3 stars in the unmasked (X, Y) area cells, subdivided into the four data samples. We removed possible members (probability > 0.1) of the open clusters tabulated by Hunt & Reffert (2023) and stars with StarHorse distance precisions worse than 10%. This left 23 551 383 stars overall: 3 262 301, 2 736 077, 5 710 148, and 11 842 857 in each of the four magnitude bins, from brightest to faintest.

Of these stars, around 60% have radial velocity measurements in Gaia DR3 (ranging from 85% in the brightest sample to 50% in the faintest). The spatial distribution of these measurements is shown in Fig. 2, showing the fraction of stars with radial velocity measurements in each of our area cells and for each magnitude sample. At all magnitudes, the stars within approximately 1 kpc are very well covered by radial velocity measurements. Beyond this distance, the radial velocity coverage decreases appreciably with brightness. For the remainder of the stars without radial velocity measurements, we used the Bayesian radial velocity predictions of Naik & Widmark (2024). These prediction posterior distributions are generated using BNNs and are trained on the subset of DR3 stars with radial velocities (Naik & Widmark 2022, 2024). In effect, these BNNs learn a latent representation of the stellar phase space distribution, convolved with any uncertainty due to a lack of training data. As a result, these prediction distributions can be rather wide, with typical widths in the range 20–30 km s⁻¹.

Fig. 2. Fraction of stars with radial velocity measurements across the 100 × 100 pc area cells used in our analysis. The four panels represent the four absolute magnitude samples, as labelled.

These velocity predictions have been validated at several stages. First, Naik & Widmark (2022) used a BNN trained on Gaia DR2 to produce a set of blind predictions for radial velocity measurements that would subsequently be released in Gaia DR3. The predictions were then compared with the measurements in Naik & Widmark (2024), and found to have been very successful: the prediction distributions were not only consistent with the measurements, but they also exhibited appropriate levels of confidence, with inflated prediction uncertainties for stars outside the footprint of the training data. Naik & Widmark (2024) then repeated the exercise and trained a new model on the DR3 data. This newer model was shown to be highly successful in validation tests with DR3 radial velocities excluded from the training data, and also with stars from other spectroscopic surveys. The model was then used to generate a publicly available catalogue of radial velocity predictions for 185 million stars still lacking radial velocity measurements in Gaia DR3. We used these tabulated predictions in this work.

Naik & Widmark (2024) suggest a natural approach for incorporating these predictive posterior distributions into a wider analysis: multiple imputation (Little & Rubin 2014; Gelman et al. 2004). Here, one draws N random realisations of the predictive posterior for each 5D star, to construct N ‘complete’ datasets (or ‘imputations’), then conducts the remainder of the analysis with each imputation separately to infer the parameters of interest. The variation in the parameters of interest across the imputations gives the uncertainty due to the missing data. This is a well-motivated technique for dealing with missing data, as it is formally equivalent to performing a single joint inference over both the missing data and the parameters of interest. In other words, ‘Bayesian shrinkage’ occurs only once for a given set of observations. In the present work, we constructed 20 such imputations (i.e. we generate 20 realisations from the prediction posterior distribution of each star) and then conducted the remainder of our analysis separately for each imputation.

For each imputation, data sample, and area cell, we converted the measured phase space coordinates (sky positions, StarHorse distances, proper motions, and radial velocities) to the heliocentric Cartesian coordinate system described above. We then fitted a parametric model that describes ${\sigma }_W^2$ as a function of Z using the data of the area cell as well as all adjacent unmasked cells (effectively smoothing our results over a 300-by-300 pc area):

$$\begin{array}{c}\hfill {\widehat{\sigma }}_W^2(Z)={\sigma }_0^2+{\sigma }_1^{2}ln[1+{\left(\frac{Z+Z_{\odot }}{H_1}\right)}^2]+{\sigma }_2^{2}ln[1+{\left(\frac{Z+Z_{\odot }}{H_2}\right)}^2],\end{array}$$(6)

where ${\sigma }_{\{0,1,2\}}^2$ and H_{1, 2} are dimensionful free parameters in the fit. The quantity Z_⊙ represents the height of the Sun with respect to the disk mid-plane; this is not taken as a free parameter but as a constant, given by the fit of the stellar number density model in that area cell. We fitted the velocity variance model (Eq. (6)) to the data, 𝒟 (comprising stellar heights, Z, vertical velocities, W, and observational vertical velocity uncertainties, σ_obs, propagated from proper motion and distance measurements) by maximising a Gaussian log-likelihood,

$$\begin{array}{c}\hfill ln\mathcal{L}(\mathcal{D}|\overline{W},{\sigma }_{\{0,1,2\}}^2,H_{\{1,2\}})=-\frac{1}{2}{\displaystyle \sum _i[\frac{{(W_i-\overline{W})}^2}{{\sigma }_i^2}+ln\left(2\pi {\sigma }_i^2\right)],}\end{array}$$(7)

where the index i runs over individual stars, and

$$\begin{array}{c}\hfill {\sigma }_i^2\equiv {\sigma }_{\mathrm{obs},i}^2+{\widehat{\sigma }}_W^2(Z_i|{\sigma }_{\{0,1,2\}}^2,H_{\{1,2\}}).\end{array}$$(8)

The mean vertical velocity, $\overline{W}$, is a sixth free parameter, in addition to the five free parameters ${\sigma }_{\{0,1,2\}}^2$, H_{1, 2} entering the dispersion model.

Figure 3 shows a few examples of fitted velocity dispersion profiles. For each of the four data samples, we chose three random area cells and plotted the best-fit model according to Eq. (6), as well as the measured velocity dispersions in adaptively spaced vertical bins. In each case, we individually plotted the data and model for each of the 20 imputations. It is clear from this figure that our model (Eq. (6)) effectively captures the overall shape of the dispersion profiles despite the observed variability in these shapes, while at the same time foregoing smaller-scale features arising from unmixed phase space substructures.

Fig. 3. Examples of vertical velocity dispersions and the best-fitting models (Eq. (6)). The four rows correspond to the four magnitude samples, while the three panels within each row are randomly chosen cells in the (X, Y) plane. The yellow points represent the measured velocity dispersions in adaptively spaced Z bins, while the solid lines show the best-fitting models. Note that the model is not fit to these binned dispersions, but to the individual stellar velocities. In each case, 20 sets of data and 20 models are shown, corresponding to the 20 random imputations of the dataset as described in the main text.

We also tried different models for the vertical velocity dispersion as a function of height. Another functional that gave rise to consistent results (in terms of the agreement between different spatial regions and between data samples) was a quartic polynomial. These fits to the velocity data were somewhat worse in terms of χ²-value and by-eye inspections, but we saw very similar results for the inferred gravitational potential, leading to the same general conclusions.

The uncertainties associated with the velocity distribution are dominant in terms of the inferred gravitational potential. These uncertainties are propagated to the global fits of the exponential function and the simple arm models (see Sect. 3).

3. Results

Our dynamical mass measurement is likely the most robust for the potential difference of roughly Φ(400 pc)−Φ(0 pc). Between these heights, the stellar number density ratio is roughly 10–15% for the three brighter data samples, and 30% for the dimmest data sample. Furthermore, the spiral arms are a thin disk phenomenon; staying closer to the Galactic plane makes the measurement more local and will better resolve the spiral arms, as the gravitational potential at greater heights is influenced by a larger area of the Galactic disk. This potential difference is a good proxy for and closely proportional to the thin disk surface density, which is discussed in Sect. 2.2.

Figure 4 shows our inferred values of Φ(400 pc)−Φ(0 pc) as a function of galactocentric radius, for all four data samples. The markers each correspond to an area cell, where the colour denotes their galactocentric angle, ϕ. Coherent structure visible in the marker colours (denoting the galactocentric angle, ϕ) reveal non-axisymmetric structure in the inferred surface density results; for example, at the solar radius (R = 8.2 kpc) the results for positive ϕ are over-dense compared to negative ϕ, in all four data samples. We also show an exponential function fitted separately to each data sample, which takes the form

$$\begin{array}{c}\hfill f(R)=A\times exp\left(\frac{R-R_{\odot }}{h_{\mathrm{L}}}\right),\end{array}$$(9)

Fig. 4. Gravitational potential, Φ(400 pc)−Φ(0 pc), as a function of galactocentric radius, for the four stellar samples (as labelled in the panels’ top-right corner). The dashed line is the exponential function of Eq. (9) fitted to data points with $\sqrt{X^2+Y^2}<2\mathrm{kpc}$, which excludes distant and dubious spatial regions. The scatter points correspond to area cells in the (X, Y) plane. Area cells that are included in the exponential function fit are coloured according to their galactocentric disk plane angle (ϕ); excluded data cells have light grey and smaller scatter points. The colour bar and axis ranges are shared between all panels.

where A is the solar radius value, and h_L is the disk scale length parameter. The data samples that are more distant ($\sqrt{X^2+Y^2}>2\mathrm{kpc}$) are likely less reliable and therefore masked in this fit. The fitted exponential function parameter values are listed in the respective panels of Fig. 4 and they agree well between the four data samples. We show only the best fit value; the statistical error is small and systematic uncertainties dominate.

The residual of the inferred Φ(400 pc)−Φ(0 pc) as compared to the fitted exponential function is shown in Fig. 5, for all four data samples. The residual of an area cell is defined according to

$$\begin{array}{c}\hfill \text{residual}=\frac{\text{measurement}-f(R)}{f(R)},\end{array}$$(10)

Fig. 5. Residual of Φ(400 pc)−Φ(0 pc) with respect to a fitted exponential, projected on the (X, Y) plane. Each panel corresponds to one of our four data samples. In each panel, dotted grey lines show galactocentric iso-radial and iso-azimuthal lines, and the black plus sign in the centre marks the solar position. The grey contour lines correspond to the spiral arms as mapped by Poggio et al. (2021), as seen in Fig. 6.

where f(R) is the fitted exponential function. The residuals of the four data samples show similar results. Most notably, they have the same over-dense structure, stretching roughly from (X, Y) = (0, 1.5) kpc to (X, Y) = (−1.5, −0.5) kpc. This feature agrees well with the Local Arm revealed in the map of young stars by Poggio et al. (2021), shown in detail in Fig. 6 and as grey contour lines in Fig. 5.

Fig. 6. Over-density of upper main sequence stars in the (X, Y) plane as a tracer of spiral arms. Results are taken directly from Poggio et al. (2021). The grey-scale contour lines are the same as in Fig. 5, corresponding to values ( − 0.5,  0,  0.5).

In order to make quantitative statements about the over-density that we identify with the Local Arm, we fitted a simple analytic model to the residual seen in Fig. 5. We modelled the ridge of the over-density as a 1D curve in the Galactic disk plane, as a logarithmic spiral of the form

$$\begin{array}{c}\hfill \stackrel{\sim }{R}(\varphi )={\stackrel{\sim }{R}}_{0}exp[-\varphi tan(\stackrel{\sim }{\alpha })],\end{array}$$(11)

where ${\stackrel{\sim }{R}}_0$ is the galactocentric radius where the arm crosses the X-axis (i.e. ϕ = 0), and $\stackrel{\sim }{\alpha }$ is the pitch angle. The relative density in this model is given by the disk plane projected distance to the ridge according to

$$\begin{array}{c}\hfill \text{residual}(R,\varphi )=\stackrel{\sim }{a}+\stackrel{\sim }{b}exp\{-\frac{{[R-\stackrel{\sim }{R}(\varphi )]}^2{cos}^2(\stackrel{\sim }{\alpha })}{2{\stackrel{\sim }{L}}^2}\},\end{array}$$(12)

where $\stackrel{\sim }{a}$ and $\stackrel{\sim }{b}$ are the residual base-line and spiral arm amplitudes, and $\stackrel{\sim }{L}$ is the spiral arm width.

The residuals have some strong outliers; for this reason, we used the 1-norm in our fit (minimising |data − model|, instead of |data − model|²), which also gives somewhat more consistent results between data samples. Additionally, we performed a joint fit of all four data samples, which is probably more robust and less prone to over-fitting. The fitted models are shown in Fig. 7, and the parameter values are listed in Table 1. The relative amplitudes of the arm over-density fall in the range 18–24%, while the joint fit gives 20%. The fitted pitch angle is roughly 50°–60°, which is large compared to typical values in the literature (Levine et al. 2006; Vallée 2005, 2017). However, our value is likely only applicable to the very local area of study (e.g. similar to the high pitch angle segment of the Sagittarius Arm analysed by Kuhn et al. 2021). The arm width of the joint fit is 0.4 kpc.

Fig. 7. Simple spiral arm model, fitted to the residual of Φ(400 pc)−Φ(0 pc), for non-masked area cells within a distance of 2 kpc.

4. Discussion

Using vertical Jeans analysis and state-of-the-art data processing, we have produced 2D maps of the thin disk surface density within a distance of 3 kpc of the Sun. Overall, the results of our four data samples agree well. We see the same general over-dense feature in all data samples, and they also agree well with spiral arm tracers based on stellar age and metallicity.

We fitted a simple axisymmetric model to our results, as defined in Eq. (9). The inferred values for normalisation constant (parameter A) agree well with previous results for the immediate solar neighbourhood (Schutz et al. 2018; Widmark et al. 2021a,c). The inferred thin disk scale lengths (parameter h_L) fall in the range 3.3–4.2 kpc. This quantity varies in previous studies: a review from 2016 by Bland-Hawthorn & Gerhard (2016) summarises 15 articles to give 2.6 ± 0.5 kpc; later analyses include 2.2 ± 0.1 kpc (Widmark et al. 2022a), 2.4 ± 0.1 kpc (Wang et al. 2022), roughly 3.9 kpc (Robin et al. 2022), and $2.{17}_{-0.08}^{+0.18}\mathrm{kpc}$ (Ibata et al. 2023). An explanation for these discrepancies could be variations in the studied tracer samples and spatial volumes; this could give rise to different results if axisymmetry is broken or if the disk scale length is not constant with respect to R. In our own results, the masked data cells at low galactocentric radii (R ≲ 7 kpc) seem to follow a steeper profile, likely better described by a broken exponential function. However, we refrain from making a more definite or quantitative statement about this, since our results are quite discrepant beyond 2 kpc and potentially affected by systematic biases, especially towards the Galactic centre. The data sample with 1 < M_G ≤ 2 seems particularly discrepant at low R; this data sample has the lowest number count, which likely makes it more prone to systematic errors (how such errors operate is very complicated considering the underlying data processing and how it is affected by low statistics, e.g. the spectro-photometric distances affected by dust extinction, which is especially severe at low R).

We have detected the Local Arm in our thin disk surface density measurements, as revealed by the residual maps in Fig. 5, and we can compare and distinguish between other recent probes of the Local Arm morphology. The Local Arm model by Reid et al. (2019), based on masers as a tracer, is somewhat discrepant from recent observations of young star tracers (see Fig. 5 in Poggio et al. 2021 for a direct comparison with Gaia Collaboration 2023a,b). In the Reid et al. (2019) model, the Local Arm barely extends into negative X; instead, the outer Perseus Arm has a low pitch angle and curves through (X, Y)≃(−2, −2) kpc. In the recent observations, the Local Arm has a higher pitch angle, at least in the local area at positive Y, and also extends into negative X and Y, passing through (X, Y)≃(−2, −2) kpc. Our own results favour the Local Arm morphology based on the recent observations, which has a higher pitch angle. However, we do not see evidence for a Local Arm over-density that continues all the way to (X, Y)≃(−2, −2) kpc, as traced by Poggio et al. (2021), although we repeat that our results could be less accurate in spatial volumes beyond roughly 2 kpc.

There are some discrepancies in the inferred surface density of the four data samples, best seen in Fig. 5, in addition to the general arm-like structure. This noise is larger than can be accounted for by statistical uncertainties, which is indicative of some uncontrolled systematic errors. For example, there are small outlier ares, for example close to (X, Y)≃(−2, −0.4) kpc for the brightest data sample, and the arm-like over-density is closer to the solar position for the dimmest data sample. Systematic errors could arise in the data or in the data processing itself; for example, there is a very complicated interplay with the dust extinction field and the spectro-astrometric distance estimates from StarHorse.

There could also be more physical reasons that could bias our results, having to do with the dynamical properties of the respective data samples. For example, the dimmest data sample has a greater spatial scale height and is kinematically hotter; around the solar position, its u and v velocity dispersions are 10–25% higher, and its vertical velocity dispersion is 20–50% higher, as compared to the other three data samples. These properties make the data sample sensitive to a larger area of the disk plane’s mass density distribution, and also affects its sensitivity to biases that can arise from a broken equilibrium assumption; for example, the response of a stellar population with respect to a disk vibrational mode can depend on the scale height and kinematic temperature of that stellar population. In fact, Widmark et al. (2022b) found evidence of a breathing mode associated with the local spiral arm (see their Figs. 9 and 10, for stars with |z|< 300 pc), seen in both the stellar number density and vertical velocity distributions; they found a contracting breathing mode (net motion towards the mid-plane) in the inner part of the arm, offset by roughly a quarter wavelength with respect to the arm’s stellar over-density, which is indicative a breathing mode that is travelling against the direction of Galactic rotation. In Appendix A, we quantify the possible bias that can arise from such time-varying phase-space structures. In particular, we investigate vertical breathing modes and the phase-space spiral. We find that the former is most significant, but that its relative bias to the inferred value of Φ(400 pc)−Φ(0 pc) is limited to a few percent. For the latter, its effect is smaller than one percent. In summary, the biases related to such time-varying structures could be sizeable, but far from large enough to explain the arm-like over-density that we observe in this work.

Despite the systematic errors discussed above, our general results are robust. A scenario where the detected Local Arm over-density is caused by data systematics seems extremely contrived, since it is observed over a range of viewing angles, distances, and across four independent data samples. Furthermore, the data samples are quite different in age and small-scale structure; especially the stars of the dimmest data sample are older, kinematically hotter and more smoothly distributed, and do not trace the spiral arms (at least not in a significant manner). The detection of the spiral arm imprint in these data samples is thus not because of the stellar number counts themselves, but because of the gravitational influence that the spiral arms have on the stellar motions.

5. Conclusion

We applied the vertical Jeans equation to the Milky Way disk, dividing the disk plane into a square grid with a grid spacing of 100 pc, in order to study the thin disk surface density and its spatial dependence along the Galactic plane. We reached distances of a few kilo-parsecs, enabled by state-of-the-art data processing, utilising spectro-photometric distance measurements, GP regression on the stellar number density field, and radial velocity predictions from BNNs. Specifically, we measured the gravitational potential difference between the disk mid-plane and a height of 400 pc, which serves as a close proxy for the thin disk surface density.

For the first time, we see dynamically measured evidence for an over-density associated with a nearby spiral arm. This over-density is detected in all four of our data samples (defined by cuts in absolute magnitude), with a relative amplitude of roughly 20%. Our results agree well with the density of the Local Arm as traced by stellar metallicity and age, especially the density measured by Poggio et al. (2021). This constrains models of spiral arms and paves the way for future dynamical mass measurements of the Milky Way spiral structure.

Acknowledgments

We thank Yassin-Rany Khalil and Giacomo Monari for useful discussions and input. APN is supported by an Early Career Fellowship from the Leverhulme Trust. AW is sponsored by the Swedish Research Council under contract 2022-04283. 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). For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

References

Appendix A: Biases from time-varying structures

It is well established that there are time-varying structures in the disk, such as the recently discovered phase-space spiral (Antoja et al. 2018). When applying the vertical Jeans equation under the assumption of equilibrium, as we do in this work, such time-varying structures can bias the inferred gravitational potential. In this section we discuss and place upper limits to the relative amplitude of such systematic biases.

In order to investigate the possible bias that arises in the presence of a breathing mode or phase-space spiral, we first assumed an underlying model for the gravitational potential and stellar tracer population, which we can then perturb. For the gravitational potential, we used the model of baryonic densities from Schutz et al. (2018), building on the model from McKee et al. (2015), which also agrees well with the results of this work and other recent dynamical mass measurements of the immediate solar neighbourhood (Buch et al. 2019; de Salas & Widmark 2021; Widmark et al. 2021c). We used a stellar tracer population that is built from a superposition of three iso-thermal components according to

$$\begin{array}{c}\hfill f(z,w)=a_i\times \frac{exp(-\frac{\mathrm{\Phi }(z)+w^2/2}{{\sigma }_{w,i}^2})}{\sqrt{2\pi {\sigma }_{w,i}^2}},\end{array}$$(A.1)

where a_i = {0.4, 0.4, 0.2} are the mid-plane amplitudes of the three components, and σ_w,i = {6,14,20} km s⁻¹ are their respective velocity dispersions. The phase-space distribution at all heights is actually fully described by a function that only depends on vertical energy: f(E_z), where E_z = Φ(z)+w²/2, the nominator in the exponential above (see Widmark 2019 for more details on this analytic description).

The phase-space distribution is shown in Fig. A.1. It is consistent with the general shape of the stellar tracer populations we used in this work, both in terms of the vertical velocity variance and stellar number density profile. The precise shape of phase-space distribution of this underlying model is not very important in terms of the arising bias. We tried a range of different tracer populations, with varying scale heights and vertical velocity variance profiles, but saw similar relative biases in the inferred value of Φ(400 pc)−Φ(0 pc).

Fig. A.1. Phase-space distribution of our non-perturbed, steady-state model. The left panel shows a 2D histogram in the (z, w) phase-space plane. The middle panel shows the stellar number density profile, normalised to unity in the mid-plane. The right panel shows the standard deviation of vertical velocity, as a function of height.

A.1. Breathing mode

First we considered a breathing mode. We applied a perturbation to the phase-space distribution as a function of vertical energy (E_z) and vertical action angle (θ), according to (Mathur 1990; Weinberg 1991; Widrow & Bonner 2015)

$$\begin{array}{c}\hfill f(E_z,\theta )=f_0(E_z)\{1+ϵE_{z}cos[2(\theta -{\theta }_0]\},\end{array}$$(A.2)

where f₀(E_z) is the unperturbed distribution, and θ₀ and ϵ are constants.

We set ϵ to a value such that Δw_{400 pc} (the difference in mean velocity 400 pc above and below the mid-plane) can be at most 1 km s⁻¹, which happens for θ₀ = π/4. In Gaia Collaboration (2023a), they compare the mean vertical velocity above and below the mid-plane, stating that |Δw| reaches values of 2–4 km s⁻¹. However, such large values are seen only in more distant regions of the disk (beyond 3 kpc; see their figure 23). For the more nearby solar neighbourhood, similar to the volume we are probing, |Δw| is limited to less than 1 km s⁻¹. Widmark et al. (2022b) did a similar analysis of the four specific data samples studied in this work, and with a more refined binning in height (see their Figs. 8, 9, and 10 and their appendix). They evaluated Δw specifically for the range |z|∈[300, 500] pc, and saw |Δw| typically below 1 km s⁻¹, although there are some smaller regions with absolute values of roughly 2 km s⁻¹, especially around (X, Y) = (0, 1.5) kpc.

The relative bias that arises from the breathing mode with respect to the non-perturbed phase-space distribution, as a function of the angle θ₀, is shown in Fig. A.2. As illustrated in the top panels, the phase-space distribution is the most spatially compressed at θ₀ = 0, which is where the bias is positive (i.e. the inferred potential is over-estimated). Conversely, it is the most spatially stretched out for θ₀ = π/2, where the bias is negative. This breathing mode is fully uniform, in the sense that particles of the same action angle are affected in the same manner. In the real Galactic disk, such an idealised mode would phase-mix. This would likely lower the dynamical mass measurement bias over time, even if the amplitude of Δw_{400 pc} is sustained.

Fig. A.2. Relative bias to the inferred value of Φ(400 pc)−Φ(0 pc) using Jeans analysis under the assumption of equilibrium, in the presence of a breathing mode as a function of θ0. The bias is shown in the bottom panel, where the breathing mode is normalised such that the difference in mean value of w between 400 pc above and below the mid-plane is Δw400 pc = 1 km s−1 when θ0 = π/4. The top panels show how the breathing mode phase, θ0, affects the phase-space distribution. The perturbation is strongly enhanced in the top panels, for better visibility.

In summary, breathing modes in the stellar disk can likely bias the inferred value of Φ(400 pc)−Φ(0 pc) by a few percent. The arm-like structure we observe in this work has a relative over-density of roughly 20%. To fully explain this as a bias arising due to a breathing mode would require amplitudes corresponding to Δw_400 pc ≳ 6 km s⁻¹, which we can confidently rule out.

A.2. Phase-space spiral

We did a similar test for a phase-space spiral. We used the model from Widmark et al. (2021b), with a relative over-density amplitude of 20%, and a time of 500 Myr since the perturbation was produced. These parameter values give rise to a phase-space spiral that is qualitatively similar to the one observed in the actual solar neighbourhood. We tried different parameter choices, giving rise to, for example, different degrees of spiral winding, and saw similar results in terms of the strength of the possible bias.

The relative bias with respect to the non-perturbed phase-space distribution is shown as a function of θ₀ in Fig. A.3. Here, θ₀ is in the range [0, 2π] since it is an asymmetric structure (as opposed to the symmetric breathing mode of Fig. A.2). Here, the bias to Φ(400 pc)−Φ(0 pc) is limited to less than a percent.

Fig. A.3. Relative bias for the inferred value of Φ(400 pc)−Φ(0 pc) using Jeans analysis under the assumption of equilibrium and in the presence of a phase-space spiral as a function of θ0. The top panels show the relative perturbation of the spiral (unlike Fig. A.2, which showed the total number densities), where red (blue) correspond to a over-density (under-density), whose relative amplitude is in the range [ − 20, 20]%.

All Tables

All Figures

Fig. 1. Stellar number density distribution for a few area cells of the data sample with 2 < MG ≤ 3. Each panel is labelled by the area cell’s mid-point (X, Y) coordinates. The function obtained by GP regression is shown in solid black; the fitted symmetric function ($\widehat{n}$) of Eq. (5) is shown in dashed red; the raw number count (N) with its associated Poisson noise uncertainty is shown as a 1σ band.

In the text

	Fig. 2. Fraction of stars with radial velocity measurements across the 100 × 100 pc area cells used in our analysis. The four panels represent the four absolute magnitude samples, as labelled.
In the text

Fig. 3. Examples of vertical velocity dispersions and the best-fitting models (Eq. (6)). The four rows correspond to the four magnitude samples, while the three panels within each row are randomly chosen cells in the (X, Y) plane. The yellow points represent the measured velocity dispersions in adaptively spaced Z bins, while the solid lines show the best-fitting models. Note that the model is not fit to these binned dispersions, but to the individual stellar velocities. In each case, 20 sets of data and 20 models are shown, corresponding to the 20 random imputations of the dataset as described in the main text.

In the text

Fig. 4. Gravitational potential, Φ(400 pc)−Φ(0 pc), as a function of galactocentric radius, for the four stellar samples (as labelled in the panels’ top-right corner). The dashed line is the exponential function of Eq. (9) fitted to data points with $\sqrt{X^2+Y^2}<2\mathrm{kpc}$, which excludes distant and dubious spatial regions. The scatter points correspond to area cells in the (X, Y) plane. Area cells that are included in the exponential function fit are coloured according to their galactocentric disk plane angle (ϕ); excluded data cells have light grey and smaller scatter points. The colour bar and axis ranges are shared between all panels.

In the text

Fig. 5. Residual of Φ(400 pc)−Φ(0 pc) with respect to a fitted exponential, projected on the (X, Y) plane. Each panel corresponds to one of our four data samples. In each panel, dotted grey lines show galactocentric iso-radial and iso-azimuthal lines, and the black plus sign in the centre marks the solar position. The grey contour lines correspond to the spiral arms as mapped by Poggio et al. (2021), as seen in Fig. 6.

In the text

	Fig. 6. Over-density of upper main sequence stars in the (X, Y) plane as a tracer of spiral arms. Results are taken directly from Poggio et al. (2021). The grey-scale contour lines are the same as in Fig. 5, corresponding to values ( − 0.5,  0,  0.5).
In the text

	Fig. 7. Simple spiral arm model, fitted to the residual of Φ(400 pc)−Φ(0 pc), for non-masked area cells within a distance of 2 kpc.
In the text

	Fig. A.1. Phase-space distribution of our non-perturbed, steady-state model. The left panel shows a 2D histogram in the (z, w) phase-space plane. The middle panel shows the stellar number density profile, normalised to unity in the mid-plane. The right panel shows the standard deviation of vertical velocity, as a function of height.
In the text

Fig. A.2. Relative bias to the inferred value of Φ(400 pc)−Φ(0 pc) using Jeans analysis under the assumption of equilibrium, in the presence of a breathing mode as a function of θ0. The bias is shown in the bottom panel, where the breathing mode is normalised such that the difference in mean value of w between 400 pc above and below the mid-plane is Δw400 pc = 1 km s−1 when θ0 = π/4. The top panels show how the breathing mode phase, θ0, affects the phase-space distribution. The perturbation is strongly enhanced in the top panels, for better visibility.

In the text

Fig. A.3. Relative bias for the inferred value of Φ(400 pc)−Φ(0 pc) using Jeans analysis under the assumption of equilibrium and in the presence of a phase-space spiral as a function of θ0. The top panels show the relative perturbation of the spiral (unlike Fig. A.2, which showed the total number densities), where red (blue) correspond to a over-density (under-density), whose relative amplitude is in the range [ − 20, 20]%.

In the text