¹²³⁴⁵⁶⁷⁸⁹¹⁰¹¹¹²¹³¹⁴¹⁵

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 , Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Mattia C. Sormani Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Jan P. Beckmann Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Eugene Vasiliev Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Simon C. O. Glover Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Ralf S. Klessen Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Juan D. Soler Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Noé Brucy Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Philipp Girichidis Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Junia Göller Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Loke Ohlin Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Robin Tress Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Sergio Molinari Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Ortwin Gerhard Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Milena Benedettini Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Rowan Smith Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Patrick Hennebelle Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential

Leonardo Testi Testing kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potentialTesting kinematic distances under a realistic Galactic potential (Accepted XXX. Received YYY; in original form ZZZ)

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, 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 investigate how these errors differ in certain regions of the Galaxy and how they relate to the underlying dynamics.

Methods. We perform 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 introduce 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. 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.:

The Galaxy – Galaxy: kinematics and dynamics – Galaxy: structure – ISM: kinematics and dynamics

\nolinenumbers

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 et al. 2023), obtaining reliable distance estimates for the gaseous component, i.e. for the various phases of interstellar medium (ISM, see e.g., Tielens, 2005; Draine, 2011) is much more challenging. Estimating distances to molecular clouds is important for understanding their properties, formation and evolution (Molinari et al., 2014) and their ability to form stars (Klessen & Glover, 2016). The same is true for the atomic and ionised components of the ISM.

Accurate 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 kpc 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 kpc 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 for their sample that 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); (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⁻¹, representative of turbulent motions in giant molecular clouds, can lead to significant (>10%) errors in the KD estimates for true distances d≲ 5 kpc, and can also lead to systematic biases of ∼ 20% despite the random motions having zero mean around the underlying circular motions. A particularly extreme example of deviations from circular velocity leading to errors in the kinematic distance estimate is reported for the Perseus Arm at Galactic longitudes $\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 + 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, 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 construct 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 run 2D hydrodynamical simulations using this potential. We then compare 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 Section 1, we introduce our new analytic description of the Milky Way potential and its various components in Section 2. We briefly describe our numerical approach and the implementation of the new potential in the Agama frameworks combined with Arepo in Section 3. Our main findings are presented in Section 4, and their implications and limitations discussed in Section 5. Finally, we summarise and conclude in Section 6.

Figure 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).

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 Figure 1, and the total midplane density profile of our potential, ρ_gal, is shown in Figure 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_{\mathrm{SgrA}^{\star}}=-\frac{GM_{\mathrm{SgrA}^{\star}}}{\sqrt{r^{2}+b^{2}}}\;,$$

(1)

where G is the gravitational constant, the mass $M_{\mathrm{SgrA}^{\star}}=4.154\times 10^{6}\,\mathrm{M}_{\odot}$ is taken from the GRAVITY Collaboration et al. (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 pc. In our model, the nuclear star cluster (NSC) follows a flattened Dehnen (1993) density profile as given in Chatzopoulos et al. (2015, see their Equation 17):

$$\rho_{\mathrm{NSC}}=\frac{(3-\gamma)M_{\mathrm{NSC}}}{4\pi q}\frac{a_{0}}{a^{\gamma}(a+a_{0})^{4-\gamma}}\;,$$

(2)

where

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

(3)

The parameters γ=0.71, q=0.73, a₀=5.9 pc and M_NSC=6.1× 10⁷M⊙ are taken from their best-fitting model. Here R refers to the Galactic radius in cylindrical coordinates. We note from our circular-velocity curve (Figure 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_{\mathrm{NSD}}=\rho_{1}\exp{\Big[ -\Big( \frac{a}{R_{1}}\Big) ^{n_{1}}\Big] }+\rho_{2}\exp{\Big[ -\Big( \frac{a}{R_{2}}\Big) ^{n_{2}}\Big] }$$

(4)

where a is as defined in Equation (3) but with q=0.37, and where n₁=0.72, n₂=0.79, R₁=5.06 pc, R₂=24.6 pc, ρ₁/ρ₂=1.311 and ρ₂=153× 10¹⁰ M⊙ kpc-3, which follows model 3 of Sormani et al. (2020). As shown in Figure 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_{\mathrm{bar}}=\underbrace{\rho_{\mathrm{bar},1}+\rho_{\mathrm{bar},2}}_{\text{bar}}+\underbrace{\rho_{\mathrm{bar},3}}_{\text{long bar}}\,.$$

(5)

The first component of Equation (5) corresponds to the X-shaped component of the observed boxy-peanut shape of the bar (Wegg & Gerhard, 2013). To describe this, we use a modified form of Equations 9 and 10 of Coleman et al. (2020) and Freudenreich (1998), respectively:

$$\rho_{\mathrm{bar},1}(x,y,z) =\rho_{1}\mathrm{sech}\left(a^{m}\right)\qquad \times\left[1+\alpha\left(e^{-{a_{+}^{n}}}+e^{-{a_{-}^{n}}}\right)\right]e^{-\left(\frac{r}{r_{\mathrm{cut}}}\right)^{2}}\,,$$

(6)

where

$$a =\left\lbrace \left[\left(\frac{|x|}{x_{1}}\right)^{c_{\perp}}+\left(\frac{|y|}{y_{1}}\right)^{c_{\perp}}\right]^{\frac{c_{\parallel}}{c_{\perp}}}+\left(\frac{|z|}{z_{1}}\right)^{c_{\parallel}}\right\rbrace ^{\frac{1}{c_{\parallel}}}\,,\qquad a_{\pm} =\left[\left(\frac{x\pm cz}{x_{c}}\right)^{2}+\left(\frac{y}{y_{c}}\right)^{2}\right]^{\frac{1}{2}}\,,\qquad r =\left(x^{2}+y^{2}+z^{2}\right)^{\frac{1}{2}}\,.$$

(7) (8) (9)

Here, α=0.626 defines the strength of the X-shape whilst c=1.342 defines the slope of the X-shape in the x-z plane. The scale lengths x₁=0.49 kpc, y₁=0.392 kpc, z₁=0.229 kpc, x_c=0.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 additonal cutoff radius at r_cut=4.37 kpc. The density profile is normalised to ρ₁=3.16× 10⁹ M⊙ kpc-3.

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

The second and third component of ρ_bar describe the ellipsoid shape of the bar, which we split into short (ρ₂) and long (ρ₃) bar components. Both components follow a modified version of Equation 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_{\mathrm{bar},i}(x,y,z) =\rho_{i}\,e^{-a_{i}^{n_{i}}}\mathrm{sech}^{2}\left(\frac{z}{z_{i}}\right)\qquad \times e^{-\left(\frac{R}{R_{\mathrm{i},out}}\right)^{n_{\mathrm{i},out}}}e^{-\left(\frac{R_{\mathrm{i},in}}{R}\right)^{n_{\mathrm{i},in}}}\,,$$

(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}}}\,,\qquad R =\left(x^{2}+y^{z}\right)^{\frac{1}{2}}\,.$$

(11) (12)

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

Table 1: Parameters for components 2 and 3 of the bar.

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 adopt a modified version of Equation 3 of McMillan (2017) using an exponential vertical profile. We obtain

$$\rho_{\mathrm{disk}}(R,z) =\frac{\Sigma_{1}}{2h_{1}}\exp\Big( -\frac{R}{R_{\mathrm{d},1}}-\frac{R_{\mathrm{cut}}}{R}-\frac{|z|}{h_{1}}\Big)\qquad {}+\frac{\Sigma_{2}}{2h_{2}}\exp\Big( -\frac{R}{R_{\mathrm{d},2}}-\frac{R_{\mathrm{cut}}}{R}-\frac{|z|}{h_{2}}\Big) \;,$$

(13)

where Σ₁=1.3719× 10³ M⊙ pc-2, R_[d,1]=2 kpc, z₁=300 pc, Σ₂=9.2391× 10² M⊙ pc-2, R_[d,2]=2.8 kpc, z₂=900 pc, and R_cut=2.4 kpc. The inner cutoff radius R_cut, scale lengths R_d and surface density normalizations Σ 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 Figure 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=4kpc)∼ 180pc 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 include two gas disks in the potential, which we take from McMillan (2017) without any further adjustments:

$$\rho_{\mathrm{gas}}(R,z) =\frac{\Sigma_{1}}{4z_{1}}\exp\bigg( -\frac{R_{m,1}}{R}-\frac{R}{R_{d,1}}\bigg) \mathrm{sech}^{2}(z/2z_{1})\qquad {}+\frac{\Sigma_{2}}{4z_{2}}\exp\bigg( -\frac{R_{m,2}}{R}-\frac{R}{R_{d,2}}\bigg) \mathrm{sech}^{2}(z/2z_{2})\;,$$

(14)

where Σ₁=53.1 M⊙ pc-2, R_[d,1]=7 kpc, z₁=85 pc, and Rₘ,₁=4 kpc represents the thick H i disk, whilst Σ₂=2.18× 10³ M⊙ pc-2, R_[d,2]=1.5 kpc, z₂=45 pc, and Rₘ,₂=12 kpc represents the thinner H2 disk. 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 introduce a perturbation to the stellar disk in the following manner:

$$\rho_{\mathrm{spiral}}(R,z,\phi)=\rho_{\mathrm{disk}}(R,z)\cdot\alpha\frac{R^{2}}{R_{0}^{2}}S(R,\phi).$$

(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 α=0.36 leading to a spiral arm strength of ∼ 17% of the disk density at the solar radius of R₀=8.179 kpc (GRAVITY Collaboration et al., 2019).² This is somewhat larger than the 10% value suggested by Eilers et al. (2020), however, we note that we measure the dynamic response in the gas rather than the field star population, which is expected to be stronger due to the more dissipative nature of the ISM. We also note that we expect the impact of the choice of α to sensitively depend on the actual implementation of the thermal physics and chemistry in the code (see also Appendix A). As outlined in Section 3, we here adopt a simple isothermal equation of state for the galactic gas. Consequently, we suggest to explore different values of α when accounting for the true multi-phase nature of the Galactic ISM and its multitude of heating and cooling processes (see e.g. Tielens 2005 or Draine 2011).

We adopt the value suggested by the GRAVITY Collaboration et al. (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₀=8.277± 0.028 kpc (see Appendix C in GRAVITY Collaboration et al., 2021).

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 take a modified form of Equation 8 of Junqueira et al. (2013):

$$S(R,\phi)=\sum_{k=1}^{2}\bigg\lbrace \exp\left(-\frac{R^{2}}{\sigma_{\mathrm{sp}}^{2}}\big[ 1-f_{m_{k},\gamma_{k}}(R,\phi)\big] \right)\quad -\exp{\left(-\frac{R^{2}}{\sigma_{\mathrm{sp}}^{2}}\right)}I_{0}\left(-\frac{R^{2}}{\sigma_{\mathrm{sp}}^{2}}\right)\bigg\rbrace ,$$

(16)

where

$$f_{m,\gamma}(R,\phi)=\cos\left(m(\phi+\gamma)-\frac{m}{\tan(i)}\ln\left(\frac{R}{R_{a}}\right)\right)\;,$$

(17)

and $i=12.5^{\circ}$ and Rₐ=9.64 kpc, m₁=m₂=2, $\gamma_{1}=139.5^{\circ}$ and $\gamma_{2}=69.75^{\circ}$. σₛₚ=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 Equation (16) is used to normalise the spiral arm potential such that the monopole component is zero. Here, I₀ is the modified Bessel function of the first kind and of zeroth order. This resulting potential is a superposition of two pairs of m=2 spiral arms with equal amplitude (Li et al., 2022). We note that the spiral arm potential does not have an m=4 pattern due to the unequal angular separation between spiral arms. We opt for this kind of pattern as it allows for the angular separation to be adjusted as needed as 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 Figure 3.

Figure 3: The 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₁,γ₁ and f_m₂,γ₂ components in blue and red respectively.

2.1.7 Dark matter halo

The dark matter halo component follows a spherical Einasto (1969) profile:

$$\rho_{\mathrm{dm}}(r)=\rho_{0}\exp\bigg[ -\bigg( \frac{r}{a}\bigg) ^{1/n}\bigg] .$$

(18)

The density normalisation, ρ₀, is determined by using the total mass of an Einasto potential:

$$M=4\pi\,\rho_{0}\,a^{3}\,n\,\Gamma(3n)\;,$$

(19)

where the total mass is M=1.1× 10¹² M⊙, the Einasto index is n=4.5 and Γ is the gamma function. The scale radius a is related to the half mass radius rₛ by a≈ rₛ (3n-1/3)⁻ⁿ. In this case, the half mass radius is rₛ=96kpc 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 (Cautun et al., 2020; Correa Magnus & Vasiliev, 2022; Vasiliev et al., 2021; Koposov et al., 2023) (see Figure 4).

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

2.2 Comparison with observations

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

2.2.1 Axisymmetric components

We begin with the Galactic circular-velocity curve, as illustrated in Figure 1 for different radial bins. The black line represents our total circular-velocity curve from the axisymmetrised potential, v_c=(R ∂Φ₀/∂ R)¹⁄², 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 Section 3.2). The resulting circular-velocity curve does not contain perturbations from the spiral arms nor contains the higher order terms needed to describe the full potential of the bar.

We make use of recent measurements of circular velocity data from Eilers et al. (2019), obtained from red giant star observed with APOGEE, WISE and Gaia, and from Mróz et al. (2019), obtained from Cepheid variable stars with Gaia. Both are in a reasonable agreement with each other and provide a coverage of Galactocentric radius of 4≲ R≲ 25 kpc. For coverage within the solar circle, R<R₀, we make of use of the terminal velocities measurements from H i and CO observations (Clemens, 1985; Fich et al., 1989; Burton & Liszt, 1993; McClure-Griffiths & Dickey, 2007). We compare the peaks of the resulting longitude-velocity (lv) diagrams from our hydrodynamical simulations with the corresponding terminal velocity measurements of the Milky Way, as discussed in detail in Section 4. Here we opt to compare terminal velocities in the lv diagram instead of the circular-velocity/rotation 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 ensure 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 Figure 5).

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

2.2.2 Non-axisymmetric components

Our potential has two rotating non-axisymmetric components, the bar and the spiral arms. We chose the pattern speeds of the bar and the spiral arms as Ω_bar=-37.5 km s-1 kpc-1 (e.g. Sormani et al., 2015b; Sanders et al., 2019; Li et al., 2022; Clarke & Gerhard, 2022) and Ωₛₚᵢᵣₐₗ=-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 Figure 6. We find that for the pattern speeds we use, the outer Lindblad resonance of the bar coincides with corotation of the spiral arms at ∼ 10.1 kpc. The outer 4:1 resonance of the bar lies close to solar circle at ∼ 8.2 kpc.

Table 2: Location of resonances of the non-axisymmetric components of the potential.

Figure 6: Frequency curve as a function of Galactic radius. The solid black line is the rotational frequency curve of the potential, whereas the dashed and dot-dashed lines are rotional frequency plus or minus 0.25 and 0.5 times the epicyclic frequency, κ₀(R)=√((2Ω/R)d(R² Ω)/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.

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 solve 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 \frac{\partial\Sigma\mathbf{v}}{\partial t}+\nabla\cdot(\Sigma\mathbf{v}\otimes\mathbf{v})=-\nabla P-\Sigma\nabla\Phi\;,$$

(20) (21)

where Σ, v, and P are gas surface density, velocity and pressure, respectively. The simulations are two-dimensional. The pressure is related to the density via the equation of state, P=cₛ² Σ, with the sound speed adopted as cₛ=10km s-1. The external potential Φ is given by our model for the Galactic potential, as explained in the next section; for ease of interpretation, we do not include the gas self-gravity, star formation or stellar feedback in our models. 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 modeling 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ₜarget=2500 M⊙). The cell will either split or merge with another with addition or removal of a mesh generating point. For further details of the code base we use, see e.g. Tress et al. (2020a, b). We achieve a minimum value of ∼ 40 pc. Figure 7 illustrates how our cellsize varies with density. The vertical spread seen at low cellsizes is a result of the minimum surface area of the cell being reached, which is set to 6 pc2. 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.

Figure 7: The relation between cellsize and density within our simulation. The cellsize in this case is the size of a square with the same area as the cell.

3.2 External Galactic potential

Arepo allows one to include an external gravitational potential in the simulation, but the Galactic potential described in Section 2 is significantly more complex than the few built-in analytic models. We calculate 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 and/or pattern speed). We created an interface between the two codes that make it possible to use any potential implemented in Agama as an external potential in Arepo (in addition to self-gravity of the simulated system, if the latter is turned on). Moreover, a very similar interface is provided for the Gadget-4 code (Springel et al., 2021), which shares a common ancestry with Arepo; both interfaces, as well as the script for generating the Galactic potential from this study, are available in the latest version of Agama.

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

Table 3: Expansion used for each component of the potential.

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 Equation (2.1.5):

$$\Sigma(R) =\Sigma_{1}\exp\bigg( -\frac{R_{m,1}}{R}-\frac{R}{R_{d,1}}\bigg)\qquad {}+\Sigma_{2}\exp\bigg( -\frac{R_{m,2}}{R}-\frac{R}{R_{d,2}}\bigg) \;,$$

(22)

where the parameters for this gas disk are the same as those used for the gas disk potential (See Section 2.1.5). We extend this gas disk to R∼ 30 kpc at which point we reduce the density significantly to prevent artifacts caused by periodic boundary conditions compromising the Galaxy itself. For simplicity, we initialise the simulation with 250000 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 2500M⊙ in which cells are refined or derefined according to the refinement/derefinement criterion mentioned previously. There is no hydrodynamics present in the meshrelax process and as such the gas is fixed until the process is complete.

The velocity of the gas is initialised to be the circular-velocity of the axisymmetric terms of the potential, v_c=(R ∂Φ₀/∂ R)¹⁄², which follows the same circular-velocity curve as Figure 1. The non-axisymmetric components of the bar and spiral arms are introduced linearly and gradually over the course of 150Myr 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₀=8.179kpc (GRAVITY Collaboration et al., 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ₓ=v_c(R₀)=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 Section 2. It should be noted that the GRAVITY Collaboration et al. (2022) has given a more recent value of R₀=8.277 kpc, however we adopt the 2019 value for the calculations in this paper. Additionally, proper motion measurements of Sgr A∗ give v_c(R₀)∼ 250 km s-1 (e.g. Reid & Brunthaler, 2020), ∼ 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 Figure 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 Figure 9.

Figure 8: Longitude-velocity maps of the simulation observed from the Sun’s position (placed at the origin of the overlayed coordinate system in Figure 9). The bar major axis is rotated by 28 degrees from the line of sight 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}$ with the left plot only taking gas into account that lies at R<6kpc from the centre and with the right one only considering gas at R>6kpc.

From Figure 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|\gt 50^{\circ}$). Towards the Galactic centre, comparing the spiral arms becomes difficult due to the perturbations generated by the bar. Some features are 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^{\circ}\lt l\lt 30^{\circ}$ and the other at ${}-20^{\circ}\lt l\lt -10^{\circ}$. 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₁ orbits, a family of orbits elongated parallel to the major axis of the bar (Contopoulos & Grosbol, 1989), until it is shocked at the end of the bar, after which it flows inwards on nearly radial orbits. Eventually, the gas stabilises onto x₂ orbits, forming a ring of material at 220 pc from the centre of the Galaxy. This ring is the equivalent of the Central Molecular Zone (CMZ) and is consistent with the larger end of estimates for the Milky Way’s CMZ (e.g. Henshaw et al., 2023).

Outside of the bar region (R>5 kpc), the gas forms a clear spiral pattern. It is rather complex and has two main components: a two-arm spiral caused by the rotation of the bar, and the four-arm structure created by the spiral component of the potential described in Section 2. These two pattern rotate at different angular speeds, Ω_bar=-37.5 km s-1 kpc-1 and Ωₛₚᵢᵣₐₗ=-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 Figure 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^{\circ}$ near the outer Lindblad resonance of the bar. We note that the pitch angle value of the bar-generated pattern depends on the sound speed of the gas, as can be understood from the dispersion relation of spiral density waves in the tight-winding limit (e.g. Binney & Tremaine, 2008). Whilst there are linear trends in the plot, there are deviations from the underlying spiral arm structure. These occur at the point where the two components intercept at R∼ 11 kpc and at around the spiral arm crossing point, R_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 Figure 9 and we use the masks to contrast the density of the simulation in the x-y projection (Figure 9, top).

Figure 9: The top plot is the x-y density projection of our simulation at t=441 Myr. Spiral arms are enhanced in this image using the mask generated by filfinder to brighten overdensities and darken other regions. The bottom figure is the polar decomposition of the density map of the top figure 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 Equation (17). The green dashed lines are the spiral arms generated by the rotation of the bar, with a pitch angle of ${}6.5^{\circ}$ around the outer Lindblad resonance, R=10.14 kpc. The light blue lines indicate the spiral arm pattern extracted with filfinder.

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_{\mathrm{los}}}\;,$$

(23)

where l is the Galactic longitude of the object, v_los is its line-of-sight velocity, and R₀ and v₀ are the Galactic radius and circular-velocity of the Sun, respectively. For consistency, we use the circular-velocity curve generated by our potential as our rotation curve, as shown in Figure 1. Because R appears on both sides of Equation (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}}\;,$$

(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_{\mathrm{los}}^{2}\gt v_{\mathrm{term}}^{2}-2v_{0}(v_{\mathrm{los}}-v_{\mathrm{term}})\sin l\;,$$

(25)

where v_term is the terminal velocity along a given Galactic longitude and is given by: (Burton & Gordon, 1978)

$$|v_{\mathrm{term}}|=v_{c}(R)-v_{0}|\sin(l)|.$$

(26)

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

Equation (24) can give two answers when observing inside the solar circle, resulting in the well-known kinematic distance ambiguity.

For the sake of simplicity, we resolve the ambiguity by selecting the kinematic distance closest to the true distance value, but we note that this may be difficult to decide in real observations. In Figure 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 Figure 10, we plot the line-of-sight velocity measured by an observer at the solar radius (left) and the KD estimates based on Equation (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 Figure 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 Equation (23) becomes ill-defined as both sin(l) and v_los tend to 0. These deviations are highlighted in the bottom row of Figure 10, where we map the difference in line-of-sight velicity and the relative error between the kinematic distance and the true distance. What stands out is the large relative error in the solar neighbourhood. For objects within 500 pc of the Sun, the kinematic distance estimate is highly unreliable and exhibits a broad error distribution. Large systematic errors also occur close to the spiral arm perturbations and at the end of the bar, where the gas flows on x₁ orbits.

Figure 10: The left column presents the LOS velocity map of the simulation and the comparison to the equivalent if the gas followes 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.

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

To give a more conclusive idea about which lines of sight to avoid, we compute the absolute kinematic distance error and plot the median and absolute median deviation (MAD) as a function of Galactic longitude in Figure 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 a line of sight of avoidance for an annulus as the line of sight where more than 20% of cells have an absolute kinematic distance error greater than 27%. Our choice of this value is motivated by the study of Wenger et al. (2018), who quantify the uncertainty in the kinematic distance inferred using the Brand & Blitz (1993) rotation curve (their method A) or the Reid et al. (2014) rotation curve (their method B) due to uncertainties in e.g. the solar Galactocentric radius and orbital velocity, the measured rotation curve, etc. The average absolute kinematic distance error they find when considering both models is 27%. For lines of sight and locations where the systematic error for most points is less than this value, the kinematic distance method should be reliable. On the other hand, if a large fraction of points have systematic errors that exceed this value, this is a good indication that the kinematic distance method will not provide reliable results. The lines of sight that should be avoided are shown in red in Figure 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 5kpc, but drops to 14% and 16% for the 5-10kpc annulus and the 10-20kpc annulus, respectively. This suggests that the kinematic distance estimate is accurate for distances beyond 5kpc from the Sun under our criterion for line of sight avoidance, but that for closer distances it should be used with great care.

Figure 11: The black line represents the median absolute relative kinematic distance error along the line-of-sight 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 line-of-sight have a relative kinematic distance error of 27% or more.

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 lines of sight for objects close to the Sun. This now poses the question: where can one reliably use kinematic distances?

Back in Figure 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, this only applies to the spiral arm features, i.e. peaks outside of the bar region, from inspection. To further analyse this, we split the Galaxy into two regions: the bar region (R<6kpc) and the disk region (R≥6kpc).

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 (Figure 10), splitting the data into overdense regions and underdense regions. We plot the PDF of each region, respectively, and compare their distributions (Figure 12).

Figure 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).

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 understand the width of the distributions. The difference between the upper and lower 20% quantiles is 0.14 for overdense region, whilst it is wider with a value of 0.22 for the underdense region. Similarly, the values are 0.29 and 0.39 for the upper and lower 10% quantiles, and 0.50 and 0.61 for the 5% quantiles respectively. This suggests that outside of the overdense regions there is a higher probability of a large error and, by consequence, an increased probability of obtaining an incorrect distance via the kinematic distance method. Going from our analysis of the error as function Galactic longitude (Figure 11), this result is not too much of a surprise given that much of the bar’s influence is lies within $l=\pm 30^{\circ}$ which is typically a line-of-sight of avoidance for all distances away from the Sun.

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=6kpc.

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

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

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

Table 4: The median and quantile differences of the distributions presented in Figures 12 and 13.

4.4 Relation to the velocity deviation

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

The deviations from rotation curve can be seen in the radial profile of the azimuthal velocity of the gas. We illustrate this in Figure 14 where in the inner most 3 kpc we observed deviations up to 70 km s-1. Between 3 kpc and 10 kpc there are small deviations from the rotation curve, on the order of few km s-1, 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.

Figure 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).

We compute the line-of-sight velocity of our simulations and map it to the face down Milky Way view at the top of the left column of Figure 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 Figure 8), with perturbations of the spiral arms appearing as displacements in the contours in comparison to the middle-left plot of Figure 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 line-of-sight velocity towards the outer Galaxy are smoother. When we subtract the two maps, we find that the largest deviations occur at the bar region with differences on the order of 100km s-1. Outside the bar region we find the line-of-sight 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 Figure 15. This can cause shocks causing the gas to get denser as it leaves the potential minima, as in the case for two of our spiral arms. Additionally, the perturbations caused by the bar can also shock the gas and causes further perturbations. In our case, at around R₀, both the bar generated spiral and two of the spiral arms overlap creating large peaks in density. However, unlike the other two spiral arms, the peak in density for these occur just before passing the potential minima of the spiral arm potential. Other more local sources of velocity perturbations come from expanding HII regions and supernova driven bubbles (see e.g. Barnes et al. 2023 or Watkins et al. 2023 for nearby galaxies, or Zucker et al. 2022 for the solar neighborhood 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 right-hand side of Figure 10). To be more quantitative, peculiar motions associated with the Galactic spiral arm structure of ∼ 6km s-1 in radial and ∼ 4km s-1 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-8km s-1 illustrated in Figure 15.

Figure 15: The top figure illustrates the Galactic azimuthal profile of the difference between azimuthal velocity, v_ϕ and circular velocity from circular-velocity curve, v_c(R), near the solar circle at R₀=8.179 kpc. The bottom figure 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.

Given how large deviations of the line-of-sight 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ᵣ and the distance error. Any deviation from circular motion is likely to occur due to motions in the radial direction. In Figure 16 we provide the distribution and correlations between absolute values of the peculiar motion |vₚ|, the deviations from circular velocity |vₗₒₛ-v_[c,los]|, and the difference between the true and kinematic distance |dₖ-dₜᵣᵤₑ|, split between arm and interarm regions. There is a slight offset between the distributions for each region, with spiral arms having lower deviations and peculiar motions. As expected, we find a clear positive correlation between the local deviations from circular velocity and the resulting distance error. It can be described by a power-law with slope α=0.94. However, we note that there are several orthogonal spurs in the 2D distribution that do not follow this trend. Additionally, we observe a power law correlation between the peculiar velocities and the difference in LOS velocity. Here, the slope is α=1.1. Despite the positive correlations with the LOS velocity deviations, there is no clear correlation between the distance deviation and the peculiar motions.

Figure 16: Distributions of the peculiar motions |vₚ|, LOS velocity difference from circular motion |vₗₒₛ-v_[c,los]|, and the difference between the kinematic and true distance |dₖ-dₜᵣᵤₑ|. The dashed vertical lines on the histograms indicate the 5%, 20%, 80% and 95% percentiles of the distributions. There is a power-law correlation between the peculiar velocity and LOS velocity deviation of slope α=1.1 (middle-left plot). Similarly, there is a power-law correlation between the LOS deviations from circular velocity and the resulting distance error (bottom-middle plot) with a slope α=0.94.

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 could be driven by stellar feedback (see also Section 5.2), and 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 generate a longitude-distance map of our kinematic distance errors in Figure 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 Section 5.2 below). The local gas flow in spiral arms is therefore approximate and we expect more stochastic behavior in more realistic simulations (see e.g. Tress et al., 2020a, 2021). As AFGL 2789 is located within the Perseus spiral arm (Oh et al., 2010) this could explain the deviation between our estimate and that of Wenger et al. (2018).

Table 5: Kinematic distances and errors of five well-studied sources in the Galactic disk.

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 Section 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 Equations (23) and (24) with the introduction of additional cos(b) terms. With the perturbations induced by the potential, the gas can also experience changes in the z-component of the velocity as it travels in and out of a spiral arm. This will impact the LOS velocity of the gas and the resulting kinematic distance. However, quantifying the size of this effect is beyond the scope of this paper.

Figure 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_{\mathrm{true}}^{\,\sf w}$ from the parallax measurements of Wenger et al. (2018).

Second, the simulations presented here are idealised with only isothermal hydrodynamics and an external galactic potential. We note that the adopted value of cₛ=10km 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 Section 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 Section 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-1 to the velocity dispersion for axisymmetric systems but can be as high as ∼ 10 km s-1 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-1 across hundreds of parsecs Lu et al. (2020). This can result in the line-of-sight 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 Figure 14 (Mróz et al., 2019; Eilers et al., 2019). There is a difference on the order of up to ∼ 10 km s-1 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. It contains density profiles for all major mass component of the Galaxy. These are 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. These are introduced and fitted to the observational constraints, such as the rotation curve and the terminal velocities, in Section 2.

We also described how the new analytic potential is implemented in the moving-mesh code Arepo using the Agama framework, as outlined in detail in Section 3. We made use of 2D hydrodynamic simulations to investigate how robustly the axisymmetric assumption holds for kinematic distance estimates. For this, we place an observer at R₀=8.179 kpc with the bar angled at 28∘, generate kinematic distance estimates to each of the gas cells present within our simulations and compute the systematic errors for each, (dₖ-dₜᵣᵤₑ)/dₜᵣᵤₑ. As discussed in Section 4 we found that the errors are high close to the Sun, with values reaching ¿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^{\circ}$ and $l=180^{\circ}$, respectively. We report pecular motions associated with the Galactic spiral arm structure of up to 6-8km s-1, 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 lines-of-sight that result in higher errors in addition to those previously mentioned for specific distance ranges. We identify these regions as zones of avoidance for the application of the kinematic distance method. We also compare our results with the distance estimates of some well-studied molecular clouds and find in general good agreement (Section 5).

In summary, the extraction of the velocity perturbations in our simulation has allowed us to determine what impact the potential has on systematic errors in the kinematic distance estimate. We find that within the spiral arms of the Galaxy, the kinematic distance errors are low as the gas lies within the local potential minima. Consequently, the line-of-sight velocity of the gas is close to what is expected for the axisymmetric version of our potential. We expect clouds within the spiral arms of the Milky Way to have low systematic kinematic distance errors. Conversely, the interarm regions present the largest deviation in both the kinematic distance and the line-of-sight velocity for a given Galactic radius. This is caused by the gas being sped up or slowed down as it travels into or out of the spiral arm. Additionally, we discovered a power law relation between the systematic kinematic distance error and the difference between the line-of-sight 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/Gadget-4 codes, along with the scripts for constructing the potential from Section 2 and for running N-body simulations with these codes in the external potential, are included in the Agama repository⁷.

https://github.com/GalacticDynamics-Oxford/Agama

Acknowledgements.

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

References

• Ahumada et al. (2020)

  Ahumada, R. et al. 2020, \hrefhttp://dx.doi.org/10.3847/1538-4365/ab929eApJS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020ApJS..249….3A249, 3

• Anderson et al. (2012)

  Anderson, L. D., Bania, T. M., Balser, D. S., & Rood, R. T. 2012, \hrefhttp://dx.doi.org/10.1088/0004-637X/754/1/62ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2012ApJ…754…62A754, 62

• Antoja et al. (2018)

  Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, \hrefhttp://dx.doi.org/10.1038/s41586-018-0510-7Nature, \hrefhttps://ui.adsabs.harvard.edu/abs/2018Natur.561..360A561, 360

• Antoja et al. (2022)

  Antoja, T., Ramos, P., López-Guitart, F., et al. 2022, \hrefhttp://dx.doi.org/10.1051/0004-6361/202244064A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2022A&A…668A..61A668, A61

• Baba et al. (2009)

  Baba, J., Asaki, Y., Makino, J., et al. 2009, \hrefhttp://dx.doi.org/10.1088/0004-637X/706/1/471ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2009ApJ…706..471B706, 471

• Balser et al. (2015)

  Balser, D. S., Wenger, T. V., Anderson, L. D., & Bania, T. M. 2015, \hrefhttp://dx.doi.org/10.1088/0004-637X/806/2/199ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2015ApJ…806..199B806, 199

• Barnes et al. (2023)

  Barnes, A. T., Watkins, E. J., Meidt, S. E., et al. 2023, \hrefhttp://dx.doi.org/10.3847/2041-8213/aca7b9ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2023ApJ…944L..22B944, L22

• Binney & Tremaine (2008)

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

• Bland-Hawthorn & Gerhard (2016)

  Bland-Hawthorn, J. & Gerhard, O. 2016, \hrefhttp://dx.doi.org/10.1146/annurev-astro-081915-023441ARA&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2016ARA&A..54..529B54, 529

• Bovy & Rix (2013)

  Bovy, J. & Rix, H.-W. 2013, \hrefhttp://dx.doi.org/10.1088/0004-637X/779/2/115ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2013ApJ…779..115B779, 115

• Brand & Blitz (1993)

  Brand, J. & Blitz, L. 1993, A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/1993A&A…275…67B275, 67

• Burton & Gordon (1978)

  Burton, W. B. & Gordon, M. A. 1978, A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/1978A&A….63….7B63, 7

• Burton & Liszt (1993)

  Burton, W. B. & Liszt, H. S. 1993, A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/1993A&A…274..765B274, 765

• Carroll & Ostlie (2017)

  Carroll, B. W. & Ostlie, D. A. 2017, An introduction to modern astrophysics, Second Edition (Cambridge University Press)

• Cautun et al. (2020)

  Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, \hrefhttp://dx.doi.org/10.1093/mnras/staa1017MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020MNRAS.494.4291C494, 4291

• Chandra et al. (2024)

  Chandra, V., Naidu, R. P., Conroy, C., et al. 2024, \hrefhttps://ui.adsabs.harvard.edu/abs/2024arXiv240601676C\hrefhttp://dx.doi.org/10.48550/arXiv.2406.01676arXiv e-prints, arXiv:2406.01676

• Chatzopoulos et al. (2015)

  Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, \hrefhttp://dx.doi.org/10.1093/mnras/stu2452MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2015MNRAS.447..948C447, 948

• Chemin et al. (2015)

  Chemin, L., Renaud, F., & Soubiran, C. 2015, \hrefhttp://dx.doi.org/10.1051/0004-6361/201526040A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2015A&A…578A..14C578, A14

• Clarke & Gerhard (2022)

  Clarke, J. P. & Gerhard, O. 2022, \hrefhttp://dx.doi.org/10.1093/mnras/stac603MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2022MNRAS.512.2171C512, 2171

• Clemens (1985)

  Clemens, D. P. 1985, \hrefhttp://dx.doi.org/10.1086/163386ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/1985ApJ…295..422C295, 422

• Coleman et al. (2020)

  Coleman, B., Paterson, D., Gordon, C., Macias, O., & Ploeg, H. 2020, \hrefhttp://dx.doi.org/10.1093/mnras/staa1281MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020MNRAS.495.3350C495, 3350

• Contopoulos & Grosbol (1989)

  Contopoulos, G. & Grosbol, P. 1989, \hrefhttp://dx.doi.org/10.1007/BF00873080A&A Rev., \hrefhttps://ui.adsabs.harvard.edu/abs/1989A&ARv…1..261C1, 261

• Correa Magnus & Vasiliev (2022)

  Correa Magnus, L. & Vasiliev, E. 2022, \hrefhttp://dx.doi.org/10.1093/mnras/stab3726MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2022MNRAS.511.2610C511, 2610

• Dehnen (1993)

  Dehnen, W. 1993, \hrefhttp://dx.doi.org/10.1093/mnras/265.1.250MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/1993MNRAS.265..250D265, 250

• Draine (2011)

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

• Eilers et al. (2020)

  Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2020, \hrefhttp://dx.doi.org/10.3847/1538-4357/abac0bApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2020ApJ…900..186E900, 186

• Eilers et al. (2019)

  Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, \hrefhttp://dx.doi.org/10.3847/1538-4357/aaf648ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2019ApJ…871..120E871, 120

• Einasto (1969)

  Einasto, J. 1969, \hrefhttp://dx.doi.org/10.1002/asna.19682910303AN, \hrefhttps://ui.adsabs.harvard.edu/abs/1969AN….291…97E291, 97

• Elia et al. (2022)

  Elia, D., Molinari, S., Schisano, E., et al. 2022, \hrefhttp://dx.doi.org/10.3847/1538-4357/aca27dApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2022ApJ…941..162E941, 162

• Fich et al. (1989)

  Fich, M., Blitz, L., & Stark, A. A. 1989, \hrefhttp://dx.doi.org/10.1086/167591ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/1989ApJ…342..272F342, 272

• Foster & MacWilliams (2006)

  Foster, T. & MacWilliams, J. 2006, \hrefhttp://dx.doi.org/10.1086/500815ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2006ApJ…644..214F644, 214

• Freudenreich (1998)

  Freudenreich, H. T. 1998, \hrefhttp://dx.doi.org/10.1086/305065ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/1998ApJ…492..495F492, 495

• Gaia Collaboration et al. (2023)

  Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, \hrefhttp://dx.doi.org/10.1051/0004-6361/202243940A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2023A&A…674A…1G674, A1

• Glover et al. (2010)

  Glover, S. C. O., Federrath, C., Mac Low, M. M., & Klessen, R. S. 2010, \hrefhttp://dx.doi.org/10.1111/j.1365-2966.2009.15718.xMNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2010MNRAS.404….2G404, 2

• Gómez (2006)

  Gómez, G. C. 2006, \hrefhttp://dx.doi.org/10.1086/508412AJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2006AJ….132.2376G132, 2376

• GRAVITY Collaboration et al. (2022)

  GRAVITY Collaboration, Abuter, R., Aimar, N., et al. 2022, \hrefhttp://dx.doi.org/10.1051/0004-6361/202142465A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2022A&A…657L..12G657, L12

• GRAVITY Collaboration et al. (2019)

  GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2019, \hrefhttp://dx.doi.org/10.1051/0004-6361/201935656A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2019A&A…625L..10G625, L10

• GRAVITY Collaboration et al. (2021)

  GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2021, \hrefhttp://dx.doi.org/10.1051/0004-6361/202040208A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2021A&A…647A..59G647, A59

• Henshaw et al. (2023)

  Henshaw, J. D., Barnes, A. T., Battersby, C., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, \hrefhttps://ui.adsabs.harvard.edu/abs/2023ASPC..534…83H83

• Junqueira et al. (2013)

  Junqueira, T. C., Lépine, J. R. D., Braga, C. A. S., & Barros, D. A. 2013, \hrefhttp://dx.doi.org/10.1051/0004-6361/201219769A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2013A&A…550A..91J550, A91

• Jurić et al. (2008)

  Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, \hrefhttp://dx.doi.org/10.1086/523619ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2008ApJ…673..864J673, 864

• Khanna et al. (2019)

  Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, \hrefhttp://dx.doi.org/10.1093/mnras/stz2462MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2019MNRAS.489.4962K489, 4962

• Klessen & Glover (2016)

  Klessen, R. S. & Glover, S. C. O. 2016, in Saas-Fee Advanced Course, Vol. 43, Saas-Fee Advanced Course, ed. Y. Revaz, P. Jablonka, R. Teyssier, & L. Mayer, \hrefhttps://ui.adsabs.harvard.edu/abs/2016SAAS…43…85K85

• Koch & Rosolowsky (2015)

  Koch, E. W. & Rosolowsky, E. W. 2015, \hrefhttp://dx.doi.org/10.1093/mnras/stv1521MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2015MNRAS.452.3435K452, 3435

• Koposov et al. (2023)

  Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, \hrefhttp://dx.doi.org/10.1093/mnras/stad551MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2023MNRAS.521.4936K521, 4936

• Krumholz & McKee (2005)

  Krumholz, M. R. & McKee, C. F. 2005, \hrefhttp://dx.doi.org/10.1086/431734ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2005ApJ…630..250K630, 250

• Lallement et al. (2019)

  Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, \hrefhttp://dx.doi.org/10.1051/0004-6361/201834695A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2019A&A…625A.135L625, A135

• Laporte et al. (2018a)

  Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V., & Garavito-Camargo, N. 2018a, \hrefhttp://dx.doi.org/10.1093/mnras/stx2146MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2018MNRAS.473.1218L473, 1218

• Laporte et al. (2018b)

  Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018b, \hrefhttp://dx.doi.org/10.1093/mnras/sty1574MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2018MNRAS.481..286L481, 286

• Leike et al. (2020)

  Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, \hrefhttp://dx.doi.org/10.1051/0004-6361/202038169A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2020A&A…639A.138L639, A138

• Levine et al. (2006)

  Levine, E. S., Blitz, L., & Heiles, C. 2006, \hrefhttp://dx.doi.org/10.1126/science.1128455Science, \hrefhttps://ui.adsabs.harvard.edu/abs/2006Sci…312.1773L312, 1773

• Li et al. (2022)

  Li, Z., Shen, J., Gerhard, O., & Clarke, J. P. 2022, \hrefhttp://dx.doi.org/10.3847/1538-4357/ac3823ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2022ApJ…925…71L925, 71

• Lu et al. (2020)

  Lu, Z.-J., Pelkonen, V.-M., Padoan, P., et al. 2020, \hrefhttp://dx.doi.org/10.3847/1538-4357/abbd8fApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2020ApJ…904…58L904, 58

• Mac Low & Klessen (2004)

  Mac Low, M.-M. & Klessen, R. S. 2004, \hrefhttp://dx.doi.org/10.1103/RevModPhys.76.125RMP, \hrefhttps://ui.adsabs.harvard.edu/abs/2004RvMP…76..125M76, 125

• Majewski et al. (2017)

  Majewski, S. R. et al. 2017, \hrefhttp://dx.doi.org/10.3847/1538-3881/aa784dAJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2017AJ….154…94M154, 94

• McClure-Griffiths & Dickey (2007)

  McClure-Griffiths, N. M. & Dickey, J. M. 2007, \hrefhttp://dx.doi.org/10.1086/522297ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2007ApJ…671..427M671, 427

• McClure-Griffiths et al. (2004)

  McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., & Green, A. J. 2004, \hrefhttp://dx.doi.org/10.1086/422031ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2004ApJ…607L.127M607, L127

• McKee & Ostriker (2007)

  McKee, C. F. & Ostriker, E. C. 2007, \hrefhttp://dx.doi.org/10.1146/annurev.astro.45.051806.110602ARA&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2007ARA&A..45..565M45, 565

• McMillan (2017)

  McMillan, P. J. 2017, \hrefhttp://dx.doi.org/10.1093/mnras/stw2759MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2017MNRAS.465…76M465, 76

• Messineo et al. (2004)

  Messineo, M., Habing, H. J., Menten, K. M., Omont, A., & Sjouwerman, L. O. 2004, \hrefhttp://dx.doi.org/10.1051/0004-6361:20034488A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2004A&A…418..103M418, 103

• Messineo et al. (2005)

  Messineo, M., Habing, H. J., Menten, K. M., et al. 2005, \hrefhttp://dx.doi.org/10.1051/0004-6361:20040533A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2005A&A…435..575M435, 575

• Messineo et al. (2002)

  Messineo, M., Habing, H. J., Sjouwerman, L. O., Omont, A., & Menten, K. M. 2002, \hrefhttp://dx.doi.org/10.1051/0004-6361:20021017A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2002A&A…393..115M393, 115

• Molinari et al. (2014)

  Molinari, S., Bally, J., Glover, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, \hrefhttps://ui.adsabs.harvard.edu/abs/2014prpl.conf..125M125–148

• Mróz et al. (2019)

  Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, \hrefhttp://dx.doi.org/10.3847/2041-8213/aaf73fApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2019ApJ…870L..10M870, L10

• Nakanishi & Sofue (2003)

  Nakanishi, H. & Sofue, Y. 2003, \hrefhttp://dx.doi.org/10.1093/pasj/55.1.191PASJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2003PASJ…55..191N55, 191

• Nakanishi & Sofue (2006)

  Nakanishi, H. & Sofue, Y. 2006, \hrefhttp://dx.doi.org/10.1093/pasj/58.5.847PASJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2006PASJ…58..847N58, 847

• Oh et al. (2010)

  Oh, C. S., Kobayashi, H., Honma, M., et al. 2010, \hrefhttp://dx.doi.org/10.1093/pasj/62.1.101PASJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2010PASJ…62..101O62, 101

• Oort et al. (1958)

  Oort, J. H., Kerr, F. J., & Westerhout, G. 1958, \hrefhttp://dx.doi.org/10.1093/mnras/118.4.379MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/1958MNRAS.118..379O118, 379

• Persic et al. (1996)

  Persic, M., Salucci, P., & Stel, F. 1996, \hrefhttp://dx.doi.org/10.1093/mnras/278.1.27MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/1996MNRAS.281…27P281, 27

• Plummer (1911)

  Plummer, H. C. 1911, \hrefhttp://dx.doi.org/10.1093/mnras/71.5.460MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/1911MNRAS..71..460P71, 460

• Pohl et al. (2008)

  Pohl, M., Englmaier, P., & Bissantz, N. 2008, \hrefhttp://dx.doi.org/10.1086/529004ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2008ApJ…677..283P677, 283

• Portail et al. (2017)

  Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, \hrefhttp://dx.doi.org/10.1093/mnras/stw2819MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2017MNRAS.465.1621P465, 1621

• Ramón-Fox & Bonnell (2018)

  Ramón-Fox, F. G. & Bonnell, I. A. 2018, \hrefhttp://dx.doi.org/10.1093/mnras/stx2866MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2018MNRAS.474.2028R474, 2028

• Reid (2022)

  Reid, M. J. 2022, \hrefhttp://dx.doi.org/10.3847/1538-3881/ac80bbAJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2022AJ….164..133R164, 133

• Reid & Brunthaler (2020)

  Reid, M. J. & Brunthaler, A. 2020, \hrefhttp://dx.doi.org/10.3847/1538-4357/ab76cdApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2020ApJ…892…39R892, 39

• Reid et al. (2016)

  Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, \hrefhttp://dx.doi.org/10.3847/0004-637X/823/2/77ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2016ApJ…823…77R823, 77

• Reid et al. (2019)

  Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, \hrefhttp://dx.doi.org/10.3847/1538-4357/ab4a11ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2019ApJ…885..131R885, 131

• Reid et al. (2014)

  Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, \hrefhttp://dx.doi.org/10.1088/0004-637X/783/2/130ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2014ApJ…783..130R783, 130

• Sanders et al. (2019)

  Sanders, J. L., Smith, L., & Evans, N. W. 2019, \hrefhttp://dx.doi.org/10.1093/mnras/stz1827MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2019MNRAS.488.4552S488, 4552

• Sofue (2011)

  Sofue, Y. 2011, \hrefhttp://dx.doi.org/10.1093/pasj/63.4.813PASJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2011PASJ…63..813S63, 813

• Sofue (2023)

  Sofue, Y. 2023, \hrefhttp://dx.doi.org/10.1007/s10509-023-04231-0Ap&SS, \hrefhttps://ui.adsabs.harvard.edu/abs/2023Ap&SS.368…74S368, 74

• Soler et al. (2022)

  Soler, J. D., Miville-Deschênes, M. A., Molinari, S., et al. 2022, \hrefhttp://dx.doi.org/10.1051/0004-6361/202243334A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2022A&A…662A..96S662, A96

• Sormani et al. (2015a)

  Sormani, M. C., Binney, J., & Magorrian, J. 2015a, \hrefhttp://dx.doi.org/10.1093/mnras/stv441MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2015MNRAS.449.2421S449, 2421

• Sormani et al. (2015b)

  Sormani, M. C., Binney, J., & Magorrian, J. 2015b, \hrefhttp://dx.doi.org/10.1093/mnras/stv2067MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2015MNRAS.454.1818S454, 1818

• Sormani et al. (2022a)

  Sormani, M. C., Gerhard, O., Portail, M., Vasiliev, E., & Clarke, J. 2022a, \hrefhttp://dx.doi.org/10.1093/mnrasl/slac046MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2022MNRAS.514L…1S514, L1

• Sormani et al. (2020)

  Sormani, M. C., Magorrian, J., Nogueras-Lara, F., et al. 2020, \hrefhttp://dx.doi.org/10.1093/mnras/staa2785MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020MNRAS.499….7S499, 7

• Sormani et al. (2022b)

  Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022b, \hrefhttp://dx.doi.org/10.1093/mnras/stac639MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2022MNRAS.512.1857S512, 1857

• Sormani et al. (2024)

  Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, \hrefhttp://dx.doi.org/10.1093/mnras/stae082MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2024MNRAS.528.5742S528, 5742

• Springel (2010)

  Springel, V. 2010, \hrefhttp://dx.doi.org/10.1111/j.1365-2966.2009.15715.xMNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2010MNRAS.401..791S401, 791

• Springel et al. (2021)

  Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, \hrefhttp://dx.doi.org/10.1093/mnras/stab1855MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2021MNRAS.506.2871S506, 2871

• Tielens (2005)

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

• Tress et al. (2020a)

  Tress, R. G., Smith, R. J., Sormani, M. C., et al. 2020a, \hrefhttp://dx.doi.org/10.1093/mnras/stz3600MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020MNRAS.492.2973T492, 2973

• Tress et al. (2024)

  Tress, R. G., Sormani, M. C., Girichidis, P., et al. 2024, \hrefhttps://ui.adsabs.harvard.edu/abs/2024arXiv240313048T\hrefhttp://dx.doi.org/10.48550/arXiv.2403.13048A&A, in press, arXiv:2403.13048

• Tress et al. (2020b)

  Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020b, \hrefhttp://dx.doi.org/10.1093/mnras/staa3120MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2020MNRAS.499.4455T499, 4455

• Tress et al. (2021)

  Tress, R. G., Sormani, M. C., Smith, R. J., et al. 2021, \hrefhttp://dx.doi.org/10.1093/mnras/stab1683MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2021MNRAS.505.5438T505, 5438

• van de Hulst et al. (1954)

  van de Hulst, H. C., Muller, C. A., & Oort, J. H. 1954, Bull. Astron. Inst. Netherlands, \hrefhttps://ui.adsabs.harvard.edu/abs/1954BAN….12..117V12, 117

• Vasiliev (2018)

  Vasiliev, E. 2018, \hrefhttps://ui.adsabs.harvard.edu/abs/2018arXiv180208255V\hrefhttp://dx.doi.org/10.48550/arXiv.1802.08255arXiv e-prints, arXiv:1802.08255

• Vasiliev (2019)

  Vasiliev, E. 2019, \hrefhttp://dx.doi.org/10.1093/mnras/sty2672MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2019MNRAS.482.1525V482, 1525

• Vasiliev et al. (2021)

  Vasiliev, E., Belokurov, V., & Erkal, D. 2021, \hrefhttp://dx.doi.org/10.1093/mnras/staa3673MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2021MNRAS.501.2279V501, 2279

• Wada et al. (2002)

  Wada, K., Meurer, G., & Norman, C. A. 2002, \hrefhttp://dx.doi.org/10.1086/342151ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2002ApJ…577..197W577, 197

• Watkins et al. (2023)

  Watkins, E. J., Kreckel, K., Groves, B., et al. 2023, \hrefhttp://dx.doi.org/10.1051/0004-6361/202346075A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2023A&A…676A..67W676, A67

• Wegg & Gerhard (2013)

  Wegg, C. & Gerhard, O. 2013, \hrefhttp://dx.doi.org/10.1093/mnras/stt1376MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2013MNRAS.435.1874W435, 1874

• Wegg et al. (2015)

  Wegg, C., Gerhard, O., & Portail, M. 2015, \hrefhttp://dx.doi.org/10.1093/mnras/stv745MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2015MNRAS.450.4050W450, 4050

• Wenger et al. (2018)

  Wenger, T. V., Balser, D. S., Anderson, L. D., & Bania, T. M. 2018, \hrefhttp://dx.doi.org/10.3847/1538-4357/aaaec8ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2018ApJ…856…52W856, 52

• Widmark et al. (2022a)

  Widmark, A., Laporte, C. F. P., & Monari, G. 2022a, \hrefhttp://dx.doi.org/10.1051/0004-6361/202142819A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2022A&A…663A..15W663, A15

• Widmark et al. (2022b)

  Widmark, A., Widrow, L. M., & Naik, A. 2022b, \hrefhttp://dx.doi.org/10.1051/0004-6361/202244453A&A, \hrefhttps://ui.adsabs.harvard.edu/abs/2022A&A…668A..95W668, A95

• Xu et al. (2006)

  Xu, Y., Reid, M. J., Zheng, X. W., & Menten, K. M. 2006, \hrefhttp://dx.doi.org/10.1126/science.1120914Science, \hrefhttps://ui.adsabs.harvard.edu/abs/2006Sci…311…54X311, 54

• Yaaqib et al. (2024)

  Yaaqib, R., Petersen, M. S., & Peñarrubia, J. 2024, \hrefhttp://dx.doi.org/10.1093/mnras/stae1363MNRAS, \hrefhttps://ui.adsabs.harvard.edu/abs/2024MNRAS.531.3524Y531, 3524

• Zucker et al. (2021)

  Zucker, C., Goodman, A., Alves, J., et al. 2021, \hrefhttp://dx.doi.org/10.3847/1538-4357/ac1f96ApJ, \hrefhttps://ui.adsabs.harvard.edu/abs/2021ApJ…919…35Z919, 35

• Zucker et al. (2022)

  Zucker, C., Goodman, A. A., Alves, J., et al. 2022, \hrefhttp://dx.doi.org/10.1038/s41586-021-04286-5Nature, \hrefhttps://ui.adsabs.harvard.edu/abs/2022Natur.601..334Z601, 334

Appendix A Spiral arm strength

We perform additional simulations to determine the impact of our choice of the spiral arm strength factor α in Equation (15). Specifically, we investigate three different values, α=0.204, 0.408 and 0.612, and find that they lead to a peak density contrast of 10%, 20% and 30% at the solar circle at R₀=8.179 kpc, respectively. We caution that Eilers et al. (2020) suggest a density contrast of ∼ 10% for the vicinity of the Sun, and note that the fiducial case of α=0.36 introduced in Section 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 Section 4.4, see also Reid et al. 2019 or Baba et al. 2009). This is the main reason for our choice.

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

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

Appendix B Parameter study

We perform a simple parameter study with the potential to understand how varying certain parameters impacts the lv diagram. This involves multiple simulations with the changed parameters. The parameters we consider are number of spiral arms (n), pitch angle of the spiral arms (i), spiral arm pattern speed (Ω_spiral) and bar pattern speed (Ω_bar). The values considered are indicated in Table 6. 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 allows 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 Figures 19 & 20. 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 $\sim-50^{\circ}$ to $\sim-130^{\circ}$ that is present in the n=4 set but not in the n=2 set. At positive galactic longitudes, we find the expected additional spiral arms with the n=4 set. We find that two of the arms in this region of the lv diagram lie very close to each other for the n=4, appearing to have split from a spiral arm in the same region in lv space in the n=2 set.

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

The spiral arm pattern speed has a larger impact as the corresponding resonances for the spiral arms end up changing with pattern speed. This in turn causes the spiral arms to become more apparent with increasing pattern speed beyond the central region of the lv diagram. Additionally, the features associated with the spiral arms tend towards 0 km s-1 in lv space. Within the inner most ${}60^{\circ}$ of the Galaxy, the structure here also moves similarly to how it does with 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 towards R=0 kpc with increasing pattern speed (see Figure 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-1 in lv space, however the broadening of the features does not happen in this case.

It should be noted for all lv diagrams presented in this section have been selected to have approximately the same phase angle between the bar and spiral arms. This means that whilst each snapshot will has similar phase angle between bar and spiral arms, they will be at different stages in evolution. We select 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 modeling that is beyond the scope of this paper.

Table 6: Parameter values considered.

Figure 19: 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.

Figure 20: Same as Figure 19 but with the n=4 spiral arms subset.