# Signatures of simulated spiral arms on radial actions

P A Palicio, A. Recio-Blanco, T Tepper-Garcia, E Poggio, Sébastien Peirani, Y Dubois, P.J Mcmillan, J Bland-Hawthorn, Katarina Kraljic, M. Barbillon

► To cite this version:

P A Palicio, A. Recio-Blanco, T Tepper-Garcia, E Poggio, Sébastien Peirani, et al. Signatures of simulated spiral arms on radial actions. Astronomy & Astrophysics - A&A, 2025, 695, pp.A193. 10.1051/0004-6361/202453611. hal-04871278

HAL Id: hal-04871278

https://hal.science/hal-04871278v1

Submitted on 26 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Signatures of simulated spiral arms on radial actions

P. A. Palicio<fnref n="1" />,*, A. Recio-Blanco<fnref n="1" />, T. Tepper-García<fnref n="2" />,<fnref n="3" />, E. Poggio<fnref n="4" />, S. Peirani<fnref n="5" />,<fnref n="6" />,<fnref n="1" />,<fnref n="7" />, Y. Dubois<fnref n="7" />, P. J. McMillan<fnref n="8" />, J. Bland-Hawthorn<fnref n="2" />,<fnref n="3" />, K. Kraljic<fnref n="9" />, and M. Barbillon<fnref n="1" />

Received 25 December 2024 / Accepted 17 February 2025

## ABSTRACT

*Context.* Spiral arms play a key role in the evolution of disc galaxies, defining their morphology, star formation, chemistry, and dynamics. Among their various implications, it has been observed in the Milky Way disc that the distribution of Gaia Data Release 3 (DR3) radial actions exhibits structures that might be related to the spiral arms.

*Aims.* Our goal is to investigate the relationship between regions of low radial action identified in simulated discs and the location of the spiral arms, such as that suggested for the Galaxy in previous studies.

*Methods.* For a sample of 23 simulated spiral galaxies, we modelled the axisymmetric component of their gravitational potential to compute the radial action of their stellar particles using the Stäckel fudge. The spatial distribution of the radial action was then compared to the location of the spiral arms, identified as overdensities in the stellar surface density using a kernel density estimator.

*Results.* Our analysis reveals a strong correlation between the radial action distribution and the spiral arms in 18 of 23 simulated galaxies. However, notable discrepancies are observed in the remaining five, since they are profoundly out-of-equilibrium systems, such as galaxies influenced by external interactions or spiral arms still in the process of winding up. Additionally, spiral arms are consistently traced across young, intermediate, and old stellar populations (≥3 Gyr) in most simulations, suggesting that they are supported by stars spanning a broad age range.

*Conclusions.* We have confirmed that, in general, there is a tendency for spatial correlation between spiral arms and stellar populations featuring low values of the radial action, as discussed in the literature using Gaia DR3 data. However, discrepancies between features in the radial action distribution and the spiral structure can be interpreted as signatures of recent disturbances, a scenario applicable to the Milky Way. Furthermore, populations at least as old as 3 Gyr trace the spiral arms with no significant misalignment across age bins, suggesting a possible theoretical interpretation of the observations obtained with Gaia data. A linear relation between the maximum value of the radial action of the spiral arms and the vertical scale-length is found, which is also satisfied by the Milky Way.

**Key words.** Galaxy: disk – Galaxy: evolution – Galaxy: kinematics and dynamics – Galaxy: structure – galaxies: kinematics and dynamics

## 1. Introduction

Spiral arms are ubiquitous structures in the Local Universe, present in about two-thirds of galaxies (Huertas-Company et al. 2011; Lintott et al. 2011; Willett et al. 2013; Kelvin et al. 2014; Hart et al. 2016, 2017; Masters et al. 2019). This prevalence has facilitated the morphological classification of spiral arms into three categories (Elmegreen & Elmegreen 1982): the ‘grand design’, typically featuring two large, continuous, well-defined arms that complete one revolution around the galaxy; flocculent spiral arms, characterised by a more irregular pattern composed of numerous fragmented segments with no apparent correlation among them (Goldreich & Lynden-Bell 1965; Binney & Tremaine 2008); and a third category comprising intermediate-scale arms, whose tracks exhibit bifurcations, branching, and are less extended than the ‘grand design’ arms.

Despite significant research over the last six decades, the formation mechanisms of spiral arms remain unclear, though several models have been proposed to explain these structures. These include the quasi-steady density wave theory (Lindblad 1963; Lin & Shu 1964; Bertin & Lin 1996; Shu 2016), the swing amplification model (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; Toomre 1981; Grand et al. 2012a,b, 2013; Dobbs & Baba 2014; Michikoshi & Kokubo 2016), transient spiral modes (Sellwood & Carlberg 2014), the dressed mass clump model (Toomre & Kalnajs 1991; D'Onghia et al. 2013), tidal interactions (Toomre & Toomre 1972; Tully 1974; Elmegreen & Elmegreen 1983; Oh et al. 2008; Dubinski et al. 2008; Dobbs et al. 2010), and recurrent cycles of groove modes (Sellwood & Carlberg 2019) among others.

As important drivers of secular evolution (Kormendy & Kennicutt 2004), spiral arms not only define the morphology of galaxies but also contribute to the redistribution of

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France</footnote>

<footnote n="2">² Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="3">³ Center of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO-3D), Australia</footnote>

<footnote n="4">⁴ INAF – Osservatorio Astrofisico di Torino, via Osservatorio 20, 10025 Pino Torinese (TO), Italy</footnote>

<footnote n="5">⁵ ILANCE, CNRS – University of Tokyo International Research Laboratory, Kashiwa, Chiba 277-8582, Japan</footnote>

<footnote n="6">⁶ Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="7">⁷ Sorbonne Université, CNRS, UMR7095, Institute d'Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France</footnote>

<footnote n="8">⁸ School of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote n="9">⁹ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote>* Corresponding author; pedro.alonso-palicio@oca.eu</footnote>
angular momentum (Lynden-Bell & Kalnajs 1972; Minchev & Famaey 2010; Minchev et al. 2011), the heating of the disc (Martinez-Medina et al. 2015; Grand et al. 2016), radial migration (Sellwood & Binney 2002; Di Matteo et al. 2013), and the thickening of the disc (Roškar et al. 2013), to cite a few examples. Therefore, the presence of spiral arms has implications for disc kinematics and dynamics. In contrast to face-on external galaxies, where the presence of the spiral structure can be detected by direct imaging, our position within the Galactic disc complicates the detection of the Milky Way's spiral arms by this technique, requiring more indirect methods such as those based on distributions of young objects (Reid et al. 2014, 2019; Xu et al. 2018, 2021a,b; Pantaleoni González et al. 2021; Zari et al. 2021; Gaia Collaboration 2023b,a; Ge et al. 2024; Rezaei et al. 2024), chemistry (Poggio et al. 2022; Spitoni et al. 2023; Hawkins 2023; Barbillon et al. 2025; Hackshaw et al. 2024; Debattista et al. 2025), kinematics (Antoja et al. 2011; Hunt et al. 2015; Martinez-Medina et al. 2022; Denyshchenko et al. 2024), and dynamics (Trick et al. 2017; Hunt et al. 2019; Sellwood et al. 2019; Palicio et al. 2023; Khalil et al. 2024). Sellwood et al. (2019) performed a detailed analysis of the Gaia Data Release 2 (DR2) dynamics (Gaia Collaboration 2018) in the very solar neighbourhood (~200 pc from the Sun) to test some of the spiral arm models mentioned above, concluding that the dressed mass model is favoured. Mata-Chávez et al. (2019) studied the relation between the physical parameters of simulated spiral arms with their formation mechanisms, although they did not perform a comparison with those of the Milky Way.

Using the full kinematics sample from Gaia DR3 (Gaia Collaboration 2023c), in Palicio et al. (2023) we mapped the radial actions, $J_R$, and identified areas in the Galactic disc characterised by low values of $J_R$ compared to the whole explored region. The location and shape of most of these areas were compatible with those believed to trace spiral arms identified using young tracers such as masers (Reid et al. 2014), upper main sequence stars (Poggio et al. 2021), and Cepheids (Lemasle et al. 2022). This correlation can be attributed to the fact that stars are believed to form on nearly circular orbits within spiral arms, resulting in low radial actions. Although satisfactory, we could not completely rule out the contribution of moving groups to the features seen in the $J_R$ maps, and some discrepancies were observed in the IV quadrant ($-90^{\circ} < \ell < 0^{\circ}$). Our results also suggested that older populations can support the structure of spiral arms, in alignment with the findings of Lin et al. (2022) and Uppal et al. (2023) concerning the distribution of red clump stars in the Local arm and the Outer arm, respectively.

In this work, we aim to confirm the relationship between spiral arms and regions of low $J_R$, originally reported in Palicio et al. (2023), through numerical simulations. We also address whether older populations can trace spiral arms. This paper is organised as follows. In Section 2, we present the set of simulations used in this work and describe the fitting procedure for the potential and acceleration field. Results and discussions are detailed in Section 3. Conclusions are summarised in Section 4.

## 2. Simulations

### 2.1. Selected sample and methodology

We analysed a sample of 23 disc galaxy simulations at redshift zero, comprising seven from the Auriga simulations (Grand et al. 2017, 2024), ten from the Illustris-TNG collaboration (Vogelsberger et al. 2014; Pillepich et al. 2019; Nelson et al. 2019a,b; Pillepich et al. 2024), one from the NewHorizon

collaboration (Dubois et al. 2021), and the g2.79e12 NIHAO-UHD simulation (Buck et al. 2020). These simulations represent Milky Way-like galaxies and were selected for the pronounced spiral arms shown in their pre-visualisation images. The galaxies provided by the Auriga, Illustris-TNG, and NIHAO-UHD teams belong to a group or cluster of galaxies and, consequently, their evolution is also affected by interactions with their neighbours, as evidenced by the visualisation of their snapshots. We noted, however, that certain relative isolation criteria were applied to consider them as Milky Way-like candidates, as reported in their respective works. In contrast, the NewHorizon volume was chosen so that most galaxies evolve in a field environment, as does the one selected in this work. Apart from these cosmological simulations, we have considered four controlled simulations of idealised galaxies, two of which were presented in Bland-Hawthorn & Tepper-García (2021) and Tepper-García et al. (2022), respectively, and two unpublished simulations created with the NEXUS framework (Tepper-García et al. 2024). These controlled simulations include the interaction with a point mass intended as a proxy for the Sagittarius dwarf galaxy (Sgr), and differ from one another in their gas content and how it is treated: NEXUS A (cf. Bland-Hawthorn & Tepper-García 2021) is a pure N-body simulation without gas, NEXUS B and NEXUS D (cf. Tepper-García et al. 2022) include an 'inert' gaseous disc component with a mass equal to 20% and 10% the total (stars and gas) disc mass, respectively. In these simulations, the gas is treated with a simple isothermal equation of state. Finally, NEXUS C adopts the same initial conditions as NEXUS B, but it includes a prescription for galaxy formation, whereby the gas is allowed to cool and to form stars, accounting for their energetic and mass feedback. We kindly refer to the above mentioned references for more details on each the simulation suites. The nominal resolutions for the baryonic particles, as reported for each family of simulations in their respective works, are: 369 pc (AURIGA), 150 pc (ILLUSTRIS), 34 pc (NEWHORIZON), 36 pc (NEXUS), and 265 pc (NIHAO-UHD). The number of stellar particles contained in each simulation, as well as other relevant details, are presented in Table 1.

### 2.2. Model of the potential

In order to model the axisymmetric component of the gravitational potential, we set the centre of mass of the main galaxy of each simulation at the origin, and rotated the whole system to align the total angular momentum with the vertical axis. For those galaxies that exhibit clumps of accretions in their disc and halo, we manually removed them as their self-gravity constitutes a significant fraction of the potential at that position. By excluding these objects, we avoided the inclusion of misleading local potential values in the fitting process, which are not representative of the global main galaxy's potential.

We considered three intervals in galactocentric distance, separated at $r = r_{min}$ and $r = r_{max}$, and performed individual fits of the potential in each of them. The primary fit was conducted in the intermediate regime $r_{min} < r < r_{max}$, where most of the features are expected. For this interval, we adopted the series proposed by Hernquist & Ostriker (1992) as fit of our potential $\Phi_{mid}$, where

$$ \Phi_{\text{mid}}(r, \theta) = \sum_{n=0}^{N} \sum_{\substack{l=0 \\ \text{even}}}^{L} A_{n,l} \Phi_{n,l}(r, \theta), \qquad (1) $$
<table><caption>**Table 1.** Parameters of the simulations considered in this work.</caption><thead><tr><th>SimID</th><th>N<sub>part</sub><br>(×10<sup>6</sup>)</th><th>r<sub>min</sub><br>(kpc)</th><th>a<br>(kpc)</th><th>r<sub>max</sub><br>(kpc)</th><th>N</th><th>L</th><th>h<br>(kpc)</th><th>H<br>(kpc)</th><th>L<sub>kpc</sub><br>(kpc)</th><th>Z<sub>lim</sub><br>(kpc)</th></tr></thead><tbody><tr><td>Auriga 2</td><td>3.78</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.5</td><td>30</td><td>2.5</td></tr><tr><td>Auriga 3</td><td>2.79</td><td>2.0</td><td>5.0</td><td>40</td><td>8</td><td>8</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Auriga 16</td><td>5.20</td><td>2.5</td><td>5.0</td><td>50</td><td>12</td><td>8</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Auriga 19</td><td>2.44</td><td>2.0</td><td>5.0</td><td>40</td><td>12</td><td>8</td><td>0.50</td><td>2.0</td><td>15</td><td>0.5</td></tr><tr><td>Auriga 23</td><td>2.90</td><td>2.5</td><td>5.0</td><td>45</td><td>10</td><td>10</td><td>0.50</td><td>2.5</td><td>20</td><td>1.0</td></tr><tr><td>Auriga 24</td><td>5.34</td><td>2.0</td><td>5.0</td><td>40</td><td>12</td><td>8</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Auriga 27</td><td>3.55</td><td>2.5</td><td>5.0</td><td>40</td><td>8</td><td>8</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Illustris 424288</td><td>2.15</td><td>1.5</td><td>5.0</td><td>40</td><td>16</td><td>12</td><td>0.75</td><td>4.0</td><td>20</td><td>3.5</td></tr><tr><td>Illustris 456326</td><td>2.47</td><td>2.0</td><td>5.0</td><td>40</td><td>16</td><td>12</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Illustris 464163</td><td>1.82</td><td>1.0</td><td>5.0</td><td>40</td><td>12</td><td>8</td><td>0.75</td><td>4.0</td><td>30</td><td>2.5</td></tr><tr><td>Illustris 479290</td><td>1.71</td><td>1.0</td><td>6.0</td><td>50</td><td>14</td><td>10</td><td>0.75</td><td>2.5</td><td>25</td><td>1.5</td></tr><tr><td>Illustris 509091</td><td>1.00</td><td>2.5</td><td>5.0</td><td>40</td><td>14</td><td>10</td><td>0.75</td><td>2.5</td><td>25</td><td>2.0</td></tr><tr><td>Illustris 532760</td><td>1.13</td><td>2.5</td><td>5.0</td><td>40</td><td>10</td><td>8</td><td>0.75</td><td>2.5</td><td>25</td><td>2.0</td></tr><tr><td>Illustris 537941</td><td>0.88</td><td>2.5</td><td>4.0</td><td>35</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>20</td><td>2.0</td></tr><tr><td>Illustris 549516</td><td>0.73</td><td>1.5</td><td>4.0</td><td>50</td><td>14</td><td>10</td><td>0.75</td><td>4.0</td><td>30</td><td>1.5</td></tr><tr><td>Illustris 552414</td><td>1.00</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>Illustris 555287</td><td>0.65</td><td>2.0</td><td>2.5</td><td>35</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>20</td><td>1.0</td></tr><tr><td>NewHorizon</td><td>8.98</td><td>0.0</td><td>3.0</td><td>40</td><td>14</td><td>14</td><td>0.50</td><td>3.0</td><td>25</td><td>1.0</td></tr><tr><td>NEXUS A</td><td>5.00</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>30</td><td>1.0</td></tr><tr><td>NEXUS B</td><td>5.00</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>30</td><td>1.0</td></tr><tr><td>NEXUS C</td><td>5.00</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>30</td><td>1.0</td></tr><tr><td>NEXUS D</td><td>5.00</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>2.0</td><td>25</td><td>1.0</td></tr><tr><td>NIHAO-UHD</td><td>8.58</td><td>2.0</td><td>5.0</td><td>40</td><td>10</td><td>10</td><td>0.50</td><td>3.0</td><td>25</td><td>1.5</td></tr></tbody></table>

**Notes.** Identifier (first column), number of stellar particles (second column), minimum radii for $\Phi_{\text{mid}}(r_{\text{min}})$, third column), scale parameter for the potential ($a$, fourth column), maximum radii for $\Phi_{\text{mid}}(r_{\text{max}})$, fifth column), number of summation terms in $\Phi_{\text{mid}}$ (sixth and seventh columns), local and global kernel bandwidths (eighth and ninth columns), $X-Y$ range shown in the plots (tenth column), and vertical distance limit (last column).

in which

$$ \Phi_{n,l}(r, \theta) = -\sqrt{2\ell + 1} \frac{(ar)^{\ell}}{(a+r)^{2\ell+1}} C_n^{(2\ell+3/2)}(\xi) P_{\ell,0}(\cos \theta), \quad (2) $$

where $\xi = (r-a)/(r+a)$, $P_{\ell,0}$ are the Legendre polynomials, and $C_n^{(\alpha)}$ denotes the Gegenbauer polynomials (Gegenbauer 1884; Abramowitz et al. 1988). The parameter $a$ is a characteristic radial scale length of the disc (McMillan 2017), with values in the range of a few kiloparsecs (~1–5 kpc), to improve the convergence of the series in Eq. (1). The restriction to even values of $\ell$ in Eq. (1) imposes a vertically symmetric fit between both hemispheres, while the value $m=0$ in the associated Legendre polynomial implies no azimuthal variations of $\Phi_{\text{mid}}$ (hence no dependence on the angle $\phi$). The coefficients $A_{n,\ell}$ were computed by a least-squared fit of the potential value reported by the simulations, performed over a subsample of randomly selected particles since considering the whole simulation in this procedure is computationally expensive.

Once $\Phi_{\text{mid}}$ is obtained, approximations in the inner ($r \le r_{\text{min}}$) and outer ($r \ge r_{\text{max}}$) region were linked to it. For the inner region, we found good fits for the potential assuming the polynomial form

$$ \Phi_{\text{in}}(r, \theta) = \Phi_{\text{mid}}(r_{\text{min}}, \theta) \cdot \left[ 1 + a_{00}u + a_{12}u \cos^2 \theta + a_{22}u^2 \cos^2 \theta + a_{32}u^3 \cos^2 \theta + a_{24}u^2 \cos^4 \theta \right], \quad (3) $$

where $u = 1 - r/r_{\text{min}}$ and the coefficients $a_{ij}$ are the free parameters set by the least-square minimisation. For the outer region,

we considered a power-law like decay of the potential with distance,

$$ \frac{\Phi_{\text{out}}(r, \theta)}{\Phi_{\text{mid}}(r_{\text{max}}, \theta)} = \left[ 1 + \frac{(r - r_{\text{max}})^2}{R_c^2} \left( 1 + \left( \frac{1 - q^2}{q^2} \right) \cos^2 \theta \right) \right]^{-\gamma/2}, \quad (4) $$

where the free parameters are the characteristic radii $R_c$, the flattening term $q$, and the index of the power law $\gamma$ such that $\Phi_{\text{out}} \sim r^{-\gamma}$ at large distances. Note that both $\Phi_{\text{in}}$ and $\Phi_{\text{out}}$ are continuous with $\Phi_{\text{mid}}$ at $r = r_{\text{min}}$ and $r = r_{\text{max}}$, respectively. The final model of the potential, $\Phi_{\text{fit}}$, is obtained by combining the fits across the three distance regimes.

The quality of the potential modelling is illustrated in Figure 1 for a galaxy from each collaboration group, while Figure A.1 shows the corresponding results for the remaining simulations considered in this work. The colour code indicates the most discrepant relative difference between the reported potential and $\Phi_{\text{fit}}$. With the exception of the galactic centre, most regions of the galaxies show good agreement between the model and the original potential, with values of $|\Phi_{\text{real}} - \Phi_{\text{fit}}| < 0.1|\Phi_{\text{real}}|$, and generally below ~2–3%.

## 2.3. Model of the accelerations

The NewHorizon simulation does not directly provide the potential of the stellar particles; instead, it provides their accelerations, $\mathbf{a}$, which relate to the potential via $\mathbf{a} = -\nabla\Phi$. If we assume the potential is axisymmetric, the azimuthal component of the acceleration, $a_\phi$, must be zero. This leaves two components of
<figure><img src="image_2.png" /><figcaption>Fig. 1. Relative residuals of the potential fits (in percentages) for one simulation from each group considered in this work. The cells show the most extreme relative discrepancy between the model and the nominal potential reported in the simulation within the range $|Z| < Z_{lim}$. The black cross in the lower-left corner indicates the scale of ±10 kpc. The corresponding maps for the other simulations can be found in Figure A.1.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Model of the accelerations of the NewHorizon simulation. Upper (lower) panel: discrepancy between the reported radial (vertical) acceleration and its fit as a function of the acceleration itself. The colour-bars indicate the percentage of stellar particles enclosed in each region. Dashed (dotted) lines represent relative errors of ±5% (±10%). All the accelerations are expressed in units of $a_{\odot} \approx 6895.6 \text{ km}^2 \text{ s}^{-2} \text{ kpc}^{-1}$.</figcaption></figure>

the acceleration to fit: the radial component, $a_R$, and the vertical component, $a_Z$. Since a single component of the acceleration does not uniquely determine the potential, we need to account for both $a_R$ and $a_Z$ to estimate the proper $\Phi(r, \theta)$. Therefore, we assumed a model of the potential given by Eq. (1) and tuned the $A_{nl}$ parameters to minimise

$$S^2 = \sum_{i=1}^{N_p} \left( a_{R,i} + \left. \frac{\partial \Phi_{\text{mid}}}{\partial R} \right|_{r_i, \theta_i} \right)^2 + \sum_{i=1}^{N_p} \left( a_{Z,i} + \left. \frac{\partial \Phi_{\text{mid}}}{\partial Z} \right|_{r_i, \theta_i} \right)^2, \quad (5)$$

where the summation is performed over $N_p$ simulation particles and $(R, Z) = (r \sin\theta, r \cos\theta)$ are the cylindrical coordinates. We noted, however, the convergence of Eq. (5) is quite slow, especially at the galactic plane (Kuijken & Dubinski 1995; Dehnen & Binney 1998), so we have added the correction term

$$\Delta\Phi_{\text{mid}}(R, Z) = \sum_{n=1}^{3} \frac{c_n}{R^n} + c_4 \log R + c_5 \exp\left(-\frac{R}{5 \text{ kpc}} - \frac{|Z|}{|c_6|}\right), \quad (6)$$

in which the free parameters $c_n$ improved the least-square minimisation of $S^2$. We also tested a correction similar to that used by Dehnen & Binney (1998) and implemented in GALPY (Bovy 2015) to accelerate the spherical harmonic decomposition; however, Eq. (6) led to lower $S^2$ values.

The resulting fit of the accelerations is illustrated in Fig. 2. As can be seen, the relative errors in the radial acceleration are mostly smaller than 5%, while for $a_Z$ these discrepancy can be larger than ±10%. We should note, however, this mismatch is

appreciated only at low vertical accelerations (e.g., close to the galactic plane), where the radial forces dominate the dynamics.

For the particular case of the NewHorizon simulation, since no improvement in the fit of $a$ was observed when $\Phi_{in}$ was considered (see Eq. (3)), we decided to retain $\Phi_{mid}$ to model the inner galaxy. On the contrary, a potential of the form $\Phi_{out}$ is used at distances $r > 40$ kpc.

## 2.4. Spiral arm overdensities

In order to identify the spiral arms, we applied a bivariate kernel density estimator (KDE) technique similar to that exploited by Poggio et al. (2021). Using this approach, at the location $(X, Y)$ we approximate the local surface density, $\Sigma(X, Y; h)$, as

$$\Sigma(X, Y; h) = \sum_{i} m_i K(x_i - X; h) K(y_i - Y; h), \quad (7)$$

where $K(x; h) \sim (1 - x^2/h^2)$ if $|x| < h$ is the Epanechnikov kernel function, which weights the contribution of a particle at distance $|x|$ to the density at origin, assuming a characteristic bandwidth $h$ (Epanechnikov 1969). For distances $|x| > h$, the Epanechnikov kernel is zero. The summation in Eq. (7) is performed over all the particles located at $(x_i, y_i)$ with mass $m_i$. Analogously, the global surface density, $\Sigma(X, Y; H)$, is estimated as its local counterpart by using a bandwidth $H > h$, so that the contrast density, $\delta_{\Sigma}$, is
<figure><img src="image_4.png" /><figcaption>Fig. 3. Maps of the density contrast (first column), its bidimensional Fourier approximation (second column), and mass-weighted median radial action (third column) for the AURIGA 2 simulation. Open circles in second and third panels denote the inferred tracks for the spiral arms. The corresponding maps for the rest of simulations considered in this work can be found in Figs. B.1–B.4.</figcaption></figure>

determined by their relative difference as

$$ \delta_{\Sigma}(X, Y) = \frac{\Sigma(X, Y; h) - \Sigma(X, Y; H)}{\Sigma(X, Y; H)}. \qquad (8) $$

The selected values for h and H for each simulation are summarised in Table 1. The parameter h is set to approximately match the width of the spiral arms, whereas H represents the typical range of large-scale density variations in the disc.

Figures 3, B.1–B.4 illustrate the maps of $\delta_{\Sigma}(X, Y)$ for each simulation in the region $|Z| < Z_{lim}$ (values in Table 1). In most simulations, the distribution of $\delta_{\Sigma}$ is influenced by short-scale distortions, particularly in the outer regions of the galactic disc, where the effects of Poissonian noise become more pronounced due to the low density in these areas. We reduced the contribution of such fluctuations by reconstructing $\delta_{\Sigma}(X, Y)$ from its truncated Fourier series approximation as

$$ \tilde{\delta}_{\Sigma}(X, Y) = \sum_{n=0}^{11} \sum_{m=0}^{11} C_{n,m} \exp\left(\frac{2\pi i}{L_{kpc}} [n \cdot X + m \cdot Y]\right), \qquad (9) $$

where $C_{n,m}$ are the coefficients of the Fourier expansion of $\delta_{\Sigma}$ and $L_{kpc}$ denotes the limits of the square maps shown in Figs. 3, B.1–B.4. Thus, by approximating $\delta_{\Sigma}$ by $\tilde{\delta}_{\Sigma}$ (second columns in Figs. 3, B.1–B.4), we are able to reduce the noise and enhance the larger scale structures in the contrast density maps.

By observing the maps of $\delta_{\Sigma}$ and $\tilde{\delta}_{\Sigma}$, we visually identified the tentative trace of the spiral arms (dotted lines in the Figures) for sake of comparison with the features in $J_R$. The criteria for considering these features as part of the spiral arms depend on: (i) their extent, where segments significantly larger than their width are readily accepted; (ii) its location, as features in the less dense outer regions require special attention due to the Poissonian noise in the $\delta_{\Sigma}$ maps or from artefacts caused by the periodic boundaries of $\tilde{\delta}_{\Sigma}(X, Y)$, a natural consequence of the Fourier reconstruction; and (iii) the degree of agreement between $\delta_{\Sigma}$ and $\tilde{\delta}_{\Sigma}$, where, if the previous criteria are not enough to accept a feature, we check whether its inclusion is consistent with the extension of other nearby structures. Since these criteria rely on visual inspection of the reported maps, some minor or ambiguous features may have been omitted. However, none of these omissions are expected to significantly affect the overall spiral pattern of the galaxy, which is dominated by the large, clearly discernible structures usually located at $R < L_{kpc}$.

## 2.5. Computation of the radial actions

For each simulation, we estimated the orbital parameters, actions and energy by using the potential $\Phi_{fit}$ (see Section 2.2) in the code used in Palicio et al. (2023) instead of the Galactic potential of McMillan (2017). Assuming an axisymmetric potential, this code implements the Stäckel-fudge approximation method (Famaey & Dejonghe 2003; Binney 2012; Sanders & Binney 2016; Mackereth & Bovy 2018; Bland-Hawthorn et al. 2019; Hadden 2024, but see also Appendix B of Palicio et al. 2023) to estimate the apocenter, pericenter, eccentricity, maximum galactic height and actions from the input position and kinematics.

The maps of the radial action, $J_R$, can be seen in the rightmost columns of Figs. 3, B.1–B.4. In contrast to Palicio et al. (2023), where the number of Gaia DR3 stars exceeded that of particles in our simulations, we had to consider an Epanechnikov kernel with $h = 0.5$ kpc to effectively increase the number of particles per cell<fnref n="1" /> in the maps of $J_R$. The contribution of each particle to the median $J_R$ within a cell is weighted by its stellar mass.

# 3. Results and discussions

## 3.1. Radial action of the spiral arms

The maps represented in the first columns of Figures 3, B.1–B.4 illustrate the contrast of stellar mass density, $\delta_{\Sigma}$, among the different regions of each galaxy. In these panels, we can discern several arc-shaped structures whose densities exceed those of their surroundings (red areas), which we identify as the loci of the spiral arms. They are oriented around the galactic centre, with some bifurcating (as those in AURIGA 2 and AURIGA 3, for example), having kicks in their pitch angles (e.g., AURIGA 16 and AURIGA 24) or showing distortions due to remnants of accreted satellites (e.g., ILLUSTRIS 532760). The contrast density maps filtered for short-scale features, $\tilde{\delta}_{\Sigma}$, are shown in the middle columns of Figs. 3, B.1–B.4, and are adopted as a visual reference for exploring the noisy regions, in which we traced, by visual inspection, those over-density structures that we attributed to the spiral arms (dotted curves in Figs. 3, B.1–B.4). Conversely, we excluded the inner galaxy from our analysis since the high mass of the bulge biases the $\delta_{\Sigma}$ estimator there, leading to unrealistic values of the density contrast (blue annular regions at origin in left and middle columns of Figs. 3, B.1–B.4). Similarly, we did not consider as representative of the main galaxy morphology

<footnote n="1">¹ For comparison, each cell in Figs. 3, B.1–B.4 has ~0.375 kpc width.</footnote>
the features in $\delta_{\Sigma}$ created by remnants of accreted satellites, as those seen, for example, in ILLUSTRIS 424288 at $(X,Y) \approx (-24.5, 9.0)$, $(5.1, -25.0)$, and $(21.0, 17.0)$ kpc.

The spatial distribution of the median radial action (right columns) also reveals arc-shape structures that trace the low $J_R$ regions. In contrast to those in $\delta_{\Sigma}$, whose width typically ranges between ~1.5 and ~2.5 kpc, the features in $J_R$ can be up to ~4–5 kpc width (e.g., AURIGA 16 and ILLUSTRIS 555287). By confronting the maps of $\delta_{\Sigma}$ and median $J_R$, we find the overdensities of stellar mass tend to be located at the low radial action areas (dotted curves in middle and right columns), although noticeable discrepancies are observed both in the external zones (e.g., AURIGA 24, ILLUSTRIS 549516, and NEXUS A) and the central regions (e.g., ILLUSTRIS 424288).

To evaluate the quality of correspondence between the high $\delta_{\Sigma}$ and the low radial action arcs, we examined the value of $J_R$ at the location of the nodes that indicate our estimation of the spiral arm tracks (i.e., the circles that form the dotted curves in Figs. 3, B.1–B.4). The separation between two consecutive nodes is 2 kpc. Subsequently, these nodes are classified into three categories based on the degree of agreement with $J_R$: category ‘G’ (good) for the cases in which the $\delta_{\Sigma}-J_R$ correlation is clear, ‘F’ (fair) if a correspondence exists between the features in $\delta_{\Sigma}$ and $J_R$ despite some mismatch, and ‘B’ (bad) if no agreement is observed. Specifically, the criteria for category ‘F’ include cases that would align after shifts of $\le 1$ kpc in their position, a natural extension of their arcs, or are located in a slightly elevated $J_R$ area (light bluish areas in the map of the median $J_R$). As an illustrative example of this intermediate category, in ILLUSTRIS 424288, six of the seven last nodes of the spiral arm that ends at $(X,Y) \approx (15.3, 0.5)$ kpc were classified as ‘F’ nodes. Another example can be found in the nodes between $(X, Y) \approx (5.0, 20.)$, and $(10.0, 23)$ kpc in ILLUSTRIS 549516, since their median $J_R$ is slightly above the border value $(0.02 L_{\odot})$ and they bridge the extremes of the arc feature they belong to.

Due to the difference in the width of the arcs seen in the $\delta_{\Sigma}$ (or $\delta_{\Sigma}$) and $J_R$ maps, as well as the Poissonian noise in the outer regions (usually at $R \gtrsim L_{kpc}$), a correlation diagram of $\delta_{\Sigma}$ vs. $J_R$ does not clearly show the connection between the overdensities and the zones of low radial action. However, we present these diagrams in Appendix C where, despite the large scatter, all the simulations show an anticorrelation between $\delta_{\Sigma}$ and $J_R$.

The resulting classification revealed the nodes tagged as ‘G’ constitute, on average, the 62.6% of the total nodes of a simulation. This percentage, however, ranges from 14.7% for the NEXUS A simulation to 88.1% for AURIGA 16, while 18 of the 21 galaxies have percentages for the ‘G’ sample $\ge 50\%^2$. If we consider the combination of categories ‘G+F’ as positive matches, the percentage of agreement rises, on average, to 80.2%; while all the simulations exceed the 70% threshold but two: NEXUS A (57%) and AURIGA 2 (67%). The percentage of agreement for each galaxy is indicated in Table 2.

## 3.2. Anomalous scenarios

Although most of the spiral arms analysed above show a correlation with the low $J_R$ regions, we have detected a few discrepancies in some cases, especially in the NEXUS family of simulations, as well as in other numerical galaxies not included above but described and discussed as follows. In the first panel of Fig. B.4, we can observe a misalignment between the spiral

<table><caption>**Table 2.** Statistics of the node classification based on the agreement between the $\delta_{\Sigma}$ overdensity and the low $J_R$ regions.</caption><thead><tr><th>SimID</th><th>N<sub>nodes</sub></th><th>G(%)</th><th>F(%)</th><th>B(%)</th></tr></thead><tbody><tr><td>Auriga 2</td><td>218</td><td>48.2</td><td>18.8</td><td>33.0</td></tr><tr><td>Auriga 3</td><td>156</td><td>66.0</td><td>17.3</td><td>16.7</td></tr><tr><td>Auriga 16</td><td>134</td><td>88.1</td><td>5.2</td><td>6.7</td></tr><tr><td>Auriga 19</td><td>81</td><td>70.4</td><td>17.3</td><td>12.3</td></tr><tr><td>Auriga 24</td><td>163</td><td>71.8</td><td>6.7</td><td>21.5</td></tr><tr><td>Auriga 27</td><td>111</td><td>72.1</td><td>14.4</td><td>13.5</td></tr><tr><td>Illustris 424288</td><td>131</td><td>75.6</td><td>11.5</td><td>13.0</td></tr><tr><td>Illustris 464163</td><td>234</td><td>80.3</td><td>6.8</td><td>12.8</td></tr><tr><td>Illustris 479290</td><td>162</td><td>55.6</td><td>15.4</td><td>29.0</td></tr><tr><td>Illustris 509091</td><td>155</td><td>86.5</td><td>11.0</td><td>2.6</td></tr><tr><td>Illustris 532760</td><td>156</td><td>75.0</td><td>7.1</td><td>17.9</td></tr><tr><td>Illustris 537941</td><td>115</td><td>67.0</td><td>9.6</td><td>23.5</td></tr><tr><td>Illustris 549516</td><td>188</td><td>60.6</td><td>11.2</td><td>28.2</td></tr><tr><td>Illustris 552414</td><td>142</td><td>69.0</td><td>19.7</td><td>11.3</td></tr><tr><td>Illustris 555287</td><td>119</td><td>73.1</td><td>16.8</td><td>10.1</td></tr><tr><td>NewHorizon</td><td>162</td><td>72.8</td><td>8.0</td><td>19.1</td></tr><tr><td>NEXUS A</td><td>300</td><td>14.7</td><td>42.7</td><td>42.7</td></tr><tr><td>NEXUS B</td><td>282</td><td>34.4</td><td>41.5</td><td>24.1</td></tr><tr><td>NEXUS C</td><td>273</td><td>25.3</td><td>46.5</td><td>28.2</td></tr><tr><td>NEXUS D</td><td>268</td><td>49.3</td><td>21.6</td><td>29.1</td></tr><tr><td>NIHAO-UHD</td><td>97</td><td>59.8</td><td>19.6</td><td>20.6</td></tr></tbody></table>

**Notes.** Identifier (first column), number of nodes on the track of the spiral arm (second column), percentage of nodes classified in the G, F, and B categories (third, fourth, and fifth columns, respectively).

arms, as traced by the overdensity in $\delta_{\Sigma}$, and the low $J_R$ regions for the NEXUS A simulation. This mismatch is more evident at the outermost parts ($R \gtrsim 20$ kpc) of the main spiral arms of NEXUS A, NEXUS C, and NEXUS D (first, third, and fourth panels in Fig. B.4, respectively); as well as in the innermost ~10 kpc where the spirals wind up bringing closer the features in $J_R$. We interpreted these mismatches as a consequence of a configuration far from the equilibrium for these simulations, as evidenced by the void in the inter-arm region at $Y \gtrsim 20$ kpc. Furthermore, we verified the radial velocity field, $V_R$, suggests an eventual wind up of the major arm there, removing the gap and leading to a more relaxed configuration. The same argument applies in the $R \le 10$ kpc regions, where the particles in the arm (inter-arm) show inwards (outwards) radial motions. In other words, the spatial shifts observed between the high density and the low $J_R$ areas are compatible with a spiral pattern that is still winding up.

Other major discrepancy concerns the spiral arm detected at the lowermost boundary of the above mentioned gap; such as that starting at $(X, Y) \approx (28, 0)$ kpc in NEXUS A but with analogous examples in all the NEXUS simulations. This spiral arms do not lie on a low radial action area but trace the envelope of the highest $J_R$ region. Although the out-of-equilibrium scenario described above can explain such discrepancy with the other simulations, there is an alternative explanation for such mismatch concerning the kernel technique described in Sect. 2.4: the abrupt transition from moderately populated (i.e., region below the dotted curve starting at $X=20$ kpc, $Y=0$ kpc) to almost empty areas (above the dotted curve) can produce misleading estimations of the density contrasts, which might be misinterpreted as spiral arms. This idea is explained in detail in Appendix E.

<footnote n="2">2 This includes the NEXUS D simulation, whose percentage for the ‘G’ subset is 49.3%.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Maps of the density contrast (upper panels) and mass-weighted median radial action (bottom panels) for the AURIGA 23 barred simulation analysed in Section 3.2. Open circles denote the inferred tracks for the spiral arms.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Absolute value of the dipolar-to-axisymmetric component of the density ratio, $|A_2/A_0|$, for all the simulations considered in this work. Large values of $A_2/A_0$ denote great contribution from elongated structures, like a bar, to the local density. For the sake of visualisation, the complete set of simulations has been divided into five panels.</figcaption></figure>

In addition to the AURIGA simulations presented in Figs. 3 and B.1, we analysed the peculiar case of AURIGA 23, whose overdensity map (upper-left panel in Fig. 4) suggests the presence of a long bar in the central regions. As commonly done in the literature (e.g., Martinez-Valpuesta et al. 2017; Thomas et al. 2023; Bland-Hawthorn et al. 2024), we used the ratio between the amplitudes of the dipolar $A_2$ and the axisymmetric $A_0$ components of the surface density to confirm the existence of the bar and quantify its strength. We computed $A_2/A_0$ over the range $R = 0.5$ to 7.5 kpc using bins of $\Delta R = 0.5$ kpc width with a 50% overlap.

As shown in Fig. 5, AURIGA 23 exhibits a large $A_2/A_0$ ratio (~0.4) compared to most of the simulations, surpassed only by AURIGA 24 ($A_2/A_0 \approx 0.54$), two earlier snapshots of AURIGA 23 from ~3 Gyr and ~1 Gyr ago ($A_2/A_0 \approx 0.53$ and ≈0.41, respectively), and ILLUSTRIS 552414 ($A_2/A_0 \approx 0.52$) with a ~2 kpc half-length bar. Despite being stronger, the bars

in AURIGA 24 and ILLUSTRIS 552414 do not seem to cause significant distortions in the distribution of radial actions. To explore this in detail, we located the resonances associated with the $m = 2$ and $m = 4$ perturbers of these galaxies and checked whether the radial action maps show any imprint of them at these positions. In particular, when two perturbers with different pattern speeds coexist, an overlap of their resonances can induce a significant exchange of angular momentum (Quillen 2003; Minchev & Famaey 2010; Minchev et al. 2011) in a non-linear fashion (i.e., the effect of the combined perturbers exceeds that produced by each of them individually) – especially if the bar’s corotation or outer Lindblad resonance is involved. This change in angular momentum would imply a variation in the radial action whenever a resonance other than corotation participates (Sellwood & Binney 2002), as suggested by the increase in the radial velocity dispersion reported by Minchev & Famaey (2010).
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Location of the resonances created by the *m* = 2 and the *m* = 4 perturbers in the AURIGA 24 and ILLUSTRIS 552414 simulations. Blue (red) vertical lines refer to the *m* = 2 (*m* = 4) perturber; while the linestyle denotes the frequency ratio involved: corotation (solid lines), first harmonics 2 : *l*<sub>b</sub> and 4 : *l*<sub>s</sub> (dashed lines), second harmonics 1 : *l*<sub>b</sub> and 2 : *l*<sub>s</sub> (dotted lines).</figcaption></figure>

Therefore, for the AURIGA 24 and ILLUSTRIS 552414 galaxies, we located the resonances of their *m* = 2 and *m* = 4 perturbers, estimating their pattern speeds through the 2D Fourier mode velocity method (Pfenniger et al. 2023). We noted that, due to the complexity of the structures in the discs of simulated galaxies—where multiple perturbers can coexist with different pattern speeds and transient behaviour—the inferred rotation Ω and epicyclic κ frequencies may not decrease monotonically with radius. Consequently, a galaxy may exhibit more than one outer Lindblad resonance, as is the case for AURIGA 24 and ILLUSTRIS 552414.

As Figure 6 illustrates, the AURIGA 24 simulation shows three resonance overlaps created by the *m* = 2 and *m* = 4 modes (hereafter denoted by the subscripts ‘b’ and ‘s’, respectively): an overlap between the corotation resonances at *R*<sub>CR,b</sub> ≈ *R*<sub>CR,s</sub> ≈ 3.0 kpc, and two overlaps involving the 2 : *l*<sub>b</sub>, 2 : *l*<sub>s</sub>, and 4 : *l*<sub>s</sub> resonances at ~5.75 kpc and 16.5 kpc. Among them, only the corotation overlap appears to have a significant impact on the radial action distribution, acting as a barrier for the innermost low-*J*<sub>R</sub> regions.

Conversely, in the ILLUSTRIS 552414 simulation, the coro-
tation radii *R*<sub>CR,b</sub> and *R*<sub>CR,s</sub> are slightly more separated, while the
2 : *l*<sub>b</sub> resonances overlap with the ultraharmonic 2 : *l*<sub>s</sub> and the
outer Lindblad resonances 4 : *l*<sub>s</sub> at 5.5, 8.5, 10.5, and 16.5 kpc.
However, we detected no distinct features in the *J*<sub>R</sub> distribu-
tion that suggest any influence from these overlaps, except for
limiting the innermost low-*J*<sub>R</sub> regions at the corotation radii,
similar to AURIGA 24. In contrast to AURIGA 24 and ILLUSTRIS
552414, the AURIGA 23 snapshots reveal two crescent-shaped
regions of low *J*<sub>R</sub> located on either side of the bar’s minor axis,
indicating a notable impact on the radial action distribution in
this area.

Finally, the analysis of the ILLUSTRIS 456326 simulation
revealed a significant mismatch between the overdensities in δΣ
and the features in the median J<sub>R</sub> (Fig. 7). The latter exhib-
ited very extended structures compared to those observed in
the other simulations. To evaluate the reasons behind these dis-
crepancies, we explored the recent history of this simulation
using earlier snapshots. 1.35 Gyr ago, the spiral arms of ILLUS-
TRIS 456326t91 showed a strong correlation with the low-J<sub>R</sub>
regions, representing one of the clearest examples in our study.
This agreement persisted until snapshot ILLUSTRIS 456326t94,
corresponding to 0.9 Gyr ago, where a ring-shaped area of low
radial action emerged at R ≈ 6.5 kpc, and slight misalignments
between the spiral arms and the arcs in J<sub>R</sub> became noticeable.

At later times, the maps from subsequent snapshots reveal
a new configuration of the overdensities in δΣ, which reorgan-
ised to eventually form two large spiral arms. Meanwhile, the

ring-shaped low-*J*<sub>R</sub> area expanded to *R* ≈ 10.8 kpc in snapshot ILLUSTRIS 456326t95 (0.7 Gyr ago). In the following snapshots, there is a complete misalignment between the overdensities and the arcs in *J*<sub>R</sub>, with no discernible patterns emerging. At the present time, the ILLUSTRIS 456326 simulation shows two banana-shaped regions characterised by low radial action, of which only the uppermost region includes a spiral arm.

We interpreted the evolution described above as the result
of an interaction between the main ILLUSTRIS 456326 galaxy
and a nearby satellite, whose mass comprises 20% of the total
system's mass. The pericenter of this encounter was located, at
most, at 44.2 kpc (snapshot ILLUSTRIS 456326t95) from the
centre of the main body, producing an overlap between the discs
of both galaxies approximately ~0.7 Gyr ago. The effects of this
interaction can be observed in Fig. A.2, where the companion
caused significant distortions and asymmetries in the potential
of the main galaxy.

3.3. Age of the spiral arm population

Apart from the low *J*<sub>R</sub>-overdensity correspondence, the results of Palicio et al. (2023) suggested the spiral arms might be also supported by populations older than expected. In that work, we were able to discern the spiral arms in the distribution of *J*<sub>R</sub> using a photometrically selected sample of Gaia DR3 giant stars. Furthermore, by reproducing the KDE technique described in Sect. 2.4, we observed overdensities of giant stars at the locations of the spiral arms reported by Poggio et al. (2021) and traced by upper main sequence stars. The idea of old populations contributing to the spiral arm structure is not new but was originally noted by Zwicky (1955) and Schweizer (1976), although for ‘grand design’ spirals (Eskridge et al. 2002). However, the spiral arms of the Milky Way are thought to be intermediate-scale or flocculent types (Binney & Tremaine 2008; Quillen et al. 2018; Martinez-Medina et al. 2022), and they are not expected to be long-lived structures (Oort 1970; Elmegreen & Elmegreen 1984).

The methodology of Palicio et al. (2023), however, offers an
alternative explanation for the features seen in their *J*<sub>R</sub> maps with
giants stars: their uncertainties in the photometric selection crite-
ria (observational errors, incorrect estimation of the extinction,
etc.) may have misclassified some young sources as giants, in
which the former would be preferentially located at the spiral
arms regions they were born. On the contrary, the giant popula-
tion would show a more uniform spatial distribution since it is
older. Although the polluting young sample were smaller in vol-
ume compared to the old one, the KDE method would enhance
the spatial asymmetries of the former revealing the spiral arms.

In order to address whether the old population can trace the
spiral arms or not, we have applied the KDE technique on three
subsets of different ranges of age, τ: the young (τ ≤ 1 Gyr), inter-
mediate (1 ≤ τ < 2 Gyr) and old (2 ≤ τ < 3 Gyr) populations.
The resulting maps for δΣ are illustrated in Figs. 8, D.1–D.3.
In general, all the spiral arms traced by the whole simulation
(dotted lines) are observed at the three age intervals, although a
few minor discrepancies are found in ILLUSTRIS 464163, where
half of the spiral arm starting at (X, Y) ≈ (1.4, −10.4) kpc is not
traced by the oldest sample; in ILLUSTRIS 479290, where the
small overdensities at the upper right corner of the maps can be
only identified in the young subset; as well as in NEWHORI-
ZON, in which the most external spiral arms are not discernible
in populations older than 1 Gyr. In contrast, the spiral arms of
the NIHAO-UHD simulation are almost exclusively supported
by the young population.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Maps of the density contrast (upper panels) and mass-weighted median radial action (bottom panels) for different snapshots of the ILLUSTRIS 456326 simulation, which is affected by a recent merger, as analysed in Section 3.2. Open circles denote the inferred tracks of the spiral arms.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Density contrast maps of the AURIGA 2 simulation for three distinct age intervals: stellar particles younger than 1 Gyr (left column), those aged between 1 and 2 Gyr (middle column), and those between 2 and 3 Gyr (right column). Open circles indicate the spiral arm tracks inferred from Fig. 3. Cross in the bottom-right corner indicates the scale of ±5 kpc.</figcaption></figure>

Our results align with the findings of Lin et al. (2022) and Uppal et al. (2023), who demonstrated that it is possible to trace the Local and Outer spiral arms, respectively, using the red clump population, whose age distribution is expected to peak at ~1–4 Gyr (Girardi & Salaris 2001; Girardi 2016). Similarly, Khanna et al. (2024) identified structures in the distribution of red clump stars consistent with the spiral arms traced by Drimmel (2000), although none correspond to the Local Arm. In a recent study, Debattista et al. (2025) identified the spiral arms of their simulated galaxy by analysing the population older than 2 Gyr.

We do not observe, however, any significant misalignment in the spiral arm traces across the three age bins (Rand & Kulkarni 1989). This contrasts with the findings of Hao et al. (2021), who reported offsets using masers and young open clusters. A key distinction is that their oldest age interval extends up to 0.1 Gyr, whereas our youngest age bin is an order of magnitude greater. Consequently, it is likely that our maps lack the resolution necessary to disentangle the offsets between populations of such small age difference.

### 3.4. Relation to the Milky Way

The arc-shaped structures identified in the maps of radial actions are not delimited by a single value of $J_R$ common to all the simulations. Instead, each galaxy requires an individual reference

value, $J_{R,w}$, that separates the high and low radial action regions (see rightmost colour-bars in Figs. 3, B.1–B.4). Motivated by the results of Jia et al. (2023) and Jia et al. (2024), who identified a correlation between the vertical scale-length, $h_Z$, and the mean radial action along the disc of the form $J_R \sim (h_Z - b)^2$, we investigated a potential connection between $J_{R,w}$ and the vertical and radial ($h_R$) scale-lengths. These parameters are represented on the horizontal axis of Fig. 9.

The scale-lengths $h_R$ and $h_Z$ were computed assuming exponential fits for the histograms of stellar masses along the radial and vertical directions, respectively. Since the discs of the analysed simulations show non-negligible deviations from exponential profiles in the radial direction, we considered the average of two estimations for $h_R$: one imposing $R < L_{kpc}$ and another with $R < L_{kpc}/2$ as the maximum galactocentric distance for the histogram. In both cases, the criteria $|Z| < Z_{lim}$ and $0.2L_{kpc} \le R$ were applied to focus on the galactic plane and exclude the bulge regions, respectively. Similarly, for $h_Z$, we averaged the estimations within $|Z| < Z_{lim}$ and $|Z| < Z_{lim}/2$, imposing $0.2L_{kpc} \le R < L_{kpc}$. The discrepancy between both estimations is denoted by the (symmetric) error bars in Fig. 9.

As shown in the upper panel of Fig. 9, no clear trend $J_{R,w}(h_R)$ can be inferred from the distribution of data points, primarily due to the large errors associated with the estimation of $h_R$. Its vertical counterpart, however, shows a positive correlation between the radial action $J_{R,w}$ and $h_Z$. To confirm this correlation, we
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Dependence of $J_{R,w}$ with the structure and kinematic parameters of the simulations: radial scale-length (upper panel) and vertical scale-length (lower panel). Coloured markers refer to simulation families: AURIGA (blue circles), ILLUSTRIS (orange triangles), NEWHORIZON (red diamond), NEXUS (purple triangles), and NIHAO-UHD (green square). Black markers denote measurements for the Milky Way reported in the literature (see main text). Open markers denote outliers excluded from the linear fit (dashed line) by the $\sigma$-clipping algorithm.</figcaption></figure>

performed a linear fit excluding outliers through a $\sigma$-clipping procedure: in each iteration, only the data points whose absolute discrepancy lay within the 75th percentile participated in the next fit (filled markers in Fig. 9). The resulting parameters of the fit confirmed our visual interpretation of the $J_{R,w}$ vs. $h_Z$ correlation, showing a positive slope of $3.69 \times 10^{-2} L_\odot/\text{kpc}$. This correlation is analogous to that found by Jia et al. (2023) at different ranges of Galactocentric radii, although we adopted a simpler functional form for the fit. Jia et al. (2023) proposed a parabolic fit for $J_{R,w}(h_Z)$, as predicted by some distribution functions for the Milky Way under the epicycle approximation (see Problem 4.21 in Binney & Tremaine 2008), while we found that a linear relation provides an adequate fit to the data points. It is important to note, however, that the parameter $J_{R,w}$ does not correspond to the mean radial action as in Jia et al. (2023), and we did not subdivide the data into radial bins as they did. Consequently, the observed trends of $J_{R,w}$ with scale lengths or velocity dispersions might appear linear rather than quadratic.

For comparison with the Milky Way, we included estimations for $h_R$ and $h_Z$ in the diagrams of Fig. 9 (black markers<fnref n="3" />). We adopted the value $J_{R,w} = (9.4 \pm 0.6) \times 10^{-3} L_\odot$ for the Milky Way from Figure 1 of Palicio et al. (2023), while those for $h_R$ were extracted from the compendium presented in Bland-Hawthorn & Gerhard (2016). For the vertical scale length, we considered $h_Z \approx 0.3$ kpc, as widely reported in the literature for decades (Gilmore & Reid 1983; Jurić et al. 2008).

The large error bars associated with the radial scale-length estimates prevent reliable comparisons, in which the Milky Way's $h_R$ remains the shortest among them. Conversely, the position of the Milky Way in the $J_{R,w}$ vs. $h_Z$ plane aligns with expectations from the linear fit, as well as with those of the NEXUS galaxies.

As noted by Jia et al. (2023), the correlation between $h_Z$ and the radial action can be explained as a consequence of the disc's underlying heating mechanisms, which broaden the distribution of radial actions and increase the velocity dispersions. This implies greater kinetic energy in the vertical direction, allowing stars to reach larger galactic heights. Since radial actions are always positive, the broadening of the $J_R$ distribution is associated with a shift of the mean and median towards higher values. This effect can be visualised assuming an exponential distribution for $J_R$, as in Fig. 3 of Jia et al. (2023), where the mean and median are proportional to the scale factor. We attribute the lack of correlation between $J_{R,w}$ and $h_R$ to their evolution through different mechanisms: the radial action is a key parameter determining the orbital eccentricity, ecc, around the Galactic centre ($ecc \propto \sqrt{J_R$ under the epicyclic approximation, as shown in Eq. 3.261 of Binney & Tremaine 2008). Thus, increasing the radial action implies larger ecc, leading to an outward displacement of the apocenter<fnref n="4" />, $r_{apo}$. In this configuration, stars can visit outer disc regions, which might suggest, erroneously, that $h_R$ has increased due to stellar mass being pushed outward. However, larger eccentricities also imply shorter galactocentric distances for the pericenter, resulting in stellar mass transport towards central regions and a reduction of the radial scale-length. Consequently, changes in radial action induced by heating mechanisms are unlikely to significantly impact $h_R$, particularly within the epicyclic approximation regime. However, when this approximation does not hold, as in the case of very eccentric orbits, a star may spend more time at galactocentric distances greater than its guiding radius than within it. This implies, on average, a net outward displacement of mass. In contrast, changes in angular momentum, $L_Z$, directly affect the guiding radii of the orbit, a proxy for the mean galactocentric distance, potentially triggering radial migration that redistributes mass outward, thereby increasing $h_R$. In alignment with this interpretation, Frankel et al. (2020) reported changes in the angular momentum of the Milky Way thin disc that are approximately 10 times larger than those in radial action. In a recent study, Hamilton et al. (2024) explored several models of spiral arms to reproduce the observed heating-to-migration ratio of ~0.1. They concluded that the non-linear horseshoe transport mechanism (Sellwood & Binney 2002) is compatible with the observations, albeit under rather specific conditions related to the spatial extent of the perturbation amplitude or the presence of greater pitch angles in the past. However, in the former scenario, the non-horseshoe scattering mechanism can also be consistent with the data.

It is interesting to compare the radial action maps presented in this work to the one obtained for the Milky Way, represented in Fig. 1 of Palicio et al. (2023). In our Galaxy, however, the comparison between the $J_R$ maps and the spiral structure is not trivial, because there is no consensus about the geometry and number of spiral arms in the Galactic disc. Low-$J_R$ features from Palicio et al. (2023) tend to be co-located with the spiral arm geometry proposed by Reid et al. (2014), with the exception of the Scutum arm and a small part of the Perseus arm in the third Galactic quadrant (see second panel of Fig. 4 in Palicio et al. 2023). The general good agreement between

<footnote n="3">³ The Milky Way data points were not included in the iterative fitting procedure.</footnote>

<footnote n="4">⁴ The eccentricity is usually approximated as $(r_{apo} - r_{peri})/(r_{apo} + r_{peri})$.</footnote>
the spiral arm model from Reid et al. (2014) and the observed low-$J_R$ features would suggest a connection between the two, as seen in the majority of the simulated galactic discs analysed in this work.

However, a larger pitch angle for the Perseus arm has been found by other works based on young stars in Gaia (Zari et al. 2021; Poggio et al. 2021; Gaia Collaboration 2023a) and Cepheids (Drimmel et al. 2025), in agreement with the model by Levine et al. (2006) based on HI observations (see their spiral arm no. 2). Such pitch angle would also be consistent with the map obtained with the giant sample in Palicio et al. (2023). Assuming the Perseus arm geometry from these works, the low-$J_R$ feature in the outer disc mapped by would not coincide with a spiral arm (see first and last panels of Fig. 4 in Palicio et al. 2023). Therefore, in this scenario, the outer disc of the Milky Way would be more similar to the out-of-equilibrium simulated galaxies analysed in this work (Antoja et al. 2018; McMillan et al. 2022; Poggio et al. 2024); although similar discrepancies have also been observed in less disturbed simulations.

Based on all spiral arm geometries explored, the structures tentatively associated with the Local and Perseus spiral arms, respectively, merge in an area of low $J_R$, showing a certain misalignment with the literature – except for the 17th group of Cepheids reported in Lemasle et al. (2022). By analysing the maps for $\delta_\Sigma$ and $J_R$ from the simulations, we identified similar structures, where the region between two nearby spiral arms also presents low radial action. This is the case for the bifurcation at (X, Y) ≈ (-8, 0) kpc in AURIGA 19 (third row in Fig. B.1), at (10, -6.6) kpc in AURIGA 27 (fifth row in Fig. B.1), at (-5.3, -7.3) kpc in ILLUSTRIS 555287 (fourth row in Fig. B.3), and the gap between two spiral arms in NEXUS D at (10.6, 0) kpc (fourth row in Fig. B.4). While reaching a definitive conclusion about the Milky Way remains challenging, there is no doubt that mapping $J_R$ features in the Galactic disc provides a unique perspective on the dynamical history of the Galaxy, as demonstrated by the simulations analysed in this work.

## 4. Conclusions

Our analysis of a set of simulations from independent groups supports the correlation between the location of the spiral arms and the low radial action features, as proposed by Palicio et al. (2023) using Gaia eDR3 and DR3 data (Gaia Collaboration 2021, 2023c). In that work, the features in $J_R$ were mainly consistent with the models of the Sagittarius, Local, and Perseus spiral arms proposed by Reid et al. (2014) in the first and second quadrants (0° ≤ $\ell$ ≤ 180°), while at negative longitudes, no preferential model emerged from the comparison with the literature. For our sample of simulations, such correspondence was observed in 20 of the 23 spiral galaxies, with an average percentage of agreement ranging from the 64% to the 88.1%, depending on the strictness of the goodness criterion. Conversely, for the remaining two simulations, the mismatch between the $\delta_\Sigma$ and the $J_R$ maps can be explained by an out of equilibrium state, such as that triggered by an interaction with a companion galaxy (ILLUSTRIS 456326, NEXUS A-D), and by border effects of the KDE technique, as these observed in the NEXUS simulations and explained in Appendix E.

Although restricted to the particular cases of AURIGA 23 and NIHAO-UHD, we found no evidence of a great influence of the bar on the distribution of the radial actions at the locations of the spiral arms. We must exclude from this observation, however, the central regions in which the bar is confined, as two arc-shaped areas of low $J_R$ have emerged near the minor axis.

In any case, all the conclusions derived for the barred galaxy scenario should be interpreted with caveat since our methodology assumes the potential is predominantly axisymmetric. Therefore, it might result unrealistic if the bar contributes significantly to the total potential, requiring more advanced approaches for the action estimation, such as that proposed by Debattista et al. (2025).

On the other hand, strong interactions with satellites can destroy the existing spiral arm pattern, leading to a new configuration in which the emerging spiral arms do not lie on low $J_R$ areas, as observed in the recent history of the ILLUSTRIS 456326 simulation. This perturbed scenario is expected to be transitory and eventually evolve to a new equilibrium state in which the $J_R-\delta_\Sigma$ anticorrelation is recovered. Further re-simulations of post-interaction galaxies are required for confirming this hypothesis.

The maps of the density contrast, $\delta_\Sigma$, computed for different bins in age revealed the population older than 1 Gyr can trace the spiral arms. Furthermore, at the spatial resolution of the considered simulations, no significant misalignment among the spiral arm tracks has been observed when compared the three age bins. This contradicts the classical picture of the spiral arms traced by the young sources recently formed in, implying differences in the mean age between the leading and trailing sides of the arms. We recall, however, this contrast in age might not be discernible in our maps due to the spatial resolution. Addressing whether the old population was formed in the spiral arms they are observed in or have been trapped later is beyond the scope of this paper.

The comparison of the limiting value for the low $J_R$ regions, $J_{R,w}$, with the vertical scale-length of the disc reveals a positive correlation between the two, which can be approximated by the linear relation presented in Fig. 9. Given that correlation does not necessarily imply causation, we interpreted that relation as consequence of disc heating, which increases the velocity dispersion raising the eccentricity of the orbits, thereby increasing the radial actions and, consequently, $J_{R,w}$. On the contrary, no clear trend is observed between the radial scale-length, $h_R$, and $J_{R,w}$, suggesting that disc heating is not the primary mechanism driving its expansion, although some influence cannot be entirely ruled out.

**Acknowledgements.** We thank the anonymous referee for their valuable comments, which have helped improve the revision of this manuscript. P. A. Palicio acknowledges the financial support from the Centre national d'études spatiales (CNES). TTG acknowledges partial financial support from the Australian Research Council (ARC) through an Australian Laureate Fellowship awarded to JBH. We have used simulations from the Auriga Project public data release (Grand et al. 2024) available at https://www.mpa.mpa-garching.mpg.de/auriga/data. The NewHorizon simulation was made possible thanks to computational resources at CINES, TGCC, and KISTI, and local clusters at the Institut d'Astrophysique de Paris and the Yonsei University. This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci, KSC-2017-G2-0003 by KISTI, and as a 'Grand Challenge' project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. Although GALPY is not explicitly used in this work, P. A. Palicio uses its source code as reference and recognises the credit for the work of Bovy (2015).

## References

* Abramowitz, M., Stegun, I. A., & Romer, R. H. 1988, Am. J. Phys., 56, 958

* Antoja, T., Figueras, F., Romero-Gómez, M., et al. 2011, MNRAS, 418, 1423

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

* Barbillon, M., Recio-Blanco, A., Poggio, E., et al. 2025, A&A, 693, A3

* Bertin, G., & Lin, C. C. 1996, Spiral Structure in Galaxies: A Density Wave Theory (Cambridge, MA: MIT Press)

* Binney, J. 2012, MNRAS, 426, 1324

* Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn.
Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bland-Hawthorn, J., & Tepper-García, T. 2021, MNRAS, 504, 3168

Bland-Hawthorn, J., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 486, 1167

Bland-Hawthorn, J., Tepper-García, T., Agertz, O., & Federrath, C. 2024, ApJ, 968, 86

Bovy, J. 2015, ApJS, 216, 29

Buck, T., Obreja, A., Macciò, A. V., et al. 2020, MNRAS, 491, 3461

Debattista, V. P., Khachaturyants, T., Amarante, J. A. S., et al. 2025, MNRAS, 537, 1620

Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429

Denyshchenko, S. I., Fedorov, P. N., Akhmetov, V. S., Velichko, A. B., & Dmytrenko, A. M. 2024, MNRAS, 527, 1472

Di Matteo, P., Haywood, M., Combes, F., Semelin, B., & Snaith, O. N. 2013, A&A, 553, A102

Dobbs, C., & Baba, J. 2014, PASA, 31, e035

Dobbs, C. L., Theis, C., Pringle, J. E., & Bate, M. R. 2010, MNRAS, 403, 625

D'Onghia, E., Vogelsberger, M., & Hernquist, L. 2013, ApJ, 766, 34

Drimmel, R. 2000, A&A, 358, L13

Drimmel, R., Khanna, S., Poggio, E., & Skowron, D. M. 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202451100

Dubinski, J., Gauthier, J. R., Widrow, L., & Nickerson, S. 2008, in ASPC, 396, Formation and Evolution of Galaxy Disks, eds. J. G. Funes, & E. M. Corsini, 321

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Elmegreen, D. M., & Elmegreen, B. G. 1982, MNRAS, 201, 1021

Elmegreen, B. G., & Elmegreen, D. M. 1983, ApJ, 267, 31

Elmegreen, D. M., & Elmegreen, B. G. 1984, ApJS, 54, 127

Epanechnikov, V. A. 1969, Theory Probab. Appl., 14, 153

Eskridge, P. B., Frogel, J. A., Pogge, R. W., et al. 2002, ApJS, 143, 73

Famaey, B., & Dejonghe, H. 2003, MNRAS, 340, 752

Frankel, N., Sanders, J., Ting, Y.-S., & Rix, H.-W. 2020, ApJ, 896, 15

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gaia Collaboration (Drimmel, R., et al.) 2023a, A&A, 674, A37

Gaia Collaboration (Recio-Blanco, A., et al.) 2023b, A&A, 674, A38

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

Ge, Q. A., Li, J. J., Hao, C. J., et al. 2024, AJ, 168, 25

Gegenbauer, L. 1884, Denkschr. Kaiserl. Akad. Wiss. Wien, Math.-Naturwiss. CL, 48, 293

Gilmore, G., & Reid, N. 1983, MNRAS, 202, 1025

Girardi, L. 2016, ARA&A, 54, 95

Girardi, L., & Salaris, M. 2001, MNRAS, 323, 109

Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125

Grand, R. J. J., Kawata, D., & Cropper, M. 2012a, MNRAS, 426, 167

Grand, R. J. J., Kawata, D., & Cropper, M. 2012b, MNRAS, 421, 1529

Grand, R. J. J., Kawata, D., & Cropper, M. 2013, A&A, 553, A77

Grand, R. J. J., Springel, V., Gómez, F. A., et al. 2016, MNRAS, 459, 199

Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179

Grand, R. J. J., Fragkoudi, F., Gómez, F. A., et al. 2024, MNRAS, 532, 1814

Hackshaw, Z., Hawkins, K., Filion, C., et al. 2024, ApJ, 977, 143

Hadden, S. 2024, ApJ, 972, 64

Hamilton, C., Modak, S., & Tremaine, S. 2024, ApJ, Submitted [arXiv:2411.08944]

Hao, C. J., Xu, Y., Hou, L. G., et al. 2021, A&A, 652, A102

Hart, R. E., Bamford, S. P., Willett, K. W., et al. 2016, MNRAS, 461, 3663

Hart, R. E., Bamford, S. P., Casteels, K. R. V., et al. 2017, MNRAS, 468, 1850

Hawkins, K. 2023, MNRAS, 525, 3318

Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375

Huertas-Company, M., Aguerri, J. A. L., Bernardi, M., Mei, S., & Sánchez Almeida, J. 2011, A&A, 525, A157

Hunt, J. A. S., Kawata, D., Grand, R. J. J., et al. 2015, MNRAS, 450, 2132

Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 490, 1026

Jia, Y., Chen, Y., Du, C., & Zhao, G. 2023, A&A, 669, A107

Jia, Y., Yang, C., Chen, Y., Du, C., & Zhao, G. 2024, A&A, 692, A167

Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810

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

Kelvin, L. S., Driver, S. P., Robotham, A. S. G., et al. 2014, MNRAS, 439, 1245

Khalil, Y. R., Famaey, B., Monari, G., et al. 2024, A&A, submitted [arXiv:2411.12800]

Khanna, S., Yu, J., Drimmel, R., et al. 2024, A&A, Submitted [arXiv:2410.22036]

Kormendy, J., & Kennicutt, Robert C., J. 2004, ARA&A, 42, 603

Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341

Lemasle, B., Lala, H. N., Kovtyukh, V., et al. 2022, A&A, 668, A40

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

Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646

Lin, Z., Xu, Y., Hou, L., et al. 2022, ApJ, 931, 72

Lindblad, B. 1963, Sto. An., 5, 5

Lintott, C., Schawinski, K., Bamford, S., et al. 2011, MNRAS, 410, 166

Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1

Mackereth, J. T., & Bovy, J. 2018, PASP, 130, 114501

Martinez-Medina, L. A., Pichardo, B., Pérez-Villegas, A., & Moreno, E. 2015, ApJ, 802, 109

Martinez-Medina, L., Pérez-Villegas, A., & Peimbert, A. 2022, MNRAS, 512, 1574

Martinez-Valpuesta, I., Aguerri, J. A. L., González-García, A. C., Dalla Vecchia, C., & Stringer, M. 2017, MNRAS, 464, 1502

Masters, K. L., Lintott, C. J., Hart, R. E., et al. 2019, MNRAS, 487, 1808

Mata-Chávez, M. D., Velázquez, H., Pichardo, B., et al. 2019, ApJ, 876, 6

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

McMillan, P. J., Petersson, J., Tepper-García, T., et al. 2022, MNRAS, 516, 4988

Michikoshi, S., & Kokubo, E. 2016, ApJ, 823, 121

Minchev, I., & Famaey, B. 2010, ApJ, 722, 112

Minchev, I., Famaey, B., Combes, F., et al. 2011, A&A, 527, A147

Nelson, D., Pillepich, A., Springel, V., et al. 2019a, MNRAS, 490, 3234

Nelson, D., Springel, V., Pillepich, A., et al. 2019b, Comput. Astrophys. Cosmol., 6, 2

Oh, S. H., Kim, W.-T., Lee, H. M., & Kim, J. 2008, ApJ, 683, 94

Oort, J. H. 1970, in IAU Symposium, 38, The Spiral Structure of our Galaxy, eds. W. Becker, & G. I. Kontopoulos, 1

Palicio, P. A., Recio-Blanco, A., Poggio, E., et al. 2023, A&A, 670, L7

Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968

Pfenniger, D., Saha, K., & Wu, Y.-T. 2023, A&A, 673, A36

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Pillepich, A., Sotillo-Ramos, D., Ramesh, R., et al. 2024, MNRAS, 535, 1721

Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104

Poggio, E., Recio-Blanco, A., Palicio, P. A., et al. 2022, A&A, 666, L4

Poggio, E., Khanna, S., Drimmel, R., et al. 2024, A&A, submitted [arXiv:2407.18659]

Quillen, A. C. 2003, AJ, 125, 785

Quillen, A. C., Carrillo, I., Anders, F., et al. 2018, MNRAS, 480, 3132

Rand, R. J., & Kulkarni, S. R. 1989, ApJ, 343, 760

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

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

Rezaei, S. K., Beuther, H., Benjamin, R. A., et al. 2021, A&A, 692, A255

Roškar, R., Debattista, V. P., & Loebman, S. R. 2013, MNRAS, 433, 976

Sanders, J. L., & Binney, J. 2016, MNRAS, 457, 2107

Schweizer, F. 1976, ApJS, 31, 313

Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785

Sellwood, J. A., & Carlberg, R. G. 2014, ApJ, 785, 137

Sellwood, J. A., & Carlberg, R. G. 2019, MNRAS, 489, 116

Sellwood, J. A., Trick, W. H., Carlberg, R. G., Coronado, J., & Rix, H.-W. 2019, MNRAS, 484, 3154

Shu, F. H. 2016, ARA&A, 54, 667

Spitoni, E., Cescutti, G., Recio-Blanco, A., et al. 2023, A&A, 680, A85

Tepper-García, T., Bland-Hawthorn, J., & Freeman, K. 2022, MNRAS, 515, 5951

Tepper-García, T., Bland-Hawthorn, J., Vasiliev, E., et al. 2024, MNRAS, 535, 187

Thomas, G. F., Famaey, B., Monari, G., et al. 2023, A&A, 678, A180

Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, eds. S. M. Fall, & D. Lynden-Bell, 111

Toomre, A., & Toomre, J. 1972, ApJ, 178, 623

Toomre, A., & Kalnajs, A. J. 1991, in Dynamics of Disc Galaxies, ed. B. Sundelius, 341

Trick, W. H., Bovy, J., D'Onghia, E., & Rix, H.-W. 2017, ApJ, 839, 61

Tully, R. B. 1974, ApJS, 27, 449

Uppal, N., Ganesh, S., & Schultheis, M. 2023, A&A, 673, A99

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Willett, K. W., Lintott, C. J., Bamford, S. P., et al. 2013, MNRAS, 435, 2835

Xu, Y., Bian, S. B., Reid, M. J., et al. 2018, A&A, 616, L15

Xu, Y., Bian, S. B., Reid, M. J., et al. 2021a, ApJS, 253, 1

Xu, Y., Hou, L. G., Bian, S. B., et al. 2021b, A&A, 645, L8

Zari, E., Rix, H. W., Frankel, N., et al. 2021, A&A, 650, A112

Zwicky, F. 1955, PASP, 67, 232
**Appendix A: Maps of the potential fit**

As mentioned in Section 2.2, for each simulation we performed a fit of the gravitational potential assuming an axisymmetric formu-
lation for it. The most deviant relative discrepancy between this fit and the potential reported in the simulations is illustrated in Figs.
A.1 and A.2.

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Relative residuals of the potential fits (in percentages) for each simulation considered in this work. Each cell shows the most extreme relative discrepancy between the model and the nominal potential reported in the simulation within $|Z| < Z_{lim}$. Black cross in the lower left corner indicates the scale of ±10 kpc.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Similar to Fig. A.1 but for the set of AURIGA 23 and ILLUSTRIS 456326 snapshots used in Section 3.2. Black cross in the lower left corner indicates the scale of ±10 kpc.</figcaption></figure>

**Appendix B: Additional maps of density contrast and radial actions**

In this section, we present the maps of the density contrasts, $\delta_{\Sigma}$ and $\tilde{\delta}_{\Sigma}$, as well as those of the radial actions for the simulations not shown in Figure 3.
<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Maps of the density contrast (first column), its bidimensional Fourier approximation (second column), and mass-weighted median radial action (third column) for a set of simulations not shown in Fig. 3.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. B.2. Continuation of Fig. B.1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. B.3. Continuation of Fig. B.2.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. B.4. Continuation of Fig. B.3.</figcaption></figure>
**Appendix C: The overdensity vs. radial action anticorrelation**

As explained in Section 3.1, the difference in the width of the structures observed in the $J_R$ maps compared to those in $\delta_\Sigma(X, Y)$ hinders a clearer discernment of their relationship through visual inspection. This is evident in Figure C.1, where the linear fit, despite significant scatter, shows negative slopes for all the simulations.

<figure><img src="image_17.png" /><figcaption>**Fig. C.1.** Distribution of particles in the correlation diagram $\delta_\Sigma$ vs. $J_R$. The solid lines represent the linear fit of the particle density shown in the background map. The slope of the linear fit is provided in the inset text.</figcaption></figure>
**Appendix D: Additional age binned maps of density contrast**

In this section, we show the density contrast maps per bin of stellar age, $\delta_{\Sigma}(X, Y)$, for those simulations not presented in Section 3.3.

<figure><img src="image_18.png" /><figcaption>**Fig. D.1.** Density contrast maps for three distinct age intervals: stellar particles younger than 1 Gyr (left column), those aged between 1 and 2 Gyr (middle column), and those between 2 and 3 Gyr (right column) for a set of galaxies not shown in Fig. 8. Open circles indicate the spiral arm tracks inferred from Figs. B.1-B.4. Cross in the bottom-right corner indicates the scale of ±5 kpc.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. D.2. Continuation of Fig. D.1.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. D.3. Continuation of Fig. D.2.</figcaption></figure>
## Appendix E: KDE technique at discontinuities

In Sect. 3.2, we explained a discrepancy between the $\delta_{\Sigma}$ and $J_R$ maps, observed in certain regions of the NEXUS simulations, as a misinterpretation of the former as a spiral arm instead of a border effect of the KDE technique. In order to illustrate the mechanisms that produces this effect, we adopted a simplified model of the scenario seen in the simulations in which a medium of constant surface density, $\Sigma_0$, limits the vacuum at $Y = 0$. We can obviate the X direction in this approximation. In the neighbourhood of the boundary between the two media, the estimation of $\delta_{\Sigma}(Y)$ depends on the distance of the kernel centre with respect to $Y = 0$ kpc, as well as on the bandwidths $h$ and $H$ used for the local and global KDEs, respectively. Consequently, there are five possible configurations for the relative location of the kernel (upper panel of Fig. E.1). In the configuration labelled as A, the support of the global kernel is completely inside the dense medium. If the support of the local kernel is completely embedded in the dense medium but not that of the global kernel we obtain the configuration B. The configuration C accounts for the partial embedding of both kernel supports, while in the configuration D only that of the global kernel is partially included in the medium. Finally, in the configuration E both kernels do not overlap with the dense medium.

From the mathematical point of view, the generalisation of the KDE technique with continuous media implies the substitution of the summation in Eq. 8 by an integral. Therefore, the local density estimated with the Epanechnikov kernel of bandwidth $h$ is

$$ \Sigma(Y; h) \sim \frac{\Sigma_0}{h} \int_{c(Y,h)}^{d(Y,h)} \left(1 - \frac{(y - Y)^2}{h^2}\right) dy, \quad (E.1) $$

where the limits of the integral depend of the configurations shown in Fig. E.1. The estimation of the global density is analogous to Eq. E.1 but replacing $h$ by $H$. The density contrast $\delta_{\Sigma}(Y)$ derived from these estimations is represented in the bottom panel of Fig. E.1 for different distances to the interface, assuming bandwidths of $h = 0.5$ kpc and $H = 2.0$ kpc. As can be seen, in configuration B the estimator $\delta_{\Sigma}$ is positively biased while at configuration C it may result either positively (peaking at ~0.5 for this example) or negatively biased. Finally, we adopted $\delta_{\Sigma} = -1$ as value for the configuration E to avoid the 0/0 indeterminate form (dashed red line in Fig. E.1).

<figure><img src="image_21.png" /><figcaption>Fig. E.1. Visualisation of the border effect produced by the KDE technique on the $\delta_{\Sigma}$ map of the NEXUS simulations. Upper panel: scheme of the different configurations of the local (small squares) and global (large squares) scanning windows with respect to the boundary between a medium of uniform density $\Sigma_0$ (grey area at negative Y) and the vacuum (white region at positive Y). Lower panel: values of $\delta_{\Sigma}$ for the mentioned configurations (red curve), assuming $h = 0.5$ kpc (vertical dotted lines) and $H = 2.0$ kpc (vertical dashed lines) bandwidths for the local and global kernels, respectively.</figcaption></figure>