© The Authors 2024

1 Introduction

Accurate distance measurements are essential for many fields of astronomy and astrophysics (e.g. Carroll & Ostlie 2017).

Whereas high-precision astrometric data are readily available within the Milky Way for the stellar component (see, e.g., Gaia data release DR3, Gaia Collaboration 2023), obtaining reliable distance estimates for the gaseous component (i.e. for the various phases of interstellar medium (ISM, see e.g., Tielens 2005; Draine 2011) is much more challenging. Estimating distances to molecular clouds is important for understanding their properties, formation, and evolution (Molinari et al. 2014) and their ability to form stars (Klessen & Glover 2016). The same is true for the atomic and ionised components of the interstellar medium (ISM).

Detailed 3D maps of the gas distribution in the solar neighbourhood have been constructed by combining precise parallax measurements of stars from Gaia and photometric measurements of reddening to the same stars (Lallement et al. 2019; Leike et al. 2020; Zucker et al. 2021). However, this approach is currently feasible only for a limited volume of a few kiloparsecs around the Sun. Reliable distances to clouds further away can be obtained from parallax measurements of molecular maser emission from high mass star-forming regions (Reid et al. 2014, 2019), but this approach is time-consuming and cannot be applied to large surveys containing thousands of clouds.

A widely used method to estimate distances to the ISM out to tens of kiloparsecs from the Sun is the kinematic distance (KD) method. This method allows one to derive the distance to a molecular cloud from its line-of-sight (LOS) velocity. Historically, it was developed by van de Hulst et al. (1954) and Oort et al. (1958), who used it to derive the first face-on maps of atomic hydrogen in the Milky Way from 21-cm spectral line observations. The same approach has since been applied several times to produce face-on maps of the neutral and molecular gas in the Milky Way (e.g. Nakanishi & Sofue 2003, 2006; Levine et al. 2006; Soler et al. 2022) and the associated star-formation rate surface density (Elia et al. 2022).

A key assumption of the KD approach is that the gas is in purely circular motion around the Galactic centre. Significant errors in the KD estimations arise if there are deviations from circular motions. Wenger et al. (2018) recently compared kinematic and parallax distances for a sample of 75 Galactic high mass starforming regions, most of which are at distances d < 10 kpc, and found that for their sample kinematic distances were overestimated by ~20% on average, with a standard deviation of ~50%, when using the KD standard approach with the Brand & Blitz (1993) rotation model for the Milky Way (their Model A). This overestimation drops below 20% when using the updated rotation model of Reid et al. (2014) (see Table 4 of Wenger et al. 2018 for a complete summary).

Errors arising from deviations from circular motions can be broadly divided into two categories: (i) random fluctuations around the average streaming motions that do not change the average velocity (e.g. a turbulent velocity dispersion) and (ii) systematic changes in the streaming velocity due to non- axisymmetric features such as spiral arms and the Galactic bar. Reid (2022) studied the effects of random motions in 3D on the KD estimates and found that a velocity dispersion of ~7 km s⁻¹, which is representative of turbulent motions in giant molecular clouds, can lead to significant (>10%) errors in the KD estimates for true distances d ≲ 5 kpc and can also lead to systematic biases of ~20% despite the random motions having zero mean around the underlying circular motions. A particularly extreme example of deviations from circular velocity leading to errors in the kinematic distance estimate is reported for the Perseus Arm at Galactic longitudes ℓ ~ 130° (see Xu et al. 2006). As pointed out by Sofue (2011) in a detailed study of the physical origins of kinematic distance uncertainties, it is expected that these deviations are more pronounced at small Galactic radii and longitudes (see also the discussion by Sofue 2023 for the very inner regions of the Milky Way).

A number of authors have investigated the effects of streaming motions due to spiral arms using simplified models of the Milky Way. For example, Gómez (2006) used 2D hydrodynamical simulations with a simple externally imposed two-armed spiral pattern to compare KD with true distances and found that errors can be large at the position of the spiral arms. Baba et al. (2009) employed self-consistent N-body plus hydrodynamical simulations with a live stellar potential and found that transient and recurring spiral arms can drive strong non-circular motions, meaning that KD estimates, whilst highly correlated with the true distance, can produce errors of 2–3 kpc. Also using hydrodynamical simulations, Ramón-Fox & Bonnell (2018) found that streaming motions can produce systematic offsets of ~1 kpc and errors of ~2 kpc and that the results are sensitive to the assumed spiral arm perturbations. Some works have tried to correct the KD method to account for the systematic non-circular motions due to spiral arms and bar, but the results are affected by large uncertainties in the gas streaming motions arising from these components (Foster & MacWilliams 2006; Pohl et al. 2008). For a summary of the key geometric and physical concepts behind these various sources of KD uncertainties, we refer to Sofue (2011) and Wenger et al. (2018).

The goal of this paper is to further characterise the KD uncertainties caused by streaming motions due to spiral arms and the bar as a function of position in the Galaxy using a much more accurate Milky Way model than previous work and therefore construct maps of the expected systematic uncertainties that can provide useful guidance as to when the KD method should be considered reliable and when it should be avoided. To do this, we constructed a realistic model of the Galactic gravitational potential that includes state-of-the-art constraints on the Galactic bar, Galactic disk, dark matter halo, and spiral arms, and we ran 2D hydrodynamical simulations using this potential. We then compared actual and kinematic distances in the model, paying particular attention to the inner regions of the Galaxy dominated by the bar and to the regions around the spiral arms.

After a brief discussion of the context of this study, in Sect. 1 we introduce our new analytic description of the Milky Way potential and its various components in Sect. 2. We briefly describe our numerical approach and the implementation of the new potential in the AGAMA frameworks combined with AREPO in Sect. 3. Our main findings are presented in Sect. 4, and their implications and limitations are discussed in Sect. 5. Finally, we summarise and conclude in Sect. 6.

2 Galactic potential

We introduce a new gravitational potential for our Milky Waylike simulation in order to investigate how non-axisymmetric perturbations affect kinematic distance estimates. The potential comprises of many components, each providing structure in different parts of the Galaxy, as detailed in the following sections. The corresponding circular-velocity curves are shown in Fig. 1, and the total midplane density profile of our potential, ρ_gal, is shown in Fig. 2.

2.1 Components of the potential

2.1.1 Supermassive black hole Sgr A^⋆

The potential of the central supermassive black hole, Sgr A^⋆, is represented by a Plummer (1911) model: $${\text{Φ}}_{{\text{SgrA}}^{\star }}=-\frac{G{M}_{{\text{SgrA}}^{\star }}}{\sqrt{r^2+b^2}},$$(1)

where G is the gravitational constant, the mass $M_{{\text{SgrA}}^{*}}=4.154\times {10}^6{\text{M}}_{\odot }$ is taken from the GRAVITY Collaboration (2019), r is the spherical Galactic radius and the scale radius b is set to 0.1 pc to avoid a singularity in the potential.

Fig. 1 Circular-velocity curve produced by our model (black line). The contribution of each component of the potential represented by the coloured lines as detailed in the legend. Observation points from Eilers et al. (2019) and Mróz et al. (2019) are also included (colored markers).

Fig. 2 Underlying density distribution within the midplane (z = 0 kpc) of the gravitational potential of the Galaxy.

2.1.2 Nuclear star cluster

The cluster of stars around Sgr A^⋆ also contributes to the potential within the core of the Galaxy, dominating in the innermost few parsecs. In our model, the nuclear star cluster (NSC) follows a flattened Dehnen (1993) density profile as given in (Chatzopoulos et al. 2015, see their Eq. (17)): $${\rho }_{\text{NSC}}=\frac{(3-\gamma )M_{\text{NSC}}}{4\pi q}\frac{a_0}{a^{\gamma }{\left(a+a_0\right)}^{4-\gamma }},$$(2)

where $$a(R,z)=\sqrt{R^2+\frac{z^2}{q^2}.}$$(3)

The parameters γ = 0.71, q = 0.73, a₀ = 5.9 pc and M_NSC = 6.1 × 10⁷M_⊙ are taken from their best-fitting model. Here R refers to the Galactic radius in cylindrical coordinates. We note from our circular-velocity curve (Fig. 1) we see little contribution to the overall potential from Sgr A^⋆ and the NSC. This is due to these components being most dominant in the inner most 100 pc of the galaxy making it difficult to compare against observational rotation curves and terminal velocities. However, we do include these components for the sake of completeness.

2.1.3 Nuclear stellar disk

For the nuclear stellar disk (NSD) surrounding the nuclear region, we adopt the parameterisation from the Jeans modelling analysis of Sormani et al. (2020) based on data from the APOGEE survey (Majewski et al. 2017; Ahumada et al. 2020) and the 86 GHz SiO maser survey of Messineo et al. (2002, 2004, 2005). The density of this component can be written as $${\rho }_{\text{NSD}}={\rho }_1\exp \left[-{\left(\frac{a}{R_1}\right)}^{n_1}\right]+{\rho }_2\exp \left[-{\left(\frac{a}{R_2}\right)}^{n_2}\right],$$(4)

where a is as defined in Eq. (3) but with q = 0.37, and where n₁ = 0.72, n₂ = 0.79, R₁ = 5.06 pc, R₂ = 24.6 pc, ρ₁/ρ₂ = 1.311 and ρ₂ = 153 × 10¹⁰ M_⊙ kpc⁻³, which follows model 3 of Sormani et al. (2020). As shown in Fig. 1, this component dominates between the inner ~20 pc and ~300 pc of the Galaxy. We opt this model over the more recent model of Sormani et al. (2022b) as the density profile of Sormani et al. (2020) is available in a closed, analytical form. This difference only affects the inner most ~300 pc of the simulation, which is only a minor impact in comparison to the larger scale of the whole Galaxy.

2.1.4 Galactic bar

The Galactic bar dominates much of the potential within the inner ~5 kpc of the Galaxy. The most realistic model for this component is the made-to-measure (m2m) model from Portail et al. (2017). It is constrained using red giant stellar density measurements and kinematics from multiple surveys across the entire bar region. Here we make use of the analytical approximation of this model presented by Sormani et al. (2022a), who provide density functions to describe the X-shaped box peanut bar and the long bar. We reiterate the functions used and their parameters here to have a complete description of our potential in this paper: $${\rho }_{\text{bar}}=\underset{\text{bar}}{\underbrace{{\rho }_{\text{bar},1}+{\rho }_{\text{bar},2}}}+\underset{\text{long bar }}{\underbrace{{\rho }_{\text{bar},3}}}.$$(5)

The first component of Eq. (5) corresponds to the X-shaped component of the observed boxy-peanut shape of the bar (Wegg & Gerhard 2013). To describe this, we used modified versions of Eqs. (9) and (10) of Coleman et al. (2020) and Freudenreich (1998), respectively: $$\begin{array}{ll}{\rho }_{\text{bar},1}(x,y,z)=\hfill & {\rho }_1\mathrm{sech}\left(a^m\right)\hfill \\ \hfill & \times \left[1+\alpha \left(e^{-a_{+}^n}+e^{-a_{-}^n}\right)\right]e^{-{\left(\frac{r}{r_{\text{cut}}}\right)}^2},\hfill \end{array}$$(6)

where $$a={\left\{{\left[{\left(\frac{\left|x\right|}{x_1}\right)}^{c_{\perp }}+{\left(\frac{\left|y\right|}{y_1}\right)}^{c_{\perp }}\right]}^{\frac{c_{\Vert }}{c_{\perp }}}+{\left(\frac{\left|z\right|}{z_1}\right)}^{c_{\Vert }}\right\}}^{\frac{1}{c_{\Vert }}},$$(7) $$a_{\pm }={\left[{\left(\frac{x\pm cz}{x_c}\right)}^2+{\left(\frac{y}{y_c}\right)}^2\right]}^{\frac{1}{2}},$$(8) $$r={\left(x^2+y^2+z^2\right)}^{\frac{1}{2}}.$$(9)

Here, α = 0.626 defines the strength of the X-shape whilst c = 1.342 defines the slope of the X-shape in the x − z plane. The scale lengths x₁ = 0.49 kpc, y₁ = 0.392 kpc, z₁ = 0.229 kpc, x_c = 0.75ɪ kpc and y_c = 0.469 kpc, shape the bar along with the shaping parameters c_⊥ = 2.232 and c_‖ = 1.991. The X- shape of the bar trails off with power-law exponents with powers m = 0.873 and n = ɪ.94 with an additonal cutoff radius at r_cut = 4.37 kpc. The density profile is normalised to ρ₁ = 3.16 × 10⁹ M_⊙ kpc⁻³.

The second and third component of ρ_bar describe the ellipsoid shape of the bar, which we split into short (ρ₂) and long (ρ₃) bar components. Both components follow a modified version of Eq. (9) of Wegg et al. (2015) with fit parameters determined by Sormani et al. (2022a) based on N-body simulations of (Portail et al. 2017): $$\begin{array}{ll}{\rho }_{\text{bar},\text{i}}(x,y,z)=\hfill & {\rho }_i{e}^{-a_i^{n_i}}{\mathrm{sech}}^2\left(\frac{z}{z_i}\right)\hfill \\ \hfill & \times e^{-{\left(\frac{R}{R_{i,\text{ out }}}\right)}^{n_{i,\text{ out }}}}{e}^{-{\left(\frac{R_{\text{i in }}}{R}\right)}^{n_{\text{i,in }}}},\hfill \end{array}$$(10)

where i = {2, 3} and $$a_i={\left[{\left(\frac{\left|x\right|}{x_i}\right)}^{c_{\perp ,i}}+{\left(\frac{\left|y\right|}{y_i}\right)}^{c_{\perp ,i}}\right]}^{\frac{1}{c_{\perp ,i}}},$$(11) $$R={\left(x^2+y^z\right)}^{\frac{1}{2}}.$$(12)

We summarise the parameters used for the components 2 and 3 of the bar in Table 1. The total mass contained within the bar is M_bar = 1.83 × 10¹⁰ M_⊙.

2.1.5 Galactic disk – axisymmetric components

The disk potential of our model takes the form of two exponential disk components with a hole in the centre, introduced to make room for the bar. We adopted a modified version of Eq. (3) of McMillan (2017) using an exponential vertical profile. We obtained $$\begin{array}{ll}{\rho }_{\text{disk }}(R,z)=\hfill & \frac{{\Sigma }_1}{2h_1}\exp \left(-\frac{R}{R_{\text{d},1}}-\frac{R_{\text{cut }}}{R}-\frac{\left|z\right|}{h_1}\right)\hfill \\ \hfill & +\frac{{\Sigma }_2}{2h_2}\exp \left(-\frac{R}{R_{\text{d},2}}-\frac{R_{\text{cut }}}{R}-\frac{\left|z\right|}{h_2}\right),\end{array}$$(13)

where Σ₁ = 1.3719 × 10³ M_⊙ pc⁻², R_d,1 = 2 kpc, z₁ = 300 pc, Σ₂ = 9.2391 × 10² M_⊙ pc⁻², R_d,2 = 2.8 kpc, z₂ = 900 pc, and R_cut = 2.4 kpc. The inner cutoff radius R_cut, scale lengths R_d and surface density normalisations Σ are obtained by fitting our model to the circular-velocity curves of Eilers et al. (2019) and Mróz et al. (2019) as shown in Fig. 1, whereas the scale heights h are fixed to the values from McMillan (2017), which are, in turn, obtained from SDSS star counts by Jurić et al. (2008). The parameterisation of the disk keeps the scale height fixed across all Galactic radii for simplicity, despite observations indicating the scale height decreases towards the Galactic centre. For example, h(R = 4 kpc) ~ 180 pc for the thick stellar disk (Wegg et al. 2015).

In order to better represent the vertical acceleration towards the midplane of the Galaxy (for z < 400 pc), we also included two gas disks in the potential, which we took from McMillan (2017) without any further adjustments: $$\begin{array}{ll}{\rho }_{\text{gas}}(R,z)\hfill & =\frac{{\Sigma }_1}{4z_1}\exp \left(-\frac{R_{m,1}}{R}-\frac{R}{R_{d,1}}\right){\mathrm{sech}}^2\left(z/2z_1\right)\hfill \\ & +\frac{{\Sigma }_2}{4z_2}\exp \left(-\frac{R_{m,2}}{R}-\frac{R}{R_{d,2}}\right){\mathrm{sech}}^2\left(z/2z_2\right),\hfill \end{array}$$(14)

where Σ₁ = 53.1 M_⊙ pc⁻², R_d,1 = 7 kpc, z₁ = 85 pc, and R_m,1 = 4 kpc represents the thick H I disk, whilst Σ₂ = 2.18 × 10³ M_⊙ pc⁻², R_d,2 = 1.5 kpc, z₂ = 45 pc, and R_m,2 = 12 kpc represents the thinner H₂ disk. We note that the gas disk in the hydrodynamical simulations in this paper is not self-gravitating and does not contribute to the potential; instead, these two gas disks are included as static components of the potential. It should also be noted that the gas disk potential does not contain a spiral perturbation, as we are interested in how the stellar potential affects the gas distribution in the simulations.

2.1.6 Galactic disk – spiral arms

To generate the spiral arms of the Galaxy, we introduced a perturbation to the stellar disk in the following manner: $${\rho }_{\text{spiral }}(R,z,\varphi )={\rho }_{\text{disk }}(R,z)\cdot \alpha \frac{R^2}{R_0^2}S(R,\varphi ).$$(15)

Here the perturbation strength increases quadratically with radius in order for the spiral arms to be strong enough in the outer regions of the Galaxy. For our fiducial model, we adopt an amplitude of a = 0.36 leading to a spiral arm strength of ~17% of the disk density at the solar radius of R₀ = 8.179 kpc (GRAVITY Collaboration 2019)¹. This is somewhat larger than the 10% value suggested by Eilers et al. (2020); however, we note that we measure the dynamic response in the gas rather than the field star population, which is expected to be stronger due to the more dissipative nature of the ISM. We also note that we expect the impact of the choice of α to sensitively depend on the actual implementation of the thermal physics and chemistry in the code (see also Appendix A). As outlined in Sect. 3, we here adopt a simple isothermal equation of state for the galactic gas. Consequently, we suggest to explore different values of α when accounting for the true multi-phase nature of the Galactic ISM and its multitude of heating and cooling processes (see e.g. Tielens 2005 or Draine 2011).

For the shaping function, S , we make use of a logarithmic spiral arm potential with the width of the arm having a Gaussian profile. We took a modified form of Eq. (8) of Junqueira et al. (2013): $$\begin{array}{l}S(R,\varphi )={\displaystyle \sum _{k=1}^2\{}\exp \left(-\frac{R^2}{{\sigma }_{\text{sp}}^2}\left[1-f_{m_k,{\gamma }_k}(R,\varphi )\right]\right)\\ \hfill -\exp \left(-\frac{R^2}{{\sigma }_{\text{sp}}^2}\right)I_0\left(-\frac{R^2}{{\sigma }_{\text{sp}}^2}\right)\},\end{array}$$(16)

where $$f_{m,\gamma }(R,\varphi )=\cos \left(m(\varphi +\gamma )-\frac{m}{\tan (i)}\ln \left(\frac{R}{R_a}\right)\right),$$(17)

and i = 12.5° and R_a = 9.64 kpc, m₁ = m₂ = 2, γ₁ = 139.5° and γ₂ = 69.75°. σ_sp = 5 kpc is the width parameter of the spiral arm, which corresponds to a physical width of 1.082 kpc perpendicular to the spiral arm. The second term in Eq. (16) is used to normalise the spiral arm potential such that the monopole component is zero. Here, I₀ is the modified Bessel function of the first kind and of zeroth order. This resulting potential is a superposition of two pairs of m = 2 spiral arms with equal amplitude (Li et al. 2022). We note that the spiral arm potential does not have an m = 4 pattern due to the unequal angular separation between spiral arms. We opted for this kind of pattern, as it allows for the angular separation to be adjusted as needed since the spiral arms of the Milky Way are not fixed to 90° in angular separation (Reid et al. 2019) The shape and intensity of the spiral arms at R = 8.179 kpc is shown in Fig. 3.

2.1.7 Dark matter halo

The dark matter halo component follows a spherical Einasto (1969) profile: $${\rho }_{\text{dm}}(r)={\rho }_0\exp \left[-{\left(\frac{r}{a}\right)}^{1/n}\right].$$(18)

The density normalisation, ρ₀ , is determined by using the total mass of an Einasto potential: $$M=4\pi {\rho }_0{a}^{3}n\Gamma (3n),$$(19)

where the total mass is M = 1.1 × 10¹² M_⊙, the Einasto index is n = 4.5 and Γ is the gamma function. The scale radius a is related to the half mass radius r_s by a ≈ r_s (3n − 1 /3)⁻ⁿ. In this case, the half mass radius is r_s = 96 kpc giving a scale radius of a = 0.88 pc These parameters are optimised to simultaneously fit the circular velocity in the inner region of the Galaxy and its mass distribution at larger distances determined from dynamical modelling of satellite galaxies and stellar streams (see our Fig. 4; Cautun et al. 2020; Correa Magnus & Vasiliev 2022; Vasiliev et al. 2021; Koposov et al. 2023).

Fig. 3 Strength of the spiral arms at solar circle, R = 8.179 kpc, as a function of azimuth. Shown is the strength of the spiral arms in black as well as the $f_{m_1,{\gamma }_1}\text{ and }{f}_{m_2,{\gamma }_2}$ components in blue and red, respectively.

Fig. 4 Enclosed mass profile of our fiducial potential (black) compared to the constraints from dynamical modelling of satellite galaxies (magenta: Correa Magnus & Vasiliev 2022, green: Cautun et al. 2020) and streams (blue: Vasiliev et al. 2021, red: Koposov et al. 2023).

2.2 Comparison with observations

The parameters of the fiducial potential were optimised to satisfy a variety of recent observational constraints, as described below.

2.2.1 Axisymmetric components

We begin with the Galactic circular-velocity curve, as illustrated in Fig. 1 for different radial bins. The black line represents our total circular-velocity curve from the axisymmetrised potential, v_c = (R ∂Φ₀/∂R)^1/2, with the colored lines indicating the contributions from individual components based on the choice of parameters outlined above. Here the axisymmetrised potential is obtained from the monopole, m = 0, component of the potential in which has been approximated by a Fourier or multipole expansion (see Sect. 3.2). The resulting circular-velocity curve does not contain perturbations from the spiral arms nor contains the higher order terms needed to describe the full potential of the bar.

We make use of recent measurements of circular velocity data from Eilers et al. (2019), obtained from red giant star observed with APOGEE, WISE and Gaia, and from Mróz et al. (2019), obtained from Cepheid variable stars with Gaia. Both are in a reasonable agreement with each other and provide a coverage of galactocentric radius of 4 ≲ R ≲ 25 kpc. For coverage within the solar circle, R < R₀, we make of use of the terminal velocities measurements from H I and CO observations (Clemens 1985; Fich et al. 1989; Burton & Liszt 1993; McClure-Griffiths & Dickey 2007). We compare the peaks of the resulting longitude-velocity (lv) diagrams from our hydrodynamical simulations with the corresponding terminal velocity measurements of the Milky Way, as discussed in detail in Sect. 4. Here, we opt to compare terminal velocities in the lv diagram instead of the circular-velocity curves within the inner most R < 4 kpc, as rotation curves obtained from observations in this region will include deviations due to the non-axisymmetric nature of the bar (Chemin et al. 2015). These deviations are not present in our axisymmetrised circular-velocity curve.

We also ensured that our potential satisfies observational constraints on the vertical acceleration at |z| = 1.1 kpc derived in Bovy & Rix (2013) from SEGUE kinematics of G-type dwarf stars, as well as the vertical acceleration at |z| = 0.4 kpc obtained in Widmark et al. (2022a) from modelling the phase spiral in Gaia EDR3 (see Fig. 5).

Fig. 5 Vertical acceleration at |z| = 1.1 kpc (solid red) and 0.4 kpc (dashed blue line). The observational constraints from Bovy & Rix (2013) in the inner Galaxy at |z| = 1.1 kpc and from Widmark et al. (2022a) in the outer Galaxy at |z| = 0.4 kpc are plotted by red and blue points.

2.2.2 Non-axisymmetric components

Our potential has two rotating non-axisymmetric components, the bar and the spiral arms. We chose the pattern speeds of the bar and the spiral arms as Ω_bar = −37.5 km s⁻¹ kpc⁻¹ (e.g. Sormani et al. 2015b; Sanders et al. 2019; Li et al. 2022; Clarke & Gerhard 2022) and Ω_spiгal = −22.5 km s⁻¹ kpc⁻¹ (Li et al. 2022), respectively. We checked the consistency of these values by running a small parameter study with our potential to generate longitude-velocity lv diagrams and comparing them to the spiral arm tracks presented in McClure-Griffiths et al. (2004), Reid et al. (2016) and Reid et al. (2019). For the sake of simplicity, we consider both non-axisymmetric components to experience solid body rotation. See Appendix B for the full details of the parameter study. The resonances for the potential can be found in Table 2 and are illustrated in the frequency curves of Fig. 6.

We find that for the pattern speeds we use, the outer Lindblad resonance of the bar coincides with corotation of the spiral arms at ~10.1 kpc. The outer 4:1 resonance of the bar lies close to solar circle at ~8.2 kpc.

Fig. 6 Frequency curve as a function of Galactic radius. The solid black line is the rotational frequency curve of the potential, whereas the dashed and dot-dashed lines are rotional frequency plus or minus 0.25 and 0.5 times the epicyclic frequency, ${\kappa }_0(R)=\sqrt{(2\text{Ω}/R)\text{d}\left(R^2\text{Ω}\right)/\text{d}R}$. The horizontal blue and red lines are the pattern speed of the bar and the spiral arms, respectively. The vertical dashed and dot-dashed lines are the corresponding resonances for bar and spiral arms.

3 Numerical simulations

Here we briefly describe the numerical methods used to simulate the dynamical evolution of the ISM in our Milky Way analog.

3.1 Numerical hydrodynamics

We solved the equations of hydrodynamics with the movingmesh code AREPO (Springel 2010). For isothermal gas in two dimensions, these are $$\frac{\partial \Sigma }{\partial t}+\nabla \cdot (\text{Σ}v)=0,$$(20) $$\frac{\partial \Sigma v}{\partial t}+\nabla \cdot (\text{Σ}v\otimes v)=-\nabla P-\Sigma \nabla \text{Φ},$$(21)

where Σ, v, and P are gas surface density, velocity and pressure, respectively. The simulations are two-dimensional. The pressure is related to the density via the equation of state, $P=c_{\text{s}}^2\text{Σ}$, with the sound speed adopted as c_s = 10 km s⁻¹. The external potential Φ is given by our model for the Galactic potential, as explained in the next section; for ease of interpretation, we do not include the gas self-gravity, star formation or stellar feedback in our models. We note that by choosing a relatively large value for c_s , we are implicitly accounting for some of the turbulent support of the gas disk, something that in reality would be provided by stellar feedback (Mac Low & Klessen 2004; Krumholz & McKee 2005; McKee & Ostriker 2007; Klessen & Glover 2016). This assumption is justified in this case as our focus is on how the large scale dynamics impact kinematic distance estimates and not on the impact by turbulent motions.

AREPO constructs a Voronoi tesselation, in which the mesh generating points are able to flow with the gas in the simulation, resulting in a quasi-Lagrangian approach to modelling the flow properties. We make use of an exact Riemann solver for isothermal flows, and the mesh can refine and derefine with the addition and removal of mesh generating points. This occurs when the mass of a given cell is a factor ∼2 larger or smaller than the cell target mass for the simulation (M_target = 2500 M_⊙). The cell will either split or merge with another with addition or removal of a mesh generating point. For further details of the code base we use, see e.g. Tress et al. (2020a,b). We achieve a minimum value of ∼40 pc. Figure 7 illustrates how our cell size varies with density. The vertical spread seen at small cell sizes is a result of the minimum surface area of the cell being reached, which is set to 6 pc² . This creates a limit on how small our cells can become by not allowing cells smaller than two times the minimum surface area to refine further.

Fig. 7 Relation between cell size and density within our simulation. The cell size in this case is the size of a square with the same area as the cell.

3.2 External Galactic potential

AREPO allows one to include an external gravitational potential in the simulation, but the Galactic potential described in Sect. 2 is significantly more complex than the few built-in analytic models. We calculated the potential using the AGAMA library for stellar dynamics (Vasiliev 2019), which, among other features, contains a powerful framework for constructing and evaluating arbitrarily complex potentials (including time-dependent features such as a rotating bar with a varying amplitude or pattern speed). We created an interface between the two codes that make it possible to use any potential implemented in AGAMA as an external potential in AREPO (in addition to self-gravity of the simulated system, if the latter is turned on). Moreover, a very similar interface is provided for the GADGET-4 code (Springel et al. 2021), which shares a common ancestry with AREPO; both interfaces, as well as the script for generating the Galactic potential from this study, are available in the latest version of AGAMA.

The Galactic potential consists of two general-purpose expansions: Multipole for spheroidal density components, and cylindrical Fourier series (CylSpline) for disk-like components. Each of the two expansions is constructed from the sum of several density components, as detailed in Table 3; the mathematical details of these potential expansions can be found in the appendix of the AGAMA documentation (Vasiliev 2018).

3.3 Initial conditions

The initial conditions of our simulations are simple. We set up the surface density of the gas following an exponential profile similar to Eq. (14): $$\begin{array}{ll}\Sigma (R)=\hfill & {\Sigma }_1\exp \left(-\frac{R_{m,1}}{R}-\frac{R}{R_{d,1}}\right)\hfill \\ & +{\Sigma }_2\exp \left(-\frac{R_{m,2}}{R}-\frac{R}{R_{d,2}}\right),\hfill \end{array}$$(22)

where the parameters for this gas disk are the same as those used for the gas disk potential (see Sect. 2.1.5). We extended this gas disk to R ∼ 30 kpc, at which point we reduced the density significantly to prevent artefacts caused by periodic boundary conditions compromising the Galaxy itself. For simplicity, we initialise the simulation with 250 000 mesh generating points distributed uniformly across a (75 kpc)² box. The mesh is then relaxed with the meshrelax method within Arepo to reach our target mass of 2500 M_⊙ in which cells are refined or derefined according to the criterion mentioned previously. There is no hydrodynamics present in the meshrelax process and as such the gas is fixed until the process is complete.

The velocity of the gas is initialised to be the circular velocity of the axisymmetric terms of the potential, v_c = (R ∂Φ₀/∂R)^1/2, which follows the same circular-velocity curve as Fig. 1. The non-axisymmetric components of the bar and spiral arms are introduced linearly and gradually over the course of 150 Myr to avoid transients, as is customary in this type of simulations (e.g. Li et al. 2022), making use of the time-dependent Evolvingpotential in AGAMA.

4 Results

4.1 The gas response

In order to test kinematic estimates properly, the simulation box needs to be rotated such that the bar is in a similar position with respect to the Sun’s position as it is for the Milky Way. For each simulation output, we rotate the system so that the angle between the bar major axis and the Sun-Galactic centre line is 28° (Bland-Hawthorn & Gerhard 2016). The Sun-Galactic centre distance is assumed to be R₀ = 8.179 kpc (GRAVITY Collaboration 2019). For each simulation we generate an lv diagram assuming the Sun moves with a velocity equal to its circular velocity in the x direction, v_x = v_c(R₀) = 229 km s⁻¹, and has no other velocity components. This value is obtained from the circular-velocity curve produced from the axisymmetric potential as detailed in Sect. 2. It should be noted that the GRAVITY Collaboration (2022) has given a more recent value of R₀ = 8.277 kpc; however we adopted the 2019 value for the calculations in this paper. Additionally, proper motion measurements of Sgr A^⋆ give v_c(R₀) ~ 250 km s⁻¹ (e.g. Reid & Brunthaler 2020), ~20 km s⁻¹ greater than the value we have derived from our circular-velocity curve. This difference can be attributed to both the peculiar velocity of the Sun relative to the local standard of rest (LSR), and the spiral perturbations, which make the LSR velocity deviate from the axisymmetric circular speed.

We compare the structures in the resulting lv diagrams with the spiral arm tracks of Reid et al. (2016, 2019) and McClure-Griffiths et al. (2004), as illustrated in Fig. 8. For the analysis presented in this section, we select the system at 441 Myr. It is very similar to the observations in the lv diagram, and the simulation at this point has been advanced for long enough that the non-axisymmetric components of the potential have had enough time to interact with the gas. The resulting density maps can be found in Fig. 9.

From Fig. 8 we find clear peaks in lv diagram that are associated with the spiral arms generated from the underlying potential. The spiral arms generally trace the spiral arm tracks of McClure-Griffiths et al. (2004), Reid et al. (2016) and Reid et al. (2019) in the regions outwith the Galactic centre (|l| > 50°). Towards the Galactic centre, comparing the spiral arms becomes difficult due to the perturbations generated by the bar. Some features match; however, there are features present in the lv diagram that do not match with any track and vice versa.

Comparing the terminal velocities we find our simulation mostly fall within what is expected for the Milky Way, with the exception of two zones; one at 10° < l < 30° and the other at −20° < l < −10°. The LOS velocities are higher than that of the terminal velocities at these angles. This is a result of the circularvelocity curve at sim3 kpc being relatively steep in comparison to similar potentials (See Li et al. 2022). This is a result of an overlap between the bar and stellar disk potentials, which causes a slight overdensity at this location.

As expected, the potential of the bar strongly influences gas dynamics in the central region of the Galaxy. The gas here follows the typical x1 orbits, a family of orbits elongated parallel to the major axis of the bar (Contopoulos & Grosbol 1989), until it is shocked at the end of the bar, after which it flows inwards on nearly radial orbits. Eventually, the gas stabilises onto x₂ orbits, forming a ring of material at 220 pc from the centre of the Galaxy. This ring is the equivalent of the Central Molecular Zone (CMZ) and is consistent with the larger end of estimates for the Milky Way’s CMZ (e.g. Henshaw et al. 2023).

Outside of the bar region (R > 5 kpc), the gas forms a clear spiral pattern. It is rather complex and has two main components: a two-arm spiral caused by the rotation of the bar, and the four-arm structure created by the spiral component of the potential described in Sect. 2. These two pattern rotate at different angular speeds, Ω_bar = −37.5 km s⁻¹ kpc⁻¹ and Ω_spiгal = −22.5 km s⁻¹ kpc⁻¹ respectively, so they periodically interfere with each other. We plot the polar decomposition of the density map in the bottom plot of Fig. 9 to better illustrate the spiral patterns. Here, a straight line would be consistent with a logarithmic spiral. We observe two gradients of spiral structure: the underlying spiral arm structure from the potential (blue dotted), and an m = 2 spiral being generated by the bar (green dashed) with an estimated pitch angle of 6.5° near the outer Lindblad resonance of the bar. We note that the pitch angle value of the bar-generated pattern depends on the sound speed of the gas, as can be understood from the dispersion relation of spiral density waves in the tight-winding limit (e.g. Binney & Tremaine 2008). Whilst there are linear trends in the plot, there are deviations from the underlying spiral arm structure. These occur at the point where the two components intercept at R ~ 11 kpc and at around the spiral arm crossing point, R_cгoss = 9.64 kpc, where the deviation is a bridging feature between the spiral arm and the bar induced spiral arm.

We extract the exact structure with the filament finding package FILFINDER (Koch & Rosolowsky 2015). This package identifies structures from a 2D image using morphological techniques. Not only does the package provide the spines of the extracted structure, it also provides the masks of the extracted regions. We overlay the extracted spines in light blue in the bottom plot of Fig. 9 and we use the masks to contrast the density of the simulation in the x − y projection (Fig. 9, top).

4.2 Kinematic distance estimates

The kinematic distance is calculated based on the assumption of purely circular orbits within the Milky Way. First, a rotation curve v_c (R) is assumed. Then the galactocentric radius of a given object is calculated by $$R=R_0\sin (l)\frac{v_c(R)}{v_0\sin (l)+v_{\text{los }}},$$(23)

where l is the Galactic longitude of the object, v_los is its LOS velocity, and R₀ and v₀ are the Galactic radius and circular velocity of the Sun, respectively. For consistency, we used the circular-velocity curve generated by our potential as our rotation curve, as shown in Fig. 1. Because R appears on both sides of Eq. (23), it is estimated through an iterative process.

The kinematic distance to the object is then obtained by $$d_k=R_0\cos l\pm \sqrt{R^2-{\left(R_0\sin l\right)}^2},$$(24)

It is possible for the kinematic distance estimate to return an undefined answer as a result of the argument inside the square root being less than zero. This occurs when $$v_{\text{los}}^2>v_{\text{term }}^2-2v_0\left(v_{\text{los}}-v_{\text{term }}\right)\sin l,$$(25)

where v_teгm is the terminal velocity along a given Galactic longitude and is given by (Burton & Gordon 1978) $$\left|v_{\text{term }}\right|=v_c(R)-v_0|\sin (l)|.$$(26)

In other words, the kinematic distance is indefinite when the observed v_los is not possible (e.g. too high) under the assumed rotation curve v_c (R). In this case, the argument of the square root of Eq. (24) is set to zero, which is equivalent to placing the object at the tangent point for a given l, where R = R₀ sin(l). This is also equivalent to setting the velocity to the terminal velocity for a given l.

Equation (24) can give two answers when observing inside the solar circle, resulting in the well-known kinematic distance ambiguity. For the sake of simplicity, we resolve the ambiguity by selecting the kinematic distance closest to the true distance value, but we note that this may be difficult to decide in real observations. In Fig. 10 we present a face-on view of the disk and compare the data in our fiducial model with the ideal case of purely rotational motions. At the left we focus on the velocities and at the right we depict the corresponding distance estimates. In the top row of Fig. 10, we plot the LOS velocity measured by an observer at the solar radius (left) and the KD estimates based on Eq. (24). The kinematic distance map is not a smooth distribution with increasing radius from the Sun. We find deviations from the true values (as illustrated in the middle row of Fig. 10) close to the perturbations caused by the spiral arms and the bar. We also observe quite large deviations close to l = 0° and l = 180°. This arises due to trigonometric effects: as l tends towards 0° or 180°, the value for the Galactic radius, R, obtained from Eq. (23) becomes ill-defined as both sin(l) and v_los tend to 0. These deviations are highlighted in the bottom row of Fig. 10, where we map the difference in LOS velicity and the relative error between the kinematic distance and the true distance. What stands out is the large relative error in the solar neighbourhood. For objects within 500 pc of the Sun, the kinematic distance estimate is highly unreliable and exhibits a broad error distribution. Large systematic errors also occur close to the spiral arm perturbations and at the end of the bar, where the gas flows on x₁ orbits.

When computing kinematic distances, observers typically avoid LOSs within ±15° of the direction towards the Galactic centre and ±20° of the anti-centre due to high errors in these LOSs. (e.g. Anderson et al. 2012; Balser et al. 2015; Wenger et al. 2018). As l tends towards 0° or 180°, the value for the Galactic radius, R, from Eq. (23) becomes more difficult to estimate as both sin(l) and v_los tend to 0. On top of this, the Galactic bar also impacts estimates towards the centre due to the high level of asymmetry in the potential.

To give a more conclusive idea about which LOSs to avoid, we computed the absolute kinematic distance error and plot the median and absolute median deviation (MAD) as a function of Galactic longitude in Fig. 11. The results are split into separate annuli around the Sun’s position, showing that the error in the kinematic distance remains high at distances closer to the Sun with some variation with Galactic longitude. We define an LOS of avoidance for an annulus as the LOS where more than 20% of cells have an absolute kinematic distance error greater than 27%. Our choice of this value is motivated by the study of Wenger et al. (2018), who quantify the uncertainty in the kinematic distance inferred using the Brand & Blitz (1993) rotation curve (their method A) or the Reid et al. (2014) rotation curve (their method B) due to uncertainties in e.g. the solar galactocentric radius and orbital velocity, the measured rotation curve, etc. The average absolute kinematic distance error they find when considering both models is 27%. For LOSs and locations where the systematic error for most points is less than this value, the kinematic distance method should be reliable. On the other hand, if a large fraction of points have systematic errors that exceed this value, this is a good indication that the kinematic distance method will not provide reliable results. The LOSs that should be avoided are shown in red in Fig. 11. For objects within 500 pc of the Sun, 99.5% of the full range of Galactic longitude should be avoided. This fraction remains above 56% out to 5 kpc, but drops to 14% and 16% for the 5–10 kpc annulus and the 10– 20 kpc annulus, respectively. This suggests that the kinematic distance estimate is accurate for distances beyond 5 kpc from the Sun under our criterion for LOS avoidance, but that for closer distances it should be used with great care.

Fig. 8 Longitude–velocity maps of the simulation observed from the Sun’s position (placed at the origin of the overlayed coordinate system in Fig. 9). The bar major axis is rotated by 28 degrees from the LOS passing through the Galactic centre. Overlayed blue dashed lines are the spiral arm tracks of Reid et al. (2016, 2019) and McClure-Griffiths et al. (2004). The green points are the terminal velocities from H I and CO observations. Top: full diagram. Bottom: zoom-in between l = −60° & +60°. The left plot only takes into account gas that lies at R < 6 kpc from the centre, and the right one only considers gas at R > 6 kpc.

Fig. 9 Spiral arm structure. The top plot is the x − y density projection of our simulation at t = 441 Myr. Spiral arms have been enhanced in this image using the mask generated by FILFINDER to brighten overdensities and darken other regions. The bottom panel is the polar decomposition of the density map of the top panel, focusing on the region with galactocentric radius 6 < R < 20 kpc. The blue dotted lines are the positions of the spiral arms according to the second half of Eq. (17). The green dashed lines are the spiral arms generated by the rotation of the bar, with a pitch angle of 6.5° around the outer Lindblad resonance, R = 10.14 kpc. The light blue lines indicate the spiral arm pattern extracted with FILFINDER.

4.3 Location of kinematic distance errors

So far we can see that velocity perturbations generated by the non-axisymmetric components of the potential can produce highly inaccurate kinematic distance estimates along most LOSs for objects close to the Sun. This now poses the question: where can one reliably use kinematic distances?

In Fig. 10, we show the map of kinematic distance errors of our simulation with the density peaks extracted with FILFINDER overlayed onto the maps. The peaks lie close to the regions of low value for the kinematic distance error. However, based on inspection, this only applies to the spiral arm features (i.e. peaks outside of the bar region). To further analyse this, we split the Galaxy into two regions: the bar region (R < 6 kpc) and the disk region (R ≥ 6 kpc).

4.3.1 Bar region

For the bar region, we employ the mask generated by FILFINDER to identify the overdensities from our simulations and applying them to the kinematic distance error map (Fig. 10), splitting the data into overdense regions and underdense regions. We plot the PDF of each region, respectively, and compare their distributions (Fig. 12).

From visual inspection, we see that both distributions peak at around 0, with the underdense regions’ distribution slightly wider than the overdense regions. Since the distribution is nonGaussian, we look at the difference between quantiles to understand the width of the distributions. The difference between the upper and lower 20% quantiles is 0.14 for overdense region, whilst it is wider with a value of 0.22 for the underdense region. Similarly, the values are 0.29 and 0.39 for the upper and lower 10% quantiles, and 0.50 and 0.61 for the 5% quantiles respectively. This suggests that outside of the overdense regions there is a higher probability of a large error and, by consequence, an increased probability of obtaining an incorrect distance via the kinematic distance method. Going from our analysis of the error as function Galactic longitude (Fig. 11), this result is not too much of a surprise given that much of the bar’s influence is lies within l = ±30°, which is typically a LOS of avoidance for all distances away from the Sun.

Fig. 10 Kinematic distance uncertainties at different positions in the disk. The left column presents the LOS velocity map of the simulation and the comparison to the equivalent if the gas follows circular orbits. The top plot illustrates LOS velocity map of the simulation. Similarly, the middle plot shows the map of the LOS velocity derived from the base axisymmetric potential. The last plot shows the difference between simulation and axisymmetric LOS velocities; top plot minus the middle plot. The right column shows the estimated and real distance maps in the simulation. The black lines indicate the spiral arm pattern extracted with FILFINDER. The top plot shows the kinematic distance maps as estimated with Eq. (24). The middle plot is the map of the true distances to the gas cells. The relative error between the kinematic and true distance is shown in the bottom plot.

Fig. 11 Kinematic distance error at different distances from the Sun. The black line represents the median absolute relative kinematic distance error along the LOS of a given Galactic longitude. Each plot represents a sampling annulus centred on the Sun. The blue shaded region is the median absolute deviation (MAD) of the error. The red shaded regions represent the Galactic longitudes where more than 20% of the cells along the LOS have a relative kinematic distance error of 27% or more.

4.3.2 Disk region

As we have done previously, we split the kinematic distance maps of the disk region into spiral arm regions and interarm regions and plot the corresponding distributions of each. In this case we apply this to everything outside of R = 6 kpc.

From Fig. 13, we find similarly shaped distributions in the errors for the Galactic disk as we did for the bar region. Once again we have the distribution peaking at 0 for both spiral and interarm regions, with a a wider distribution in the interarm regions. Looking at the ranges between quantiles, the ranges are more similar between the spiral arm and interarm region. The difference between the upper and lower 20% quantiles is 0.1 for both regions. For the difference in the 10% quantiles the ranges are 0.17 and 0.2 for the spiral arm and interarm regions respectively. Similarly it is 0.26 and 0.39 for the 5% qunatile difference. Although the difference in the quantiles are similar between the regions, the interarm regions has a wider distribution. This indicates that the interarm regions of the Galaxy have a higher probability of an incorrect kinematic distance, similar to that for the underdense regions of the bar though with a lower probability of a larger associated error. The narrower distribution of errors of the spiral arm region is a rather reassuring fact, because it implies that the kinematic distance method is more reliable in the regions of the Galaxy where most of the dense gas and subsequent star-forming regions are found. We look into the dynamics causing this result in the subsequent subsection.

It should be noted that in both overdense regions there is not much skewness in the distribution and, as such, there is no clear way to indicate whether if the kinematic distance would be under- or overestimated. We summarise the statistics generated from the systematic kinematic distance error distributions in Table 4.

4.4 Relation to the velocity deviation

The key parameters that are needed for kinematic distance methods are the LOS velocity of the object and a rotation curve for the Milky Way. As mentioned previously, the Milky Way is not axisymmetric and as such there are deviations away from the rotation curve velocities. Quantifying the correlation between these deviations and the systematic kinematic distance errors can give an insight into how the velocity impacts the kinematic distance estimates.

The deviations from rotation curve can be seen in the radial profile of the azimuthal velocity of the gas. We illustrate this in Fig. 14 where in the inner most 3 kpc we observed deviations up to 70 km s⁻¹ . Between 3 kpc and 10 kpc there are small deviations from the rotation curve, on the order of a few kilometers per second, due to perturbations of the spiral arms in our system. We include the rotation curve from Brand & Blitz (1993) and the universal rotation curve in physical units of Persic et al. (1996) with the updated parameters from Reid et al. (2019) as comparison to other rotation curves used in kinematic distance estimates.

We compute the LOS velocity of our simulations and map it to the face down Milky Way view at the top of the left column of Fig. 10. The non-axisymmetric perturbations of the potential are apparent here. Towards the bar region, we observe a sharp transition across the Galactic centre (as expected from the lv diagram of Fig. 8), with perturbations of the spiral arms appearing as displacements in the contours in comparison to the middleleft plot of Fig. 10; the LOS velocity maps for an axisymmetric potential. When we look at what we expect from an axisymmetric potential, the transition towards the Galactic centre is not as strong and the contours of the LOS velocity towards the outer Galaxy are smoother. When we subtract the two maps, we find that the largest deviations occur at the bar region with differences on the order of 100 km s⁻¹. Outside the bar region we find the LOS velocity difference to be close to zero along the spiral arms of the Galaxy. Gas experiences an acceleration or deceleration as it flows into or out of the spiral arm, respectively as shown in Fig. 15. This can cause shocks causing the gas to get denser as it leaves the potential minima, as in the case for two of our spiral arms. Additionally, the perturbations caused by the bar can also shock the gas and causes further perturbations. In our case, at around R₀ , both the bar generated spiral and two of the spiral arms overlap creating large peaks in density. However, unlike the other two spiral arms, the peak in density for these occur just before passing the potential minima of the spiral arm potential. Other more local sources of velocity perturbations come from expanding HII regions and supernova driven bubbles (see e.g. Barnes et al. 2023 or Watkins et al. 2023 for nearby galaxies, or Zucker et al. 2022 for the solar neighbourhood in the Milky Way).

All of this results in the largest deviations from the rotation curve to occur in the interarm regions, which in turn causes shifts in the estimated kinematic distances, increasing their systematic error within these regions (see the right-hand side of Fig. 10). To be more quantitative, peculiar motions associated with the Galactic spiral arm structure of ∼6 km s⁻¹ in radial and ∼4 km s⁻¹ in tangential direction are reported by Reid et al. (2019). We note that this is very much in line with the streaming motions of up to 6–8 km s⁻¹ illustrated in Fig. 15.

Given how large deviations of the LOS velocity occur in the same regions where the systematic kinematic distance errors are highest, we look at the correlation between the galactocen- tric radial velocity v_r and the distance error. Any deviation from circular motion is likely to occur due to motions in the radial direction. In Fig. 16 we provide the distribution and correlations between absolute values of the peculiar motion |v_p|, the deviations from circular velocity |v_los − v_c,los|, and the difference between the true and kinematic distance |d_k − d_true|, split between arm and interarm regions. There is a slight offset between the distributions for each region, with spiral arms having lower deviations and peculiar motions. As expected, we find a clear positive correlation between the local deviations from circular velocity and the resulting distance error. It can be described by a power law with slope α = 0.94. However, we note that there are several orthogonal spurs in the 2D distribution that do not follow this trend. Additionally, we observe a power-law correlation between the peculiar velocities and the difference in LOS velocity. Here, the slope is α = 1.1. Despite the positive correlations with the LOS velocity deviations, there is no clear correlation between the distance deviation and the peculiar motions.

Fig. 12 Probability distributions of the systematic kinematic distance errors within the overdensities in the bar region (blue) and the underdensities in the bar region (red).

Fig. 13 Probability distributions of the systematic kinematic distance errors within the spiral arm region (blue) and the interarm regions (red).

Fig. 14 Comparison of rotation curve against the average azimuthal velocity (black solid and dashed line respectively). The red and blue lines are the rotation curve from Brand & Blitz (1993) and Reid et al. (2019).

Fig. 15 Deviations from circular velocity and relation to column density at the solar circle. The top panel illustrates the Galactic azimuthal profile of the difference between azimuthal velocity, vϕ, and circular velocity, vc(R), at R0 = 8.179 kpc. The bottom panel shows the azimuthal surface density profile. The vertical dashed lines are the potential minima of the underlying spiral arm potential. Rotation of the system is from left to right.

5 Discussion

5.1 Implications

From our analysis, we find that kinematic distance estimates are most unreliable close to the Sun and along Galactic longitudes towards the Galactic (anti-)centre. Additionally, there is a higher deviation within interarm regions than within spiral arms. This implies that clouds within underdense regions of the galaxy are more likely to have an incorrect distance estimate from the kinematic distance method.

An observational study carried out by Wenger et al. (2018) compared the distance estimates obtained by the parallax method to those obtained with the kinematic distance method. By treating the parallax distance as the true distance, they found distance deviations of ±40% for their entire sample (see their Table 4). This is in line with the findings reported by Sofue (2011) and our analysis here. As a cautionary note, we reiterate that our simulations do not contain any additional velocity perturbations, such as those that could be driven by stellar feedback (see also Sect. 5.2), so our numbers are likely to be lower bounds. However, to quantify this effect requires more detailed simulations that include all relevant physical processes self-consistently (for a preliminary discussion see e.g. Tress et al. 2024).

To help give an idea as to what implications this has for the observations, we generated a longitude-distance map of our kinematic distance errors in Fig. 17. We overplot some of the sources listed by Wenger et al. (2018) to indicate where real objects would lie on the map, taking the values derived from maser parallax measurements as true distances. We provide a quantitative comparison in Table 5 and find that in some instances the errors are similar, within a factor of two or less. However, in others they are quite different. For example, the error associated with AFGL 2789 has an error on the order of −8% in our system but is around ∼50% from observations. We note that our hydrodynamic model for the response of the gas to the Galactic potential is chosen to be highly idealised, and for example include neither star formation nor stellar feedback (see the discussion in Sect. 5.2). The local gas flow in spiral arms is therefore approximate and we expect more stochastic behaviour in more realistic simulations (see e.g. Tress et al. 2020a, 2021). As AFGL 2789 is located within the Perseus spiral arm (Oh et al. 2010) this could explain the deviation between our estimate and that of Wenger et al. (2018).

Fig. 16 Distributions of the peculiar motions |vp|, LOS velocity difference from circular motion |vlos − vc,losI, and the difference between the kinematic and true distance |dk − dtrue|. The dashed vertical lines on the histograms indicate the 5%, 20%, 80% and 95% percentiles of the distributions. There is a power-law correlation between the peculiar velocity and LOS velocity deviation of slope α = 1.1 (middle-left plot). Similarly, there is a power-law correlation between the LOS deviations from circular velocity and the resulting distance error (bottom-middle plot) with a slope α = 0.94.

Fig. 17 Map of the relative kinematic distance error as a function of Galactic longitude and true distance. The contours are placed at 0, ±0.1, ±0.3, ±0.5, ±1, ±1.5, and ±2. The green lines indicate the spiral arms extracted with FILFINDER. Overplotted in orange are a selection of known sources for which we compute the true distance $d_{\text{true }}^{\text{w}}$ from the parallax measurements of Wenger et al. (2018).

5.2 Caveats

There are a few limitations to bear in mind when considering the maps of kinematic distance errors derived from our simulations. First, as mentioned in Sect. 3, the simulations performed here are 2D dimensional. This corresponds to the assumption that the gas in our simulation is integrated vertically, along the z-axis. The acceleration of the gas due to the potential is computed as if the gas lies in the midplane of the Galaxy (z = 0). This completely neglects the 3D structure of the Galaxy and vertical motions present within the gas. This additional component will impose changes to Eqs. (23) and (24) with the introduction of additional cos (b) terms. With the perturbations induced by the potential, the gas can also experience changes in the z-component of the velocity as it travels in and out of a spiral arm. This will impact the LOS velocity of the gas and the resulting kinematic distance. However, quantifying the size of this effect is beyond the scope of this paper.

Second, the simulations presented here are idealised with only isothermal hydrodynamics and an external galactic potential. We note that the adopted value of c_s = 10 km s⁻¹ implies that we use an effective sound speed, which includes a strong turbulent component providing additional support of the gas disk (Mac Low & Klessen 2004; Krumholz & McKee 2005; McKee & Ostriker 2007; Klessen & Glover 2016). Besides a self-consistent treatment of stellar feedback, more realistic models should include time-dependent chemistry (e.g. Glover et al. 2010). Altogether, we would expect the sound speed to vary across the Galaxy, causing the structure of the Galaxy to alter. Indeed, parameter studies with different sound speeds have shown that spiral arms tend to get wider with increasing cs (Li et al. 2022), whereas the size of the CMZ becomes smaller (Sormani et al. 2015a, 2024).

The spiral arm potential we present in this paper is not an exact match to the spiral arm pattern of the Milky Way. It is well noted that the spiral arm shape of the Milky Way is not regular, with differences in phase angles between spiral arms and pitch angle changes along the spiral arm (see Reid et al. 2016, 2019). These irregularities are hard to model when constructing the potential so approximations are required. Here, we keep the pitch angle of the spiral arms the same and try to account for the change in phase angle with the two pairs of Gaussian shaped spiral arms (see Sect. 2.1.6). Additionally, perturbations generated by the interaction with orbiting satellite galaxies can impact the velocity space of the Galaxy, creating wave-like oscillations throughout the Galaxy (Khanna et al. 2019). These differences in structure will make some of the kinematic distance estimates in our simulation very different from what they would be for the Milky Way, as illustrated in Sect. 5.1.

Turbulent motions induced by physics such as self-gravity, stellar and supernova feedback would contribute to the velocity dispersion of the system. The effects of self-gravity can add an additional ∼2–5 km s⁻¹ to the velocity dispersion for axisymmetric systems but can be as high as ∼10 km s⁻¹ for non- axisymmetric systems such as the one presented in this paper (Wada et al. 2002). Despite this, supernova feedback is believed to give the largest contributions to the velocity dispersion on large scales, potentially producing a velocity dispersion of as much as ∼10 km s⁻¹ across hundreds of parsecs (Lu et al. 2020). This can result in the LOS velocity deviating from the values derived here by a similar amount, causing kinematic distance estimates to deviate further from the true value.

Kinematics distance are normally computed with one of two Galaxy rotation models, the rotation curve of Brand & Blitz (1993) and the universal rotation curve of Persic et al. (1996) with updated parameters from Reid et al. (2014). Both of these rotation curves have been obtained from the gas within the Milky Way; the former making use of HII regions and HI tangent point data, whilst the latter makes use of maser parallaxes. We do see differences between the rotation curves presented in these papers and our values, since our potential is modelled based on the rotation curves from stellar data, as shown in Fig. 14 (Mróz et al. 2019; Eilers et al. 2019). There is a difference on the order of up to ∼10 km s⁻¹ between the rotation curve within the disk of the Galaxy, and larger deviations within the inner most 3 kpc. An investigation into how these differences in the standard rotation curves can impact kinematic distance is beyond the scope of this paper, but does warrant future investigation.

In addition, we note that the simple assumption of the Milky Way being close to dynamical equilibrium requires revisions. Besides the local velocity fluctuations due to stellar feedback, as discussed above, there is increasing evidence that our Galaxy is significantly perturbed on global scales. One of the main results of the Gaia satellite has been to reveal that the Milky Way’s disk is much farther from equilibrium than previously expected. It is subject to strong perturbations from satellite galaxies, such as the Magellanic Clouds and the Sagittarius dSph, that cause it to depart significantly from axisymmetry (see e.g. Antoja et al. 2018, 2022; Widmark et al. 2022b), and generate warps in the outer disk or even lead to a displacement of the baryonic component with respect to the global dark matter distribution (e.g. Laporte et al. 2018b,a; Yaaqib et al. 2024; Chandra et al. 2024). Future estimates of local distance estimates, especially for the outer disk, need to take this into consideration.

6 Conclusions

In this paper, we have presented a realistic analytic potential for the Milky Way, and it contains density profiles for all major mass components of the Galaxy. These include the supermas- sive black hole in the very centre; the nuclear stellar cluster and nuclear stellar disk; the Galactic bar; the Galactic disk, which we split into axisymmetric components for field stars and gas and a spiral arm component for the field stars only; and finally an extended dark matter halo that dominates the potential at large distances. We introduced these components and fitted them to the observational constraints, such as the rotation curve and the terminal velocities, in Sect. 2.

We also described how the new analytic potential is implemented in the moving-mesh code AREPO using the AGAMA framework, outlined in detail in Sect. 3. We made use of 2D hydrodynamic simulations to investigate how robustly the axisymmetric assumption holds for kinematic distance estimates. For this, we placed an observer at R₀ = 8.179 kpc with the bar angled at 28°, generated kinematic distance estimates to each of the gas cells present within our simulations, and computed the systematic errors for each, (d_k − d_true)/d_true. As discussed in Sect. 4, we find that the errors are high close to the Sun, with values reaching greater than 50% on average for any sources with 1 kpc. Along with proximity, we found that errors also reach these values when viewing towards the Galactic centre and anti-centre, l = 0° and l = 180°, respectively. We report peculiar motions associated with the Galactic spiral arm structure of up to 6–8 km s⁻¹, which is in line with the values inferred from maser observations in the Milky Way.

When considering both Galactic longitude and distance, there are certain LOSs that result in higher errors in addition to those previously mentioned for specific distance ranges. We identified these regions as zones of avoidance for the application of the kinematic distance method. We also compared our results with the distance estimates of some well-studied molecular clouds and found generally good agreement (Sect. 5).

In summary, the extraction of the velocity perturbations in our simulation allowed us to determine what impact the potential has on systematic errors in the kinematic distance estimate. We find that within the spiral arms of the Galaxy, the kinematic distance errors are low, as the gas lies within the local potential minima. Consequently, the LOS velocity of the gas is close to what is expected for the axisymmetric version of our potential. We expect clouds within the spiral arms of the Milky Way to have low systematic kinematic distance errors. Conversely, the interarm regions present the largest deviation in both the kinematic distance and the LOS velocity for a given Galactic radius. This is caused by the gas being sped up or slowed down as it travels into or out of the spiral arm. Additionally, we discovered a power-law relation between the systematic kinematic distance error and the difference between the LOS velocity and the projected circular velocity.

We conclude that the assumption of axisymmetry for the kinematic distance method can result in large systematic deviations depending on where a source is situated within the Milky Way. These deviations can alter derived values that depend on distance and, as such, the corresponding systematic errors should be accounted for in these derived values. The interface between AGAMA and AREPO orGADGET-4, along with the scripts for constructing the potential from Sect. 2 and for running N- body simulations with these codes in the external potential, are included in the AGAMA repository².

Acknowledgements

We thank the referee for insightful commments that helped to improve the clarity of this paper. This investigation is funded by the European Research Council under ERC Synergy Grant ECOGAL (grant 855130). The authors gratefully acknowledge computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG. They also acknowledge data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. RSK furthermore acknowledges financial support from the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany’s Excellence Strategy (EXC-2181/1 – 390900948), and from the German Ministry for Economic Affairs and Climate Action in project MAINN (funding ID 50OO2206). The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR project a104bc. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG) – 440719683. MCS acknowledges financial support from the European Research Council under the ERC Starting Grant “GalFlow” (grant 101116226) and from the Royal Society (URF\R1\221118). EV acknowledges support from an STFC Ernest Rutherford fellowship (ST/X004066/1). NB acknowledges support from the ANR BRIDGES grant (ANR-23-CE31-0005). RSK also thanks the Harvard- Smithsonian Center for Astrophysics and the Radcliffe Institute for Advanced Studies for their hospitality during his sabbatical, and the 2024/25 Class of Radcliffe Fellows for highly interesting and stimulating discussions. JG is a member of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD).

Appendix A Spiral arm strength

We perform additional simulations to determine the impact of our choice of the spiral arm strength factor α in Eq. (15). Specifically, we investigate three different values, α = 0.204, 0.408 and 0.612, and find that they lead to a peak density contrast of 10%, 20% and 30% at the solar circle at R₀ = 8.179 kpc, respectively. We caution that Eilers et al. (2020) suggest a density contrast of ~10% for the vicinity of the Sun, and we note that the fiducial case of α = 0.36 introduced in Sect. 2.1.6 gives a density contrast of 17%. However, we also point out that it leads to peculiar gas motions consistent with the maser data (as discussed in Sect. 4.4, see also Reid et al. 2019 or Baba et al. 2009). This is the main reason for our choice.

We also note that the dynamic response ISM to spiral potential perturbations strongly depends on the thermodynamic properties of the medium. Our fiducial value of α = 0.36 gives good results for our highly simplified description of the Galactic ISM, where we adopted an isothermal equation of state, ignored self-gravity in the gas (and consequently neither treated star formation nor stellar feedback), and neglected the potential impact magnetic fields or cosmic ray physics (for details of the numerical set-up, see Sect. 3.1). This put us in the position to study the gas response to the external potential without the need to correct for additional physical processes. We note that the inclusion of appropriate heating and cooling processes as well as star formation and stellar feedback is expected to increase the local deviations from circular velocity and hence the overall velocity dispersion. We speculate that simulations that remain close to the observed level of peculiar gas motions (e.g. Baba et al. 2009; Reid et al. 2019) are likely to require a smaller α and reach spiral arm density contrasts more in line with the values proposed by Eilers et al. (2020). A detailed analysis of these effects will be subject to a follow-up investigation.

Fig. A.1 Density maps of simulations of differing spiral arm strength at t = 441 Myr. The strength of the spiral arm pertibation increase from 10% to 30% stellar disk density at R = 8.179 kpc from top to bottom.

Appendix B Parameter study

We perform a simple parameter study with the potential to understand how varying certain parameters impacts the lv diagram. This involves multiple simulations with the changed parameters. The parameters we consider are number of spiral arms (n), pitch angle of the spiral arms (i), spiral arm pattern speed (Ω_spiral) and bar pattern speed (Ω_bar). The values considered are indicated in Table B.1. We consider two sets of parameter studies, one for each spiral arm number, in which we vary one of the other parameters and keep all other fixed to the fiducial values. Our fiducial parameters are based on the same parameters in Li et al. (2022). This allowed us to investigate how each individual parameter alters the lv diagram and later to vary the parameters to obtain an appropriate approximate match between the observed spiral arms and those in our simulations in lv space.

We present the lv diagrams of our parameter study in Figs. B.1 and B.2. The first thing that stands out between the figures is the number of filamentary structures in the lv diagram away from the centre of the galaxy, most notably the feature that goes between two spiral arms in the region of ~ − 50° to ~ − 130° that is present in the n = 4 set but not in the n = 2 set. At positive galactic longitudes, we find the expected additional spiral arms with the n = 4 set. We find that two of the arms in this region of the lv diagram lie very close to each other for the n = 4, appearing to have split from a spiral arm in the same region in lv space in the n = 2 set.

From the first column of both figures we see that pitch angle has little impact on the shape of the spiral arms in lv space outside the central region. However, within the central 60°, differences can be seen between lv diagrams. Here some features move towards the Galactic centre in lv space with increasing pitch angle though not all features, with those associated with the bar’s rotation remaining fixed in position.

The spiral arm pattern speed has a larger impact as the corresponding resonances for the spiral arms end up changing with pattern speed. This in turn causes the spiral arms to become more apparent with increasing pattern speed beyond the central region of the lv diagram. Additionally, the features associated with the spiral arms tend towards 0 km s⁻¹ in lv space. Within the inner most 60° of the Galaxy, the structure here also moves similarly to how it does with the pitch angle, moving towards the Galactic centre with increasing pattern speed; however, the features here that move are different to those with increasing pitch angle suggesting these are resonance features from the spiral arms.

Similar to the spiral arm pattern speed, changes to bar pattern speed gives arise to different positions for the resonances of the bar, with them moving inward with increasing pattern speed (see Fig. 6). In the lv diagram we see an effect similar to that of the spiral arm pattern speed. Here we see the spiral arms generated by the bar moving towards 0 km s⁻¹ in lv space; however, the broadening of the features does not happen in this case.

It should be noted for all lv diagrams presented in this section have been selected to have approximately the same phase angle between the bar and spiral arms. This means that whilst each snapshot will has similar phase angle between bar and spiral arms, they will be at different stages in evolution. We selected snapshots as close as possible in time late into the systems’ evolution (t > 381 Myr). However, there is a range of 147 Myr within the snapshots selected.

In all cases, the resulting lv diagrams are similar enough that altering one parameter within the constraints of our parameter range does not induce large deviations. However, an exact match to the spiral arm models of McClure-Griffiths et al. (2004), Reid et al. (2014) and Reid et al. (2019) requires more sophisticated modelling that is beyond the scope of this paper.

Fig. B.1 Longitude velocity diagrams of our parameter study with n = 2 spiral arms with similar phase angle between the spiral arm potential and the bar. Left column is where we vary the pitch angle. The middle and right columns are the variation in spiral arm pattern speed and bar pattern speed respectively.

Fig. B.2 Same as Fig. B.1 but with the n = 4 spiral arms subset.

References

All Tables

All Figures

	Fig. 1 Circular-velocity curve produced by our model (black line). The contribution of each component of the potential represented by the coloured lines as detailed in the legend. Observation points from Eilers et al. (2019) and Mróz et al. (2019) are also included (colored markers).
In the text

	Fig. 2 Underlying density distribution within the midplane (z = 0 kpc) of the gravitational potential of the Galaxy.
In the text

	Fig. 3 Strength of the spiral arms at solar circle, R = 8.179 kpc, as a function of azimuth. Shown is the strength of the spiral arms in black as well as the $f_{m_1,{\gamma }_1}\text{ and }{f}_{m_2,{\gamma }_2}$ components in blue and red, respectively.
In the text

	Fig. 4 Enclosed mass profile of our fiducial potential (black) compared to the constraints from dynamical modelling of satellite galaxies (magenta: Correa Magnus & Vasiliev 2022, green: Cautun et al. 2020) and streams (blue: Vasiliev et al. 2021, red: Koposov et al. 2023).
In the text

	Fig. 5 Vertical acceleration at |z| = 1.1 kpc (solid red) and 0.4 kpc (dashed blue line). The observational constraints from Bovy & Rix (2013) in the inner Galaxy at |z| = 1.1 kpc and from Widmark et al. (2022a) in the outer Galaxy at |z| = 0.4 kpc are plotted by red and blue points.
In the text

Fig. 6 Frequency curve as a function of Galactic radius. The solid black line is the rotational frequency curve of the potential, whereas the dashed and dot-dashed lines are rotional frequency plus or minus 0.25 and 0.5 times the epicyclic frequency, ${\kappa }_0(R)=\sqrt{(2\text{Ω}/R)\text{d}\left(R^2\text{Ω}\right)/\text{d}R}$. The horizontal blue and red lines are the pattern speed of the bar and the spiral arms, respectively. The vertical dashed and dot-dashed lines are the corresponding resonances for bar and spiral arms.

In the text

	Fig. 7 Relation between cell size and density within our simulation. The cell size in this case is the size of a square with the same area as the cell.
In the text

Fig. 8 Longitude–velocity maps of the simulation observed from the Sun’s position (placed at the origin of the overlayed coordinate system in Fig. 9). The bar major axis is rotated by 28 degrees from the LOS passing through the Galactic centre. Overlayed blue dashed lines are the spiral arm tracks of Reid et al. (2016, 2019) and McClure-Griffiths et al. (2004). The green points are the terminal velocities from H I and CO observations. Top: full diagram. Bottom: zoom-in between l = −60° & +60°. The left plot only takes into account gas that lies at R < 6 kpc from the centre, and the right one only considers gas at R > 6 kpc.

In the text

Fig. 9 Spiral arm structure. The top plot is the x − y density projection of our simulation at t = 441 Myr. Spiral arms have been enhanced in this image using the mask generated by FILFINDER to brighten overdensities and darken other regions. The bottom panel is the polar decomposition of the density map of the top panel, focusing on the region with galactocentric radius 6 < R < 20 kpc. The blue dotted lines are the positions of the spiral arms according to the second half of Eq. (17). The green dashed lines are the spiral arms generated by the rotation of the bar, with a pitch angle of 6.5° around the outer Lindblad resonance, R = 10.14 kpc. The light blue lines indicate the spiral arm pattern extracted with FILFINDER.

In the text

Fig. 10 Kinematic distance uncertainties at different positions in the disk. The left column presents the LOS velocity map of the simulation and the comparison to the equivalent if the gas follows circular orbits. The top plot illustrates LOS velocity map of the simulation. Similarly, the middle plot shows the map of the LOS velocity derived from the base axisymmetric potential. The last plot shows the difference between simulation and axisymmetric LOS velocities; top plot minus the middle plot. The right column shows the estimated and real distance maps in the simulation. The black lines indicate the spiral arm pattern extracted with FILFINDER. The top plot shows the kinematic distance maps as estimated with Eq. (24). The middle plot is the map of the true distances to the gas cells. The relative error between the kinematic and true distance is shown in the bottom plot.

In the text

Fig. 11 Kinematic distance error at different distances from the Sun. The black line represents the median absolute relative kinematic distance error along the LOS of a given Galactic longitude. Each plot represents a sampling annulus centred on the Sun. The blue shaded region is the median absolute deviation (MAD) of the error. The red shaded regions represent the Galactic longitudes where more than 20% of the cells along the LOS have a relative kinematic distance error of 27% or more.

In the text

	Fig. 12 Probability distributions of the systematic kinematic distance errors within the overdensities in the bar region (blue) and the underdensities in the bar region (red).
In the text

	Fig. 13 Probability distributions of the systematic kinematic distance errors within the spiral arm region (blue) and the interarm regions (red).
In the text

	Fig. 14 Comparison of rotation curve against the average azimuthal velocity (black solid and dashed line respectively). The red and blue lines are the rotation curve from Brand & Blitz (1993) and Reid et al. (2019).
In the text

Fig. 15 Deviations from circular velocity and relation to column density at the solar circle. The top panel illustrates the Galactic azimuthal profile of the difference between azimuthal velocity, vϕ, and circular velocity, vc(R), at R0 = 8.179 kpc. The bottom panel shows the azimuthal surface density profile. The vertical dashed lines are the potential minima of the underlying spiral arm potential. Rotation of the system is from left to right.

In the text

Fig. 16 Distributions of the peculiar motions |vp|, LOS velocity difference from circular motion |vlos − vc,losI, and the difference between the kinematic and true distance |dk − dtrue|. The dashed vertical lines on the histograms indicate the 5%, 20%, 80% and 95% percentiles of the distributions. There is a power-law correlation between the peculiar velocity and LOS velocity deviation of slope α = 1.1 (middle-left plot). Similarly, there is a power-law correlation between the LOS deviations from circular velocity and the resulting distance error (bottom-middle plot) with a slope α = 0.94.

In the text

Fig. 17 Map of the relative kinematic distance error as a function of Galactic longitude and true distance. The contours are placed at 0, ±0.1, ±0.3, ±0.5, ±1, ±1.5, and ±2. The green lines indicate the spiral arms extracted with FILFINDER. Overplotted in orange are a selection of known sources for which we compute the true distance $d_{\text{true }}^{\text{w}}$ from the parallax measurements of Wenger et al. (2018).

In the text

	Fig. A.1 Density maps of simulations of differing spiral arm strength at t = 441 Myr. The strength of the spiral arm pertibation increase from 10% to 30% stellar disk density at R = 8.179 kpc from top to bottom.
In the text

	Fig. B.1 Longitude velocity diagrams of our parameter study with n = 2 spiral arms with similar phase angle between the spiral arm potential and the bar. Left column is where we vary the pitch angle. The middle and right columns are the variation in spiral arm pattern speed and bar pattern speed respectively.
In the text

	Fig. B.2 Same as Fig. B.1 but with the n = 4 spiral arms subset.
In the text