Orbital Torus Imaging: Acceleration, Density, and Dark Matter in the Galactic Disk Measured with Element Abundance Gradients

Under the assumption of a simple and time-invariant gravitational potential, many Galactic dynamics techniques infer the milky Way’s mass and dark matter distributions from stellar kinematic observations. These methods typically rely on parameterized potential models of the Galaxy and must take into account nontrivial survey selection effects, because they make use of the density of stars in phase space. Large-scale spectroscopic surveys now supply information beyond kinematics in the form of precise stellar label measurements (especially element abundances). These element abundances are known to correlate with orbital actions or other dynamical invariants. Here, we use the Orbital Torus Imaging framework that uses abundance gradients in phase space to map orbits. In many cases these gradients can be measured without detailed knowledge of the selection function. We use stellar surface abundances from the Apache Point Observatory Galactic Evolution Experiment survey combined with kinematic data from the Gaia mission. Our method reveals the vertical (z-direction) orbital structure in the Galaxy and enables empirical measurements of the vertical acceleration field and orbital frequencies in the disk. From these measurements, we infer the total surface mass density, Σ, and midplane volume density, ρ 0, as a function of Galactocentric radius and height. Around the Sun, we find Σ⊙(z=1.1kpc)=72−9+6M⊙pc−2 and ρ⊙(z=0)=0.081−0.009+0.015M⊙pc−3 using the most constraining abundance ratio, [Mg/Fe]. This corresponds to a dark matter contribution in surface density of Σ⊙,DM(z = 1.1 kpc) = 24 ± 4 M ⊙ pc−2, and in total volume mass density of ρ ⊙,DM(z = 0) = 0.011 ± 0.002 M ⊙ pc−3. Moreover, using these mass density values we estimate the scale length of the low-α disk to be h R = 2.24 ± 0.06 kpc.


INTRODUCTION
Galaxies are dynamic and evolving structures driven by the interactions and assembly of their constituent stars, gas, and dark matter (White & Frenk 1991).The resulting observed structure of galaxies provides important information about the properties and behavior of dark matter, and the processes involved in galaxy formation.Unfortunately, we only Corresponding author: Danny Horta dhortadarrington@flatironinstitute.org observe galaxies at (effectively) a single snapshot in time.
To turn this snapshot of stellar and gas kinematics into inferences about galaxy formation processes and dark matter, the study of galactic dynamics relies on theoretical frameworks and methods built on statistical physics (Binney & Tremaine 2008).Over the last century, this endeavor has led to important measurements of the structure of mass and dark matter in our galaxy, the Milky Way.
A key quantity that has been used historically to summarize the mass distribution in the Milky Way is the total surface mass density near the sun, Σ ⊙ (e.g., Oort 1960;Kuijken & Gilmore 1989, 1991;Creze et al. 1998;Holmberg & Flynn 2000).The total surface mass density is typically measured from stellar kinematics using methods that rely on strong assumptions about the stellar distribution function (DF) and form of the gravitational potential (e.g., Kuijken & Gilmore 1991;Flynn et al. 2006;McMillan 2011;Bovy & Rix 2013;Zhang et al. 2013).In all cases, measurements of the total surface mass density together with estimates of the surface density from stars and gas provides a means to infer the local contribution of dark matter to the mass distribution in the galaxy (e.g., McKee et al. 2015;Bland-Hawthorn & Gerhard 2016).
One prominent approach is "Jeans modeling" (e.g., Jeans 1919;Binney & Tremaine 2008, and references therein), which employs a set of equations (the Jeans equations) that relate spatial derivatives of second-order velocity moments to the gravitational potential.More general approaches instead attempt to model the DF explicitly (e.g., Kuijken & Gilmore 1991;Bovy & Rix 2013;Piffl et al. 2014;Widmark et al. 2022b;Li & Widrow 2021, 2023;Cheng et al. 2023).However, these models become computationally expensive for large data sets or are unable to incorporate flexible model forms.Moreover, these methods rely on an accurate accounting of the survey's selection function, making them challenging to apply precisely.
The advent of large-scale stellar surveys -especially the Gaia mission (Gaia Collaboration et al. 2016, 2022) -has revolutionised Galactic astrophysics.We are now able to measure the phase-space coordinates of over a billion stars throughout the Galaxy.Owing to large spectroscopic surveys, it is also now possible to obtain high-quality measurements of stellar parameters/labels (temperatures, surface gravities, metallicities, and other element abundances) for over a million stars in the Milky Way (e.g., APOGEE : Majewski et al. 2017;GALAH : Martell et al. 2017;LAMOST : Cui et al. 2012;amongst others).This unprecedented amount of highquality information is supplying the necessary means to further develop dynamical inference frameworks.It also provides a window of opportunity to study the relationship between stellar labels, that should be approximately invariant over long timescales, with the kinematics and orbits of stars in the Galaxy.
Recent results have suggested that combining spectroscopic and astrometric data for large samples can provide precise measurements of the Galactic mass distribution (e.g., Sanders & Binney 2015;Price-Whelan et al. 2021;Binney & Vasiliev 2023).One example of such techniques is Orbital Torus Imaging (OTI; Price-Whelan et al. 2021), which utilises chemical abundance information to directly map the orbit structure of stars in the Galactic disk.It uses the idea that, in a well-mixed (steady-state) population, abundance moments (like mean abundances) can depend only on ac-tions, or dynamical invariants, or integrals of the motion; they cannot depend on orbital phases.OTI was originally developed and used in the context of modeling the vertical kinematics of stars (i.e., position z and velocity v z perpendicular to the Galactic disk).While this framework is still new, it is a promising avenue toward developing models of our Galaxy that relaxes many assumptions that have been necessary in past modeling efforts.In particular, provided that the stars are selected in an abundance-neutral way, the method does not require much (or any) knowledge of the survey selection function.
In the first of a series of papers focusing on this topic, here we set out to combine stellar labels (namely, mean element abundances ratios) obtained from spectroscopic measurements using the Apache Point Observatory Galactic Evolution Experiment (APOGEE Majewski et al. 2017) survey, with kinematic observations from the latest Gaia data release 3 (Gaia Collaboration et al. 2022), in order to measure the surface mass density and vertical acceleration field across the disk of the Milky Way.We do this by constructing a model under the Orbital Torus Imaging framework, building on the method described in (Price-Whelan et al, in prep).
In Section 2 we describe the data and sample used.In Section 3 we show and discuss the abundance gradients in phase space.In Section 4 we describe the method.Section 5 contains our results, which are then discussed in Section 6.We finalise by providing our concluding thoughts in Section 7.

DATA
We use a cross-matched catalogue of the latest spectroscopic data from the APOGEE survey (DR17 Majewski et al. 2017;Abdurro'uf et al. 2022) and astrometric data from the Gaia mission data release 3 (Gaia DR3 Gaia Collaboration et al. 2022).The APOGEE data are based on observations collected by two high-resolution, multi-fibre spectrographs (Wilson et al. 2019) attached to the 2.5m Sloan telescope at Apache Point Observatory (Gunn et al. 2006) and the du Pont 2.5 m telescope at Las Campanas Observatory (Bowen & Vaughan 1973), respectively.Element abundances are derived using the ASPCAP pipeline (García Pérez et al. 2016) based on the FERRE code (Allende Prieto et al. 2006) and the line lists from Cunha et al. (2017) and Smith et al. (2021).The spectra themselves were reduced by a customized pipeline (Nidever et al. 2015).For details on target selection criteria, see Zasowski et al. (2013) for APOGEE, Zasowski et al. (2017) for APOGEE -2, Beaton et al. (2021) for APOGEE north, and Santana et al. (2021) for APOGEE south.
The Gaia mission/survey (Gaia Collaboration et al. 2016) delivers detailed sky positions, proper motion, and parallax measurements for ∼ 2 billion stars, limited only by their apparent magnitudes (Gaia G ≲ 20.7).Here we use only astro- Together, the APOGEE and Gaia data supply full 6D phase space information from which kinematic and orbital properties can be derived.The APOGEE data additionally provide detailed chemical abundance measurements for up to ∼ 20 different elemental species spanning the α, odd-Z, iron-peak, and s-process nucleosynthetic channels.

Parent Sample
The parent sample used in this work is comprised of stars that satisfy the following selection criteria:

Red Giant Branch stars:
APOGEE -determined atmospheric parameters, effective temperature and surface gravity 3500 < T eff < 5500 K and 1.5 < log g < 3.5,  Note that we zero-out the z-component of the solar velocity, v z,⊙ , and set the solar height above the galactic midplane to z ⊙ = 0 as we measure these quantities below.Lastly, for our modelling we only use stars within |z| < 2 kpc from the solar position, and have vertical velocities with respect to the Sun smaller than |v z | < 80 km s −1 .This removes likely halo contaminants and/or stars not part of the Milky Way disk.We additionally select stars that have small eccentricities (i.e. on near-circular orbits) by enforcing |v R | < 50 km s −1 and |∆R| < 1 kpc, where R is the present-day cylindrical radius of each star, and ∆R = R − R g , where R g = L z /v circ is the approximate guidingcenter radius of each star computed using the rotation curve from Eilers et al. (2019).We chose a cut of |v R | < 50 km s −1 as after inspecting the ∆R-v R plane, where we find that the approximate maximum radial velocity a star with |∆R| < 1 kpc could reach is v R ∼ 50 km s −1 .Our final working sample is comprised of 94,685 stars.

ELEMENT ABUNDANCE GRADIENTS IN VERTICAL KINEMATICS
The stars in the Milky Way disk have a well-known dependence on position, velocity, and chemical abundance within the Galaxy.For example, there is a known interrelation between Galactic position, iron abundance [Fe/H], and age for stars in the Galactic disk (e.g., Frankel et al. 2018).The Milky Way disk has also been shown to be homogeneous in terms of its star formation environment for most elements (e.g., Ness et al. 2019Ness et al. , 2022;;Horta et al. 2022), and to show a continuity in mono-abundance populations (e.g., Bovy et al. 2016;Mackereth et al. 2017).These element abundance gradients also exist with respect to the vertical kinematics of stars: the mean abundance (in any abundance ratio) of stars in the Galactic disk typically depends on both the vertical height z and vertical velocity v z of the stars (e.g., Price-Whelan et al. 2021).We leverage these gradients by developing a model in the Orbital Torus Imaging framework to fit mono-abundance contours in the vertical phase space (z, v z ) for different element abundance ratios.
In this section, we first visualize the available abundance data as a function of vertical phase space coordinates (Section 3.1) and show how these abundance gradients contain information about the orbital structure of the Galactic disk (Section 3.2).We then show that some element abundance ratios have a stronger dependence on orbital properties, which should make them more informative for our analysis (Section 3.3).Finally, we explore the dependence of the vertical abundance distribution as a function of spatial location in the disk for our most promising abundance ratio, [Mg/Fe] (Section 3.4).

Exploring different element abundances
Figure 2 shows mean element abundances, ⟨[X/Y]⟩, binned by vertical phase space coordinates, (z, v z ), for a subset of our parent sample positioned around the Sun that have low eccentricities.Each pixel has a size of δz = 0.04 kpc and δv z = 1.76 km s −1 .Specifically, we show six element abundance ratios provided by the APOGEE survey that represent the α elements (O, Mg, Si), odd-Z elements (Al), and iron peak elements (Fe and Ni).We provide equivalent plots for all elements measured by APOGEE in Appendix A. The six elements shown in Figure 2 are, in theory, the most reliable in the APOGEE survey (Jönsson et al. 2020).With these six elements, we are sensitive to differences in contributions from different supernova (SN) types (SN II and SN Ia) for stars with different orbital properties.
Upon inspection of ⟨[X/Y]⟩ in (z, v z ), it is clear is that there are gradients with respect to both height z and vertical velocity v z .For elements synthesised predominantly in the explosions of SN II (like α and odd-Z), the gradient in ⟨[X/Y]⟩ is positive with increasing z, v z .Conversely, for heavier elements like [Fe/H], synthesized in both SN II and SN Ia, the gradient in ⟨[X/Y]⟩ is in the opposite direction, growing in magnitude with decreasing z, v z .We also note that the range, [min ⟨[X/Y]⟩ ,max ⟨[X/Y]⟩ ], in the gradient for each element abundance is different.We find that [Fe/H] has the largest range, spanning from ∼ −0.3 − 0.1 dex, and [Ni/Fe] has the smallest range, spanning just ≈ 0.06 dexthe fact that such a small gradient is detectable speaks to the high quality element abundance measurements produced by the APOGEE survey.
Another interesting property that is apparent in the element abundance gradients in z,v z space is the shape of the monoabundance contours.Figure 2 reveals that both the shape and element abundance gradient are more pronounced/steeper for element abundances synthesised via SN II (for example, Mg and Si) as compared with other elements (e.g., Ni).These small yet noticeable differences are likely caused by different elements being synthesised at different rates, leading to stellar populations manifesting different mono-abundance contour shapes in phase space.Despite such differences, the fact that we observe mono-abundance z, v z contours for all elements is proof that this chemical-dynamical data is rich in information, and in principle should retain vital clues about the structure of the Galactic disk.

Mapping orbits with mono-abundance contours
As seen in Figure 2, there are clear gradients of ⟨[X/Y]⟩ in z, v z space.Stars at small heights z and small vertical velocities v z have (on average) different mean element abundance ratios than stars at larger z or v z .In fact, these ⟨[X/Y]⟩ gradients are related with respect to phase space position in z, v z .Stars that have higher vertical velocities will reach higher maximum heights above the disk, z max , as they orbit, and vice versa.As stars orbit in this phase plane, their trajectories roughly follow elliptical shapes.Those elliptical trajectories are the shapes we see in Figure 2 for pixels of constant ⟨[X/Y]⟩.
Figure 3 illustrates that orbital trajectories in z, v z space are identifiable using curves of ⟨[X/Y]⟩ gradients seen in Figure 2. The left panel shows eight orbital trajectories in z, v z with different values of the vertical action, J z , obtained from using a toy model of the Galactic potential from Price-Whelan et al., (in prep).We note that this is our only usage of a parameterized Galactic mass model in this work, and this is only for demonstration purposes.The orbits of stars with lower J z values trace out nearly elliptical trajectories that reach smaller heights above the Galactic plane, z max ,

Identifying informative element abundance gradients
The shapes of contours of constant ⟨[X/Y]⟩ in Figure 2 nearly correspond to orbits in the vertical phase space (Fig-ure 3).However, our ability to identify contours of constant ⟨[X/Y]⟩ will depend on the steepness and scatter of the gradients of ⟨[X/Y]⟩ vs. vertical action J z .Thus, in order to assess which element abundance ratio we expect to give the best constraining power, we investigate the trends between our set of well-measured APOGEE abundances as a function of a proxy for the maximum z height a star reaches, ∼ z max .To compute z max , which is related to the vertical action J z , one has to assume a form for the Galactic mass distribution, which we are trying to avoid.We instead here sub-select our data to have |v z | < 10 km s −1 so that those stars' instantaneous vertical positions, z, can serve as a proxy for z maxi.e.z |vz|<10 km s −1 ∼ z max .Figure 4 shows the resulting de- Right: We paint element abundances onto the simulated stars with a linear dependence on vertical action Jz with a slope similar to that of [Mg/Fe] (see Figure 2).The filled contours then show curves of constant element abundance ratio for a simulated population of orbits in the vertical phase space.See Price-Whelan et al. (in prep) for more details.
pendence of each element abundance ratio as a function of our potential-independent proxy for the maximum z height, ∼ z max .The pixels show the column-normalized densities of stars in this space for stars with orbits similar to the Sun, R g −R g,⊙ < 0.5 kpc, where R g is the guiding-centre radius.All of our selected element abundance ratios show a trend in their mean value (red dashed line) below z ≲ 1 kpc.Above this height, z ≳ 1 kpc, the clearest trends are shown by the elements in the α group, with [Mg/Fe] being the clearest.
To quantify the power an element abundance ratio has for inferring the orbital structure of the Milky Way disk, in Figure 5 we show the signal-to-noise value for every element abundance ratio supplied by APOGEE.This value is computed by taking the ratio between the slope of ⟨[X/Y]⟩ with ∼ z max with the average variance around the median line, σ 2 ⟨[X/Y]⟩ .Figure 5 confirms our qualitative analysis from Figure 4, and shows that [Mg/Fe] is the element abundance ratio with the most constraining power.Given this test, we will use [Mg/Fe] as our representative element abundance ratio for our model fitting (Section 4).

Changes in vertical mono-abundance contours across the Galactic disk
Figure 2 shows element abundance trends with vertical kinematics for stars around the Sun.However, our sample contains stars with a wide range of orbits in the Galactic disk, owing to the exquisite volume and precision of data from Gaia.In Figure 6, we therefore examine how the ⟨[X/Y]⟩ contours in vertical kinematics for [Mg/Fe] vary as a function of orbital location across the Galactic disk.We illustrate this using the same visualization of ⟨[Mg/Fe]⟩ in z, v z , but now divide our sample into bins of similar guiding-center radius, R g , from top left to bottom right going from orbits in the inner to outer Galactic disk.Only stars with 5 < R g < 15 kpc are shown, binned in 2 kpc-wide bins for illustration purposes.
Figure 6 shows that the morphology of constant ⟨[Mg/Fe]⟩ contours varies as a function of guiding-center radius.Stars in the inner Galaxy (top left) show a more oblate monoabundance contour than stars in the outer Galaxy (bottom right).This "squashing" is around z for different v z , and becomes less pronounced as one increases in R g , turning from an oblate shape to a more circular and then prolate shape.As we will see in Section 5, the change in mono-abundance contour shapes is related to the surface mass density in the disk.Towards the inner disk, where the surface mass density within a given column is higher, stars in the Galactic disk feel a stronger gravitational influence as they oscillate above and below the Galactic midplane.This leads to stars of a given element abundance ratio not being able to reach larger vertical heights for a given vertical velocity value, which leads to flattened mono-abundance contours in z, v z space.Conversely, towards the outer regions of the Galaxy, where the surface mass density in the Galactic disk is lower, and the influence from the gravitational pull of the dark halo is stronger, stars are able to reach higher vertical heights given a particular vertical velocity magnitude.This leads to a more prolate and then diamond-shaped profile for stars with a given element abundance ratio.
The fact that we observe such changes in the morphology of constant ⟨[Mg/Fe]⟩ in z, v z highlights that the disk is rich in chemical-dynamical structure that is useful for inferring properties of the Galaxy.From the contours of constant ⟨[Mg/Fe]⟩, we can already see, conceptually, how these data can be used to infer mass density: The aspect ratio of an ellipse in these panels has units of frequency (i.e.v z /z), and frequency is related to density through Poisson's equation, (1) In the following section, we describe a method for empirically modelling the mono-abundance z, v z contours, following the method from Price-Whelan et al. (in prep.).

ORBITAL TORUS IMAGING: MODELING THE ORBITAL STRUCTURE AND MASS OF THE DISK WITH ELEMENT ABUNDANCES
Orbital Torus Imaging (OTI) is a framework for exploiting the idea that, in a steady-state galaxy, the distribution Trends of [X/Fe] abundance with vertical orbit amplitude, ∼ z max  , for all element abundance ratios supplied by APOGEE.The higher the signal-to-noise ratio, the more constraining power an element abundance ratio has for modelling the Milky Way disc in z, vz space.Quantitatively, [Mg/Fe] is the element abundance ratio with the highest constraining power.
of intrinsic stellar invariants of a tracer population can only depend on dynamical invariants and not on time or orbital phase.In a 2D slice of phase-space, for example the vertical kinematics w = (z, v z ) for a small region of the Galactic disk, an implication of this is that contours of constant mean element abundance should correspond to orbital trajectories.With orbital trajectories, we can measure properties of the acceleration field and mass distirbution of the Galaxy.A full description of our methodology is described in a companion paper (Price-Whelan et al., in prep.), but here we briefly summarize the framework we use to infer the disk properties from the combined vertical kinematics and element abundance data.Our model for fitting ⟨[X/Y]⟩ works only in z, v z , but we will explore how the properties of the disk vary as a function of radius by binning the data.Our framework is summarized in the following steps: 1. We bin the the parent sample data into 25 equally-sized (2 kpc-wide) but overlapping bins of guiding-center radius, R g , ranging from 5 < R g < 13 kpc.For each bin of R g , we compute ⟨[X/Y]⟩ and the median absolute deviation (MAD, defined as median(|X i − X|), where X is ⟨[X/Y]⟩), for each of our element abundance ratios in small bins (δ z = 0.04 kpc, δ vz = 1.76 km s −1 ) of z, v z .
2. We fit a flexible, empirical model of the element abundance trends to the pixelized using the method described in Price-Whelan et al. (in prep.)2 .We use the empirical scatter of the abundances in each pixel (the MAD values) scaled by 1/ √ N , where N is the number of stars in a pixel, as our uncertainty when fitting this model.
3. We optimize the following likelihood function, where for each pixel in w = (z, v z ) space.Here, ξ = 1.5 × MAD(w) data / N stars (w), and θ is the set of model parameters we aim to optimize (described in detail at the end of this Section).
4. We use the best-fit model and interpret the level sets of the fitted function as orbits, allowing us to estimate the vertical acceleration a z (z), surface mass density Σ(z), and midplane volume density ρ ⊙ (z = 0 kpc) for each bin of guiding radius and for each element abundance ratio.To compute the uncertainties on these measured quantities, we bootstrap sample with replacement, instead of using the standard error estimate.
In the following, we briefly outline the Orbital Torus Imaging approach and assumptions that go into this method and discuss how we map orbits from element abundance gradients (Section 4.1), before presenting the method for determining the vertical acceleration and surface mass density from our model (Section 4.2).

Modeling approach and underlying assumptions
The core idea of Orbital Torus Imaging (OTI) is that any stellar label (e.g., element abundance ratios) that correlates intrinsically with the orbital properties of stars in the Galaxy can be used in place of dynamical quantities to "image" the orbits.The element abundance distribution should be constant along an orbit, and stars at different phases along an orbit should have the same element abundance distribution.Stellar labels can therefore be used to trace the orbits by finding curves of constant label value (see Price-Whelan et al. 2021 for more details).Once the shapes of orbits are traced using this method, we can use standard dynamical arguments to extract quantities of interest like the acceleration, surface mass density, and midplane mass density.The advantage of this approach is that our model is very flexible: we never parameterize the global form of the gravitational potential and instead only locally (in guiding radius) parameterize the shapes of orbital trajectories in the vertical phase space.Our model form and methodology is described in detail in Price-Whelan et al. (in prep), which is based on ideas presented previously (Price-Whelan et al. 2021).
Our current implementation of OTI relies on a set of assumptions described below: Axisymmetric and separable: The gravitational potential that the stars orbit in is axisymmetric, smooth, and separable in cylindrical radius R and vertical position z: . This allows us to work with just one position and velocity component (namely, z and v z ), but this could equally be applied to radial kinematics R and v R .In future work, we will explore models that simultaneously model the radial and vertical phase space.
Near circular orbits: For the previous assumption to be a valid simplification, we also require that stars have negligible eccentricity (or close to zero radial action J R = 0) and have the same z-component of the angular momentum L z (or azimuthal action J ϕ ).In order to get close to this, we have restricted our disk sample to |R − R g | < 1 kpc in any bin of guiding radius R g , and have ensured the stars have low radial velocities (or radial action), |v R | < 50 km s −1 .
Phase-mixed: The stellar distribution function in vertical kinematics, f (z, v z ), in any bin of guiding radius is phase mixed.This assumption is the key assumption; it is what ensures that contours of constant element abundance will correspond to curves of constant orbital actions (or other invariants).
Negligible kinematic measurement uncertainties: When fitting our model below, we always bin our stellar data into small bins of vertical phase-space coordinates (δz, δv z ).We assume that most measurements of vertical position z and velocity v z for stars in our sample are more precise than our adopted bin sizes, so we ignore measurement uncertainties on the kinematic quantities.
For a slice of phase space at fixed values of the other actions (namely, J R and J ϕ ), under our assumptions, contours of constant phase space ⟨[X/Y]⟩ delineate orbital trajectories (Price-Whelan et al. 2021).This concept is powerful and motivates a path towards empirically mapping the orbital structure in vertical kinematics directly from the observed element abundance ratios in a slice of phase space.
Our method works by modeling the shapes of level sets of ⟨[X/Y]⟩ in z, v z space.We do this by parameterizing the shapes of the contours as a low-order Fourier distortion away from ellipses of constant axis ratio (i.e.constant frequency).Very near the Galactic midplane, where the density distribution is approximately uniform, we expect the orbits (and therefore contours of constant ⟨[X/Y]⟩) to be very nearly elliptical.For orbits that reach larger heights above the plane, we expect their frequencies to be smaller (i.e. they have longer oscillation periods).For orbits that reach many scale heights above the midplane (i.e.z max ≳ 1 kpc), the influence of the Milky Way's dark matter halo will further distort the shapes of the orbits to be more diamond-shaped and less elliptical.
We therefore parameterize the shapes of these level sets using a mixture of m = 2 and m = 4 Fourier distortions away from a constant axis ratio ellipse.That is, we define an elliptical (euclidian) distance rz and angle θz in the vertical phase space, where and z 0 , v z,0 are parameters that set the solar position and velocity with respect to the midplane, and Ω 0 is a parameter that sets the frequency of an orbit with zero rz .We parameterize the level sets of the element abundances as curves of constant distorted distance r z , defined as In the above expression, the functions e 2 (r z ) and e 4 (r z ) are yet-to-be specified amplitude functions that set the oblateness of the contours and diamond-shaped distortion due to the dark matter halo.We require that these functions go to zero as rz → 0 so that Ω 0 can be interpreted as the asymptotic midplane orbital frequency.As we expect the orbits to get more prolate with increasing J z , we require that the m = 2 coefficient function (e 2 ) is monotonic and strictly positive.
We also expect the orbits to become more "diamond-shaped" at large J z and not more square, so we require the m = 4 coefficient function to be monotonic and negative.To enforce these things within a fitting procedure, we represent the functions e 2 (r z ) and e 4 (r z ) as monotonic quadratic splines.We fix the knot locations of the splines to be equally spaced in r2 z , and treat the function derivative values at the fixed knot locations as free parameters.
With the above, we have a formalism for labeling (and foliating) curves in the vertical phase space using the distorted elliptical distance r z (Equation 8).To fit ⟨[X/Y]⟩ as a function of z, v z , we then must also model the dependence of ⟨[X/Y]⟩ as a function of this distorted radius, i.e. ⟨[X/Y]⟩(r z ).We represent this function also using a monotonic quadratic spline.The knot locations of the splines are also set to be equally spaced in r2 z .Our model is a generative model of ⟨[X/Y]⟩ at any location in z, v z .We fit the model parameters by maximizing the log-likelihood of the data (Equation 2) -the value of the ⟨[X/Y]⟩ at locations of z, v z space -by treating the MAD value in each pixel of z, v z as our uncertainty, and ignoring any uncertainties in the phase space coordinates.Our implementation is built on the JAX framework (Bradbury et al. 2018) such that we can easily auto-differentiate through our model, allowing us to efficiently use a gradient-based optimizer (BFGS; Fletcher 1987) to fit our model parameters.For a given element abundance ratio, for each bin of guiding radius, our model has eighteen free parameters: the midplane frequency Ω 0 , the solar position and velocity z 0 and v z,0 , the spline function derivative values (10) for the e 2 and e 4 coefficient functions, and the spline function values (5) for the abundance function ⟨[X/Y]⟩(r z ).

Measuring the vertical acceleration field and surface mass density
Any function of the phase-space coordinates, F , that is time-invariant and constant along an orbit (i.e. a constant of the motion), F (z(t), v z (t)) = const., provides a means to measure the acceleration: For a given setting of our parameters (see Section 4.1), we can compute the vertical acceleration inferred by our model using the chain rule applied to our orbit label function, the distorted elliptical radius r z , which, by construction, is a constant of the motion.As we are interested in evaluating the acceleration as a function of z and not v z , we take the limit of this expression as v z → 0 and θz → π 2 : This already reveals a subtlety in our parameterization and modeling method: in a gravitating system, the acceleration should not be allowed to be a function of the velocity v z , but our model permits this; we return to this point in the discussion (Section 6).If we substitute our expression for r z (Equation 8) into the above limit, we find that for m = {2, 4}, and the whole expression is evaluated at rz = |∆z| √ Ω 0 (see Price-Whelan et al., in prep., for a full derivation).We note that the factor −Ω 2 0 z is equivalent to the acceleration in a simple harmonic oscillator potential, but for our parameterization, this is made anharmonic through the Fourier coefficient functions e 2 and e 4 .
At this point, it is worth asking "How reliably can we empirically measure a z ?"In a companion paper, Price-Whelan et al., (in prep), we focus on the method and test the Orbital Torus Imaging framework described here in a number of different scenarios to assess the accuracy of the a z measurement.Specifically, Price-Whelan et al., (in prep) test: i) a simple harmonic oscillator; ii) a quasi-isothermal disk DF simulation; iii) a quasi-isothermal disk DF simulation with a selection function imprinted resembling that of the APOGEE survey; iv) a perturbed N-body disk simulation.In all cases, the empirical measurement of the a z matches well the true value for measurements below the range we explore in this work (|z| < 1.2 kpc, see Section 5).
Once the vertical acceleration is determined, we then compute the volume and surface mass density using Poisson's equation.We assume that the total mass distribution is axisymmetric such that the volume density only depends on cylindrical radius and vertical position, ρ = ρ(R, z).In an axisymmetric and separable mass distribution, Poisson's equation is where a R and a z are the radial and vertical acceleration, respectively, and the second line replaces the radial terms with an expression involving the circular velocity, v c (R).We note that the midplane density, at any radius, ρ 0 (R), can be estimated using the fact that the z-derivative of Equation 14evaluated at z = 0 is simply related to our frequency parameter Ω 0 , where we have neglected the term involving the radial gradient of the circular velocity.We tested the impact of including this radial dependence in our results, and found a maximum impact on the order of ∼ 10%, which is well within our measurement uncertainties.
Integrating the expression in Equation 16over z, we obtain an expression for the surface mass density where we emphasize that Σ(z) is the total mass density within a vertical column (stars, gas, and dark matter).In the following section, we utilize this method to estimate key quantities about the mass distribution of the Galactic disk.We note that with our current method, we do not measure the circular velocity curve v c (R), so we adopt a linear fit to the circular velocity measured in Eilers et al. (2019) over the radius range considered below (6 < R < 13 kpc).

Measuring the contribution of dark matter to the surface and volume mass density
If knowledge of the baryonic contribution to the column mass density, Σ bary , is available, it is then straightforward to determine the contribution from dark matter to the total surface mass density, Σ tot = Σ bary + Σ DM .Similarly, under the assumption that the density of dark matter, ρ DM , is constant below a given vertical height, |z|, it is also possible to directly infer ρ DM using In the following Section, we will use estimates of Σ bary from the literature (e.g., McKee et al. 2015) to compute Σ DM and ρ DM around the Sun, employing estimates of the total surface mass density computed using element abundances.

Measuring the disc scale length
With measurements of the total volume mass density at the midplane, ρ(z = 0), across radius, it is also possible to estimate the disk's exponential scale length, h R , and the total volume mass density at the midplane at the solar radius, ρ(z = 0) ⊙ , via the following relation (Kuijken & Gilmore 1989) We will use this equation in Section 5.5 to estimate h R and ρ(z = 0) ⊙ , respectively.

Modelling the vertical element abundance trends across the Galactic disk
As an initial demonstration of our method, we use the modeling framework described above (Section 4 to fit the vertical element abundance gradient of ⟨[Mg/Fe]⟩ in bins of guiding radius for our APOGEE -Gaia parent sample.Figure 7 illustrates the result of optimizing our model applied to this data in three bins of guiding radius, 8 < R g < 9 kpc, 9 < R g < 10 kpc, and 10 < R g < 11 kpc, where each row in the figure corresponds to one guiding radius bin.The panels in each row of Figure 7 show the data, our optimized model, and the residual, from left to right.The data panel (left) shows ⟨[Mg/Fe]⟩ in bins of z, v z coordinates, where each pixel has a size of δz = 0.04 kpc and δv z = 1.76 km s −1 .The model panel (middle) shows our optimized model, fit to the data in the left column, evaluated for all pixel locations in the left panel.The residual panel (right) shows the difference between the data and the optimized model.
Comparing the top row of Figure 7 to the bottom row (i.e.stars in 8 < R g < 9 kpc versus 10 < R g < 11 kpc), we note that the optimized model prefers a more oblate contour shape at a given value of ⟨[Mg/Fe]⟩ (e.g., ⟨[Mg/Fe]⟩ ≈ 0.05) for the inner radial bin.Conversely, the outer radial bin (i.e. the bottom row) shows a more oblate/diamond-like contour shape far from the Galactic midplane.The fact that the residuals are generally small suggests that the functional form used is flexible enough to capture the overall structure of the data, at least in the regions of vertical phase space that are well-represented by the data.However, we note that the residuals do show a clear spiral pattern, which is most visible in the outer radial bin (i.e. the lower right panel).This closely resembles the morphology of the Gaia "phase spiral" (Antoja et al. 2018) in phase-space density, suggesting that this feature also manifests in chemical abundance space (see also Gaia Collaboration et al. 2023, , Frankel et al. in prep).This feature is further evidence of the complex chemicaldynamical structure of the Galactic disk, but we reserve any study of the chemical z, v z spiral for future work.
In the subsections that follow, we show summary properties of the mass distribution of the Milky Way disk as in-ferred from the most constraining element abundance ratio, [Mg/Fe], in all bins of guiding-center radius.

Measuring the vertical acceleration field and surface mass density across the Galactic disk
In this Section, we measure the vertical acceleration, a z (z), and surface mass density, Σ(z), as a function of height above the Galactic plane in bins of guiding-center radius across the disk for all six element abundance ratios considered above.In detail, we use the abundance ratios shown in Figures 2 and 4 and use vertical kinematic data for stars in overlapping, 2 kpc-wide bins of guiding radius, R g , spaced by 0.25 kpc from 6 to 12 kpc.We evaluate the quantities (a z and Σ) as a function of vertical position between z = 0−1.2kpc in steps of 50 pc.We truncate at z = 1.2 kpc as this is the approximate height to which our model represents well the data (see the residuals panels in Figure 7 and Price-Whelan et al., in prep).We use the bootstrap resamplings to estimate our measurement error on the resulting acceleration and surface density values; we compute the 50 th and [16 th ,84 th ] percentiles and use these as our representative median and uncertainty values for all measurables.
The left panel of Figure 8 shows our inferred mean vertical acceleration values when using our best element abundance ratio ([Mg/Fe]).From this, it is clear evident that a z varies with both R g and z, revealing beautiful trends across radius and height above the plane.Close to the Galactic midplane (z ∼ 0.1 kpc), the vertical acceleration is almost zero for all guiding-centre radii, as is to be expected for any model.
In a similar vein, the right panel of Figure 8 shows estimates of the surface mass density given our measurements of a z and Equation 19 for a range of z and R g (solid lines).As a comparison, we also show the Σ(z) values (dashed lines) computed using the MilkyWayPotential2022 in gala (Price-Whelan 2017a).We find that the inferred vertical surface mass density profile of the disk is generally consistent with this potential model at low z.At R g = 6 kpc and z < 0.2 kpc, we estimate a surface mass density Σ(z < 0.2 kpc) ∼ 80 M ⊙ pc −2 .For similar heights but larger guiding-center radius (R g = 10 kpc), we find a much lower surface mass density estimate (Σ(z < 0.2 kpc) ∼ 20 M ⊙ pc −2 ).Around the Sun (R g ∼ 9 kpc), we estimate a value of the total surface mass density of Σ ⊙ (z) ∼ 70 M ⊙ pc −2 at z = 1.1 kpc.
In the following subsections, we set out to constrain this quantity in more detail and compare it to estimates from the literature.We will also estimate the contribution of dark matter to the surface mass density around the Sun by using our estimates from Figure 8 and subtracting from it values of the stellar (and stellar remnants) surface mass density at around the Solar Neighborhood, Σ ⊙,bary (z).

The total and dark matter contribution to the surface mass density in the solar neighborhood
In this section, we focus on the surface mass density around the solar neighborhood and estimate the contribution to this value from dark matter.To do this, we select stars with precise abundance measurements within a narrow range of radius around the Solar neighbourhood (i.e., |R − R ⊙ | < 1 kpc for R ⊙ = 2.75 kpc).This cut yields 37,729 stars.We then bootstrap sample with replacement 25 trials, and compute the 50 th and [16 th ,84 th ] percentiles and use them as our representative median and 1σ measurement and uncertainty.Upon obtaining our model fits, we repeat the procedure from Section 4.2, and estimate the vertical acceleration as a function of vertical height, binning the data in 50 pc bins spanning from z = 0 kpc to z = 1.2 kpc.As a last step, we convert the vertical acceleration into an estimate of the total surface mass density.
A commonly referenced value in the literature is Σ ⊙ (z = 1.1 kpc), the total surface mass density around the Solar neighborhood at a height of z = 1.1 kpc.The reason for the choice of 1.1 kpc is historic (Kuijken & Gilmore 1991), and is hypothesized to be the height where the dark halo starts to dominate the surface mass density of the disk.For this reason, in this Section we will use this value to compare to previous measurements.

The total volume density at the midplane across the Galaxy disk
Another important quantity in the Milky Way is the volume mass density throughout the Galactic disk.Here, we use our inferred asymptotic midplane frequency values, Ω 0 , to measure the total volume mass density at the Galactic midplane (Equation 17) for [Mg/Fe].Figure 10 shows our Figure 7.A demonstration of our modeling procedure (Section 4) that fits ⟨[Mg/Fe]⟩ trends in vertical kinematics in three bins of guiding radius.Each row of panels corresponds to a given guiding-center radius bin: From top to bottom, 8 < Rg < 9 kpc, 9 < Rg < 10 kpc, and 10 < Rg < 11 kpc.The left panels show the mean [Mg/Fe] abundance in bins of (z, vz).The middle panels show the optimized OTI models.
The right panels show the residuals (i.e. the difference between the data and the optimized model, in units of [Mg/Fe] dex).The model captures the general trends of chemical gradients despite the obvious sampling limitations, but a clear spiral pattern is apparent in the residual plots.This feature illustrates that there is a chemical phase-space spiral in the Galactic disk, which is likely related to the kinematic phase spiral (see, e.g., Antoja et al. 2018).
resulting measurements of the total volume density at the midplane, ρ(R, z = 0) as a function of guiding-center radius, R g .The median and [16 th ,84 th ] percentiles from 25 bootstrap with replacement samples are shown as the solid line and shaded regions.The general trend shows a declining value of ρ(z = 0) with increasing guiding-center radius, as expected.At R g = 6 kpc, we determine a value of ρ(z = 0) ∼ 0.3 M ⊙ pc −3 .This value drops to ρ(z = 0) ∼ 0.08 M ⊙ pc −3 at R g = 9 kpc, and ρ(z = 0) ∼ 0.02 M ⊙ pc −3 at R g = 12 kpc.[6,7,8,9,10,11,12] kpc).Also shown here as dashed lines are the analytical values computed using the MilkyWayPotential2022 in Gala (Price-Whelan 2017a).Using orbital torus imaging, we empirically map the acceleration field and surface mass density across the Galaxy disk with element abundances.

The scale radius
Table 1.From left to right, the work or element abundance ratio used to compute the the total surface mass density at the Solar neighbourhood and z = 1.1 kpc, the stellar (and stellar remnants) surface mass density at the Solar neighbourhood and z = 1.1 kpc, and the dark matter surface mass density at the Solar neighbourhood and z = 1.1 kpc.Our estimates of Σ⊙,DM, (z = 1.1 kpc) are computed assuming Σ ⊙,bary (z = 1.1 kpc) = 47 ± 3 M⊙ pc −2 (Bland-Hawthorn & Gerhard 2016).For reference, we also list the estimates from previous studies (Kuijken & Gilmore 1991;Bovy & Rix 2013;Zhang et al. 2013;Bienaymé et al. 2014;Bland-Hawthorn & Gerhard 2016;Binney & Vasiliev 2023).for Mg agree well with previous estimates of Σ⊙(z = 1.1 kpc) around the solar neighbourhood computed using other methods.We note that previous works determined these values with smaller samples of stars, on the order of Nstars ≲ 10, 000, and the uncertainties are likely underestimated.Our measurement is computed using a sample of 37,729 stars; the uncertainties on our measurement will only become smaller with the advent of more data (e.g., SDSS-V/MWM Kollmeier et al. 2017).
Using these values and equation 21, we optimize 3 the exponential model to find the best fitting parameters representing the scale radius, h R , and the midplane total volume mass density at the solar radius, ρ(z = 0) ⊙ .We obtain a value (for the low-α disc) of h R = 2.24 ± 0.06 kpc, and ρ(z = 0) ⊙ = 0.081 +0.015 −0.009 M ⊙ pc −3 .The scale radius estimate we obtain is in good agreement with previous results (h R ∼ 2 − 3; e.g., Bovy & Rix 2013).5.6.The total and dark matter contribution to the total volume mass density in the solar neighbourhood 3 To optimize equation 21 we used the curve-fit function in scipy.In Figure 11, we show our measurement of ρ ⊙ (z = 0) in the solar neighbourhood.For reference, we also show the value from previous work (Holmberg & Flynn 2000;Garbari et al. 2012;McKee et al. 2015) for both ρ ⊙ (z = 0) and ρ ⊙,DM (z = 0).We estimate a value for the total midplane volume mass density at the Sun of ρ ⊙ (z = 0) = 0.081 +0.015 −0.009 M ⊙ pc −3 .This is in good agreement with the value reported in recent studies (ρ ⊙ (z = 0) ∼ 0.08 − 0.1 M ⊙ pc −3 ; e.g., Holmberg & Flynn 2000;Garbari et al. 2012;McKee et al. 2015).
In addition, using Equation 20, we also compute the value of the local volume mass density contributed by dark matter to be ρ ⊙,DM (z = 0) = 0.011 ± 0.002 M ⊙ pc −3 , using the value of Σ ⊙,bary (z = 1.1 kpc) M ⊙ pc −2 from McKee et al. (2015).Our estimate of ρ ⊙,DM (z = 0) is also in good agreement with the wealth of previous measurements (de Salas & Widmark 2021, and references therein).Our results are summarised in Table 2.

The midplane position and vertical velocity at the Solar Neighborhood
In addition to the main results in this paper, other quantities that we can directly measure from our fitting procedure are the midplane position, vertical velocity, and vertical frequency around the Solar neighbourhood.To do so, we take the same stars from Section 4.2 and measure z ⊙ , v z,⊙ , and Ω z,⊙ around the Sun.We perform this exercise using  3 shows the median and [16 th ,84 th ] percentile values from fitting the data using each element abundance ratio and bootstrapping with replacement 25 times.We find that our fits return measurements for z ⊙ and v z,⊙ that are close to previous estimates (e.g., z ⊙ = 29.5 ± 1 pc and v z,⊙ = 8.66 ± 0.05 km s −1 , see Schönrich et al. 2010;Bennett & Bovy 2019).This leads to an estimated value for the Solar vertical frequency of Ω z,⊙ = 0.30 ± 0.01 Myr −1 .Our measurements are in good agreement with previous estimates (see Table 3).
6. DISCUSSION 6.1.First empirical estimates of the acceleration field and mass density across the disk from stellar labels We believe that this article presents the first dynamical measurement of the Milky Way disk's profile empirically from stellar labels.The ability to infer dynamical properties of the disk from kinematics via contours of constant density was proposed approximately four decades ago (Kuijken & Gilmore 1989).However, until the advent of Gaia, insufficient numbers of stars have prohibited such analysis.In this work, we have inferred dynamical properties of the Milky Way disk by modelling the mono-abundance contours, which reveals the vertical orbital structure of stars across the Galaxy.We have shown that this unique synergy between stellar labels and kinematics enables measuring the vertical acceleration field, and therefore the surface mass density profile of the Galactic disk without the assumption of a functional form for the gravitational potential.Our method is also much less affected by spatial selection function effects when compared to previous techniques (e.g., Kuijken & Gilmore 1991;Bovy & Rix 2013;Zhang et al. 2013;McMillan & Binney 2013;Sanders & Binney 2015;Li & Widrow 2021), as we are much less sensitive to spatial selection functions the survey data we use.

The impact of survey selection effects on our measurements
In detail, the effective selection function for our APOGEE -Gaia sample depends on position, magnitude, and color of the stars.However, in our analysis below, we only model the element abundance trends as a function of orbital properties and not the absolute distribution of orbital properties represented in our sample (i.e.we model the element abundance conditioned on the orbital actions).Because of this, if the effective selection function factorizes in terms of its dependence on position and its dependence on element abundances, our analysis will not be sensitive to selection biases.This is especially important for cross-matches to APOGEE because of the pencil-beam nature of the survey (see for example Lane et al. 2022).We do expect the effective selection function of our sample to approximately factorize in this way because of the uniform selection criteria APOGEE used to define its core samples of targets, which only depends on J − K s color and H-band magnitude (Zasowski et al. 2013(Zasowski et al. , 2017;;Beaton et al. 2021;Santana et al. 2021).However, even with simple selection criteria, stars at different evolutionary stages are represented differently as a function of sky position and distance.This could introduce correlations in the selection function as the reported chemical abundances depend on stellar parameters in complex ways (Eilers et al. 2019(Eilers et al. , 2022)).To mitigate this issue, we select stars in a narrow range of surface gravity log g along the red giant branch, 1.5 < log g < 3.5.Moreover, in the low-α disc, stars of constant mean element abundance and orbit (e.g., R g ) present similar age distributions (Ness et al. 2019).In summary, the advantage of OTI over other methods is that it does not require detailed knowledge of the positional selection function, provided that the sample selection does not strongly impact the chemical abundance ratios as a function of phase space selections.For the case of our APOGEE disk sample, we argue that this condition is sufficiently satisfied.

Orbital torus imaging in the context of chemical-dynamical models of the Galaxy
In this work we have built on the Orbital Torus Imaging (OTI) framework (Price-Whelan et al. 2021) to map the orbital structure of disk stars using their element abundance ratios.Our framework effectively models the conditional distribution over some stellar label or element abundance ratio, [X/Y], at (conditioned on) a given location in phasespace, i.e. p([X/Y] | w) (where we use w to represent the phase-space coordinates; here, these would be w = (z, v z )).Other techniques have constructed models that unify stellar labels with the distribution function of stars to create chemo-dynamical models of the Galaxy.One example of this is the dynamical modeling of "mono-abundance popu-lations" (MAPs; Bovy et al. 2016;Mackereth et al. 2017), which constrains the stellar distribution function conditioned on element abundances, or p(w | [X/Y]).Another example is the "extended distribution function" (eDF) framework used (Sanders & Binney 2015;Binney & Vasiliev 2023), where here the model uses and constrains the joint distribution over kinematics and element abundances, p(w, [X/Y]).In both of these cases, the phase-space coordinates appear on the left of the conditional, in the distribution being modeled.The power of the OTI framework is that we model the element abundance (or stellar label) distribution conditioned on phasespace coordinates, which, as mentioned above, requires much less knowledge of the survey selection functions that supply our data.

The impact of disequilibrium and phase spirals
As a result of our method, in the residual plots of the 10<R g <11 kpc (bottom row of Figure 7, but also other R g bins not shown) we have detected a spiral pattern that matches the phase space spiral discovered by Antoja et al. (2018) (but also hinted at by Widrow et al. 2012).This result highlights that in addition to the over-arching smooth distribution, there is substructure in the chemical gradients of disk populations in z, v z space.The properties and origin of the phase space spiral are currently a debated topic in the field of Galactic dynamics using Gaia (e.g., Binney & Schönrich 2018;Li & Shen 2020;Bland-Hawthorn & Tepper-García 2021;Hunt et al. 2022;Widmark et al. 2022b,c,a;Antoja et al. 2023;Frankel et al. 2023;Tremaine et al. 2023) and numerical/cosmological simulations (Laporte et al. 2019;Khoperskov et al. 2019;Hunt et al. 2021;Grand et al. 2022).Although we have not explored the properties of this chemical phase spiral in this work, it would be interesting to model the shape, size, and amplitude of the chemical spiral across the Galaxy disk, and compare its properties to the density phase spiral (see Frankel et al. in prep for an example).We reserve this exploration for future work.
In principle, the effect of disequilibrium in the Milky Way disc, and the phase spirals that arise from it, could bias our measurements of the midplane vertical position and velocity.However, we choose to ignore this effect as we argue that it will be small (our model is able to capture well the distribution of mean element abundance ratio in vertical kinematic space.If the effect of the chemical spiral was large, this would not be the case).In the future, it would be worth-while accounting for the effect of the chemical spiral in our method, using a model that can also varies with vertical phase, θ z .

Future prospects of Orbital Torus Imaging
In this paper we have used the Orbital Torus Imaging (OTI) framework to directly map orbits and mass in the Galactic disc using one element abundance ratio with the high-est signal to noise gradient, [Mg/Fe] (see Figure 5).However, there are clear observed element abundance gradients in many other element abundance ratios as well (see Figure 12).In principle, all element abundance ratios should provide constraints on the orbit shapes and mass distribution, as done here.However, most element abundance ratios are not independent (e.g., Ness et al. 2022), and therefore cannot be naively combined in a likelihood sense.In future work, it will be interesting to identify and use independent combinations of abundance ratios as our stellar labels, especially those originating from different nucleosynthetic channels (i.e., SN II, SN Ia, AGB winds, neutron-star mergers), which can be separated by studying the element abundance distributions of stellar populations (e.g., Weinberg et al. 2022;Griffith et al. 2023).
This work has made use of spectroscopic data from the APOGEE survey, but many other large spectroscopic surveys are operating or coming online soon (e.g., MWM/SDSS-V: Kollmeier et al. 2017;WEAVE: Dalton et al. 2012;4MOST: de Jong et al. 2012).The next generation of surveys will deliver precise stellar label (element abundance) measurements for many millions of stars across large swaths of the Milky Way.This will extend the reach of the OTI framework to new parts of the Galaxy and enable even more precise constraints on the Galactic mass distribution.

CONCLUSIONS
In this paper we use the latest APOGEE and Gaia data to measure the bulk properties of the mass distribution around the Galactic disk.We use a novel modeling framework (Orbital Torus Imaging; Price-Whelan et al. 2021, Price-Whelan et al. in prep.) that leverages correlations between stellar element abundance measurements (provided by APOGEE ) with stellar kinematics (provided by Gaia) to infer the orbit structure of stars in the vertical phase-space, z, v z .We never use a parameterized model of the global Galactic potential and instead partition the data into bins of guiding-center radius to perform our measurements with the vertical phase-space coordinates.Our main results and findings are summarized below: • The stellar label (i.e., mean element abundances) gradients in the vertical (z, v z ) phase space can be well mapped by an empirical model that includes a loworder Fourier expansion distortion to an ellipsoidal radius (Figure 7).From this mapping, it is possible to estimate the vertical acceleration field and surface mass density across the low-α disk (Figure 8).
• At the Solar radius and at a vertical height of z = 1.1 kpc from the midplane, we estimate an average value for the total surface mass density of Σ ⊙ (z = 1.1 kpc) = 72 +6 −9 M ⊙ pc −2 using our the element abundance ratio with highest constraining power, [Mg/Fe] (Figure 9).Using the value from McKee et al. (2015) for the contribution from baryons, we estimate a contribution of dark matter of Σ ⊙,DM (z = 1.1 kpc) = 24 ± 4 M ⊙ pc −2 .
• We fit an exponential profile to the ρ(z = 0)-R g relation (Fig 10), and obtain a scale length of the low-α disc of h R = 2.24 ± 0.06 kpc.
• Our model directly fits the vertical position and velocity at the Sun.Using [Mg/Fe] as our most reliable element abundance ratio, we find a midplane vertical position of z ⊙ = 29 ± 1 pc, a vertical velocity relative to the midplane of v z,⊙ = 8.66±0.05km s −1 , and a midplane vertical frequency Ω z,⊙ = 0.30 ± 0.01 Myr −1 .
• In the residual from our best model fits (right panels of Figure 7), we find a clear imprint of the vertical phasespace spiral, but now detected in element abundances.This "chemical spiral" matches in position and structure the recently-discovered Gaia phase space spiral (Antoja et al. 2018).
Orbital Torus Imaging (OTI) is a promising new method for inferring the orbital structure and mass distribution of the Galaxy that does not rely on parameterised functions of the Galactic potential and is much less affected by the spatial selection function.Instead, it relies on the ability to image orbital tori from gradients of stellar labels in phase space.In this work, we have used this technique to map the acceleration field and surface mass density of a region of the Galaxy disc in vertical phase space.To do so, we have had to assume the potential can be decomposed in R and z.However, the next steps should involve combining the radial dependence in our model.Along those lines, insofar we have only been able to map the orbital structure of the Galaxy disc with a limited sample.The advent of the Galactic Genesis Milky Way Mapper survey (SDSS-V; Kollmeier et al. 2017) will deliver over ≳ 5 million all-sky spectra with APOGEE precision, that will be perfect for measuring the structure of the Milky Way with OTI.
and Debra for their constant support.The Flatiron Institute is funded by the Simons Foundation.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.

DATA AVAILABILITY
All APOGEE and Gaia data used in this study are publicly available and can be downloaded directly from https://www.sdss4.org/dr17/and https://gea.esac.esa.int/archive/,respectively.
Figure 1.Our parent sample of red giant branch stars in the low-α disk sequence used in this work.The left panel shows a 2D histogram of the positions of stars in our sample projected onto the Galactic plane (i.e. in Galactocentric x, y coordinates); The position of the Sun is marked with the ⊙ symbol.The middle panel shows the selection used to define the low-α sequence using a cut in element abundance space, i.e., the [Mg/Fe] vs. [Fe/H] plane.This chemical cut was performed in order to restrict our sample to stars that are likely part of the Galaxy's thin disk.The right panel shows the spectroscopic stellar parameters, effective temperature T eff and surface gravity log g.In the middle and right panels, the full APOGEE -Gaia sample is shown as the faint background histogram.metric (parallax and proper motion) measurements released in Gaia DR3 (Gaia Collaboration et al. 2022).Together, the APOGEE and Gaia data supply full 6D phase space information from which kinematic and orbital properties can be derived.The APOGEE data additionally provide detailed chemical abundance measurements for up to ∼ 20 different elemental species spanning the α, odd-Z, iron-peak, and s-process nucleosynthetic channels.
ure 1), and thus are likely part of the Milky Way's thin disk.Our selection yields a total of 172,656 stars in our APOGEE -Gaia parent data sample; These stars are shown in Figure1compared to the full APOGEE -Gaia cross-matched catalog in spatial location within the Galaxy (left panel), element abundance measurements ([Mg/Fe] vs. [Fe/H]; middle panel) and spectroscopic parameters (right panel).We use the 6D phase space information 1 to compute Galactocentric Cartesian coordinates assuming a solar velocity v ⊙ = (8.4,251.8, 0) km s −1 and distance to the Galactic center R 0 = 8.275 kpc(GRAVITY Collaboration et al. 2022).

Figure 2 .
Figure 2.Each panel shows the mean element abundance, ⟨[X/Y]⟩, of stars in bins of vertical phase space coordinates, z, vz, for the abundance ratios stated in the panel titles.All stars are red giant branch stars in the low-α disk sequence (see Section 2) around the Sun (i.e., |R − R⊙| < 0.5 kpc) with close-to-circular orbits (|R − Rg| < 1 kpc and |vR| < 50 km s −1 ).The typical spread in abundance values in each pixel (as estimated with the mean absolute deviation) is ≲ 0.02 dex, except for [Fe/H] where it is ≲ 0.08 dex.We only show the abundances that have the smallest scatter in our low-α sample and/or are best determined by ASPCAP/APOGEE.than stars with larger J z values.The right panel shows the same vertical phase plane but now shows contours of constant ⟨[X/Y]⟩, computed by painting on chemical abundance values to stars based on a linear relation between ⟨[Mg/Fe]⟩ and the vertical action J z (see Figure 4).Under assumptions enumerated below (see Section 4.1), contours of constant ⟨[X/Y]⟩ in z, v z delineate orbits.In what follows, we use this link between element abundances and orbits to model the mass distribution across the Galactic disk.

Figure 3 .
Figure3.This figure demonstrates, with simulated data, that contours of constant element abundance ratio (right panel) nearly delineate orbits in z, vz space (left panel).Left: This panel shows eight orbits with different values of the vertical action Jz computed in a simple Milky Way mass model.For all of these orbits, the value of the other actions are set such that the orbits have zero radial action, JR = 0, and a constant azimuthal action set to the solar value, J ϕ = J ϕ,⊙ .Orbits with lower Jz values are more elliptical and reach smaller maximum heights from the Galactic plane, zmax.We note that this is our only usage of a parameterized Galactic mass model in this work, and this is only for demonstration purposes.

Figure 4 .
Figure 4. Column-normalized histograms of element abundance values as a function of approximate maximum vertical height, ∼ zmax, for stars with orbits similar to the Sun (i.e., Rg − Rg,⊙ < 0.5 kpc) with low vertical velocities, |vz| < 10 km s −1 , from our parent sample.Each panel shows a different element abundance ratio indicated on the vertical axis.The dashed (red) line shows the running mean as a function of ∼ zmax.A steeper gradient in the mean abundance of a given panel indicates that a given abundance ratio will lead to better constraints on the orbit structure.The abundance ratio with the steepest and smoothest trend between element abundance and vertical kinematics is [Mg/Fe].

Figure 5 .
Figure 5. Signal-to-noise measurement, computed by taking the ratio of the slope, ∆⟨[X/Y]⟩/∆ ∼ zmax, with the average variance around the median line, σ 2⟨[X/Y]⟩ , for all element abundance ratios supplied by APOGEE.The higher the signal-to-noise ratio, the more constraining power an element abundance ratio has for modelling the Milky Way disc in z, vz space.Quantitatively, [Mg/Fe] is the element abundance ratio with the highest constraining power.
Figure 6.The same as Figure 2 but now showing only ⟨[Mg/Fe]⟩ and binned into eight overlapping bins of guiding-center radius, Rg.Each panel shows a 2 kpc-wide bin of Rg, increasing from the top left (inner galaxy) to the bottom right (outer galaxy), and the Sun would appear in the top right panel.The aspect ratio of each panel is kept constant.Stars orbiting in the inner galaxy have more oblate contours of constant abundance as compared to stars further out, which is especially apparent in for the low-[Mg/Fe] contour values.For example, compare the darker colors in the two panels in the far left column.This is a result of the surface mass density, as stars in the inner galaxy feel a stronger gravitational pull and will therefore reach smaller vertical heights for a given vertical velocity.

FigFigure 8 .
Fig 10 shows the total volume mass density estimates we obtain at the midplane (z = 0) across guiding-center radii.

Figure 9 .
Figure 9.Estimated median and [16 th ,84 th ] percentile values of the total (circles) and dark matter (squares) surface mass density at z = 1.1 kpc for stars in the solar neighbourhood, computed by modelling 25 bootstrapped samples with replacement of the ⟨[Mg/Fe]⟩ distribution in z, vz space, using the most constraining element abundance ratio, [Mg/Fe].For comparison, we also show previous estimates of this quantity from previous work (see Table 1), including the values of the surface mass density from baryonic material (diamonds).The values of the dark matter surface mass density are computed by subtracting the Σ ⊙,bary (z = 1.1 kpc) = 47 ± 3 value from McKee et al. (2015) to our estimates of the total surface mass density.Our estimates of the total and dark matter Σ⊙(z = 1.1 kpc)for Mg agree well with previous estimates of Σ⊙(z = 1.1 kpc) around the solar neighbourhood computed using other methods.We note that previous works determined these values with smaller samples of stars, on the order of Nstars ≲ 10, 000, and the uncertainties are likely underestimated.Our measurement is computed using a sample of 37,729 stars; the uncertainties on our measurement will only become smaller with the advent of more data (e.g., SDSS-V/MWMKollmeier et al. 2017).
Figure 10.Total volume density at the midplane ρ(z = 0) computed using [Mg/Fe] as our element abundance tracer as a function of guiding-center radius.The solid line indicates the median of 25 bootstrap with resampling, whilst the shaded regions show the [16 th ,84 th ] percentile ranges.There is a declining density at the midplane with increasing guiding-center radius, which is best fitted by a scale radius of hR = 2.24 ± 0.06 kpc.The dashed line shows an exponential profile assuming a value for the scale radius of hR = 2.24.The vertical dotted line corresponds to the Rg value of the Sun.

Figure 11 .
Figure11.Estimated median and [16 th ,84 th ] percentile values of the total volume mass density in the Solar Neighbourhood, computed by modelling the 25 bootstrapped samples with replacement of the ⟨[Mg/Fe]⟩ distribution in z, vz space.For comparison, we also show estimates fromHolmberg & Flynn (2000),Garbari et al. (2012), andMcKee et al. (2015).See de Salas & Widmark (2021) for a more extensive review of ρ⊙,DM.Our estimates of ρ⊙ and ρ⊙,DM are in good agreement with the latest previous measurements.

Table 3 .
From left to right, the work or element abundance ratio to compute z⊙, vz,⊙, and Ωz,⊙, the Solar midplane position, the Solar vertical velocity, and the Solar vertical frequency.