# Testing kinematic distances under a realistic Galactic potential: Investigating systematic errors in the kinematic distance method arising from a non-axisymmetric potential

Glen H. Hunter, Mattia C. Sormani, Jan P. Beckmann, Eugene Vasiliev, Simon C. O. Glover, Ralf S. Klessen, Juan D. Soler, Noé Brucy, Philipp Girichidis, Junia Göller, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Glen H. Hunter, Mattia C. Sormani, Jan P. Beckmann, Eugene Vasiliev, Simon C. O. Glover, et al. Testing kinematic distances under a realistic Galactic potential: Investigating systematic errors in the kinematic distance method arising from a non-axisymmetric potential. Astronomy & Astrophysics - A&A, 2024, 692, 10.1051/0004-6361/202450000. insu-04858881

HAL Id: insu-04858881

https://insu.hal.science/insu-04858881v1

Submitted on 30 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Testing kinematic distances under a realistic Galactic potential

## Investigating systematic errors in the kinematic distance method arising from a non-axisymmetric potential

Glen H. Hunter<fnref n="1" />, Mattia C. Sormani<fnref n="2" />,<fnref n="3" />, Jan P. Beckmann<fnref n="1" />, Eugene Vasiliev<fnref n="4" />,<fnref n="3" />, Simon C. O. Glover<fnref n="1" />, Ralf S. Klessen<fnref n="1" />,<fnref n="5" />,<fnref n="6" />,*,* *<fnref n="1" />, Juan D. Soler<fnref n="7" />, Noé Brucy<fnref n="1" />,<fnref n="8" />, Philipp Girichidis<fnref n="1" />, Junia Göller<fnref n="1" />, Loke Ohlin<fnref n="1" />, Robin Tress<fnref n="9" />, Sergio Molinari<fnref n="7" />, Ortwin Gerhard<fnref n="10" />, Milena Benedettini<fnref n="7" />, Rowan Smith<fnref n="11" />, Patrick Hennebelle<fnref n="12" />, and Leonardo Testi<fnref n="13" />,<fnref n="14" />

Received 15 March 2024 / Accepted 26 October 2024

## ABSTRACT

**Context.** Obtaining reliable distance estimates to gas clouds within the Milky Way is challenging in the absence of certain tracers. The kinematic distance approach has been used as an alternative, and it is derived from the assumption of circular trajectories around the Galactic centre. Consequently, significant errors are expected in regions where gas flow deviates from purely circular motions.

**Aims.** We aim to quantify the systematic errors that arise from the kinematic distance method in the presence of a Galactic potential that is non-axisymmetric. We investigated how these errors differ in certain regions of the Galaxy and how they relate to the underlying dynamics.

**Methods.** We performed 2D isothermal hydrodynamical simulation of the gas disk with the moving-mesh code AREPO, adding the capability of using an external potential provided by the AGAMA library for galactic dynamics. We introduced a new analytic potential of the Milky Way, taking elements from existing models and adjusting parameters to match recent observational constraints.

**Results.** In line with results of previous studies, we report significant errors in the kinematic distance estimate for gas close to the Sun along sight lines towards the Galactic centre and anti-centre and associated with the Galactic bar. Kinematic distance errors are low within the spiral arms, as gas resides close to local potential minima and the resulting line-of-sight velocity is similar to what is expected for an axisymmetric potential. Interarm regions exhibit large deviations at any given Galactic radius, and this is caused by the gas being sped up or slowed down as it travels into or out of spiral arms. In addition, we identify 'zones of avoidance' in the *lv*-diagram, where the kinematic distance method is particularly unreliable and should only be used with caution, and we find a power-law relation between the kinematic distance error and the deviation of the projected line-of-sight velocity from circular motion.

**Key words.** ISM: general – ISM: kinematics and dynamics – ISM: structure – Galaxy: kinematics and dynamics – Galaxy: structure

## 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)

<footnote n="1">¹ Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany</footnote>

<footnote n="2">² Università dell'Insubria, via Valleggio 11, 22100 Como, Italy</footnote>

<footnote n="3">³ Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote n="4">⁴ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, U.K.</footnote>

<footnote n="5">⁵ Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="7">⁷ INAF – Istituto di Astrofisica e Planetologia Spaziali, Via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote n="8">⁸ Centre de Recherche Astrophysique de Lyon UMR5574, ENS de Lyon, Univ. Lyon1, CNRS, Université de Lyon, 69007 Lyon, France</footnote>

<footnote n="9">⁹ Institute of Physics, Laboratory for galaxy evolution and spectral modelling, EPFL, Observatoire de Sauverny, Chemin Pegais 51, 1290 Versoix, Switzerland</footnote>

<footnote n="10">¹⁰ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse, 85748 Garching, Germany</footnote>

<footnote n="11">¹¹ School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK</footnote>

<footnote n="12">¹² Laboratoire AIM, Paris-Saclay, CEA/IRFU/SAp – CNRS – Université Paris Diderot. 91191, Gif-sur-Yvette Cedex, France</footnote>

<footnote n="13">¹³ Dipartimento di Fisica e Astronomia “Augusto Righi”, Viale Berti Pichat 6/2, Bologna, Italy</footnote>

<footnote n="14">¹⁴ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote>* Corresponding author; klessen@uni-heidelberg.de</footnote>

<footnote>** Elizabeth S. and Richard M. Cashin Fellow, the Radcliffe Institute for Advanced Studies at Harvard University, 10 Garden Street, Cambridge, MA 02138, USA</footnote>
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 star-forming 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 \lesssim 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 $\ell \sim 130^\circ$ (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 Way-like 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, $\rho_{\text{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:

$$ \Phi_{\text{SgrA}^*} = -\frac{GM_{\text{SgrA}^*}}{\sqrt{r^2 + b^2}}, \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Underlying density distribution within the midplane ($z=0$ kpc) of the gravitational potential of the Galaxy.</figcaption></figure>

where $G$ is the gravitational constant, the mass $M_{SgrA^*} = 4.154 \times 10^6 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.

### 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}(a + a_0)^{4-\gamma}}, \quad (2) $$

where

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

The parameters $\gamma = 0.71$, $q = 0.73$, $a_0 = 5.9$ pc and $M_{\text{NSC}} = 6.1 \times 10^7 M_{\odot}$ 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], \quad (4) $$

where $a$ is as defined in Eq. (3) but with $q = 0.37$, and where $n_1 = 0.72$, $n_2 = 0.79$, $R_1 = 5.06$ pc, $R_2 = 24.6$ pc, $\rho_1/\rho_2 = 1.311$ and $\rho_2 = 153 \times 10^{10} M_{\odot} \text{ kpc}^{-3}$, 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}} = \underbrace{\rho_{\text{bar},1} + \rho_{\text{bar},2}}_{\text{bar}} + \underbrace{\rho_{\text{bar},3}}_{\text{long bar}}. \quad (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:

$$ \rho_{\text{bar},1}(x, y, z) = \rho_1 \text{sech}(a^m) \times \left[1 + \alpha \left(e^{-a_+^m} + e^{-a_-^m}\right)\right] e^{-\left(\frac{r}{r_{\text{cut}}}\right)^2}, \quad (6) $$
<table><caption>**Table 1.** Parameters for components 2 and 3 of the bar.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Value</th></tr><tr><th>Component 2</th><th>Component 3</th></tr></thead><tbody><tr><td>ρ<sub>i</sub> (M<sub>☉</sub> kpc<sup>-3</sup>)</td><td>0.5 × 10<sup>9</sup></td><td>1.743 × 10<sup>13</sup></td></tr><tr><td>x<sub>i</sub> (kpc)</td><td>5.364</td><td>0.478</td></tr><tr><td>y<sub>i</sub> (kpc)</td><td>0.959</td><td>0.297</td></tr><tr><td>z<sub>i</sub> (kpc)</td><td>0.611</td><td>0.252</td></tr><tr><td>R<sub>i,in</sub> (kpc)</td><td>0.558</td><td>7.607</td></tr><tr><td>R<sub>i,out</sub> (kpc)</td><td>3.19</td><td>2.204</td></tr><tr><td>c<sub>⊥,i</sub></td><td>0.97</td><td>1.879</td></tr><tr><td>n<sub>i,in</sub></td><td>3.196</td><td>1.63</td></tr><tr><td>n<sub>i,out</sub></td><td>16.731</td><td>-27.291</td></tr></tbody></table>

Notes. Parameter definitions are found in the main text.

where

$$a = \left\{ \left[ \left( \frac{|x|}{x_1} \right)^{c_{\perp}} + \left( \frac{|y|}{y_1} \right)^{c_{\perp}} + \left( \frac{|z|}{z_1} \right)^{c_{\parallel}} \right]^{\frac{1}{c_{\parallel}}} \right\}, \quad (7)$$

$$a_{\pm} = \left[ \left( \frac{x \pm cz}{x_c} \right)^2 + \left( \frac{y}{y_c} \right)^2 \right]^{\frac{1}{2}}, \quad (8)$$

$$r = (x^2 + y^2 + z^2)^{\frac{1}{2}}. \quad (9)$$

Here, $\alpha = 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_1 = 0.49$ kpc, $y_1 = 0.392$ kpc, $z_1 = 0.229$ kpc, $x_c = 0.751$ kpc and $y_c = 0.469$ kpc, shape the bar along with the shaping parameters $c_{\perp} = 2.232$ and $c_{\parallel} = 1.991$. The X-shape of the bar trails off with power-law exponents with powers $m = 0.873$ and $n = 1.94$ with an additional cutoff radius at $r_{\text{cut}} = 4.37$ kpc. The density profile is normalised to $\rho_1 = 3.16 \times 10^9$ M$_{\odot}$ kpc$^{-3}$.

The second and third component of $\rho_{\text{bar}}$ describe the ellipsoid shape of the bar, which we split into short ($\rho_2$) and long ($\rho_3$) 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):

$$\rho_{\text{bar},i}(x, y, z) = \rho_i e^{-a_i^{n_i}} \operatorname{sech}^2\left(\frac{z}{z_i}\right) \times e^{-\left(\frac{R}{R_{\text{out}}}\right)^{n_i, \text{out}}} e^{-\left(\frac{R_{\text{in}}}{R}\right)^{n_i, \text{in}}}, \quad (10)$$

where $i = \{2, 3\}$ and

$$a_i = \left[ \left( \frac{|x|}{x_i} \right)^{c_{\perp,i}} + \left( \frac{|y|}{y_i} \right)^{c_{\perp,i}} \right]^{\frac{1}{c_{\perp,i}}}, \quad (11)$$

$$R = (x^2 + y^2)^{\frac{1}{2}}. \quad (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_{\text{bar}} = 1.83 \times 10^{10}$ M$_{\odot}$.

## 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

$$\rho_{\text{disk}}(R, z) = \frac{\Sigma_1}{2h_1} \exp\left(-\frac{R}{R_{d,1}} - \frac{R_{\text{cut}}}{R} - \frac{|z|}{h_1}\right) + \frac{\Sigma_2}{2h_2} \exp\left(-\frac{R}{R_{d,2}} - \frac{R_{\text{cut}}}{R} - \frac{|z|}{h_2}\right), \quad (13)$$

where $\Sigma_1 = 1.3719 \times 10^3$ M$_{\odot}$ pc$^{-2}$, $R_{d,1} = 2$ kpc, $z_1 = 300$ pc, $\Sigma_2 = 9.2391 \times 10^2$ M$_{\odot}$ pc$^{-2}$, $R_{d,2} = 2.8$ kpc, $z_2 = 900$ pc, and $R_{\text{cut}} = 2.4$ kpc. The inner cutoff radius $R_{\text{cut}}$, scale lengths $R_d$ and surface density normalisations $\Sigma$ 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:

$$\rho_{\text{gas}}(R, z) = \frac{\Sigma_1}{4z_1} \exp\left(-\frac{R_{m,1}}{R} - \frac{R}{R_{d,1}}\right) \operatorname{sech}^2\left(\frac{z}{2z_1}\right) + \frac{\Sigma_2}{4z_2} \exp\left(-\frac{R_{m,2}}{R} - \frac{R}{R_{d,2}}\right) \operatorname{sech}^2\left(\frac{z}{2z_2}\right), \quad (14)$$

where $\Sigma_1 = 53.1$ M$_{\odot}$ pc$^{-2}$, $R_{d,1} = 7$ kpc, $z_1 = 85$ pc, and $R_{m,1} = 4$ kpc represents the thick H I disk, whilst $\Sigma_2 = 2.18 \times 10^3$ M$_{\odot}$ pc$^{-2}$, $R_{d,2} = 1.5$ kpc, $z_2 = 45$ pc, and $R_{m,2} = 12$ kpc represents the thinner H$_2$ 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, \phi) = \rho_{\text{disk}}(R, z) \cdot \alpha \frac{R^2}{R_0^2} S(R, \phi). \quad (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 $\alpha = 0.36$ leading to a spiral arm strength of ~17% of the disk density at the solar radius of $R_0 = 8.179$ kpc (GRAVITY Collaboration 2019)<fnref n="1" />. This is somewhat larger than the 10% value suggested by Eilers et al. (2020); however, we

<footnote n="1">¹ We adopt the value suggested by the GRAVITY Collaboration (2019) for comparison with other work in the literature, but we note that there is still some discussion about the exact value with a more recent estimate being $R_0 = 8.277 \pm 0.028$ kpc (see Appendix C in GRAVITY Collaboration 2021).</footnote>
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 $\alpha$ 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 $\alpha$ 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):

$$S(R, \phi) = \sum_{k=1}^{2} \left\{ \exp\left(-\frac{R^2}{\sigma_{\text{sp}}^2}[1 - f_{m_k, \gamma_k}(R, \phi)]\right) - \exp\left(-\frac{R^2}{\sigma_{\text{sp}}^2}I_0\left(-\frac{R^2}{\sigma_{\text{sp}}^2}\right)\right) \right\}, \quad (16)$$

where

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

and $i = 12.5^{\circ}$ and $R_a = 9.64$ kpc, $m_1 = m_2 = 2$, $\gamma_1 = 139.5^{\circ}$ and $\gamma_2 = 69.75^{\circ}$. $\sigma_{\text{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_0$ 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$^{\circ}$ 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]. \quad (18)$$

The density normalisation, $\rho_0$, is determined by using the total mass of an Einasto potential:

$$M = 4\pi\rho_0 a^3 n \Gamma(3n), \quad (19)$$

where the total mass is $M = 1.1 \times 10^{12}$ M$_{\odot}$, the Einasto index is $n = 4.5$ and $\Gamma$ is the gamma function. The scale radius $a$ is related to the half mass radius $r_s$ by $a \approx r_s (3n - 1/3)^{-n}$. 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).

<figure><img src="image_5.png" /><figcaption>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}$ and $f_{m_2, \gamma_2}$ components in blue and red, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>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).</figcaption></figure>

### 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 \partial\Phi_0/\partial 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.
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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 \le R \le 25$ kpc. For coverage within the solar circle, $R < R_0$, 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).

### 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 $\Omega_{bar} = -37.5$ km s$^{-1}$ kpc$^{-1}$ (e.g. Sormani et al. 2015b; Sanders et al. 2019; Li et al. 2022; Clarke & Gerhard 2022) and $\Omega_{spiral} = -22.5$ km s$^{-1}$ kpc$^{-1}$ (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.

<table><caption>**Table 2.** Location of resonances of the non-axisymmetric components of the potential.</caption><thead><tr><th>Resonance</th><th>ILR<br>(kpc)</th><th>In 4:1<br>(kpc)</th><th>CR<br>(kpc)</th><th>Out 4:1<br>(kpc)</th><th>OLR<br>(kpc)</th></tr></thead><tbody><tr><td>Bar</td><td>1.02</td><td>3.45</td><td>6.08</td><td>8.24</td><td>10.14</td></tr><tr><td>Spiral arms</td><td>2.36</td><td>6.50</td><td>10.07</td><td>13.06</td><td>15.89</td></tr></tbody></table>

**Notes.** Notation of resonances: ILR = Inner Lindblad, In 4:1 = inner 4:1, CR = corotation, Out 4:1 = outer 4:1, OLR = outer Lindblad.

<figure><img src="image_8.png" /><figcaption>**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 rotational frequency plus or minus 0.25 and 0.5 times the epicyclic frequency, $\kappa_0(R) = \sqrt{(2\Omega/R)d(R^2\Omega)/dR}$. 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.</figcaption></figure>

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.

## 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 moving-mesh code AREPO (Springel 2010). For isothermal gas in two dimensions, these are

$$ \frac{\partial \Sigma}{\partial t} + \nabla \cdot (\Sigma \mathbf{v}) = 0, \qquad (20) $$

$$ \frac{\partial \Sigma \mathbf{v}}{\partial t} + \nabla \cdot (\Sigma \mathbf{v} \otimes \mathbf{v}) = -\nabla P - \Sigma \nabla \Phi, \qquad (21) $$

where $\Sigma$, $\mathbf{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_s^2 \Sigma$, with the sound speed adopted as $c_s = 10$ km s$^{-1}$. The external potential $\Phi$ is given by our model for the Galactic potential, as
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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_{\text{target}} = 2500 \text{ M}_\odot$). 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$^2$. 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.

## 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

<table><caption>**Table 3.** Expansion used for each component of the potential.</caption><thead><tr><th>Component</th><th>Expansion type</th></tr></thead><tbody><tr><td>Sgr A*</td><td>Multipole</td></tr><tr><td>NSC</td><td>Multipole</td></tr><tr><td>NSD</td><td>Multipole</td></tr><tr><td>Bar</td><td>CylSpline (1)</td></tr><tr><td>Stellar disk</td><td>CylSpline (1)</td></tr><tr><td>Gas disk</td><td>CylSpline (1)</td></tr><tr><td>Spiral arms</td><td>CylSpline (2)</td></tr><tr><td>DM halo</td><td>Multipole</td></tr></tbody></table>

**Notes.** When present, spiral arms are represented by a separate CylSpline potential since they rotate with a different pattern speed than the bar.

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):

$$ \Sigma(R) = \Sigma_1 \exp\left(-\frac{R_{m,1}}{R} - \frac{R}{R_{d,1}}\right) + \Sigma_2 \exp\left(-\frac{R_{m,2}}{R} - \frac{R}{R_{d,2}}\right), \quad (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 \sim 30 \text{ 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)$^2$ box. The mesh is then relaxed with the meshrelax method within Arepo to reach our target mass of 2500 M$_\odot$ 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 \partial\Phi_0/\partial 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 Evolving potential 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_0 = 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_0) = 229$ km s$^{-1}$, 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_0 = 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_0) \sim 250$ km s$^{-1}$ (e.g. Reid & Brunthaler 2020), $\sim 20$ km s$^{-1}$ 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 circular-velocity 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 $x_1$ 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_2$ 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, $\Omega_{bar} = -37.5$ km s$^{-1}$ kpc$^{-1}$ and $\Omega_{spiral} = -22.5$ km s$^{-1}$ kpc$^{-1}$ 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 \sim 11$ kpc and at around the spiral arm crossing point, $R_{cross} = 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_{los}}, \qquad (23)$$

where $l$ is the Galactic longitude of the object, $v_{los}$ is its LOS velocity, and $R_0$ and $v_0$ 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 - (R_0 \sin l)^2}, \qquad (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_{los}^2 > v_{term}^2 - 2v_0(v_{los} - v_{term}) \sin l, \qquad (25)$$

where $v_{term}$ is the terminal velocity along a given Galactic longitude and is given by (Burton & Gordon 1978)

$$|v_{term}| = v_c(R) - v_0 |\sin(l)|. \qquad (26)$$
<figure><img src="image_10.png" /><figcaption>**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^\circ$ & $+60^\circ$. 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.</figcaption></figure>

In other words, the kinematic distance is indefinite when the observed $v_{\text{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_0 \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^\circ$ and $l = 180^\circ$. This arises due to trigonometric effects: as $l$ tends towards $0^\circ$ or $180^\circ$, the value for the Galactic radius, $R$, obtained from Eq. (23) becomes ill-defined as both $\sin(l)$ and $v_{\text{los}}$ tend to 0. These deviations are highlighted in the bottom row of Fig. 10, where we map the difference in LOS velocity 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_1$ orbits.
<figure><img src="image_11.png" /></figure>

**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.

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*<sub>los</sub> 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 dis-
tance 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 con-
sidering both models is 27%. For LOSs and locations where the
systematic error for most points is less than this value, the kine-
matic 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.

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 non-
Gaussian, we look at the difference between quantiles to under-
stand 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 respec-
tively. 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.
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

### 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 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% quantile 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.

<table><caption>**Table 4.** Median and quantile differences of the distributions presented in Figs. 12 and 13.</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">Median</th><th colspan="3">Quantile difference</th></tr><tr><th>20%</th><th>10%</th><th>5%</th></tr></thead><tbody><tr><td>Overdense bar</td><td>0.001</td><td>0.14</td><td>0.3</td><td>0.5</td></tr><tr><td>Underdense bar</td><td>-0.018</td><td>0.22</td><td>0.39</td><td>0.61</td></tr><tr><td>Spiral arms</td><td>0.001</td><td>0.1</td><td>0.17</td><td>0.26</td></tr><tr><td>Interarm</td><td>0.004</td><td>0.1</td><td>0.2</td><td>0.39</td></tr></tbody></table>

**Notes.** The quantile difference indicates the difference between upper and lower percent quantile, for example 20% means the difference between the upper and lower 20% quantile.

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
<figure><img src="image_14.png" /><figcaption>**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).</figcaption></figure>

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 middle-left 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

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Probability distributions of the systematic kinematic distance errors within the spiral arm region (blue) and the interarm regions (red).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**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).</figcaption></figure>

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
<figure><img src="image_17.png" /><figcaption>**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_\phi$, and circular velocity, $v_c(R)$, at $R_0 = 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.</figcaption></figure>

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_0$, 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 galactocentric 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_{\text{los}} - v_{\text{c,los}}|$, and the difference between the true and kinematic distance $|d_k - d_{\text{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 $\alpha = 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 $\alpha = 1.1$. Despite the positive correlations with the LOS velocity deviations, there is no clear correlation between the distance deviation and the peculiar motions.

# 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
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Distributions of the peculiar motions $|v_p|$, LOS velocity difference from circular motion $|v_{los} - v_{c,los}|$, and the difference between the kinematic and true distance $|d_k - d_{true}|$. 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 $\alpha = 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 $\alpha = 0.94$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**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_{true}^w$ from the parallax measurements of Wenger et al. (2018).</figcaption></figure>

explain the deviation between our estimate and that 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,
<table><caption>**Table 5.** Kinematic distances and errors of five well-studied sources in the Galactic disk.</caption><thead><tr><th>Source</th><th><i>l</i><br>(degrees)</th><th><i>d</i><sub>true</sub><sup>w</sup><br>(kpc)</th><th><i>d</i><sub>k</sub><sup>w</sup><br>(kpc)</th><th><i>d</i><sub>k</sub><br>(kpc)</th><th>(<i>d</i><sub>k</sub><sup>w</sup> - <i>d</i><sub>true</sub><sup>w</sup>) / <i>d</i><sub>true</sub><sup>w</sup></th><th>(<i>d</i><sub>k</sub> - <i>d</i><sub>true</sub><sup>w</sup>) / <i>d</i><sub>true</sub><sup>w</sup></th></tr></thead><tbody><tr><td>M17</td><td>15.03</td><td>1.97</td><td>2.33</td><td>1.10</td><td>0.18</td><td>-0.44</td></tr><tr><td>W49N</td><td>44.20</td><td>10.93</td><td>11.52</td><td>11.15</td><td>0.05</td><td>0.02</td></tr><tr><td>NML Cyg</td><td>80.80</td><td>1.6</td><td>1.33</td><td>1.32</td><td>-0.17</td><td>-0.18</td></tr><tr><td>AFGL 2789</td><td>94.60</td><td>3.49</td><td>5.48</td><td>3.21</td><td>0.57</td><td>-0.08</td></tr><tr><td>G240.31+00.07</td><td>-163.84</td><td>7.11</td><td>5.75</td><td>9.17</td><td>-0.19</td><td>0.29</td></tr></tbody></table>

**Notes.** For five Galactic sources (Col. 1), we list the longitude *l* (Col. 2) and the true distance *d*<sub>true</sub><sup>w</sup> (Col. 3) derived from the maser parallax measurements reported by Wenger et al. (2018) as well as their KD estimate *d*<sub>k</sub><sup>w</sup> (Col. 4) for comparison with the value *d*<sub>k</sub> from our model. We also provide the corresponding relative KD errors in columns 6 and 7. Despite the simplicity of the physical processes included in our simulation, we find that our approach reproduces the observational data reasonably well.

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 \text{ km s}^{-1}$ 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 $c_s$ (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<sup>-1</sup> to the velocity dispersion for axisymmetric systems but can be as high as ~10 km s<sup>-1</sup> 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<sup>-1</sup> 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<sup>-1</sup> 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 supermassive 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 imple-
mented 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<sub>k</sub> − d<sub>true</sub>)/d<sub>true</sub>. 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<sup>-1</sup>, 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 poten-
tial 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 poten-
tial minima. Consequently, the LOS velocity of the gas is close
to what is expected for the axisymmetric version of our poten-
tial. 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 dis-
covered 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 or GADGET-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<fnref n="2" />.

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 Ger-
many's Excellence Strategy (EXC-2181/1 – 390900948), and from the Ger-
man Ministry for Economic Affairs and Climate Action in project MAINN
(funding ID 50002206). The authors gratefully acknowledge the scientific
support and HPC resources provided by the Erlangen National High Perfor-
mance 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 par-
tially 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 (URFR1V221118). 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 Rad-
cliffe 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).

References

Ahumada, R. et al. 2020, ApJS, 249, 3

Anderson, L. D., Bania, T. M., Balser, D. S., & Rood, R. T. 2012, ApJ, 754, 62

Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360

Antoja, T., Ramos, P., López-Guitart, F., et al. 2022, A&A, 668, A61

Baba, J., Asaki, Y., Makino, J., et al. 2009, ApJ, 706, 471

Balser, D. S., Wenger, T. V., Anderson, L. D., & Bania, T. M. 2015, ApJ, 806, 199

Barnes, A. T., Watkins, E. J., Meidt, S. E., et al. 2023, ApJ, 944, L22

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)

Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529

Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115

Brand, J., & Blitz, L. 1993, A&A, 275, 67

Burton, W. B., & Gordon, M. A. 1978, A&A, 63, 7

Burton, W. B., & Liszt, H. S. 1993, A&A, 274, 765

Carroll, B. W., & Ostlie, D. A. 2017, An Introduction to Modern Astrophysics, 2nd edn. (Cambridge University Press)

Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291

Chandra, V., Naidu, R. P., Conroy, C., et al. 2024, arXiv e-prints
[arXiv:2406.01676]

Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948

Chemin, L., Renaud, F., & Soubiran, C. 2015, A&A, 578, A14

Clarke, J. P., & Gerhard, O. 2022, MNRAS, 512, 2171

Clemens, D. P. 1985, ApJ, 295, 422

Coleman, B., Paterson, D., Gordon, C., Macias, O., & Ploeg, H. 2020, MNRAS, 495, 3350

Contopoulos, G., & Grosbol, P. 1989, A&A Rev., 1, 261

Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610

Dehnen, W. 1993, MNRAS, 265, 250

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)

Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120

Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2020, ApJ, 900, 186

Einasto, J. 1969, AN, 291, 97

Elia, D., Molinari, S., Schisano, E., et al. 2022, ApJ, 941, 162

Fich, M., Blitz, L., & Stark, A. A. 1989, ApJ, 342, 272

Foster, T., & MacWilliams, J. 2006, ApJ, 644, 214

Freudenreich, H. T. 1998, ApJ, 492, 495

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Glover, S. C. O., Federrath, C., Mac Low, M. M., & Klessen, R. S. 2010, MNRAS, 404, 2

Gómez, G. C. 2006, AJ, 132, 2376

GRAVITY Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10

GRAVITY Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59

GRAVITY Collaboration (Abuter, R., et al.) 2022, A&A, 657, L12

Henshaw, J. D., Barnes, A. T., Battersby, C., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 83

Junqueira, T. C., Lépine, J. R. D., Braga, C. A. S., & Barros, D. A. 2013, A&A, 550, A91

Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864

Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 489, 4962

Klessen, R. S., & Glover, S. C. O. 2016, in Saas-Fee Advanced Course, 43, Saas-Fee Advanced Course, eds. Y. Revaz, P. Jablonka, R. Teyssier, & L. Mayer, 85

<footnote n="2">² https://github.com/GalacticDynamics-Oxford/Agama</footnote>
Koch, E. W., & Rosolowsky, E. W. 2015, MNRAS, 452, 3435

Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135

Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V., & Garavito-Camargo, N. 2018a, MNRAS, 473, 1218

Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018b, MNRAS, 481, 286

Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138

Levine, E. S., Blitz, L., & Heiles, C. 2006, Science, 312, 1773

Li, Z., Shen, J., Gerhard, O., & Clarke, J. P. 2022, ApJ, 925, 71

Lu, Z.-J., Pelkonen, V.-M., Padoan, P., et al. 2020, ApJ, 904, 58

Mac Low, M.-M., & Klessen, R. S. 2004, RMP, 76, 125

Majewski, S. R., et al. 2017, AJ, 154, 94

McClure-Griffiths, N. M., & Dickey, J. M. 2007, ApJ, 671, 427

McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., & Green, A. J. 2004, ApJ, 607, L127

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

McMillan, P. J. 2017, MNRAS, 465, 76

Messineo, M., Habing, H. J., Sjouwerman, L. O., Omont, A., & Menten, K. M. 2002, A&A, 393, 115

Messineo, M., Habing, H. J., Menten, K. M., Omont, A., & Sjouwerman, L. O. 2004, A&A, 418, 103

Messineo, M., Habing, H. J., Menten, K. M., et al. 2005, A&A, 435, 575

Molinari, S., Bally, J., Glover, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 125

Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10

Nakanishi, H., & Sofue, Y. 2003, PASJ, 55, 191

Nakanishi, H., & Sofue, Y. 2006, PASJ, 58, 847

Oh, C. S., Kobayashi, H., Honma, M., et al. 2010, PASJ, 62, 101

Oort, J. H., Kerr, F. J., & Westerhout, G. 1958, MNRAS, 118, 379

Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27

Plummer, H. C. 1911, MNRAS, 71, 460

Pohl, M., Englmaier, P., & Bissantz, N. 2008, ApJ, 677, 283

Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621

Ramón-Fox, F. G., & Bonnell, I. A. 2018, MNRAS, 474, 2028

Reid, M. J. 2022, AJ, 164, 133

Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131

Sanders, J. L., Smith, L., & Evans, N. W. 2019, MNRAS, 488, 4552

Sofue, Y. 2011, PASJ, 63, 813

Sofue, Y. 2023, Ap&SS, 368, 74

Soler, J. D., Miville-Deschênes, M. A., Molinari, S., et al. 2022, A&A, 662, A96

Sormani, M. C., Binney, J., & Magorrian, J. 2015a, MNRAS, 449, 2421

Sormani, M. C., Binney, J., & Magorrian, J. 2015b, MNRAS, 454, 1818

Sormani, M. C., Magorrian, J., Nogueras-Lara, F., et al. 2020, MNRAS, 499, 7

Sormani, M. C., Gerhard, O., Portail, M., Vasiliev, E., & Clarke, J. 2022a, MNRAS, 514, L1

Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022b, MNRAS, 512, 1857

Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, MNRAS, 528, 5742

Springel, V. 2010, MNRAS, 401, 791

Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871

Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge University Press)

Tress, R. G., Smith, R. J., Sormani, M. C., et al. 2020a, MNRAS, 492, 2973

Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020b, MNRAS, 499, 4455

Tress, R. G., Sormani, M. C., Smith, R. J., et al. 2021, MNRAS, 505, 5438

Tress, R. G., Sormani, M. C., Girichidis, P., et al. 2024, A&A, 691, A303

van de Hulst, H. C., Muller, C. A., & Oort, J. H. 1954, Bull. Astron. Inst. Netherlands, 12, 117

Vasiliev, E. 2018, arXiv e-prints [arXiv:1802.08255]

Vasiliev, E. 2019, MNRAS, 482, 1525

Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279

Wada, K., Meurer, G., & Norman, C. A. 2002, ApJ, 577, 197

Watkins, E. J., Kreckel, K., Groves, B., et al. 2023, A&A, 676, A67

Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874

Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050

Wenger, T. V., Balser, D. S., Anderson, L. D., & Bania, T. M. 2018, ApJ, 856, 52

Widmark, A., Laporte, C. F. P., & Monari, G. 2022a, A&A, 663, A15

Widmark, A., Widrow, L. M., & Naik, A. 2022b, A&A, 668, A95

Xu, Y., Reid, M. J., Zheng, X. W., & Menten, K. M. 2006, Science, 311, 54

Yaaqib, R., Petersen, M. S., & Peñarrubia, J. 2024, MNRAS, 531, 3524

Zucker, C., Goodman, A., Alves, J., et al. 2021, ApJ, 919, 35

Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334
# Appendix A: Spiral arm strength

We perform additional simulations to determine the impact of our choice of the spiral arm strength factor $\alpha$ in Eq. (15). Specifically, we investigate three different values, $\alpha = 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_0 = 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 $\alpha = 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 $\alpha = 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 $\alpha$ 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.

# 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 ($\Omega_{\text{spiral}}$) and bar pattern speed ($\Omega_{\text{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

<figure><img src="image_20.png" /><figcaption>**Fig. A.1.** Density maps of simulations of differing spiral arm strength at *t* = 441 Myr. The strength of the spiral arm perturbation increase from 10% to 30% stellar disk density at *R* = 8.179 kpc from top to bottom.</figcaption></figure>

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
<table><caption>**Table B.1.** Parameter values considered.</caption><thead><tr><th>Parameter</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td><i>n</i></td><td>2, 4</td><td></td></tr><tr><td><i>i</i></td><td>10, 12.5*, 15, 17.5</td><td>[°]</td></tr><tr><td>Ω<sub>spiral</sub></td><td>-17.5, -20, -22.5*, -25</td><td>[km s<sup>-1</sup> kpc<sup>-1</sup>]</td></tr><tr><td>Ω<sub>bar</sub></td><td>-36.25, -37.5*, -38.75, -40</td><td>[km s<sup>-1</sup> kpc<sup>-1</sup>]</td></tr></tbody></table>

**Notes.** Parameters as defined in the main text, with * indicating fiducial values.

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 innermost 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 have 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.
<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. B.2. Same as Fig. B.1 but with the *n* = 4 spiral arms subset.</figcaption></figure>