# Azimuthal metallicity variations, spiral structure, and the failure of radial actions based on assuming axisymmetry

Victor Debattista, Tigran Khachaturyants, João Amarante, Christopher Carr, Leandro Beraldo E Silva, Chervin Laporte

► To cite this version:

Victor Debattista, Tigran Khachaturyants, João Amarante, Christopher Carr, Leandro Beraldo E Silva, et al.. Azimuthal metallicity variations, spiral structure, and the failure of radial actions based on assuming axisymmetry. Monthly Notices of the Royal Astronomical Society, 2025, 537 (2), pp.1620-1645. 10.1093/mnras/staf035. obspm-05266963

HAL Id: obspm-05266963

https://hal-obspm.ccsd.cnrs.fr/obspm-05266963v1

Submitted on 19 Sep 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
# Azimuthal metallicity variations, spiral structure, and the failure of radial actions based on assuming axisymmetry

Victor P. Debattista<sup>1,★</sup> Tigran Khachaturyants,<sup>2</sup> João A. S. Amarante<sup>1,2,3★</sup>† Christopher Carr<sup>1,4</sup>,<sup>†</sup> Leandro Beraldo e Silva<sup>5,6</sup> and Chervin F. P. Laporte<sup>3,7,8</sup>

Accepted 2024 December 22. Received 2024 November 18; in original form 2024 February 5

## ABSTRACT

We study azimuthal variations in the mean stellar metallicity, $\langle[\text{Fe/H}]\rangle$, in a self-consistent, isolated simulation in which all stars form out of gas. We find $\langle[\text{Fe/H}]\rangle$ variations comparable to those observed in the Milky Way and which are coincident with the spiral density waves. The azimuthal variations are present in young and old stars and therefore are not a result of recently formed stars. Similar variations are present in the mean age and $\alpha$-abundance. We measure the pattern speeds of the $\langle[\text{Fe/H}]\rangle$-variations and find that they match those of the spirals, indicating that spirals are the cause of the metallicity patterns. Because younger stellar populations are not just more [Fe/H]-rich and $\alpha$-poor but also dynamically cooler, we expect them to more strongly support spirals, which is indeed the case in the simulation. However, if we measure the radial action, $J_R$, using the Stäckel axisymmetric approximation, we find that the spiral ridges are traced by regions of high $J_R$, contrary to expectations. Assuming that the passage of stars through the spirals leads to unphysical variations in the measured $J_R$, we obtain an improved estimate of $J_R$ by averaging over a 1 Gyr time interval. This time-averaged $J_R$ is a much better tracer of the spiral structure, with minima at the spiral ridges. We conclude that the errors incurred by the axisymmetric approximation introduce correlated deviations large enough to render the instantaneous radial actions inadequate for tracing spirals.

**Key words:** Galaxy: abundances – Galaxy: disc – Galaxy: evolution – Galaxy: kinematics and dynamics – galaxies: abundances.

## 1 INTRODUCTION

Observations of Milky Way (MW) disc stars have been finding azimuthal variations in the mean stellar metallicity at fixed radius. Early work using small samples of Cepheids, which are less than 1 Gyr old, found an azimuthal variation in the metallicity at the Solar radius (Luck, Kovtyukh & Andrievsky 2006). On a larger scale, Pedicelli et al. (2009) showed that the relative abundance of metal-rich and metal-poor Cepheids varies across the different quadrants. While Luck et al. (2011) were unable to detect any azimuthal metallicity variations in a sample of 101 Cepheids, Lépine et al. (2011) found azimuthal variations in both Cepheids and open clusters comparable to the radial gradient, and attributed them to the spiral structure.

Large surveys, such as APOGEE and the Gaia data releases (Majewski et al. 2016; Gaia Collaboration 2018) have significantly

accelerated the characterization of azimuthal stellar metallicity variations. Bovy et al. (2014) used an $\sim 10^4$ red clump star sample from APOGEE, with an age distribution favouring 1–4 Gyr, to search for azimuthal metallicity variations. They found azimuthal variations to be smaller than 0.02 dex within an azimuthal volume spanning 45°. Poggio et al. (2022) used Gaia DR3 data to map metallicity variations within 4 kpc of the Sun. They showed that metal-rich stars trace the spiral structure. They also found that these variations are stronger for young stars, but are also present in old stellar populations. Hawkins (2023) used a sample of OBAF stars from LAMOST, and a separate sample of giant stars from Gaia DR3 to measure radial and vertical metallicity gradients. He also found azimuthal metallicity variations of the order of $\sim 0.1$ dex, on average, which are co-located with the spiral structure in the Gaia data (but not in the LAMOST data). Imig et al. (2023) used 66K APOGEE DR17 red giants to suggest that the azimuthal metallicity variations in the Solar Neighbourhood do not seem to track the spiral structure. Using APOGEE DR17, Hackshaw et al. (2024) recently found azimuthal variations in $[\alpha/\text{Fe}]$ of the order of 0.1 dex; they also showed that azimuthal variations seem to be stronger in older populations.

© 2025 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, P.R. China</footnote>

<footnote n="3"><sup>3</sup>Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain</footnote>

<footnote n="4"><sup>4</sup>Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Astronomy, University of Michigan, 500 Church Street, Ann Arbor, MI 48109, USA</footnote>

<footnote n="6"><sup>6</sup>Department of Astronomy & Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="7"><sup>7</sup>LIRA, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, F-92195 Meudon, France</footnote>

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

<footnote>* E-mail: vpdebattista@gmail.com (VPD); joaoant@gmail.com (JASA)</footnote>

<footnote>† Visiting Fellow at UCLan.</footnote>
Theoretical studies of azimuthal metallicity variations attribute them to spiral structure, via a number of different mechanisms. The simplest idea is that the metallicity variations reflect the higher metallicity of recently formed stars, since the gas metallicity itself has an azimuthal variation of ~0.1 dex (Wenger et al. 2019). Comparable variations in gas metallicity have also been found in external galaxies (e.g. Sánchez-Menguiano et al. 2016; Vogt et al. 2017; Ho et al. 2018; Hwang et al. 2019; Kreckel et al. 2019). Spitoni et al. (2019) developed a chemical evolution model of the Galactic disc perturbed by a spiral. They found large metallicity variations near the corotation of the spiral. Solar, Tissera & Hernandez-Jimenez (2020) measured the metallicity of star-forming regions in EAGLE (Schaye et al. 2015) simulated galaxies finding variations of ~0.12 dex. Bellardini et al. (2022) used a suite of FIRE-2 cosmological simulations (Hopkins et al. 2018) to show that the azimuthal variation in the stellar metallicity *at birth* decreases with time, mirroring the evolution of the gas-phase metallicity (Bellardini et al. 2021).

Alternatively, azimuthal variations of the mean stellar metallicity have been viewed as a signature of radial migration. Using an *N*-body simulation of a barred spiral galaxy, in which star particles had metallicity painted on based on their initial radius, Di Matteo et al. (2013) proposed that azimuthal metallicity variations are produced by the radial migration induced by the non-axisymmetric structure. Meanwhile, in a zoom-in cosmological simulation, Grand et al. (2016b) found that the kinematic motions induced by the spirals drive an azimuthal variation in the metallicity, with metal-rich stars on the trailing side and metal-poor ones on the leading side, which they interpreted as a signature of migration of stars from a stellar disc with a radial metallicity gradient. Carr et al. (2022) noted that, in high-resolution *N*-body simulations (Laporte et al. 2018) of the Sagittarius dwarf's encounter with the MW, azimuthal metallicity variations appear in the disc's impulsive response to the passage through the disc, during which migration is also enhanced, and then continuing within the spirals in the periods between passages. Sánchez-Menguiano et al. (2016) observed an offset between the metallicity of the gas and the location of the spiral, which they interpreted as the result of gas migration.

A third hypothesis on the origin of azimuthal metallicity variations posits that they are produced by the different reactions of different stellar populations to a spiral perturbation. Khoperskov et al. (2018) used pure *N*-body (no gas or star formation) simulations of compound stellar discs to show that azimuthal variations can be produced by the varying response of cool, warm and hot stellar populations to spirals. Even in the absence of a metallicity gradient in the disc, if the different populations have different metallicities, azimuthal variations result. They showed that, when migration dominates the azimuthal variations (because all populations had the same radial metallicity profile), the metallicity and density peaks are offset, whereas when the variations are driven by the different responses of different populations then the density and metallicity peaks tend to align.

In this paper, we explore the link between spiral structure and azimuthal metallicity variations. In agreement with Khoperskov et al. (2018), we show, using a model with star formation and self-consistent chemical evolution, that the main driver of these variations are differences in how strongly a stellar population can trace the spiral structure based on its radial action. In Section 2, we describe the simulation used in this paper. In Section 3, we describe the azimuthal chemical variations, both for [Fe/H] and [O/Fe]. We find that variations in the average [Fe/H] closely track the spiral structure. Therefore, in Section 4, we construct a method to track the pattern speeds of the average metallicity variations, taking care

to ensure that the average metallicity is not weighted by the density. Thus, the pattern speeds of the metallicity variations we obtain are independent from those of the density. We show that there is such a close correspondence between the pattern speeds of the average metallicity and of the density in the disc region that it must be the spirals that are driving the metallicity variations. Section 5 explores the idea that the metallicity variations are due to young stars but finds that the azimuthal metallicity variations are present in all age populations, both very young and those older than 6 Gyr. Then, in Section 6, we explore the dependence on the radial action. We find that the radial action provides a very robust explanation for the azimuthal metallicity variations, provided that the radial action is time averaged to reduce the errors introduced by the assumption of axisymmetry. We show that instantaneous values of the radial action instead provide almost exactly the wrong answer, in that stars with seemingly large radial action trace the spiral. In Section 7, we explore whether migration is driving the azimuthal metallicity variations in the model. Multiple lines of evidence suggest that migration is not the dominant factor in setting up the azimuthal variations. In Section 8, we explore whether the temporal changes in the radial action are due to libration about resonances. We show that, even when there are signs of libration, there is also a large jitter in the time-evolution of the instantaneous radial action computed under the assumption of axisymmetry. We present our conclusions in Section 9. The appendices present comparable analyses at other times, explore using medians instead of means, and track a sample of stars in action space chosen to favour those at resonances.

## 2 SIMULATION DESCRIPTION

The model we use here was first described in Fiteni et al. (2021), who referred to it as model M1_c.b. The model evolves via the cooling of a hot gas corona, in equilibrium within a spherical dark matter halo. The dark matter halo has a virial mass $M_{200} = 10^{12}$ M$_{\odot}$ and a virial radius $r_{200} \simeq 200$ kpc, making it comparable to the MW (see the review of Bland-Hawthorn & Gerhard 2016, and references therein). The gas corona has the same profile but with an overall mass that is 10 per cent that of the dark matter. The initial gas corona, as well as the dark matter halo, is comprised of 5 million particles, which means gas particles have an initial mass of $2.8 \times 10^4$ M$_{\odot}$. The dark matter particles come in two mass species, one of mass $2 \times 10^5$ M$_{\odot}$ inside $r = 200$ kpc and the other of mass $7.1 \times 10^5$ M$_{\odot}$ outside. Gas particles are softened with $\epsilon = 50$ pc, while the dark matter particles have $\epsilon = 100$ pc. A cylindrical spin is imparted to the gas particles such that the gas corona has $\lambda \equiv L|E|^{1/2}/(GM_{\text{vir}}^{5/2}) = 0.065$, where $L$ and $E$ are the total angular momentum and energy of the gas, respectively, and $G$ is the gravitational constant (Peebles 1969). After setting up the system, we evolve it adiabatically for 5 Gyr without gas cooling or star formation to ensure that it is fully relaxed.

We then evolve these initial conditions with gas cooling and star formation. As the gas cools, it settles into a rotationally supported disc and ignites star formation. Thus, all stars form directly out of gas. Individual star particles inherit their softening, $\epsilon = 50$ pc, and chemistry from the gas particle from which they form. Star formation, with an efficiency of 5 per cent, requires a gas particle to have cooled below $15\,000$ K, exceeded a density of $0.1$ a.m.u cm$^{-3}$ and to be part of a convergent flow. Star particles are all born with the same mass, $9.3 \times 10^3$ M$_{\odot}$, which is $1/3$ of the (initial) gas particle mass. Gas particles continue to form stars until their mass falls below $5.6 \times 10^3$ M$_{\odot}$ (20 per cent of the initial gas particle mass), at which point a gas particle is removed and its mass distributed to the nearest gas particles. After 13 Gyr, roughly 2.6 million gas particles remain.
Once star formation begins, feedback from supernovae Ia and II, and from AGB winds, couples $4 \times 10^{51}$ erg per supernova to the interstellar medium using the blastwave prescription of Stinson et al. (2006). Thermal energy and chemical elements are mixed amongst gas particles using the turbulent diffusion prescription of Shen, Wadsley & Stinson (2010) with a mixing coefficient of 0.03.

The model is evolved for 13 Gyr with GASOLINE (Wadsley, Stadel & Quinn 2004; Wadsley, Keller & Quinn 2017), an N-body+SPH code.¹Gravity is solved using a KD-tree; we employ a tree opening angle of $\theta = 0.7$ with a base time-step of $\Delta t = 5$ Myr. Time-steps of individual particles are then refined by multiples of 2 until they satisfy the condition $\delta t \equiv \Delta t/2^n < \eta\sqrt{\epsilon/a_g}$, where $a_g$ is the gravitational acceleration of a particle, and we set the refinement parameter $\eta = 0.175$. Time-steps for gas particles are also refined in a similar manner to satisfy the additional condition $\delta t_{gas} < \eta_{coulant}h/[(1+\alpha)c + \beta\mu_{max}]$, where $\eta_{coulant} = 0.4$, $h$ is the SPH smoothing length set over the nearest 32 particles, $\alpha = 1$ is the shear coefficient, $\beta = 2$ is the viscosity coefficient, $c$ is the sound speed and $\mu_{max}$ is the maximum viscous force between gas particles (Wadsley et al. 2004). With these time-stepping recipes, typically 10 rungs (maximum $n = 9$, corresponding to $\delta t = 9766$ yr) are required to move all the particles.

We save snapshots of the simulation every 5 Myr. This high cadence gives a Nyquist frequency of $\sim 300$ km s$^{-1}$ kpc$^{-1}$ (for $m = 2$ perturbations), which is well above any frequency of interest in our MW-sized galaxy. By $t = 13$ Gyr, which is the main snapshot we consider here, the simulation has formed 11 587 120 star particles (2.9M in the radial range $5 \le R/\text{kpc} \le 10$, where we focus our attention), allowing us to dissect the model by stellar population properties while ensuring that we have enough particles for reliable study.

Fiteni et al. (2021) showed that this model experiences an episode of mild clump formation lasting to $\sim 2$ Gyr, which gives rise to a small population of retrograde stars in the Solar Neighbourhood. A bar also forms at $t \simeq 4$ Gyr, which subsequently weakens at $t \simeq 7$ Gyr before recovering. Meanwhile, Khachaturyants et al. (2022a) showed that weak bending waves propagate through the disc, with peak amplitudes $\lesssim 50$ pc in the Solar Neighbourhood. Khachaturyants et al. (2022b) computed the pattern speeds of the density perturbations (bar+spirals), as well as those of the breathing and bending waves in this model. They showed that the breathing waves have the same pattern speeds as the spirals and bar, but very different pattern speeds from the bending waves, from which they concluded that breathing waves are driven by density perturbations. Finally, Ghosh, Debattista & Khachaturyants (2022) used this model at 12 Gyr to show that breathing motions induced by spirals get larger with distance from the mid-plane, and are stronger for younger stars than for older stars.

We focus most of our analysis on the single snapshot of this model at $t = 13$ Gyr. We have verified that results are broadly similar at other times, and in Appendix A we present the azimuthal metallicity variations at other times.

# 3 CHEMICAL AZIMUTHAL VARIATIONS

## 3.1 The azimuthal variation of $\langle[\text{Fe/H}]\rangle$

In Fig. 1, we present the mean metallicity, $\langle[\text{Fe/H}]\rangle$, of the model at 13 Gyr in both Cartesian (top left panel) and cylindrical (top right panel) coordinates. Clear azimuthal variations in $\langle[\text{Fe/H}]\rangle$ are

present. The density contours show a weak bar extending to $R \sim 2$ kpc and open spirals beyond. In the disc region ($R > 5$ kpc), the peak $\langle[\text{Fe/H}]\rangle$ coincides with regions of high density on the spirals.

The left panel of Fig. 2 compares the azimuthal variations of the stellar surface density, $\delta\Sigma/\langle\Sigma\rangle$, with those of the mean metallicity, $\delta\langle[\text{Fe/H}]\rangle$, in the disc region. At these radii ($5.75 \le R/\text{kpc} \le 10.25$), where the spirals dominate, the variations in $\langle[\text{Fe/H}]\rangle$ closely follow those of the density, with peaks and troughs at the same locations. This panel also shows that the $\delta\langle[\text{Fe/H}]\rangle$ of stars older than 2 Gyr behaves very similarly to that of the full distribution, with only slightly weaker amplitudes, indicating that recent metal-rich star formation is not driving the $\langle[\text{Fe/H}]\rangle$ variations. The main variation in both the density and $\langle[\text{Fe/H}]\rangle$ has an $m = 2$ multiplicity. The maximum to minimum variations in $\langle[\text{Fe/H}]\rangle$ are of the order of 0.05–0.1 dex outside the bar radius. In the MW, Hawkins (2023) found metallicity variations of the same order. The top panels of Figs B1 and B2 show maps and azimuthal profiles of the median $[\text{Fe/H}]$, finding a comparable behaviour as for the means.

Since the azimuthal variations in $\langle[\text{Fe/H}]\rangle$ are not limited to the young stars, we proceed by exploring the $\langle[\text{Fe/H}]\rangle$-azimuthal variations of different age populations. Fig. 3 shows $\langle[\text{Fe/H}]\rangle$ in cylindrical coordinates for four age groups. A clear azimuthal variation is present in the very young ($\tau \le 200$ Myr) stars, as seen in the MW (Luck et al. 2006; Lemasle et al. 2008; Pedicelli et al. 2009). A similar pattern is seen in slightly older stars, up to 1 Gyr old. This behaviour might be expected if the spirals are the locii of metal-rich star formation. But a similar azimuthal variation in $\langle[\text{Fe/H}]\rangle$ is evident in the intermediate-age stars ($1\langle\tau/\text{Gyr}\rangle < 6$) and even, to a lesser extent, in the old stars ($6 < \tau/\text{Gyr} < 13$). Thus, the azimuthal variation of $\langle[\text{Fe/H}]\rangle$ is not a result of recent star formation but must be inherent to spiral structure.

In Fig. 4, we compare the azimuthal variations of $\langle[\text{Fe/H}]\rangle$ and the density for each of the age populations and with the full density distribution. The density variations in each population (solid blue lines) track those of the full distribution (dashed black lines); the comparison is noisy in the youngest population (there is a factor of 100 more stars in the oldest population than in the youngest) but seems to hold for this population too. While the amplitude of the density variations becomes smaller with age, we find no significant evidence of a phase shift between the azimuthal variation of the density in populations of different ages. In general, the peak $\langle[\text{Fe/H}]\rangle$ also seems to trace the peak density, particularly at large radii, but the density minima are not always traced by prominent minima in $\langle[\text{Fe/H}]\rangle$.

From the observations that the $\langle[\text{Fe/H}]\rangle$ peaks are radially extended but coincident with the density peaks, we conclude that the cause of the $\langle[\text{Fe/H}]\rangle$ variations is unlikely to be a resonant phenomenon, since azimuthal variations due to stars librating about a resonance might be expected to have $[\text{Fe/H}]$ peaks azimuthally displaced relative to the perturbation (e.g. Grand et al. 2016b; Khoperskov et al. 2018).

## 3.2 The azimuthal variation of $\langle[\text{O/Fe}]\rangle$

The simulation also tracks the evolution of $[\text{O/Fe}]$; the bottom left panel of Fig. 1 shows the azimuthal variation of $\langle[\text{O/Fe}]\rangle$ at 13 Gyr. We find very similar behaviour as for $\langle[\text{Fe/H}]\rangle$, but since $[\text{Fe/H}]$ and $[\text{O/Fe}]$ are anticorrelated, the minima in $\langle[\text{O/Fe}]\rangle$ are at the density peaks.

The right panel of Fig. 2 compares the azimuthal variations of $\langle[\text{O/Fe}]\rangle$ with those of the density (with a minus sign on $\langle[\text{O/Fe}]\rangle$ so maxima and minima are aligned). The peak-to-trough variation is of the order of 0.03 dex at small radii, and declining at larger radii. Overall a match in location of the peaks and troughs is present but

<footnote n="1">¹GASOLINE is available at https://www.gasoline-code.com.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. 〈[Fe/H]〉 in Cartesian (top left panel) and cylindrical (top right panel) coordinates. The bottom left panel shows 〈[O/Fe]〉 in cylindrical coordinates while the bottom right panel shows 〈[Fe/H]〉 excluding stars older than 9 Gyr. In all panels, the contours correspond to the total stellar surface density. The disc is rotating in a counter-clockwise sense, i.e. rotation is in the direction of increasing φ. The model is shown at 13 Gyr.</figcaption></figure>

the correspondence is somewhat less exact than between 〈[Fe/H]〉 and δΣ/⟨Σ⟩.

## 4 PATTERN SPEEDS

Having shown the close correspondence between spirals and 〈[Fe/H]〉 variations at one time, we explore their co-evolution across time by treating the 〈[Fe/H]〉 variations as waves, and comparing their pattern speeds with those of the spirals. We use the spectral analysis code of Khachaturyants et al. (2022b), which uses binned data to measure pattern speeds. We bin the mass and [Fe/H] data in polar coordinates with R ∈ [0, 12] kpc and ΔR = 0.5 kpc, and φ ∈ [0, 360]° and Δφ = 10°. Thus, our binning has N_R × N_φ = 24 × 36 bins. In order to ensure that our pattern speeds for [Fe/H] are not produced by young stars recently born in the spiral arms, which would trivially result in the same pattern speed as the spirals, we exclude all stars younger than 2 Gyr from all the pattern speed analysis. We compute the stellar mass and 〈[Fe/H]〉 in the binned space and produce a 2D polar array at each snapshot. The resulting 2D arrays are M(Rᵢ, φⱼ, tₖ) and fᵢ(Rᵢ, φⱼ, tₖ), where Rᵢ and φⱼ are the coordinates

of the centres of each bin in the array, and tₖ is the time of the snapshot, which are spaced at intervals of Δt = 5 Myr. Following the notation of Khachaturyants et al. (2022a, b), we define the density Fourier coefficients as

$$c_m(R_i, t_k) = \frac{1}{\sum_j M(R_i, \phi_j, t_k)} \sum_j [M(R_i, \phi_j, t_k)e^{im\phi_j}], \quad (1)$$

and similarly for 〈[Fe/H]〉:

$$\gamma_m(R_i, t_k) = N_\phi^{-1} \sum_j [f_e(R_i, \phi_j, t_k)e^{im\phi_j}], \quad (2)$$

where m = 1–4 is the azimuthal multiplicity of the Fourier term. We then use a discrete Fourier transform

$$C_{m,n}(R_i) = \sum_{k=0}^{S-1} c_m(R_i, t_k) w_k e^{-2\pi i k n / S}, \quad (3)$$

with n = 0, ..., S − 1, where S is the number of snapshots in each 1 Gyr baseline (we use S = 200), and the wₖ are the weights of a window function, which we set to Gaussian:

$$w_k = e^{-(k-S/2)^2/(S/4)^2}. \quad (4)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The variation of 〈[Fe/H]〉 (left panel) and − 〈[O/Fe]〉 (right panel) compared with the density variation at a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. The average values for each radius are subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled by the indicated factors for ease of comparison. The dotted horizontal lines indicate the zero for each radius.</figcaption></figure>

The frequency associated with each *n* is given by $\Omega_n = 2\pi n/(mS\Delta t)$. Thus, given we are using time baselines of 1 Gyr $\equiv S\Delta t$, the frequency resolution of our analysis is $\Delta\Omega = 0.98 \times 2\pi/m$ km s$^{-1}$ kpc$^{-1}$ (where the factor of 0.98 is the conversion factor to km s$^{-1}$ kpc$^{-1}$). The power spectrum, defined as the power in each radial bin as a function of frequency, is given by

$$P_m(R_i, \Omega_n) = W^{-1} \left[ |C_{m,n}(R_i)|^2 + |C_{m,S-n}(R_i)|^2 \right], \quad (5)$$

where $W = S \sum_{k=0}^{S-1} w_k^2$.

In order to compute the total (radially integrated) power spectrum for a given *m*, from which we recover the most important pattern speeds, we compute the sum of the mass-squared-weighted power in each radial bin:

$$P_R(\Omega_n) = \frac{\sum_i M^2(R_i) P_m(R_i, \Omega_n)}{\sum_i M^2(R_i)} \quad (6)$$

(where we ignore the multiplicity index *m* on the notation since this will be clear). This mass weighting ensures that the power better reflects the overall importance of a wave, rather than having large relative amplitude perturbations in low-density regions dominating the signal. As we did in Khachaturyants et al. (2022b), we focus on the positive pattern speeds and use the method described in Roškar et al. (2012) to obtain the pattern speeds by fitting Gaussians to $P_R(\Omega_n)$, and subtracting successive peaks.

We adopt a similar analysis for the 〈[Fe/H]〉 pattern speeds, replacing $c_m(R_i, t_k)$ by $\gamma_m(R_i, t_k)$. The spectral analysis of the density (i.e. the measurement of the pattern speeds of the bar and spirals) has already been presented by Khachaturyants et al. (2022a) (see their figs 6 and 7). They found that the density Fourier expansion is dominated by $m=2$, with $m=1$ significantly weaker. In Fig. 5, we plot spectrograms of the density (top row) and 〈[Fe/H]〉 (bottom row) distributions for the time interval $t = 12-13$ Gyr, for $m=1$ to 4. Each block contains a spectrogram (left panel) and the radially

integrated total power, $P_R(\Omega)$ (right panel) of a given wave. For $m \ge 2$, the spectrograms and total power spectra of the density and 〈[Fe/H]〉 waves are in good agreement with each other, with the strongest power in $m=2$ (note the different scales of $P_R$ for different $m$ in Fig. 5). We therefore restrict our attention to $m=2$ from here on.

Fig. 6 therefore shows spectrograms of $m=2$ for the density (left column) and 〈[Fe/H]〉 (right column) at different times in the evolution of the model. A clear correspondence between the pattern speeds in the density and in 〈[Fe/H]〉 is evident at each time. When a peak in the power is present in the density we generally find a corresponding peak for 〈[Fe/H]〉. At the earlier times, the disc hosts a larger number of pattern speeds; the number declines to about three pattern speeds at later times. The growth and slowdown of the bar from 6 to ~9 Gyr and its subsequent weakening by 12 Gyr are also evident in the strong feature at 60–80 km s$^{-1}$ kpc$^{-1}$.

The temporal evolution of the $m=2$ pattern speeds in the density (black crosses) and 〈[Fe/H]〉 (red circles) can be seen in Fig. 7. The most prominent density perturbations, including the bar, but also the slower spirals, have pattern speeds that are well matched by pattern speeds in the metallicity. Conversely, in almost all cases, the pattern speeds identified in 〈[Fe/H]〉 are matched by pattern speeds in the density. The size of the markers represents the power normalized by the maximum power value each pattern reaches over the entire time interval of this analysis, i.e. 5–13 Gyr. The density and 〈[Fe/H]〉 patterns fail to match only for the weakest waves (smallest markers).

Khachaturyants et al. (2022b) showed that there is a strong connection between spirals and breathing waves by showing that the pattern speeds of strong breathing waves are always matched by pattern speeds of density waves, and vice versa. They argued therefore that the breathing motions must be part of the nature of spirals. As in Khachaturyants et al. (2022b), we quantify the significance of the match between the density and the 〈[Fe/H]〉 waves
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The distribution of stellar $\langle [\text{Fe/H}] \rangle$ at 13 Gyr split by age, as indicated. White represents bins with $\langle [\text{Fe/H}] \rangle$ off the scale (inside ~ 2 kpc) or with less than 10 particles (elsewhere). In all panels, the identical contours correspond to the total surface density. The disc is rotating in a counter-clockwise sense, i.e. rotation is in the direction of increasing $\phi$.</figcaption></figure>

in Fig. 8 by showing the fraction of identified density pattern speeds that are matched by pattern speeds in $\langle [\text{Fe/H}] \rangle$ (solid lines). The match fraction is defined as the ratio of the number of density pattern speeds with metallicity matches to the total number of density pattern speeds, i.e.

$$f = \frac{N(\Omega_{\text{dens+met}})}{N(\Omega_{\text{dens}})}, \qquad (7)$$

where $N(\Omega_{\text{dens}})$ ($N(\Omega_{\text{dens+met}})$) is the number of density pattern speeds above the threshold (the number of density pattern speeds above the threshold that are matched by a metallicity pattern speed) for the full time interval 5–13 Gyr (not just a single time). In order to avoid including the bar, we consider only $\Omega_p < 49 \text{ km s}^{-1} \text{ kpc}^{-1}$ for each multiplicity and baseline. This fraction is shown for different minimum power thresholds of the density pattern speeds, $\min(P_\Sigma)$. The different panels show the matches for different wave multiplicities. Better matches between the density and $\langle [\text{Fe/H}] \rangle$ waves occur for $m > 1$. We check if these match fractions are compatible with random matches via a bootstrap algorithm where we select 2500 iterations of random pattern speeds, with a time distribution exactly matching those of the $\langle [\text{Fe/H}] \rangle$ waves, and compute the match

fraction with the density pattern speeds as before. Fig. 8 presents the 5–100 per cent probability intervals (shaded regions) of the match fraction for these random samples. Only in the case of the $m = 1$ waves are the match fractions compatible with the random matches, while $m > 1$ match fractions are clearly higher than random at all $\min(P_\Sigma)$. Since we have binned the $\langle [\text{Fe/H}] \rangle$ field, and computed the Fourier time series only of the binned values, the metallicity time series is not contaminated by the density. As a result, the pattern speeds of the two would be independent of each other if the two were unrelated. The high match probability we find therefore indicates that the $\langle [\text{Fe/H}] \rangle$ variations are due to the spirals. Thus, we conclude that the $\langle [\text{Fe/H}] \rangle$ waves are driven by density waves and, in the disc, we must identify them with spirals. It remains to be understood how spirals cause the $\langle [\text{Fe/H}] \rangle$ variations.

## 5 DISSECTION BY AGE

Having shown that the $\langle [\text{Fe/H}] \rangle$ variations are part of the spirals, and do not result from recent star formation, we next explore the physical mechanism of the $\langle [\text{Fe/H}] \rangle$ variations.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The variation of 〈[Fe/H]〉 compared with the density variation, split by age for the same populations as in Fig. 3, at a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. In each panel, the solid lines refer to 〈[Fe/H]〉 (red) and population density (blue) while the dashed solid lines refer to the density of all stars. The average values for each radius are subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled by the same factor of 0.2 for ease of comparison. The dotted horizontal lines indicate the zero for each radius.</figcaption></figure>

First, we consider whether 〈[Fe/H]〉 variations reflect variations in the mean age. Fig. 9 shows the age–metallicity relation (AMR) of stars between 5 and 10 kpc. A quite shallow relation is evident for stars younger than ~ 8 Gyr, getting steeper for older stars. The scatter in [Fe/H] increases with age. We have also checked that the AMR

of different subsections of this radial range look reasonably similar. The AMR of the model for stars younger than ~ 5 Gyr is shallow but not perfectly flat. The AMR in the Solar Neighbourhood is perhaps similar (Feuillet et al. 2019; Anders et al. 2023; Johnson et al. 2024), although some studies find a flatter AMR (Mackereth et al. 2017;
<figure><img src="image_6.png" /></figure>

**Figure 5.** Power spectra (spectrograms) for density (top panel) and 〈[Fe/H]〉 (bottom panel) waves of different multiplicities (columns) at *t* = 12 – 13 Gyr. The solid lines show the rotation curve, Ω(*R*), evaluated at the central time of the interval using an interpolated potential. Dashed lines indicate Ω ± *κ*/m. To the right of each spectrogram, we show the radius-integrated power. The colour scale limits are the same for each wave type over all multiplicities and models. The pattern speeds of *m* > 1 density and 〈[Fe/H]〉 waves match each other. Only stars older than 2 Gyr are considered.

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

**Figure 6.** Similar to Fig. 5, but for the power spectra of density (left panel) and 〈[Fe/H]〉 (right panel) *m* = 2 waves at the different time intervals (rows) indicated at the top right of each panel. Only stars older than 2 Gyr are considered.

Delgado Mena et al. 2019). The existence of an AMR implies that, if populations of different age have spirals of different amplitude, then azimuthal variations in 〈[Fe/H]〉 would arise from the azimuthal variation of the spiral density (e.g. Khoperskov et al. 2018). The

top panel of Fig. 10 plots the radial profiles of the *m* = 2 amplitude of the density; it shows that spirals at *R* > 5 kpc are present in all age populations, although they are weaker in older populations, as seen already in Fig. 3 (see also Ghosh et al. 2022). The bottom panel shows that the older populations support weaker spirals because their radial velocity dispersion, σ<sub>R</sub>, is higher.

Fig. 11 shows maps of the mean age, 〈τ〉. The top panel, which includes all stars, shows age variations, which have similar spatial distribution as the metallicity variations of Fig. 1. The bottom panel shows 〈τ〉 excluding stars younger than 2 Gyr. The variations in 〈τ〉 retain a similar character to those of the full population, showing again that recent star formation is not the cause of the azimuthal variations. As expected, the spiral ridges are regions of younger stars, which follows because these stars have lower radial velocity dispersions and thus host stronger spirals (Fig. 10).

Fig. 12 compares the azimuthal variation of 〈τ〉 and the density at different radii. Age variations of the order of 0.5 – 1 Gyr are present at all radii. The green lines in Fig. 12 show the azimuthal profiles of 〈τ〉 when stars younger than 2 Gyr are excluded. The azimuthal variations are now weaker, but still evident. Moreover, these profiles largely track those of all stars where the deviation from the radial average is significant. Note that the density peaks correspond to locations where the stars are, on average, younger.

While there is a reasonably good correspondence between the variation of the density and mean age, comparison with Fig. 2 shows that 〈[Fe/H]〉 traces the density variations slightly better than 〈τ〉 does. This is particularly evident at larger radii. 〈τ〉 therefore is probably not directly the driver of the 〈[Fe/H]〉 variations. In any case, the variations in 〈τ〉 are too small, given the AMR, to be directly responsible for the variations in 〈[Fe/H]〉.

# 6 DEPENDENCE ON THE RADIAL ACTION

We have shown above that older populations with larger σR have weaker spirals. The radial motions of individual stars are best quantified by the radial action, J<sub>R</sub>, which is an integral of motion when a system evolves adiabatically. We therefore expect to see strong correlations between the azimuthal variations of 〈[Fe/H]〉
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The most prominent pattern speeds for the *m* = 2 density (black crosses) and 〈[Fe/H]〉 (red circles) waves. The marker size is indicative of the radius-integrated power at a given pattern speed. The peaks are scaled by the highest peak in the particular type of wave (density or 〈[Fe/H]〉) over all baselines, *P*₂/max(*P*₂) with the legend showing the relationship between the marker size and *P*₂/max(*P*₂). Only pattern speeds with *P*₂/max(*P*₂) ≥ 0.01 are shown. We observe a correspondence between the *m* = 2 density and 〈[Fe/H]〉 perturbations: for the most part, the pattern speeds of the two match, except for the weakest patterns. The horizontal dashed line shows Ω = 49 km s⁻¹ kpc⁻¹; the region above this frequency is dominated by the bar. Each marker is shown at the mid-point of the corresponding Fourier time series; since these are 1 Gyr long, consecutive points are correlated (data points in the second half of each time series are the first half of the next time series). Only stars older than 2 Gyr are considered.</figcaption></figure>

and those of 〈*J*<sub>R</sub>〉. To test this idea, we compute the actions of the model using the axisymmetric Stäckel fudge of Binney (2012), as implemented in the AGAMA package (Vasiliev 2019).

## 6.1 Instantaneous radial actions

The top panel of Fig. 13 shows the *J*<sub>R</sub>-[Fe/H] relation. Across the range in *J*<sub>R</sub> the mean metallicity varies significantly, and nearly linearly, although the scatter about the mean is large. The bottom panel shows the *τ*-*J*<sub>R</sub> relation. The average *J*<sub>R</sub> grows linearly with age, comparable to the trend seen in the MW (Beane, Ness & Bedell 2018; Gandhi & Ness 2019), with 〈*J*<sub>R</sub>〉 at 10 Gyr roughly double that measured in the MW. These relations suggest that the driver of the azimuthal variations of stellar populations may be *J*<sub>R</sub>, which both determines how strong the spiral structure is in a given population, and what the mean metallicity is.

The top panel of Fig. 14 maps the average radial action, 〈*J*<sub>R</sub>〉, excluding stars younger than 2 Gyr. Significant variations are evident, with 〈*J*<sub>R</sub>〉 maxima near the locations of density peaks. However, at some radii (e.g. 6.5 kpc), the location of the 〈*J*<sub>R</sub>〉 peaks are clearly offset from the location of the density peaks. The left panel of Fig. 15, which compares the azimuthal profiles of the density and of 〈*J*<sub>R</sub>〉, also shows that the peaks and troughs are offset relative to each other, with the offset varying with radius.

Figs 14 (top panel) and 15 (left panel) are puzzling because they imply that, broadly speaking, the spiral density peaks are better traced

by stars with large *J*<sub>R</sub>, which should be older, and therefore represent [Fe/H]-poor stars. However, we have seen that the 〈[Fe/H]〉-peaks are located at the density peaks. Moreover, in Fig. 10, we showed that older populations are hotter and host weaker spirals, contradicting the behaviour seen in *J*<sub>R</sub>.

We resolve this difficulty in the next section by showing that the inherent assumption of axisymmetry in computing *J*<sub>R</sub> gives rise to correlated errors.

## 6.2 Time-averaged radial actions

Because AGAMA assumes axisymmetry when computing actions, subtle but spurious correlations between the azimuthal phase of a star with respect to a spiral and its computed *J*<sub>R</sub> may arise. In order to mitigate, somewhat, such effects, we also compute the time average of *J*<sub>R</sub>, which we denote 〈*J*<sub>R</sub>〉<sub>t</sub>. We carry out these time averages over a 1 Gyr time interval. This is long enough for stars to have drifted across spirals several times. For those stars trapped at the corotation resonance, 1 Gyr is long compared to the lifetimes of transient spirals (see Roškar et al. 2012), so that they will not have been trapped throughout this time. We compute the actions at 21 snapshots, between 12 and 13 Gyr, inclusive, and average *J*<sub>R</sub> to obtain 〈*J*<sub>R</sub>〉<sub>t</sub>.

We compare the time-averaged radial action to the instantaneous one (at 13 Gyr) in the top panel of Fig. 16. The mean of the distribution (solid line) is slightly offset towards *J*<sub>R</sub> > 〈*J*<sub>R</sub>〉<sub>t</sub>; this is
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Match fraction between the density and ⟨[Fe/H]⟩ waves, given by equation (7), versus the minimum power of the density waves, min(*P*<sub>Σ</sub>). The shaded region indicates the 5–100 per cent distribution of the match fractions resulting from a bootstrap resampling of random pattern speeds of ([Fe/H]) waves. All random bootstrap samples have match fraction lower than in the simulation for *m* > 1, implying the matches are not random. For this analysis, we only consider pattern speeds Ω ≤ 49 km s<sup>-1</sup> kpc<sup>-1</sup> (the dashed horizontal line in Fig. 7), to exclude the bar. We do not show min(*P*<sub>Σ</sub>) with less than 10 identified pattern speeds across all times. For *m* > 1, we note an increase in match fraction with increasing min(*P*<sub>Σ</sub>) that substantially exceeds the random bootstrap-sampled match fractions. Only stars older than 2 Gyr are considered in this analysis.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The age–metallicity relation for star particles at $5 \le R$/ kpc $\le 10$ at 13 Gyr. The solid line shows ([Fe/H]) while the dashed lines show the 16 and 84 per cent distributions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Top panel: the radial profiles of the amplitude of the *m* = 2 density Fourier moments for different age populations, as indicated in the legend. Outside 5 kpc, spiral structure is weaker in old populations than in younger ones. Bottom panel: the radial velocity dispersion, $\sigma_R$, profiles of the same stellar populations. Stars younger than 1 Gyr are not plotted.</figcaption></figure>

possibly the effect of stellar orbit heating during this time interval. We can indeed see, in the bottom panel of Fig. 10, that stellar populations are slowly heating. The distribution is broadened, spreading across the 1:1 line, suggesting that the instantaneous action is varying significantly over the 1 Gyr time interval. The most likely cause of this variation is that the assumption of axisymmetry is introducing errors in the measurements.

The distribution in the top panel of Fig. 16 appears featureless. In contrast, the comparable plot for $J_\phi - \langle J_\phi \rangle_t$ (bottom panel) shows clear features over the same time interval. These features are likely produced by the radial migration driven by transient spirals (Sellwood & Binney 2002).

To confirm the hypothesis that the instantaneous value of $J_R$ is being perturbed by the spirals, we measure the overdensity of individual particles, $\Sigma/\langle\Sigma\rangle$, where $\langle\Sigma\rangle$ is the azimuthally averaged density. We measure the local density at each particle by constructing a KD tree of the star particle distribution based on the Euclidean 2D mid-plane distances using the SCIPY.SPATIAL.KDTREE package. For each particle, we measure the distance to the nearest 32 particles and compute the density over these particles, and then compute the ratio of the particle's density to the azimuthal average. We plot the $\langle J_R \rangle_t - J_R$ distribution coloured by this density excess in Fig. 17. The deviation of the instantaneous $J_R$ from the diagonal is clearly very well correlated with the excess density, with stars having $J_R > \langle J_R \rangle_t$ ($J_R \langle\langle J_R \rangle\rangle_t$) in high (low) density regions, indicating that $J_R$ has correlated errors.

In order to compare the measurement uncertainty with the heating rate, we first bin particles by $J_R$ at 13 Gyr. For each bin, we estimate the uncertainty in the measurement of $J_R$ as the standard deviation in $\langle J_R \rangle_t$ at a given $J_R$, which we term $\sigma_{J_R}$. We estimate the heating rate as the difference between $J_R$ and the average $\langle J_R \rangle_t$ in each bin, i.e. $\Delta J_R = \langle J_R - \langle J_R \rangle_t \rangle$. Fig. 18 shows our measurements of $\Delta J_R$.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Mean age, $\langle\tau\rangle$, in cylindrical coordinates. Top panel: for all stars. Bottom panel: for stars older than 2 Gyr. The contours correspond to the surface density. The disc is rotating in a counter-clockwise sense, i.e. rotation is in the direction of increasing $\phi$. The model is shown at 13 Gyr.</figcaption></figure>

and $\sigma_{J_R}$. Over this 1 Gyr time interval, $\sigma_{J_R}$ is smaller than $\Delta J_R$, i.e. the uncertainty in measuring the instantaneous $J_R$ introduced by the assumption of axisymmetry is larger than the radial heating. In reality, we have overestimated the heating, since this is based on the assumption that all the difference between $\langle J_R \rangle$ and $\langle\langle J_R \rangle\rangle$ (corresponding to the difference between the 1 : 1 (dashed) line and the mean $J_R$ (solid) line in the top panel of Fig. 16) is due entirely to heating. However, the dotted line in the same panel shows the relation between the same $\langle J_R \rangle$ and $J_R$ now computed at 12 Gyr; this line is also below the diagonal. The offset between the diagonal and the solid line, therefore, is partly due to the intrinsically skewed distribution of $\langle J_R \rangle$, which means that we have overestimated the heating rate. We have also estimated the heating rate by averaging $J_R$ between 12 and 12.5 Gyr and between 12.5 and 13 Gyr, which also shows that $\Delta J_R$ overestimates the heating rate. Crucially then, we are not mistaking actual (physical) heating, from whatever source, for the (unphysical) scatter arising from the assumption of axisymmetry.

The middle panel of Fig. 13 shows the $\langle J_R \rangle$-[Fe/H] relation; this is very similar to the $J_R$-[Fe/H] relation, since the errors in

<figure><img src="image_13.png" /><figcaption>**Figure 12.** The variation of $\langle\tau\rangle$ compared with the density variation at a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. (We use $-\langle\tau\rangle$ to permit easier comparison with the density.) The average value at each radius is subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled by the factor 3.5 for ease of comparison. The dotted horizontal lines indicate the zero for each radius. The red (green) lines show $\langle\tau\rangle$ profiles of all stars (stars older than 2 Gyr).</figcaption></figure>

$J_R$ should be independent of [Fe/H]. Thus, the azimuthal variations of $\langle[\text{Fe/H}]\rangle$ may be driven by those of $\langle J_R \rangle$. The bottom panel of Fig. 14 maps $\langle\langle J_R \rangle\rangle$ in cylindrical coordinates. The spiral arms are very well delineated by $\langle\langle J_R \rangle\rangle$, much better than by $\langle J_R \rangle$ (top panel). The differences between the two maps are quite striking with the peak densities being associated with minima in $\langle\langle J_R \rangle\rangle$ but closer to maxima in $\langle J_R \rangle$. That low $\langle J_R \rangle$ should track the spirals is natural since radially cool stars support the spiral structure better. At the same time, the larger radial perturbations these stars experience means that when the actions are computed they will appear to be radially hotter, leading to apparently larger $J_R$.

The right panel of Fig. 15 plots the azimuthal variations of $\langle\langle J_R \rangle\rangle$ and compares them with those of the density. Compared with $\langle J_R \rangle$ in the left panel, a number of important differences are evident. First, we confirm that, while the variations in $\langle J_R \rangle$ somewhat correlate with the density variations, those in $\langle\langle J_R \rangle\rangle$ anticorrelate with the density variations, as expected. Furthermore, while the correlation with $\langle J_R \rangle$ is relatively poor, that with $-\langle\langle J_R \rangle\rangle$ is quite strong. Combined with the $\langle J_R \rangle$-[Fe/H] relation of Fig. 13 $\langle J_R \rangle$ is able to explain the strong correlation between the density and $\langle[\text{Fe/H}]\rangle$ variations.

Fig. 19 shows the Fourier amplitudes of $m=2, 4$ and $6$ of the stellar surface density at $R=6$ kpc (where the spirals are strong) separated into 3 roughly equal-sized populations of low, medium and high $J_R$ and $\langle J_R \rangle$. The low $\langle J_R \rangle$ population has the strongest amplitudes, including at $m=6$, indicating that this population is able to support sharp features in the spiral structure. The mid and high $\langle J_R \rangle$ have progressively weaker amplitudes. Instead, in $J_R$, the strongest amplitudes are in the high-$J_R$ population, although these amplitudes are smaller than the ones in the low-$\langle J_R \rangle$ population. The amplitudes decrease as $J_R$ decreases, the opposite of what is expected. Moreover, the $m=6$ amplitude is nearly zero for all $J_R$.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** The $J_R$–[Fe/H] (top panel), $\langle J_R \rangle_t$–[Fe/H] (middle panel), and $\tau - J_R$ relations (bottom panel) for star particles at $5 \le R/\text{kpc} \le 10$ at 13 Gyr. The solid lines show $\langle [\text{Fe/H}] \rangle$ (top and middle panels) and $\langle J_R \rangle$ (bottom panel). In all panels, the dashed curves show the 16 and 84 per cent distributions. The middle panel excludes stars younger than 2 Gyr, which results in the paucity of low $\langle J_R \rangle_t$ stars.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** The mean instantaneous radial action, $\langle J_R \rangle$, (top panel) and the mean time-averaged radial action, $\langle \langle J_R \rangle_t \rangle$, (bottom panel) in cylindrical coordinates. The contours correspond to the surface density. Only stars older than 2 Gyr have been considered. The disc is rotating in a counter-clockwise sense, i.e. rotation is in the direction of increasing $\phi$. The model is shown at 13 Gyr.</figcaption></figure>

populations, indicating that sharp features in physical space cannot be represented by $J_R$ populations. Clearly, $J_R$ does a poorer job of tracing the spiral structure than does $\langle J_R \rangle_t$. The middle and bottom panels of Fig. B1 show that the median $\langle J_R \rangle_t$ also traces the spirals better than the median $\langle J_R \rangle$. The bottom panel of Fig. B2 shows that the azimuthal variations of the median $\langle J_R \rangle_t$ track those of the density very well.

Thus, while the time-averaged actions cannot be the true actions at any given time, they are a sufficiently close approximation to be useful for studying the effects of spiral structure.

## 6.3 Correlations

We compute the correlations between the variations in the density, and those of the population properties for stars in the disc region. In other words we quantify the correlations seen in Fig. 2, and others like it. From the cylindrical maps (in the radial range $5 \le R/\text{kpc} \le 10$, which covers $N_\phi = 101$ and $N_R = 51$ bins, giving $N_\phi \times N_R = 5151$ pixels), we subtract the radial profile of the azimuthal average in each
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The variation of $\langle J_R \rangle$ (left panel) and $-\langle\langle J_R \rangle_t\rangle$ (note the sign difference) compared with the density variation at a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. The average value at each radius is subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled by the factor 50 for ease of comparison. The dotted horizontal lines indicate the zero for each radius. The red lines show $\langle J_R \rangle$ profiles of all stars (left panel), while the green lines show stars older than 2 Gyr only (both panels).</figcaption></figure>

of the quantities of interest ($\Sigma$, [Fe/H], [O/Fe], $\tau$, $J_R$, and $\langle J_R \rangle_t$). We then compute the correlation coefficient between all the pixel values for the variations in the density and each of the other variables in turn.

Fig. 20 shows how stellar population properties vary with density variations. The spiral ridges correspond to high $\delta\langle[\text{Fe/H}]\rangle$, and low $\delta\langle[\text{O/Fe}]\rangle$ and $\delta\langle\tau\rangle$. It also shows that the correlation between density variations and $\delta\langle\langle J_R \rangle_t\rangle$ is stronger than that between density variations and $\delta\langle J_R \rangle$. This latter relation is, in fact, even weaker than the correlation between the density variations and $\delta\langle[\text{Fe/H}]\rangle$ and $\delta\langle\tau\rangle$. Even more strikingly, the correlation between the density variations and $\delta\langle\langle J_R \rangle_t\rangle$ is negative (i.e. high-density regions have small time-averaged radial action) while the correlation is positive for the instantaneous one. Fig. 21 shows how the variations in the stellar population properties (age, metallicity, and $\alpha$-abundance) correlate with the variations in $J_R$ and $\langle J_R \rangle_t$. The stellar population variations are poorly correlated with $\delta\langle J_R \rangle$ but strongly correlated with $\delta\langle\langle J_R \rangle_t\rangle$. Together, Figs 20 and 21 demonstrate that the errors associated with the axisymmetric approximation in computing the actions are sufficiently large to mask the actual dependence of the population variations on the radial action.

We have also explored the vertical action, $J_z$ and $\langle J_z \rangle_t$ (defined similar to $\langle J_R \rangle_t$). We find that $J_z$ behaves more similar to $\langle J_z \rangle_t$, with spiral structure becoming progressively weaker as $J_z$ increases. This indicates that the vertical motion is less affected by the spirals than the radial motion, consistent with the idea that in-plane and vertical motions are only weakly coupled.

## 6.4 Synthesis

In summary, the azimuthal $\langle[\text{Fe/H}]\rangle$ variations appear to be the result of the spirals being more strongly supported by stars with low intrinsic $J_R$. Such stars are, on average, metal-rich, and relatively young, which means that the locations of the spiral density peaks (troughs)

are also where these stars are over (under) represented, giving rise to the mean metallicity peaks (troughs). However, recovering this behaviour requires that the radial actions are corrected for the error from the assumption of axisymmetry; here, we have partially done this by time-averaging $J_R$.

## 7 THE ROLE OF MIGRATION

Trapping at the corotation resonance of transient spirals drives stellar migration (Sellwood & Binney 2002; Grand, Kawata & Cropper 2012; Roškar et al. 2012). Azimuthal variations in $\langle[\text{Fe/H}]\rangle$ have been proposed to result from radially migrating populations. We show that migration is not a leading cause of the azimuthal $\langle[\text{Fe/H}]\rangle$ variations by showing that populations with a stronger radial gradient do not produce an azimuthal offset, and that non-migrating populations also exhibit azimuthal variations. Additional support comes from considering the phase difference between the $\langle[\text{Fe/H}]\rangle$ and density waves.

### 7.1 Populations with strong radial gradients

Khoperskov et al. (2018) emphasize that migration-dominated $\langle[\text{Fe/H}]\rangle$ variations require a radial metallicity gradient. Fig. 22 shows the $\langle[\text{Fe/H}]\rangle$ profiles of different age populations in the model. This shows the usual influence of radial migration, in that metallicity profiles become flatter for older populations, which is the opposite trend expected for a disc growing inside-out. A comparable behaviour has been observed in the MW (Nordström et al. 2004; Casagrande et al. 2011; Bergemann et al. 2014; Vickers, Shen & Li 2021; Anders et al. 2023; Willett et al. 2023), including the slightly positive gradient of the oldest stars (Nordström et al. 2004; Casagrande et al. 2011; Anders et al. 2023). The overall $\langle[\text{Fe/H}]\rangle$ profile (dashed cyan line) is quite flat over $5 \le R/\text{kpc} \le 10$, which might explain why migration
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparison of $\langle J_R \rangle_t$ and $J_R$ (top panel) and $\langle J_\phi \rangle_t$ and $J_\phi$ (bottom panel) for stars which are at $5 \le R$/kpc $\le 10$ at 13 Gyr. The dashed lines show the diagonal, $J_x = \langle J_x \rangle_t$, while the solid lines show the mean relation. The radial action shows the effect of modest heating while the angular momentum shows the effect of migration via resonances. Stars younger than 1 Gyr are not included. In the top panel, the dotted line shows the mean $J_R$ if we instead use $J_R$ at 12 Gyr on the horizontal axis.</figcaption></figure>

has not seemed to matter thus far, as implied by the coincidence of the $\delta\langle[\text{Fe/H}]\rangle$ and density peaks. However, if we consider only stars younger than 9 Gyr, to match the bulk of the MW's thin disc, we find a metallicity gradient of $\nabla\langle[\text{Fe/H}]\rangle \simeq 0.03 \text{ dex kpc}^{-1}$ (dashed black line). This gradient, although roughly half that of the MW's thin disc (e.g. Vickers et al. 2021; Imig et al. 2023; Hackshaw et al. 2024), allows us to test whether stellar migration of populations with a radial metallicity gradient may dominate azimuthal $\langle[\text{Fe/H}]\rangle$ variations. In the bottom right panel of Fig. 1, we map $\langle[\text{Fe/H}]\rangle$ for stars younger than 9 Gyr, superimposing the density map of all stars. No obvious offsets between the peak $\langle[\text{Fe/H}]\rangle$ and peak density are evident, and the overall metallicity map looks broadly consistent with that for all stars, apart from the different $\langle[\text{Fe/H}]\rangle$ values. The same result can be seen comparing the cyan line with the red or green lines in the left panel of Fig. 2, suggesting that migration does not dominate the azimuthal $\langle[\text{Fe/H}]\rangle$ variations.

<figure><img src="image_18.png" /><figcaption>**Figure 17.** The density variation of individual stars at $t = 13$ Gyr in $(J_R, \langle J_R \rangle_t)$-space. The deviation of $J_R$ from the diagonal line is highly correlated with the overdensity. Only bins with more than 20 particles are shown.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** The heating, $\Delta J_R = \langle J_R - \langle J_R \rangle_t \rangle$ (blue), and the instantaneous radial action measurement uncertainty, $\sigma_{J_R}$, computed as the standard deviation of $(J_R - \langle J_R \rangle_t)$ (red). Only stars older than 1 Gyr in the disc are considered. See the text for why $\Delta J_R$ overestimates the heating.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** The $m = 2, 4$, and 6 Fourier amplitudes of the density distribution at $5.5 < R$/kpc $< 6.5$ separated into equal-sized low, medium and high populations of $\langle J_R \rangle_t$ (solid lines) and $J_R$ (dashed lines).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 20.** The dependence of azimuthal variations in [Fe/H] (top row), [O/Fe] (second row), τ (third row), *J*<sub>R</sub> (fourth row), and ⟨*J*<sub>R</sub>⟩<sub>t</sub> (bottom row) on the azimuthal variations in density. The Pearson *r* value of each correlation is indicated in each panel. The strongest correlation is between the density variations and those in the ⟨*J*<sub>R</sub>⟩<sub>t</sub>. Note the change in the sign of *r* for the correlation between the density variations with δ⟨*J*<sub>R</sub>⟩ and with δ⟨⟨*J*<sub>R</sub>⟩<sub>t</sub>⟩. Populations younger than 2 Gyr are excluded from this plot. The units of δ⟨τ⟩ are Gyr, while the units of δ⟨*J*<sub>R</sub>⟩ and δ⟨⟨*J*<sub>R</sub>⟩<sub>t</sub>⟩ are kpc km s<sup>-1</sup>.</figcaption></figure>

## 7.2 Non-migrating stars

We further explore whether migration plays a leading role in the metallicity variations by studying a non-migrating population over the time interval 12.65–13 Gyr, during which a strong spiral forms. We compute the change in angular momentum, $\Delta J_{\phi} = J_{\phi}(12.65 \text{ Gyr}) - J_{\phi}(13 \text{ Gyr})$. Fig. 23 plots the distribution of $\Delta J_{\phi}$ versus $J_{\phi}$. This exhibits the usual stripes associated with spiral corotation migration (Sellwood & Binney 2002; Roškar et al. 2012). The dotted horizontal lines delineate $|\Delta J_{\phi}| < 100 \text{ kpc km s}^{-1}$, which we take as a non-migrating population (constituting ~70 per cent of all stars).

Fig. 24 shows azimuthal profiles of $\delta([\text{Fe/H}])$ for all stars, and for the non-migrating population. We present this twice, once for all ages (other than those that form after 12.65 Gyr for which we cannot compute $\Delta J_{\phi}$), which we show with dashed lines, and excluding stars older than 9 Gyr (solid lines). The $\langle [\text{Fe/H}] \rangle$ profiles of all stars (shown in red) follow the density variations quite well, regardless of whether stars older than 9 Gyr are excluded (solid lines) or included (dashed lines). When non-migrating stars are considered, the profiles still somewhat track those of the density but there is now considerably more noise in the profiles, as attested by the point-to-point scatter. However, for the most part, the blue solid line follows the density variations, most notably peaking at the same locations. If we include stars older than 9 Gyr, the non-migrating population (blue-dashed line) still generally peaks near the density peaks.

## 7.3 Evolution of the phase difference

Our final evidence that migration does not dominate the azimuthal variations comes from considering the time-evolution of the relative phase of the density and metallicity. We compute the phases of the waves using the Fourier analysis of Section 4. Since star formation occurs preferentially on the leading (trailing) side of spirals outside (inside) corotation, which leads to phase offsets for young stars, as in Section 4, we exclude stars younger than 2 Gyr in this analysis. From the Fourier moments at each timestep, we compute the phase shift between the density and $\langle [\text{Fe/H}] \rangle$ waves, $\Delta\phi_m = \phi_{[\text{Fe/H}],m} - \phi_{\Sigma,m}$. Because single phases are noisy, particularly when the spirals are weak, or when the two spirals which overlap in the disc have similar amplitudes, we average over 1 Gyr time intervals. Fig. 25 shows $\Delta\phi_m$ for $1 \le m \le 4$. The phase shift for $m=1$ perturbations is large, but the density waves of this multiplicity are weak. For $m=2$, which has the largest amplitudes, $|\Delta\phi_2|$ is everywhere less than $10^\circ$. Fig. 25 includes $\Delta\phi_m$ for $R=1$ kpc, where the bar is resident, to show that even in this case there is some deviation from $\Delta\phi_m = 0^\circ$, which therefore indicates some of the uncertainty in the measurement method. No sustained offset, $|\Delta\phi_2|$, is evident; the same is true for the $m=3$ and 4 moments. Thus, we find no evidence that the $\langle [\text{Fe/H}] \rangle$ waves have a consistent phase shift relative to the spirals, even though the pattern speeds of the main spirals are only slowly evolving over this time. We have repeated this analysis excluding stars younger than 2 Gyr and those born in the first 4 Gyr (to ensure a significant disc [Fe/H] radial gradient) and still find no consistent phase shifts.

We conclude therefore that, at least for this simulation, there is no support for the idea that migration dominates the azimuthal [Fe/H] variations. We discuss this result, and reasons why the model might be underestimating the role of migration relative to the effect of varying radial actions, in Section 9.

# 8 THE VIEW IN ACTION SPACE

Resonant trapping leads to $J_R$ and $J_\phi$ non-conservation, with stars librating around the resonance when viewed in the space of axisymmetric actions. In this section, we examine whether the $J_R$ variations are due to such librations, rather than errors arising from the axisymmetric approximation. We consider the joint time-evolution of $J_R$ and $J_\phi$, in order to show that the variations in $J_R$ we have found are not due to resonances.

Recall that when a resonance is present, stars are either trapped by it, in which case their actions librate back and forth across the resonant action, or they are untrapped, in which case their actions vary periodically over their orbits. The variations of $J_R$ and $J_\phi$ of both trapped and untrapped stars can be related by noting that the Jacobi
<figure><img src="image_22.png" /><figcaption>**Figure 21.** The dependence of variations in stellar population properties, $\delta([\text{Fe/H}])$ (top panel), $\delta([\text{O/Fe}])$ (middle panel), and $\delta(\tau)$ (bottom panel) on the azimuthal variations of the radial action computed instantaneously, $\delta(J_R)$ (left panel) and time-averaged, $\delta\langle(J_R)_t\rangle$ (right panel). The Pearson $r$ value of each correlation is indicated in each panel. The population variations are strongly correlated with those of $\langle J_R \rangle$, but very poorly correlated with those of $J_R$.</figcaption></figure>

energy, $E_J \equiv E - \Omega_p J_\phi$, is a conserved quantity (e.g. Binney & Tremaine 2008, chapter 3.3.2), where $E$ is a star's energy and $\Omega_p$ the pattern speed of the perturbation, provided that the system is nearly stationary in the rotating frame of the perturber. Thus, in the $(J_\phi, J_R)$-plane, resonances cause stars to oscillate tangential to the contours of constant $E_J$, even if they are not resonantly trapped.

Furthermore, a star is at a planar resonance when $N_R \Omega_R + N_\phi(\Omega_\phi - \Omega_p) = 0$, where $N_R$ and $N_\phi$ are integers and $\Omega_\phi$ and $\Omega_R$ are

a star's azimuthal and radial frequencies. The principal resonances are the corotation (CR; $N_R = 0$), and the outer and inner Lindblad resonances (OLR and ILR; $N_R = 1$, $N_\phi = \pm 2$). The conservation of $E_J$ leads to the condition

$$ \Delta J_R = \frac{N_R}{N_\phi} \Delta J_\phi \qquad (8) $$
<figure><img src="image_23.png" /><figcaption>**Figure 22.** 〈[Fe/H]〉 profiles at *t* = 13 Gyr for different age populations as indicated. Outside the inner 5 kpc, the profiles become progressively flatter for older populations. The overall 〈[Fe/H]〉 profile (dashed cyan line) is flat. Excluding populations older than 9 Gyr (dashed black line) gives a gradient ∇〈[Fe/H]〉 ≃ 0.03 dex kpc⁻¹.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 23.** The distribution of stars in the Δ*J*<sub>Φ</sub> − *J*<sub>Φ</sub> plane. Prominent ridges are evident, indicative of migration driven by spirals. The horizontal lines delineate the non-migrators, defined as those having |Δ*J*<sub>Φ</sub>| < 100 kpc km s<sup>−1</sup>. The horizontal axis shows *J*<sub>Φ</sub> at 12.65 Gyr while the vertical axis shows *J*<sub>Φ</sub>(13 Gyr) − *J*<sub>Φ</sub>(12.65 Gyr). Only stars younger than 9 Gyr are included.</figcaption></figure>

(e.g. Sellwood & Binney 2002; Chiba, Friske & Schönrich 2021). Therefore, for the ILR, CR, and OLR, the directions of oscillation of resonant orbits have slopes −1/2, 0 and +1/2, respectively, in the (*J*<sub>Φ</sub>, *J*<sub>R</sub>)-plane. We can therefore test whether the variations in *J*<sub>R</sub> and *J*<sub>Φ</sub> follow these conditions by plotting their evolution in the (*J*<sub>Φ</sub>, *J*<sub>R</sub>)-plane.

We carry out this test for a subsample of stars chosen at random. We begin by selecting stars older than 2 Gyr in the radial region $7 \le R/\text{kpc} \le 9$ at 13 Gyr. From the 940 800 stars that satisfy these conditions, we select those with small vertical excursion through $\langle J_z \rangle_t < 5 \text{ kpc km s}^{-1}$ (corresponding to an RMS height of 0.17 kpc), which leaves 441 931 stars. Fig. 26 shows the orbital evolution of a sample of 54 such stars chosen at random in the ($J_\phi$, $J_R$)-plane

<figure><img src="image_25.png" /><figcaption>**Figure 24.** The effect of radial migration on 〈[Fe/H]〉 variations. The variations are shown for a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. The average values for each radius are subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled by the factor of 0.2, for ease of comparison. All profiles exclude stars formed after 12.65 Gyr, for which we cannot compute Δ*J*<sub>Φ</sub>. The red lines show the variations of stellar populations with no cuts on Δ*J*<sub>Φ</sub>, while the blue lines show the variations for the non-migrating populations, i.e. those with |Δ*J*<sub>Φ</sub>| < 100 kpc km s<sup>−1</sup>. Dashed lines show all ages while solid lines show only populations younger than 9 Gyr. Note that the density, shown by black lines, include both a solid line (|Δ*J*<sub>Φ</sub>| < 100 kpc km s<sup>−1</sup>) and a dashed line (no |Δ*J*<sub>Φ</sub>| cut); other than at 6 kpc, these overlap. The dotted horizontal lines indicate the zero for each radius. We show the 13 Gyr distribution.</figcaption></figure>

over the time interval 12 – 13 Gyr.<fnref n="2" /> The age of the selected stars decreases from top to bottom. The three columns show *E*<sub>J</sub> contours for the three main pattern speeds at this time (including the bar’s). It is immediately obvious that the computed actions of all the stars are varying. While some stars do appear to be at resonances [e.g. ILR: (a) cyan and blue; CR: (g) black; OLR: (a) green] for the most part the stars do not appear to be librating about resonances. Nor do they faithfully follow the contours of *E*<sub>J</sub> for any of the main pattern speeds. In all cases, we find considerable jitter in the tracks, which we associate with measurement errors.

We conclude that the variations we find are not predominantly the result of resonances. This does not mean of course that the actions are not changing as a result of resonances, but the actions computed assuming axisymmetry are not able to properly follow this behaviour. While we have shown this for only a tiny fraction of stars, we have verified that this behaviour is generic to other stars. In Appendix C,

<footnote n="2">²We show a random distribution where we only avoid orbits that overlay the resonance arrows in panel (e).</footnote>
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Moving average over $\Delta t = 1$ Gyr intervals, of the phase shift, $\Delta\phi_m$, between the density and $\langle[\text{Fe/H}]\rangle$ waves at different radii, as indicated by the colours, and harmonics (indicated at top left of each panel). $\Delta\phi_2$ appears to fluctuate about 0° throughout the 7 Gyr evolution with $\max(\Delta\phi_2) \le 10°$. Only stars older than 2 Gyr are considered in this analysis.</figcaption></figure>

we present a choice of particles which favours stars trapped at ILRs, and find the same pronounced jitter in the instantaneous actions.

Nor does it seem likely that resonance overlap is driving the changes in the actions, since such scattering leads to heating overall, which we earlier found to be quite modest. We conclude that the variations in $J_R$ are largely the result of errors introduced by the assumption of axisymmetry.

# 9 DISCUSSION

## 9.1 The cause of azimuthal metallicity variations

We have shown that the azimuthal metallicity variations are a result of spiral structure. The peaks and troughs of the $\langle[\text{Fe/H}]\rangle$ waves match those of the density. Their pattern speeds also match those of the spirals. The possibility that the $\langle[\text{Fe/H}]\rangle$ variations are just a result of young, recently formed stars is excluded by the simulation, and by observations, which both show that $\langle[\text{Fe/H}]\rangle$ variations occur also in stars older than 2 Gyr (e.g. Bovy et al. 2014), which will have had sufficient time to move out of their natal spiral arms. Moreover, Roškar et al. (2012) found (their fig. 6) that the lifetimes of spirals is typically a few 100 Myr. We have shown that azimuthal $\langle[\text{Fe/H}]\rangle$ variations are caused by the stronger response to a spiral perturbation of kinematically cooler (therefore younger and metal-rich) stars compared with kinematically hotter (older and metal poorer) stars. As a result, the metal-rich stars trace the density ridge of the spirals. For the same reason, the spirals are traced by minima in [O/Fe] and age (even when young stars are excluded). Our conclusion with a fully self-consistent simulation agrees with that of Khoperskov et al. (2018). We do not find the offsets between the peak metallicity and peak density characteristic of metallicity variations produced primarily by migrating stars (Grand et al. 2016a; Khoperskov et al.

2018, 2023). We find that non-migrating populations also follow the overall trend of the [Fe/H] variation. Because the population that is migrating at any one time is a relatively small fraction of the total (consisting of those currently resonantly trapped at the corotation of any given spiral) this is not an efficient way to generate azimuthal variations, because the majority of the stars do not partake in the migration. Thus, we propose that the azimuthal $\langle[\text{Fe/H}]\rangle$ variations are dominated by the weaker spiral response of populations with increasing radial action and decreasing metallicity. Our conclusion agrees with models 1 and 2 of Khoperskov et al. (2018), where metallicity is a function of both radius (in model 2) and stellar population temperature (both models), but differs from that of model 3, in which, by design, the metallicity is a function only of radius, and independent of $J_R$. As a result, azimuthal variations in model 3 can only be driven by radial migration. Since the age of stars is related to both the metallicity and their dynamical temperature, the conditions of model 3 are unlikely to arise in nature.

The reason for the disagreement with the result of Grand et al. (2016b) is less clear. These authors only considered stars younger than 3 Gyr; it may be that these populations were too similar in $J_R$ to generate much of a different spiral response. Another possible cause for a strong effect of radial migration is that the spiral structure in their model is corotating at all radii (Grand et al. 2012), which may enhance migration-driven azimuthal variations by migrating many stars over very large radial ranges. Alternatively, migration is enhanced by satellite interactions (Grand et al. 2016a; Carr et al. 2022) in their simulation, which are absent in our model.

On the other hand, some important caveats need to be acknowledged about our model, which may lead to underestimating the effect of radial migration. The model may have a steeper AMR in the disc than found in some studies of the MW (e.g. Mackereth et al. 2017; Delgado Mena et al. 2019); this would imply that migration is weaker in the model than in the MW. However, other studies find a steeper AMR in the MW (e.g. Feuillet et al. 2019; Anders et al. 2023; Johnson et al. 2024), consistent with the model. The AMR therefore does not provide an unambiguous test of the model’s migration history.

The model’s metallicity gradient for stars younger than 9 Gyr is lower than the MW’s thin disc’s by a factor of about 2. Since the metallicity gradient depends on both the birth metallicity gradient across time, as well as migration, it also does not test the model’s rate of migration. However, the lower gradient implies that migrating populations generate a weaker azimuthal metallicity variation. Moreover, the $J_R$–[Fe/H] relation in the model is steeper than in the MW, which means that the spirals are bound to produce stronger azimuthal variations because of $J_R$. All these suggest that the relative effect of migration may be underestimated in this model compared with the effect of varying $J_R$. Nevertheless, the fact that only stars trapped at the spirals’ corotation resonances can produce the migration-induced azimuthal variations implies that this is a less efficient mechanism than the $J_R$ variations.

Since, at present, the MW’s spiral structure is not fully known (see, for instance, Shen & Zheng 2020), the azimuthal metallicity variations themselves may therefore be used to probe the spiral structure. Unlike molecular gas clouds, the distance to stars within ~5 kpc can now be determined accurately with Gaia; furthermore, metallicity variations can be mapped away from the dusty mid-plane.

## 9.2 Errors in the radial action

That the axisymmetric approximation needed in computing actions using the Stäckel approximation does not hold has long been obvious. The fact that in the Stäckel approximation itself the errors incurred
<figure><img src="image_27.png" /><figcaption>**Figure 26.** The evolution of a random set of 54 particles, selected as described in the text, in the action space ($J_\phi, J_R$). In each panel, we present six stars with different colours over a 1 Gyr time interval plotting every 50 Myr. The plotted particles are ordered by age, with the oldest in the top row (9.2 $\le \tau$/ Gyr $\le 11.7$), the intermediate-age ones in the middle row (4.7 $\le \tau$/ Gyr $\le 8.8$) and the youngest in the bottom row (2.2 $\le \tau$/ Gyr $\le 4.2$). The dashed yellow lines show contours of constant $E_J$ for $\Omega_p = 15.4 \text{ km s}^{-1} \text{ kpc}^{-1}$ (left panel), 40.0 km s$^{-1}$ kpc$^{-1}$ (centre panel), and 61.5 km s$^{-1}$ kpc$^{-1}$ (right panel). The white lines in panel (e) indicate the libration vectors at the OLR (upper line) and the ILR (lower line). The background shows a log scale representation of the number density of stars.</figcaption></figure>

are less than 10 per cent (Binney 2012) has given rise to the expectation that other sources of error are equally small. Since results based on actions of Solar Neighbourhood disc stars computed using the axisymmetric Stäckel approximation give physically reasonable results (e.g. Monari et al. 2019; Trick, Coronado & Rix 2019; Frankel et al. 2020; Kawata et al. 2021; Drimmel et al. 2023), this has perhaps allayed concerns that the local spiral structure disturbs the computed actions much. We have shown (Fig. 18) that the presence of spiral structure causes errors of the order of $< 100 \text{ kpc km s}^{-1}$ in $J_R$ through most of the disc. The error increases slowly with $J_R$ such that the fractional error is larger at small $J_R$.

More worryingly, $J_R$ computed using the axisymmetric approximation must have *correlated* errors large enough to render them compromised when employed to study spiral structure. In principle, this difficulty can be sufficiently circumvented in simulations by replacing $J_R$ by the time-averaged radial actions, $\langle J_R \rangle_t$, over a time-

scale sufficiently short to limit the effect of resonant heating but longer than the finite lifetime of spirals, to ensure any resonant stars are not trapped for the whole time interval. Naturally, this is not possible in the MW. It may be possible to account for the perturbations from spirals to correct for the radial action errors; the torus-mapping technique (Binney & McMillan 2011; Binney 2016, 2018, 2020) is well suited to this problem (see also Al Kazwini et al. 2022). However, attempting to map the spiral structure using just the radial actions based on an axisymmetric Stäckel approximation gives the wrong result. For instance, Palicio et al. (2023) found that the spirals traced by $J_R$ computed this way did not fully match with any of the current models of the MW’s spiral structure, although it should also be noted that none of the models agree fully with any other.

Despite these errors, we note that the density distribution in Fig. 26, shown by the grey background, still exhibits multiple ridges in the space of instantaneous actions ($J_\phi, J_R$). Similar ridges have been
observed in the MW (Monari et al. 2019; Trick et al. 2019). Thus, the errors are not so large, or the spirals so strong in the simulation, that this important observable is smeared out.

## 9.3 Summary

Our key results are the following:

(i) We have shown that an isolated simulation develops azimuthal $\langle[\text{Fe/H}]\rangle$ variations with peaks at the locations of the density peaks, corresponding to the spiral ridges. These azimuthal variations occur also when stars younger than 2 Gyr are excluded, indicating that high-metallicity star formation is not driving the $\langle[\text{Fe/H}]\rangle$ variations. Indeed, the $\langle[\text{Fe/H}]\rangle$ variations occur in stars of all ages. The peak to trough variations of $\langle[\text{Fe/H}]\rangle$ are of the order of 0.1 dex (see Section 3).

(ii) We bin [Fe/H] on a cylindrical grid and use this binning to compute the pattern speeds of the $\langle[\text{Fe/H}]\rangle$ variations. We compare these pattern speeds to those of the spiral density waves binned in the same way. Because we bin the average [Fe/H], the density is factored out of the $\langle[\text{Fe/H}]\rangle$ pattern speeds. We find that the pattern speeds of $\langle[\text{Fe/H}]\rangle$ match those of the density waves, indicating that the $\langle[\text{Fe/H}]\rangle$ waves are merely a manifestation of the spiral density waves (see Section 4).

(iii) Azimuthal variations occur also in $\langle[\text{O/Fe}]\rangle$ (peak to trough ~0.025 dex) and age, $\langle\tau\rangle$, (peak to trough ~1 Gyr). The variations in $\langle\tau\rangle$ and $\langle[\text{O/Fe}]\rangle$ track density variations slightly less well than does $\langle[\text{Fe/H}]\rangle$, and both are anticorrelated with the density, as would be expected since [Fe/H]-rich stars should be younger and $\alpha$-poor (see Sections 3.2 and 5).

(iv) The radial velocity dispersions, $\sigma_R$, of stars in the simulation increases with age, $\tau$. As a consequence, the spirals are weaker in the old stars than in the younger stars. This variation in spiral strength with $\tau$ is why the [Fe/H]-rich stars are concentrated on the spiral density ridges. However, when we compute the radial action, $J_R$, and compare the azimuthal variation of $\langle J_R \rangle$, we find that they generally correlate with those of the density, the opposite of what is expected (see Sections 5 and 6.1).

(v) We demonstrate that the axisymmetry approximation is introducing correlated errors in the computation of $J_R$. In order to explore this further, we measure $J_R$ over 1 Gyr, and compute the time-averaged radial action, $\langle J_R \rangle_i$. We note both a heating of the stellar populations, and a comparable (but larger) scatter in $J_R$. We show that the azimuthal variations of $\langle J_R \rangle_i$ are anticorrelated with the variations of the density, as expected if the spirals are stronger in cooler, younger, metal-rich populations. We further show that the variations of metallicity, age and $\alpha$-abundance are poorly correlated with those of $\langle J_R \rangle_i$ but very well correlated with those of $\langle J_R \rangle_i$ (see Section 6).

(vi) Thus, the axisymmetric Stäckel approximation overestimates the radial action in high-density regions and underestimates it in low-density regions. For studies that depend critically on the behaviour of one $J_R$ population versus that of another, it is important to correct for this effect. We have done this by time-averaging $J_R$ over 1 Gyr, during which time star particles will have drifted in and out of spiral arms. Such an approach is of course not possible in the MW; other ways of correcting for these errors are needed (see Section 6).

## ACKNOWLEDGEMENTS

We thank the anonymous referee whose comments helped improve this paper, and Pedro Palicio for useful discussions. TK acknowl-

edges support from the National Science Foundation of China (NSFC) under grant no. 12303013 and support from the China Postdoctoral Science Foundation under grant no. 2023M732250. JASA and CFPL acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 852839). JA also acknowledges support by the National Natural Science Foundation of China under grant no. 12122301, 1223001, by the '111' project of the Ministry of Education under grant no. B20019, and the sponsorship from Yangyang Development Fund. L.BeS. acknowledges the support provided by the Heising Simons Foundation through the Barbara Pichardo Future Faculty Fellowship from grant no 2022-3927, and of the NASA-ATP award 80NSSC20K0509 and U.S. National Science Foundation AAG grant AST-2009122. The simulation was run at the DiRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/J005673/1, STFC capital grant ST/H008586/1 and STFC DiRAC Operations grant vST/K00333X/1. DiRAC is part of the National E-Infrastructure. The analysis used the PYTHON library PYNBODY (Pontzen et al. 2013). The pattern speed analysis was performed on Stardynamics, which was funded through Newton Advanced Fellowship NA150272 awarded by the Royal Society and the Newton Fund.

## DATA AVAILABILITY

The simulation snapshot in this paper can be shared on reasonable request.

## REFERENCES

* Al Kazwini H. et al., 2022, A&A, 658, A50

* Anders F. et al., 2023, A&A, 678, A158

* Beane A., Ness M. K., Bedell M., 2018, ApJ, 867, 31

* Bellardini M. A., Wetzel A., Loebman S. R., Faucher-Giguère C.-A., Ma X., Feldmann R., 2021, MNRAS, 505, 4586

* Bellardini M. A., Wetzel A., Loebman S. R., Bailin J., 2022, MNRAS, 514, 4270

* Bergemann M. et al., 2014, A&A, 565, A89

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

* Binney J., 2016, MNRAS, 462, 2792

* Binney J., 2018, MNRAS, 474, 2706

* Binney J., 2020, MNRAS, 495, 886

* Binney J., McMillan P., 2011, MNRAS, 413, 1889

* Binney J., Tremaine S., 2008, in Binney J., Tremaine S., eds, Galactic Dynamics: 2nd edn. Princeton Univ. Press, Princeton, NJ

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

* Bovy J. et al., 2014, ApJ, 790, 127

* Carr C., Johnston K. V., Laporte C. F. P., Ness M. K., 2022, MNRAS, 516, 5067

* Casagrande L., Schönrich R., Asplund M., Cassisi S., Ramírez I., Meléndez J., Bensby T., Feltzing S., 2011, A&A, 530, A138

* Chiba R., Friske J. K. S., Schönrich R., 2021, MNRAS, 500, 4710

* Delgado Mena E. et al., 2019, A&A, 624, A78

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

* Drimmel R. et al., 2023, A&A, 670, A10

* Feuillet D. K., Frankel N., Lind K., Frinchaboy P. M., García-Hernández D. A., Lane R. R., Nitschelm C., Roman-Lopes A., 2019, MNRAS, 489, 1742

* Fiteni K., Caruana J., Amarante J. A. S., Debattista V. P., Beraldo e Silva L., 2021, MNRAS, 503, 1418
Frankel N., Sanders J., Ting Y.-S., Rix H.-W., 2020, ApJ, 896, 15

Gaia Collaboration, 2018, A&A, 616, A11

Gandhi S. S., Ness M. K., 2019, ApJ, 880, 134

Ghosh S., Debattista V. P., Khachaturyants T., 2022, MNRAS, 511, 784

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

Grand R. J. J., Springel V., Gómez F. A., Marinacci F., Pakmor R., Campbell
D. J. R., Jenkins A., 2016a, MNRAS, 459, 199

Grand R. J. J. et al., 2016b, MNRAS, 460, L94

Hackshaw Z., Hawkins K., Filion C., Horta D., Laporte C. F. P., Carr C.,
Price-Whelan A. M., 2024, ApJ, 977, 15

Hawkins K., 2023, MNRAS, 525, 3318

Ho I. T. et al., 2018, A&A, 618, A64

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Hwang H.-C. et al., 2019, ApJ, 872, 144

Imig J. et al., 2023, ApJ, 954, 124

Johnson J. W. et al., 2024, preprint (arXiv:2410.13256)

Kawata D., Baba J., Hunt J. A. S., Schönrich R., Ciucă I., Friske J., Seabroke
G., Cropper M., 2021, MNRAS, 508, 728

Khachaturyants T., Beraldo e Silva L., Debattista V. P., Daniel K. J., 2022a,
MNRAS, 512, 3500

Khachaturyants T., Debattista V. P., Ghosh S., Beraldo e Silva L., Daniel K.
J., 2022b, MNRAS, 517, L55

Khoperskov S., Di Matteo P., Haywood M., Combes F., 2018, A&A,
611, L2

Khoperskov S., Sivkova E., Saburova A., Vasiliev E., Shustov B., Minchev
I., Walcher C. J., 2023, A&A, 671, A56

Kreckel K. et al., 2019, ApJ, 887, 80

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

Lemasle B., François P., Piersimoni A., Pedicelli S., Bono G., Laney C. D.,
Primas F., Romaniello M., 2008, A&A, 490, 613

Lépine J. R. D. et al., 2011, MNRAS, 417, 698

Luck R. E., Kovtyukh V. V., Andrievsky S. M., 2006, AJ, 132, 902

Luck R. E., Andrievsky S. M., Kovtyukh V. V., Gieren W., Graczyk D., 2011,
AJ, 142, 51

Mackereth J. T. et al., 2017, MNRAS, 471, 3057

Majewski S. R. et al., 2016, Astron. Nachr., 337, 863

Monari G., Famaey B., Siebert A., Wegg C., Gerhard O., 2019, A&A, 626,
A41

Nordström B. et al., 2004, A&A, 418, 989

Palicio P. A., Recio-Blanco A., Poggio E., Antoja T., McMillan P. J., Spitoni
E., 2023, A&A, 670, L7

Pedicelli S. et al., 2009, A&A, 504, 81

Peebles P. J. E., 1969, ApJ, 155, 393

Poggio E. et al., 2022, A&A, 666, L4

Pontzen A., Roškar R., Stinson G. S., Woods R., Reed D. M., Coles J., Quinn
T. R., 2013, synbody: Astrophysics Simulation Analysis for Python

Roškar R., Debattista V. P., Quinn T. R., Wadsley J., 2012, MNRAS, 426,
2089

Sánchez-Menguiano L. et al., 2016, ApJ, 830, L40

Schaye J. et al., 2015, MNRAS, 446, 521

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

Shen J., Zheng X.-W., 2020, Res. Astron. Astrophys., 20, 159

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Solar M., Tissera P. B., Hernandez-Jimenez J. A., 2020, MNRAS, 491, 4894

Spitoni E., Cescutti G., Minchev I., Matteucci F., Silva Aguirre V., Martig
M., Bono G., Chiappini C., 2019, A&A, 628, A38

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006,
MNRAS, 373, 1074

Trick W. H., Coronado J., Rix H.-W., 2019, MNRAS, 484, 3291

Vasiliev E., 2019, MNRAS, 482, 1525

Vickers J. J., Shen J., Li Z.-Y., 2021, ApJ, 922, 189

Vogt F. P. A., Pérez E., Dopita M. A., Verdes-Montenegro L., Borthakur S.,
2017, A&A, 601, A61

Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137

Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS, 471, 2357

Wenger T. V., Balser D. S., Anderson L. D., Bania T. M., 2019, ApJ, 887, 114

Willett E. et al., 2023, MNRAS, 526, 2141

APPENDIX A: OTHER TIMESTEPS

We present ⟨[Fe/H]⟩ maps (Fig. A1) and azimuthal profiles of Σ*
and ⟨[Fe/H]⟩ (Fig. A2) of the model at *t* = 8 to 13 Gyr. The general
behaviour described in the main text holds at earlier times.
<figure><img src="image_28.png" /><figcaption>**Figure A1.** ([Fe/H]) in cylindrical coordinates at *t* = 8–13 Gyr, as labelled. Only stars older than 2 Gyr are considered (compare Fig. 1). The contours correspond to the full surface density. Note the changing colour scale.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure A2. Same as the left panel of Fig. 2 at t = 8–13 Gyr, as labelled.</figcaption></figure>
## APPENDIX B: MEDIAN [Fe/H] AND RADIAL ACTION

Since distributions of [Fe/H] and $J_R$ are generally skewed, we demonstrate that the results in the main text would still hold using the medians instead of the means. Fig. B1 shows maps of medians, while Fig. B2 shows azimuthal profiles of medians. Despite the ~0.1 dex difference in [Fe/H] between the mean and the median, the trends are comparable. A similar result holds for $J_R$.

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

**Figure B1.** Top panel: median [Fe/H]. Compare with the middle panel of Fig. 1. Middle and bottom panels: median $J_R$ and $\langle J_R \rangle_t$. Compare with Fig. 14. In all panels, the contours correspond to the surface density. The disc is rotating in a counter-clockwise sense, i.e. rotation is in the direction of increasing $\phi$. The model is shown at 13 Gyr.
APPENDIX C: OTHER PARTICLES IN ACTION
SPACE

We consider a different selection of particles to track in the ($J_\phi$, $J_R$)-space. We now add the additional selection that the variation of $J_R$ is large by requiring that its dispersion for values computed over 1 Gyr is large: $25 < \sigma_{J_R}/\text{kpc km s}^{-1} < 50$. The results are shown in Fig. C1. The new selection favours stars trapped at the Lindblad resonances and we see that many stars appear to be trapped at the ILR of the slowest spiral. Nevertheless even for these particles many of the tracks show considerable jitter, as was the case in Fig. 26, reinforcing the conclusion that measuring actions assuming axisymmetry introduces significant errors.

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

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

**Figure B2.** The variation of median [Fe/H] (top panel) and −median $\langle J_R \rangle_t$ (bottom panel) compared with the density variation at a series of 500 pc-wide annuli centred at 6 to 10 kpc, as indicated. The average value at each radius is subtracted from each azimuthal profile and the profiles are then vertically offset by a fixed amount to show the variation. The azimuthal profiles of the density have been scaled for ease of comparison. The dotted horizontal lines indicate the zero for each radius. In the top panel, the red lines show profiles of median [Fe/H] for all stars, while the green lines show stars older than 2 Gyr only (both panels). Compare these with Fig. 2 and the right panel of Fig. 15.
<figure><img src="image_33.png" /></figure>

**Figure C1.** The evolution of a random set of 54 particles in the action space ($J_\phi, J_R$). These particles have been selected as in Fig. C1 with the extra condition that the spread in $J_R$ of each particle is large. In each panel, we present six stars with different colours over a 1 Gyr time interval, plotting every 50 Myr. The plotted particles are ordered by age, with the oldest in the top row ($8.9 \le \tau$/ Gyr $\le 12.7$), the intermediate-age ones in the middle row ($5.8 \le \tau$/ Gyr $\le 8.9$) and the youngest in the bottom row ($2.4 \le \tau$/ Gyr $\le 5.7$). The dashed yellow lines show contours of constant $E_J$ for $\Omega_p = 15.4 \text{ km s}^{-1} \text{ kpc}^{-1}$ (left panel), $40.0 \text{ km s}^{-1} \text{ kpc}^{-1}$ (centre panel), and $61.5 \text{ km s}^{-1} \text{ kpc}^{-1}$ (right panel). The white lines in panel (e) indicate the libration vectors at the OLR (upper line) and the ILR (lower line). The background shows a log scale representation of the number density of stars.

This paper has been typeset from a TEX/LATEX file prepared by the author.