The Implications of Local Fluctuations in the Galactic Midplane for Dynamical Analysis in the Gaia Era

TL;DR

Local midplane fluctuations in the Milky Way induce significant systematic errors in stellar action calculations, affecting dynamical analyses in the Gaia era, and must be accounted for in large-scale Galactic studies.

Contribution

This paper identifies and quantifies the systematic errors caused by local midplane variations on stellar action calculations, highlighting their importance for Galactic dynamics analysis.

Findings

Midplane variations induce up to 25% error in vertical actions.

Simulations show midplane fluctuations of ~100 pc at the solar circle.

Observed vertical velocity variations suggest midplane shifts of 60-170 pc.

Abstract

Orbital properties of stars, computed from their six-dimensional phase space measurements and an assumed Galactic potential, are used to understand the structure and evolution of the Galaxy. Stellar actions, computed from orbits, have the attractive quality of being invariant under certain assumptions and are therefore used as quantitative labels of a star's orbit. We report a subtle but important systematic error that is induced in the actions as a consequence of local midplane variations expected for the Milky Way. This error is difficult to model because it is non-Gaussian and bimodal, with neither mode peaking on the null value. An offset in the vertical position of the Galactic midplane of $\sim15\,\text{pc}$ for a thin disk-like orbit or $\sim 120\,\text{pc}$ for a thick disk-like orbit induces a $25\%$ systematic error in the vertical action $J_z$. In FIRE simulations of Milky Way-mass galaxies, these variations are on the order of $\sim100\,\text{pc}$ at the solar circle. From observations of the mean vertical velocity variation of $\sim5\text{--}10\,\text{km}\,\text{s}^{-1}$ with radius, we estimate that the Milky Way midplane variations are $\sim60\text{--}170\,\text{pc}$, consistent with three-dimensional dust maps. Action calculations and orbit integrations, which assume the global and local midplanes are identical, are likely to include this induced error, depending on the volume considered. Variation in the local standard of rest or distance to the Galactic center causes similar issues. The variation of the midplane must be taken into account when performing dynamical analysis across the large regions of the disk accessible to Gaia and future missions.

Figures

Tables

Table 1. Table 1: Names and properties of the three orbits considered in this work, where z max subscript 𝑧 max z_{\text{max}} is the maximum height of the orbit, 1 2 ​ ( R max − R min ) 1 2 subscript 𝑅 max subscript 𝑅 min \frac{1}{2}(R_{\text{max}}-R_{\text{min}}) is the magnitude of the radial excursions of the orbit, and κ 𝜅 \kappa and ν 𝜈 \nu are the radial/epicyclic and vertical frequencies of the orbit. In the epicyclic approximation, A z = z max subscript 𝐴 𝑧 subscript 𝑧 max A_{z}=z_{\text{max}} and A R = 1 2 ​ ( R max − R min ) subscript 𝐴 𝑅 1 2 subscript 𝑅 max subscript 𝑅 min A_{R}=\frac{1}{2}(R_{\text{max}}-R_{\text{min}}) .

name	initial position	initial velocity	$J_R$	$J_{\varphi }$	$J_z$	$z_{\text{max}}$	$\frac{1}{2}\left(R_{\text{max}}-R_{\text{min}}\right)$	$\kappa$	$\nu$
	$\left(\mathrm{kpc}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{kpc}\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{kpc}\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{kpc}\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{kpc}\right)$	$\left(\mathrm{kpc}\right)$	$\left({\mathrm{Myr}}^{-1}\right)$	$\left({\mathrm{Myr}}^{-1}\right)$
thin-disk	(8, 0, 0)	(0, -190, 10)	40	-1500	0.69	0.12	1.3	0.049	0.093
thick-disk	(8, 0, 0)	(0, -190, 50)	33	-1500	23	0.85	1.2	0.048	0.061
halo	(8, 0, 0)	(0, -190, 190)	33	-1500	530	6.2	2.3	0.033	0.025

Table 2. Table 2: Stellar and gas disk scale heights of the Milky Way and the FIRE galaxies considered in this work (described in Section  3.1 ). For comparison, we also give the median softening lengths for the FIRE galaxies, computed for cold gas ( T < 1000 ​ K 𝑇 1000 K T<1000\,\text{K} ) and stars with | R − R 0 | < 0.5 ​ kpc 𝑅 subscript 𝑅 0 0.5 kpc \left|R-R_{0}\right|<0.5\,\text{kpc} and | z | < 1 ​ kpc 𝑧 1 kpc \left|z\right|<1\,\text{kpc} . We have assumed that R 0 = 8.2 ​ kpc subscript 𝑅 0 8.2 kpc R_{0}=8.2\,\text{kpc} .

galaxy	coldaa gas disk scale height	thin disk scale height	thick disk scale height	coldaa gas softening length	stellar softening length
	(pc)	(pc)	(pc)	(pc)	(pc)
Milky WaybbJurić et al. (2008); Bland-Hawthorn & Gerhard (2016)	40	300	900	$\cdots$	$\cdots$
m12iccSanderson et al. (2018)	800ddThe azimuthally averaged gas vertical density profile in m12i is nearly constant to this height, though individual regions show smaller scale heights and dense clouds.	480	2000	53.4	11.2
m12fccSanderson et al. (2018)	360	440	1280	57.2	11.2
m12mccSanderson et al. (2018)	250	290	1030	60.1	11.2

Equations

$$J_i\equiv \frac{1}{2\pi }{\oint }_{\text{orbit}}{p}_i\text{d}{x}_i\text{,}$$

$$\begin{array}{rl}R(t)& =R_g+A_R\sin (\kappa t+\alpha )\\ \varphi (t)& ={\Omega }_{c}t\\ z(t)& =A_z\sin (\nu t+\beta )\text{,}\end{array}$$

$$\begin{array}{rl}{v}_R(t)& =\kappa A_R\cos (\kappa t+\alpha )\\ v_{\varphi }(t)& =v_c\\ v_z(t)& =\nu A_z\cos (\nu t+\beta )\text{.}\end{array}$$

$$J_{\varphi }=R_g{v}_c\text{,}$$

$$\begin{array}{rl}{J}_R& =E_R/\kappa \\ J_z& =E_z/\nu \text{,}\end{array}$$

$$\begin{array}{rl}{J}_R& =\frac{v_R^2+{\kappa }^2(R-R_g{)}^2}{2\kappa }\\ J_z& =\frac{v_z^2+{\nu }^2{z}^2}{2\nu }\text{.}\end{array}$$

$$\begin{array}{rl}{J}_R& =\frac{\kappa A_R^2}{2}=\frac{v_{R,\text{max}}^2}{2\kappa }\\ J_z& =\frac{\nu A_z^2}{2}=\frac{v_{z,\text{max}}^2}{2\nu }\text{,}\end{array}$$

$$\frac{\Delta J_R}{J_R}=\frac{2(R-R_g)}{A_R^2}\Delta R+\frac{2v_R}{v_{R,\text{max}}^2}\Delta v_R\text{.}$$

$$\frac{\Delta J_{\varphi }}{J_{\varphi }}=\frac{\Delta R}{R_g}+\frac{\Delta v_{\varphi }}{v_c}\text{.}$$

$$\frac{\Delta J_z}{J_z}=\frac{2z}{A_z^2}\Delta z+\frac{2v_z}{v_{z,\text{max}}^2}\Delta v_z\text{.}$$

$$\begin{array}{rl}\frac{\Delta J_R}{J_R}& =\frac{2\Delta R}{A_R}\\ \frac{\Delta J_{\varphi }}{J_{\varphi }}& =\frac{\Delta R}{R_g}+\frac{\Delta v_{\varphi }}{v_c}\\ \frac{\Delta J_z}{J_z}& =\frac{2\Delta z}{A_z}\text{,}\end{array}$$

$$\frac{\Delta E}{E}\approx \frac{\Delta J_R}{J_R}+2\frac{\Delta J_{\varphi }}{J_{\varphi }}+\frac{\Delta J_z}{J_z}\text{.}$$

$$A\sin \left(\varphi +B\right)+C\text{,}$$

$$(x+\Delta x{)}^2+y^2=R_{\text{err}}^2\text{.}$$

$$\begin{array}{rl}{R}_{\text{err}}^2& =R^2-2R\cos \varphi \Delta x\\ ⟹\Delta R& \equiv \left|R_{\text{err}}-R\right|=\left|\cos \varphi \right|\Delta x\text{.}\end{array}$$

$$⟨\Delta R⟩=\frac{2}{\pi }\Delta x\text{.}$$

Full text

Milky Way disk – Milky Way dynamics – Milky Way evolution – Galaxy structure – Orbits – Stellar dynamics

††journal: ApJ††software: This work made use of the following software: astropy (Astropy Collaboration et al., 2013; Price-Whelan et al., 2018), gala (Price-Whelan, 2017; Price-Whelan et al., 2019), matplotlib (Hunter, 2007), numpy (Van Der Walt et al., 2011), scipy (Jones et al., 2001–), and tqdm (da Costa-Luis, 2019).

1 Introduction

Our understanding of the Milky Way is currently undergoing a revolution as a result of Gaia Data Release 2 (DR2). Recent major discoveries include the affirmation of remnants of a major merger (Koppelman et al., 2018; Belokurov et al., 2018; Helmi et al., 2018; Lancaster et al., 2019; Mackereth et al., 2019) hinted at in pre-Gaia work (e.g., Meza et al., 2005; Navarro et al., 2011), a phase-space “spiral” in the solar neighborhood (Antoja et al., 2018) possibly indicating local substructure infall (Binney & Schönrich, 2018; Laporte et al., 2019b), and a gap suggestive of perturbation by a dark matter subtructure in the tidal stream GD1 (Price-Whelan & Bonaca, 2018; Bonaca et al., 2018). These discoveries all indicate that the Milky Way’s stellar distribution, which demonstrably departs significantly from axisymmetry, is undergoing phase mixing and dynamical interactions across a range of spatial and temporal scales.

The assumption of a global, axisymmetric Galactocentric coordinate system (Binney & Tremaine, 2008) underlies much of the quantitative analysis of the mechanisms that give rise to these signatures. In order to construct such a system, the Sun’s relative position and velocity must be defined and measured both precisely and accurately. This involves determining the angular position of and distance to the Galactic center, the orientation of and distance to the Galactic midplane, and the local standard of rest (LSR). We review and discuss the observational efforts to measure these parameters in Section 4.2.

Once a Galactocentric coordinate system has been established and a six-dimensional (6D) phase space measurement of a star has been made, it is often desirable to convert this measurement into action space to concisely summarize its projected orbit, model the stellar distribution function, or find stars with similar dynamical properties. Actions are conserved quantities that describe the orbit of a star under the assumption of regular, bound orbits in a system where the equations of motion are separable in a particular coordinate system. They are the cyclical integral of the canonical momentum over its conjugate position:

[TABLE]

where $p_{i}$ are the conjugate momenta. Under the assumption of axisymmetry, $i=R,\phi,z$ are the radial, azimuthal, and vertical coordinates respectively in a cylindrical coordinate system. In a slowly-evolving axisymmetric potential, these actions are invariant and $J_{\phi}\equiv L_{z}$, where $L_{z}$ is the $z$-component of the angular momentum (Binney & Tremaine, 2008; Sellwood, 2014).

With the advent of 6D phase-space measurements over a relatively large ($\gtrsim 2$ kpc) volume from the Gaia satellite, the study of stellar actions has gained new popularity. One reason is dimensionality reduction — an individual stellar orbit is concisely described by three actions, as opposed to six phase space coordinates. Second, under the assumption of a phase-mixed system, the dynamical properties of a population of stars should be uniquely a function of their actions and independent of the conjugate angles. This allows one to use actions to study the relationship between orbital properties of stars and other intrinsic, and, at least partially, invariant properties such as age or metallicity (Beane et al., 2018; Ting & Rix, 2018; Sanders & Das, 2018; Gandhi & Ness, 2019; Das et al., 2019; Bland-Hawthorn et al., 2019). Actions also provide a convenient basis for constructing models of the stellar distribution function (e.g., Jeans, 1915; Villumsen & Binney, 1985; Trick et al., 2017), or for associating stars with similar dynamical properties, e.g., to potentially determine membership in moving groups.

If the system being considered departs from axisymmetry in a significant and/or non-adiabatic way, the actions computed using an axisymmetric approximation to the true potential can exhibit cyclic dependence on the orbital phase (or time at which the star’s position and velocity are observed), large-scale migration, or diffusion from their initial values. In the Milky Way, stellar actions are expected to diffuse on short time scales due to scattering with gas clouds and to evolve on longer time scales in the case of orbits near resonances with spiral arms, bar(s), and other large scale perturbations (Sellwood, 2014). For this reason, actions have been used to study stellar scattering in the Milky Way disk using the improved astrometry of Gaia DR2 and various age catalogues (Beane et al., 2018; Ting & Rix, 2018). Actions have also been used to study different models of spiral structure in the Milky Way (Sellwood et al., 2019). Characteristics of the distribution of stars in the extended solar neighborhood in action space are discussed in Trick et al. (2019).

The true Galactic potential is not strictly axisymmetric, beyond even the typically-quoted gas clouds, spiral arms, and bar(s). The presence of disk oscillations has been known since its observation in H i by Kerr (1957). The presence of a North-South asymmetry in stellar number density and velocity in the solar neighborhood has been found in SDSS data (Widrow et al., 2012), and in velocity in the LAMOST survey (Carlin et al., 2013; Xu et al., 2015), and RAVE data (Williams et al., 2013). Beyond the solar neighborhood, low-latitude overdensities have been observed, e.g., the Monoceros Ring (Newberg et al., 2002; Ibata et al., 2003; Crane et al., 2003; Slater et al., 2014; Morganson et al., 2016), an overdensity in the direction of the Triangulum and Andromeda galaxies (Martin et al., 2007; Sheffield et al., 2014; Xu et al., 2015; Price-Whelan et al., 2015) and A13 (Sharma et al., 2010; Li et al., 2017), and in the direction of the Galactic Center (Feast et al., 2014). There is evidence that these stellar populations originated from the disk (Price-Whelan et al., 2015; Sheffield et al., 2018; Bergemann et al., 2018). See, e.g., Laporte et al. (2018) for a discussion of these observations.

With the vast improvement in the quality of phase-space measurements due to Gaia the assumption of axisymmetry is increasingly inadequate (e.g., Antoja et al., 2018; Laporte et al., 2019b). The presence of the phase-space spiral itself is evidence of disk oscillations (Antoja et al., 2018), but the fact that its strength changes across the disk is further evidence (Laporte et al., 2019b; Bland-Hawthorn et al., 2019). Even if the axisymmetric assumption were close enough for many purposes, the parameters used in axisymmetric models of the Galactic potential may be biased by the non-axisymmetries in the disk.

Such variations in the density of gas and stars cause the local midplane position to vary as a function of radius and azimuth. Stars far from the Sun have a local midplane that differs from our local midplane extrapolated onto their position. Converting positions and velocities of more distant stars from a heliocentric to a Galactocentric coordinate system thus introduces a systematic bias in the $z$ coordinate. We show that this bias induces non-Gaussian errors in the actions computed for these stars. The further from the solar neighborhood the target star is, the more likely the mismatch will result in large systematic uncertainty, especially in the vertical action $J_{z}$. A similar argument applies to any remaining uncertainty in measurements of the Galactocentric radius of the Sun, and to variations in the LSR.

In Section 2, we describe the general impact coordinate system errors have on the measured actions. In Section 3, we examine the azimuthal variations of the midplane itself in examples from two classes of simulations: cosmological, hydrodynamical, zoom-in simulations of isolated Milky Way-mass galaxies from the Feedback in Realistic Environments (FIRE) collaboration111https://fire.northwestern.edu (Hopkins et al., 2014; Wetzel et al., 2016; Hopkins et al., 2018), and a controlled N-body simulation of a Sagittarius encounter with a galaxy otherwise tailored to the stellar mass, scale length, and scale height of the Milky Way (Laporte et al., 2018). In Section 4, we discuss the implications of midplane variations, and the resulting systematic uncertainty in the vertical action, for action-space analyses. We also estimate the expected midplane variations of the Milky Way based on the observed velocity variations and three-dimensional dust maps. We summarize our main results and conclude in Section 5.

2 Motivation

We first demonstrate the significance to action computations of a systematic offset in the determination of the Galactic midplane, distance to the Galactic center, or LSR. We will find that such offsets are especially important for disk-like orbits. The consequences we explore here may also arise from various other systematic errors. For instance, the axisymmetric Galactic potential model used in many works to compute actions may not be a good description of the true potential — or the parameters used may yield a potential that is systematically incorrect outside an original fitted region. In this work, we assume that the Galaxy is perfectly described by our model axisymmetric potential, and simply explore the consequences of offsets in the Galactocentric coordinate system.

2.1 Effect of Midplane Offset

We present an illustration of an orbit in Figure 1 to show how an inaccurate or erroneous determination of the midplane leads to a dependence on orbital phase of the value of the actions calculated from a point in phase space and an assumed potential model. The $x$-axis corresponds to orbital phase and the $y$-axis to vertical height. The solid gray curve indicates the true orbit of the star as it oscillates around the true midplane. The dashed gray line, offset from the true midplane, is the midplane location used by an observer to integrate the orbit of the star and estimate its actions. The model potential is otherwise identical to the one in which the star is actually moving.

Now suppose this observer makes a measurement of the star’s position and velocity at the teal point or the orange point (i.e., at two different orbital phases). Then, the teal and orange curves correspond to the orbits that the observer would compute for each point based on the potential model with the offset midplane. In action space, this would correspond to a different value of $J_{z}$ for the teal and orange points. In this way, assuming the wrong coordinate system induces a phase-dependence in the actions estimated for the star, which in the correct potential (in this example, the one with the correct midplane) should be phase independent.

This example uses an offset in $z$, but analagous effects occur from offsets in other coordinates, such as the distance to the Galactic center or the LSR. A similar effect in which actions gain time-dependence due to a time-varying potential was pointed out by Buist & Helmi (2015).

In reality one only ever measures the position and velocity of a star at a single instance in time. However, there are at least two possibilities for going beyond a single temporal snapshot of the phase-space distribution. If more than one star on the same or similar orbit could be identified (e.g., by identifying the remnants of a disrupted open cluster), then it would be possible to constrain the disk potential in analogy to how stellar streams are used to constrain the Galactic potential (e.g. Johnston et al., 1999; Ibata et al., 2001; Helmi, 2004; Sanderson et al., 2015; Küpper et al., 2015; Laporte et al., 2019a). This would be challenging since, in principle, one would want to use the stellar actions or other conserved quantities, but as we discuss in this work they are subject to substantial errors. Another possibility would be to use chemical tagging to identify a disrupted open cluster (Freeman & Bland-Hawthorn, 2002), but this may be difficult due to the large number of stars that are chemically indistinguishable though presumably born in different clusters (Ness et al., 2018). A second way to improve on a single temporal snapshot is to measure directly the accelerations of stars (Silverwood & Easther, 2018). While this would allow a direct measurement of the aliasing effect in Figure 1, the observational challenges will likely limit this type of measurement to a small number of local stars for the near future.

2.2 Epicyclic Approximation

Before turning to numerical methods, we derive analytic expressions for the systematic error in the actions induced from offsets in the position (e.g., the midplane or Galactic center distance) or velocity (e.g., the LSR) of the assumed coordinate system’s origin.

We use the epicyclic approximation, which assumes that the motion in the $z$ and $R$ components of the orbit are decoupled and follow simple harmonic motion about a circular and planar guiding orbit (Binney & Tremaine, 2008, Section 3.2 and references therein). The radius of this orbit is referred to as the guiding radius $R_{g}$. This approximation is an excellent description of the thin disk and a good description of the thick disk in an axisymmetric potential that ignores the influence of the Galactic bar and spiral arms. We also make the assumption of a perfectly flat circular velocity curve with $v_{c}(R)=v_{c}$, a good approximation near the solar circle (e.g., McMillan, 2017).

Under this approximation, we can write down the cylindrical components of the orbits as

[TABLE]

where $\kappa$ and $\nu$ are the radial/epicyclic and vertical frequencies, $\Omega_{c}\equiv v_{c}/R_{g}$ is the orbital frequency of the guiding center, $A_{R}$ and $A_{z}$ are the amplitudes in the radial and vertical coordinates, and $\alpha$ and $\beta$ are the initial orbital phases. Similarly, the velocities of the orbit are given by:

[TABLE]

In this case, the azimuthal action is (Binney & Tremaine, 2008, Section 3.5.3b):

[TABLE]

and the radial and vertical actions are

[TABLE]

where $E_{R}$ and $E_{z}$ are the energy per unit mass in the radial and vertical coordinates, respectively. Therefore,

[TABLE]

Using Equations (2) and (3), we can simplify this:

[TABLE]

where the last equality in each line comes from the fact that $v_{R,\text{max}}=\kappa A_{R}$ and $v_{z,\text{max}}=\nu A_{z}$.

Notice that while the value for each of $J_{\phi}$, $J_{R}$, and $J_{z}$ is phase independent, the contribution from the kinetic and potential terms in Equation (6) is phase dependent. Now assume that the coordinates $(R,z,v_{\phi},v_{R},v_{z})$ are offset by $(\Delta R,\Delta z,\Delta v_{\phi},\Delta v_{R},\Delta v_{z})$. We can then apply the standard propagation of errors formula to Equation (6) to determine the error in each of the actions. For $J_{R}$, the induced error is:

[TABLE]

For $J_{\phi}$, the induced error is:

[TABLE]

For $J_{z}$, the induced error is:

[TABLE]

We have ignored second order contributions.

Since most of the time stars will be at maximum amplitude (i.e., turnaround) in both $R$ and $z$, we can approximate the order of magnitude of the systematic error in the actions by

[TABLE]

where we have again ignored second order terms.

Because the energies in the radial and vertical coordinates are related to the radial and vertical actions by a constant (Equation (5)), the fractional error in these two energies is the same as the actions.222This ignores any potential error in the epicyclic and vertical frequencies. The energy in the azimuthal coordinate can be readily computed ($E_{\phi}=\frac{1}{2}v_{\phi}^{2}$). Using standard propagation of errors, we see that $\Delta E_{\phi}/E_{\phi}$ and $\Delta J_{\phi}/J_{\phi}$ (Equation (9)) differ by a $\Delta R/R_{g}$ term. Since $R_{g}\approx 8\,\text{kpc}$, we expect this term to be negligible and we therefore conclude that:

[TABLE]

In the remainder of this section, we compare our analytic estimates of the effect of a midplane offset on actions against numerical calculations. A numerical evaluation of the effect of velocity offsets on actions is deferred to future work, as we discuss in Section 3.4.

2.3 Numerical Methods

We now quantify the argument made in Section 2.1 using numerical computations of the actions for a range of orbits in a model Galactic potential. We compute actions as in Beane et al. (2018), using the code gala v0.3 to perform orbit integrations and conversion to action space (Price-Whelan, 2017; Price-Whelan et al., 2019). To compute actions we use the torus-mapping technique first presented by McGill & Binney (1990) and adapted by Sanders & Binney (2014) to calculate actions for an orbital time-series starting from a phase-space position $(x,v)$ and integrated in a potential $\Phi$. For our Galactic potential we use MWPotential, based on the Milky Way potential available in galpy (Bovy, 2015), which includes a Hernquist bulge and nucleus (Hernquist, 1990), a Miyamoto–Nagai disk (Miyamoto & Nagai, 1975), and a Navarro, Frenk, & White 1997 halo, and is fit to empirically match some observations. We use the Dormand-Prince 8(5,3) integration scheme (Dormand & Prince, 1980) with a timestep of $1\,\text{Myr}$ and integrate for $5\,\text{Gyr}$, corresponding to $\sim 20$ orbits for a Sun-like star.

We assume the Sun is located at $(8.2,0,0)\,\text{kpc}$. None of our orbit integrations depend on the value of the LSR in this toy potential (though this is important when using real data, since the conversion from heliocentric to Galactocentric coordinates depends on the LSR). In this potential, we have that the circular velocity $v_{\text{circ}}$ is $231\,\text{km}\,\text{s}^{-1}$ at the solar circle.

Other methods for computing actions are used in the literature. For example, the Stäckel Fudge method (Sanders & Binney, 2016), which uses a single Stäckel potential (with analytic actions) to approximate the Galactic potential (de Zeeuw, 1985; Binney, 2012), was used in many recent works exploring actions in the Galactic disk (e.g., Trick et al., 2019; Sanders & Das, 2018; Ting & Rix, 2018). For disk-like orbits, existing implementations of the Stäckel Fudge method are of acceptable accuracy, but since we also consider halo-like orbits in this work (where the Stäckel Fudge method is inaccurate) we choose to use orbit integration and torus mapping throughout (Sanders & Binney, 2016).

2.4 Quantification of the Midplane Effect

We now quantify how a systematic error in the Galactocentric coordinate system induces phase-dependence in the actions calculated from the observed position and velocity of a star. We consider three orbits in the model potential described in Section 2.3 that are typical of stars in the thin disk, thick disk, and halo. We summarize their initial positions in phase space, the actions computed by integrating their orbits in the correct potential, and other properties in Table 1. Each orbit, integrated without systematic coordinate errors, is plotted in Appendix A. We will refer to these orbits by their names (thin-disk, thick-disk, halo) henceforth.

We begin with thick-disk. Consider an observer who can measure the orbit’s phase-space position at many different times (and hence different orbital phases), but does so using a coordinate system in which the midplane is systematically offset in height by $100\,\text{pc}$ from its actual location. To model this we subtract the vector $(0,0,100)\,\text{pc}$ from each position in the orbit. This corresponds to an observer physically located at, e.g., the position $(8,0,0)\,\text{kpc}$ in the coordinate system of the true potential, but erroneously thinking they are located at $(8,0,0.1)\,\text{kpc}$.

We consider the observer making a measurement, integrating an orbit, and computing actions every megayear using the prescription above. However, we specify the star’s starting position using the systematically offset coordinate system. Essentially we are shifting and then reintegrating at each point along the original orbit. The actions computed using the offset coordinate system for each phase-space starting point are shown for the first gigayear of the orbit in the upper panels of Figure 2.333Occasionally the numerical scheme fails and very large actions are reported by gala — we perform a $4\sigma$ clip on each action to exclude such orbits, but this only excludes a total of $5$ orbits out of the $1000$ considered for Figure 2. Some numerical artifacts remain, but the vast majority of orbits are computed properly.

We also perform the same procedure in the lower panels but assuming an $x$ component offset of $100\,\text{pc}$, i.e., subtracting the vector $(100,0,0)\,\text{pc}$. This is equivalent to a measurement error in the distance from the Sun to the Galactic center.

Figure 2 shows that the actions computed in the offset coordinate systems oscillate as a function of the time/orbital phase at which the star’s phase-space position is observed. This time dependence comes even though the observer is using the correctly constructed, best-fit, static, axisymmetric potential. The relative size of the phase variation in each action depends on the direction of the systematic offset as well as the true values of the actions (i.e. the type of orbit). In reality, we will have one measurement of the phase-space position to work with, in which case the determination of the orbital phase in $R$ or $z$ is degenerate with the degree of systematic offset in that coordinate (see Figure 1). In the following we therefore quote percentile ranges for the possible values computed for each action as a proxy for the effect of these systematic errors in the coordinate system.

For a systematic offset in $z$ (upper panels), the 95th minus 5th percentiles are $2.2$ and $6.2\,\text{kpc}\,\text{km}\,\text{s}^{-1}$ for $J_{R}$ and $J_{z}$, respectively. As a fraction of the true values, these are $5.7\%$ and $86\%$, respectively. The spread induced in $J_{\phi}$ is negligible, as expected since $J_{\phi}$ only depends on the $x$- and $y$-components of the position and velocity of the stars.444In practice, however, $J_{\phi}$ is computed as part of the torus-fitting method. It is worth pointing out that a $100\,\text{pc}$ offset in an orbit with $z_{\text{max}}=850\,\text{pc}$ — a $12\%$ error — induced an $86\%$ spread in the computation of $J_{z}$.

For a systematic offset in $x$ (or distance to the Galactic center), the 95th minus 5th percentiles are $6.9$, $47$, and $0.71\,\text{kpc}\,\text{km}\,\text{s}^{-1}$ for $J_{R}$, $J_{\phi}$ and $J_{z}$, respectively. These are fractionally $21\%$, $3.1\%$ and $3.1\%$ in these actions, respectively, despite only a $1.2\%$ error in the distance to the Galactic center.

In Figure 3, we plot a histogram of the values of $J_{z}$ computed at different orbital phases for thick-disk (top panel) and thin-disk (bottom panel), assuming a $z$ offset of $100\,\text{pc}$ (as in the upper right panel of Figure 2). The true value is plotted as a vertical dashed line. The systematic error in $J_{z}$ induced by a systematic offset in $z$ is non-Gaussian and bimodal; neither of the modes is centered on the null value. In the case of thin-disk (bottom panel), we see that, in addition to the prior complications, the distribution is not even centered on the true value. This comes about when the midplane error is $\gtrsim z_{\text{max}}$, where $z_{\text{max}}$ is the maximum height of the orbit (equivalent to $A_{z}$ in the epicyclic approximation, see Section 2.2).

In Appendix C we plot the same histogram as in Figure 3, but for the distributions of $J_{R}$ induced by a $z$ offset (upper left panel of Figure 4) and the distributions of $J_{R}$ and $J_{\phi}$ errors induced by an $x$ offset (lower left and lower center panels of Figure 4, respectively). We find similar error distributions as in Figure 3, with the exception that the computed $J_{R}$ distribution induced by an $x$ offset more closely resembles a Gaussian distribution.

We now suggest a heuristic explanation for the shape of Figure 3. Consider first the thick-disk (top panel), where the offset in $z$ is much less than the $z_{\text{max}}$ of the orbit. The peaks in the distribution correspond to the turning points of the orbit (or points of maximum vertical amplitude), where $v_{z}\sim 0$ and where the star is on most of it orbit. This is why the distribution, which is calculated at evenly spaced time intervals, peaks at these values. For thin-disk (bottom panel), the offset in $z$ is comparable to $z_{\text{max}}$. Now, there will be some points in the orbit where $v_{z}=0$ and $z=0$ (in the erroneous coordinate system). At these points, the computed $J_{z}$ will vanish. The asymmetry and systematic offset then comes about because of the constraint that $J_{z}\geq 0$.555This argument is similar to ones given in cosmology for why gravity produces non-Gaussianity in the density field, since the density cannot become negative but it can grow arbitrarily large.

Gaussian summary statistics are clearly insufficient to describe the distribution shown in Figure 3. We therefore elect to measure this error by computing one half the 95th percentile minus the 5th percentile of the distribution of action values. We refer to this quantity as $\Delta J_{i}$ for each action and plot it in Figure 3 as a horizontal arrow anchored on the true action value. Because of the bimodality of the error distribution, this quantity roughly measures the distance from the true action value to the peak of one of the modes. Furthermore, this bimodality also implies that $\Delta J_{i}$ is not very sensitive to the exact percentiles used. This summary statistic does not reflect the bias induced when the midplane error is $\gtrsim z_{\text{max}}$.

We now repeat the same procedure as in Figure 2 but for systematic offsets between [math] and $500\,\text{pc}$ in the $z$ and $x$ components. In Figure 4, we report $\Delta J_{i}/J_{i}$ for the three different fiducial orbits in Table 1. The upper panels of Figure 4 shows the spread induced in each action for an offset in the $z$-component. In the lower panels we consider offsets in the $x$ component (i.e. the solar radius). The left, center, and right columns show the fractional spread in the values computed for $J_{R}$, $J_{\phi}$, and $J_{z}$, respectively.

In the upper middle panel of Figure 4, there is essentially no spread in the determination of $J_{\phi}$. This is expected since $J_{\phi}$ is independent of $z$ and is thus unaffected by offsets in $z$, as discussed earlier. Indeed, the result we found in Figure 2 for thick-disk holds for all orbit types. This is also a demonstration of the robustness of the integration and action calculation methods we use.

The upper right panel of Figure 4 shows that the fractional error in $J_{z}$ is more exaggerated for more planar (disk-like) orbits. For thin-disk, a systematic offset of $15\,\text{pc}$ in the $z$-coordinate results in a $25\%$ deviation in $J_{z}$, while a $120\,\text{pc}$ offset results in the same deviation for thick-disk. We find that halo is relatively insensitive to errors in the midplane, with only $\sim 15\%$ error in $J_{z}$ out to an offset of $500\,\text{pc}$.

For the offset in the solar radius (lower panels), the error is largest for $J_{R}$, with some deviations resulting in $J_{\phi}$ and relatively small deviations in $J_{z}$. In the lower center and lower right panels all three lines nearly overlap.

In each panel of Figure 4, where relevant, we include the estimation of the action errors derived under the epicyclic approximation from Equation (11), with $\Delta v_{\phi}=0$, as dashed lines in the color of each orbit. This equation is relevant since during most of the orbit the star will be close to maximum radial and vertical amplitude. Note that we consider an error in the $x$-coordinate $\Delta x$, which is not exactly the same as $\Delta R$. For observations of stars close to us, we have that $\Delta x\sim\Delta R$, but for the experiment performed in this section we consider observations of the star throughout its entire orbit. This introduces a factor of $2/\pi$ when converting from $\Delta x$ to $\Delta R$, which we derive in Appendix B.

The epicyclic approximation is a good predictor of $\Delta J_{z}$, even for halo. It performs similarly for $\Delta J_{R}$, now underpredicting for halo and slightly overpredicting for thin-disk. Note that for the particular orbits we chose, thin-disk has slightly larger $A_{R}$ than thick-disk, and so we actually expect the epicyclic approximation to perform slightly worse for thin-disk in this case. The epicyclic approximation underpredicts $\Delta J_{\phi}$ for all orbits.

To further understand the effect of the midplane error, we also plot the fractional error in $J_{z}$ as a function of $J_{z}$ for $z$ offsets of $10$, $50$, and $100\,\text{pc}$ (orange, teal, and green, respectively) in Figure 5. For each orbit, we set the initial position to be $(8,0,0)\,\text{kpc}$ and the initial velocity to be $(0,-190,v_{z})\,\text{km}\,\text{s}^{-1}$, where we vary $v_{z}$.666One can recover thin-disk, thick-disk, and halo by setting $v_{z}=10$, $50$, and $190$, respectively, giving $J_{z}\simeq 0.7,20,$ and $500\,\text{kpc}\,\text{km}\,\text{s}^{-1}$ (see Table 1). For a thin-disk-like orbit ($J_{z}\sim 0.7\,\text{kpc}\,\text{km}\,\text{s}^{-1}$), even a $10\,\text{pc}$ offset in $z$ is enough to induce a $\sim 20\%$ error in $J_{z}$. For larger values of $J_{z}$, the fractional errors are suppressed, but the induced error can still be large depending on how great the $z$ offset is. We also plot the epicylic prediction for $\Delta J_{z}/J_{z}$ from Equation (11) as dashed lines for each $z$ offset. We find that the epicyclic approximation matches the numerical estimate quite well.

3 Azimuthal Midplane Variations

The stellar midplane of the Galaxy should vary as a function of azimuth and Galactocentric radius due to small, local variations in the stellar density. Hints of this variation as a function of Galactocentric radius have been noted through their impact on the stellar velocity distribution pre-Gaia by Widrow et al. (2012), Carlin et al. (2013), and Williams et al. (2013) and recently post-Gaia by Friske & Schönrich (2019). As pointed out by, e.g., Goodman et al. (2014) and Anderson et al. (2019), among many others, the gas distribution in the Galaxy also shows significant density variation across the disk.

The local Galactocentric coordinate system is defined based on the location of the Sun relative to the midplane. Extending this coordinate system to a global one therefore introduces systematic errors in the $z$ components of stellar positions. As discussed in Section 2, this systematic error introduces errors in integrating orbits and computing actions.

We specifically consider azimuthal variations in the midplane at the solar circle, as defined by the stellar mass density. Since, to our knowledge, there are no direct empirical measurements of these variations in the Milky Way, we use example simulations from two classes of simulations to estimate the size of this effect.

One set are three zoom-in, cosmological hydrodynamical simulations of isolated Milky Way-mass galaxies from the FIRE collaboration, described briefly in Section 3.1. These include stars, gas, and dark matter in a fully cosmological setting but are not tailored to specific properties of the Milky Way (such as the scale height or scale length, or the details of the accretion history). We use these simulations to span the range of possibilities for azimuthal midplane variations.

The other set of simulations are isolated N-body simulations of interactions between the Milky Way and a Sagittarius-like dwarf galaxy, described briefly in Section 3.2. These include dark matter and stars and are tailored to existing measurements of the structure of the Milky Way’s disk and of the orbit and properties of the Sagittarius dwarf galaxy. Comparing the azimuthal midplane variations in the host galaxy of these simulations before and after the interaction with the Sagittarius-like object gives an idea of the effect of one minor merger whose properties are relatively well measured. Azimuthal variations of the mean vertical height of stars has been explicitly pointed out in a different simulation of a Sagittarius-like encounter by Gómez et al. (2013).

3.1 Description of FIRE Simulations

The FIRE cosmological hydrodynamic simulations (Hopkins et al., 2014, 2018) use the zoom-in technique (e.g., Katz & White, 1993; Oñorbe et al., 2014) to model the formation of a small group of galaxies at high resolution in a full cosmological context. Feedback from supernovae, stellar winds, and radiation from massive stars is implemented at the scale of star forming regions following stellar population synthesis models, generating galactic winds self-consistently (Muratov et al., 2015; Anglés-Alcázar et al., 2017) while reproducing many observed galaxy properties, including stellar masses, star formation histories, metallicities, and morphologies and kinematics of thin and thick disks (Hopkins et al., 2014; Ma et al., 2016, 2017; Wetzel et al., 2016; Garrison-Kimmel et al., 2018; Hopkins et al., 2018).

For this work, we focus on the three Milky Way-mass zoom-ins considered in Sanderson et al. (2018), which were simulated as part of the Latte suite and show broad agreement of many of their global properties with observations of the Milky Way (Wetzel et al., 2016; Garrison-Kimmel et al., 2018). The $z_{r}=0$ snapshots777In this work, to avoid confusion with the vertical height $z$, we refer to cosmological redshift as $z_{r}$. of these three simulations, The snapshots of these three simulations at cosmological redshift $z_{r}=0$, named m12i, m12f, and m12m, are publicly available alongside associated mock Gaia DR2 catalogues generated from them.888http://ananke.hub.yt

These simulations contain dark matter particles of mass $\sim 35,000\,\text{M}_{\odot}$, gas particles of mass $\sim 7000$ to $20,000\,\text{M}_{\odot}$, and star particles of mass $\sim 5000$ to $7000\,\text{M}_{\odot}$, with the lower end coming from stellar evolution (Sanderson et al., 2018). Softening lengths for dark matter and star particles are fixed at $112\,\text{pc}$ and $11.2\,\text{pc}$, respectively.999This is $2.8$ times the often-quoted Plummer-equivalent. The gas softening length is adaptive, but at $z_{r}=0$ the median softening length for cold ($T<100\,\text{K}$) gas particles around roughly solar positions (with galactocentric cylindrical radii within $0.5\,\text{kpc}$ of $8.2\,\text{kpc}$ and $\left|z\right|<1\,\text{kpc}$) is $53.4$, $57.2$, and $60.1\,\text{pc}$ for m12i, m12f, and m12m, respectively. These values are summarized in Table 2, along with measurements of the stellar and gas disk scale heights.

The softening lengths used in the simulations can affect the ability to resolve the very thinnest planar structures, which in turn can affect how much the density-based midplane varies as a function of azimuth. The Milky Way’s dense, star-forming gas disk is thought to have a scale height of about $40\,\text{pc}$, on the order of the cold gas softening length (Anderson et al., 2019). The thin stellar disk has a scale height of about $300\,\text{pc}$, $\sim 30$ times the stellar softening length (Jurić et al., 2008). We therefore expect that resolution effects are still affecting the scale heights of these components in the simulations, especially the cold gas. Indeed, the stellar scale heights of the simulated galaxies are equal to or larger than the Milky Way’s while the gas scale heights are significantly larger (although the proper basis comparison is less clear in the case of the gas; the quoted value for the Milky Way comes from studies of high-mass star-forming regions). The midplanes defined by gas and stars can be tilted with respect to one another as well, precluding extending the precision of the gas midplane definition to the stellar component.

Cosmological simulations of Milky Way-mass galaxies are not perfect representations of the true Milky Way in other ways as well, as discussed in Sanderson et al. (2018). The failure of cosmological simulations to exactly reproduce the Milky Way is not necessarily due to limitations of the numerical model. Candidate Milky Way-like galaxies are chosen solely on their mass and isolation, for which there are a wide variety of possible galaxies with qualitatively different properties. For example, the velocity structure of m12i is closer to M31’s than the Milky Way’s (S. Loebman et al., in preparation).

However, in this work we are most interested in the global properties of the potential, and specifically in deviations from axisymmetry. From this perspective, the simulated galaxies are actually more axisymmetric than we might expect of the Milky Way. While they have prominent spiral arms, none has as strong a bar as the Milky Way does at present day, and none has a nearby companion like the Large Magellanic Cloud. One of the three we consider (m12f) does have an ongoing interaction with a satellite galaxy similar in mass to Sagittarius, which has punched through the Galactic disk outside the solar circle, leaving behind some of its stars and inducing warping in the disk.

In this work, we take the galactocentric coordinate system described in Section 3 of Sanderson et al. (2018) as our fiducial coordinate system for each galaxy. In short, the center of the galaxy is found iteratively. The center of mass velocity is then determined by all star particles within $15\,\text{kpc}$ of this center. The galaxy is then rotated onto a principal axis frame determined by stars younger than $1\,\text{Gyr}$ inside of the fiducial solar radius $R_{0}=8.2\,\text{kpc}$, such that the disk plane is the $x$–$y$ plane.

3.2 Description of Milky Way-Sagittarius Interaction Simulation

In addition to the cosmological zoom-ins, we will also briefly consider results from a live N-body simulation of a Sagittarius-like encounter. This simulation offers us the ability to see how the midplane varies in a more controlled environment. The simulation is described by Laporte et al. (2018), but we briefly summarize the most relevant details here.

For the Milky Way, the dark halo is modeled as a Hernquist sphere of mass $10^{12}\,\text{M}_{\odot}$ and scale length of $52\,\text{kpc}$ (Hernquist, 1990), the disk is modeled as an exponential disk with a scale radius of $3.5\,\text{kpc}$, scale height $0.53\,\text{kpc}$, and mass $6\times 10^{10}\,\text{M}_{\odot}$, and the bulge as a Hernquist sphere of mass $10^{10}\,\text{M}_{\odot}$ and scale length $0.7\,\text{kpc}$. The Sagittarius dwarf is modeled with two components: a dark matter Hernquist sphere of mass $8\times 10^{10}\,\text{M}_{\odot}$ and scale length $8\,\text{kpc}$, and a stellar component modeled as a Hernquist sphere of mass $6.4\times 10^{8}\,\text{M}_{\odot}$ and scale length $0.85\,\text{kpc}$. All components are realized with distributions of live N-body particles; the Milky Way and Sagittarius are each initialized to be in equilibrium in isolation.

The mass resolution of the simulation is $2.6\times 10^{4}$, $1.2\times 10^{4}$, and $1.0\times 10^{4}\,\text{M}_{\odot}$ for the dark matter, disk, and bulge components, respectively. For the disk and bulge components, a softening length of $30\,\text{pc}$ is used whereas for the halo a softening length of $60\,\text{pc}$ is used. For Sagittarius, the softening length for the dark matter and the stars is $60$ and $40\,\text{pc}$, respectively.

The fiducial coordinate system for these N-body simulations is the rest frame of the aligned host galaxy at the beginning of the simulation.

3.3 The Local Midplane

Using the two sets of simulations, we determine the local midplane as a function of azimuth at the solar circle that an observer might measure if they were situated in each of these galaxies. Starting from the coordinate system described in the previous section, which is aligned so that the $z$-coordinate is approximately perpendicular to the disk plane at the solar circle, we place our imaginary observer at $z=0$ and a galactocentric cylindrical radius of $8.2\,\text{kpc}$ and vary the azimuth between $0<\phi<2\pi$. At equally spaced values of $\phi$ we then compute the median $z$ for stars within a cylinder of radius $0.5\,\text{kpc}$ and height $1\,\text{kpc}$ perpendicular to the fiducial disk and centered on it. We choose to use 50 bins in azimuth, sufficiently few that no cylinder shares stars with its neighbors. We then re-define the new midplane of the cylinder to be the median $z$, re-select stars, and iterate until the median $z$ value converges. We find that only $10$ iterations of this procedure are necessary for convergence. The resulting median $z$ is taken to be what our observer would measure as the local galactic midplane at each $\phi$.

This procedure assumes perfect density estimation, and therefore perfect corrections for extinction within the cylinder defining the “solar neighborhood.” Imperfect extinction correction is likely to increase the amplitude of the estimated fluctuations in $z$.

To account for the effect of particle noise, we bootstrap resample stars within a cylinder of height $2\,\text{kpc}$ and the same radius $1000$ times and determine the $1\,\sigma$ error bars by repeating the midplane determination with that reselection.

To allow for potential small inaccuracies in the determination of the original fiducial coordinate system, we also subtract the best fit curve of the form

[TABLE]

from the midplane as a function of azimuth to account for an overall tilt of the midplane (a simplified version of the strategy described in Anderson et al. 2019). For simulations (m12i, m12f, m12m), the best fit values are $A=(-170,45,8.8)\,\text{pc}$, $B=(38,-5.0,1.8)\,\text{deg}$, and $C=(-69,19,-18)\,\text{pc}$. For the assumed solar radius of $8.2\,\text{kpc}$, we can approximate the angle offset $\Delta\theta$ for the $z$-axis from the values of $A$. We find for the same simulations $\Delta\theta=(1.15,0.31,0.062)\,\text{deg}$. These angle offsets are consistent with the values given in Sanderson et al. (2018) for the difference between the $z$-axis as defined by the gas and stars.

Figure 6 shows the relative $z$ location of the inferred midplane an observer would determine as a function of azimuth for each galaxy, using their local solar neighborhood (the cylinder defined above). The $1\,\sigma$ error from the bootstrap procedure is shown as the dashed-line error bars. The middle $90\%$ of midplane values across the solar circle spans $(190,160,84)\,\text{pc}$ for these simulations. In two of the three cases the midplane therefore varies by more than $\pm 100\,\text{pc}$ depending on the azimuth along the solar circle; in the third (m12m, which has the thinnest “thin disk” of stars, but the largest stellar mass) the variation is closer to $\pm 50\,\text{pc}$.

We compute the same midplane variation in Figure 7, but for four succcessive timesteps of the live N-body simulation of a Sagittarius encounter (Laporte et al., 2018). Again we have subtracted a best fit curve of the form given in Equation (13), with the values at times (2.0, 4.0, 6.0, 6.9) Gyr being $A$ = (9.5, 2.5, -210, -390) pc, $B$ = (0.074, 0.039, -36, -57) deg, and $C$ = (2.6,-7.8, -65, -53) pc. The middle $90\%$ of midplane values across the solar circle spans (49, 62, 140, 120) pc.

These values for the midplane variation are consistent with the azimuthal midplane variations seen by Gómez et al. (2013). However, they only saw significant variations in their Heavy but not their Light Sagittarius model (virial masses of $10^{11}\,\text{M}_{\odot}$ and $\sim 3\times 10^{10}\,\text{M}_{\odot}$, respectively). The model we used (L2 from Laporte et al. (2018)) has a virial mass of $6\times 10^{10}\,\text{M}_{\odot}$, intermediate between their two models.

In the upper panels, we see that the midplane is relatively flat in the inner galaxy, but additional encounters drive strong midplane variation. In the lower left panel, we see a strong $m=2$ mode develop, consistent with the $R=8\,\text{kpc}$ panel of Figure 17 in Laporte et al. (2018) ($m=0$ and $m=1$ modes are stronger, but these are removed in our sine-curve subtraction). The lower right panel, which shows the galaxy at $t=6.9\,\text{Gyr}$ when some relaxation has occured, is qualitatively similar to the midplane variations we saw in the FIRE simulations (Figure 6), evidence that they are at least partially driven by mergers.

3.4 Velocity Variations

We also expect that the LSR should vary as a function of azimuth. We perform this calculation in Appendix D to estimate the components of the LSR as a function of azimuth, but performing a best-fit subtraction to correct for misalignment of the original coordinate system (as in the previous section) is more involved. Since we find that the variation in the LSR is less pronounced than for the midplane, and since offsets in velocity only contribute to second order to $\Delta J_{R}$ and $\Delta J_{z}$ when a star is at maximum amplitude in $R$ or $z$ (where the majority of the orbit is, see Section 2.2), we defer this calculation to future work.

4 Discussion

We have used high-resolution simulations to illustrate why we expect the local midplane defined by stellar density to vary with azimuth by up to $\pm 100\,\text{pc}$ as a natural consequence of the non-axisymmetry of the Galactic disk at small scales. While this is not in itself surprising or new, we have also argued that the discrepancy between our local midplane and that of distant stars introduces a systematic error in the $z$ component when converting from heliocentric to Galactocentric coordinates. This systematic error introduces a non-Gaussian error in the vertical action, $J_{z}$, when starting from the present-day positions and velocities of stars as measured by, e.g., Gaia.

These systematic errors are most important for stars on thin disk-like orbits, where they can be large enough to yield actions representative of orbits in the thick disk. This effect is entirely due to the extension of a local to a global coordinate system, and is separate from real diffusion in stellar integrals of motion caused by interactions with these same deviations from axisymmetry, such as resonant perturbation by spiral arms or scattering from molecular clouds (Sellwood, 2014).

4.1 Estimates of Milky Way Midplane Offsets

Systematic variations in $v_{z}$ and number density were first noted as asymmetries in the local velocity distribution towards the North and South Galactic Caps from the radial velocity surveys of the Sloan Digital Sky Surveys (Widrow et al., 2012) and RAdial Velocity Experiment (Williams et al., 2013). Subsequently, Carlin et al. (2013) pointed out suggestions of an oscillation in average vertical velocities of order $5\,\text{km}\,\text{s}^{-1}$ on roughly kiloparsec scales looking toward the Galactic anticenter.

Work by the Gaia collaboration confirmed these preliminary results on the velocity and spatial scales of oscillation with clear spatial maps made using DR2 data of median $v_{z}$ over a significant Galactic volume (Gaia Collaboration et al., 2018a; Friske & Schönrich, 2019), which can be explained with models of Sagittarius-like encounters (Gómez et al., 2013; Laporte et al., 2018, 2019b). We also see in the FIRE simulation that the vertical velocity variation as a function of azimuth is $\sim 5\textup{--}10\,\text{km}\,\text{s}^{-1}$ (Figure 10), consistent with these observations.

The vertical frequency of thin-disk and thick-disk are $\sim 0.09\,\text{Myr}^{-1}$ and $\sim 0.06\,\text{Myr}^{-1}$, respectively (Table 1). By dimensional analysis, and assuming a vertical velocity variation of $5\textup{--}10\,\text{km}\,\text{s}^{-1}$, we therefore expect the midplane offsets to be $\sim 57\textup{--}170\,\text{pc}$. We stress that this is a rough calculation.

Three-dimensional dust maps also offer a view into the expected variation of the stellar disk, since dust should trace regions of massive star formation. Figure 9 of Chen et al. (2019), Figure 1 of Leike & Enßlin (2019), and Figure 2 of Green et al. (2019) all show that the midplane varies by $\sim 10\degree$ at a distance of $\sim 0.75\,\text{kpc}$, corresponding to a physical vertical variation of $\sim 130\,\text{pc}$.

Already we see evidence in the data from velocities and dust maps for midplane offsets on the order of what we saw in both sets of simulations.

4.2 Uncertainties in the Solar Position and Velocity

Uncertainties in measurements of the position and velocity of the Sun relative to the Galactic center can also contribute to systematic error in the actions, since converting from heliocentric to Galactocentric coordinates relies on these measurements. Therefore, errors in their values will induce a systematic offset in the Galactocentric phase-space position of any observed star. Considerable effort has been placed on each of these measurements, but uncertainties remain, and detailed modeling across the disk — particularly for dynamically cold stars — may have to take them into account. Here we briefly review the current measurements of the four relevant quantities, their uncertainties, and the implications for the calculation of actions.

4.2.1 Galactic Center Position

First, one must define the center of the Galaxy. This is usually taken to be the location of the central supermassive black hole, Sagittarius A* (Sgr A*, e.g., Reid & Brunthaler, 2004). From stellar motions near Sgr A*, the distance from the Sun to Sgr A*, $R_{0}$, can be precisely measured (Gillessen et al., 2009; Gravity Collaboration et al., 2018). A recent measurement using near-infrared interferometry places $R_{0}$ at $8.178\pm 0.035$ kpc (Abuter et al., 2019), or a $0.4\%$ uncertainty.

However, the location of Sgr A* may not be equivalent to the location of the dynamical Galactic center, the point in three-dimensional space about which the stars in the solar neighborhood are orbiting. This assumption, although sensible and frequently made, has not yet been justified.

If the dynamical Galactic center is offset from Sgr A* by $100\,\text{pc}$, only a $1.2\%$ difference, then this induces a $\sim 15\%$ error in $J_{R}$ for the disk-like orbits we considered (see Section 2.4). The reason such a large error in $J_{R}$ can be generated by a small error in $R_{0}$ can be understood from the epicyclic approximation (Equation (11)), which states that $\Delta J_{R}/J_{R}=2\Delta R/A_{R}$. The fractional error in $J_{R}$ is related to the error in $R_{0}$ as a fraction of the radial amplitude of the orbit, which is much smaller than $R_{0}$ ($\sim 1.2\,\text{kpc}$ for thin-disk and thick-disk). This also implies the very precise $0.4\%$ measurement of $R_{0}$ still translates to a $\sim 6\%$ uncertainty in $J_{R}$.

The assumption that the dynamical Galactic center and Sgr A* are colocated is tested in any construction of a dynamical model where $R_{0}$ is a free parameter. For example, Küpper et al. (2015) measured $R_{0}$ while modeling the dynamics of the stream Palomar 5. Many other dynamical measurements of $R_{0}$ have been made (Bland-Hawthorn & Gerhard 2016 summarize many pre-Gaia results), but none have yet achieved a precision comparable to that of the distance to Sgr A*.

We did not consider in this work the effect of the angular position of the dynamical Galactic center being offset from Sgr A*.

4.2.2 Galactic Orientation

Second, one must define the angular orientation of the Galaxy. This was defined in 1958 by the IAU subcomission 33b (Blaauw et al., 1960) by defining the coordinates of the Galactic center in B1950 coordinates as (17:42:26.6, -28:55:00) and the North Galactic pole as (12:49:00, +27:24:00). These two quantities, together with $R_{0}$, define the orientation of the Galactic plane. However, there is growing evidence that the stellar midplane is tilted relative to this coordinate system (Goodman et al., 2014; Bland-Hawthorn & Gerhard, 2016), though not the H ii midplane (Anderson et al., 2019).

This tilt will contribute a systematic offset in $z$, with the exact magnitude depending on the position of the observed star. For instance, Goodman et al. (2014) quote a $\sim 0.4\degree$ tilt at $3.1\,\text{kpc}$, corresponding to a vertical height of $\sim 22\,\text{pc}$. This corresponds to a $37\%$ error in $J_{z}$ for thin-disk and a 5% error for thick-disk.

4.2.3 Solar Height

Third, one must define the Sun’s vertical distance from the Galactic midplane, which can be determined by identifying where the stellar density and velocities reach a maximum (effectively the median height of all disk stars). The solar height is usually taken to be $\sim 25\,\text{pc}$ (Chen et al., 2001), with a more recent measurement from Gaia DR2 placing it at $20.8\pm 0.3\,\text{pc}$ (Bennett & Bovy, 2019). Another strategy is to use the cold gas or H ii regions in the disk to define the Galactic midplane, leading to slightly different values (by $\sim 5\,\text{pc}$) for the Sun’s relative height (e.g., Anderson et al., 2019). A pre-Gaia review of these measurements is given by Bland-Hawthorn & Gerhard (2016). The discrepancy between gas-based and stellar-based determinations of the solar height is small, and thus only likely to dominate over intrinsic midplane variations on small scales, but will be relevant for detailed modeling of young and kinematically cold stars. For instance, it will induce a $\sim 10\%$ error in $J_{z}$ for an orbit with $z_{\text{max}}\sim 100\,\text{pc}$.

4.2.4 Local Standard of Rest

Finally, one must define the LSR, or mean velocity of stars near the Sun relative to the Galactic center (which is defined to have zero velocity), and the velocity of the Sun relative to the LSR. The radial ($U_{\odot}$) and vertical ($W_{\odot}$) components are computed by taking the mean motions of different stellar groups (e.g., Schönrich, 2012). The azimuthal component ($V_{\odot}$) is more difficult to measure, but can be modeled using the asymmetric drift relation (Binney & Tremaine, 2008). The values of the components of the LSR are usually taken from Schönrich et al. (2010). Their uncertainties should also lead to systematic errors in the actions, as given in Equations (9)–(10). For example, the value of the circular velocity is taken to be $\sim 220\,\text{km}\,\text{s}^{-1}$ (e.g., Bovy et al., 2012) with roughly $10\%$ uncertainty. We expect this to translate to at least a $10\%$ systematic error in $J_{\phi}$.

4.3 Orbit Integration

We have mainly been concerned with actions, since they provide a convenient way to quantify different types of orbits. However, all of our conclusions also apply to studies that simply rely on orbit integrations, since the two are equivalent. For instance, computing orbital properties of open or globular clusters (e.g., Cantat-Gaudin et al., 2016, 2018; Gaia Collaboration et al., 2018b) should ideally take the midplane variation into account. Orbit integrations of nearby systems over short time periods (e.g., Mamajek & Bell, 2014; Bailer-Jones et al., 2018) are unlikely to be impacted. It should also be unimportant for halo applications, e.g., in modeling of stellar streams (e.g., Bovy, 2014) or the substructure potentially responsible for the gap in GD1 (Bonaca et al., 2018).

5 Conclusions

Determining the orbital properties of stars is important for understanding the structure and evolution of the Galaxy. Actions have been argued to be excellent orbit labels. If the Galaxy can be well approximated as axisymmetric and 6D phase space positions can be measured accurately and precisely enough, then the computed actions are invariant with orbital phase. However, we have shown that the fact that the Galactic midplane is not constant across the disk presents a significant complication to computed actions actually being invariant. Our main conclusions are:

• •

Inaccuracy in the Galactocentric coordinate system induces orbital phase dependence in the actions calculated from the observed positions and velocities of stars (Figures 1 and 2). Since stars’ instantaneous phase-space positions are measured without prior knowledge of their orbital phases, this results in systematic error in the computed actions (Figure 4).

• •

Inaccuracy in the midplane location most severely affects computation of the vertical action $J_{z}$. A midplane offset of $\sim 15\,\text{pc}$ for a typical thin disk orbit results in a $25\%$ error in $J_{z}$, and even for a thick disk orbit a $120\,\text{pc}$ offset will result in the same size error. The fractional error is significantly less for halo orbits.

• •

The distribution of systematic errors in the actions induced by a coordinate system offset is highly non-Gaussian. The distribution is bimodal with neither mode at null. As a result, error propagation of coordinate system offsets is complex when considering actions, and is likely to significantly deform the action-space distribution function.

• •

Dynamical modeling across large regions of the disk, over which the midplane location varies by more than the limits discussed above, is susceptible to this type of systematic error, since the assumption that our local Galactic midplane is the global Galactic midplane is not true a priori. A violation of this assumption (by, e.g., intrinsic midplane variations) leads to a systematic error in $z$ which generates the large errors in actions summarized above.

• •

We show that such midplane variation is likely by measuring the local galactic midplane along the solar circle in three different high-resolution, zoom-in simulations of Milky Way mass galaxies from the FIRE collaboration, as well as a controlled simulation of the interaction of the Milky Way with Sagittarius. We found that the midplane varies as a function of azimuth at the solar circle by $60\textup{--}185\,\text{pc}$ in these simulations.

• •

Assuming a vertical velocity variation of the Milky Way of $\sim 5\textup{--}10\,\text{km}\,\text{s}^{-1}$, consistent with recent results from Gaia and our results from the FIRE simulations (Figure 10), we estimated that the corresponding midplane offsets are $\sim 60\textup{--}170\,\text{pc}$ by dimensional analysis using the vertical frequencies of disk-like orbits. This range of values is consistent with the variations seen in the simulations. Similar offsets are seen in three-dimensional dust maps.

• •

Inaccuracies in the parameters of the currently adopted Galactocentric coordinate system are likely important for some applications. In particular, it is imperative to test the assumption that the dynamical Galactic center is colocated with Sgr A*. We discuss how to do this in Section 4.2.

• •

This work underlines the importance of combining chemistry and dynamics. Since chemical tagging (Freeman & Bland-Hawthorn, 2002) is not subject to the same systematic errors discussed in this work, it should be used to confirm dynamical associations and to offset the effect of these systematic errors on the action-space distribution function.

• •

While in this work we have focused on systematic errors in action computation, all of our conclusions also extend to studies of stars that simply rely on orbit integration, since the computation of actions and orbit integrations are equivalent.

Our main point is that the local midplane varies between different points in the Galaxy, and that this variation can lead to significant systematic errors in the computation of actions under the assumption of a global axisymmetric potential. Current observations from Gaia should soon permit a measurement of the real azimuthal dependence of the midplane location. For some applications, such as those using actions as labels to group stars on similar orbits, using such a measurement to shift stars to a consistent midplane height as a function of azimuth before using a global axisymmetric approximation to the potential may be sufficient, although this ignores the dynamical implications of shifts in the midplane height (which result from fluctuations in the local density). However, for other applications, such as the study of action diffusion, a more extensive perturbative approach is likely needed. We plan to explore the mitigation of these effects in future work.

Supplementary Material

All code used in this work is available at https://github.com/gusbeane/actions_systematic.

We thank the anonymous referee for providing helpful comments. We would like to thank Megan Bedell, Robert A. Benjamin, Maria Bergemann, Joss Bland-Hawthorn, Tobias Buck, Elena D’Onghia, Benoit Famaey, and Adrian Price-Whelan for helpful discussions. A.B. would like to thank Todd Phillips for helpful discussions. This work uses data hosted by the Flatiron Institute’s FIRE data hub. The Flatiron Institute is supported by the Simons Foundation. This project was developed in part at the 2019 Santa Barbara Gaia Sprint, hosted by the Kavli Institute for Theoretical Physics (KITP) at the University of California, Santa Barbara. This research was supported in part at KITP by the Heising-Simons Foundation and the National Science Foundation (NSF) under Grant No. PHY-1748958. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF Grant No. OCI-1053575. A.B. was supported in part by the Roy & Diana Vagelos Program in the Molecular Life Sciences and the Roy & Diana Vagelos Challenge Award. K.V.J.’s contributions were supported in part by the NSF under Grant No. AST-1715582. M.-M.M.L. was partly supported by the NSF under Grant No. AST-1815461.

Appendix A Orbits

We plot the three orbits considered throughout the work (Table 1) in Figure 8.

Appendix B $\Delta R$-$\Delta x$ Relation

We considered the effect on actions of an inaccuracy in the distance from the Sun to the Galactic center, which introduces an offset in the $x$ coordinate, $\Delta x$, of each star when converting to a Galactocentric coordinate system. In observations of nearby stars, we have that $\Delta x\sim\Delta R$. However, for the experiment we performed in Section 2.4 we considered observations of a star throughout its entire orbit. Therefore, we must average $\Delta R$ over the course of the orbit. We derive this relation now.

An offset $\Delta x$ results in an erroneous radius $R_{\text{err}}$ related by the formula,

[TABLE]

Keeping only terms to first order in $\Delta x$, we have that,

[TABLE]

Averaging over the circle, we therefore have that,

[TABLE]

Appendix C $J_{R}$ and $J_{\phi}$ Distributions

In Figure 9 we plot the distribution of $J_{R}$ as a function of orbital phase induced by an offset in $x$ and $z$ and the distribution of $J_{\phi}$ for an offset in $x$. We plot the distributions for thick-disk (upper panels) and thin-disk (lower panels). We find that the $J_{R}$ distribution induced by an offset in $x$ more closely resembles a Gaussian distribution, while the $J_{R}$ distribution induced by an offset in $z$ and the $J_{\phi}$ distribution induced by an offset in $x$ are both similar to the $J_{z}$ distribution induced by an offset in $z$ (see Figure 3).

Appendix D LSR Variations

We consider the variations of the LSR as a function of azimuth at the fiducial solar circle ($R_{0}=8.2\,\text{kpc}$) in Figure 10. At each azimuth, $\phi$, we take the median velocity in cylindrical coordinates of all stars within $200\,\text{pc}$ of the position, following Sanderson et al. (2018). No best-fit subtraction was performed as in Figure 6.