Radial phase spirals in the Solar neighbourhood

The second data release of ESA's $Gaia$ mission revealed numerous signatures of disequilibrium in the Milky Way's disc. These signatures are seen in the planar kinematics of stars, which manifest as ridges and ripples in $R-v_{\phi}$, and in vertical kinematics, where a prominent spiral is seen in the $z-v_z$ phase space. In this work, we show an equivalent $\Delta R-v_{\mathrm{R}}$ phase spiral forms following a perturbation to the disc. We demonstrate the behaviour of the $\Delta R-v_{\mathrm{R}}$ phase spirals in both a toy model and a high resolution $N$-body simulation of a satellite interaction. We then confront these models with the data, where we find partial $\Delta R-v_{\mathrm{R}}$ phase spirals in the Solar neighborhood using the most recent data from $Gaia$ DR3. This structure indicates ongoing radial phase mixing in the Galactic disc, suggesting a history of recent perturbations, either through internal or external (e.g., satellite) processes. Future work modelling the $z-v_z$ and $\Delta R-v_{\mathrm{R}}$ phase spirals in tandem may help break degeneracy's between possible origins of the perturbation.


INTRODUCTION
The European Space Agency's Gaia mission (Gaia Collaboration et al. 2016) has revolutionised studies of stellar dynamics in the Milky Way by measuring the position and motion of ∼ 1.8 × 10 9 stars in our Galaxy.In turn, this detailed catalogue has revealed the degree of disequilibrium present in the Galactic disc, motivating a need to develop modelling and analysis techniques capable of quantifying and interpreting disequilibrium dynamics.
One of the most striking disequilibrium features is the z − vz phase spiral discovered by Antoja et al. (2018) in Gaia DR2 (Gaia Collaboration et al. 2018), which showed that the local Milky Way disc is phase mixing vertically following some perturbation.Later studies have shown such phase spirals are present over several kpc in the Galactic disc (e.g.Hunt et al. 2022;Antoja et al. 2023;Frankel et al. 2023;Alinder et al. 2023), and simulations have shown that we expect them to occur across the entire disc (e.g.Bland-⋆ E-mail: j.a.hunt@surrey.ac.uk Hawthorn & Tepper-García 2021; Hunt et al. 2021).There are currently many proposed explanations for the origins of the z − vz phase spirals: The Sagittarius dwarf galaxy, (e.g.Antoja et al. 2018;Laporte et al. 2018a;Hunt et al. 2021), bar buckling (Khoperskov et al. 2019), torque from a triaxial halo (Grand et al. 2023), stochastic perturbation from a dark matter subhalo population (Tremaine et al. 2023), or transient spiral structure (Hunt et al. 2022).
Regardless of whether the z − vz phase spirals are instigated by one of the above processes or a complex combination of them (e.g.García-Conde et al. 2023), the same perturbing influence is likely to impart a force to the radial and azimuthal components of a star's orbit.There has long been evidence that the Milky Way is 'ringing' following some perturbation (Minchev et al. 2009), and Gaia has now revealed large scale departures from equilibrium in radial and azimuthal stellar kinematics such as the ridges in R − v ϕ (e.g.Antoja et al. 2018;Kawata et al. 2018), which can also be created by the Galactic bar (e.g.Fragkoudi et al. 2019), spiral arms (e.g.Hunt et al. 2018Hunt et al. , 2019a) ) or interac-tion with a satellite (e.g.Laporte et al. 2018a;Khanna et al. 2019;Hunt et al. 2021).
What has currently not been explored is whether a phase spiral also exists in the radial position and motion of stars in the disc following some perturbation, analogous to the z − vz spiral in the vertical motion.Belokurov et al. (2023) recently discovered such radial phase mixing is ongoing in the accreted part of the stellar halo.They find phase mixing 'chevrons' in r − vr which are likely tidal debris of the ancient merger event known as Gaia-Sausage/Enceladus (Belokurov et al. 2018;Helmi et al. 2018).While Belokurov et al. (2023) work in spherical coordinates, r and vr, as is appropriate for the stellar halo, similar structures should also be present in the cylindrical radius and radial motion R and vR of stars in the disc.We would expect a radial phase spiral to be a natural consequence of radial phase mixing following some more recent perturbations to the in-situ disc stars, e.g. from the Sagittarius dwarf galaxy, or the Milky Way's bar and spiral arms.
The R − vR plane for disc stars has been explored previously (e.g.Khanna et al. 2019, who find 'cone like' structure as an indication of phase mixing), but owing to high dust extinction in the Galactic mid-plane, it is challenging to resolve such structure further than ∼ 2 kpc, and the distribution is dominated by stars close to the Sun.However, the equivalent plane to z − vz is not R − vR, but rather ∆R − vR, where ∆R = R−RG, e.g.capturing the amplitude and phase of the radial epicyclic oscillation around the Guiding radius, RG. ∆R − vR has the advantage that stars near the Sun have a range of radial eccentricities and phases, and we can probe a range of ∆R even with a local sample.
In this work we reveal ∆R − vR phase spirals in simulations and in the Galactic disc in Gaia DR3 (Gaia Collaboration et al. 2023).In Section 2 we illustrate the shape of orbits in the ∆R − vR plane compared to the z − vz plane.In Section 3 we describe a toy model (3.1) and the N -body simulation used in this paper (3.2), and in Section 3.2.1 we show the spirals and substructure in the ∆R − vR plane both locally and across the disc of the simulation.In Section 4 we show the ∆R − vR plane data from Gaia DR3.We describe our treatment of the Gaia data in Section 4.1, and in Section 4.2 we show the ∆R−vR plane in the Solar neighborhood and over 2 kpc.In Section 4.2.1 we decompose the 2 kpc sample as a function of angular momentum, and the describe the introduced biases in Section 4.2.2.In Section 4.2.3 we discuss what we can learn from this projection, and in Section 5 we summarise our conclusions.

ORBITS AND EPICYCLES IN AN AXISYMMETRIC POTENTIAL
In this section we introduce the ∆R − vR phase plane and illustrate how individual disc stars evolve along their orbits in the absence of any external perturbations.Stars on sufficiently circular orbits move along nearly closed curves in the z − vz plane as they oscillate about the Galactic midplane.Similarly, stars evolve along nearly closed curves in ∆R − vR as they progress along their radial epicycles centered on their guiding radius, RG.In both cases, as the orbits become less circular (higher action), these tra- Four example orbits, each with the same J ϕ and J R = 0 and increasing Jz (dark to light).Right: Six orbits with identical J ϕ and Jz = 0, but increasing J R (dark to light).Note that in both cases for orbits with higher action, the trajectories remain nearly closed but increasingly diverge from ellipses in a different way.
jectories in phase space become less like simple harmonic motion ellipses and distort into more complex shapes.
Figure 1 shows example orbits calculated in the MilkyWayPotential2022 potential from gala (Price-Whelan 2017).Each orbit has the same J ϕ =1900 km s −1 kpc −1 .The left panel shows orbits with zero radial action, JR = 0, and increasing vertical action (dark to light) in the z − vz plane.The right panel shows orbits with zero vertical action, Jz = 0, and increasing radial action (dark to light) in the ∆R − vR plane.
Note the difference in the shape of the distortion for orbits in the ∆R − vR and z − vz planes.Since the effective potential is symmetric about the midplane, the z − vz trajectories remain symmetric even as they diverge from simple ellipses.However, since the effective potential decreases with increasing galactic radius, it is not symmetric about RG and consequently the expected trajectories are themselves not symmetric.Instead, stars in the ∆R − vR plane move along "arrowhead"-like trajectories, reaching a larger ∆R when further from the galactic center than when nearer it.
Regardless, in both cases, when considering a phase mixed system of many stars, these shapes can be traced by any label invariant to orbital phase, such as the local stellar density (assuming no selection effects) or stellar abundances (Price-Whelan et al. 2021;Horta et al. 2023).

SIMULATED DATA
In this section we explore the creation of ∆R − vR phase spirals in a toy model, and a high resolution N -body simulation to aid with later comparison and interpretation of the Gaia data.

Toy model
As an initial illustration we present a toy model of particles evolved in a fixed Milky Way-like potential following an artificial kick.

Phase mixing following a kick
We apply three different impulsive 'kicks' by offsetting the radial velocities from the equilibrium distribution by 25, 50 and 75 km s −1 , and then integrate the sample within the MilkyWayPotential2022 for 1 Gyr.
Figure 2 shows the initial distribution of the sample in ∆R − vR (left column), the sample after 500 Myr (second column) and after 1 Gyr (third column), for the sample with the 25 km s −1 (upper row), 50 km s −1 (second row), and 75 km s −1 (third row) kicks.The range of the left three columns is chosen for later comparison with the data, and the right hand column shows the full distribution for completeness.
The lower row of Figure 2 shows the same sample kicked by 20 km s −1 in the vertical direction, for comparison.
Figure 2 shows that phase spirals do form in ∆R − vR as the samples phase mix back towards equilibrium following the 'kick', as expected.The trend of increasingly strong ∆R−vR phase spirals with increasingly strong kicks is clear, and unsurprising.While the amplitude of the spiral changes, the general morphology does not, as it depends on the shape of the galactic potential and the direction of the 'kick', which remain fixed between samples.The ∆R − vR phase spirals in the stronger cases are quite striking, but this toy example shows only a limited sample without selection effects, and without self-gravity, which is known to be important in phase spiral formation and morphology (Widrow 2023).
Note that the distribution in ∆R − vR is asymmetric in the radial dimension in both the initial sample, (which matches the shape of the orbits in the right panel of Figure 1) and following the perturbation, owing to the radially decreasing nature of the galactic potential.This is in contrast to the symmetric nature of the initial z − vz plane shown in the lower row (although the kick and subsequent phase mixing makes the z − vz distribution asymmetric until it phase mixes), owing to the symmetric nature of the vertical potential.This illustrates that while phase mixing and thus phase spirals occur both radially and vertically, we do not expect them to produce phase spirals with the same morphology.
Also, in this simple toy model a vertical kick has almost no discernible effect on the morphology of the ∆R−vR spiral, and similarly a radial kick has almost no discernible effect on the morphology of the z − vz spiral.In reality we would expect a perturber such as a satellite galaxy to impart a radial and vertical kick together, and we would expect there to be some coupling in the response in a self-gravitating disc.However, this simple toy is enough for now to illustrate the expected shape and properties of a R − vR phase spiral, and we defer a more thorough exploration of coupling to a future work.

The N -body simulation
To explore ∆R−vR phase spiral formation in a more realistic model, we make use of the N -body simulation M1 from Hunt et al. (2021), which consists of a merger of a dwarf galaxy into a highly stable disc.Full details of the simulation setup and evolution are available in Hunt et al. (2021), but in brief; the host galaxy is based on the MWb model from Widrow & Dubinski (2005) with total mass ∼ 6 × 10 11 M⊙, consisting of a 8.8 × 10 8 particle live NFW halo (Navarro et al. 1997), a 2.2 × 10 7 particle Hernquist bulge (Hernquist 1990) and a 2.2 × 10 8 particle exponential disc with a high Toomre parameter Q = σRκ/3.36GΣ= 2.2 (Toomre 1964) which is stable against bar and spiral formation for several Gyr without the addition of the satellite perturber (Hunt et al. 2021).The satellite is the L2 model from Laporte et al. (2018b), which is comprised of two Hernquist spheres (Hernquist 1990), with a virial mass of ∼ 6 × 10 10 M⊙, making the merger approximately 1-10 mass ratio.
The combined model is evolved with the N -body tree code Bonsai (Bédorf et al. 2012(Bédorf et al. , 2014) ) for 8.292 Gyr, during which period the satellite merges into the host galaxy.In this instance we choose to work with both an 'early time' snapshot 454 at t = 4.45 Gyr, which has experi-enced a single pericentric passage of the satellite 2 Gyr previously, and a 'late time' snapshot 720, at t = 7.05 Gyr owing to the presence of both a tightly wound z − vz phase spiral such as seen in the Gaia data, and the presence of a ∆R − vR phase spiral in the 'Solar neighborhood' when assuming the observer to be located at (x, y, z)⊙ = (−8.2,0, 0.02) kpc.In the 'late time' snapshot the dwarf galaxy is located at (x, y, z, vx, vy, vz) = ( −8.1, 1.2, 24.8, −117.4, 8.6, −30.2) kpc, km s −1 compared to the estimated position and motion of the Sagittarius remnant (x, y, z, vx, vy, vz) = (9.5, −0.4,−6.7, 115.7, 5.9, 311.9) kpc, km s −1 (Vasiliev & Belokurov 2020).Thus, the dwarf galaxy in the Simulation is not at the same phase as Sagittarius, and we are not attempting to reproduce the conditions in the Milky Way, merely illustrate the phase spirals following multiple perturbations.
In general, we stress that this model is not designed to specifically reproduce the Milky Way Sagittarius interaction (for which we refer the reader to the models presented in Bennett et al. 2022;Stelea et al. 2023), but rather be a general laboratory for investigating merger dynamics in an otherwise quiet disc.The satellite in model M1 has experienced seven disc crossings by the 'late time' snapshot and is about to cross again in ∼ 150 Myr (see Hunt et al. 2021, for the satellite orbit and mass loss history).While the exact number of disc crossing experienced by the Sagittarius dwarf is not fully known (although it can be constrained by changes in star formation and disk velocity dispersion; e.g.Ruiz-Lara et al. 2020;Das et al. 2024) it has certainly experienced multiple passages since first infall.
To calculate action-angle variables for the simulation we first use Agama (Vasiliev 2019) to reconstruct the simulation potential as a composite of two multipole expansions for the stellar bulge and dark matter halo, and an axisymmetric CylSpline expansion for the disc.We then use Agama's ActionFinder routine to calculate actions, angles and frequencies in the reconstructed potential.(lower left) for disc stars in the 'late time' snapshot within a local volume of 2 kpc selected around (x, y, z) = (7.2,0, 0) kpc, and with Guiding radii 7 < RG < 7.4, i.e. within 0.2 kpc of the 'local' radius.This is an arbitrary choice in the simulation, chosen merely as an illustration of a volume which contains two clear phase spirals.
For visualisation purposes we have subtracted a Gaussian smoothed background distribution, such that the left column shows the difference in density between the data and the smoothed version.The smoothing, is different in z − vz (20% of axis range) and ∆R − vR (30% of axis range), although all subsequent ∆R − vR panels have the same smoothing.The z − vz phase spiral in the top left panel is qualitatively similar to the local Gaia data, it has been reproduced in many simulations and has been discussed in numerous other works (e.g.Antoja et al. 2018;Laporte et al. 2018a;Hunt et al. 2021;Bland-Hawthorn & Tepper-García 2021).
Any perturbation to the disc should phase mix away vertically, azimuthally and radially, and thus such a ∆R − vR phase spiral should be expected, as shown in the toy model in Section 3.1.However, the z − vz phase spiral is especially striking because the vertical frequency of stars, Ωz, is a strong function of vertical action, Jz, as shown in the upper right panel of Figure 3 (where Jz is approximately the distance from the centre of the panel).This causes a perturbation to rapidly 'wind up' into a spiral pattern.The radial frequency of stars, ΩR is a much weaker function of the radial action, JR, as shown in the lower right panel of Figure 3.
To further illustrate this, the upper row of Figure 4 shows the vertical (left) and radial (right) action as a function of angular momentum, J ϕ , colored by the vertical and radial frequencies Ωz and ΩR, respectively, for disc stars in the simulation M1.The upper left panel shows that Ωz is a strong function of both J ϕ and Jz.In contrast, the upper right panel shows that while ΩR is also a strong function of J ϕ , the gradient of ΩR with respect to JR is much shallower.However, while shallow, there is a gradient, which is required for phase spiral formation.
The lower row of Figure 4 shows the difference between the vertical (left) and radial (right) frequencies of the disc stars and the asymptotic frequency, ν = Ωz|Jz = 0 (km s −1 kpc −1 ) and κ = ΩR|JR = 0 (km s −1 kpc −1 ), respectively, for a given J ϕ .The lower left panel shows that the spread between the asymptotic vertical frequency, ν, as a function of action is much larger than the spread between the asymptotic radial frequency, κ, at a given J ϕ .
Thus, the timescale for mixing in this space is much longer than the vertical phase mixing in the z − vz plane, resulting in a much more open pitch angle of the spiral in ∆R − vR compared to z − vz(as also shown when comparing the toy models for ∆R − vR and z − vz in Figure 2).Given the long time scale we would expect the presence of bar resonances and spiral arms to readily interfere with older patterns, but perhaps we can still exploit any structure in ∆R − vR much as we have done with the phase spirals in z − vz.
In order to test whether such ∆R − vR phase spirals are common within the simulation we relax the local selection.Following Hunt et al. (2021), Figure 5 shows the ∆R − vR distribution in 30x30 1 kpc bins across the face of the disc for the late time snapshot at 7.05 Gyr. Figure 5 shows a large variation in the morphology of the individual ∆R − vR planes.Some show spiral-like patterns (most commonly around a radius of ∼ 7 kpc), while others are chaotic, or simply arcs (more commonly in the outer disc around ∼ 13 kpc).So, 'clean' ∆R − vR phase spirals such as are shown in Figure 3 are not 'common' within the simulation snapshot, but neither are they particularly rare, or unique.When examining the time evolution of the local R − vR plane3 , the Solar neighborhood does not always contain ∆R−vR phase spirals.Some snapshots contain complex structure in the ∆R − vR plane but no phase spirals, as we would expect if they are disrupted by secular dynamics such as the bar resonances or the passage of a spiral arm.
In Hunt et al. (2021) we showed that the z−vz phase spirals became increasingly messy following repeated satellite impacts, while earlier snapshots showed well defined spirals after the first impact alone (see Figures 7 & 9 in Hunt et al. 2021).As such, Figure 6 shows the same plane of ∆R − vR distributions at an earlier time (t = 4.49 Gyr).This snapshot is 2 Gyr after the first impact and shortly before the second interaction.Around the 'Solar radius' (R ∼ 8 kpc) many weak partial phase spirals have formed, while the outer disc is dominated by arcs of varying completeness, resembling "pac-men", or "fortune cookies", instead of spirals.As expected, the shallow gradient of ΩR as a function of JR (shown in Figure 3) has resulted in only partial phase mixing (and minimal spiral formation) especially in the outer galaxy where stellar frequencies are lower, even over 2 Gyr.
While the ∆R − vR planes do not show 'clean' phase spirals such as are seen in the z − vz direction (see Figures 7 and 9 of Hunt et al. 2021, or the animated evolution of both side by side here 3 ), it remains clear that the passage of the satellite galaxy imparts coherent substructure in the radial oscillation of stars as well as the vertical, as expected.Note that the individual ∆R − vR planes trace out galaxy spanning spiral patterns, similar to the z − vz planes shown in Hunt et al. (2021).
Most importantly, these ∆R − vR features can perhaps be exploited in a similar fashion as is being done in studies of the z − vz spirals to infer information about the Galactic potential, or the nature of the perturbation (e.g.Widmark et al. 2022a;Frankel et al. 2023;Darragh-Ford et al. 2023;Guo et al. 2023) if such structures are also seen in the data.

GAIA DR3 DATA
In this section we describe our treatment of the Gaia DR3 data, and the ∆R − vR substructure in the Solar neighborhood and beyond.

∆R − vR structures in the Gaia data
In this section we show the dynamical structure present in ∆R − vR in the Solar neighborhood, and as a function of angular momentum.
As an initial illustration, the left panel of Figure 7 shows the local vR − v ϕ plane for all stars within 200 pc.The classical moving groups originally discovered in a series of papers by Olin Eggen (see Eggen 1996, and references within for a summary) are clearly visible, and have been marked.The right panel shows the ∆R − vR plane (axes swapped to match vR) for the same 200 pc sample, in which the structure is clearly the same.This is to be expected within a very local volume, because we see only a small section of these dynamical structures even though we know they span several kpc across the Galactic disc (e.g.Kawata et al. 2018;Gaia Collaboration et al. 2021).The ∆R − vR panel is effectively an orbit map, yet within a very local volume stars on the same orbit must have the same kinematics, and very little is learned beyond the classic vR − v ϕ plane.For example, at a single radius there is a clear relation between v ϕ and ∆R.Only a single v ϕ will lead to a circular orbit (∆R = 0), and a star with higher or lower v ϕ will have ∆R negative or positive respectively, while on the inner or outer part of their radial epicycles.However, when analysing a larger volume, orbit maps become more informative than kinematic maps (see e.g.Hunt et al. 2020).For example, when exploring a large range of R, there are many v ϕ which can correspond to ∆R = 0 for stars at different radii (assuming a non-flat rotation curve), and v ϕ is no longer a direct map to ∆R (for any radial eccentricity).The effect over a small radial range is minimal, but as we expand our sample to increasingly large distances ∆R will increasingly diverge from v ϕ .
As such, the left panel of Figure 8 shows the ∆R − vR plane for stars within 2 kpc from the Sun.We now recover arcs with higher | ∆R |, but still we do not find a clean phase spiral, but more of a 'rose-like' structure.In addition, the imprint of the well known local moving groups; Hyades, Pleiades, Hercules, etc. remain visible with comparatively low radial eccentricity in the centre of the panel, and have been marked again for comparison with Figure 7 (note that many of the unmarked features are also previously discovered kinematic groups; see e.g.Ramos et al. 2018, but a detailed census of kinematic substructure is not the purpose of this work).To illustrate the utility of the ∆R − vR plane over the R − vR plane, the right panel of Figure 8 shows the R − vR plane for the same sample of stars.Here, some of the substructure is faintly visible as chevrons (as shown in r − vr for the halo in Belokurov et al. 2023) but it is not possible to resolve any full wraps within a local sample.
While the left panel is not a perfect phase spiral either, we do also expect the ∆R−vR plane to be a superposition of many different ∆R − vR structures as a function of angular momentum, Lz (or RG), and θ ϕ as has been shown in several works for the z − vz phase spiral (Li 2021;Hunt et al. 2021;Gandhi et al. 2022;Darragh-Ford et al. 2023).

As a function of angular momentum, J ϕ
In order to test how the structures vary with angular momentum, Figure 9 shows the ∆R − vR plane for stars with d < 2 kpc, in nine guiding radius bins from 6.2 < RG < 10.2 kpc, each selected such that | RG − RC |< 0.5 kpc (where RC is the centre of the bin).We recover several fragments of arcs and ridges, but no complete phase spirals, owing to the selection function of the sample.
From these partial fragments in phase and angular momentum it's difficult to conclusively state whether the complete ∆R − vR planes would form phase spirals, or whether they instead consist solely of arcs and ridges of stars following the same orbital tracks around their radial epicycle.Both cases are present in the late time simulation snapshot presented in Figure 5, which contains many ∆R − vR planes that are qualitatively similar to the structure in Figure 9 (although the features are finer in the data).Conversely, the early time simulation snapshot in Figure 6 predominantly presents single arcs of varying lengths in the outer disc or faint partially wound spirals in the inner disc, neither of which are a good match to the data.This is as unsurprising, because the early time snapshot in Figure 6 has only experienced one single satellite impact approximately 2 Gyr previously, and the 'hotness' of the disc has prevented additional secular perturbations such as from a bar or internally excited spiral arms.In contrast, the late time snapshot in Figure 5 has experienced several satellite impacts and the disc contains significant non-axisymmetric structure (even if the disc remains hot).The Milky Way has experienced multiple pericentric passages from the Sagittarius dwarf galaxy, and contains a central bar and spiral arms.Thus, we would expect the response in the ∆R − vR plane in the data to closer resemble the late time snapshot, which it does.

Selection function and biases
Note that Figure 9 also has an implicit bias in radial action, JR, and conjugate radial angle, θR, because stars can either be close to their guiding radius by having low radial action (and are 'always' close) while stars with high JR will only appear if they are at the appropriate part of their radial epicycle.The physical selection limits mean that for stars with low Lz we only see them close to apocentre (e.g. the top left panel), and for stars with high Lz we only see them close to their pericentre (e.g. the lower right panel).Thus across the range of guiding radius bins we probe all radial phases, but for different stars with different JR and θR and different ∆R − vR planes.
To examine the selection effects and JR, θR biases, we first apply the same selection (d < 2 kpc, and nine bins in Guiding radius, where | RG −RC |< 0.5 kpc) to the idealised toy model with a 50 km s −1 kick from Section 3.1, as shown in the middle row of Figure 2. We resample 10 7 particles within a range of 1144 < J ϕ < 2738 kpc km s −1 , to cover a larger range of RG (4.8 ≲ RG ≲ 11.4 kpc), but we do not include the effect of dust extinction or Gaia observational uncertainties.
Figure 10 shows the same angular momentum split as Figure 9 but for the toy model.The similarity for higher JR stars is striking, in that the arcs in the middle and lower rows closely resemble the data, with the middle rows capturing only the arcs at high vR and low | ∆R |, for stars close to their Guiding radius, while the lower rows capture the arcs at negative ∆R, from stars at pericentre.The top row's of Figures 3 and 9 showing the stars at apocentre are also similar, with the flattening of the distributions towards a point at higher ∆R, as expected from Figure 2, although the model contains more arcs at high JR than are recovered in the data.
We then repeat the exercise with the N -body model, while also including the effect of dust extinction and Gaia observational uncertainties.We use the stellar population synthesis code Snapdragons4 (Hunt et al. 2015(Hunt et al. , 2019b) ) to construct mock Gaia data from the N -body particles.Because M1 is a pure N -body simulation and the particles do not have a self-consistent age or metallicity, we first assign a simple random age and metallicity to each particle, 1 < τ < 12 Gyr and 0.0019 < Z < 0.019.This will not result in a realistic spread of stellar populations, nor a realistic age or metallicity gradient across the galaxy.However, it should be sufficient for this simple illustration of the effects of extinction and error.
We then use Snapdragons to sample stars from the Nbody particles using the Padova isochrones (Marigo et al. 2008), and a Kroupa IMF (Kroupa 2001), and convolve them with extinction calculated from a 3D version of the extinction map from Schlegel et al. (1998) and Gaia errors based on 34 months of observations (see Hunt et al. 2015;Grand et al. 2018, for a full explanation of the mock data synthesis process).We then select stars with a fractional parallax error of less than 10 percent, and a line-of-sight velocity error of less than 5 km s −1 to match our Gaia sample.This leaves us with a sample of 17,790,885 stars drawn from the simulation, compared to the 16,673,097 which were selected from the Gaia data.Figure 11 shows the same as Figure 9 but for the mock data (selected again with d < 2 kpc, and nine bins in Guiding radius, where | RG − RC |< 0.5 kpc) generated from model M1 with Snapdragons.We would not expect it to quantitatively reproduce the Gaia data, as again, this is not a model of the Milky Way, and has not experienced the same perturbation history.However, we recover a clear phase spiral in the upper right and middle left angular momentum bins, and fragments of spirals, in the other bins, qualitatively similar to the Gaia data (although the lower row, with RG outside the Solar radius only contain arcs at negative vR, unlike the data).
We further illustrate the phase dependence and bias with the diagram reproduced from Figure 6 of Darragh-Ford et al. (2023).Figure 12 shows an illustration of three stars that are currently in the observed sample (inner black box) at physical positions (R, ϕ) in the disk (star symbols), with their radial epicycles (large circles) and guiding centres (RG, θ ϕ ) shown (dots).The 'red' star with low angular momentum, J ϕ , must be on the outer part of its epicycle to appear in the local volume.Conversely, the 'blue' star with high J ϕ must be on the inner part of its epicycle to appear in the local volume.In addition, stars with their guiding centres further ahead (blue) or behind (red) the Solar neighborhood must have higher radial action, JR, and a specific radial phase angle, θR, to appear in the local volume.A star with its guiding centre within the local volume (black) must have sufficiently small JR to be physically located in the local volume (see Figure 6 and ascociated discussion in Darragh-Ford et al. 2023).
While the selection effects prevent direct observation of a complete ∆R − vR phase spiral in the data, the fact that the arc fragments are similar in Figures 9, 10 and 11, and we know that the features in the toy model in Figure 10 are parts of coherent ∆R − vR phase spirals as shown in Figure 2, we consider it likely that the arcs in the data are part of partially observed ∆R − vR phase spirals, arising from ongoing radial phase mixing, that will be fully resolvable with future surveys.However, note that the data shows much finer substructure than either the N -body simulation or the toy model, especially at low JR, where Figure 9 and Figure 10 are significantly different.This is owing to a combination of factors.Firstly, both the simulation and the toy model are limited in resolution which makes it difficult to resolve fine grained phase space structure.Secondly, as mentioned above, the disc in Hunt et al. ( 2021) is 'hot' with a Toomre Q parameter of 2.2 which suppresses substructure formation while the toy model lacks the self gravity to form self-consistent substructure at all.We would expect bar resonances and spiral arms to produce kinematic substructure at low JR (such as the classical moving groups, see e.g.Trick et al. 2021;Sellwood et al. 2019;Hunt et al. 2019a) which are not included in Figure 10 (note that bar resonances also induce kinematic structure at high JR, e.g.Kawata et al. 2021).Thirdly, the simulation is a pure N -body model without gas or star formation, which are also not included in the toy model.The real Milky Way has continued to form stars at low JR, replenishing the inner parts of the ∆R − vR plane after earlier perturbations excite older stars to higher JR (although M1 does contain some low velocity moving groups in Figure 11

Utility and context of this space
As discussed above, the z − vz phase spirals are of interest because they allow us to directly see dynamics at work locally in the Galactic disc.Whether phase mixing results from a single perturbation such as Sagittarius, the response to a triaxial dark matter halo, or stochastic scattering, the orbital tracks in z − vz allow us to 'see' the potential in a similar fashion to how streams from dissolving clusters or dwarf galaxies shows us their orbits in the Milky Way.If we can 'see' orbits we can constrain the Galactic potential (e.g.Orbital Torus Imaging; Price-Whelan et al. 2021;Horta et al. 2023).
Whether the ∆R − vR planes in the data are showing partial phase spirals or simply arcs, they still illustrate the orbital tracks and contain information on the radial gradient of the Galactic potential.For example, the smooth background in ∆R − vR can be exploited with a variation of Orbital Torus Imaging (Price-Whelan et al. 2021) to measure the radial profile of the Galactic potential, while the pitch angle and amplitude of the ∆R − vR phase spiral fragments can be used to 'time' the perturbation in the same fashion as being done in studies of the z − vz phase spirals (e.g.Antoja et al. 2023;Frankel et al. 2023;Darragh-Ford et al. 2023).
If the z − vz phase spirals and the ∆R − vR phase spirals originate from the same perturbation, combining the information from both should be more discriminatory than either dimension alone.

SUMMARY
In this work we explore the radial ∆R − vR equivalent of the well known z − vz phase spirals.Our conclusions are as follows: • Using a toy model and a simulation of a dwarf galaxy merging into a disc galaxy (Hunt et al. 2021), we show that radial phase mixing of disc stars following a perturbation induces spirals and arc-like substructure in the ∆R − vR plane, analogous to the well known z − vz phase spirals.The radial frequency of stars, ΩR, is a much weaker function of radial action, JR, than the vertical frequency Ωz is as a function of vertical action, Jz.The ∆R − vR phase spirals therefore take much longer to phase mix away (see Figure 3).
• We present similar substructure in the Gaia DR3 data.In a sample of stars in the very local solar neighbourhood (i.e.within d < 200 pc of the Sun), the ∆R − vR distribution is dominated by the solar neighborhood moving groups, which are local manifestations of extended dynamical structures.However, when exploring a large volume with discrete angular momentum bins we recover numerous arcs which resemble the substructure in the simulation and the toy model.
• Owing to observational selection effects we recover only partial ∆R − vR planes for varying angular momentum bins, making it unclear if they are partial recovery of coherent phase spirals or independent arcs.However, applying the same selection to the toy model of ∆R − vR phase spirals qualitatively matches the structure observed in the data.
• Regardless, these arcs and ridges follow the approximate orbital paths of stars on their radial epicycles provid-ing us an opportunity to directly observe dynamics at work in the Solar neighborhood, and constrain the Galactic potential in the radial direction (e.g.Price-Whelan et al. 2021).
• While we discuss the results in the context of a satellite induced perturbation, any radial perturbation which takes the stellar distribution out of equilibrium will of course result in radial phase mixing, and an ∆R − vR phase spiral.As discussed in the introduction, there are many currently proposed theories which can explain the z −vz phase spirals; e.g. the Sagittarius dwarf galaxy, buckling of the Milky Way bar, transient spiral structure, a triaxial/lumpy dark matter halo, or stochastic perturbation of the disc by a population of dark matter subhaloes.Such explanations will cause different features in the radial kinematics of stars, and if we can build models of such phase mixing which are coupled in z − vz and ∆R − vR, we may be able to break the degeneracy between competing theories, but we defer this to future work.
Future data releases from Gaia, complemented by spectroscopic surveys such as SDSS-V Milky Way Mapper (Kollmeier et al. 2019) will increase the number of stars for which we have 6 dimensional phase space information, and extend the coverage of this data across a larger fraction of the Galactic disc.This may enable us to map full ∆R − vR phase spirals, and apply methods to 'unwind' them similar to recent work on the z − vz phase spirals (Antoja et al. 2023;Frankel et al. 2023;Darragh-Ford et al. 2023).For the moment, it may still be possible to fit models to these partial ∆R − vR phase spiral fragments (e.g. as done for z − vz in Widmark et al. 2022b), or to use the background ∆R − vR distribution to constrain the radial component of the Galactic potential (similar to Orbital Torus Imaging as presented in Price-Whelan et al. 2021;Horta et al. 2023).
Figure 1.Left: Four example orbits, each with the same J ϕ and J R = 0 and increasing Jz (dark to light).Right: Six orbits with identical J ϕ and Jz = 0, but increasing J R (dark to light).Note that in both cases for orbits with higher action, the trajectories remain nearly closed but increasingly diverge from ellipses in a different way.

Figure 2 .
Figure 2. Toy model of 2 × 10 7 particles in a fixed potential showing ∆R − v R phase spirals that form following an impulsive kick of 25 (top row), 50 (second row) and 75 (third row) km s −1 .The left column shows the initial equilibrium distribution, the second and third columns show the ∆R − v R phase spiral after 0.5 and 1 Gyr respectively, over the same range as Figure 9 for later comparison.The rightmost column shows the full radial extent of the ∆R − v R spiral after 1 Gyr.The bottom row shows the same sample kicked vertically by 20 km s −1 , for comparison.
The full time evolution of simulation M1 from Hunt et al. (2021) in 10 Myr snapshots is now freely available at SciServer 1 (Taghizadeh-Popp et al. 2020), along with the rest of the 'SMUDGE' (Satellite Mergers Usher Disc Galaxy Evolution) simulation suite M2, D1 & D2.The simulations can be accessed and analysed within the SciServer Compute app and the package for interfacing with the simulations is introduced and explained in a provided example Jupyter Notebook, available both on SciServer and GitHub 2 .

NFigure 4 .
Figure 3. Upper row: z −vz phase spiral selected within a local volume in the 'late time' simulation snapshot, shown in density (left) and as a function of Ωz (right).Lower row: ∆R−v R phase spiral within the same volume in the simulation, as a function of density (left) and Ω R (right).The z − vz phase spiral winds up more rapidly than the ∆R − v R spiral as the vertical frequency gradient with respect to Jz is much steeper than the radial frequency gradient with respect to J R , as illustrated in the right column, where Ωz and Ω R are shown on the same scale.

Figure 5 .
Figure 5. ∆R − v R planes for 900 (30 by 30) 1 kpc-square regions in the disc at snapshot 720 (t = 7.05 Gyr) from simulation M1 from Hunt et al. (2021).The colormap in each subpanel corresponds to the number of particles in each bin of the ∆R − v R phase space, with darker color corresponding to more particles.The structure of the perturbed phase-space distribution is more irregular at larger radius because the self-gravity of the disc is weaker.The time evolution of ∆R − v R and z − vz can be seen side by side in M1-vertical-radial at https://zenodo.org/records/8402668.

Figure 6 .
Figure 6.Same as Figure 5 but occurring earlier in the simulation at time 4.45 Gyr, 2 Gyr after the first passage of the satellite galaxy, but before the second passage (see Hunt et al. 2021, for the orbit of the satellite).
Figure 7. v R − v ϕ (left) and v R − ∆R (right) for a 200 pc volume around the Sun.Note that over such a local volume, they are almost identical.The classical moving groups are marked in the left panel.

Figure 8 .
Figure 8. Left: ∆R − v R plane stars with physical distance d < 2 kpc from the Sun.As with Figure 3, the colormap shows the residual phase-space density of stars after unsharp-masking the data (i.e. after subtracting a smoothed version of the phasespace density from itself).The classical moving groups, Sirius (blue dot), Coma Berenices (green plus), Pleiades (orange cross), Hyades (yellow tri) and Hercules (red dash) are marked matching Figure 7. Right: R − v R plane for the same stars, illustrating the difference in R − v R and ∆R − v R for a local volume Figure 9. Substructure in the radial phase space, ∆R − v R , for Gaia data within d < 2 kpc of the Sun and in nine bins of Guiding radius.In each bin, we select stars with | R G − R C |< 0.5 kpc.Note that this figure has an implicit radial action and radial angle dependence, because stars will only appear if they are at the appropriate part of their radial epicycle.The physical selection limits mean that for stars with low Lz we only see them close to apocentre (e.g., the top left panel), and for stars with high Lz we only see them close to their pericentre (e.g., the lower right panel).

Figure 10 .
Figure10.Same as Figure9but for the toy model shown in Figure2with a 50 km s −1 kick, after 1 Gyr.Note that the toy does not contain observational uncertainties, selection effects from extinction, or secular effects from a galactic bar or spiral arms.

N
body simulation with extinction and Gaia observational uncertainties

Figure 11 .
Figure 11.Same as Figure 9 but for the N -body simulation M1, after splitting the N -body particles into stars and adding extinction and Gaia observational uncertainties with the population synthesis code Snapdragons.

Figure 12 .
Figure 12.Illustration of three stars that are currently in an observed sample (black box) at physical positions (R, ϕ) in the disk (star), with their radial epicycles and guiding centres (R G , θ ϕ ) shown (dots).The diagram shows that there are biases in radial action, J R , and radial phase, θ R as a function of azimuthal action, J ϕ , and phase θ ϕ , when selecting stars in a physically local sample (Reproduced from figure 6 of Darragh-Ford et al. 2023).
, as also shown in vR − v ϕ in Hunt et al. 2021).