© The Authors 2023

1. Introduction

The Gaia satellite (Gaia Collaboration 2016, 2018, 2022c; de Bruijne 2012) constitutes the most advanced astrometric mission to date. After its launch in 2013, it has been providing positions, parallaxes, proper motions, and line-of-sight (LoS) velocities for an increasing number of sources in subsequent data releases (Katz et al. 2004, 2019, 2022; Cropper et al. 2018). This exquisite astrometry has improved our understanding of Galactic structures that are already known, such as the spiral arms and the warp (Antoja et al. 2016; Poggio et al. 2018, 2021; Chrobáková et al. 2022; Gaia Collaboration 2022b), in addition to revealing a complex formation for the Milky Way and strong interactions with other galaxies (Belokurov et al. 2018; Myeong et al. 2018, 2019; Helmi et al. 2018; Koppelman et al. 2019; Helmi 2020). In this context, many asymmetries in different parameter spaces have been interpreted as a consequence of this scenario, including velocities (Antoja et al. 2017), distribution of proper motions (Palicio et al. 2020), ridges in projected velocities (Ramos et al. 2018; Fragkoudi et al. 2019; Khoperskov & Gerhard 2022; Gaia Collaboration 2022b,a; McMillan et al. 2022), distribution of actions (Hunt et al. 2019; Sellwood et al. 2019; Bland-Hawthorn et al. 2019; Trick et al. 2019, 2021; Trick 2022), and distribution of metallicity (Poggio et al. 2022). In this work, we report on structures in the Galactic plane revealed by the distribution of the radial action, J_R, computed with the Gaia DR3 astrometry and LoS velocities (Gaia Collaboration 2022c).

This Letter is organised as follows. In Sect. 2, we explain the selection criteria applied to the Gaia data of our sample. In Sect. 3, we describe the model and the performance adopted for computing the orbital parameters and actions from the input data. Results are shown and discussed in Sects. 4 and 5, respectively. The conclusions can be found in Sect. 6. Finally, in Appendices A and B, we specify the data query performed on the Gaia archive¹ and the detailed procedure for the estimation of the actions and orbital parameters, respectively. In Appendix C, we reproduce our analysis with a subsample of giants stars selected by photometry.

2. Gaia data and selection criteria

We made use of all the Gaia DR3 stars with full astrometric information available (parallaxes, positions, proper motions and line of sight velocities) and select those with non-null geometric distance estimation (Bailer-Jones et al. 2021). The corresponding ADQL query can be found in Appendix A. This sample totals 33 653 049 million sources – of these, we select the ones with good kinematic measurements by imposing a maximum LoS velocity error of 5 km s⁻¹ and a relative error in proper motion lower than 15%. For the heliocentric distance, we imposed a maximum relative error of 20%. Since we have focussed our study on the disc, we excluded those stars whose maximum distance from the Galactic plane is larger than 500 pc (see Sect. 3). The resulting sample size is ∼12.4 million sources. We corrected the LoS velocities and proper motions, assuming (U_⊙, V_⊙, W_⊙)=(9.5, 250.7, 8.56) km s⁻¹ for the solar motion (GRAVITY Collaboration 2021; Reid & Brunthaler 2020) and R_⊙ = 8.249 kpc for the Galactocentric distance of the Sun (GRAVITY Collaboration 2021).

In order to propagate the errors, we considered the correlations between the astrometric parameters. Motivated by the approach of Gaia Collaboration (2022b) and Kordopatis et al. (2022), we modelled the distribution of the errors of the geometric distances with a broken Gaussian distribution parameterised by the input confidence intervals (i.e. r_lo and r_hi in Bailer-Jones et al. 2021).

3. Orbital parameters and actions

We modeled the forces of the Milky Way with a rescaled version of the potential of McMillan (2017) such that the circular velocity at R_⊙ = 8.249 kpc is V_⊙ = 238.5 km s⁻¹, consistent with our assumed Solar motion and the velocity of the Sun with respect to the local standard of rest taken from Schönrich et al. (2010). This potential is fully axisymmetric and models the contribution of the halo, bulge, and thin and thick stellar discs as well as the H I and H II gas discs. We estimate the orbital parameters (apocenter, r_apo, pericenter, r_peri, and maximum orbital distance to the galactic plane, Z_max) and the non-trivial actions, J_R and J_Z, by using own implementation of the Stäckel-Fudge approximation (Binney 2012; Sanders & Binney 2016; Mackereth & Bovy 2018). We refer to Appendix B for a detailed description of this procedure. Apart from these parameters, the vertical component of the angular momentum, L_z, and the total energy, E, are obtained as output. The actions, J_R, J_Z, and the angular momentum, L_z, presented in this work are expressed in units of L_⊙ = R_⊙V_⊙.

4. Results

In this section, we explore the map of the distribution of the radial action J_R in the Galactic Plane. Figure 1 shows the spatial distribution of the median J_R, while each panel in Fig. 2 refers to other percentiles to illustrate the variation of the distribution of J_R across the Galactic plane. Due to the variations in the observed trends as a function of the considered percentile, the colorbar is tuned to enhance the contrast between the high and low J_R regions in each panel. We identify three main structures in the low J_R regions for the first four percentiles shown in the figure (labelled as A, B, and C in Fig. 1), whereas for the 94th percentile, they are highly distorted. We observe an additional feature (labelled D) in the outer part of the disc (10 kpc ≲ R ≲ 11 kpc) characterised by high J_R values. The innermost structure, labelled A, extends from R ≈ 6.0 kpc at (X, Y) ≈ (−2, −5.5) kpc to R ∼ 7 kpc at the solar azimuth (X = 0 kpc direction), while for X <  0 kpc it shows an almost constant radii of R ≈ 7 − 7.2 kpc. This results in a longitudinally asymmetric arc-shape structure of variable pitch angle.

Fig. 1. Distribution of the median (equivalent to the P50 percentile) of JR on the Galactic Plane (|Zmax|< 0.5 kpc). The solid black circle denotes the solar position. Features labelled A to D are discussed in the text.

Fig. 2. Distribution of JR percentiles on the Galactic Plane (|Zmax|< 0.5 kpc). The x percentile is denoted by Px. The solid black circle denotes the solar position.

Structure B also shows significant variations with longitude²: for X <  0 we observe a well-defined low J_R area that extends from (X, Y) ≈ (−4, −6.5) kpc to (0, 8.5) kpc, embedding the solar neighbourhood. However, its prolongation at negative X is highly distorted, resulting in a wide area of low J_R between structure A and the ℓ = −90° direction.

In contrast to the previous features, structure C is sharply defined at negative X, where it extends from (X, Y) ≈ (4, −9) kpc to (−2, −10) kpc, although it is possible to discern a tail of relatively low J_R at X <  −2 kpc. At positive X, this structure is connected with one of the extensions of the feature B, located in the large low J_R area found between A and B (X >  0 kpc and −8 ≲ Y ≲ −7 kpc), and creating a gap of high J_R with B. As can be seen in Fig. 2, this feature extends towards outer radii for percentiles larger than P₈₃, and constitutes the only low J_R structure at large percentile (P₉₄).

The outermost feature (D) is a high J_R region with an arc-shape of almost constant radii of 10.5 kpc and ∼1.0 kpc width. It remains almost unchanged for percentiles lower than P₇₇ and becomes blurred for higher values. Apart from the main features, it is worthwhile mentioning the bifurcation in structure A at (X ≲ −2 kpc, Y ≈ −6 kpc) towards positive X, although it gets distorted in the maps for the large percentiles (P₆₆ and above). Finally, we can discern a subtle arc-shape structure between A and B, with very low median radial action (P₅₀ <  0.008) from (X, Y)≈(0, −7.7) to (2, −7.7).

In order to check whether the features described above are a consequence of the distribution of stars, we represent in Fig. 3 the density map of the selected sample. As it can be seen in the left panel, the density map cannot explain all the structures identified in Figs. 1 and 2. The density map peaks at the solar position and decreases with the heliocentric distance as fainter stars are excluded, showing no correspondence with the arc-shaped structures in the J_R distribution.

Fig. 3. Density map of the selected sample in the Cartesian plane (X, Y) on the left. Distribution of median errors in JR on the Galactic plane on the right.

We verify the significance of the features in J_R with the observational errors by evaluating the map of the median error of the radial action, δJ_R, estimated from 25 realisations of the input data (right panel in Fig. 3). Although it is possible to distinguish some selection effects in an annular region centred at the Sun, the structure associated with them does not correspond to that reported in Figs. 1 and 2. Furthermore, in the vast majority of the plane, the errors of J_R are at least 3.5 times smaller than the median J_R, supporting the robustness of the features found in the percentile distributions.

5. Discussion

In this section, we discuss three possible scenarios to explain the observed features in J_R.

5.1. Spiral arms

The spatial distribution and shape of the structures reported above suggest a connection with the spiral arms. To explore this hypothesis, we compare these structures with the fit of the spiral arms inferred from the kinematics of one hundred masers (Reid et al. 2014) from the distribution of Cepheids (Lemasle et al. 2022) and from the distribution of Gaia EDR3 upper main sequence stars (UMS stars, Poggio et al. 2021), which considers the same astrometric measurements (but for a different sample) as this work. We complement these references with the overdensity map of our subsample of giant stars (see Appendix C). Following the procedure described in Poggio et al. (2021), we compute the local (average) density using an Epanechnikov kernel (Epanechnikov 1969) of bandwidth 0.3 kpc (2.0 kpc). This kernel assumes that the contribution of a star to the density at a reference point is weighted by a term ∝max{1 − x²/h², 0}, where h is the bandwidth and x is the separation between the star and the reference point. For sake of visualisation, the references of spiral arms described above are shown in individual panels in Fig. 4.

Fig. 4. Comparison of the maps of P50(JR) with the spiral arms reported in literature. First panel: contour lines enclose the overdensities found in the subsample of giants. Second panel: solid lines represent the Scutum (cyan), Sagittarius (yellow), Local (blue), and Perseus (black) spiral arms of Reid et al. (2014), while the dotted lines correspond to their extrapolation in azimuth. Third panel: Solid lines represent the segments of spiral arms of Lemasle et al. (2022), in which their same naming convention is used, while the colorcode results from a visual comparison with these of Reid et al. (2014). The additional structures are indicated by red and pink lines for description convenience. Fourth panel: contour lines illustrate the overdensities reported by Poggio et al. (2021). Background image: reproduction of Fig. 1 using a grey color-scale to increase the contrast between the coloured lines and the background map. The solid white circle denotes the solar position.

The first panel in Fig. 4 illustrates the overdensities in the distribution of our sample of giants. We find a correspondence between the overdensities in this sample and these reported by Poggio et al. (2021) for the younger UMS population (fourth panel); although some discrepancies are observed at (X, Y)≈(−1, −9) kpc and (−2,−6.5) kpc. The presence of the spiral arms traced by an old population has been recently proposed by Lin et al. (2022), who identified the Local Arm in a sample of 87 000 Gaia EDR3-2MASS (Gaia Collaboration 2021; Skrutskie et al. 2006), Red Clump stars (RC) and could be related to the metallicity asymmetry in Sample B and C in Poggio et al. (2022, see their Fig. 1).

In general terms, we find a good agreement between the low J_R areas and the spiral arms, especially in the innermost regions, where the distribution of giant stars (first panel) reveals an overdensity consistent with Structure A. Furthermore, the bifurcation observed in A can be explained by segments 16 and 22 of Lemasle et al. (2022), likely to be part of Sagittarius and Scutum, respectively. On the contrary, we find a shift of ∼0.5 kpc between structure A and the location of segment 9, where the extrapolation of Reid et al. (2014) fits the lowest J_R region of A (dotted lines) perfectly. Compared to Poggio et al. (2021), we can identify most of Structure A in the innermost overdensity of UMS stars, although no bifurcation is observed at X ≈ −2 kpc. The extension of this overdensity, however, is compatible with the area of low J_R that connects the structures A and B in Fig. 1.

As mentioned in Sect. 4, we find a subtle arc-shape structure between A and B close to the solar neighbourhood. This feature has no counterpart in the spiral arms of Reid et al. (2014) and Poggio et al. (2021) or in our distribution of giant stars, but it is located at the same position as the segment 18 (pink line) of Lemasle et al. (2022), being a potential continuation of the Sagittarius arm. As the percentile increases (Fig. 2), this small area of low J_R becomes more evident (a gap with A emerges) and consistent with the segment 18 and its extension towards positive X.

The part of structure B located at negative X is compatible with the fit of the Local spiral arm of Reid et al. (2014), the segment 23 of Lemasle et al. (2022) and the overdense regions found in the UMS and giant population. On the contrary, at positive X only the Local arm of Reid et al. (2014) might provide a good explanation for structure B, but only if a shift of ∼0.5 kpc is considered. It is worthwhile mentioning the significant differences among the references for that part of the Local arm: assuming the segment labelled 17 is part of the Local arm, it implies a pitch angle of opposite sign compared to that in Poggio et al. (2021); whereas according to Reid et al. (2014), the Local arm is more tangential. This variety of observations suggests a complex definition of the extension and limits of the Local spiral arm despite its proximity to the Sun.

One major discrepancy is found in the solar neighbourhood: according to our maps, the Sun is embedded in the intersection of the Local and the Sagittarius spiral arm, while the predictions of all three spiral arms maps report a solar location in the inner boundary of the Local Arm.

The J_R maps suggest a connection between the Perseus and the Local Arm (Structures C and B, respectively). However, the spatial geometry of the spiral arms from UMS stars (Poggio et al. 2021), Red Clump stars (Lin et al. 2022), and the giant sample does not coincide with the observed features in J_R in this region.

The comparison of structure C reveals a good agreement with the Perseus spiral arm of Reid et al. (2014) for X ≲ 0 kpc and the segment 12 of Lemasle et al. (2022) within |X|≲1 kpc. However, at positive X, structure C exhibits a different pitch angle, as compared to both Reid et al. (2014) and Lemasle et al. (2022).

It is worth mentioning that the spiral structure of the Milky Way might be different depending on the considered stellar population. Here, the contrast in J_R is observed mainly in the giant old population (see Appendix C for the specific analysis of the giant stars), even though the stars in the spiral arms tend to be young and, through the age–velocity dispersion relation, exhibit lower values of J_R. For instance, the referred spiral arms have been traced by selecting masers, Cepheids and UMS stars, namely, the young population. Thus, the dynamics of the old stars seem to be in agreement with the spatial distribution of the young population in some regions, but present some discrepancies in others. Such discrepancies can be either due to the fact that the geometry of the spiral arms might be different for different stellar populations or that the dynamical nature of the spiral arms somehow leads to the observed features.

5.2. Moving groups

We also explored the possible origin of the reported structures in the moving groups. As Ramos et al. (2018) have shown, it is possible to identify the moving groups as stripes in the azimuthal velocity, V_ϕ vs. R, diagram. Figure 5 represents the distribution of the median J_R in the (R,V_ϕ) plane, including some of the moving groups reported by Ramos et al. (2018) as reference (yellow dashed lines). For sake of visualisation, we focus on the range 220 <  V_ϕ <  250 km s⁻¹ and use a logarithmic colorscale for median(J_R) to enhance the features. As expected, the values of J_R tend to increase as V_ϕ differs from the rotation curve. As Fig. 5 shows, the Dehnen98-6, Hyades and Sirius moving groups are predominantly located in areas of relatively high J_R in the (R,V_ϕ) plane, in contrast to the low J_R values that characterise the features described in Sect. 4. On the contrary, Coma Berenices lies close to a transition from low to high J_R. The Hercules and most of the Horn-Dehnen98 moving groups lie in the region of high J_R (blue saturated region) and we do not see any clear correspondence for the Arch1-Hat moving group at this point. In any case, these groups could be related to the structures of J_R in the R − V_ϕ projection that extend to higher J_R (not seen in our figures due to the colour range).

Fig. 5. Azimuthal velocity Vϕ vs. R diagram colorcoded with the median JR. The colorbar has been intentionally set in logarithmic scale to cover a wide range of values in JR. The moving groups (dashed yellow lines) are displayed from the bottom left to the upper right corner as follows: Hercules, Dehnen98-6, Horn-Dehnen98, Hyades, Coma Berenices, Sirius, and Arch1-Hat. Black ellipses enclose the two selected areas (see the text), while the Sun is denoted by the solid black circle.

Apart from the ridges, we can identify two interesting areas of low J_R: one located between the Dehnen98-6 and the Hyades moving groups (as a prolongation of Horn-Dehnen98 at inner radii) and another more extended low J_R area close to Coma Berenices. In order to evaluate the contribution of this potential members of moving groups, we exclude the stars within these regions (black dashed ellipses in Fig. 5). We have verified the exclusion of the stars close to Coma Berenices raises the median values of J_R at ∼(2, −8) kpc, improving the separation between the A and B structures at positive X. The exclusion of the other selection, however, leads to an annular distortion at R ∼ 7.7 kpc that increases the gap between the A and B structures, especially at negative X. This distortion, however, is more likely to be an artefact caused by the exclusion of a significant number of sources within 7.6 ≲ R ≲ 8.0 kpc, rather than a true contribution of the Horn-Dehnen moving group.

Based on our tests, the features observed in the V_ϕ vs. R plane are not as clear as those found in the maps of J_R, although an apparent relation between the high J_R values and the position of some ridges can be inferred. A deeper analysis of this relation is needed to evaluate the contribution of the moving groups to the features in J_R(X, Y) and (potentially) its connection with the spiral arms. Such an analysis is beyond the scope of this Letter and will be explored in a future work.

5.3. Galactic bar

Apart from the spiral arms, the location and shape of the high J_R region at R ∼ 10.5 kpc is consistent with some values reported for the Outer Lindblad Resonance (OLR; Liu et al. 2012; Portail et al. 2017; Pérez-Villegas et al. 2017). In order to evaluate this possible connection, we must verify whether the OLR corresponds to a region of increasing radial action. Under the epicyclic approximation (see Binney & Tremaine 2008), the Galactocentric distance R(t) of a star trapped by a resonance varies with time as:

$$\begin{array}{c}\hfill R(t)=R_{\mathrm{g}}-C_{2}cos(2\mathrm{\Delta }\mathrm{\Omega }t)(\mathrm{\Delta }\mathrm{\Omega }\equiv \mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{p}}),\end{array}$$(1)

where the factor of 2 in the cosine comes from the assumption of a dipolar disturbance of the potential (m = 2), R_g is the guiding radius, Ω is the circular frequency at R_g, Ω_p is the pattern speed, and C₂ is a constant that depends on the bar potential Φ_b(R, ϕ, t) as:

$$\begin{array}{c}\hfill C_2=\frac{1}{{\kappa }^2-4\mathrm{\Delta }{\mathrm{\Omega }}^2}\times {\left(\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{b}}}{\mathrm{d}R}+\frac{2\mathrm{\Omega }{\mathrm{\Phi }}_{\mathrm{b}}}{R\mathrm{\Delta }\mathrm{\Omega }})\right|}_{R_{\mathrm{g}}},\end{array}$$(2)

with κ the epicyclic frequency at R = R_g. Differentiating Eq. (1) with respect to the time and substituting in the integral for J_R (see Eq. (B.11)) we have:

$$\begin{array}{c}\hfill J_{\mathrm{R}}=\frac{2\mathrm{\Delta }\mathrm{\Omega }{C}_2}{\pi }{\displaystyle {\int }_{R_{\mathrm{g}}-C_2}^{R_{\mathrm{g}}+C_2}}sin(2\mathrm{\Delta }\mathrm{\Omega }t)\mathrm{d}R=C_2^2·\mathrm{\Delta }\mathrm{\Omega },\end{array}$$(3)

where the upper and lower limits correspond to the cases in which cos(ΔΩt) = −1 (apocenter) and +1 (pericenter), respectively. Thus, Eq. (3) diverges in the Outer Lindblad Resonance (κ + 2ΔΩ = 0). Although Eq. (1) assumes a small deviation in azimuth with respect to the circular orbit defined by the guiding radius, which is not true in the resonance regime, it is enough to demonstrate the radial action increases towards the resonances (Chiba & Schönrich 2021). A more detailed analysis, such as that described for the corotation in Sect. 3.3b of Binney & Tremaine (2008), would predict a large but finite action. However, the calculus of this more general case is not straightforward (Goldreich & Tremaine 1981).

According to Sellwood & Binney (2002) and Sellwood (2010), it is not only the OLR, but also the inner Lindblad resonance (ILR), that should be characterised by an increment in J_R due to the outward flow of angular momentum (Lynden-Bell & Kalnajs 1972). Assuming a pattern speed for the bar between 34 and 47 km s⁻¹ kpc⁻¹ (Bland-Hawthorn & Gerhard 2016, and references therein), the ILR is expected to be located out of our region of study.

6. Conclusions

The statistics of the radial actions reveal arc-shape structures in the Galactic disc. These structures are characterised by a predominance of more circular orbits that contrasts to the high radial action feature found at R ∼ 10.5 kpc. The analysis of the errors in J_R confirms the reported structures are not spurious but robust from the statistical point of view. Furthermore, they cannot be explained by the selection effects inherent in Gaia.

The characteristic arc shape of the structures in J_R motivates the comparison with the Milky Way spiral arms, whose fit parameters have been reported in previous studies. We find that in the innermost region, structure A clearly defines the Sagittarius arm, with its upper boundary is delimited by the Scutum arm. At larger Galactocentric radii, structure B tracks the Local Arm at negative X, while no clear correspondence with literature is found at X >  0, where the variety of models suggests a complex definition for this arm. On the contrary, for the Perseus Arm we observe a good concordance with the spatial distribution of young stellar population for X ∈ ( − 2, 0) kpc, while at positive X, the orientation of the J_R feature has a different pitch angle compared to all the considered models. Our results suggest that the Perseus Arm in the J_R map is connected to the Local Arm at ∼3.6 kpc from the Sun, in the direction of ℓ ≈ −100°. This would result in a mismatch with some geometries of the spiral arms from young stellar populations, which will be studied in the future. We observe a correspondence between the segment 18 in Lemasle et al. (2022) and a region of very low J_R between structures A and B that has not clear spiral arm assignation.

We also explored the moving groups as a possible explanation for the features. The J_R arc-shape structures in the (X, Y) plane are likely to be related to the structures in J_R in the R − V_ϕ plane but mapped onto different projections of phase space, in particular, showing also their complex dependency with position (e.g., azimuth) in the (X, Y) case. We observe some features in the V_ϕ vs. R plane that might be anti-correlated with some known moving groups. However, this connection between the moving groups and the J_R features in the Galactic plane, if present, is not obvious and should be explored in future studies.

We identify an area of high radial action centered at ∼10.5 kpc, where the outer Lindblad resonance (OLR) caused by the bar is expected. Apart from the features in the maps of the radial action, we find the distribution of the giant stars in the disc is consistent with the spiral arms traced by younger populations, in particular: the upper main sequence stars.

The analysis presented in this work indicate that multiple agents might be causing the structures found in the distribution of J_R. Although the spiral arms account for most of the features reported in this work, there are still many discrepancies that must be addressed. In this context, further studies with numerical simulations and analytical models are required to explain these differences and shed light on the Galactic dynamics.

Acknowledgments

The authors acknowledge J. Bland-Hawthorn for his constructive contribution to this work as referee. We thank P. de Laverny for his useful comments. P. A. Palicio acknowledges the financial support from the Centre national d’études spatiales (CNES). E. Spitoni and A. Recio-Blanco received funding from the European Union’s Horizon 2020 research and innovation program under SPACE-H2020 grant agreement number 101004214 (EXPLORE project). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N. 101063193. TA acknowledges the grant RYC2018-025968-I funded by MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your future”. This work was (partially) funded by the Spanish MICIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” by the “European Union” through grant RTI2018-095076-B-C21, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. PJM acknowledges project grants from the Swedish Research Council (Vetenskaprådet, Reg: 2017- 03721; 2021-04153). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Although GALPY is not explicitly used in this work, P. A. Palicio uses its source code as reference and recognises the credit for the work of Bovy (2015).

References

Appendix A: ADQL query

SELECT source_id, ra, dec, pmra, pmdec,
 radial_velocity, parallax, ruwe, ra_error,
 dec_error, pmra_error, pmdec_error,
 radial_velocity_error, parallax_error,
 ra_dec_corr, ra_pmra_corr, ra_pmdec_corr,
 dec_pmra_corr, dec_pmdec_corr, pmra_pmdec_corr,
 grvs_mag,r_med_geo, r_lo_geo, r_hi_geo
FROM user_dr3int6.gaia_source INNER JOIN
 external.gaiaedr3_distance USING(source_id)
WHERE (radial_velocity is not NULL) and (pmra is not
 NULL) and (pmdec is not NULL)

Listing 1. ADQL query for the Gaia DR3 considered in this work.

Appendix B: Stäckel-Fudge approximation

Within this approach, the orbital parameters can be computed assuming the considered Galactic potential Φ(R, z) satisfies some properties of the so-called Stäckel potentials. Given an axisymmetric oblate distribution of mass, its potential Φ(R, z) is said to be a Stäckel potential if there are two single-variable functions U(u) and V(v) such that:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_S(u,v)=\frac{U(u)-V(v)}{{sinh}^{2}u+{sin}^{2}v},\end{array}$$(B.1)

where (u, v) are the ellipsoidal coordinates (de Zeeuw 1985) related to (R, z) through the transformation:

$$\begin{array}{c}\hfill R=\mathrm{\Delta }sinhusinvz=\mathrm{\Delta }coshucosv,\end{array}$$(B.2)

with Δ as the focal length of the elliptical (hyperbolic) curves of constant u (v). Since the Galactic potential is known to be oblate, we do not describe the prolate case (for the prolate case see de Zeeuw 1985). By differentiating both sides of Eq. B.2 with respect to time, the transformation of the momentum between (R, z) and the (u, v) coordinate systems is as follows:

$$\begin{array}{cc}\hfill p_u& =p_{\mathrm{R}}\mathrm{\Delta }coshusinv+p_z\mathrm{\Delta }sinhucosv,\hfill \\ \hfill p_v& =p_{\mathrm{R}}\mathrm{\Delta }sinhucosv-p_z\mathrm{\Delta }coshusinv,\hfill \end{array}$$(B.3)

where p_i is the momentum associated with the coordinate i ∈ {R, z, u, v}. The Hamiltonian constructed with the momenta of Eq.B.3 and the potential Φ_S(u, v) results in a expression that can be separated into two single variable terms:

$$\begin{array}{cc}\hfill E{sinh}^{2}u& =\frac{p_u^2}{2{\mathrm{\Delta }}^2}+U(u)+I_3+\frac{L_z^2}{2{\mathrm{\Delta }}^2{sinh}^{2}u},\hfill \\ \hfill E{sin}^{2}v& =\frac{p_v^2}{2{\mathrm{\Delta }}^2}-V(v)-I_3+\frac{L_z^2}{2{\mathrm{\Delta }}^2{sin}^{2}v},\hfill \end{array}$$(B.4)

in which E is the total energy of the system (since the Hamiltonian does not depend explicitly on time), L_z is the vertical component of the angular momentum, and I₃ is the third integral of motion.

For a reference point with coordinates (u, v)=(u₀, π/2), the expression for u in Eq. B.4 is:

$$\begin{array}{cc}\hfill E{sinh}^2{u}_0& =\frac{{p_u^2}_0}{2{\mathrm{\Delta }}^2}+U(u_0)+I_3+\frac{{L_z}^2}{2{\mathrm{\Delta }}^2{sinh}^2{u}_0},\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill E& =\frac{p_{0.5\pi }^2}{2{\mathrm{\Delta }}^2}-V(\pi /2)-I_3+\frac{{L_z}^2}{2{\mathrm{\Delta }}^2},\hfill \end{array}$$(B.6)

where the choice for u₀ is discussed later. Subtracting Eq. B.4 from B.5 and solving for p_u we find:

$$\begin{array}{c}\hfill \frac{p_u^2}{2{\mathrm{\Delta }}^2}=\frac{{p_u}_0^2}{2{\mathrm{\Delta }}^2}+E({sinh}^{2}u-{sinh}^2{u}_0)-U(u)+U(u_0)\\ \hfill -\frac{{L_z}^2}{2{\mathrm{\Delta }}^2}(\frac{1}{{sinh}^{2}u}-\frac{1}{{sinh}^2{u}_0}),\end{array}$$(B.7)

where the term δU ≡ U(u)−U(u₀) can be approximated using the definition of the Stäckel potential (Eq. B.1) as:

$$\begin{array}{c}\hfill \delta U\equiv U(u)-U(u_0)\approx ({sinh}^{2}u+{sin}^{2}v)\mathrm{\Phi }(u,v)\\ \hfill left({sinh}^2{u}_0+{sin}^{2}v\mathrm{\Phi }(u_0,v).\end{array}$$(B.8)

Similarly, we can define δV = V(v)−V(π/2) such that:

$$\begin{array}{ccc}\hfill \delta V& \equiv \hfill & \hfill {cosh}^{2}u\mathrm{\Phi }(u,\pi /2)-\\ \hfill & -\hfill & \hfill ({sinh}^{2}u+{sin}^{2}v)\mathrm{\Phi }(u,v),\end{array}$$(B.9)

then p_v can be written as

$$\begin{array}{c}\hfill \frac{p_v^2}{2{\mathrm{\Delta }}^2}=E{sin}^2(v)+I_3+V(\pi /2)+\delta V-\frac{L_z^2}{2{\mathrm{\Delta }}^2{sin}^{2}v}.\end{array}$$(B.10)

Expressions B.7 and B.10 for p_u and p_v, respectively, can be substituted in the integrals for the definition of the actions (see Eq. 6 in Binney 2012):

$$\begin{array}{cc}\hfill J_u& =\frac{1}{2\pi }{\displaystyle \oint p_u}du=\frac{1}{\pi }{\displaystyle {\int }_{u_{\mathit{min}}}^{u_{\mathit{max}}}}{p}_{u}du,\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill J_v& =\frac{1}{2\pi }{\displaystyle \oint p_v}dv=\frac{2}{\pi }{\displaystyle {\int }_{v_{\mathit{min}}}^{\pi /2}}{p}_{v}dv,\hfill \end{array}$$(B.12)

Both the limits and the integrals of Eq. B.11 have to be computed numerically. The limits of integration correspond to the roots of Eq. B.7 and B.10 (therefore, the actions are always real). In our case, we computed these roots using the bisection method, while the numerical integration is performed by Gaussian Quadrature with ten nodes. We approximated the radial and vertical actions as J_R ≈ J_u and J_z ≈ J_v, respectively, since the R (z) coordinate varies more with u (v), as Fig B.1 illustrates. The choice of u₀ is rather arbitrary (see Section 2 in Binney (2012) for the discussion), so we use the coordinate u given by the input value (R, z) of the star.

Fig. B.1. Example of an orbit in the (R, Z) plane (blue curve) with the lines of constant u (dot-dashed ellipses) and v (dashed hyperbolas) in the background. The units of the axis are arbitrary.

In order to account for the error propagation, we performed 25 random realisations of the input data and compute the median values and the 16th and 84th percentiles of the output.

Appendix C: Results with the subsample of giants

Using a photometric selection of stars, we demonstrate that the dynamical pattern reported here is mainly supported by the old giant population, although the stars in spiral arms tend to be younger than average and, therefore, they have lower values of J_R according to the age-velocity dispersion relation. We reproduce the percentiles of J_R (shown in Fig. 2) by applying the following photometric criteria:

$$\begin{array}{c}\hfill G_{\mathit{RVS}}+5-5{log}_{10}{d}_{\mathit{pc}}<G_{\mathit{BP}}-G_{\mathit{RP}}-1,\end{array}$$(C.1)

which assumes no extinction as a first approximation to restrict the sample to the giants (hereafter, we refer to this subset as ’giant subsample’). The expression in Eq. C.1 visually separates the red giant branch (RGB) from the main sequence stars in the Hertzsprung–Russell (HR) diagram, using the red clump as reference for the boundary (Fig. C.1). By selecting giants, we are able to retain stars intrinsically brighter and reduce the effect of the selection function as well as the contribution from the faint dwarf stars that dominate the sample in the Solar neighbourhood (Gaia Collaboration 2022b).

Figure C.2 illustrates the distribution of the percentiles of J_R for the ’giant’ subsample. In general, the features found in the whole sample are observed in the giant subsample, with the exception of the high J_R region between the Local and Perseus (B-C) arm which is more distorted. Similarly, the high radial action region near the Sun disappears. In contrast to Fig. 3, for the giant subsample, the highest density area corresponds to the innermost low J_R region, although no evidence of the other structures are observed.

Fig. C.1. Distribution of stars in the Hertzsprung–Russell diagram for the Galactic disc sample (|Zmax|< 0.5 kpc) shown in the upper panel. The dashed black line represents the boundary condition considered to separate Main Sequence and RGB stars (Eq. C.1) neglecting extinction. Density map of the giant sample in the (X, Y) plane in the bottom panel. The solar position is denoted by a solid black circle.

Fig. C.2. Distribution of JR percentiles for the giant subsample. Details are the same as those in Fig. 2. Upper-left panel corresponds to the median (equivalent to Fig. 1), while two additional percentiles are shown in the central and right upper panels.

We performed an additional test to address the age of the tracers of the arc-shaped structures. We compared the kinematics of our giant sample with that of the Cepheids in Gaia DR3 (Ripepi et al. 2022) to get a proxy of the relative age of both populations. In Fig. C.3, we show the distributions of azimuthal (V_ϕ) and vertical velocities (V_Z) for the Cepheid and giant subsamples. As we can see, the distribution of V_ϕ for the Cepheids is more peaked than that of the giants, as expected for cooler (and younger) populations, with a median absolute deviation³ of σ_C = 36.05 km/s (σ_G = 20.41 km/s) for the Cepheid (giant) samples. Similarly, the distribution of V_Z observed for the giant subsample is consistent with a hotter and older population compared to the Cepheids, with σ_G = 22.56 km/s and σ_C = 8.56 km/s, respectively. Therefore, the distributions presented in Fig. C.3 indicate that the giant subsample is dominated by stars typically older than the Cepheids.

It is interesting to note that our sample is typically old, whereas the spiral arms feature is most often seen in young stars (classical Cepheids, masers, and UMS stars). When comparing the features in J_R to the spiral arms, we should therefore bear in mind that intrinsic differences might be present because they are two different populations.

Fig. C.3. Probability distribution function of the azimuthal (left panel) and vertical (right panel) velocities for the subsample of giant stars (red bars) and Cepheids (blue bars). The median values (P50) and median absolute deviations (σ) for both samples are indicated in the inset. Vertical dashed lines denote median values.

All Figures

	Fig. 1. Distribution of the median (equivalent to the P50 percentile) of JR on the Galactic Plane (|Zmax|< 0.5 kpc). The solid black circle denotes the solar position. Features labelled A to D are discussed in the text.
In the text

	Fig. 2. Distribution of JR percentiles on the Galactic Plane (|Zmax|< 0.5 kpc). The x percentile is denoted by Px. The solid black circle denotes the solar position.
In the text

	Fig. 3. Density map of the selected sample in the Cartesian plane (X, Y) on the left. Distribution of median errors in JR on the Galactic plane on the right.
In the text

Fig. 4. Comparison of the maps of P50(JR) with the spiral arms reported in literature. First panel: contour lines enclose the overdensities found in the subsample of giants. Second panel: solid lines represent the Scutum (cyan), Sagittarius (yellow), Local (blue), and Perseus (black) spiral arms of Reid et al. (2014), while the dotted lines correspond to their extrapolation in azimuth. Third panel: Solid lines represent the segments of spiral arms of Lemasle et al. (2022), in which their same naming convention is used, while the colorcode results from a visual comparison with these of Reid et al. (2014). The additional structures are indicated by red and pink lines for description convenience. Fourth panel: contour lines illustrate the overdensities reported by Poggio et al. (2021). Background image: reproduction of Fig. 1 using a grey color-scale to increase the contrast between the coloured lines and the background map. The solid white circle denotes the solar position.

In the text

Fig. 5. Azimuthal velocity Vϕ vs. R diagram colorcoded with the median JR. The colorbar has been intentionally set in logarithmic scale to cover a wide range of values in JR. The moving groups (dashed yellow lines) are displayed from the bottom left to the upper right corner as follows: Hercules, Dehnen98-6, Horn-Dehnen98, Hyades, Coma Berenices, Sirius, and Arch1-Hat. Black ellipses enclose the two selected areas (see the text), while the Sun is denoted by the solid black circle.

In the text

	Fig. B.1. Example of an orbit in the (R, Z) plane (blue curve) with the lines of constant u (dot-dashed ellipses) and v (dashed hyperbolas) in the background. The units of the axis are arbitrary.
In the text

Fig. C.1. Distribution of stars in the Hertzsprung–Russell diagram for the Galactic disc sample (|Zmax|< 0.5 kpc) shown in the upper panel. The dashed black line represents the boundary condition considered to separate Main Sequence and RGB stars (Eq. C.1) neglecting extinction. Density map of the giant sample in the (X, Y) plane in the bottom panel. The solar position is denoted by a solid black circle.

In the text

	Fig. C.2. Distribution of JR percentiles for the giant subsample. Details are the same as those in Fig. 2. Upper-left panel corresponds to the median (equivalent to Fig. 1), while two additional percentiles are shown in the central and right upper panels.
In the text

	Fig. C.3. Probability distribution function of the azimuthal (left panel) and vertical (right panel) velocities for the subsample of giant stars (red bars) and Cepheids (blue bars). The median values (P50) and median absolute deviations (σ) for both samples are indicated in the inset. Vertical dashed lines denote median values.
In the text