© The Authors 2025

1 Introduction

Understanding our home galaxy is one of the greatest challenges in modern astrophysics. In the Galactic disc, where most of the stars are contained, a great variety of physical mechanisms are expected to act in concert. For instance, it has been known for a long time that the disc contains large-scale non-axisymmetric features, including a central bar (Okuda et al. 1977), spiral arms (Georgelin & Georgelin 1976), and a warp (Burke 1957; Kerr 1957) in the outer parts. These features can influence the evolution of the Galactic disc and leave their footprints in the observable properties of the stars. Additionally, vertical and inplane disturbances in the Galactic disc can be induced by the satellite galaxies that orbit the Milky Way. The co-existence and interplay of both internal (e.g. the influence of the bar and spiral arms) and external (e.g. accretion events) physical mechanisms determine how the entire Galactic system evolves over time and, consequently, the properties of the stars that can be observed today (e.g. their spatial distribution, their motions, the amount of metals in their atmospheres, etc.).

In the era of large astronomical surveys, the physical mechanisms at work in our Galaxy can be investigated in great detail. The unparalleled wealth of data from the Gaia satellite (Gaia Collaboration 2016; Gaia Collaboration 2018a, 2023c) and ground-based surveys revealed the richness and complexity of the stellar motions in the Galactic disc, triggering new questions and opening new scenarios. Using astrometry from the first Gaia Data Release, the kinematic signature of the warp was revealed by Schönrich & Dehnen (2018) as a monotonic increase in V_Z as a function of angular momentum L_Z, on top of which an additional ‘wave-like pattern’ was also noted (See also Huang et al. 2018). Using line-of-sight velocities V_los ¹ from Gaia Data Release 2 (hereafter, Gaia DR2), Friske & Schönrich (2019) and Khanna et al. (2019) both noted large scale undulations in V_R as a function of angular momentum, L_Z, or energy. Poggio et al. (2018) photometrically selected both giants and upper main sequence stars from Gaia DR2 with the help of 2MASS photometry, mapped over a large extent of the Galactic disc the kinematic signature of the warp in the vertical velocities, V_Z. Romero-Gómez et al. (2019) found a high degree of complexity in the vertical distribution and velocities of OB and red giant branch (RGB) stars, triggering the need for complex kinematic models flexible enough to combine both wave-like patterns and an S-shaped lopsided warp. Not long after, it was shown that the astrometric data for the giants could be modelled with a precessing warp (Poggio et al. 2020; Cheng et al. 2020). Chrobáková & López-Corredoira (2021) also tried to estimate the warp precession rate, but due to their large uncertainties, their result was both consistent with a precessing and a static warp (which they favoured, due to Occam’s razor). Using accurate distances and kinematics of classical Cepheids, Dehnen et al. (2023) and Cabrera-Gadea et al. (2024) mapped the Galactic disc out to ~15 kpc of the Sun, finding not only that the warp is precessing, but also that its precession rate gradually declines with the galactocentric radius. Recently, using stars with measured line-of-sight velocities in Gaia Data Release 3 (hereafter DR3), Jónsson & McMillan (2024) also found that the warp appears to be rapidly precessing, but additional kinematic disturbances would be needed to match the data.

In the vicinity of the Sun, a Galactic North-South asymmetry was discovered by Widrow et al. (2012) by analysing the number density and bulk velocity of main sequence stars from the SDSS-DR8 (Sloan Digital Sky Survey Data Release 8, Aihara et al. 2011) and SEGUE (Sloan Extension for Galactic Understanding and Exploration, Yanny et al. 2009). Not long after, Yanny & Gardner (2013) confirmed and expanded the discovery of Widrow et al. (2012), detecting a significant Galactic NorthSouth asymmetry in the number density of K and M dwarf stars from the SDSS-DR9 (Ahn et al. 2012). Using LAMOST spectroscopic velocities (Cui et al. 2012; Zhao et al. 2012) and PPMXL proper motions (Roeser et al. 2010), the study of Carlin et al. (2013) reported sub-structures in bulk velocities of disc stars near the Galactic anti-centre. Williams et al. (2013) analysed the vertical velocity field of red clump giant stars in RAVE (Siebert et al. 2011) and detected a rarefaction and compression pattern, suggestive of density wave-like behaviour.

The hypothesis that a passing satellite or dark matter subhalo can excite coherent oscillations of the Galactic stellar disc has been extensively studied in the literature. It has been shown that encounters with satellite galaxies can excite bending and breathing modes in the Galactic disc (Toomre & Toomre 1972; Mathur 1990; Weinberg 1991). Using numerical simulations, Minchev et al. (2009) showed that an initial energy kick approximating a massive Galactic merger can induce waves in the stellar velocity distribution of an axisymmetric galactic disc. Gómez et al. (2013) found that a satellite similar to the Sagittarius dwarf galaxy (Sgr) can produce north-south asymmetries and vertical wave-like behaviour in the Galactic disc that are in qualitatively agreement with the observational results from Widrow et al. (2012). Toy-model calculations and simulations of disc-satellite interactions from Widrow et al. (2014) showed that the response of the disc depends on the relative velocity of the satellite. Using cosmological simulations, Gómez et al. (2017) detected ‘integral sign’ warps and vertical waves in the stellar discs of Milky Way-sized galaxies.

The Gaia phase spiral, discovered by Antoja et al. (2018) and subsequently studied by several works (e.g. Binney & Schönrich 2018; Bland-Hawthorn et al. 2019; Laporte et al. 2019; Bland-Hawthorn & Tepper-García 2021; Khoperskov et al. 2019; Hunt et al. 2022; Grand et al. 2023; Alinder et al. 2024; Frankel et al. 2023, 2024; Tremaine et al. 2023) suggests that the Galactic disc is somewhat out-of-equilibrium and is still recovering from a recent perturbation. Several studies indicate that the phase spiral can be generated by a perturbation induced by the passage of the Sagittarius dwarf galaxy (Binney & Schönrich 2018; Bland-Hawthorn et al. 2019; Bland-Hawthorn & Tepper-García 2021) and this is currently the favoured explanation overall. However, Bennett et al. (2022) were unable to reproduce the amplitude and wavelength of the vertical perturbation observed in the Milky Way disc testing a range of possible models for Sgr, triggering the need for a more complex solution. Mapping the phase spiral with Gaia DR3, Hunt et al. (2022) discovered a transition from one-armed ‘bending mode’ spirals in the Solar neighbourhood to two-armed ‘breathing spirals’ in the inner Galaxy, containing signatures of multiple perturbations.

In the anti-centre region, an apparent ring of stars was discovered at low galactic latitudes at 18 kpc from the Galactic Centre (GC, Newberg et al. 2002; Yanny et al. 2003). This feature has been called the Monoceros ring, Monoceros stream, or Galactic anti-centre stellar stream (Crane et al. 2003; Rocha-Pinto et al. 2003). Using data from SDSS between Galactic longitudes 110° < l < 229°, Xu et al. (2015) detected an asymmetry in the number counts of main-sequence stars, which they explained as roughly concentric disc oscillations, opening in the direction of the Milky Way’s spiral arms. On the other hand, using the Pan-STARRS1 survey, Morganson et al. (2016) mapped the three-dimensional (3D) structure of the Monoceros ring, finding that it is composed of two roughly concentric arcs, which do not appear to align with the Milky Way spiral arms. Even further from the GC than the Monoceros ring, additional sub-structures have been detected, such as the Triangulum Andromeda stream (Majewski et al. 2004; Rocha-Pinto et al. 2003; Martin et al. 2007; Price-Whelan et al. 2015) or the Pisces globular cluster stream (Bonaca et al. 2012; Martin et al. 2014). Vertical corrugations have also been observed in external galaxies (e.g. Matthews & Uson 2008a,b).

Using high-resolution N-body simulations, Bland-Hawthorn & Tepper-García (2021) found that a satellite galaxy similar to Sgr can generate a corrugated bending wave. Tepper-García et al. (2022) found that the satellite-induced corrugation is present both in the stars and gas; the corrugation appears to be initially in phase, the two components move apart after a few rotation periods (500-700 Myr).

In this contribution, we explore the vertical structure and kinematics of the Galactic disc, looking for vertical perturbations on top of the warp of the Milky Way. We used two samples of young stellar populations: (i) a sample of young giant stars and (ii) the Cepheids catalogue by Skowron et al. (2024, S24 hereafter). Specific details on our datasets and on their selection are given in Section 2, as well as a description of the derivation of their distances. In Section 3, we give an overview of the young giant sample. For an overview of the Cepheids sample, we refer to the catalogue featured in S24, as well as a recent application of their data to map the large-scale spiral structure of the Milky Way (Drimmel et al. 2024). We analyse the vertical distribution of our two samples in Section 4 and their kinematics in Section 5. We discuss the obtained results in Section 6 and present our summary and conclusions in Section 7.

2 Data

This work is based on two different datasets, which are described below.

Young giant sample². This dataset has been selected from the Gaia DR3 catalogue (Gaia Collaboration 2023c) with the specific objective to study young populations in the Galactic disc using both spatial and kinematic information. As illustrated in detail below (see Section 2.1), we wanted to select stars with typical effective temperatures of T_eff ~ 4500-5000 K and surface gravities log(g) of about ~1 dex. This corresponds to the blue loop evolutionary stage (see the massive sample in Gaia Collaboration 2023b) and, specifically, to its cold phases (see sample A in Poggio et al. 2022). These tracers have been found to map the segments of the nearest spiral arms in the Galaxy, as expected for young populations in the Galactic disc. Indeed, a comparison between theoretical isochrones in the range of typical disc metallicities indicates that this portion of the T_eff-log(g) diagram is expected to be populated by young stars; that is, approximately less than 100 Myr old (see Figure 1, more details in Section 2.1). Those stars are also expected to be typically bright (with absolute magnitudes between −3 and −7 magnitudes in the G band), making them visible out to relatively large distances. To select our sample, we first perform a preliminary selection based on the astrophysical parameters from Andrae et al. (2023), as described in Section 2.1. Then, we inferred Bayesian distances to each star using a prior specifically constructed to be self-consistent with the adopted selection criteria and compared the obtained distances with other values found in the literature (Section 2.2). Finally, we refined our selection with additional cleaning criteria (Section 2.3). We end up with a catalogue of 17 030 young giant stars, which sample the Galactic disc out to approximately 7 kpc in heliocentric distance. Gaia DR3 line-of-sight velocities are available for 14 219 young giant stars, which is the great majority of the sample (≃83%).

Cepheid sample. This dataset is based on the classical Cepheids catalogue recently presented by S24, which contains new distances based on mid-infrared photometry for 3425 stars, out to approximately 15 kpc in terms of heliocentric distance. From their catalogue, we use the ‘d_av’ and ‘dd_av’ columns as the heliocentric distance and uncertainty, respectively. Typically, the distance uncertainties are less than 5%, however, S24 suggested the uncertainties could be up to 13% at most. Following Drimmel et al. (2024), we removed the stars whose distances are clearly inconsistent with their astrometry, applying a cut Q<5 (where the Q parameter is defined in Equation 14 of S24). After applying this cut, we obtained 3357 stars, of which 2034 also have available line-of-sight velocities and proper motions from the Gaia DR3 catalogue, allowing us to calculate the full 3D kinematic information for each star. For additional details on the Cepheids sample, we refer to Skowron et al. (2024) and Drimmel et al. (2024).

Our two datasets are clearly complementary, that is, the Cepheid sample reaches very large distances, but with only three thousand stars; on the other hand, the young giant sample contains approximately five times the number of stars in a relatively smaller region (out to a distance of ~6-7 kpc) and, thus, with a denser coverage. It should be noted, however, that our two samples are both young, and therefore are expected to trace similar sub-structures, providing two different views of the young Galactic disc, from independent sub-samples of the same young stellar population.

Fig. 1 Comparison between the distribution of XGBoost stars (Andrae et al. 2023) in the Kiel diagram and the PARSEC isochrones for different metallicities. The portion of the Kiel diagram selected in this work is shown by a dark blue dashed line.

2.1 Preliminary selection of the young giant sample

At the time of writing, the largest publicly available catalogue of astrophysical parameters comes from Andrae et al. (2023), who derived the stellar parameters (T_eff,log(g),[M/H]) using Gaia DR3 XP spectra (De Angeli et al. 2023; Montegriffo et al. 2023), combined with CatWISE photometry (Marocco et al. 2021), for a total of ~175 million stars. We refer to it hereafter as the XGBoost sample, after the algorithm they employed. To select young giants in Gaia Data Release 3, we selected sources that fall in the region of the Kiel diagram with: $$\begin{array}{r}2<\text{log}(g)/\text{dex}<0.5,\\ \text{log}(g)>(C\cdot \text{T}+I_L),\\ \text{log}(g)<(C\cdot \text{T}+I_R),\end{array}$$(1)

where the coefficient C=0.00192 dex K⁻¹ approximately follows the natural slope of the RGB branch (using the inclination adopted in Gaia Collaboration 2023b, to select their massive sample, but with a different coverage of the Kiel diagram), while the intercepts I_L = −8.3 dex and I_R = −7.3 determine the position of two inclined lines, respectively, at the left and right borders of the selected area in Figure 1. The cut in log(g) helps us to select part of the blue loop stars, keeping only values lower than the red clump, which contains less massive and typically older stars. This inclined portion of the Kiel diagram follows the selection criteria applied in our previous work (Poggio et al. 2022), where young giants were selected based on Gaia GSP-spec data (Recio-Blanco et al. 2023). In this work, we aim to perform a similar selection, but going to fainter magnitudes, to map a larger portion of the Galactic disc. From the above-selected objects, we removed stars having apparent magnitude G > 15, where we found that the sample to be incomplete. Finally, we removed possible globular cluster members using the catalogue of Hunt & Reffert (2023). After these cuts, we obtained a total of 138 829 stars.

2.2 Distances of the young giants sample

For each star selected in Section 2.1, we infer the probability density function (PDF) of the heliocentric distance d for the star’s galactic coordinates (l, b) via Bayes’ theorem, P(d | l, b, ϖ_corr, σ_ϖ) ∝ P(ϖ_corr | d, σ_ϖ)P(d | l, b), where ϖ_corr is the observed parallax corrected for the parallax zero-point offset (see below) and σ_ϖ is the parallax uncertainty. We assume a Gaussian likelihood P(ϖ | d, σϖ), and a prior similar to Astraatmadja & Bailer-Jones (2016), as per their Equation (7), the Milky Way prior: $$P(d|l,b)\propto \rho (l,b,d)\mathcal{F}(l,b,d),$$(2)

where ρ(l, b, d) is the density distribution of stars in the Galactic disc and F(l, b, d) is the fraction of observable objects (based on the selection function of the catalogue, as explained below). Here, we assume a simple model for the disc density ρ(l, b, d), composed of an exponential disc with a galactocentric radius, R, with radial scale length of LR = 2.6 kpc (Bland-Hawthorn & Gerhard 2016) and vertical height, Z, with a vertical scale height of h_z = 150 pc, following Poggio et al. (2018). The assumed distance of the Sun from the GC is R = 8.277 kpc GRAVITY Collaboration (2022). It is important to note that the simple density prior assumed here does not include non-axisymmetric features (e.g. warp, spiral arms, etc.). This is because non-axisymmetric structures represent one of the main features of interest of this study and we want to avoid introducing artefacts due to our assumptions. Therefore, any deviation from symmetry observed in the 3D distribution of stars is not expected to be attributable to the assumed prior. The impact of various assumptions incorporated in our prior is discussed in the following sections.

The term F(l, b, d) is calculated following the approach of Astraatmadja & Bailer-Jones (2016) and Poggio et al. (2018), but now re-adapted to the specific case of the young giants selected in this work. First, we constructed the luminosity function Φ(M_G) in the G band by applying the same cuts specified in Section 2.1 to the Kiel diagram based on the PARSEC isochrones (Bressan et al. 2012; Chen et al. 2014, 2015; Pastorelli et al. 2019), assuming the canonical two-part power law IMF corrected for unresolved binaries Kroupa (2001, 2002), a star formation rate constant with time (which might be realistic given the typical age of the target stars of ~ 100 Myr), and solar metallicity.

The fraction of observable objects, F, can then be estimated by integrating the luminosity function for all possible distances along the line-of-sight (after including the completeness of the sample at each apparent magnitude G; see details below), down to the faintest observable apparent magnitude (see Equations (4)–(6) in Astraatmadja & Bailer-Jones 2016), after taking extinction into account. Following the procedure laid out in Khanna et al. (2024), we estimate the extinction by combining publicly available dust maps (Schlegel et al. 1998; Green et al. 2019; Lallement et al. 2019, etc.). Then, we have $$\mathcal{F}(l,b,d)={\int }_{G_{min}}^{G_{max}}S(l,b,G)\mathrm{\Phi }(M_G)dG$$(3)

where S (l, b, G) is the selection function of our selected sample and M_G = G - 5 log(d) + 5 + A_G(l, b, d) for d in parsecs. To estimate S(l, b, G), we follow the approach laid out in Castro-Ginard et al. (2023), computing the ratio of stars that end up in the catalogue of Andrae et al. (2023), to the entire Gaia catalogue (n) in bins of (l, b, G), between 4 < G < 15. On the Gaia archive (see Appendix B), we computed this ratio in bins of G (1 mag wide) and sky position (l, b) tiled at HEALPix level 5 (Górski et al. 2005). Using these number counts, we then take $S(l,b,G)=\frac{k+1}{n+2}$ as the mean probability for a source to end up in our sub-sample (with k stars) as proposed by Castro-Ginard et al. (2023). The completeness at the edge of the brightness limit (14 < G < 15) for the young giants is shown in Figure B.1.

After cleaning the young giant sample (see Section 2.3 below), we found that 593 Cepheids are also included in this sample. This is a consequence of the instability strip partially overlying the area of the Kiel diagram used to define the initial selection of these stars. We took advantage of this overlap by using these common stars to check their derived distances against their photometric Cepheid distances, d_S24, taken from Skowron et al. (2024). The first step of this process is to make a direct comparison of the measured Gaia parallaxes with the photometric parallax (1/d_S24). Figure 2 shows the differences between the measured and photometric parallax with no zeropoint correction to the measured parallaxes (black points) and after applying the zero-point correction of Lindegren et al. (2021) (red points). The corresponding running medians (black and red curves) show substantial improvement after the zeropoint correction is applied; however, a mean offset of 7 μas is still present, which is nonetheless in agreement with the offset found by S24. Such an offset will have an impact on the inferred distances for distant objects (while it is totally irrelevant for nearby stars). Therefore, we can apply it on top of Lindegren et al. (2021)’s predicted zero-point offset zp_L21, that is, ϖ_corr = ϖ - zp_L21 - 0.007 mas, where ϖ is the parallax from the Gaia archive. The blue points in Fig. 2 (left panel) show the results using the final corrected parallaxes, with a running median of the final parallax difference (blue curve) tending to coincide with 0. The right panel of Figure 2 shows the comparison between the distances calculated as prescribed above, with the photometric distances, d_S24.

To test the reliability of the distances obtained in this Section, we also cross-matched our young giants sample with other catalogues available in the literature (e.g. Hunt & Reffert 2023; Bailer-Jones et al. 2021). Within 8 kpc in heliocentric distance (where our sample has the best coverage), 93% or more of the stars deviate less than 20% from the distance estimates taken from the literature. We also note that 80% of the young giant stars have measured parallaxes, ϖ, with a relative uncertainty of less than 20% (i.e. satisfying ϖ/σ_ϖ > 5). Finally, ~90% of the stars in our catalogue have a relative distance error σ_d/d < 20%, where σ_d is the average between σ_up = d₈₄ - d₅₀ and σ_down = d₅₀ - d₁₆, where d₈₄, d₅₀, d₁₆ are, respectively, the 84^th, 50^th and 16^th percentiles of the posterior distribution of the distance for each single star. More details are given in Appendix A.

2.3 Cleaning of the young giant sample

The selection described in Section 2.1, which is based on T_eff-log(g) parameters, is not expected to lead to a perfectly pure sample of young giant stars. This is because, in this region of the Kiel diagram, the isochrones of a given age tend to move to the left (i.e. larger T_eff values) for lower values of metallicity, as shown in Fig. 1. This is a consequence of the decreasing opacity of the atmosphere for stars of lower metallicity. Therefore, it is possible that our selected box in the Kiel diagram does not only contain young stars of disc-like metallicity, but also some metalpoor contaminants of significantly larger age. Additionally, we note that random errors on T_eff-log(g) can display a fraction of the red giant branch (RGB) stars as potential contaminants that fall into the selected region. Andrae et al. (2023) did not provide individual uncertainties for the single sources, but they reported the median uncertainties of the sample, namely, ~50 K in T_eff and 0.08 dex in log(g), respectively.

To clean our sample, we removed stars whose metallicities were lower than [M/H]_lim = −0.04 · R - 1.25, where R is the galactocentric radius and −0.04 dex/kpc follows the typical values of the metallicity gradient of the Galactic disc (e.g. Gaia Collaboration 2023a). The intercept of our [M/H]_lim cut has been tuned to take into account the metallicity dispersion at a constant R and, at the same time, to try to minimise the contamination from metal-poor stars (see Figure A.2). In the inner parts of the Galaxy, we see a significant contribution from metalpoor stars with high velocity dispersions; therefore, we applied a cut at R < 5 kpc. Since the above-mentioned metallicity cut is very broad, we refined it by removing the low-angular momentum stars with V_φ < 180 km/s (assuming, for the solar velocity, the values reported in Drimmel & Poggio 2018), for the stars that have available radial velocities (>90% of the sample). Inside R < R, we remove the stars that have both [M/H] < −1 dex and |z| > 0.5 kpc. We also applied a cut of |b| < 30°, which removes most of the Magellanic Clouds. To more efficiently remove stars in the Large Margellanic region on the sky, we also cut stars with b < −28^° between galactic longitudes 275^° and 285^°. As a test, we also tried to remove the open cluster members from Hunt & Reffert (2023) to ensure that our maps (e.g. stellar velocity fields, see the sections below) were typically mapping field stars and not dominated by individual clusters. Since we obtained very similar results by including or not open cluster members, we decided to keep them to have larger statistics. We also removed six stars that are expected to be members of the globular cluster NGC 6656, as explained in Appendix A. The apparent magnitude distribution of the final catalogue is shown in Figure 3.

Fig. 2 Left panel: comparison between the photometric parallaxes (obtained from Cepheids distance estimates from Skowron et al. 2024) and trigonometric parallaxes without zero-point correction (black), with Lindegren et al. (2021) correction (red), and Lindegren et al. (2021) + additional offset correction (blue). See main text for a detailed explanation. Right panel: derived distances for the 593 young giants (dYG), which were also classified as Cepheids, compared to their photometric (Leavitt Law) based distances (dS 24).

Fig. 3 Apparent G magnitude distribution of the young giant sample.

3 Overview of the young giant sample

In this section, we analyse the 3D distribution of the samples presented in Section 2. A first overview of the young giants sample is given in Figure 4, where the distribution of stars in the vertical direction Z as a function of galactocentric radius, R, is shown (for stars at all azimuthal angles in the Galactic disc). As we can see, the stars of our sample are located close to the Galactic midplane, as expected from young disc stellar populations. However, we note that the distribution of stars in the disc is not continuous: clumps in the stellar density are apparent. Such “blobs” are vertical sections of the Galactic spiral arms. Indeed, the inplane distribution of the stars shows segments of the nearest spiral arms in the Galaxy, as shown in Figure 5. This is a further indication that our selected stars are expected to be typically young. The position and orientation of the observed segments are in agreement with the spiral structure found by previous works (Zari et al. 2021; Poggio et al. 2021; Gaia Collaboration 2023a). A more detailed analysis of the in-plane distribution of this sample will be presented in a separate paper.

In the outer regions of the disc (R ≳ 10 kpc), we can observe a displacement of the Galactic mid-plane with respect to Z=0 kpc, together with an increase of the vertical dispersion as a function of R. This is a 2D projection of the Galactic warp and flare, which both act on the outer disc. To get a 3D view of the vertical displacement of the Galactic disc, we can have a look at Figure 6, which is similar to Fig. 4, but split in different overlapping slices of Galactic azimuth. (In this work, we define the Galactic azimuth φ as increasing in the direction of Galactic rotation, and equal to 0° toward the Galactic anti-centre.) Figure 6 shows that the outer disc bends downwards for φ ≲ 0°, and in the opposite direction for positive φ ≳ 0°. This large-scale distortion of the outer Galactic disc is the well-known Galactic warp, which was first detected in the HI data (Kerr 1957; Burke 1957; Levine et al. 2006) and later confirmed in different stellar populations in numerous previous works (e.g. López-Corredoira et al. 2002; Reylé et al. 2009; Amôres et al. 2017; Uppal et al. 2024, and others), including classical Cepheids (Skowron et al. 2019a,b; Chen et al. 2019; Lemasle et al. 2022; Dehnen et al. 2023; Cabrera-Gadea et al. 2024). For our specific sample of Cepheids, the presence of the warp has been illustrated in Fig. 19 (left panel) of Skowron et al. (2024).

In the following sections, we describe how we inferred the large-scale 3D shape of the warp using our two samples of young giants and Cepheids (Section 4.1). Next, we examine the residuals of the inference (Section 4.2). As shown below, a classical description of the warp is not sufficient to explain the observed vertical distribution of stars.

Fig. 4 Vertical section of the Galactic disc as a function of galacto-centric radius, R. Over-dense regions are the spiral arms of the Galaxy seen edge-on. Stars are colour-coded by the probability density function at each bin, calculated as the number of stars in each bin over the total number of stars per bin area, where the bin area is 60 pc (in R) times 30 pc (in Z).

Fig. 5 Distribution of the young giant sample in the Galactic plane. The Sun’s position is shown by the symbol. The GC is to the right, at (X,Y)=(R, 0 kpc), with R = 8.277 kpc (GRAVITY Collaboration 2022). Galactic rotation is clockwise.

4 Analysis of the vertical distribution

4.1 Inference of the warp shape

In this section, we infer the shape of the Galactic warp using a classical m = 1 description (similar to Skowron et al. 2019a; Chen et al. 2019, and numerous other works in literature). For a given dataset, we performed the inference using different approaches, to test the robustness of our results. The same procedures were applied to our two samples of young giants and Cepheids, followed by a comparison of the results.

We started with a classical ‘tilted rings’ approach, which we first applied to our young giants sample. We split our dataset into non-overlapping galactocentric annulii of 0.5 kpc in width. Each annulus typically contains between 900 and 1200 stars. As a test, annulii of different width (ranging from 0.1 kpc to 1.5 kpc) were also considered, to check that the final results of this paper do not depend on the chosen bin size. Stars inside R = R were not considered here, given that the amplitude of the warp has been found to be very small or even null inside that radius by several works in the literature. Young giant stars with R > 14.5 kpc were ignored, due to the low number statistics in those regions. For the i-th ring, we fit a simple m = 1 distortion as $$Z_{w,i}=h_{w,i}\sin (\varphi -{\varphi }_{LON,i}),$$(4)

where h_w,i and φ_LON,i are, respectively, the amplitude and phase angle of the line-of-nodes at the radial position covered by the i-th ring. For each star, uncertainties in the Z-coordinate were computed by converting the entire distribution in heliocentric distance (using 1000 mock resamples for each star), and taken as the 16th and 84th percentiles of the posterior distribution described in Section 2.2, to a distribution in Z, assuming no errors in the Galactic latitude, b. For each ring, we performed the inference using the emcee code (Foreman-Mackey et al. 2013), using a broad uninformative prior for h_w,i, and imposing that the phase angle φ_LON,i must be within the range [−180°, 180°]. The likelihood function used to perform the inference is the product of individual Gaussian likelihoods (i.e. sum of log-likelihoods) obtained for each star, assuming a random Gaussian distribution around the model prediction Z_w,i from Equation (4), namey, $$\mathcal{L}=\frac{1}{\sqrt{2\pi }{\sigma }_z}\exp [-{\mathrm{\Delta }}_Z^2/(2{\sigma }_z^2)],$$(5)

where Δ_Z = Z_w,i - Z_i, and Z_i is the individual vertical coordinate of each star. To describe the total dispersion in the Z-coordinate σ_z, individual uncertainties in Z are summed in quadrature to the intrinsic thickness of the (young) stellar disc, using 150 pc as a reference value (broadly consistent with the vertical distribution of the sample), although other values have also been tested (see below). The 16th and 84th percentiles of the obtained marginalised posterior distributions on h_w,i and φ_LON,i for each ring are shown as brown error bars in Figure 7, while the bestfit results for each ring are shown as brown points, for the warp amplitude (left panel) and phase angle (right panel) as a function of galactocentric radius, R. As pointed out by several works (e.g. Dehnen et al. 2023; Cabrera-Gadea et al. 2024), the tilted ring approach has the advantage that it does not assume a parametric functional form of the warp parameters (amplitude, shape of the line-of-nodes) as a function for the galactocentric radius.

As we can see, the warp amplitude systematically increases with R. This is expected from previous observations of the largescale warp. The phase angle, however, varies abruptly from one ring to another inside approximately R = 11 kpc, showing large uncertainties. Such a behaviour is relatively anticipated from geometrical considerations: the warp amplitude is very small at these radii (see discussion in Dehnen et al. 2023), making the measurement of the line-of-nodes more challenging. We notice, however, that there is a net increase of the phase angle in the direction of Galactic rotation for R ≳ 12 kpc, in agreement with Dehnen et al. (2023) and Cabrera-Gadea et al. (2024).

A complementary approach is to infer the warp shape using the individual sources without binning in R, but assuming a parametric shape for the amplitude and phase angle as a function of galactocentric radius. We refer to this approach as the ‘global’ parametric fit. For the amplitude, we assume the typical functional form of $$h_w(R)=h_{w,0}(R-R_w{)}^{{\alpha }_w}$$(6)

This expression is valid for R > R_w, and zero otherwise. We first considered a scenario with a curved (twisted) line-of-nodes (LON), using a functional form that follows the Briggs rule for warps in external galaxies: $${\varphi }_{LON}(R)=\{\begin{array}{ll}{\varphi }_{LON,0}& \text{for R <= R\_T},\\ {\varphi }_{LON,0}+{\beta }_w(R-R_T)& \text{for R > R\_T},\end{array}$$

according to which the phase angle is constant in the inner parts of the disc, and curves outside a certain radius, R_T (which ism by definition, greater than the starting radius of the warp R_W). The parameter β_w defines the rate at which the line-of-nodes changes as a function of R (for R > R_T). The final modelled distortion of the Galactic disc along the Z-coordinate is $$Z_w(R,\varphi )=h_w(R)\sin (\varphi -{\varphi }_{LON}(R))\cdot ,$$(7)

which includes the amplitude h_w (R) and the phase angle φ_LON described above. The results of the statistical inference are presented in Table 1. Systematic uncertainties on each parameter were estimated by testing five different values for the thickness of the disc, based on the typical thickness of young stellar populations, namely 100, 150, 200, 250, and 300 pc. The range of the derived warp parameters assuming different values of disc thickness is adopted as the systematic uncertainty on each specific parameter.

We note that based on this approach, the starting radius of the warp, R_W, was treated as a free parameter, and we obtained a very low value of R_w. We should, however, avoid over-interpreting such a result: it is clearly evident from Figure 7 (left panel) that the warp amplitude at R=9 kpc is smaller than the typical vertical scale-height of the sample. Different mathematical prescriptions with different values of R_w can lead to very similar warp amplitudes in the inner disc (e.g. R ≲ 8-10 kpc), whose differences are presumably too small to be detected, taking into account the intrinsic thickness of the Galactic disc. We note that, for the young giants sample, the tilted rings method and the parametric global fit give very similar results, as shown in Figure 7.

We note that if instead a straight LON (i.e. with a constant phase angle) is assumed, so that Equation (7) becomes φ_LON = Φ_LON,0 at all radii, the resulting warp amplitude of the young giants is similar to the twisted LON scenario, but now with a phase angle φ_LON = 9.9° ± 0.6°.

The analysis performed with the young giants catalogue was also repeated with our Cepheids sample. Figure 7 shows the results obtained with the m=1 non-overlapping tilted rings method (blue circles). With this sample, due to the relatively low-number statistics, we increased the width of the rings to 0.6 kpc in terms of the galactocentric radius. Each ring typically contains between 100 and 230 stars. Figure 7 shows the results for R > 9 kpc, where the warp signal is more apparent. We do not consider stars outside R ≳ 16.5 kpc and/or on the other side of the Galaxy, because the number of objects is too small to give reliable results. We note that this sample reaches much larger distances than the young giant sample. For this sample, we also apply the parametric global fit with a twisted LON, obtaining the results presented in Table 2. For the R_w parameter, the same considerations made for the young giant sample are also valid for the Cepheid sample. We also test a straight LON model, obtaining a phase angle φ_LON = 14.0° ± 0.6°.

The results obtained for the warp structure with our two samples are in good agreement with each other. This is somehow expected, because they are both tracing the young stellar populations in the Galactic disc. This is also in very good agreement with previous works based on young stars. For instance, Dehnen et al. (2023) find an inclination angle of 3° at R ≥ 14 kpc, which corresponds to an amplitude of 0.7 kpc at R = 14 kpc, in perfect agreement with the left panel of Figure 7. Moreover, we note that the general behaviour of the line-of-nodes in our two samples is in agreement with the results obtained in the literature with classical Cepheids (Dehnen et al. 2023; Cabrera-Gadea et al. 2024).

Fig. 6 Same as Figure 4, but now in overlapping bins in Galactic azimuth φ. Azimuthal ranges are chosen to cover the regions where our sample has more data. Isodensity contours are show as black solid lines.

Fig. 7 Inferred amplitude (left panel) and phase angle (right panel) for the m = 1 classical warp model using the young giants sample (orange/brown markers) and the Cepheids sample (blue/light blue markers) based on different approaches. Individual points show the results of the m=1 tilted ring model, while solid lines correspond to the parametric fit, performed on a star-by-star basis, of the twisted line-of-nodes for the young giant (in brown) and Cepheid (in blue) samples. Orange (light blue) transparent lines represent 5 000 random realisations of the MCMC chain for the fits of the young giant (Cepheid) sets.

4.2 Analysis of the residuals

It is natural to wonder whether the classical m = 1 description adopted in the previous Section is sufficient to explain the 3D distribution observed in our datasets. Several works in the literature have already pointed out that a simple m = 1 model is not enough to explain the warp in our Galaxy, suggesting that it is lopsided (Romero-Gómez et al. 2019), or even fitting a Fourier decomposition with m = 0, 1, 2 terms (Levine et al. 2006; Cabrera-Gadea et al. 2024). The approach presented here is complementary to previous works: we fit a simple m = 1 term, and then investigate the residuals, with the aim of determining which additional ingredient is needed to explain the observed vertical structure of the young Galactic disc.

Figure 8 shows the comparison between the median Z in the Galactic plane observed in the young giant dataset (upper left panel), the prediction from the tilted ring model (left panel, middle row), straight line-of-nodes (LON) model (middle panel, middle row) and the twisted LON model (right panel, middle row) described in Section 4 (middle panel), and the corresponding residuals (lower panels). As we can see, the residuals of the fit show a clear coherent pattern rather than random noise.

The pattern is similar for all three m=1 models considered here. Specifically, we can identify a specific region of the outer Galactic disc where the distribution of stars is systematically shifted toward positive Z compared to the global best-fit model (i.e. red region labelled as Feature 1), extending from approximately (X, Y) ≃ (−4 kpc, 4 kpc) to (X, Y) ≃ (−2 kpc, −4 kpc). The morphology of the observed feature is consistent with a nearly radial ripple or corrugation, which cannot be explained by a simple m = 1 warp. The observed corrugation reaches about 150-200 pc in vertical height, and extends over the entire volume covered by the young giant sample for about 10 kpc (corresponding to about 50 degrees in Galactic azimuth). On the other hand, no major vertical displacement is detected within the galactocentric radii of 8-9 kpc (see also Xu et al. 2015; Hunt et al. 2022); however, it is possible that another ripple with very small vertical amplitude is present close to the Solar circle (we defer the study of its possible existence to future studies).

In the lower left corner of the residual maps (lower panels), we can see that stars are systematically located lower than the prediction from our models (blue region labelled as Feature 2). It is impossible not to wonder whether the observed blue region is the ‘negative counterpart’ of the red region observed at inner Galactic radii, being both part of the same oscillatory pattern. We note, however, that the blue region is at the border of our map, and therefore we treat it with caution.

To test the robustness of the obtained maps, we re-performed the warp inference using only the sub-set of stars with high quality distances (i.e. based on a good parallax signal-to-noise ratio, ϖ/σ_ϖ > 5), or using other distance estimates than those calculated in Section 2.2 (for instance from Bailer-Jones et al. 2021; Anders et al. 2022, or selecting stars with ϖ/σ_ϖ > 5 and calculating the distance as 1/parallax), always obtaining residuals maps with features similar to those shown in Figure 8.

Figure 9 shows again a comparison between the data (upper left panel), the best-fit predictions from our models (middle panels) and the corresponding residuals (lower panels), similar to Figure 8, but now for our sample of Cepheids. We note that the lower panels in Figure 9 exhibit some similarities and differences (and there are several reasons for this being so, as discussed in the following). The residuals of the tilted ring model (lower left panel) and the twisted LON model (lower right panel) exhibit a region of systematically positive residuals (red area) in the lower parts of the plot (Y ⪅ 4 kpc), for galactocentric radii between 10 and 14 kpc. In the residuals of the straight LON model (lower central panel), the observed feature is even more prominent, with systematically positive residuals that extend from Y ~ 10 kpc to Y ~ −10 kpc in the range of galactocentric radii 10-14 kpc. For all three considered models, the red region detected in the residuals of the Cepheids sample might correspond to the corrugation (labelled as ‘Feature 1’) found with the young giant sample, suggesting a similar behaviour.

To enable a proper comparison between our two samples, the observed residuals are plotted on the same XY-range in Figure 10. As we can see, the red region of systematically positive residuals discussed above is located at similar positions in the Galactic plane for our two samples. To better explore the geometry of the observed residuals, our maps are then compared to spiral arms models available in literature. For the Cepheids sample, we overplot the model of Drimmel et al. (2024) (coloured lines in the right panel of Figure 10), which is derived from the same dataset used here. For both the young giant and Cepheids sample, we overplot the spiral arm model from Reid et al. (2019) (black solid lines), which is based on masers. Finally, for the young giant sample, we also overplot the overdensity contours of the upper main sequence stars from Poggio et al. (2021) (grey shaded areas in the left panel of Figure 10), which we found to be in good agreement with the spiral structure observed in the young giant sample (Figure 5). The general picture emerging from Figure 10 is that the vertical residuals observed in this work do not exhibit an obvious correlation with the different spiral arm geometries considered here: unlike the red region in the vertical residuals, the spiral arms tend to have a trailing geometry (i.e. positive pitch angle), so that they are opening in the direction opposite to Galactic rotation. We note, however, that the outer arm from the model of Reid et al. (2019) (outermost black solid line in Figure 10) has a very small pitch angle (specifically, 3° in the region of overlap with our maps), which implies that that arm is almost a radial feature. While there is a partial overlap in some regions (approximately −2.5 kpc < Y < 5 kpc), the orientation of the red feature in the vertical residuals is different, because it is more similar to a leading (i.e. negative pitch angle) rather than a trailing (i.e. positive pitch angle) spiral. We also compared our maps to the spiral arm model from Hou & Han (2014) (here not shown for brevity), reaching conclusions similar to those obtained with the spiral models discussed above.

It is already evident from the lower panels of Figures 8 and 9 that the corrugation is not exactly radial: in the lower parts of the plot (Y ⪅ −2 kpc) it tends to be within 10 and 12 kpc in terms of the galactocentric radius, R, while it is located between 12 and 14 kpc in the upper parts (Y ⪆ −2 kpc). For visual reference, Figure 11 shows a comparison between our maps and two leading spiral curves with pitch angles of −15° (dashed-dotted black curve) and −7.5° (solid black curve), which cross the radius R=12.3 kpc at Galactic azimuth φ = 0° (toward the Galactic anticentre). The black dotted line shows a curve at constant radius R=12.3 kpc, shown for reference.

Fig. 8 Comparison between the maps of median Z-coordinate from observations, best-fit model and residuals. Top-left: median Z-coordinate of the young giant sample. The position of the Sun is shown by the black cross at (X,Y)=(0 kpc, 0 kpc). Middle: prediction from the m=1 tilted rings model (left), straight line-of-nodes parametric model (centre) and twisted line-of-nodes parametric model (right). Bottom: corresponding vertical residuals for the m=1 tilted rings model (left), straight line-of-nodes parametric model (centre), and twisted line-of-nodes parametric model (right). Dotted lines in the three lower panels show rings at constant radii of R=10, 12, and 14 kpc from the GC. An interactive 3D view of this figure is available online.

Fig. 9 Same as Figure 8, but using the Cepheids sample. The black cross shows the position of the Sun, while the black dot shows the position of the GC. Due to the different spatial coverage, the range in Z is different compared to Figure 8.

Fig. 10 Vertical residuals, ΔZ, of the m=1 straight line-of-nodes model using young giants (left panel) and Cepheids (right panel), compared to spiral arm geometries taken from the literature: Reid et al. (2019) (black solid lines, both panels), Poggio et al. (2021) (grey shaded contours, left panel), and Drimmel et al. (2024) (coloured lines, right panel). For the Drimmel et al. (2024) model, different colours indicate the Perseus arm (grey), Local/Orion arm (orange), Sagittarius-Carina arm (blue), and Scutum arm (red).

5 Kinematics

5.1 Radial motions

To learn more about the residuals observed in the vertical distribution of stars, we now compare them to the stellar velocity field, looking for possible correlations. In the left panel of Figure 12, we can see the same residuals shown in Figure 8 (right panel) for the young giant sample, now overplotted with two dashed lines that roughly delimit the region where the majority of the positive residuals in Z are located (i.e. Feature 1 discussed above), indicated simply for the purpose of comparison with the velocity maps of the young giants (as well as the Cepheid sample).

The middle panel of Figure 12 shows the median radial component of the velocity, V_R, in terms of the galactocentric cylindrical coordinates for the young giant sample. We note that the V_R map shown here does not depend on any warp model. This is because the warp is mainly a vertical distortion. In our coordinate system, positive (negative) values of V_R correspond to velocities directed towards the outer (inner) parts of the Galactic disc. Surprisingly (at least to us), the stars located in Feature 1 (i.e. between the dashed lines) are typically moving outwards, with velocities of about 10-15 km/s. Roughly outside the dashed lines, the typical orientation of the velocities changes in sign, with stars that are typically moving inwards.

A similar pattern is also observed in the median VR of our Cepheids sample, shown in the middle panel of Figure 13. As we can see, also with this sample the positive VR region coincides with the positive crest identified with the young giants (delimited by the dashed lines). Furthermore, the Cepheids sample indicates that such features are even more extended, reaching about 20 kpc in length.

To better quantify the statistical correlation between the Δ_Z and V_R maps, we dissected the XY-plane in large (nonoverlapping) pixels of 1 kpc (2 kpc) width for the young giants (Cepheids) sample. The median Δ_Z and V_R have been calculated for the pixels that contain at least 10 stars from the young giant sample (five stars from the Cepheids catalogue), and only those with R>10.5 kpc (where the above-discussed features are located). To quantify the correlation between these two quantities, we calculated Kendall’s correlation coefficient (Kendall 1938), obtaining values between τ_b = 0.33 and 0.37 for the young giants sample, and between 0.3 and 0.46 for the Cepheids sample, depending on which of the three different m = 1 warp models was considered (with the straight LON model showing larger correlation). In all three cases, and for both samples, the obtained τ_b indicates a strong positive correlation³, suggesting that regions of the disc with systematically positive residuals in Z (i.e. shifted upwards with respect to a simple m = 1 warp model) statistically tend to exhibit positive values of V_R, and vice versa.

The observed statistical correlation might suggest that the quantities Δ_Z and V_R are linked by a common physical origin (see discussion in Section 6). However, correlation doesn’t necessarily imply causation. It is therefore important to consider alternative scenarios, like, for instance, the possible impact of the spiral structure on the in-plane motions. Figure 14 (right panel) shows the comparison between the V_R map and the spiral model from Drimmel et al. (2024). As we can see, the observed V_R features do not exhibit an obvious alignment with this spiral model. Specifically, the observed large-scale positive V_R feature dominating the outer (>R_/odot) disc is crossing the Perseus arm (grey curve) and the Local Arm (orange curve). A similar consideration can be done while comparing the V_R map of the young giant sample with the overdensity contours from Poggio et al. (2021) (which, from left to right, show segments of the Perseus, the Local and the Sagittarius-Carina arm) in the left panel of Figure 14. Again, a large-scale positive V_R feature in the outer disc is connecting the Perseus and the Local Arm. The overdensity in the upper part of the Perseus arm (the Cassiopeia region, at approximately (X,Y)~(−2 kpc, 3 kpc)) systematically moves inwards (negative V_R), showing a different behaviour compared to the rest of the Perseus arm. Similarly, the Local Arm exhibit alternate regions of inwards/outwards motions. On the other hand, the Sagittarius-Carina arm is mostly dominated by inward motions, especially in the lower parts of the plot). The spiral arms from Reid et al. (2019) (solid black lines) have a different orientation with respect to the other previously considered spiral arms geometries. In the outer disc, the positive V_R feature does not exhibit an obvious alignment with the Perseus arm from this model, while there is an overlap with the outer arm (left black curve), though with a slightly different orientation: in the upper part of the plot, the outer arm seems to be located in the portion of the red V_R feature which is closer to the GC, while in the lower parts of the plot the outer arm is at the outer edge (i.e. more distant parts from the GC) of the positive V_R feature.

Fig. 11 Same as Figure 10, but now compared to different curves, plotted for visual reference. The dashed-dotted and solid black lines show a leading spiral curve with pitch angles −15° and −7.5°, respectively (see text for more details). The black dotted line shows a circle (pitch angle 0°) at constant radius R=12.3 kpc. In the right panel, the Perseus arm from Drimmel et al. (2024) (grey line) is also shown for reference.

Fig. 12 Comparison between the residuals in the vertical coordinate ΔZ = Z - ZWARP (left panel), the median radial motions VR (middle panel) and the residuals in the vertical velocity ΔVZ = Vz — VZ,WARP for the young giant sample. Dashed lines give a rough indication of the region where positive residuals ΔZ are obtained with this sample. Dotted lines are the same, but shifted by half the amplitude of the region (see text).

Fig. 13 Same as Figure 12, but now for the Cepheids sample.

5.2 Vertical motions

If the residuals discussed in Section 4.2 are related to a vertical wave in the Galactic disc, we should see a trace in the vertical velocities V_Z of the stars of our samples. However, the vertical velocity field of stars in the Galactic disc is expected to be dominated by the signature of the large-scale warp (e.g. Schönrich & Dehnen 2018; Poggio et al. 2018; Romero-Gómez et al. 2019). Therefore, the natural approach would be to remove the signature expected from the Galactic warp, and then investigate the residuals (similar to the approach performed with the Z-coordinate).

However, modelling the warp signature is not trivial. Several works have found that, in addition to the spatial parameters that describe the 3D shape of the warp, at least one kinematic parameter is needed (Poggio et al. 2020; Cheng et al. 2020; Dehnen et al. 2023; Jónsson & McMillan 2024; Cabrera-Gadea et al. 2024; Zhou et al. 2024). Using classical Cepheids Dehnen et al. (2023) and Cabrera-Gadea et al. (2024) found a prograde precession that decreases as a function of the galactocentric radius, R.

Figure 15 shows the results obtained for the young giants (red points) and Cepheids sample (black points) after dividing our samples in radial bins, and fitting for the warp precession ω_p separately (adopting the warp best-fit models discussed in Section 4). With the young giant sample, the bins are of 1.5 kpc in width, and typically contain between 2500 and 3000 stars. For the Cepheid sample, due to the low number of stars with available line-of-sight measurements, we chose a relatively large bin size of 4 kpc, allowing us to have approximately 500-800 stars for each bin. We infer the warp precession using Equation (9) in Cheng et al. (2020), which also includes the mean radial motions V_R (although we obtain very similar results also if the V_R term is not included). As we can see, our results are in good agreement with other works from the literature, especially for R ⪆ 11 kpc, where the warp amplitude is larger (see previous sections).

After subtracting the prediction from the above-described kinematic model, the residual vertical velocities Δ_Vz are shown in the right panels of Figures 12 and 13, respectively for the young giants and the classical Cepheids. In both plots, we see a feature with systematically positive Δ_Vz. This feature is similar to the one observed in the V_R maps and Δ_Z, but now slightly shifted towards the outer parts of the disc. For comparison, we plot the same lines shown in the previous plots, as well as two additional lines, now shifted by 1.4 kpc (which is the half distance between the two dashed lines), shown as dotted lines. It is interesting to note that the new dotted lines roughly contain systematically positive Δ_VZ.

To test the robustness of the obtained results against different warp kinematic models, we tested different values of ω_p=12, 10, 8 and 3 km/s/kpc, and checked the corresponding Δ_VZ maps. As expected, we found that the Δ_VZ structure associated with the corrugation was always present, though with a different magnitude in Δ_VZ.

Fig. 14 Face-on view of the median radial motions, VR, for the young giant sample (left panel) and the Cepheids sample (right panel), compared to different spiral arm geometries from literature: Reid et al. (2019) (black solid lines, both panels), Poggio et al. (2021) (grey shaded contours, left panel) and Drimmel et al. (2024) (coloured lines, right panel). For the Drimmel et al. (2024) model, different colours indicate the Perseus arm (grey), Local/Orion arm (orange), Sagittarius-Carina arm (blue), and Scutum arm (red).

Fig. 15 Comparison between the obtained values of warp precession with our young giants (red points) and Cepheids sample (black points), together with other values from the literature.

6 Discussion

To facilitate the interpretation of the obtained results, we compared the observed maps to a simple toy model describing a vertical corrugation propagating in the Galactic stellar disc. The toy model, based on the 0^th moment of the Collisionless Boltzmann Equation (see details in Appendix D), is not designed to provide a sophisticated and realistic description of the Milky Way, but is simply a first step towards understanding the results obtained in this work. Figure 17 (left panel) shows the prediction for a corrugation propagating towards the outer disc, which would manifest itself with a vertical velocity pattern (red curve) shifted in radius R by about half the radial width of the corrugation towards the outer disc with respect to the vertical displacement in Z (blue curve). The opposite behaviour would be expected for a corrugation moving towards the inner Galaxy (i.e. with V_Z features shifted towards the inner disc with respect to the Z pattern, see the right panel of Figure 17). It is important to note that the toy model is simply describing the expected kinematic signature for a given type of corrugation propagating in the Galactic disc, and does not explain in any way which phenomenon is at the origin of the observed perturbation. Thus, based on this simple model, a possible interpretation for our observed maps is that a vertical wave extending over a large portion of the Galactic disc is present, and it is propagating outwards. The spatial manifestation of such a wave would be a radial corrugation on top of the large-scale Galactic warp, which we detect in the vertical residuals Δ_Z. Future works based on more realistic models will help improve our understanding of the propagation of waves in the Milky Way disc.

The pattern found in our datasets is in good agreement with the numerical simulations analysed in Gómez et al. (2013), which exhibit vertical perturbations of the Milky Way disc excited by a Sagittarius-like satellite galaxy. In their Figure 6 (Panel E), they find that where the mean 〈Z〉 of the Galactic disc takes an extrema, 〈V_Z〉 ~ 0 (similar to the toy model in Figure 17). They interpret this oscillatory behaviour as the signature of a vertical wave.

Using N-body simulations of a Milky Way-like galaxy perturbed by a dwarf galaxy, Price-Whelan et al. (2015) showed that, in the right circumstances, concentric rings propagating outwards from the GC can plausibly produce vertical structures similar to the Triangulum-Andromeda stellar clouds observed in the real Milky Way. On the other hand, using cosmological high-resolution hydrodynamical simulations, Gómez et al. (2016) found a well defined and strong vertical pattern in the disc of a Milky Way-like galaxy, which was the result of a satellite-host halo-disc interaction and reproduced qualitatively the observable properties of the Monoceros ring.

It is interesting to compare the 3D shape of the corrugation detected in this work to previously known corrugations in the Galactic disc. For instance, Morganson et al. (2016) mapped the Monoceros Ring in three dimensions using PAN-STARRS1. Specifically, they map an overdensity of stars below the Galactic plane at about R ~ 14 kpc (see their Fig. 16, lower left panel), and another one above the Galactic plane at about R ~ 17 kpc (lower right panel in their Fig. 16). The features they detect are not well fitted by a galactocentric ring, and do not align with the spiral arms. Moreover, the observed features are more distant from the GC in the north than in the south (i.e. leading spiral geometry). These results are in perfect agreement with the properties of the corrugation detected in this work. Given their respective locations in the Galactic disc, and given they share a similar orientation, it is possible that the corrugation presented in this work is an inner ripple with respect to Monoceros ring. We note, however, that previous works (to our knowledge) mapping vertical corrugations in the Galactic disc do not take the warp into account. While this might not be crucial in the regions close to the warp’s line-of-nodes (i.e. approximately towards the Galactic anti-centre), where the warp amplitude is small, other directions can potentially be biased from the assumption of a perfectly axisymmetric disc, which is clearly not the case of our Milky Way.

In addition to the vertical structure and motions in the Galactic disc, we also explored the radial component of the individual stellar velocities. The obtained maps indicate that stars in the corrugation are systematically moving outwards (V_R > 0, see also Figures 16, 12 and 13). It is interesting to note the agreement between the direction of the mean radial motions V_R in the corrugation (pointing toward the outer disc) and the direction of propagation of the wave that would be suggested from the observed oscillatory behaviour Δ_Z-ΔV_Z if the above wave-like interpretation is assumed.

As mentioned in Section 1, previous studies have suggested that there is evidence of vertical waves present in the disc. In particular, Schönrich & Dehnen (2018) saw a kinematic signature of the warp as well as an additional wave-like pattern in vertical velocity vs. angular momentum space, with an amplitude of ~1 km/s and a scale-length of ~2.5 kpc in guiding the centre radius, R_g. Friske & Schönrich (2019) analysed the Gaia DR2 RV sample, covering approximately a volume of ~2-2.5 kpc from the Sun, and found a wave-like pattern in the mean galactocentric radial velocity vs. R_g space, with a short-wave pattern of order of 1.2 kpc in R_g. They found that the wave-peaks in the vertical velocity line up with the outer troughs in the radial velocity wave with a strong correlation. This is in good agreement with the results presented in this work, which are based on both the 3D distribution of stars and their kinematics on a large scale. While the typical age of the tracers is different (our datasets are mostly young), and therefore some differences can be expected, it is possible that our results are connected, or even that they are part of the same oscillatory pattern. In this work we only clearly see about a half oscillation, and the radial extent of the detected pattern suggests a wavelength of at least 4 kpc.

It is generally assumed that vertical waves in the Galactic disc do not have associated radial motions, as radial and vertical motions are usually assumed to be decoupled. The main reasons for this are as follows. Firstly, the gravitational potential in a galaxy can be approximated as separable into components that affect the radial and vertical motions independently. Vertical motion is influenced by the potential perpendicular to the disc, while radial motion is influenced by the potential within the plane of the disc. Secondly, vertical oscillations typically have a higher frequency compared to radial oscillations. In local regions of the galactic disc, stars are assumed to follow approximately epicyclic orbits where vertical oscillations with respect to the disc plane are largely independent and out of phase of their radial motions.

However, our observational results suggest a different scenario. Our maps suggest that the vertical wave detected here in the Galactic disc also exhibits an associated net radial motion. If this interpretation is correct, we can draw an intriguing analogy with ocean waves, which propagate through water as a result of a circular motion of volume elements that includes both radial and vertical components, that is, with a component of motion transverse to the direction of wave propagation, and a ‘longitudinal’ component in the direction of wave propagation. These types of waves are referred to as Rayleigh waves. While it is true that in the context of the Galactic disc, we are clearly not seeing a surface wave, the coupling of radial and vertical motions here nevertheless suggests that the stars experience both upward and outward movements as the wave propagates.

Our results also indicate that, at least for young stars, the V_R velocity field in the outer disc is dominated by an alternating pattern, whose geometry is correlated with a corrugation in the spatial vertical direction. This is not the first time that an alternate pattern in V_R has been detected. For instance, Gaia Collaboration (2018b, 2023a), Friske & Schönrich (2019); Khanna et al. (2023) and others mapped streaming motions in the Galactic disc, usually with stars typically older than the samples studied here. Gaia Collaboration (2023a) mapped the velocity field of the OB stars, but only reached 1-2 kpc in heliocentric distance due to limited availability of the line-of-sight velocities for these stars. Eilers et al. (2020) mapped a V_R radial wave using upper Red Giant Branch (RGB) stars, which they interpreted as the signature of the Galactic spiral arms. The radial wave mapped by Eilers et al. (2020) might be the RGB counterpart of the V_R feature mapped here with young stellar populations. In this work, the V_R feature was discussed in relation to a vertical corrugation, but also compared to spiral arm geometries available in literature, to take into considerations different interpretative scenarios. Similarly, Palicio et al. (2023) detected a pattern in the radial action (tracing the absolute value of V_R) of disc stars. Recently, Stelea et al. (2024) used a suite of N-body simulations to study the response of the Galactic disc to the Large Magellanic Clouds and the Sagittarius dwarf galaxy (Sgr), and found that the corrugations induced by Sgr can also reproduce the large radial velocity wave seen in the data by Eilers et al. (2020). Simulations (Monari et al. 2016; Vislosky et al. 2024, etc.) have also shown that coupling between gravitational perturbations due to the Galactic bar and spiral density modes (e.g., Lin & Shu 1964; Lin et al. 1969), can generate non-linear effects on the velocity field, and could offer other mechanisms of generating large scale waves in the disc.

Finally, it is important to note that our sample is very young, and thus may be revealing the bulk motions of the gas from which they were born. We therefore cannot expect these young stars to behave as a kinematically relaxed population. Instead they may be illuminating large-scale non-circular motions in the gas disc that could have been excited by a perturbation.

The maps obtained here also imply that the disc of the Milky Way is very active, and exhibit, at the present day, compelling evidence of past (and/or present) perturbations. This has been predicted by numerical simulations (D’Onghia et al. 2016; Laporte et al. 2018; Bland-Hawthorn & Tepper-García 2021; Tepper-García et al. 2022) and has already been shown with real data by a great variety of physical phenomena, like the phase space spirals (Antoja et al. 2018; Binney & Schönrich 2018; Bland-Hawthorn et al. 2019; Laporte et al. 2019; McMillan et al. 2022; Antoja et al. 2023; Alinder et al. 2024), local vertical waves (Widrow et al. 2012, 2014; Bennett & Bovy 2019), substructures in velocity space and ridges in V_φ vs. R space (Kawata et al. 2018; Antoja et al. 2018; Ramos et al. 2018; Bernet et al. 2022; Antoja et al. 2022; Lucchini et al. 2023). All these phenomena are showing that the disc is a complex system, and it is not always trivial to link the observational features to their possible generating mechanisms (Cao et al. 2024).

In this work, we are adding a new element to the (already complicated) picture of the disc of our Galaxy, by presenting the detection of a vertical corrugation on top of the Galactic warp, whose velocities are consistent with a wave propagating towards the outer Milky Way disc. Several theoretical works predicted the presence of similar features in the Galactic disc, triggered by a passing satellite (D’Onghia et al. 2016; Laporte et al. 2018; Bland-Hawthorn & Tepper-García 2021; Tepper-García et al. 2022). However, the specific origin of the feature detected is this work is to date unknown.

The influence of a satellite galaxy has also been proposed by Thulasidharan et al. (2022) to explain the Radcliffe Wave (hereafter RW, Alves et al. 2020), a coherent structure of dense gas clouds in the Solar Neighbourhood, extending for 2.7 kpc in length. Alternative explanations proposed in the literature include the possibility that the RW is the result from a Kelvin-Helmholtz instability (Fleck 2020). Moreover, Swiggum et al. (2022) proposed the scenario that the RW constitutes the gas reservoir of the Local (Orion) Arm, based on the observed offset between those two features. Using the proper motions of young stellar objects anchored inside the clouds, Li & Chen (2022) identified a damped wave-like pattern in the vertical velocities, with a phase difference between the spatial and vertical velocity oscillations of about 2π/3. Recently, Konietzka et al. (2024) suggested that the RW is oscillating, with a phase difference of π/2, while also drifting away from the GC with a velocity of about 5 km/s. They found that, even though a gravitational interaction with a perturber seems a natural possible origin, the RW’s stellar velocities are not fully consistent with models of a perturbed-based scenario.

It is interesting to compare the RW to the wave detected in this work. Figure 18 shows that the two features are located in different regions of the Galactic disc: the RW is a filament located relatively close to the Sun (250 pc at the closest point), while the wave mapped here is more distant (at galactocentric radii of R ⪆ 10 kpc, heliocentric distances d_HEL ⪆ 2 kpc at the closest point) and spans over a much larger area of the outer Galactic disc. Moreover, the mean characteristic wavelength of the RW is about ~2 kpc (Alves et al. 2020; Konietzka et al. 2024), while that of the feature mapped here is of at least 4 kpc. Their different properties and locations in the Galactic disc suggest that they should be treated as two distinct features, but it is possible that the two waves are somehow related. For instance, they might be two different fluctuations of a common oscillatory pattern, which permeates the Galactic disc. In this context, it is worth noting that the RW has both radial and vertical motions (Konietzka et al. 2024), with a behaviour that resembles that of the wave studied in this work. Moreover, the RW and the wave mapped here are both based on young tracers. The wave mapped here might be an extension of the RW in the outer disc. However, it is beyond the scope of this paper to assess whether this scenario is realistic. Future works will shed further light on whether the two waves are connected or not.

Fig. 16 Edge-on view of the detected corrugation using the young giant sample. Left columns : vertical spatial residuals ΔZ as a function of galactocentric radius, R, for different slices in Galactic azimuth φ, as indicated by the individual labels. Points are colour-coded by the residuals in the vertical velocity ΔVZ. Black arrows show the direction and magnitude of the median residual vertical velocity ΔVZ. Right columns : same as left, but colour-coded by the median radial velocity VR. Black arrows show the direction and magnitude of the median radial velocity VR.

Fig. 17 Predicted Z coordinate (blue) and VZ pattern (red) as a function of Galactictocentric radius R based on the simple toy model described in Appendix D (see Equation (D.2)). The curves show the case of a corrugation propagating towards the outer parts of the disc at vc=10 km/s (left panel) or the inner parts of the disc at vc = −10 km/s (right panel).

Fig. 18 Comparison between the feature studied in this work and the Radcliffe wave (Alves et al. 2020), based on the model of Konietzka et al. (2024).

7 Summary and conclusions

In this paper, we analyse the vertical structure and kinematics of the Galactic disc using two samples of young giant stars and Classical Cepheids. The two samples are complementary: Cepheids can be seen to typically large distances (~ 15 kpc), but are relatively rare (~3000 objects); on the other hand, the young giant sample is mostly confined to within 6-7 kpc in heliocentric distance, but contains significantly more stars (~17 000 stars). Notwithstanding their different properties, we found that the two samples offer a self-consistent view of the Galactic disc, as expected for two samples with similar young ages. Our main results are summarised in the following.

The warp’s three-dimensional shape. The results of our inference based on a classic m = 1 warp shape are in good agreement with previous results based on young stellar populations. We obtained a warp amplitude of ~0.45 kpc at R ~ 12 kpc and ~0.7 kpc at R ~ 14 kpc, in perfect agreement with previous estimates from Skowron et al. (2019a); Chen et al. (2019); Dehnen et al. (2023); Cabrera-Gadea et al. (2024). Both samples show that the line-of-nodes begins increasing in the direction of Galactic rotation starting at about R⪆12 kpc, in agreement with previous works (Chen et al. 2019; Dehnen et al. 2023; Cabrera-Gadea et al. 2024).

Vertical corrugation on top of the warp. The analysis of the fit residuals, Δ_Z, reveals a coherent vertical feature, which is systematically shifted upwards by about 150-200 pc with respect to the best-fit m = 1 warp shape, located approximately at galac-tocentric radiii of R - 10-14 kpc. Using the sample of young giants, which has the better statistical signal, we find that this feature is not exactly radial; rather, the feature is slightly shifted at different radii depending on the considered azimuthal position (at least, in the region covered by the young giants sample).

Corrugation interpreted as a wave propagating outwards. To study the vertical kinematics, we first infer the warp precession rate using our two samples. We obtain a prograde warp precession (in agreement with Poggio et al. 2020; Cheng et al. 2020; Jónsson & McMillan 2024), which is gradually declining with R, confirming previous works by Dehnen et al. (2023) and Cabrera-Gadea et al. (2024). After subtracting the warp model in the vertical velocities, V_Z, we investigated the residuals and found that the corrugation detected in the vertical coordinate Δ_Z (see point above) has a kinematic self-consistent counterpart. In the same region, but slightly shifted toward the outer disc, we see a similar feature with positive vertical velocities, V_Z (after the warp kinematic signature has been subtracted). The shift approximately corresponds to about half of the radial width of the corrugation observed in the Z residuals, so that the lines of the maxima of the vertical velocities roughly coincides with the lines where the Z residuals are zero, indicating an oscillatory behaviour. The comparison of the observed shift with a simple toy model suggests that this is consistent with a wave propagating towards the outer disc.

The stars located in the corrugation exhibit systematic radial motions. The feature in the outer disc with Δ_Z > 0 observed in vertical fit residuals interestingly appears in the same position of systematically positive radial motions. Stars that are located in the positive crest of the corrugation today exhibit systematic outwards motions of about 10-15 km/s, while those located outside the crest are moving towards the inner parts of the disc.

Taken together, these findings lead us to explore the hypothesis that there is a vertical wave extending over a large portion of the outer disc (R > 10 kpc) that is moving away from the GC. This movement could possibly involve both vertical and radial motions, which appear to be statistically correlated. This wave, detected in young stellar populations, could primarily be part of the gaseous component of the Galactic disc, revealed by the kinematics of the young stars which have inherited the bulk motions of the gas from which they were born. Future works based on other datasets, as well as different galaxy formation models, will reveal further details on the formation and evolution of the observed feature, as well as a more detailed dynamical description of how waves propagate in the disc of the Milky Way.

Data availability

The young giant sample is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A199

The interactive 3D view of Fig. 8 is available at https://www.aanda.org

Acknowledgements

EP thanks Dorota Skowron, Joss Bland-Hawthorn, Ralph Schoenrich and Walter Dehnen for useful comments and suggestions. EP also thanks Walter Dehnen for sharing his data. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 101063193. RD and EP are supported in part by the Italian Space Agency (ASI) through contract 2018-24-HH.0 and its addendum 2018-24-HH.1-2022 to the National Institute for Astrophysics (INAF). SK & RD acknowledge support from the European Union’s Horizon 2020 research and innovation program under the GaiaUn-limited project (grant agreement No. 101004110). SK acknowledges use of the INAF PLEIADI@IRA computing resources. 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.

Appendix A Additional information on the young giant sample selection and distance calculation

Figure A.1 shows a comparison between the distances obtained in this work for the young giant sample (see Section 2.2) and other values taken from the literature. It is well known that heliocentric distances can be studied with great accuracy for open clusters in the Galactic disc (e.g. Cantat-Gaudin et al. 2020; Spina et al. 2022; Cavallo et al. 2024). For this reason, we crossmatched our sample with the open cluster members from Hunt & Reffert (2023), finding 307 stars in common. In the left panel of Figure A.1, we show the comparison between our distances and those given by Hunt & Reffert (2023). The middle and right panels of Figure A.1 show, respectively, the comparison between the photogeometric and geometric distances from Bailer-Jones et al. (2021). 95 % of our distance estimates agree within 25 % with the values taken from the literature (for stars within 8 kpc from the Sun).

Figure A.2 shows the distribution of our sample in the R-Z plane for different slices in Galactic azimuth φ before applying the cleaning described in Section 2.3. Stars are colour-coded by metallicity [M/H] from the XGBoost catalogue Andrae et al. (2023). As we can see, the selection criteria described in Section 2.1 are not enough to select young giants in the Galactic disc: we can see a significant fraction of metal-poor contaminants at large distances from the Galactic plane (which, in contrast, is typically populated by stars that are more metal-rich).

To further refine our sample from potential contaminants, we look for potential clustered stars in velocity space. We apply the python scikit-learn implementation of the HDBSCAN algorithm (Campello et al. 2013) to our dataset. Based on the output, we find 6 tightly clustered stars in the V_R, Vφ and V_Z space, which appear to be located at the same (l, b) coordinates, at the same distance. We compare the (l, b) position of the 6 stars to the globular clusters catalogue from Harris (1996), and find that they perfectly coincide with the sky coordinates of the globular cluster NGC 6656. Those 6 stars are therefore potential contaminants, and for this reason we removed them from the final catalogue.

Appendix B Sub-selection function

We estimated the sub-selection function using the Gaia query, run on the archive, in 1 magnitude wide bins in G, and at healpix level 5 for bins in l, b. The XGBoost catalogue by Andrae et al. (2023) is based on sources in Gaia DR3 for which the XP spectra were published. Furthermore, they restrict to those sources with CatWISE photometry. At the time of writing, the Gaia archive only had an official crossmatch provided with AllWISE which is shallower than CatWISE. However, we checked that our catalogue of young giants, owing to their relative brightness, have an entry both for AllWISE & CatWISE. Thus, we can still proceed with using the ratio-ing method directly on the archive but using the counts based on AllWISE instead, as shown below:

SELECT healpix_, phot_g_mean_mag_, COUNT(*) AS n, SUM(selection) AS k

FROM (SELECT

to_integer(GAIA_HEALPIX_INDEX(5,source_id)) AS healpix_, to_integer(floor((phot_g_mean_mag - 4)/1.)) AS phot_g_mean_mag_, to_integer(

IF_THEN_ELSE(has_xp_continuous=’true’,1.0,0.0) )

AS selection FROM gaiadr3.gaia_source as gdr3 left outer join gaiadr3.allwise_best_neighbour as wisenb USING (source_id)

WHERE phot_g_mean_mag > 4 AND phot_g_mean_mag < 15 and wisenb.source_id is not null )

AS subquery GROUP BY healpix_, phot_g_mean_mag_

Listing 1. ADQL query for the Gaia DR3 data considered in the sub-selection function.

A sky map of the resulting selection function S(l, b, G) for the XGBoost catalogue in the magnitude bin 14 < G < 15 is shown in Fig B.1. Maps of S(l, b, G) at other apparent magnitudes (not shown) are very similar.

Appendix C Corner plots

Figure C.1 and C.2 show, respectively, the corner plots of the global parametric fit for the young giant sample and the Cepheid sample.

Appendix D Toy model

Here, we present a simple toy model for a radial corrugation propagating in the Galactic disc, and derive the expected mean vertical velocities $\overline{V_Z}$ based on the 0^th moment of the collisionless Boltzmann Equation (CBE)⁴.

We stress that the goal of this section is not to give a sophisticated and realistic description of the Galaxy, but simply to make a first step towards the understanding of how waves propagate in the disc, and give a possible framework to interpret the observational results presented in this work. We consider a simple double-exponential disc, whose density distribution is modelled by the equation $$\rho (R,\varphi ,z,t)={\rho }_0{e}^{-\frac{(R-R_{\odot })}{L_R}}{e}^{-\frac{|z-Z_C(R,\varphi ,t)|}{h_z}},$$(D.1)

where L_R and h_z are the scale-length and scale-height of the disc, respectively, and the coordinate Z_C(R, φ, t) is the time-dependent vertical displacement of the stellar disc associated with the corrugation. For the sake of brevity, here we analyse the case of a m = 0 corrugation (i.e. radial ripples), which has no azimuthal dependency. For simplicity, we assume a corrugation described by a sinusoidal behaviour $$Z_C(R,t)=A_C\sin (\frac{2\pi }{{\lambda }_C}(R-R_{C0}-v_{C}t)),$$(D.2)

where A_C is the amplitude of the vertical corrugation, λ_C is the wavelength of the radial oscillation, R_C0 is the radial offset with respect to the GC, t is time and v_C is the propagation velocity of the corrugation. For positive (negative) values of v_C, the corrugation moves towards the outer (inner) Galactic disc. Moreover, we assume for the moment that the mean radial motions $\overline{V_R}=0$ (noting that we relax this assumption later) and that the mean vertical velocity $\overline{V_Z}$ does not depend on Z (but only on (R, φ)). These assumptions help us to keep the problem simple, and build upon that by adding new terms, if necessary. Taking the above considerations into account, the CBE in terms of galactocentric cylindrical coordinates becomes $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\overline{V_Z}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }Z}=0,$$(D.3)

and we can solve it for the vertical velocity $$\overline{V_Z}=-2\pi \frac{A_C}{{\lambda }_C}{v}_c\cos (\frac{2\pi }{{\lambda }_C}(R-R_{CO})),$$(D.4)

which is evaluated at t = 0 (today). The parameters A_C and λ_C can, in general, be function of R (see for example Xu et al. 2015), but here for simplicity we show the case where they are constant. Figure 17 shows the displacement in Z (blue curve) and in $\overline{V_Z}$ (red curve) caused by a vertical corrugation Z_C with A_C = 0.2 kpc, λ_c = 4 kpc, R_C0 = 10 kpc, and propagating toward the outer parts of the disc with a velocity v_c=10 km/s (left panel), or towards the inner parts of the disc at v_c=−10 km/s (right panel).

Fig. A.1 Comparison between different distance estimates. Left panel: Subsample of young giants classified as members of open clusters by Hunt & Reffert (2023). The comparison between the distances dYG calculated in Section 2.2 and the distances to each cluster from Hunt & Reffert (2023) is shown. The red dashed line shows the line of equality. Middle panel: Same as left panel, but now with the distances dYG compared to the photogeo distances from Bailer-Jones et al. (2021). Right panel: Same as middle panel, but now compared to the geometric distances from Bailer-Jones et al. (2021).

Fig. A.2 Spatial distribution of the young giant candidates before applying the cuts from Section 2.3. Stars are colour-coded by metallicity to better highlight the (metal-poor) contaminants at high Z at different slices in Galactic azimuth φ. The location of the Sun is indicated by the black dot.

As we can see, in this simple experiment, the vertical velocity pattern is slightly shifted outwards (inwards) with respect to the vertical displacement for a corrugation moving towards the outer (inner) parts of the disc. Similar considerations remain valid when less simplifying assumptions are made to Equation D.3, for instance including a non-zero $\overline{V_Z}$ pattern similar to the one observed in the data, or an azimuthal dependency in the shape of the corrugation (e.g. assuming a pitch angle of ±10°).

Fig. B.1 Mollweide projection of the selection function S(l, b, 14 < G < 15), for the XGBoost sample in Galactic coordinates. Here we show the magnitude bin 14 < G < 15, at HEALpix level 5. The yellow regions indicate the highest completeness, while blue regions are incomplete.

Fig. C.1 Corner plot for the global parametric fit with a twisted line-of-nodes using the young giant catalogue.

Fig. C.2 Corner plot for the global parametric fit with a twisted line-of-nodes using the Cepheid catalogue.

References

All Tables

All Figures

	Fig. 1 Comparison between the distribution of XGBoost stars (Andrae et al. 2023) in the Kiel diagram and the PARSEC isochrones for different metallicities. The portion of the Kiel diagram selected in this work is shown by a dark blue dashed line.
In the text

Fig. 2 Left panel: comparison between the photometric parallaxes (obtained from Cepheids distance estimates from Skowron et al. 2024) and trigonometric parallaxes without zero-point correction (black), with Lindegren et al. (2021) correction (red), and Lindegren et al. (2021) + additional offset correction (blue). See main text for a detailed explanation. Right panel: derived distances for the 593 young giants (dYG), which were also classified as Cepheids, compared to their photometric (Leavitt Law) based distances (dS 24).

In the text

	Fig. 3 Apparent G magnitude distribution of the young giant sample.
In the text

	Fig. 4 Vertical section of the Galactic disc as a function of galacto-centric radius, R. Over-dense regions are the spiral arms of the Galaxy seen edge-on. Stars are colour-coded by the probability density function at each bin, calculated as the number of stars in each bin over the total number of stars per bin area, where the bin area is 60 pc (in R) times 30 pc (in Z).
In the text

	Fig. 5 Distribution of the young giant sample in the Galactic plane. The Sun’s position is shown by the symbol. The GC is to the right, at (X,Y)=(R, 0 kpc), with R = 8.277 kpc (GRAVITY Collaboration 2022). Galactic rotation is clockwise.
In the text

	Fig. 6 Same as Figure 4, but now in overlapping bins in Galactic azimuth φ. Azimuthal ranges are chosen to cover the regions where our sample has more data. Isodensity contours are show as black solid lines.
In the text

Fig. 7 Inferred amplitude (left panel) and phase angle (right panel) for the m = 1 classical warp model using the young giants sample (orange/brown markers) and the Cepheids sample (blue/light blue markers) based on different approaches. Individual points show the results of the m=1 tilted ring model, while solid lines correspond to the parametric fit, performed on a star-by-star basis, of the twisted line-of-nodes for the young giant (in brown) and Cepheid (in blue) samples. Orange (light blue) transparent lines represent 5 000 random realisations of the MCMC chain for the fits of the young giant (Cepheid) sets.

In the text

Fig. 8 Comparison between the maps of median Z-coordinate from observations, best-fit model and residuals. Top-left: median Z-coordinate of the young giant sample. The position of the Sun is shown by the black cross at (X,Y)=(0 kpc, 0 kpc). Middle: prediction from the m=1 tilted rings model (left), straight line-of-nodes parametric model (centre) and twisted line-of-nodes parametric model (right). Bottom: corresponding vertical residuals for the m=1 tilted rings model (left), straight line-of-nodes parametric model (centre), and twisted line-of-nodes parametric model (right). Dotted lines in the three lower panels show rings at constant radii of R=10, 12, and 14 kpc from the GC. An interactive 3D view of this figure is available online.

In the text

	Fig. 9 Same as Figure 8, but using the Cepheids sample. The black cross shows the position of the Sun, while the black dot shows the position of the GC. Due to the different spatial coverage, the range in Z is different compared to Figure 8.
In the text

Fig. 10 Vertical residuals, ΔZ, of the m=1 straight line-of-nodes model using young giants (left panel) and Cepheids (right panel), compared to spiral arm geometries taken from the literature: Reid et al. (2019) (black solid lines, both panels), Poggio et al. (2021) (grey shaded contours, left panel), and Drimmel et al. (2024) (coloured lines, right panel). For the Drimmel et al. (2024) model, different colours indicate the Perseus arm (grey), Local/Orion arm (orange), Sagittarius-Carina arm (blue), and Scutum arm (red).

In the text

Fig. 11 Same as Figure 10, but now compared to different curves, plotted for visual reference. The dashed-dotted and solid black lines show a leading spiral curve with pitch angles −15° and −7.5°, respectively (see text for more details). The black dotted line shows a circle (pitch angle 0°) at constant radius R=12.3 kpc. In the right panel, the Perseus arm from Drimmel et al. (2024) (grey line) is also shown for reference.

In the text

Fig. 12 Comparison between the residuals in the vertical coordinate ΔZ = Z - ZWARP (left panel), the median radial motions VR (middle panel) and the residuals in the vertical velocity ΔVZ = Vz — VZ,WARP for the young giant sample. Dashed lines give a rough indication of the region where positive residuals ΔZ are obtained with this sample. Dotted lines are the same, but shifted by half the amplitude of the region (see text).

In the text

	Fig. 13 Same as Figure 12, but now for the Cepheids sample.
In the text

Fig. 14 Face-on view of the median radial motions, VR, for the young giant sample (left panel) and the Cepheids sample (right panel), compared to different spiral arm geometries from literature: Reid et al. (2019) (black solid lines, both panels), Poggio et al. (2021) (grey shaded contours, left panel) and Drimmel et al. (2024) (coloured lines, right panel). For the Drimmel et al. (2024) model, different colours indicate the Perseus arm (grey), Local/Orion arm (orange), Sagittarius-Carina arm (blue), and Scutum arm (red).

In the text

	Fig. 15 Comparison between the obtained values of warp precession with our young giants (red points) and Cepheids sample (black points), together with other values from the literature.
In the text

Fig. 16 Edge-on view of the detected corrugation using the young giant sample. Left columns : vertical spatial residuals ΔZ as a function of galactocentric radius, R, for different slices in Galactic azimuth φ, as indicated by the individual labels. Points are colour-coded by the residuals in the vertical velocity ΔVZ. Black arrows show the direction and magnitude of the median residual vertical velocity ΔVZ. Right columns : same as left, but colour-coded by the median radial velocity VR. Black arrows show the direction and magnitude of the median radial velocity VR.

In the text

	Fig. 17 Predicted Z coordinate (blue) and VZ pattern (red) as a function of Galactictocentric radius R based on the simple toy model described in Appendix D (see Equation (D.2)). The curves show the case of a corrugation propagating towards the outer parts of the disc at vc=10 km/s (left panel) or the inner parts of the disc at vc = −10 km/s (right panel).
In the text

	Fig. 18 Comparison between the feature studied in this work and the Radcliffe wave (Alves et al. 2020), based on the model of Konietzka et al. (2024).
In the text

Fig. A.1 Comparison between different distance estimates. Left panel: Subsample of young giants classified as members of open clusters by Hunt & Reffert (2023). The comparison between the distances dYG calculated in Section 2.2 and the distances to each cluster from Hunt & Reffert (2023) is shown. The red dashed line shows the line of equality. Middle panel: Same as left panel, but now with the distances dYG compared to the photogeo distances from Bailer-Jones et al. (2021). Right panel: Same as middle panel, but now compared to the geometric distances from Bailer-Jones et al. (2021).

In the text

	Fig. A.2 Spatial distribution of the young giant candidates before applying the cuts from Section 2.3. Stars are colour-coded by metallicity to better highlight the (metal-poor) contaminants at high Z at different slices in Galactic azimuth φ. The location of the Sun is indicated by the black dot.
In the text

	Fig. B.1 Mollweide projection of the selection function S(l, b, 14 < G < 15), for the XGBoost sample in Galactic coordinates. Here we show the magnitude bin 14 < G < 15, at HEALpix level 5. The yellow regions indicate the highest completeness, while blue regions are incomplete.
In the text

	Fig. C.1 Corner plot for the global parametric fit with a twisted line-of-nodes using the young giant catalogue.
In the text

	Fig. C.2 Corner plot for the global parametric fit with a twisted line-of-nodes using the Cepheid catalogue.
In the text