Action-based Dynamical Modeling for the Milky Way Disk: The Influence of Spiral Arms

TL;DR

This study demonstrates that the RoadMapping dynamical modeling method can reliably infer the Milky Way's gravitational potential from star data, even in the presence of spiral arms, using simulated galaxy data.

Contribution

We show that RoadMapping remains robust in modeling the MW potential with spiral arms, providing reliable constraints from various survey volumes and informing future Gaia data analyses.

Findings

Potential constraints are robust even with spiral arms.

Models from large survey volumes can be extrapolated reliably.

Biases occur when data is dominated by spiral arms.

Abstract

RoadMapping is a dynamical modeling machinery developed to constrain the Milky Way's (MW) gravitational potential by simultaneously fitting an axisymmetric parametrized potential and an action-based orbit distribution function (DF) to discrete 6D phase-space measurements of stars in the Galactic disk. In this work we demonstrate RoadMapping's robustness in the presence of spiral arms by modeling data drawn from an N-body simulation snapshot of a disk-dominated galaxy of MW mass with strong spiral arms (but no bar), exploring survey volumes with radii 500pc<=r_max<=5kpc. The potential constraints are very robust, even though we use a simple action-based DF, the quasi-isothermal DF (qDF). The best-fit RoadMapping model always recovers the correct gravitational forces where most of the stars that entered the analysis are located, even for small volumes. For data from large survey volumes, RoadMapping finds axisymmetric models that average well over the spiral arms. Unsurprisingly, the models are slightly biased by the excess of stars in the spiral arms. Gravitational potential models derived from survey volumes with at least r_max=3kpc can be reliably extrapolated to larger volumes. However, a large radial survey extent, r_max~5kpc, is needed to correctly recover the halo scale length. In general, the recovery and extrapolability of potentials inferred from data sets which were drawn from inter-arm regions appear to be better than those of data sets drawn from spiral arms. Our analysis implies that building axisymmetric models for the Galaxy with upcoming Gaia data will lead to sensible and robust approximations of the MW's potential.

Tables

Table 1. Table 1 Vantage points within the galaxy simulation snapshot around which we center survey volumes of radius r max subscript 𝑟 max r_{\text{max}} .

name	position	$R_0$ [kpc]	${\varphi }_0$ [degrees]	legend
S8	on spiral arm	8	5
I8	in inter-arm region	8	-15
S5	on spiral arm	5	60
I5	in inter-arm region	5	0

Table 2. Table 2 Best fit parameters of the DEHH-Pot .

circular velocity	$v_{\text{circ}}\left(R_{\text{☼}}\right)$	$222{\text{km s}}^{-1}$
disk scale length	$R_{\text{s}}$	$2.5\text{kpc}$
disk scale height	$z_{\text{s}}$	$0.17\text{kpc}$
halo fraction	$f_{\text{halo}}$	$0.54$
halo scale length	$a_{\text{halo}}$	$29\text{kpc}$
bulge mass	$M_{\text{bulge}}$	$0.95\times {10}^{10}{\text{M}}_{\odot }$
bulge scale length	$a_{\text{bulge}}$	$0.25\text{kpc}$

Table 3. Table 3 Best fit MNHH-Pot and qDF parameters as recovered from the RoadMapping analysis of a survey volume with r max = 4 ​ kpc subscript 𝑟 max 4 kpc r_{\text{max}}=4~{}\text{kpc} centered on a spiral arm at R 0 = 8 ​ kpc subscript 𝑅 0 8 kpc R_{0}=8~{}\text{kpc} (position S8 ).

circular velocity at $R_{\text{☼}}=8\text{kpc}$	$v_{\text{circ}}\left(R_{\text{☼}}\right)$	$\left(223.0\pm 0.1\right){\text{km s}}^{-1}$
Miyamoto-Nagai disk scale length	$a_{\text{disk}}$	$\left({3.62}_{-0.05}^{+0.06}\right)\text{kpc}$
Miyamoto-Nagai disk scale height	$b_{\text{disk}}$	$\left(0.26\pm 0.02\right)\text{kpc}$
halo fraction at $R_{\text{☼}}=8\text{kpc}$	$f_{\text{halo}}$	$\left(0.53\pm 0.02\right)$
halo scale length	$a_{\text{halo}}$	$\left(21\pm 2\right)\text{kpc}$
bulge mass	$M_{\text{bulge}}$	$0.95\times {10}^{10}{\text{M}}_{\odot }$ (fixed)
bulge scale length	$a_{\text{bulge}}$	$0.25\text{kpc}$ (fixed)
qDF tracer scale length	$h_R$	$\left({3.34}_{-0.04}^{+0.05}\right)\text{kpc}$
qDF radial velocity dispersion	${\sigma }_{R,0}$	$\left(15.91\pm 0.08\right){\text{km s}}^{-1}$
qDF vertical velocity dispersion	${\sigma }_{z,0}$	$\left({14.0}_{-0.1}^{+0.2}\right){\text{km s}}^{-1}$
qDF radial velocity dispersion scale length	$h_{\sigma ,R}$	$\left(4.6\pm 0.5\right)\text{kpc}$
qDF vertical velocity dispersion scale length	$h_{\sigma ,z}$	$\left({5.65}_{-0.07}^{+0.06}\right)\text{kpc}$

Equations

$${\rho }_{\text{dm}}(r)=\frac{M_{\text{dm}}}{2\pi }\frac{a_{\text{dm}}}{r(r+a_{\text{dm}}{)}^3}$$

$${\rho }_{\ast }(R,z)=\frac{M_{\ast }}{4\pi z_{\text{s,init}}{R}_{\text{s}}^2}{\text{sech}}^2\left(\frac{z}{z_{\text{s,init}}}\right)\exp \left(-\frac{R}{R_{\text{s}}}\right),$$

$${\Delta }_{\text{Spiral}}(x_k,y_k)\equiv \frac{{\Sigma }_{\text{1.5kpc,disk},T}(x_k,y_k)}{{\Sigma }_{\text{1.5kpc,disk},S}(x_k,y_k)}-1$$

$${\Sigma }_{\text{1.5kpc,disk},\alpha }(x_k,y_k)\equiv {\int }_{-1.5\text{kpc}}^{1.5\text{kpc}}{\rho }_{\text{disk},\alpha }(x_k,y_k,z)\mathrm{d}z.$$

$$⟨{\Delta }_{\text{Spiral}}⟩$$

$${\sigma }_{\Delta \text{Spiral}}$$

$$(x_k-x_c{)}^2+(y_k-y_c{)}^2\le r_{\text{max}}^2.$$

$$\text{SF}(\mathbf{x})\equiv \{\begin{array}{ll}1& \text{mbox}if\left|\mathbf{x}-{\mathbf{x}}_0\right|\le r_{\text{max}}\\ 0& \text{mbox}otherwise\end{array},$$

$${\mathbf{J}}_i\equiv \mathbf{J}[{\mathbf{x}}_i,{\mathbf{v}}_i\mid p_{\Phi }],$$

$$\text{DF}({\mathbf{J}}_i\mid p_{\text{DF}})\equiv \text{DF}({\mathbf{x}}_i,{\mathbf{v}}_i\mid p_{\Phi },p_{\text{DF}}).$$

$${\mathcal{L}}_i$$

$${\mathcal{L}}_i⟶\max \left({\mathcal{L}}_i,ϵ\times \text{median}(\mathcal{L})\right),$$

$$n(R_g\mid p_{\text{DF}})$$

$${\sigma }_R(R_g\mid p_{\text{DF}})$$

$${\sigma }_z(R_g\mid p_{\text{DF}})$$

$$p_{\text{DF}}\equiv \left\{\ln h_R,\ln {\sigma }_{R,0},\ln {\sigma }_{z,0},\ln h_{\sigma ,R},\ln h_{\sigma ,z}\right\}.$$

$${\Phi }_{\text{disk}}(R,z)=-\frac{GM}{\sqrt{R^2+(a_{\text{disk}}+\sqrt{z^2+b_{\text{disk}}^2}{)}^2}},$$

$$f_{\text{halo}}\equiv {\frac{F_{R,\text{halo}}}{F_{R,\text{disk}}+F_{R,\text{halo}}}|}_{\stackrel{R=R_{\text{sun}}}{z=0}}.$$

$$v_{\text{circ}}(R_{\text{sun}}=8\text{kpc})\equiv {\sqrt{R\frac{\partial \Phi }{\partial R}}|}_{\stackrel{R=R_{\text{sun}}}{z=0}}.$$

$$p_{\Phi }\equiv \left\{v_{\text{circ}}(R_{\text{sun}}),a_{\text{disk}},b_{\text{disk}},a_{\text{halo}},f_{\text{halo}}\right\}.$$

$$F_{R,\text{typ}}(R)$$

$$F_{z,\text{typ}}(R)$$

$$\Delta F_R({\mathbf{x}}_{\ast ,i})\equiv \frac{|F_{R,M}(R_i,z_i)|-|F_{R,T}(x_i,y_i,z_i)|}{|F_{R,\text{typ}}(R_i)|}$$

$$\Delta F_z({\mathbf{x}}_{\ast ,i})\equiv \frac{|F_{z,M}(R_i,z_i)|-|F_{z,T}(x_i,y_i,z_i)|}{|F_{z,\text{typ}}(R_i)|}$$

$$\Delta F_R({\mathbf{x}}_{g,j})\equiv \frac{|F_{R,M}(R_j,z_j)|-|F_{R,T}(x_j,y_j,z_j)|}{|F_{R,\text{typ}}(R_j)|}$$

$$\Delta F_z({\mathbf{x}}_{g,j})\equiv \frac{|F_{z,M}(R_j,z_j)|-|F_{z,T}(x_j,y_j,z_j)|}{|F_{z,\text{typ}}(R_j)|},$$

Full text

††software: GADGET-3 (Springel, 2005), galpy (Bovy, 2015), emcee (Foreman-Mackey et al., 2013), Matplotlib (Hunter, 2007)

I. Introduction

An important basis for learning more about the Milky Way’s (MW) overall gravitational potential and orbit distribution function (DF) is to find the “best possible” axisymmetric model for the Galaxy. Given such a model the identification and characterization of non-axisymmetries like spiral arms or stellar streams in stellar phase-space (and chemical abundance) data would then become more straightforward.

Several approaches to constrain an axisymmetric potential and/or orbit DF have recently been put forward: Bovy & Rix (2013) and Piffl et al. (2014) fitted potential and DF simultaneously to stellar kinematics in the disk and got precise constraints on the overall potential; Sanders & Binney (2015) and Das & Binney (2016) investigated extended DFs for the disk and halo respectively (given a fiducial potential), that included the metallicity of each star, in addition to the distribution in orbit space.

In this work we will continue our investigation of the RoadMapping approach (“Recovery of the Orbit Action Distribution of Mono-Abundance Populations and Potential INference for our Galaxy”). The first application of RoadMapping was performed by Bovy & Rix (2013). Trick et al. (2016), hereafter Paper I, subsequently performed a detailed analysis of the strengths and limitations of the approach. RoadMapping presumes that simple stellar populations in the MW disk—be it mono-abundance populations (MAPs), i.e., stars with the same $[\mathrm{Fe}/\mathrm{H}]$ and $[\alpha/\mathrm{Fe}]$ (Bovy et al., 2012a, c, b, 2016), or maybe also mono-age populations (Bird et al., 2013; Martig et al., 2014, 2016; Minchev et al., 2014; Ness et al., 2016; Minchev et al., 2017)—follow simple orbit DFs, like, e.g., the quasi-isothermal DF (qDF) by Binney & McMillan (2011). That MAPs in the MW are well described by the qDF was first shown by Ting et al. (2013). The qDF is expressed in terms of the orbital actions $\boldsymbol{J}=(J_{R},J_{\phi}=L_{z},J_{z})$, which are integrals of motion, quantify the amount of the orbit’s oscillation in each of the coordinate directions $(R,\phi,z)$, and are therefore excellent orbit labels. Given an assumed gravitational potential one can calculate the orbital actions from the stars’ current phase-space positions $(\boldsymbol{x},\boldsymbol{v})$ (Binney 2012; Sanders & Binney 2016; see also Bovy (2014) for a general method to compute $(\boldsymbol{x},\boldsymbol{v})\longrightarrow(\boldsymbol{J},\boldsymbol{\theta})$). Only if this assumed gravitational potential is close to the true potential, the action distribution of the stellar MAP in question will follow an orbit DF of qDF-shape. This is the idea on which RoadMapping builds, and which allows us to simultaneously fit potential and orbit DF to observations.

Bovy & Rix (2013) employed this approach to measure the Milky Way’s surface density profile within $1.1~{}\text{kpc}$ using 43 MAPs in the Galactic disk from the SDSS/SEGUE survey (Yanny et al., 2009). To avoid spiral arm effects, they did not use in-plane motions. Their potential model had only two free parameters (disk scale length and relative halo-to-disk contribution to the radial force at the solar radius). To account for missing model flexibility they constrained the surface density for each MAP only at one best radius. The profile they derived in this fashion had a scale length of $R_{\text{s}}=2.5~{}\text{kpc}$ and was—in the regime $R>6.6~{}\text{kpc}$—later confirmed by Piffl et al. (2014) using a different action-based procedure.

Given the success of this first application and in anticipation of the upcoming data releases from Gaia in 2016-2022 (Eyer et al., 2013), Paper I improved the RoadMapping machinery and studied its strengths and breakdowns in detail, by investigating a large suite of mock data sets. Under the prerequisite of axisymmetric data and model, we found that RoadMapping’s modeling success is stable against minor misjudgments of DF or selection function, and that—if the true potential is not contained in the proposed family of model potentials—one can still find a good fit that returns the correct forces, given the limitations of the model. Paper I also found that measurement uncertainties of the order of those by the final Gaia data release should be good enough (within $3~{}\text{kpc}$ from the Sun) to allow for precise and unbiased modeling results.

The MW is, however, not axisymmetric. The bulge contains a strong bar (Liszt & Burton, 1980; Blitz & Spergel, 1991; Binney et al., 1991, 1997; Hammersley et al., 2000; Wegg & Gerhard, 2013) and the disk itself is threaded by spiral arms (Oort et al., 1958; Georgelin & Georgelin, 1976; Churchwell et al., 2009; Reid et al., 2014), and (ring-like) overdensities (Newberg et al., 2002; Jurić et al., 2008; Xu et al., 2015), which induce non-circular motions and asymmetries in stellar number counts. There is also kinematic evidence in the disk for moving groups (Dehnen, 1998; Famaey et al., 2005; Bovy et al., 2009; Bovy & Hogg, 2010) and streaming motions (in 21 cm or velocities) (Bovy et al., 2015; Williams et al., 2013; Siebert et al., 2012), both of which are likely caused by non-axisymmetric perturbations to the gravitational potential.

As RoadMapping and related approaches can only build axisymmetric models, this is an important breakdown of modeling assumptions which was not investigated in Paper I. In this paper we want to understand in which respects RoadMapping will still give reliable constraints on the MW’s gravitational potential in the presence of spiral arms.

Our investigation makes use of an $N$-body simulation snapshot of a spiral galaxy with strong spiral arms presented in D’Onghia et al. (2013) (see their Figure 8, top left panel). From this snapshot we draw mock data in regions with different spiral arm strengths. We then apply the RoadMapping machinery to these data sets and test how well we recover the local and overall gravitational potential.

In Paper I we confirmed and tested separately the robustness of RoadMapping in the case that the data came from a different model family—for either the potential or DF—than assumed in the dynamical modeling. What would happen if both potential and DF model families were slightly wrong at the same time? The set-up of this study will automatically cover this important test case. The potential and orbit DF model that we are using were picked as a pragmatic compromise between (i) being a reasonable choice given the initial axisymmetric set-up of the galaxy simulation, and (ii) because of their simplicity, computational advantages, and—in case of the qDF—because that’s what we are planning to use in the MW. Given that the simulation has evolved away from its axisymmetric beginnings, we expect our chosen model to be reasonable, but not particularly well-suited to model this galaxy.

Spiral arms introduce another—but possibly minor—breakdown of the modeling assumptions: In a non-axisymmetric gravitational potential the three actions will not be strict integrals of motions anymore (Binney & Tremaine, 2008; Minchev et al., 2011, 2012; Solway et al., 2012; Grand et al., 2012; Vera-Ciro & D’Onghia, 2016). It will be interesting to see if the action-based DF in RoadMapping modeling is still informative, even if it only uses the approximate actions estimated in an axisymmetric potential.

In the MW we expect the central bar to introduce additional non-axisymmetries in the Galactic disk—but presumably not stronger ones than the spiral arms. As the galaxy simulation in this work does not have a central bar, we do not investigate specific bar effects here.

Though non-axisymmetry could be a severe problem for RoadMapping, we show in this paper that RoadMapping potential estimates are still surprisingly accurate, which makes us optimistic that they will also be so for the MW.

This paper is organized as follows. Section II describes the $N$-body simulation snapshot of a spiral galaxy that we are going to model in this study, explains how we extract 6D stellar phase-space data from it, and how we quantify the spiral arm strength. There we also review similarities and differences between the simulation in this work and what we know about the MW. Section III summarizes the RoadMapping dynamical modeling framework, and introduces the DF and potential model that we will fit to the data. Section IV is dedicated to presenting the results: In Section IV.1 we discuss in detail the RoadMapping modeling results derived from one data set within a survey volume with radius $r_{\text{max}}=4~{}\text{kpc}$ around the Sun. Section IV.2 then investigates a whole suite of RoadMapping analyses, corresponding to survey volumes of different sizes and different positions within the galaxy and with respect to the spiral arms. In Section V we discuss the results and give an outlook to the application of RoadMapping to Gaia data. We conclude in Section VI.

II. Data from a galaxy simulation

RoadMapping requires 6D phase-space coordinates $(\boldsymbol{x}_{i},\boldsymbol{v}_{i})$ for a large set of stars that move independently in a collisionless galactic potential. If we want to test RoadMapping on a simulated galaxy, it is most convenient to apply it to an $N$-body simulation with a huge number of low-mass “star” particles. In that way, we can directly take the positions and the motions of individual particles as independent tracers of the potential, just as with the stars in the MW, without having to use an error-prone prescription to turn a single particle into many stars. The high-resolution simulations with its millions of particles by D’Onghia et al. (2013) satisfy this requirement.

II.1. Description of the galaxy simulation snapshot

The high-resolution $N$-body simulation snapshot of a disk galaxy by D’Onghia et al. (2013), which we use in this work, was carried out with the GADGET-3 code, and set up in the manner described in Springel et al. (2005). In this simulation, overdensities with properties similar to giant molecular clouds induced four prominent spiral arms—and therefore a non-axisymmetric sub-structure—via the swing amplification mechanism. This galaxy simulation was also investigated by D’Onghia et al. (2013) (their Figure 8, top left panel) and D’Onghia (2015) (their Figure 2, top left panel). For details see D’Onghia et al. (2013), here we summarize the essential characteristics.

The simulation has a gravitationally evolving stellar disk within a static/rigid analytic dark matter (DM) halo.

The analytic halo follows a Hernquist (1990) profile

[TABLE]

with total halo mass $M_{\text{dm}}=9.5\times 10^{11}~{}\text{M}_{\odot}$ and scale length $a_{\text{dm}}=29~{}\text{kpc}$.

The disk consists of $10^{8}$ “disk star” particles, each having a mass of $\sim 370~{}\text{M}_{\odot}$, and 1000 “giant molecular cloud” particles with mass $\sim 9.5\times 10^{5}~{}\text{M}_{\odot}$. The initial vertical mass distribution of the stars in the disk is specified by the profile of an isothermal sheet with a radially constant scale height $z_{\text{s,init}}$, i.e.,

[TABLE]

with total disk mass $M_{*}=0.04M_{\text{dm}}=3.8\times 10^{10}~{}\text{M}_{\odot}$. The scale-length $R_{\text{s}}$ is assumed to be $2.5~{}\text{kpc}$ and $z_{\text{s,init}}=0.1R_{\text{s}}$. In this model the disk fraction within $2.2$ scale lengths is 50% of the total mass, leading to a formation of approximately four arms (D’Onghia, 2015) (see Figure 2(b)).

The bulge consists of $10^{7}$ “bulge star” particles with mass $\sim 950~{}\text{M}_{\odot}$ and they are distributed following a spherical Hernquist profile analogous to Equation (1), with total mass $M_{\text{bulge}}=0.01M_{\text{dm}}=9.5\times 10^{9}~{}\text{M}_{\odot}$ and scale length $a_{\text{bulge}}=0.1R_{\text{s}}=0.25~{}\text{kpc}$.

The initial velocity setup of the “disk star” particles assumes for simplicity Gaussian velocity dispersion profiles (Springel et al., 2005).

The simulation snapshot which we are using in this work has evolved from these initial conditions in isolation for $\sim 250~{}\text{Myr}$, which corresponds to approximately one orbital period at $R\sim 8~{}\text{kpc}$. The mass density of simulation particles (without the DM halo) at this snapshot time is shown in Figure 1. The “molecular cloud perturbers”, which caused the formation of the four pronounced spiral arms, can be seen in Figure 1 as small overdensities in the disk. The spherical bulge and very flattened disk are shown in Figure 1(b).

We have confirmed that the gravitational center of the particles corresponds to the coordinate origin.

II.2. Data selection and survey volume

The selection function of all-sky surveys like Gaia, that are only limited by the brightness of the tracers, are contiguous and—when ignoring anisotropic effects like dust obscuration—spherical in shape. For simplicity we will use spherical survey volumes centered on different vantage points, and with sharp edges at a distance $r_{\text{max}}$ around it (see Equation (8)), which corresponds to a magnitude cut for stellar tracers all having the same luminosity.

Figure 1 illustrates the different survey volume positions analyzed in this study. We selected volumes with $r_{\text{max}}=[0.5,1,2,3,4,5]~{}\text{kpc}$ centered on a spiral arm (S) and on an inter-arm region (I) at both the equivalent of the solar radius, $R_{0}=8~{}\text{kpc}$, and at $R_{0}=5~{}\text{kpc}$, where the disk strongly dominates (see Figure 3), and the spiral arms are more pronounced than at $R_{0}=8~{}\text{kpc}$ (see Figure 2). The exact positions of the vantage points S8, I8, S5, and I5 are summarized in Table 1.

From within each volume we drew $N_{*}=20,000$ random “disk star” particles, and used their phase-space positions $(\boldsymbol{x}_{i},\boldsymbol{v}_{i})$ within the simulated galaxy’s rest-frame as data.

To make the data sample more realistic, one would actually have to add measurement uncertainties, especially to the distances from the survey volume’s central vantage point and the proper motions measured from there. We decided not to include measurement uncertainties: Firstly, their effect on RoadMapping modeling has been already investigated in Paper I, and we found that the measurement uncertainties of the last data release of Gaia should be small enough to not significantly disturb the modeling. Secondly, in this study we want to isolate and investigate the deviations of the data from axisymmetry and the assumed potential and DF model independently of other effects.

II.3. True symmetrized potential

For a galaxy with pronounced spiral arms, an axisymmetric model matter distribution per se cannot reproduce the true matter distribution globally. We therefore obtain an “overall best fit symmetrized” potential model from the distribution of particles to be able (i) to quantify the non-axisymmetries in the simulation snapshot better and (ii) to compare how close our axisymmetric RoadMapping results can get to it.

We derive this model by fitting axisymmetric analytical functions to the density distribution of each of the galaxy components’ particles. The bulge and halo follow Hernquist profiles by construction (see Section II.1).

The disk in this simulation snapshot deviates from its initial conditions in Equation (2): after $250~{}\text{Myr}$ pronounced spiral arms have formed, also causing some in-plane heating. Except of an overdensity around $R\sim 6~{}\text{kpc}$ (see Figure 6(d) in Section IV.1.2), the overall radial surface density profile (i.e., the azimuthal average) did not change by much, and appears smooth and exponential. We therefore chose a double exponential disk model to fit the particle distribution in the disk. The fit assumes the total disk mass to be known and be equal to the total mass of all disk particles. The best fit double exponential disk profile is found by maximizing the likelihood for all disk particles to be drawn from this axisymmetric density profile. In this way the fit does not depend on binning choices and is driven by the number and location of the stars—analogous to our RoadMapping procedure.

The best fit parameters for this reference potential, to which we will refer as the DEHH-Pot (Double-Exponential disk + Hernquist halo + Hernquist bulge) in the remainder of this work, are given in Table 2. As can be seen in Figure 6 in Section IV.1.2 below, the DEHH-Pot fits the overall true density distribution very well, with the exception of $z\sim 0$ where the particle distribution is not as cuspy as the exponential disk. Figure 3 shows the circular velocity curve of the DEHH-Pot, and its decomposition into disk, halo, and bulge contribution. The disk clearly dominates between $R\sim 2~{}\text{kpc}$ and $R\sim 7~{}\text{kpc}$.

II.4. Quantifying the strength of spiral arms

Depending on the size and position of the survey volume, spiral arms and inter-arm regions dominate the stellar distribution within the volume to different degrees. To quantify the strength of the spiral arms, we introduce the quantity

[TABLE]

where $\Sigma_{\text{1.5kpc,disk},\alpha}$ is the true surface density of the disk component of the simulation snapshot ($\alpha=T$ for “true”), or of the symmetrized snapshot model DEHH-Pot in Section II.3 ($\alpha=S$ for “symmetrized”),

[TABLE]

$(x_{k},y_{k})$ are the coordinates of regular grid points with spacing $\delta=0.25~{}\text{kpc}$.111We average the particle surface density of the true simulation potential over area element sizes of $\delta\times\delta$ around $(x_{k},y_{k})$, when calculating $\Sigma_{\text{1.5kpc,disk},T}(x_{k},y_{k})$. $(x_{c}=R_{0,c}\cdot\cos\phi_{0,c},y_{c}=R_{0,c}\cdot\sin\phi_{0,c},z_{c}=0)$ is the position of the survey volume’s center within the simulation, with $c\in\left\{\texttt{S8},\texttt{I8},\texttt{S5},\texttt{I5}\right\}$ and $(R_{0,c},\phi_{0,c})$ given in Table 1. We consider all $n\simeq\pi r_{\text{max}}^{2}/\delta^{2}$ values of $\Delta_{\text{Spiral}}(x_{k},y_{k})$ inside a given survey volume of radius $r_{\text{max}}$ around position $c$ and calculate the mean and standard deviation,

[TABLE]

These quantities tell us if and how much a spiral arm or an inter-arm region dominates the survey volume ($\langle\Delta_{\text{Spiral}}\rangle>0$ for spiral arms, $\langle\Delta_{\text{Spiral}}\rangle<0$ for inter-arm regions) and how large the relative contrast between spiral arms and inter-arm regions is ($\sigma_{\Delta\text{Spiral}}$). For example, volumes will have a smaller relative spiral contrast $\sigma_{\Delta\text{Spiral}}$, if they are either small and sitting completely within an inter-arm region, or if they are large volumes that contain—in addition to some spiral arms and depleted inter-arm regions—large areas of unperturbed disk.

Figure 4 shows $\Delta_{\text{Spiral}}$ as function of $(x_{k},y_{k})$, a histogram over all $\Delta_{\text{Spiral},k}$ within the galaxy, and $\langle\Delta_{\text{Spiral}}\rangle$ and $\sigma_{\Delta\text{Spiral}}$ calculated for all test survey volumes in this work (see Section II.2), depending on position and size.222When considering the whole galaxy or a large survey volume, $\langle\Delta_{\text{Spiral}}\rangle$ in Figures 4(b) and 4(c) is not exactly at [math], but slightly larger ($<0.05$). We account this small bias to the different functional forms of the DEHH-Pot and the initial axisymmetric disk in the simulation, Equation (2). This bias will however not affect our results.

II.5. Comparison of the $N$-body simulation to the MW

The $N$-body simulation in this work is not a perfect match to the MW. But to be suitable for this study it is most important that it satisfies our main requirements: it is a disk-dominated spiral galaxy with a very high number of star particles and strong spiral arms. To set the context, we discuss in the following the differences between the simulation at hand and what we know about the MW.

Bar.

The MW has a central bar; our simulation does not. The Galactic bar can introduce non-axisymmetries around the co-rotation and Lindblad resonances (Dehnen 2000; Fux 2001; Quillen 2003; Quillen & Minchev 2005; Sellwood 2010; Minchev & Famaey 2010; Monari et al. 2017; see also review by Gerhard 2011). Bar effects are consequently not part of this study and remain a possible source of uncertainty in RoadMapping. We suspect however that its influence is modest outside of the bar region, except close to the outer Lindblad resonance.

Mass components.

The MW’s dark halo is estimated to have a mass $M_{\text{dm},200,\text{MW}}\approx 10^{12}~{}\text{M}_{\odot}$ (Bland-Hawthorn & Gerhard, 2016); the simulation’s halo is only slightly less massive with $M_{\text{dm},200}\approx 7\times 10^{11}~{}\text{M}_{\odot}$. The stellar bulge mass of $M_{\text{bulge}}=9.5\times 10^{9}~{}\text{M}_{\odot}$ in this simulation is a bit smaller than the estimated stellar mass of the MW bulge with $M_{\text{bulge,MW}}=(1.4-1.7)\times 10^{10}~{}\text{M}_{\odot}$ (Portail et al., 2015). The total stellar mass of the MW is estimated to be $M_{\text{stars,MW}}=(5\pm 1)\times 10^{10}~{}\text{M}_{\odot}$ (Bland-Hawthorn & Gerhard, 2016; Bovy & Rix, 2013), consistent with the total baryonic mass in the simulation of $M_{\text{stars}}=4.75\times 10^{10}~{}\text{M}_{\odot}$. The fraction of bulge mass to total stellar mass of the MW is $M_{\text{bulge,MW}}/M_{\text{stars,MW}}=0.3\pm 0.06$; in our simulation it is $0.25$.

Disk.

The best estimate for the MW’s thin disk scale length from combining several measurements in the literature is $R_{\text{s}}=2.6\pm 0.5~{}\text{kpc}$ (Bland-Hawthorn & Gerhard, 2016). Bovy & Rix (2013), for example, found a stellar disk scale length of $R_{\text{s}}=2.15\pm 0.14~{}\text{kpc}$. The disk of this $N$-body simulation has a similar scale length, $R_{\text{s}}=2.5~{}\text{kpc}$. The stellar disk is however thinner than in the MW. Jurić et al. (2008) found scale heights $300~{}\text{pc}$ and $900~{}\text{pc}$ (with $20\%$ uncertainty) for the thin and thick disk of the MW, respectively; Bovy et al. (2012c), who considered the disk as a superposition of many exponential MAPs, measured scale heights from $\approx 200~{}\text{pc}$ up to $1~{}\text{kpc}$, continuously increasing with the age of the sub-population. In our simulation there is only a very thin stellar disk component with scale height $z_{\text{s}}=170~{}\text{pc}$, and no gas and thick disk component as compared to the MW. This discrepancy does, however, not affect the objective of this study. The thick disk has a much higher velocity dispersion and is less prone to spiral perturbations. When we will apply RoadMapping to real MW data, the gas disk can be included as an additional component in the mass model, analogous to Bovy & Rix (2013), and if needed also the thick disk. By slicing data according to MAPs, the thick disk will be accounted for implicitly in the tracer selection.

Disk fraction.

The strength of the perturbations in the disk depends on the disk fraction (e.g., D’Onghia 2015). The total disk mass in the simulation might be only 4% of the halo mass, but within $2.2R_{\text{s}}$ it is already 50% of the total mass. It is still under debate if the MW disk is maximal (i.e., rotational support of the disk at $2.2R_{\text{s}}$ is $\sim 55-90\%$; Sackett 1997; Binney & Tremaine 2008, §6.3.3). Bovy & Rix (2013) found, for example, a maximum disk with rotational support $\left(v_{\text{circ,disk}}/v_{\text{circ,total}}\right)^{2}=(69\pm 6)\%$ at $2.2R_{\text{s}}$ with $R_{\text{s}}=2.15\pm 0.14~{}\text{kpc}$. The DEHH-Pot (and therefore our $N$-body model) has $\left(v_{\text{circ,disk}}/v_{\text{circ,total}}\right)^{2}\approx 47\%$ at $2.2R_{\text{s}}$ and is therefore slightly sub-maximal (see Figure 3). The decomposed rotation curve in Figure 3 is qualitatively similar to the MW models by McMillan (2011) (their Figure 5) and by Barros et al. (2016) (their Figure 5, left panels, model MI). Their disks dominate between $R\approx 2.5-8.5~{}\text{kpc}$ and $R\approx 5-9.5~{}\text{kpc}$, respectively, while our simulation’s disk dominates between $R\approx 2-7~{}\text{kpc}$. We account for this slight difference by also drawing mock data sets from regions around $R_{0}=5~{}\text{kpc}$ (see Section II.2).

Number of spiral arms.

The exact number of spiral arms in the MW is still under debate. The distribution of star-forming regions traced, e.g., by HII regions (Georgelin & Georgelin, 1976), maser sources (Reid et al., 2009, 2014), and young massive stars (Urquhart et al., 2014), suggest that the MW has four major spiral arms (see also Vallée (2008, 2014) and references therein). Observations in the infrared, e.g., of old red-clump giant stars in the Spitzer/GLIMPSE survey (Churchwell et al., 2009) or of stellar NIR-emission by the COBE satellite (Drimmel & Spergel, 2001) indicate that the MW is a grand-design two-armed spiral. One hypothesis is that the four-armed spiral observed in young stars and gas is due to the response of the gas to the two-armed spiral in old stars (Drimmel, 2000; Martos et al., 2004). Our simulation is overall a four-armed spiral galaxy, with the mode $m=4$ dominating between $R=4-7~{}\text{kpc}$, $m=2$ at smaller and $m\geq 6$ at larger radii (see Figure 2(b)).

Strength of spiral arms.

As demonstrated in Figure 2(a) the spiral arms introduce strong peak-to-peak differences in the stellar surface density (e.g., $\sim 200\%$ at $R=5~{}\text{kpc}$). The disk dominates inside $R=8~{}\text{kpc}$ (see Figure 3), and in Section IV.1.2, Figure 7, we will see that the spiral arms introduce relative perturbations in the total gravitational forces of up to $30\%$. Figure 5(c), which we will discuss in Section IV.1.1, suggests an excess of stars with radial velocities up to $50~{}\text{km s}^{-1}$ in our simulated galaxy as compared to an axisymmetric model. This can be compared to Reid et al. (2014), who measured, for example, that typical peculiar non-circular motions in the MW spiral arms were around $10-20~{}\text{km s}^{-1}$ for $R\gtrsim 4~{}\text{kpc}$. Siebert et al. (2012) and Bovy et al. (2015) found velocity fluctuations of the same order in the solar neighborhood. We therefore expect the strength of the perturbation to the total potential due to the spiral arms in this simulation—especially inside $R=8~{}\text{kpc}$ where most of our mock data is drawn from—to be similar or even stronger than those in the MW.

III. RoadMapping modeling

In this section we summarize the mathematical ingredients of RoadMapping, and motivate the DF and potential model that we are going to fit to the data. RoadMapping makes extensive use of the galpy python library by Bovy (2015)333The galpy python package by Bovy (2015) can be downloaded from http://github.com/jobovy/galpy.. For full details on the RoadMapping machinery see also Paper I.

III.1. Likelihood

As already laid out in Section II.2, we use as data the 6D $(\boldsymbol{x}_{i},\boldsymbol{v}_{i})$ coordinates of $N_{*}$ stars within a spherical survey volume. The corresponding, purely spatial selection function $\text{SF}(\boldsymbol{x})$ is

[TABLE]

with $\boldsymbol{x}_{0}=(R_{0},\phi_{0},z_{0}=0)$ from Table 1.

Given a parametrized axisymmetric potential model $\Phi(R,z)$ with parameters $p_{\Phi}$, the probability that the $i$-th star is on an orbit with the actions

[TABLE]

is proportional to the given orbit distribution function $\text{DF}(\boldsymbol{J})$ with parameters $p_{\text{DF}}$,

[TABLE]

The joint likelihood of a star being within the survey volume and on a given orbit is therefore

[TABLE]

The details of how we numerically evaluate the likelihood normalization to a sufficiently high precision are discussed in Paper I.444In the terminology of Paper I we use the high numerical accuracy $N_{x}=20$, $N_{v}=28$, $n_{\sigma}=5.5$ to calculate the likelihood normalization, or in other words, to evaluate the spatial and velocity integrals over the qDF within the survey volume.

In the scenario considered in this paper it can happen that there are a few ($\sim 1$ in 20,000) stars entering the catalog that are for some reason on rather extreme orbits, e.g., moving radially directly towards the center. These kinds of orbits do not belong to the set of orbits that we classically expect to make up an overall smooth galactic disk. To avoid that such single stars with very low likelihood have a strong impact on the modeling we employ here a simple strategy to ensure a robust likelihood,

[TABLE]

where $\epsilon=0.001$ for $N_{*}=20,000$ stars and $\text{median}(\mathscr{L})$ is the median of all the $N_{*}$ stellar likelihoods $\mathscr{L}_{i}$ with the given $p_{\Phi}$ and $p_{\text{DF}}$. This robust likelihood method was not used in Paper I. For future applications to real MW data an outlier model similar to the one in Bovy & Rix (2013) could be added: a stellar outlier distribution with constant spatial number density and velocities following a broad Gaussian.

Following Paper I, we assume for now uninformative flat priors on the model parameters $p_{\Phi}$ and $p_{\text{DF}}$ and find the maximum and width of the posterior probability function $pdf(p_{\Phi},p_{\text{DF}}\mid\text{data})\propto\prod_{i=1}^{N_{*}}\mathscr{L}_{i}\times prior(p_{\Phi},p_{\text{DF}})$ using a nested-grid approach and then explore the full shape of the $pdf$ using a Monte Carlo Markov Chain (MCMC)555We use the MCMC software emcee by Foreman-Mackey et al. (2013).. Full details on this procedure are given in Paper I.

III.2. Distribution function model

The most simple action-based orbit DF is the quasi-isothermal DF (qDF) introduced by Binney (2010) and Binney & McMillan (2011), which has been a successful ingredient in Paper I and many disk modeling approaches (Bovy & Rix, 2013; Piffl et al., 2014; Sanders & Binney, 2015). The exact functional form of the qDF$(J_{R},L_{z},J_{z}\mid p_{\text{DF}})$ is given, for example, in Binney & McMillan (2011), or in Equations (2)-(4) of Paper I.

The qDF is expressed in terms of actions, frequencies, and scaling profiles for the radial stellar tracer density $n(R_{g})$, and velocity dispersion profiles $\sigma_{z}(R_{g})$ and $\sigma_{R}(R_{g})$. The latter are functions of the guiding-center radius $R_{g}$, i.e., the radius of a circular orbit with given angular momentum $L_{z}$ in a given potential. We set the scaling profiles to

[TABLE]

The free model parameters of the qDF are

[TABLE]

In an axisymmetric potential superimposed with non-axisymmetric perturbations the actions are not exact integrals of motion. But if one simply considers action-angles as phase-space coordinates, a DF that is a function of the actions only might still be a good model for $\text{DF}(\boldsymbol{x},\boldsymbol{v},t)=\text{DF}(\boldsymbol{J},\boldsymbol{\theta},t)$ at a given point in time and if the system is well-mixed in phase.

We motivate the use of the qDF as specific action-based DF for the simulation snapshot in this work as follows: There is no stellar abundance or age information in the simulation. We therefore cannot define stellar sub-populations, as we normally would for the MW (see Paper I and Bovy & Rix 2013). However, the disk of the galaxy simulation was originally set up as a single axisymmetric flattened particle population whose density decreases exponentially with radius (see Section II.1). This is actually very similar to the stellar distribution generated by a single qDF (see, e.g., Ting et al. 2013). As all particles in the disk have evolved for the same $\sim 250~{}\text{Myr}$ since its axisymmetric set-up, we can consider the disk essentially as a mono-age population. All of this motivates us therefore to use one single qDF to model the whole disk.

Locally, the current particle distribution in the snapshot at hand might be dominated by non-axisymmetries which evolved later in the simulation. We have no indication if for small survey volumes the qDF is still a good model for the data. We will use it anyway—to see how far we can get with the simplest model possible and to test if actions are still informative in this case.

III.3. Potential model

In all RoadMapping analyses in this work we will fit an axisymmetric gravitational potential model consisting of a (fixed and known) Hernquist bulge, a free Hernquist halo and a free Miyamoto-Nagai disk (Miyamoto & Nagai, 1975),

[TABLE]

where $a_{\text{disk}}$ and $b_{\text{disk}}$ are the equivalents of a disk scale length and scale height. Using Hernquist profiles for halo and bulge is motivated by our knowledge of the snapshot galaxy, and we fix the bulge’s total mass and scale length to the true values (see Section II.1 and Table 3). As the bulge contribution to the total radial force at $R_{\text{\sun}}\equiv 8~{}\text{kpc}$ is only $\sim 9-10\%$, this will not give the modeling an unfair advantage. The free model parameters of the halo are the halo scale length $a_{\text{halo}}$ and the halo fraction, i.e., the relative halo-to-disk contribution to the radial force at $R_{\text{\sun}}$, defined as

[TABLE]

As a parameter that scales the total mass of the galaxy model we use the circular velocity at the “solar” radius $R_{\text{\sun}}$,

[TABLE]

The total set of free potential model parameters is therefore

[TABLE]

We will call this potential model the MNHH-Pot (Miyamoto-Nagai disk + Hernquist halo + Hernquist bulge) in the remainder of this work.

To estimate the stellar actions $\boldsymbol{J}=(J_{R},L_{z},J_{z})$ in the axisymmetric MNHH-Pot, we use the Stäckel fudge algorithm by Binney (2012) with fixed focal length $\Delta=0.45$, and interpolate the actions on a grid (Binney, 2012; Bovy, 2015). We made sure that the accuracy of the parameter estimates are not degraded by interpolation errors.666For the action interpolation grid following Bovy (2015), we use $R_{\text{max}}=40~{}\text{kpc}$, $n_{E}=70$, $n_{\psi}=40$, $n_{L_{z}}=50$ in their terminology.

Galaxy disks in general, as well as the simulated disk in this work, have exponential radial density profiles. A single Miyamoto-Nagai disk is more massive at large radii than an exponential disk (see e.g., Smith et al. 2015). By construction, the DEHH-Pot introduced in Section II.3 is therefore better suited to reproduce the overall density distribution in the simulation than the MNHH-Pot (as we will see in Figure 6 in the next section). However, the closed form expression of the Miyamoto-Nagai potential in Equation (17) has the crucial advantage of allowing much faster force and therefore action calculations. In addition, by using a potential model where we already know that it is not be the optimal model for the galaxy’s disk, we challenge RoadMapping even further.

IV. Results

At the core of this work is a suite of 22 data sets consisting of the phase-space coordinates of stellar tracer particles, drawn from the spiral galaxy simulation snapshot introduced in Section II.1. Each data set comes from a different survey volume within the galaxy’s disk (see Section II.2). We modeled all data sets with RoadMapping as described in Section III, by fitting to it a single qDF (see Section III.2) and the potential model MNHH-Pot (introduced in Section III.3). This resulted in 22 independent measurements of the simulated galaxy’s potential and DF.

We present our results in two steps. In Section IV.1 we look at one of these RoadMapping models in detail. In Section IV.2 we then compare all 22 RoadMapping results, and discuss their differences in the context of spiral arms.

IV.1. An axisymmetric galaxy model from RoadMapping

In this section we will discuss all aspects of a RoadMapping model for one single data set. This data set has $N_{*}=20,000$ stars that were drawn from the spherical volume with $r_{\text{max}}=4~{}\text{kpc}$ centered on a spiral arm at the “solar” radius $R_{0}=8~{}\text{kpc}$. This volume is shown in orange in Figure 1 (position S8). We chose this volume because of its position centered on a smaller spiral arm, similar to our Sun being located in the Orion spiral arm. It is a bit larger than our conservative guess for the survey volume size for which we currently expect unbiased potential estimates from the final Gaia data release (DR), $r_{\text{max}}=3~{}\text{kpc}$ from the Sun (see Paper I and Section V.3). But improvements in RoadMapping with respect to the treatment of measurement errors might ultimately allow to also model larger volumes. This example survey volume ranging from $R=4-12~{}\text{kpc}$ also has the advantage of being crossed by several spiral arms of different spiral strength (see Figures 1(a) and 2). In this section our goal is now to investigate the ability of the best fit RoadMapping model to serve as an overall axisymmetric model for the galaxy.

The parameters of the best fit MNHH-Pot and qDF recovered with RoadMapping from this data set are summarized in Table 3. The circular velocity $v_{\text{circ}}(R_{\text{\sun}})$ and halo fraction $f_{\text{halo}}$ are especially well-recovered (compare to Table 2).

IV.1.1 Recovering the stellar distribution

The RoadMapping fit itself takes place in action space. However, an important sanity check to decide if the fit was successful, is to test if the best fit RoadMapping model (i.e., best fit action-based DF in best fit potential and in given selection function) generates a stellar distribution that reproduces the distribution of data points in observable phase-space, $(\boldsymbol{x},\boldsymbol{v})$. This comparison is shown in Figure 5.

We note that the spiral arms introduce very strong non-axisymmetries in the data, both in the spatial and the velocity distribution (especially in $v_{R}$, where a significantly larger number of stars move outward than inward as compared to an axisymmetric model). We therefore compare the data and fit separately for different spatial regions, $R>8~{}\text{kpc}$ and $R<8~{}\text{kpc}$, and $\phi<5^{\circ}$ and $\phi>5^{\circ}$. In the region where the spiral arm dominates (blue in Figure 5) the best fit RoadMapping model is actually a very poor model. However, what the model underestimates in the spiral arm, it slightly overestimates in the other regions and is therefore indeed something like a good average model for the overall distribution. The region at $R>8~{}\text{kpc}$ and $\phi<5^{\circ}$ (yellow in Figure 5), where neither the spiral arm nor the inter-arm regions dominate strongly, is especially well-described by the model.

Overall the qDF appears to be a good model for unperturbed regions of the disk and averages over spiral arms.

IV.1.2 Recovering the gravitational potential

As shown in the previous section the best fit RoadMapping model seems to reproduce the average stellar phase-space distribution quite well. But is the corresponding potential close to the true potential?

Figure 6 compares the true potential from the simulation snapshot (symmetrized by averaging over the whole $\Delta\phi=2\pi$) and the axisymmetric reference DEHH-Pot from Table 2 with the best fit MNHH-Pot from the RoadMapping analysis. In particular, Figure 6 illustrates the overall matter density distribution, the rotation curve and the surface density profile. Figure 7 compares the true and recovered (median) gravitational forces at the position of each star in the data set.

The recovery of density, surface density and circular velocity curve is especially good in the region where most of the stars are located, around $R\sim 6~{}\text{kpc}$ and in the plane of the disk. In large regions inside the survey volume, and even outside, the density is recovered to within 15%. The circular velocity curve is recovered to within 5%, which is also approximately the extent of perturbation that the spiral arms cause with respect to a smooth rotation curve. There is however a very small ($<1.5\%$) underestimation of $v_{\text{circ}}$ at larger radii. We suspect that this bias is introduced by the spiral arms (see discussion in Sections IV.2.6 and IV.1.4). The overall surface density profile is a bit overestimated ($\sim 15\%$) at smaller radii; this is clearly due to the local spiral arm at $R\sim 6~{}\text{kpc}$ with its higher surface density and many stars entering the analysis, which bias the result and which was to be expected. At larger radii ($R\sim 10-12~{}\text{kpc}$) the fit of the local density in the disk and surface density profile starts to flare due to the choice of the Miyamoto-Nagai disk family with its shallow profile (see Section IV.1.4). But again, where most of the stars are located, our RoadMapping model is a very good average model for the true galaxy.

The aspect of the potential to which the stellar orbits are actually sensitive is the gravitational forces. In Figure 7 we therefore compare the true force that each star in the data set feels (i.e., the radial force, $F_{R,T}(\boldsymbol{x}_{*,i})$, and vertical force, $F_{z,T}(\boldsymbol{x}_{*,i})$, calculated as the sum of the individual contributions by each particle in the simulation and the analytic DM halo at the position of each star $\boldsymbol{x}_{*,i}\equiv(x_{i},y_{i},z_{i})$) with the force that the RoadMapping median model predicts for each star ($F_{R,M}(\boldsymbol{x}_{*,i})$ and $F_{z,M}(\boldsymbol{x}_{*,i})$; $M$ for “median model”). We scale the difference between truth and model by a typical radial or vertical force at the given radius for which we use

[TABLE]

where $v_{\text{circ},S}$ and $F_{z,S}$ are the circular velocity and vertical force evaluated in the “true symmetric” reference potential DEHH-Pot in Table 2. As the typical vertical force at a given radius we use $F_{z,S}$ evaluated at the scale height $z_{\text{s}}=0.17~{}\text{kpc}$ of the disk. This is of the same order as the true vertical force averaged over all stars at this radius. Figure 7 shows therefore

[TABLE]

for each star ($R_{i}=\sqrt{x_{i}^{2}+y_{i}^{2}}$). The recovery is as expected: The true vertical force is stronger in the spiral arms (due to the higher surface density) and weaker in the inter-arm regions as compared to the axisymmetric best fit model. The radial force is well recovered where the majority of the stars are located, i.e., in the wide inter-arm regions and in the peaks of the spiral arms. Misjudgments happen in the wings of the spiral arms: The true radial force (i.e., the pull towards the galactic center) is stronger at the outer edge/leading side of the spiral arm because of the additional gravitational pull towards the massive spiral arm, and for the same reason weaker at the inner edge/trailing side. Overall the recovered RoadMapping model appears to be a good mean model, averaging over spiral arms and inter-arm regions.

IV.1.3 Recovering the action distribution

In Sections IV.1.1 and IV.1.2 we have demonstrated the goodness of the fit in the configuration space of the data, and of the recovered gravitational potential. What RoadMapping is actually fitting, however, is the distribution in action space. Figure 8 compares the data and the model action distribution (generated by the best fit qDF) given the best fit median MNHH-Pot in Table 3. (We use this axisymmetric potential to calculate the actions which lead to the best fit model, and do not attempt to estimate the true actions in the true potential.)

We note that the radial and vertical action distribution fits quite well; the axisymmetric model however contains many more stars on close-to-circular orbits $(J_{R}\sim 0,J_{z}\sim 0)$ than the simulation. In the data set there is an excess of stars in the galactic plane $(J_{z}\sim 0)$ that have more eccentric orbits than the axisymmetric model would predict. In Figure 5(a) we have marked the radial extent of the stronger spiral arm with dotted lines ($R_{\text{spiral}}\in[5.6,6.8]~{}\text{kpc}$), and overplotted the corresponding angular momenta $L_{z}=R_{\text{spiral}}\times v_{\text{circ}}(R_{\text{spiral}})$ in Figure 8. This serves as a rough estimate for the region in action space where we expect the stars of this spiral arm to be located. It is again obvious that this spiral arm contains (i) more stars in general and (ii) more stars with eccentric orbits $(J_{R}>0)$ which are (iii) mostly located close to the plane $(J_{z}\sim 0)$, as compared to the axisymmetric model. All of this confirms our expectations for orbits in a spiral arm.

One of the open tasks that the Galactic dynamical modeling community faces is the description of the orbit distribution of spiral arms. The above exercise of comparing the data and the model actions in a best-fit axisymmetric potential should therefore be performed for any future application to data in the Milky Way as well. It could help to learn more about the approximate orbits that stars move on in real spiral arms, and how spiral arms perturb axisymmetric action DFs.

IV.1.4 Calibrating the method by modeling a snapshot without spiral arms

To better understand the effect of spiral arms on the RoadMapping modeling in the previous sections, we also performed a calibration test run. For that we applied RoadMapping to a mock data set drawn from the same spherical survey volume as in Section IV.1 (centered at S8 with maximum radius $r_{\text{max}}=4~{}\text{kpc}$), but from the initial axisymmetric snapshot of the galaxy simulation ($t=0~{}\text{Myr}$), in which no spiral arms had evolved yet.

The spatial tracer distribution was very well recovered by RoadMapping, as well as the distribution of tracers in $v_{R}$ and $v_{z}$. There were some deviations in the $v_{T}$ distribution between data and best fit model, however, with our model predicting more asymmetric drift. We attribute this to the simple initial setup of the disk’s stellar velocities, which does not follow a physical distribution function but only assumes a triaxial Gaussian velocity distribution (Springel et al., 2005). We expect the modeling in this calibration run to be slightly biased by this. Other than that, this confirms our expectation that the qDF is indeed a good model for the unperturbed disk, and that the deviations between data and model in Figure 5 were purely due to the spiral arms.

The surface density profile generated by the best fit potential was a perfect fit inside of $R\sim 7~{}\text{kpc}$, but it started flaring further out. This was expected, because the Miyamoto-Nagai disk is known to be more massive at larger radii than an exponential disk (Smith et al., 2015). Also, if the fit is driven by the majority of stars, the fit is expected to be better at smaller radii with its higher tracer number density. The same flaring in the surface density also showed up in Figure 6. Overall our chosen potential model will therefore systematically bias vertical forces at large radii to be too strong as compared to the truth.

The recovered potential parameters for the initial snapshot, $a_{\text{disk}}=(3.73\pm 0.04)~{}\text{kpc}$, $b_{\text{disk}}=(0.34\pm 0.03)~{}\text{kpc}$, $f_{\text{halo}}=(0.52\pm 0.03)$, and $a_{\text{halo}}=(24\pm 2)~{}\text{kpc}$, are consistent with the spiral arm-affected measurements in Table 3 to within 2, 4, 1, and 1.5 times the statistical error, respectively. While the halo scale length in this calibration run is also consistent with the truth within 3 times the error, we expect it to be underestimated to at least partly account for the Miyamoto-Nagai disk being too massive at large radii. As we will see later (in Section IV.2.3 and Figure 12) we seem to need an even larger survey volume to have enough radial coverage to constrain the halo scale length properly.

Compared to the truth, the circular velocity curve derived from the initial snapshot was underestimated by $\sim 2.5\%$; $v_{\text{circ}}(R_{\text{\sun}})=(219.36\pm 0.08)~{}\text{km s}^{-1}$. For the snapshot with spiral arms (see Figure 6) we observe an underestimation of only $\sim 1.5\%$, so the bias in the initial snapshot might be related to its unphysical setup of the $v_{T}$ velocities. One the one hand, because this underestimation is present for both snapshots, it could also be a systematic bias due to the chosen Miyamoto-Nagai disk model. On the other hand, the bias also shows up in the DEHH-Pot’s circular velocity curve in Figure 6. The DEHH-Pot has a more realistic shape and was found as a direct fit to the spiral arm affected disk particle distribution. The $\sim 1.5\%$ underestimate is therefore most likely introduced by the spiral arms (see Section IV.2.6).

Because the simulation immediately develops spiral arms after the initial snapshot, there was no snapshot that was still axisymmetric, yet already in a dynamical steady state.

IV.2. The influence of spiral arms in RoadMapping modeling

In the previous section we showed that for a large survey volume ($r_{\text{max}}=4~{}\text{kpc}$) RoadMapping can construct a good average axisymmetric potential (and DF) model for a galaxy with spiral arms. In the following we want to investigate how this modeling success depends on the position and the size of the survey volume within the galaxy and with respect to the spiral arms.

IV.2.1 A suite of data sets drawn from spiral arms and inter-arm regions

To investigate a range of data sets affected in different proportions by spiral arms, we center our test survey volumes at the positions marked in Figures 1 and 4 (see also Table 1) and consider volume sizes with $r_{\text{max}}\in[0.5,1,2,3,4,5]~{}\text{kpc}$ for $R_{0}=8~{}\text{kpc}$ and $r_{\text{max}}\in[0.5,1,2,3,4]~{}\text{kpc}$ for $R_{0}=5~{}\text{kpc}$ (to avoid the galactic center). As demonstrated in Figure 4 the spiral arm strength is very different in these test volumes. Each data set that we draw from the simulation contains a random selection of $N_{*}=20,000$ stars inside the given spherical volume and we fit a single qDF and MNHH-Pot to it.

IV.2.2 Recovering the circular velocity curves and surface density profiles

It turns out that RoadMapping is successful in finding reasonable and even very good best fit potential models for each one of the 22 test data sets independent of size and location—given the data and limitations of the model. To illustrate this and to make this encouraging result immediately obvious, we show the circular velocity curves and surface density profiles of all analyses in Figures 9-11.

In contrast to Figure 6, these Figures only show the true profiles for the region within the galaxy where the data comes from: averaged over the angular wedge covering the radial extent of the survey volume, $\Delta\phi=\phi_{0}\pm\arcsin(r_{\text{max}}/R_{0})$, and within $|z|\leq 2z_{\text{s}}=0.34~{}\text{kpc}$, i.e., twice the scale height of the disk, which contains most of the disk mass. This is the matter distribution in which the stars are currently moving, and therefore the potential to which the modeling should be sensitive. In Figures 9-11 we also mark the survey volume and the radial bins of size $\Delta R=200~{}\text{pc}$ with the highest number of stars.

Even though the curves vary extremely between the individual data sets, it becomes very obvious that it is indeed the regions in which the majority of stars is located that drives the RoadMapping fit. Furthermore, whether this region is dominated by a spiral arm or an inter-arm region, and even if this region is only as small as $r_{\text{max}}=500~{}\text{pc}$, RoadMapping indeed constrains the local potential where most of the stars of the data set are located. Also, the constraints are not only most accurate but also most precise in these regions.

Only in the two volumes with $r_{\text{max}}=[0.5,1]~{}\text{kpc}$ at position S8 RoadMapping has some difficulties fitting the circular velocity curve; the model expects a flat or falling rotation curve and is presented with a steeply rising rotation curve due to the spiral arm dominating the region. But given that RoadMapping recovers a good average surface density profile and the circular velocity at least at the center of the small volume, the fit is still quite successful.

In an application to real data in the MW we would also have the possibility to impose some informative prior information on the potential shape (e.g., on the rotation curve), to avoid very unrealistic results (see also discussion in Section V.4).

IV.2.3 Discussion of the model parameter recovery

Figures 9-11 in the previous Section have illustrated how well the potential is recovered by RoadMapping. Figure 12 compares the potential and qDF parameters found with RoadMapping to the parameters of the reference DEHH-Pot from Table 2. And at first glance there appear to be several discrepancies. In the following we will discuss the deviations and explain why each set of parameters still corresponds to a good fit to the data.

Overall the statistical random errors on the parameter recovery are very small for $N*=20,000$ and possible systematic errors dominate. There are only a few exceptions ($r_{\text{max}}=500~{}\text{pc}$, $r_{\text{max}}=1~{}\text{kpc}$ at I5, $r_{\text{max}}=2~{}\text{kpc}$ at I8), which we will discuss later.

We will first consider the parameters of the gravitational potential (left column in Figure 12): All volumes recover $v_{\text{circ}}(R_{\text{\sun}})$ within a few $\text{km s}^{-1}$; in the largest volumes—where the circular velocity curve is probed over several kpc—the estimate is the most accurate. The halo fraction $f_{\text{halo}}$ of the radial force at the “solar” radius $R_{\text{\sun}}$ is very well recovered, especially for $r_{\text{max}}\gtrsim 2~{}\text{kpc}$. The estimate that we get for the best fit Miyamoto-Nagai disk scale height $b_{\text{disk}}$ seems to be also approximately independent of the size of the volume. We can even recover the true halo scale length $a_{\text{halo}}$, however only for a volume as large as $r_{\text{max}}=5~{}\text{kpc}$. The models at $r_{\text{max}}=500~{}\text{pc}$ appear to be too small to constrain the halo at all, and the MCMCs diverged completely for this parameter. Smaller volumes that underestimate $a_{\text{halo}}$ get slightly larger estimates for the disk scale length $a_{\text{disk}}$ and the overall radial density slope is then probably closer to the truth, even if the individual parameters are not. Outliers can often be explained by having a look at the data: The large disk scale length recovered from the $r_{\text{max}}=2~{}\text{kpc}$ volume at S5, for example, mirrors the comparably flat matter distribution caused by two spiral arms close together and dominating the volume (see Figure 9(a) and the large $\sigma_{\Delta\text{Spiral}}$ for this analysis in Figure 4(d)).

The right column of Figure 12 compares the recovered qDF parameters for the different survey volumes with the qDF parameters we got from fixing the potential model to the DEHH-Pot and fitting the qDF only in a $r_{\text{max}}=5~{}\text{kpc}$ volume at S8. Even though the qDF parameters for small volumes are widely different for different positions within the galaxy, they all approach the values recovered with the DEHH-Pot for larger volumes. There seems, therefore, to be an overall best-fit qDF describing the average tracer distribution in the galaxy’s disk. The only difference is in the $h_{\sigma,z}$ parameter, where the models fitting a MNHH-Pot recover a slightly larger value than the models using the known DEHH-Pot. The suspected reason is that the Miyamoto-Nagai disk flares at larger radii as compared to the double exponential-disk (see Figure 6), which leads to a less-steep radial decline in the vertical forces, and therefore mean vertical orbital energies $\langle E_{z}\rangle\sim\nu\times J_{z}$, and therefore to a slightly longer $h_{\sigma,z}$ scale length. In general, volumes centered on spiral arms have larger velocity dispersion parameters $\sigma_{R,0}$ and $\sigma_{z,0}$ as compared to volumes at the same radius $R_{0}$ but centered on an inter-arm region. And the volumes at $R_{0}=5~{}\text{kpc}$ with their stronger spiral arms have larger velocity dispersions than those at $R_{0}=8~{}\text{kpc}$—which is what we expect. Most volumes recover similar tracer scale lengths $h_{R}\sim 2.5\pm 0.5~{}\text{kpc}$ close to the known disk scale length $R_{\text{s}}$. Only the volumes centered on the inter-arm region at $R_{0}=8~{}\text{kpc}$ (position I8) recover much longer $h_{R}$. This might be related to the fact that volumes at I8 are dominated by an especially extended inter-arm region. The volumes at I5 with $r_{\text{max}}=[0.5,1]~{}\text{kpc}$ were not able to constrain the tracer scale length at all because of the unfortunate position between the rising density wings of two strong spiral arms (see Figure 10).

There are a few survey volumes for which the recovered parameters show some peculiarities: The models from volumes with $r_{\text{max}}=1~{}\text{kpc}$ at I5 and $r_{\text{max}}=2~{}\text{kpc}$ at I8 reject the DM halo completely, i.e., $f_{\text{halo}}=0$. The corresponding halo scale lengths $a_{\text{halo}}$ are therefore unconstrained,777For these analyses and the $r_{\text{max}}=500~{}\text{pc}$ analyses the fit could not constrain $a_{\text{halo}}$ and the MCMC was diverging. We had to stop the MCMC after some time, so the $a_{\text{halo}}$ might in truth be even less constrained than shown in Figure 12. while the corresponding disk scale heights $b_{\text{disk}}$ are grossly overestimated to account for the missing contribution of the spherical halo. We have investigated the reason for this fitting result and found that for the way in which the spiral arms affect the circular velocity curve in these volumes, the recovered models with unusual radial profiles are indeed a better description for the data (see Figure 9(a) and 10(b)). Also, while most analyses average the vertical forces radially over the spiral arms (see Figures 7, lower right panel), for these analyses the averaging happens vertically, i.e., at approximately one scale height above the plane where the model’s vertical forces are equally well recovered at all radii (in spiral arms and between), while at small and large $|z|$ the model is bad. RoadMapping therefore also found a good average fit model for the stars in these volumes.

Overall we find that if the volume is large enough to average over several spiral arms and inter-arm regions, an unlucky positioning with respect to the spiral arms does not lead to strong biases in the parameter recovery. We stress again that for particularly large volumes, $r_{\text{max}}=5~{}\text{kpc}$, we were able to recover all model parameters, including the halo scale length $a_{\text{halo}}$.

IV.2.4 Recovering the local gravitational forces

In the previous section we found that the potential and qDF parameters recovered from different survey volumes can be quite different. While the differences can be explained qualitatively, it is not yet clear how good the corresponding potential constraints actually are in a quantitative sense. To test this we calculate again $\Delta F_{R}(\boldsymbol{x}_{*,i})$ and $\Delta F_{z}(\boldsymbol{x}_{*,i})$ from Equations (23)-(24) at the position of each star $\boldsymbol{x}_{*,i}$ in each data set (analogous to Figure 7). From the corresponding histograms of number of stars vs. $\Delta F$ we derive the median and the $16th$ and $84th$ percentiles ($1\sigma$ range) and show them in Figure 13. We chose this diagnostic because the forces at the positions of the stars are the quantities of the potential to which our modeling is sensitive.

The important key result from Figure 13 is that we get very close to recovering the true forces $\Delta F(\boldsymbol{x}_{*,i})\lesssim 10\%$ at the positions of the majority of stars in the survey volume, no matter how large or small the survey volume is. On average, the force recovery is also unbiased888An exception are the radial forces for small volumes strongly dominated by spiral arms (e.g., at S5). The small systematic bias in $F_{R}$ is discussed in detail in Section IV.2.6. for the ensemble of stars.

Figure 13 also contains some subtle clues that suggest that the quality of the force recovery could be correlated with the position of the data set with respect to the spiral arms. We investigate this further by relating in Figure 14 the local force recovery, i.e., the distribution of $\Delta F_{R}(\boldsymbol{x}_{*,i})$ and $\Delta F_{z}(\boldsymbol{x}_{*,i})$ for each data set, to the relative spiral contrast within the respective survey volume, $\sigma_{\Delta\text{Spiral}}$ (Equation (6); see also Figure 4(d)). Figure 14 shows that the average fraction of stars for which the recovery of the radial or vertical force is bad (i.e., larger than 10%) increases with increasing spiral contrast $\sigma_{\Delta\text{Spiral}}$. This is as expected: Volumes in which the steep gradient in surface density around a strong spiral arm ($|\Delta_{\text{Spiral},k}|>0$) is not balanced by larger areas with less perturbations ($\Delta_{\text{Spiral},k}\sim 0$) have (i) a large relative spiral contrast $\sigma_{\Delta\text{Spiral}}$, and (ii) a large relative number of stars affected by the non-axisymmetric kinematics of the spiral arms. And for these individual stars, the axisymmetric RoadMapping model is less successful in recovering the correct forces. (But as we saw in Figure 13, the ensemble average is even in these cases unbiased.)

Interestingly, and even though there is some scatter, the force recovery at a given $\sigma_{\Delta\text{Spiral}}$ is on average very similar for the radial and vertical forces (compare the linear fits in Figure 14). This means that RoadMapping attempts to fit both the radial and vertical forces at the positions of the stars, and is not particularly sensitive to just one of them.

As we saw in Figure 4(d), the spiral contrast $\sigma_{\Delta\text{Spiral}}$ increases for the different test volume positions approximately in this order: I8 $\longrightarrow$ S8 $\longrightarrow$ I5 $\longrightarrow$ S5. From Figure 14 it follows that this is also the order in which the accuracy of the force recovery decreases. (We did not include this piece of additional information in Figure 14, but it can be seen in Figure 13, especially for the smaller volumes.)

IV.2.5 Extrapolating the gravitational potential model

Also, it is interesting to see how well the extrapolation of a recovered potential describes the overall gravitational potential of the galaxy. To investigate the extrapolability, we introduce another diagnostic which uses a cylindrical grid centered on the respective positions in Table 1, always having a radius of $r_{\text{max}}=5~{}\text{kpc}$ and a height of $z=1.5~{}\text{kpc}$ both above and below the plane. In the $(x,y)$ plane the regular grid points have a distance of $0.25~{}\text{kpc}$ and in $z$ they have a distance of $0.125~{}\text{kpc}$ to better sample the thin disk (we throw out grid points close to the galactic center with $R<0.125~{}\text{kpc}$, however). We then evaluate at the position $\boldsymbol{x}_{g,j}\equiv(x_{j},y_{j},z_{j})$ of each regular grid point the force residuals

[TABLE]

analogous to Equations (21)-(24). The two panels in Figure 15 show the [$16th$,$84th$] percentile range and the median of the grid points’ distribution in $\Delta F_{R}(\boldsymbol{x}_{g,j})$ and $\Delta F_{z}(\boldsymbol{x}_{g,j})$.

First, we find that the radial forces are overall very well predicted, especially when derived from large survey volumes. There is however an overestimation of $\sim 5-20\%$ in the vertical forces (depending on volume size and position) which is induced by the spiral arms and is partly also due to a systematic error caused by the choice of potential model (see explanation below in Section IV.2.6).

Second, the constraints we get on the spatially averaged forces inside $r<5~{}\text{kpc}$ are almost as good when derived from a survey volume of $r_{\text{max}}=3~{}\text{kpc}$ as compared to survey volumes of $r_{\text{max}}=4$ or $5~{}\text{kpc}$. If we had to decide between a $r_{\text{max}}=3~{}\text{kpc}$ volume with good data quality and a larger volume with worse data quality, we would lose nothing in terms of extrapolability when using the smaller volume (only the halo scale length might not be as well constrained, see Figure 12).

Third, there are further indications in Figure 15 that the position of the survey volume with respect to the spiral arms matters for the force recovery.

In Figure 16 we stress that even more by relating the extrapolability, i.e., the volume-averaged distribution of force residuals, $\Delta F_{R}(\boldsymbol{x}_{g,j})$ and $\Delta F_{z}(\boldsymbol{x}_{g,j})$, to the dominance of the spiral arm $\langle\Delta_{\text{Spiral}}\rangle$ in the survey volume (Equation (5); see also Figure 4(c)). We derive the fraction of grid points $\boldsymbol{x}_{g,j}$ in the reference cylinder for each data set/RoadMapping model with forces that are misjudged by more than 10%, and plot it against $\langle\Delta_{\text{Spiral}}\rangle$. Positive or negative $\langle\Delta_{\text{Spiral}}\rangle$ quantify how much the spiral arms or the inter-arm regions dominate the corresponding survey volume, respectively. A low fraction of grid points with $|\Delta F|>10\%$ means that the extrapolation of the potential model works well.

Firstly, we note that there is a clear trend that the extrapolability gets worse if a spiral arm strongly dominated the survey volume from which the potential constraint was derived ($\langle\Delta_{\text{Spiral}}\rangle\gg 0$). The same trend can be seen for the dominance of inter-arm regions ($\langle\Delta_{\text{Spiral}}\rangle\ll 0$), but it is weaker and less clear.

Secondly, we note again that $F_{z}(\boldsymbol{x}_{g,j})$ is predicted less good than $F_{R}(\boldsymbol{x}_{g,j})$. The reason for this is laid out in Section IV.2.6.

The main result of Figure 16 is, however, the following: The extrapolability of models derived from data sets drawn from survey volumes centered on inter-arm regions appears to be in general better than that of data sets centered on spiral arms. We suspect that the reason for this is that the stellar distribution between spiral arms is smoother, more extended, and closer to the overall axisymmetric average model, such that the potentials recovered from these volumes have real predictive power for a much larger volume.

IV.2.6 Biases in the potential recovery caused by the spiral arms

What are the reasons for the biases that we observe in Figure 15?

The peak of the distribution in $\Delta F_{R}(\boldsymbol{x}_{g,j})$ (and also in $\Delta F_{R}(\boldsymbol{x}_{*,i})$ in Figure 13) is slightly ($\sim 1.5\%$) biased towards an underestimation of $|F_{R}|$ in our RoadMapping models. This bias already showed up in the circular velocity curve in Figure 6 and also in the reference potential model DEHH-Pot. We therefore suspect that this bias is caused by the spiral arms. The line of argument goes like this: Spiral arms are very thin. If a spiral arm crosses the observation volume, both its leading side (at large radii) and its trailing side (at small radii) are also in the volume. Stars on the trailing side feel a lower gravitational pull towards the galaxy center than they would if there was no spiral arm. Because there are in general more stars at smaller radii, the RoadMapping fit is slightly biased to reproduce in general slightly weaker radial forces.999In the special case that the survey volume coincidentally only contains part of a spiral arm—as was the case with the analyses for $r_{\text{max}}=[0.5,1]~{}\text{kpc}$ at position I5 (see Figure 10, right panels)—the fitting behaves differently anyway, as was already discussed in Section IV.2.3.

The peak in the distribution of $\Delta F_{z}(\boldsymbol{x}_{g,j})$ is strongly biased towards an overestimation of $|F_{z}|$ by the RoadMapping model. We illustrate the reason for this in Figure 17, where we show how $\Delta F_{z}(\boldsymbol{x}_{g,j})$ varies as a function of $R$ for three example data sets with $r_{\text{max}}=2~{}\text{kpc}$. One reason for this bias is a relic of the assumed MNHH-Pot disk model family, which flares outside of $R\sim 8~{}\text{kpc}$ (see Section III.3). The vertical forces in this region are therefore much higher in the model than in the true galaxy with its exponential disk. The main reason for this overestimation of $F_{z}$ comes however directly from the spiral arms: There are much more stars in the spiral arms than in the inter-arm regions, and the stars in the spiral arm feel stronger vertical forces because of the higher surface mass density. The RoadMapping fit is driven by the majority of stars and the best fit model therefore predicts in general higher vertical forces. As expected, the overestimation of $F_{z}$ in Figure 15 is especially strong ($\sim 20\%$) for small survey volumes dominated by spiral arms, while small volumes dominated by an inter-arm region result in much better estimates for the spatially averaged $F_{z}(\boldsymbol{x}_{g,j})$ ($\sim 5\%$ bias). Large volumes lie somewhere in between (bias of $\sim 10\%$).

Why do the biases show up in different strength in Figures 13 and 15?

The stellar number asymmetry in the trailing vs. leading sides of spiral arms is much smaller than the stellar number asymmetry in the spiral arm vs. the inter-arm region. The bias is therefore visible in the distribution of $\Delta F_{R}(\boldsymbol{x}_{*,i})$ (because the $F_{R}$ recovery is biased only by a few stars, which leads to a bias that is visible for the majority of stars) and not in $\Delta F_{z}(\boldsymbol{x}_{*,i})$ (because the majority of stars bias the fit and we therefore also recover $F_{z}$ for the majority of stars). The bias becomes particularly pronounced for $\Delta F_{z}(\boldsymbol{x}_{g,j})$ (because the inter-arm regions dominate when averaging spatially which leads to a large average overestimation of $F_{z}$) and stays small for $\Delta F_{R}(\boldsymbol{x}_{g,j})$ (because trailing and leading sides of spiral arms are similarly important when averaging spatially so it becomes visible that the bias is actually not that big).

V. Discussion and Outlook

V.1. On the informativeness of an orbit distribution function

The qDF appears to be very informative. We did expect it to be at least a reasonable model for the overall symmetrized disk of the galaxy simulation, considering its initial set-up as an axisymmetric, exponentially decreasing particle distribution that subsequently evolved as a mono-age population (see Sections II.1 and III.2). In Sections IV.1.1, IV.1.3, and IV.1.4 we demonstrated that the qDF is indeed a good average model for the tracer distribution in a large survey volume—even though the spiral arms did introduce considerable deviations.

We had, however, no indications beforehand of how well the axisymmetric qDF would perform in a small survey volume completely dominated by non-axisymmetric spiral arms. It would not have been surprising if RoadMapping had failed. But in Section IV.2.2 it turned out that the potential measurements were reliable even in most of the small volumes with $r_{\text{max}}=[0.5,1]~{}\text{kpc}$. And the corresponding qDF parameters were tightly constrained and reasonable as well.

We deduce that the qDF is indeed flexible and robust enough to work with data affected by non-axisymmetries.

That the corresponding potential constraints were reliable as well leads to the following conclusion: A potential model that does not fit the gravitational forces acting on the stars appears to lead to such an unrealistic orbit/action distribution, that a fit with even such a simple orbit DF as the qDF is impossible. This demonstrates once more how powerful the concept of an orbit DF is.

V.2. On the restrictiveness of the parametrized potential model

How much does the choice of potential model matter for the success of the modeling?

We used, on the one hand, a bulge and halo model that reproduces the true bulge and halo better than we can hope to use in reality for the MW. The fact that for all except the largest volumes the true halo scale length is not remotely recovered (and some small volumes even have $f_{\text{halo}}=0$), and that the contribution of the bulge to the overall potential is small (i.e., the bulge contribution to the total radial force at $R=8~{}\text{kpc}$ is only $\sim 9-10\%$) remedies this apparent advantage.

On the other hand, we used a disk model, the Miyamoto-Nagai disk, that we chose purely for its convenient parametric form and of which we know that it is not a good model for the simulation snapshot; especially not for the radial density profile at large radii (Smith et al., 2015). As we saw in most figures in this paper this lead to biases in predicting the potential at radii where we have only a few or no stars. But because the spiral arms are such strong perturbations in the overall potential, a better disk model would not give much better results.

It appears that a potential model with a reasonable shape and flexibility (here: disk+bulge+halo structure with 5 free parameters) can do well enough in finding a good fit, both locally for small volumes and overall for large volumes.

This is in agreement with one of our key results of Paper I, where we managed to successfully fit data from a MW-like (but axisymmetric) galaxy model with a bulgeless potential of a restrictive Stäckel form. This was illustrated in Figure 16 of Paper I. Considering that we used there the same number of stars, $N_{*}=20,000$, the potential uncertainties were much greater than in the analogous figure of this work, Figure 6(a). We believe that this is how RoadMapping accounts for an inconvenient potential parametrization—by increasing the uncertainty of the model estimate—which is exactly as it should be.

V.3. Gaia measurement errors and choosing the survey volume size

Considering measurement uncertainties of distances and proper motions, we found in Paper I that for a survey volume with $r_{\text{max}}=3~{}\text{kpc}$, distance uncertainties of $<10\%$ and proper motion uncertainties of less than $3~{}\text{mas yr}^{-1}$ RoadMapping still gives unbiased parameter results. Even if the proper motion errors are not perfectly known.

The measurement uncertainties of Gaia in proper motions (already in the first data release $\delta\mu\sim 1~{}\text{mas yr}^{-1}$; Lindegren et al. 2016) and in distances (at least for the final data release within $3~{}\text{kpc}$ and for bright stars; de Bruijne et al. 2014) lie below these limits.

In paper I we focused on recovering completely unbiased model parameters and found that RoadMapping is robust against moderate deviations of the model assumptions. In this work we released the condition that the model parameters themselves had to be recovered accurately, but allowed RoadMapping to simply find an overall best fit for the data strongly affected by spiral arms—which was surprisingly successful in recovering the local potential even if the model parameters were not recovered.

We therefore presume that in reality we probably have an even larger margin of error than we found in Paper I, and before the measurement uncertainties noticeably muddle the constraints.

In addition, we found in this work that a volume of $r_{\text{max}}=3~{}\text{kpc}$ should be already big enough to find an overall best fit axisymmetric model for the Galaxy. At larger distances dust starts affecting the measurements. And inside of $R=3-4~{}\text{kpc}$ the stellar motions become increasingly non-axisymmetric, possibly because of the Galactic bar (e.g., Reid et al. 2014; Bovy et al. 2015, and others, see Introduction in Section I).

Overall we should therefore be very well-off by applying RoadMapping to the final Gaia data set within $r_{\text{max}}=3~{}\text{kpc}$ only.

How well we can do with the first few Gaia data releases remains to be seen. The Gaia DR1 from September 2016 has parallax measurements of $\sim 16\%$ for red clump giants at a distance of $\sim 500~{}\text{pc}$ from the Sun (de Bruijne et al., 2014; Michalik et al., 2015), which is not precise enough for RoadMapping. We could, however, use photometric distances, which should be precise enough for red clump stars and even extend over a larger volume. The Gaia DR2 in April 2018 might however already cover $\sim 2~{}\text{kpc}$ with good parallaxes for fainter stars as well.

V.4. Spiral arms in the solar neighbourhood

The Sun is located in one of the smaller spiral arms of the MW, the local Orion spur/arm (Morgan et al., 1953). Two of the MW’s major spiral arms pass by the Sun within a few kpc: The Perseus arm is $\sim 2~{}\text{kpc}$ from the Sun (towards the outer MW) (Xu et al., 2006), and the Sagittarius arm at $\sim 1~{}\text{kpc}$ (Sato et al., 2010) (towards the Galactic center). It is, however, still under dispute which arms are actually major arms of the MW (Blaauw, 1985; Xu et al., 2013; Zhang et al., 2013).

How reliable the RoadMapping results from the Gaia DR1 will be depends on the strength of the local Orion arm, which will dominate the survey volume. RoadMapping had, for example, some difficulties recovering the circular velocity curve for small volumes (see Figure 10).

However, recent measurements of the MW’s rotation curve (Bovy et al., 2012d; Reid et al., 2014) confirm again that it is flat. We could impose this condition as a prior constraint in RoadMapping (or fix the rotation curve slope as Bovy & Rix (2013) did in their RoadMapping analysis). Based on independent information on the location and strength of the MW spiral arms, we could also use the present approach to estimate systematic uncertainties.

For later Gaia data releases, where the tracers extend further into the Galaxy ($\sim 2~{}\text{kpc}$ from the Sun), the Sagittarius and Perseus arms could also play a role. In general, RoadMapping should do much better with larger volumes, and if several spiral arms and inter-arm regions within the survey volume average each other out. This averaging should work especially well if the MW is—like the simulation in this study—a four-armed spiral. The exact number of spiral arms is however still disputed (see Section II.5). If Gaia should show that the MW is a two-armed grand design spiral instead, having a very large survey volume would become even more important to achieve a good axisymmetric average RoadMapping model for the MW. Local measurements of the potential as in Figure 10 should however still be possible.

V.5. Absence of a central bar in the simulation

The simulation we analyzed in this work does not have a prominent bar, and so we have not explicitly explored the impact of such a feature. Bars can play an important role in the dynamics when very small volumes near a resonance are considered (see, e.g., Dehnen 2000, and references in Section II.5). When considering volumes of $\gtrsim 1~{}\text{kpc}$, we have no reason to believe that this should severely affect the robustness of such an analysis.

V.6. Interpreting RoadMapping results

The two main findings of this work give a clear directive regarding how we should deal with any future result about the MW’s potential derived with RoadMapping. (i) If the data spans a large volume, a significant proportion of the disk (at least $R\sim 5-11~{}\text{kpc}$), and averages over several spiral arms and inter-arm regions, we can trust and use the resulting model as an overall axisymmetric potential for the MW. For $R\sim 3-13~{}\text{kpc}$ we should even be able to make definite statements about the DM distribution. (ii) If the data do not span such a large volume, we can still believe the local constraints. In particular, the surface density within 1-2 disk scale heights, the circular velocity within the survey volume, and the average gravitational forces where the majority of stars is located.

Fortunately, this is consistent with the procedure by Bovy & Rix (2013), who used RoadMapping to constrain the vertical force $F_{z}$ for each MAP at only one radius that corresponds to a typical radius. It will be interesting to see whether RoadMapping potential constraints from Gaia data will agree with their findings.

In addition, one should always compare the distribution of the data and the recovered model in configuration space $(\boldsymbol{x},\boldsymbol{v})$ and in action space, as we did in Section IV.1. This is not only a sanity check to confirm the goodness of the fit, but it might also reveal some substructure in the data that only becomes visible when comparing it to an axisymmetric smooth model.

V.7. Recovering non-axisymmetric structures in the potential from modeling small volumes

As we saw in Sections IV.2.2 and IV.2.4, RoadMapping also makes a very good attempt at constraining the local potential for volumes as small as $r_{\text{max}}=500~{}\text{pc}$ or $1~{}\text{kpc}$. It recovers, for example, the higher surface density in volumes dominated by spiral arms, or correctly estimates the circular velocity at the median radial position of the stars (see Figure 10).

At a time after the final Gaia data release, when we have data of high accuracy covering a large proportion of the disk, we could make use of this interesting property of RoadMapping modeling. We could split the data set not only into different MAPs analogous to Bovy & Rix (2013), but also into different spatial bins in the $(x,y)$ plane of the MW and model each of the smaller volumes separately. This approach would only probe the local potential, even when using an axisymmetric potential model, and should be sensitive to the overdensities induced by the spiral arms. In this way it should be possible to build up a non-axisymmetric map of the MW potential—albeit with very large spatial pixels—with constraints from dynamical modeling only.

VI. Conclusion

RoadMapping is a well-tested axisymmetric dynamical modeling machinery that simultaneously fits an action-based DF and gravitational potential to the individual 6D phase-space coordinates of stellar populations in the MW disk. RoadMapping builds on previous work by Binney & McMillan (2011), Binney (2012), Bovy (2015), and was first applied by Bovy & Rix (2013). RoadMapping was improved and tested in detail against the breakdowns of its modeling assumptions by Trick et al. (2016).

In this paper we investigated the robustness of RoadMapping when modeling a non-axisymmetric system. We explore this for the first time explicitly, by modeling a simulated spiral galaxy from D’Onghia et al. (2013), and by comparing the results to the true potential. This simulation has stronger spiral arms than we expect in the MW, and—except of the absence of a bar and thick and gas disk components—it has matter components similar to the MW, as discussed in Section II.5. We find that RoadMapping-like action-based dynamical modeling is very robust against perturbations of spiral arms in this simulation, especially if the survey volume is large enough to encompass both spiral and inter-arm regions. In Section IV.1 we demonstrated this in detail for a single RoadMapping analysis of a data set with a spatial coverage of radius $r_{\text{max}}=4~{}\text{kpc}$ around a position equivalent to that of the Sun.

In Section IV.2 we have investigated the role of survey volumes differing in size and position with respect to the spiral arms and inter-arm regions within the simulated galaxy. We find that the gravitational forces are mostly well-recovered at the locations of the stars that entered the analysis. For survey volume sizes $r_{\text{max}}\geq 3~{}\text{kpc}$ the recovered potential model already becomes a good average potential model for a large portion of the galaxy. For some positions of the survey volume center, e.g., in a smooth and not-too-depleted inter-arm region, smaller volumes can also give good overall constraints. If a small volume is dominated by a very strong spiral arm the constraints become less reliable, as expected. The correct DM halo scale length was, however, only recovered for a survey volume as large as $r_{\text{max}}=5~{}\text{kpc}$.

This overall robustness of RoadMapping is particularly notable, as the breakdown of the assumption of axisymmetry implies a breakdown of several model assumptions simultaneously: (i) orbital actions are not fully conserved anymore, (ii) the true potential is not spanned by the family of model potentials, (iii) the quasi-isothermal DF need not, or will not, describe the orbit distribution within spiral arms. However, the qDF seems to be informative enough to guide the fit to potential shapes that correctly measure the average surface density (within $\sim 2\times$ the disk scale height) and the circular velocity where most of the stars that entered the analysis are located—even for small volumes with $r_{\text{max}}=500~{}\text{pc}$ dominated by spiral arms.

The results of this paper imply that RoadMapping should be well-suited to making new measurements of the MW’s gravitational potential with the upcoming Gaia data releases. It might even potentially work with the Gaia DR1 with its smaller coverage of the disk ($r_{\text{max}}\sim 1~{}\text{kpc}$) because the local Orion arm, in which the Sun is located, is thought to be only a minor spiral arm in the MW and should not significantly disturb the modeling.

VII. Acknowledgments

W.H.T. thanks Stephen Pardy for helpful advice on handling $N$-body simulations. W.H.T. and H.-W.R. acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP 7) ERC Grant Agreement n. ${\rm[321035]}$. J.B. received support from the Natural Sciences and Engineering Research Council of Canada and the Alfred P. Sloan Foundation. E.D. gratefully acknowledges the support of the Alfred P. Sloan Foundation and ATP NASA Grant No NNX144AP53G, and support by Sonderforschungsbereich SFB 881 ”The Milky Way System” (A3) of the German Research Foundation (DFG).