© The Authors 2023

1. Introduction

The European Space Agency Gaia mission has transformed our understanding of the Milky Way, thanks to its ability to trace the motion of stars in the sky (Gaia Collaboration 2023a). The observation of these movements has allowed us to see the Galaxy as an evolving system. Components that were previously thought to be distinct (the thin disc in the Galactic plane with ongoing star formation, the more diffuse and older thick disc, the central bulge, and the extended stellar halo) now appear to be interlinked formation phases of a system in clear interaction with its environment. In particular, studies of stellar orbits and kinematics have uncovered a considerable proportion of merger debris in the halo (e.g. Helmi et al. 2018; Belokurov et al. 2018; Malhan et al. 2018; Myeong et al. 2019; Helmi 2020, and references therein) and the Galactic disc (e.g. Sestito et al. 2020; Re Fiorentin et al. 2021). Additionally, investigations of stellar motions have revealed the Galactic disc phase mixing process, which is subsequent to a mildly disturbed state (Antoja et al. 2018). A massive disc-crossing perturber (e.g. Binney & Schönrich 2018) – possibly akin to the Sagittarius dwarf galaxy (Laporte et al. 2019a; Bland-Hawthorn et al. 2019) – or a strong buckling of the stellar bar (e.g. Khoperskov et al. 2020) are the most likely interpretations of this phenomenon. A recent or ongoing encounter with a satellite galaxy seems also to be responsible for the rapidly precessing disc warp (Poggio et al. 2020, 2021a). In summary, the picture of a ‘living and breathing Galaxy’ has emerged thanks to Gaia data (Belokurov 2019; Brown 2021).

Despite the above-mentioned transformational results, stellar motions alone do not allow a complete reconstruction of the intricate puzzle of Galactic history. The orbit of a star evolves in response to fluctuations in the Galaxy’s gravitational field (e.g. Sellwood & Binney 2002). As a consequence, the reconstruction of the history of the Milky Way based on the interpretation of stellar motions in terms of evolutionary processes is hampered by degenerate explanations and the complex interplay of different physical mechanisms.

Indeed, understanding how galaxies like the Milky Way form and evolve remains a challenge. In the cold dark matter Universe, galaxies grow through a sequence of merger and accretion events. However, the impact of these events on the evolution of a galaxy is extremely difficult to predict because of the complex physics of baryons. As a consequence, studying the chemo-physical properties of matter is essential to comprehend the Galaxy in which we live. Fortunately, we have a powerful tool at our disposal: stellar spectroscopy.

The study of stellar spectra gave rise to the science of astrophysics in the 19th century (e.g. Huggins & Miller 1864; Huggins & Huggins 1899) and, since then, the varying characteristics of spectral absorption lines have been used by researchers to decipher the physical properties of stars (Maury & Pickering 1897; Cannon & Pickering 1918). Stellar parametrisation became a powerful decryption tool, allowing us to unveil the chemical composition of stellar atmospheres (Payne 1925a,b), and to provide observational evidence of the stellar nucleosynthesis theory (Burbidge et al. 1957). Stars form during the collapse of molecular clouds of gas and dust. Like alchemists, stars of different masses synthesise all chemical elements except hydrogen¹; they partially return them in the later stages of their life into the interstellar medium, from which new stars are born. As a consequence, the stellar chemical composition evolves from one generation to the next, and reflects the gas conditions at the time and place of the formation of a star. Moreover, contrary to stellar motions, the chemical abundances of a star’s atmosphere are conserved² from its birth, and can therefore be used to trace its origin. Chemical abundances break degenerated dynamical scenarios with a variety of conserved parameters (e.g. they play a key role in merger debris studies; e.g. Helmi 2020). Therefore, stellar atmospheres record the past in their chemical abundances, allowing a look-back time that varies between a few hundred million years and the age of the first stars in the Universe.

In this framework, previous intermediate Gaia data releases had to be complemented with chemical data from ground-based observations. However, ground-based spectroscopic surveys like GALAH (e.g. Buder et al. 2021), APOGEE (Abdurro’uf et al. 2022), Gaia-ESO Survey (Gilmore et al. 2022; Randich et al. 2022), and RAVE (Steinmetz et al. 2020), despite the recent improvement of multiplex capabilities, are still hampered by spatially biased samples. In addition, the inhomogeneity induced by different analysis procedures, targeted stellar types, and spectral configurations blur the collected chemical information. Moreover, ground-based spectroscopy suffers from time-dependent effects such as the Earth’s atmospheric absorption and instrumental systematic effects, which are difficult to model with discontinuous data collections.

Fortunately, the context is now evolving favourably. Gaia Data Release 3 (DR3; Gaia Collaboration 2023c) opens a new era of all-sky spectroscopy and chemo-physical characterisation of Galactic stellar populations, and includes a new transformational data set that confirms Gaia’s leading role in the golden age of Galactic archaeology: the largest homogeneous spectral analysis performed so far with a total of 5 594 205 stars observed by the Radial Velocity Spectrometer (RVS; Cropper et al. 2018; Katz et al. 2023) and parametrised by the General Stellar Parametriser - spectroscopy (GSP-Spec; Recio-Blanco et al. 2023). With continuous data collection for 34 months outside the Earth’s atmosphere, and a large volume coverage reaching distances of about 8 kpc from the Sun (thanks to the population of giant stars), the Gaia DR3 spectroscopic survey provides stellar parameters and chemical abundances in all major Galactic populations, sharpening our global view of the Milky Way. In addition to the sky coverage advantage, it is worth comparing this Gaia DR3 GSP-Spec catalogue with high-resolution ground-based surveys in other crucial characteristics for Milky Way studies, as the number of analysed stars, the limiting magnitude, and the explored chemical diagnostics. For magnitudes brighter than G = 13.6³, there are more stars in the Gaia DR3 GSP-Spec catalogue than in any other ground-based survey (with both GALAH and APOGEE representing only ∼8% of Gaia GSP-Spec). For magnitudes fainter than 13.6, Gaia GSP-Spec has about 61 000 stars (reaching G = 16.1 mag), GALAH has about 130 000 stars (20% of the survey, reaching G = 18 mag) and APOGEE about 314 000 stars (43% of the survey, reaching G ≃ 20 mag). Concerning the nucleosynthesis diagnostics⁴ through individual abundance estimates, Gaia DR3 GSP-Spec explores five different nucleosynthetic channels with 13 chemical elements, while GALAH covers seven nucleosynthetic channels with 21 elements, and APOGEE six channels with 24 elements.

The goal of this paper is to demonstrate the scientific quality of Gaia’s chemical cartography through a chemo-dynamical analysis of disc and halo populations. To this purpose, Sect. 2 presents the data that are used, including (i) DR3 atmospheric parameters and chemical abundances (Sect. 2.1), (ii) DR3 radial velocities (Sect. 2.2), (iii) EDR3 astrometric data and distances (Sect. 2.3), (iv) a set of stellar velocities and orbits specifically derived for this work (Sect. 2.4), and (v) the definition of working subsamples (Sect. 2.5) allowing us to optimise the scientific analysis, and illustrating the use of quality flags defined in Recio-Blanco et al. (2023).

In Sect. 3 we present the global chemical properties of the Milky Way through sky and Galactic maps (Sect. 3.1) and explore selection function effects (Sect. 3.2). Section 4 presents the radial and vertical chemical gradients of field stellar populations. In Sect. 5 we present our analysis of large-scale chemo-kinematical correlations, while in Sect. 6 we explore the relation between the orbital parameters and stellar chemistry. Subsequently, Sect. 7 is dedicated to chemo-dynamical relations in solar cylinder populations using individual element abundances, and in Sect. 8 we use the open clusters population to study chemo-kinematical correlations and the temporal evolution of disc radial chemical gradients.

Finally, Sect. 9 summarises the results of our Galactic chemo-dynamical analysis using the Gaia RVS GSP-Spec database. In particular, we discuss the observed chemical markers of Milky Way structure (Sect. 9.1), disc kinematic disturbances (Sect. 9.2), and satellite accretion (Sect. 9.3). This is completed with the examination of the detected chemo-dynamical trends of the last billion years (Sect. 9.4) and, finally, with a discussion of the Sun’s chemo-dynamical properties in the context of its Galactic environment (Sect. 9.5). Our overall conclusions are presented in Sect. 10.

2. Data

2.1. Stellar atmospheric parameters and chemical data

This work makes use of the stellar physical parameters and chemical abundances derived from the Gaia RVS spectra by the GSP-Spec module and available through the astrophysical_parameters table of Gaia DR3. It is worth mentioning that the present work does not use the global metallicity [M/H] derived from BPRP spectra by the General Stellar Parametrizer from Photometry (GSP-Phot) and published in the GaiaSource and AstrophysicalParameters tables (mhgspphot field). Although GSP-Phot metallicities are suitable for different scientific purposes, their application to large-scale Galactic chemo-dynamical studies requires a calibration that at the time of writing this paper was not available.

GSP-Spec estimates the main atmospheric parameters (effective temperature T_eff, stellar surface gravity log(g)⁵, global metallicity [M/H]⁶, and the global abundance of α-elements⁷ with respect to iron [α/Fe]), together with the individual abundances of 12 different chemical elements from RVS spectra of single stars. The RVS wavelength range is [845 − 872] nm, and its medium resolving power is R = λ/Δλ ∼ 11 500 (Cropper et al. 2018). This spectral parametrisation is based on the MatisseGauguin GSP-Spec workflow and is described in detail in the GSP-Spec processing paper (Recio-Blanco et al. 2023). It is worth recalling that the GSP-Spec [M/H] estimation considers all the non-α metallicity indicators in the observed spectra thanks to a four-dimensional synthetic spectra grid including not only the [M/H] dimension but also the [α/Fe] one. Non-α indicators are dominated by Fe lines in the RVS domain. As a consequence, the estimated [M/H] value follows the [Fe/H] abundance with a tight correlation.

In the following sections, T_eff is taken from the teffteffgspspec field; log(g) comes from the logggspspec field; [M/H] is taken from mhgspspec; and [α/Fe] corresponds to alphafegspspec with a calibration⁸ imposed that requires a zero average value for [α/Fe] in the solar neighbourhood for any gravity (see Recio-Blanco et al. 2023, Table 4).

In a similar way, all the stellar individual chemical abundances come from the GSP-Spec estimates. In particular, this paper makes use of [N/Fe] (nfegspspec), [Mg/Fe] (mgfegspspec), [Si/Fe] (sifegspspec), [S/Fe] (sfegspspec), [Ca/Fe] (cafegspspec), [Ti/Fe] (tifegspspec), [Cr/Fe] (crfegspspec), [Fe/M] (femgspspec), [Ni/Fe] (nifegspspec), and [Ce/Fe] (cefegspspec). As for the [α/Fe] estimates, these individual abundances were calibrated following the prescriptions indicated in Table 4 of Recio-Blanco et al. (2023). It is important to note here that GSP-spec assumes the reference solar abundances of Grevesse et al. (2007).

In addition, we make use of the GSP-Spec quality flags reported in the flagsgspspec string chain (which consists of 41 characters) defined in Recio-Blanco et al. (2023). For example, we make use of the first three characters in this chain (that is, flagsgspspec[0], flagsgspspec[1], and flagsgspspec[2], reporting on the degree of parameter biases from line broadening) and called vbroadT, vbroadG, and vbroadM, respectively (see naming convention in Recio-Blanco et al. 2023).

Finally, the uncertainty on any derived parameter (Paramunc) or abundance ([X/Fe]unc) is defined as half of the difference between its upper and lower confidence levels (e.g. Teffunc = [teffgspspecupper − teffgspspeclower]/2).

2.2. Radial velocities

The complete GSP-Spec sample contains 5 594 205 stars in total (based on the query flagsgspspec is not null). DR3 provides radial velocities (radialvelocity, V_Rad) for 33 834 834 stars (Katz et al. 2023) through the gaia_source table. The GSP-Spec sample is a subset of the V_Rad sample because the GSP-Spec sample was selected based on the signal-to-noise ratio (S/N) of the mean RVS spectra. An unpublished GSP-Spec S/N > 20 was used (Recio-Blanco et al. 2023). This was found to be very similar to rvexpectedsigtonoise⁹ in the gaia_source table. V_Rad is used to Doppler shift RVS CCD spectra to rest before combining them into a mean RVS spectrum (Seabroke et al., in prep.). The sensitivity of GSP-Spec parametrisation to spectra that are not perfectly corrected for their radial velocity shift is flagged in the GSP-Spec flagsgspspec string. In particular, characters 3 to 5 (flagsgspspec[3], flagsgspspec[4] and flagsgspspec[5]), called vradT, vradG, and vradM respectively, report on the degree of parameter biases from V_Rad uncertainties (see Recio-Blanco et al. 2023, for more details on these flags).

2.3. Astrometric data and distances

High-precision astrometric parameters (α, δ, μ_α*, μ_δ, ϖ) from Gaia EDR3 (Gaia Collaboration 2021a) complement the above-described radial velocities, allowing three-dimensional study of stellar velocity. Thanks to the bright magnitude limit of the spectroscopic sample, the median parallax uncertainty is better than 20 μas for most of our targets and the median uncertainty increases up to a maximum of ∼40 μas for the faintest stars with magnitude G ≃ 14.

Based on these Gaia EDR3 astrometric data, Bailer-Jones et al. (2021) geometric distances, r_Geo, are adopted in this study. To test the implications on our analysis of the Galaxy prior (a 3D model of the stellar distribution and extinction) adopted by Bailer-Jones et al. (2021), we defined an astrometric control sample (ACS, Appendix A). This sample is composed of stellar distances, r_ϖ, exclusively based on the observed Gaia trigonometric parallaxes after removing the individual parallax zero points, Δ_ϖ, following Lindegren et al. (2021). In addition, the ACS selection fulfills the following criteria: (i) a parallax uncertainty of less than 20% derived from (ii) five-parameter astrometric solutions with renormalised unit weight error (RUWE) < 1.4, (iii) sources not flagged as duplicated by the Gaia cross-matching algorithm (Torra et al. 2021), and (iv) with full photometry (G, G_BP, G_RP) and radial velocity available. To compare the two distance estimates, we consider the fractional residual (r_Geo − r_ϖ)/r_ϖ for objects in common. The median value of the residuals is −0.04% and its median absolute deviation is 0.2%, showing very good overall agreement between both estimates. As shown in Appendix A, the main discrepancy occurs in the Galactic plane, where Bayesian distances are greater than the pure trigonometric parallax-based distances by ≲4%. These very small discrepancies have a negligible impact on the results of this work and, as a consequence, r_Geo distances are adopted in the following.

2.4. Galactocentric positions, velocities, and orbits

For the computation of the stellar positions (Galactocentric Cartesian X, Y, and Z positions, and cylindrical radius R) and the Galactocentric cylindrical velocities (V_R, V_ϕ, and V_Z in a right-handed frame with positive V_ϕ for most of the disc stars), we adopted the Sun’s Galactocentric position (R, Z)_⊙ = (8.249, 0.0208) kpc (Gravity Collaboration 2021; Bennett & Bovy 2019) and Galactocentric cylindrical velocities (V_R, V_ϕ, V_Z)_⊙ = ( − 9.5, 250.7, 8.56) km s⁻¹ (Reid & Brunthaler 2020; Gravity Collaboration 2021), as well as Gaia’s α, δ, line-of-sight velocity, and Bailer-Jones et al. (2021) geometric distances for each target. The orbital parameters (actions J_R, J_ϕ, J_Z, apocenter and pericenter distance R_a, R_p, maximum distance from the plane reached during the star’s orbit Z_max) are computed using the Stäckel fudge method Binney (2012) implemented in the Galpy code (Bovy 2015), combined with a rescaled version of the McMillan (2017) axisymmetric Galactic potential. In this model, the Galactocentric distance of the Sun is R_⊙ = 8.249 kpc while the circular velocity of the Local Standard of Rest is set to 238.5 km s⁻¹ (Gravity Collaboration 2021; Schönrich et al. 2010). The lower and upper confidence limits for each of these parameters are obtained by drawing Monte Carlo realisations of μ_α*, μ_δ, V_Rad, and distance, and propagating them through the positions, velocities, and orbits. We assumed Gaussian symmetric distributions for the proper motions and line-of-sight velocities, and a broken Gaussian distribution for distance (based on the lower and upper confidence limits of Bailer-Jones et al. 2021). The data on Galactocentric positions, velocities and orbits used in this paper can be retrieved via the Gaia archive table gaiadr3.chemical_cartography.

2.5. Working samples

To optimise the analysis of this work, different data selections and quality cuts are required. To this purpose, we defined several data samples with the following characteristics. The precise conditions used for these sample selections are presented in Appendix B, including the corresponding catalogue queries.

General sample (5 527 090 stars; Sects. 3 and 6): From the complete GSP-Spec sample (containing 5 594 205 stars in total), we excluded about 56 000 stars located near the cool end of the GSP-Spec parameter space and suffering from biases in different parameters (Recio-Blanco et al. 2023). The remaining objects without r_Geo distances were removed. For this sample, the median uncertainty in [M/H] is 0.075 dex and the median uncertainty in [α/Fe] is 0.05 dex.

Medium quality sample (4 140 759 stars; Sect. 3, Sect. 5, Sect. 6): This results from additional quality criteria rejecting stars with grid border effects and large parameter uncertainties. For this sample, the median uncertainty in [M/H] is 0.06 dex and the median uncertainty in [α/Fe] is 0.04 dex.

Gradient analysis sample (2 762 809 stars, Sect. 4): This optimises the quality of the [M/H] and [α/Fe] parameters and the astrometric quality to ensure reliable distance estimates, with the goal of estimating radial and vertical gradients. For this sample, the median uncertainty in [M/H] is 0.05 dex and the median uncertainty in [α/Fe] is 0.035 dex.

High-quality sample (2 218 573 stars, Sect. 5): This selects the most precise stellar parameters and chemical abundances. We note, in particular, that the choice of the fluxNoise = 0 filter flag imposes a selection of stars with very low parameter uncertainties. For this sample, the median uncertainty in [M/H] is 0.03 dex and the median uncertainty in [α/Fe] is 0.015 dex.

Individual abundances samples (Sect. 3 and 7): these have been defined to optimise the individual abundance distributions of different chemical elements. The number of selected stars depends on the chemical element. Specifically, as expected, the precision of abundance measurements is a function of different conditions, such as the number of analysed lines, their equivalent widths, and possible blends, which vary with effective temperature, surface gravity, mean metallicity, and other individual abundances, and is different from one element to another.

3. Global chemical properties of the Milky Way

The ∼5.5 million star database of spectroscopic stellar parameters and chemical abundances (see Sect. 2.1), as part of the Gaia all-sky data collection, provides significant radial and vertical coverage of our Galaxy, but also probes a considerable range in azimuth. In this section, we first analyse the global spatial distribution of the stars in the Milky Way and their chemical properties, and then we illustrate the selection function in radial and vertical bins. Here, we show the chemical properties and spatial coverage of three representative stellar populations selected from the Kiel diagram of the near-plane objects: massive young stars, bright red giant branch (RGB) objects and hot turn-off (TO) stars.

3.1. Spatial and chemical trends

Sky maps in Galactic coordinates. The distribution of all the stars in the General sample is presented in Galactic coordinate maps colour coded according to distance (Fig. 1), metallicity (Fig. 2), and [α/Fe] abundance (Fig. 3).

We first note that the ecliptic pole scanning law (EPSL; Boubert & Everall 2020) pattern is visible in Fig. 3 with slightly higher [α/Fe] values. Between 25 July and 21 August 2014 (the first weeks of Gaia’s nominal mission), Gaia operated in EPSL mode (Gaia Collaboration 2016). During this period, stars close to the ecliptic poles and along great circles connecting them were observed a greater number of times than during Gaia’s nominal scanning law. This leads to spectra with higher S\N towards the EPSL footprint. As Gaia continues to observe in nominal scanning law, the S/N will homogenise over the sky for future data releases. The higher S/N of the spectra along the EPSL pattern allows the parameterisation of fainter (and hence more distant) stars in these regions, leading to a slightly lower mean metallicity and especially larger mean [α/Fe] values per pixel.

The metal-rich and relatively α-low thin disc is easily identified in Figs. 1–3. In the direction of the Galactic bulge, outside the Galactic plane, one sees a dominant metal-poor population. Indeed, due to the fact that the bulge is more radially concentrated, many thick disc stars are present in these regions, decreasing the median metallicity. The regions located far from the plane are filled with a mixture of nearby metal-rich α-low disc stars and more distant metal-poor α-rich halo stars (see also Sect. 4).

To complement this first global view, we cut Fig. 2 into different distance intervals from the Sun (see Fig. 4). Similar figures for the [α/Fe] distribution are provided in the Appendix (Fig. C.1). The scanning law features described above are visible in the closest distance bin (D < 0.5 kpc, upper left panel). Figures 4 and C.1 show that, whereas the sky is filled by metal-rich and α-low thin disc stars within 500 pc of the Sun, the Galactic discs start to progressively emerge with increasing distance. The more metal-poor and α-rich halo dominates beyond ∼4 kpc. In the more distant panel, beyond 6 kpc, the thick disc is well seen on both sides of the Galactic centre and the thin disc has almost disappeared because of interstellar extinction. An asymmetry can also be seen in the metallicity map close to the Galactic centre (with the central region being more metal-rich at positive longitudes than at negative ones), which may result from the presence of the central bar and an inhomogeneous distance distribution around the Galactic centre in this more distant bin due to extinction. Finally, the Small and Large Magellanic Clouds (SMC and LMC, respectively) are also clearly visible in this panel, both being rather α poor. We note that the LMC appears rather metal rich. Our LMC sample is indeed dominated by bright massive and young stars (referred to as blue loop stars by Gaia Collaboration 2021b), which likely sample the metal-rich LMC populations.

Galactic maps in Cartesian coordinates. Figure 5 (top panels) shows three maps in (X–Y) Cartesian coordinates of all the General sample stars superimposed on the Milky Way sketch of Hurt & Benjamin (Churchwell et al. 2009). From left to right, the maps are colour coded according to stellar counts, mean metallicity, and enrichment in α-elements with respect to iron. We note that, in the panels colour coded according to [M/H] and [α/Fe], the apparent bulge asymmetry described in Figs. 4 and C.1 is again seen between the Sun and the Bar.

The bottom row of Fig. 5 shows the distribution of the same stars in the (R, Z) plane, where R is the Galactic centre distance and Z the distance from the Galactic plane. It can be seen that the spatial coverage of this sample with measured [M/H] and [α/Fe] is very large, spanning several kiloparsec. Close to the Galactic plane and towards the inner regions, a thin stripe of metal-rich low-α stars is clearly seen. Its vertical extension is about 300 pc below and above the plane, which is the typical scale height of the thin disc. This stripe is therefore the thin disc seen edge-on. It is not detected for R ≲ 3.5 kpc because of interstellar extinction and the resultant lack of sufficiently bright stars for analysis by GSP-Spec at smaller Galactic radii. Towards the external regions near the plane, this thin disc is also seen up to ∼14 kpc from the Galactic centre. Its flare in the outer regions is traced by the progressive increase of high-metallicity low-[α/Fe] stars farther away from the plane, as we move towards higher R values (see Sect. 5 for a chemo-kinematical perspective of the disc flare). Moreover, the stars located close to the Galactic plane but in the external parts of the disc appear more metal poor and slightly more α rich: a clear gradient in [M/H] and [α/Fe] versus R is thus revealed (see also Sect. 4). As a consequence, at large heights from the Galactic plane, internal regions are dominated by low-metallicity thick disc and halo stars, while for R > 12 kpc, flared disc stars of higher metallicity and lower [α/Fe] values dominate.

3.2. Exploration of the selection function

The large spatial coverage of the Gaia chemical database requires a global presentation of the underlying parameter distribution across the Galaxy. In addition, the large variety of stellar types included in the catalogue demands that we understand the underlying spatial distribution across the Kiel diagram. Without entering into a detailed characterisation of the selection function, in this section we explore its impact so as to help DR3 users and to support the chemo-dynamical analysis of this paper.

3.2.1. Stellar parameter distribution across the Galaxy

To dissect the covered spatial volume, we divided the (R, Z) plane into bins of 2 × 2 kpc². The resulting [R, Z] grid of Kiel diagrams colour coded according to [M/H] is shown in Fig. 6, and that of [α/Fe] versus [M/H] distributions is presented in Fig. 7. In both figures, the global (R, Z) density contour plot is shown in the background. For the purpose of clarity, only the stars in the Medium Quality sample (4 140 759 objects) are included. As anticipated from the previous section, about 90% of this sample is found between 6 and 10 kpc in R and |Z|< 2 kpc.

The majority of the stars are metal-rich and alpha-poor dwarfs on the main sequence and belonging to the thin disc. Close to the Sun, all evolutionary stellar stages from main sequence stars up to the tip of the asymptotic giant branch (AGB) are encountered with a large mixture of [M/H] and [α/Fe] combinations. Moving away from the solar location in both directions, the selection function progressively favours intrinsically brighter stars with lower log(g) values. The radial and vertical metallicity gradients can also be qualitatively appreciated (see Sect. 4 for quantification of the gradients). For |Z|> 4 kpc, only metal-poor RGB and AGB stars are observed and most of them have metallicities of between ∼ − 2.5 and ∼ − 0.5 dex with high α-enrichment. Nevertheless, a few more metal-rich stars from the flared disc are visible in the outer regions. We also note a high symmetry with respect to the Galactic plane for the Kiel and [α/Fe] versus [M/H] diagrams. As illustrated in Fig. 7, the global standard decreasing trend of [α/Fe] with metallicity can be seen across the Galaxy.

3.2.2. Kiel diagram selections near the Galactic plane

To further illustrate the selection function close to the Galactic plane, where stellar type variations are more important, we first show in Fig. 8 the Kiel diagram of the Medium Quality stars with |Z|< 1 kpc. We then define three subsamples by selecting: (i) the hottest dwarfs of the main sequence around the turn-off, (ii) the bright part of the RGB above the clump, and (iii) the massive bright stars hotter than the RGB, including blue loop stars and variables such as Cepheids. This last selection, which is biased towards young stars, is indeed dominated by stars located very close to the Galactic plane and an additional filter keeping stars at |Z|< 400 pc is applied to avoid contaminants (mainly from the RGB population that is nearby in the Kiel diagram). Hereafter, these subsamples are respectively labelled hot turn-off (HotTO), RGB, and Massive as indicated in Fig. 8.

The (X,Y) spatial distribution and chemical properties of the three populations are shown in the left and central panels of Figs. 9–11, respectively. The metallicity colour scale is the same in the three figures for comparison purposes. In addition, the right top panels in the three figures show the evolution of the calcium abundance ratio with respect to iron, [Ca/Fe], as a function of metallicity for the three populations selected in the Kiel diagram of Fig. 8. The abundance of this chemical element is the easiest to derive from RVS spectra thanks to the presence of the strong absorption Ca II IR triplet lines. Finally, the abundance trends of additional elements are shown in the right bottom panels as a function of stellar type and therefore the presence or absence of the lines of an element in the related spectra. Sulfur lines are only detected for hot dwarfs (because of its rather high excitation energy) and its abundance is illustrated in Fig. 9. The only nickel line selected in the RVS range provides reliable estimates for bright cool giants and its corresponding abundance is shown in Fig. 10. Finally, nitrogen abundances are nicely illustrated using the hot massive population in Fig. 11.

To optimise the individual chemical abundance plots, the Individual abundance sample defined in Sect. 2.5 is used. We note that, even after applying rather strict selection criteria (cf. Sect. 2.5), the number of reported abundances is huge (from ∼10⁴ to a couple of ∼10⁶, depending on the chemical species) and generally much larger than any other published abundance catalogue. As shown below, the behaviour of these chemical species is in perfect agreement with previously published catalogues of stellar abundances and with Galactic chemical evolution models (see e.g. Prantzos et al. 2018; Kobayashi et al. 2020; Matteucci 2021, and references therein for observational datasets), confirming their high quality. Finally, we note here that the individual abundance calibrations proposed by Recio-Blanco et al. (2023) have been applied.

Spatial distributions in the disc plane. First, the left panels of Figs. 9–11 reveal that almost 680 000 HotTO stars are found within about 1 kpc of the Sun, because they are not bright enough to sample more distant regions. In contrast, the RGB sample contains somewhat fewer stars (about 580 000) but allows us to explore much more extended regions up to 5–6 kpc from the Sun. Similarly, Massive stars are distributed up to about 3 kpc from the Sun. Their spatial distribution is not homogeneous and nicely reveals the closest spiral arms, that is, the Sagittarius/Carina, Local, and Persus arms, in agreement with the spatial maps derived by Poggio et al. (2021b) from Gaia EDR3 astrometry and photometry. These massive stars also show the enhanced equivalent width of the diffuse interstellar band at 862 nm, confirming their nature as tracers of the spiral arms (Gaia Collaboration 2023b).

Metallicity distributions in the Galactic plane. The metallicity distributions in the (X,Y) plane of the three samples are shown in the central panels of Figs. 9–11. Most strikingly, we see that a large fraction of the HotTO stars are metal rich (up to [M/H] ∼ +0.5 dex). Indeed, the selection of the HotTO sample favours the younger and therefore more chemically enriched TO stars because of the temperature cut at T_eff > 6000 K. Finally, the slightly more metal-poor yellowish feature centred on the solar position comes from an inhomogeneous Z distribution of these stars. In fact, near the Sun, the sample reaches slightly greater distances from the plane, with median values of around 0.2 kpc, while outside the solar vicinity, and probably driven by the younger objects in the sample, most of the HotTO stars are slightly hotter (in order to be brighter) and therefore closer to the plane at |Z|≲0.1 kpc. Therefore, the vertical metallicity gradient (see also Sect. 4) seems to explain the central metallicity feature. Nevertheless, the solar metallicity (represented by the colour of the central starred marker) is in agreement with the global [M/H] distribution of HotTO stars, as expected.

Regarding the RGB sample (cf. Fig 10), stars are on average more metal-poor than the HotTO ones, with values around [M/H] < –0.4 dex in most parts of the (X,Y) plane. However, there is a metallicity gradient that is visible for X ≳ 5 kpc, with [M/H] values decreasing with distance from the Galactic centre, as expected (see also Sect. 4). In this framework, the Sun’s metallicity is clearly higher that the average observed one in the solar surroundings. This is probably due to the underlying age distribution of RGB stars. In fact, the selected locus of the Kiel diagram is populated by a large range of stellar masses, corresponding to a wider age distribution of these objects. In addition, the Z distribution of this RGB sample includes stars farther away from the plane that reduce the median metallicity.

Interestingly, we clearly see that the distribution is not azimuthally symmetric. Stars with X in the range between 5 and 7 kpc and Y within [−2.0,+4.0] kpc seem to define a metal-rich sustructure, which probably coincides with the Sagittarius arm position. Moreover, stars with X in the range between 7 and 10 kpc show an azimuthal metallicity trend, with objects in azimuthal angles closer to the major axis of the Galactic bar (lower Y values in the figure) having higher metallicities. This could also be the consequence of the presence of the Local Arm and generally of the pitch angle of the spiral arms. For the innermost regions of the Galaxy, it could also partially result from a spatially inhomogeneous extinction distribution. In any case, these azimuthal dependencies reveal an inhomogeneous metallicity distribution, even for this upper RGB population for R ≳ 5–6 kpc. In this sense, it is worth noting that the innermost regions (X < 5 kpc) in this figure appear metal poor and therefore do not follow the radial metallicity gradient. However, as is visible in the bottom panels of Fig. 5, the inner Galaxy sample is dominated by metal-poor thick-disc stars, particularly at |Z|> 0.5 kpc.

Finally, the central panel of Fig. 11 shows that the arms internal to the position of the Sun and mapped by the Massive stars are much more metal-rich than the external ones, following again the expected radial gradient. It can also be noted that, as already discussed for the RGB sample, Massive stars show azimuthal dependencies in their metallicity distribution. In addition, it can be seen that the Sun lies between two distinct spiral arms and that its mean metallicity is higher than that of the surrounding younger Massive stars. This seems to indicate that the younger populations are chemically impoverished, which is in contradiction with a continuous chemical evolution where metallicity increases with time. It is also in contrast with the older but more metal-rich turn-off stars discussed above. This point is discussed in more detail below, when looking at other chemical diagnostics. It is worth noting that although the parametrisation of young stars from high-resolution spectra (R ≳ 100 000) can be difficult because of a combination of intrinsic factors (activity, fast rotation, magnetic fields; see for instance Zhang et al. 2021; Spina et al. 2020), this is mainly affecting particularly young stars (age ≲200 Myr) and does not alter the RVS parametrisation, which is performed at medium resolution (R ∼ 11 500). Visual inspection of the GSP-Spec parameter and individual abundance solution and synthetic spectral fitting for one typical star of the Massive population can be found in Recio-Blanco et al. (2023). No particular parametrisation problems are reported for these stars; they have generally very high S/N values. Additional support to the derived metallicities for this population comes from the validation by Clementini et al. (2023) of RR Lyrae metallicities. The authors show that spectroscopic GSP-Spec metallicities are in good agreement with those derived from the period–luminosity relation. RR Lyrae stars are indeed a subsample of our Massive population.

[Ca/Fe] distributions. We now look at the [Ca/Fe] abundances presented in the upper right panels of Figs. 9–11 (individual abundances samples). Concerning HotTO stars, Fig. 9 presents the [Ca/Fe] trend with metallicity. As expected for this faint nearby population, the distribution is dominated by a low-α thin-disc sequence, with [α/Fe] decreasing as metallicity increases and reaching [M/H] values as high as +0.5 dex. The relatively large metallicity range of this sample is therefore probably due to stars of different origins visiting the solar neighbourhood as already observed for ground-based surveys of dwarf stars. We emphasise that the decrease in [Ca/Fe] with [M/H] is continuous even for super-solar metallicity stars, in perfect agreement with Galactic evolution chemical models. This is seen for the first time with Gaia DR3 data thanks to an optimised spectral normalisation of metal-rich stars (Santos-Peral et al. 2020).

[Ca/Fe] abundances are presented for the RGB population in the right upper panel of Fig. 10. Two sequences are easily identifiable, with the thick disc population being more numerous for the RGB selection than for the HotTO one. As already mentioned, the metallicity distribution of the RGB is shifted to more metal-poor values than that of HotTO stars. We also note that the continuous decline in [Ca/Fe] abundance with metallicity seen for the HotTO sample is also visible here; however it is less clearly observable due to the presence of [Ca/Fe]-poor RGBs at metallicities larger than ∼ − 0.3 dex coming from the Massive sample, which is contiguous in our Kiel diagram selection. The [Ca/Fe] decline with [M/H] in the thin-disc sequence can also be seen by examining the upper envelope of the [Ca/Fe] distribution, which reaches [Ca/Fe] ∼ +0.0 dex at [M/H] ∼ +0.25 dex.

Finally, Fig. 11 illustrates the [Ca/Fe] abundance for about 7400 Massive stars. The corresponding distribution for the stars in the RGB population is shown in grey for reference. Surprisingly, most stars in the Massive sample are Ca poor, presenting [Ca/Fe] values down to ∼–0.3 dex with [M/H] varying between −0.5 and +0.0 dex. This confirms the apparent chemical impoverishment of this young population discussed above. We investigated whether these peculiarly low abundances could be due to a zero-point offset for the Massive stars or an artifact of the abundance calibration. The different tests that we performed, which are presented in Appendix D, allow us to conclude that, for all the combinations of calibration flavours including no calibration, the Massive sample is impoverished in [α/Fe] with respect to the thin disc sequence. In addition, the maximum metallicity reached by the Massive sample is also lower than that of the RGB and HotTO samples, independently of the chosen calibration flavour. It is also interesting to consider the fact that this type of young cool massive star in the spiral arms is rare (because of its shorter lifetime). Therefore, only a high number statistics survey, as in the Gaia/RVS one, can observe enough stars of this type to be statistically relevant in the different abundance distributions.

As noted above, the observed chemical impoverishment of the Massive population is in contradiction with a continuous chemical evolution of the thin disc. Moreover, because of the decreasing stellar density with Galactic radius, the mechanism whereby inward migration favours lower metallicities (as required to explain the chemical pattern of the Massive population) should not dominate. Chemical evolution models and/or simulations, which are beyond the scope of this paper, will be necessary to develop the interpretation of the chemical pattern seen in Massive stars.

Other individual abundance distributions. We now discuss the additional individual chemical abundance distributions presented in the right lower panels of Figs. 9–11. For HotTO stars, Fig. 9 shows the variation of sulfur abundance with metallicity for more than 36 000 stars. The behaviour of sulfur with metallicity, although exhibiting a larger spread than the [Ca/Fe], is in agreement with the observed global [α/Fe] trend, again showing a continuous decrease with metallicity. Indeed, the ratio [S/Ca] (and the [S/α] one) is almost equal to zero (with a median value of 0.02 dex) over the whole metallicity range, with a rather small median absolute deviation (MAD = 0.06 dex). This confirms the α-like nature of sulfur in agreement with previous studies (see Perdigon et al. 2021, for a recent review and consistent results derived from high-resolution spectra).

Regarding the RGB sample, about 28 000 nickel abundances are reported for the RGB sample in Fig. 10 (we adopt here [Ni/Fe]_unc < 0.15 dex to increase the size of the sample). [Ni/Fe] stays close to zero for all metallicities, in agreement with the iron-peak nature for this chemical element. Regarding the Massive population in Fig. 11, the approximately 4000 calcium-poor stars located in the spiral arms have nitrogen abundances close to zero (with a very small dispersion) for metallicities varying between −0.5 and 0.0 dex.

Finally, we show in Fig. 12 the behaviour of cerium with metallicity for the coolest part of the Massive stars plus the whole RGB sample. This allows us to select a disc stellar sample with a wide age distribution, from the younger populations in the spiral arms to the older stars in the RGB sample. The two populations show a large, and very similar spatial coverage. We remind the reader that cerium is a neutron-capture element that is predominantly produced by slow-neutron captures and belongs to the second peak of the s-elements. We find about 7000 stars with Ce estimates in these two subsamples, relaxing the Ce upper limit flag to allow values equal to 2. The upper panel of Fig. 12, which is colour coded according to stellar density, shows that the [Ce/Fe] abundance decreases with mean metallicity in agreement with previous studies (see also Sect. 7) and Galactic chemical evolution models. Moreover, there is a group of stars with very negative [Ce/Fe] for [M/H] between about −0.5 and −0.25 dex. As shown in the lower panel of Fig. 12, in which the [Ce/Fe] versus [M/H] distribution is colour coded according to [Ca/Fe] abundance, these low-[Ce/Fe] stars correspond to the Massive Ca-poor population located in the Galactic arms. Therefore, Ce abundances confirm the chemical impoverishment seen in metallicity and other chemical elements for this Massive star sample.

4. Radial and vertical gradients

In this section, we investigate the overall radial and vertical gradients for metallicity and calibrated [α/Fe] (cf. Sect. 2.1 and Recio-Blanco et al. 2023), linking these results to the ones previously found in the literature. For this purpose, we use the Gradient analysis sample defined in Sect. 2.5.

Figure 13 shows the median metallicity (upper panel) and [α/Fe] (lower panel) radial values for five different distances below (dashed lines) and above (full lines) the Galactic midplane (i.e. in a total of ten different bins in Z). This is indeed achievable with good statistics, given the large number of stars for which we have parameters; however, we note that we do not separate our sample into chemically thin or thick disc (i.e. based on their [α/Fe]-[M/H] sequence). It is interesting to note the limited range in [α/Fe] enhancement in this figure, peaking at about ∼0.2 dex. This is explained by the fact that GSP-Spec [α/Fe] abundances are governed by calcium indicators dominating the RVS α-element information in many pixels. The calcium abundance range in the thin- and thick-disc distributions is more limited than in other elements such as magnesium, as already observed in the data from ground-based surveys such as the Gaia-ESO Survey (Mikolaitis et al. 2014) and APOGEE (Abdurro’uf et al. 2022). Figure 14 is similar to Fig. 13, showing this time the vertical gradients in metallicity (upper panel) and [α/Fe] (lower panel) for four different R bins.

It is remarkable to find that both metallicity and [α/Fe] trends show similar behaviours above and below the plane, when the same |Z| range is considered. More specifically, we find that both metallicity and [α/Fe] gradients have a break in their slope at galactocentric radii R ∼ 7 kpc, which is particularly visible for |Z|≲1 kpc. This behaviour is in agreement with results from Haywood et al. (2019) using α-low stars from the APOGEE-DR14 catalogue, Kordopatis et al. (2020) from a compilation of spectroscopic catalogues of field stars (LAMOST, RAVE, GALAH-DR2, and APOGEE-DR14), and Katz et al. (2021) using both α-high and α-low stars from the APOGEE-DR16 catalogue. More specifically, we find that the median metallicity either flattens or even decreases towards the inner Galaxy, while, at the same R-position, the median [α/Fe] increases (see lower plot of Fig. 14). We note that this change of regime stays roughly similar when one selects stars based on Z_max and guiding radius (see Fig. 15).

Table 1 reports the linear [M/H] and [α/Fe] trends as a function of R, in the form of ax + b (and their associated uncertainties e_a and e_b, respectively), only considering the galactocentric radii larger than 7 kpc. We find that the metallicity slope for the closest bin to the plane is ≈ − 0.056 ± 0.007 dex kpc⁻¹. This is slightly different from but still consistent with the metallicity gradients found using Cepheids (e.g. −0.060 ± 0.002 dex kpc⁻¹Genovali et al. 2014, and −0.045 ± 0.07 dex kpc⁻¹, Lemasle et al. 2018), open clusters (e.g. −0.076 ± 0.09 dex kpc⁻¹Spina et al. 2021, and −0.068 ± 0.001 dex kpc⁻¹, Donor et al. 2020), or thin disc (i.e. low α-sequence) field stars (e.g. ∼ − 0.053 to ∼ − 0.068 dex kpc⁻¹Bergemann et al. 2014; Recio-Blanco et al. 2014; Hayden et al. 2015; Anders et al. 2017). The trends become flatter as one moves to larger |Z|, eventually becoming null or slightly positive for bins farther than 1.5 kpc from the plane (−0.005 ± 0.003 at 1.5 < Z < 2 kpc and +0.002 ± 0.004 at 1.5 < Z < 2 kpc). This flattening is again qualitatively consistent with previous results using field stars and interpreted as a smooth transition from a thin disc to more radially concentrated thick disc (see e.g. Boeche et al. 2014; Kordopatis et al. 2020; Nandakumar et al. 2022). We note that the metallicity gradient when considering the Massive stars sample (see the selection in Fig. 8; the gradient close to the plane is flatter than the one found using the RGB sample (−0.036 ± 0.002 dex kpc⁻¹). This difference might be due to the difference in age between these two samples, and is consistent with the gradient found using the youngest sample of open clusters (see Sect. 8).

As far as the ∂[α/Fe]/∂R gradients are concerned (see Table 1), for R ≳ 7 kpc, we find a value of +0.007 ± 0.001 dex kpc⁻¹ close to the plane, in agreement with the values published by Carrera & Pancino (2011), Yong et al. (2012), Mikolaitis et al. (2014), Genovali et al. (2015), Reddy et al. (2016), Casamiquela et al. (2019), Donor et al. (2020). The trend becomes steeper as one moves farther away from the Galactic plane, reaching eventually −0.016 ± 0.002 dex kpc⁻¹ for 2 < |Z|< 2.5 kpc. Despite not reaching typical thick disc [α/Fe] enhancements at |Z| where the thick disc dominates¹⁰ (> 1.5 kpc; see e.g. Hayden et al. 2015; Kordopatis et al. 2015), the qualitative trends that we find show indeed that one moves from a thin-disc-dominated population to a thick-disc-dominated population as |Z| increases. This is also seen from the bottom plot of Fig. 14, for the curves with R < 11 kpc.

The points above, especially the breaks in ∂[M/H]/∂R and ∂[α/Fe]/∂R, are more clearly illustrated by Fig. 14, which shows the vertical gradients as a function of galactocentric radius. Only four R bins are investigated, each of 2 kpc in width, for which a sufficient number of stars are available at all Z (see Fig. 5). As anticipated from the shift in the slopes of Fig. 13, we find non-zero vertical gradients at all R. Considering stars at |Z|< 1 kpc, we find that the inner disc shows steeper vertical metallicity and [α/Fe] gradients. Eventually, at the outermost radii (R = [11 − 13] kpc), we see only a mild vertical metallicity gradient and a null [α/Fe] gradient. This is in agreement with Bensby et al. (2011), Hayden et al. (2015), Kordopatis et al. (2015), Mackereth et al. (2019), Katz et al. (2021), for example, who found that the thick disc is more radially concentrated than the thin disc.

Finally, it is worth mentioning that we also find some wiggles in the radial metallicity trends for |Z|< 0.5 kpc and R ∼ 8.5 kpc – both above and below the plane – that are not related to statistical noise (these bins contain a lot of stars, as they are relatively close to the Sun). While we cannot exclude the possibility that these wiggles could partly be due to dynamical effects (e.g. resonances Antoja et al. 2018; Kawata et al. 2018), it is clear that our geometrical selection biases play an important role in producing them. Indeed, hot TO stars, that is, relatively young and super-solar metallicity stars, can only be seen relatively close to the solar neighbourhood (see Fig. 9). This introduces an age (and metallicity) spatial bias in the sense that the local median metallicity is higher than farther away, simply because TO stars cannot be seen at large distances. Similarly, red clump stars are also lost in the most distant regions, creating a secondary but non-negligible spatial bias. Selecting only giant stars above the red clump (log(g) ≤2), which can be seen in a much more unbiased way at all radii, leads to trends without metallicity wiggles (see Fig. 16) with gradients that are similar (but slightly flatter by ∼0.01 dex kpc⁻¹) than the ones quoted above (see Table 2). This indicates that the selection function plays an important role in the study of the metallicity gradients.

5. Chemo-kinematical analysis

In this section, we present the large-scale chemo-kinematical properties of the Medium Quality and High Quality samples presented in Sect. 2.5. We imposed additional conditions on the quality of the astrometric parameters (see also A.1) to ensure the use of reliable kinematical velocities. In particular, we selected stars with both five- and six-parameter astrometric solutions and a quality index RUWE < 1.4. In addition, duplicated sources have been rejected. We also performed this analysis by means of the trigonometric distances, r_ϖ, from the ACS (see Appendix A). No significant difference was detected when using the trigonometric distances with respect to the chemo-kinematical distributions presented in the following sections.

The considered spatial range in this section covers the inner disc up to R = 16 kpc, as well as vertical distances from the Galactic plane in the range −6 < Z < +6 kpc.

5.1. Toomre diagram and Galactic populations

Figure 17 shows the Toomre diagrams colour coded according to the median [M/H], for 3 468 272 stars with the Medium Quality sample within R < 16 kpc and |Z|< 6 kpc. Each panel presents the distribution versus V_ϕ within bins of ΔR = 2 kpc and ΔZ = 2 kpc. The circular dashed line defines the region where the disc stars are dominant:

(1)

Although the mean disc rotation varies as a function of R and Z, here we adopt a constant value, V₀ = 238.5 km s⁻¹, corresponding to the velocity of the LSR at the Sun position (see Sect. 2.4).

We note that the chemo-kinematical distributions appear quite symmetrical above and below the plane. Focusing on the annulus, 4 < R < 10 kpc, close to the Galactic plane, |Z|< 2 kpc, we clearly identify the stellar populations belonging to the thin disc with typical solar metallicity [M/H] ≈ −0.1 dex (red colour) and the thick disc with [M/H] ≈ −0.6 dex (green colour). The low metallicity Galactic halo with [M/H] < −1 dex is also present (blue colours). As expected, the fraction of disc stars decreases significantly above the plane for |Z|> 2 kpc. At |Z|> 4 kpc, the main stellar component is the halo formed by Gaia-Enceladus-Sausage (Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Di Matteo et al. 2019; Gallart et al. 2019) and heated thick disc stars.

In the outer disc, 10 < R < 16 kpc, the mean metallicity of the disc stars decreases to [M/H] ∼ −0.5 dex (see Sect. 4). In the last panels, 14 < R < 16 kpc, the thick disc disappears and only thin disc and halo stars are present. Moreover, we still find fast rotating V_ϕ ≈ 200 km s⁻¹ thin disc stars possibly associated to the flaring of the disc (e.g. López-Corredoira et al. 2020; Ding et al. 2021) up to |Z|∼4 kpc (cf. also Fig. E.1).

In the inner disc, 2 < R < 4 kpc, the high-metallicity disc population seems to decrease. In reality, this is a selection effect produced by the interstellar extinction that reduces the completeness of the sample close to the Galactic plane inside a given R − Z bin.

5.2. Velocity distribution V_Z versus V_ϕ

Figure 18 shows the velocity distribution V_Z versus V_ϕ of giant stars (log(g) < 3.5) with |Z|< 1 kpc selected from the Medium Quality sample (with the above described additional conditions on the astrometric parameters). This log(g) cut allows us to select a tracer population with a more homogeneous spatial coverage (cf. Sect. 3), thus removing dwarf stars to avoid selection function artifacts. However, contrary to the previously defined RGB sample, this selection includes red clump and bottom RGB stars to increase the number statistics, while losing some spatial homogeneity. The three panels represent the distributions for the subsets within the Galactocentric distance ranges 6 < R < 8 kpc (left), 8 < R < 10 kpc (middle), and 10 < R < 12 kpc (right).

No significant difference is found between the normalised density distributions above and below the Galactic plane, except for the outer annulus with 10 < R < 12 kpc. Here we find about 80 000 sources with Z > 0 kpc, compared to the approximately 59 000 sources with Z < 0 kpc. Such an effect, which is also present in the General sample, is possibly the result of greater extinction at negative latitudes than at positive ones. As for the velocity distributions shown in the upper panels of Fig. 18, we notice that they are quite smooth and do not display evidence of substructures, such as those found for R > 10 kpc towards the Galactic anticentre by Gaia Collaboration (2021c). In terms of (normalised) density distribution, such a difference seems associated with the fact that ∼80% of the stars in our outer region are in the range 10 < R < 11 kpc, while the asymmetries found by Gaia Collaboration (2021c) are more evident for R > 11 kpc (see below for more on this point).

The velocity distributions, V_Z versus V_ϕ, colour coded according to the median value [M/H], are shown in the bottom panels of Fig. 18. As expected, the mean metallicity decreases as a function of Galactocentric distance R (cf. Sect. 4). Also, stars with higher moduli of vertical velocity, |V_Z|, show typically lower [M/H] than stars with smaller vertical velocity dispersions. This effect is consistent with the vertical metallicity gradient.

Moreover, on average, stars with lower metallicity are rotating faster than those with higher metallicity. This is also an expected result that comes from the negative rotation–metallicity gradient of the thin disc stellar population (Kordopatis et al. 2017; Re Fiorentin et al. 2019). This effect is also investigated in Sect. 5.5.

Finally, we return to the issue of possible kinematical substructures, such as those found for R > 10 kpc towards the Galactic anticentre in Gaia Collaboration (2021c). Figure 19 shows the V_Z versus V_ϕ distribution for the Medium Quality sample with 11 < R < 12.5 kpc and |Z|< 1 kpc. From top to bottom, the panels refer to three Galactocentric bins, each of 0.5 kpc in width. Stars above and below the plane are colour-coded according to their density (first two columns) and their metallicity (last two columns), respectively.

We observe a significant difference in the velocity-space density distribution either side of the Galaxy’s midplane. In particular, below the plane (second column), we see two stellar populations that are similar to those found by these latter authors for 11 < R < 13 kpc, as seen in their Fig. 13. Furthermore, we notice a first group of stars with V_Z ≃ 0 and V_ϕ ≃ 220 km s⁻¹ that is evident in the range 11 < R < 11.5 kpc but decreases in the outer bins, R = 11.5 − 12 kpc and R = 12.0 − 12.5 kpc, where a second population of stars with V_Z > 0 and faster rotation, V_ϕ ≃ 240 km s⁻¹, is more significant. This is clear evidence of the bimodal distribution found by Gaia Collaboration (2021c) and confirmed by McMillan et al. (in prep.).

In addition, the chemical information (last column) clearly shows that the population with lower V_ϕ is more metal-rich than the second group of stars with higher V_ϕ and positive V_Z. This is an interesting new result that requires further investigation. This chemo-kinematical signature is confirmed at larger radii, as shown in Fig. E.1 which shows the V_Z versus V_ϕ distribution for a set of radial bins from R = 10 kpc to 14 kpc.

5.3. Vertical velocity oscillations

We focus on the Z versus V_Z distribution, where vertical wave-like motions have been detected in the form of a snail-like pattern especially present when stars are colour-coded according to their azimuthal velocities (Antoja et al. 2018). Here, we analyse the Gaia DR3 GSP-Spec database selecting 3 229 518 stars from the Medium Quality sample (with the above-described additional conditions on the astrometric parameters), requiring that the stars be within |Z|< 1 kpc, |V_Z|< 60 km s⁻¹, and Galactocentric distance range 5.25 < R < 11.25 kpc. This region is approximately centred in the smaller volumes of the solar neighbourhood studied previously by Antoja et al. (2018) and other authors such as Bland-Hawthorn et al. (2019).

Figure 20 shows the projection of phase space in vertical position and velocity, Z versus V_Z. The panels refer to three Galactocentric rings, each of 2 kpc in width, with central radius at R = 6.25, 8.25, 10.25 kpc, respectively. The stars are colour-coded according to their density (top panels), radial and azimuthal velocities (second and third panels from the top, respectively), and residual metallicity (bottom panels).

The density distributions exhibit structures that are reminiscent of the effects of incomplete phase-mixing. Also, the V_R colour-coded distributions show strong spiral substructures over the explored Galactocentric range, especially for R > 7 kpc.

As in Antoja et al. (2018), the distribution colour coded according to V_ϕ for the stars near the Sun (central panel of Fig. 20) follows a curled, spiral-shaped distribution particularly pronounced up to |V_Z|≃40 km s⁻¹. In the left panel, that is, 5.25 < R < 7.25 kpc, we still clearly see a snail-like structure of fast rotating stars, although with a smaller extension in height (|Z|< 0.3 kpc) and a vertical velocity amplitude which remains below ∼30 km s⁻¹ for V_Z > 0 km s⁻¹ and reaches ∼ − 40 km s⁻¹ for V_Z < 0 km s⁻¹.

On the right panel, that is 9.25 < R < 11.25 kpc, stars with V_ϕ ∼ 230 km s⁻¹ follow a spiral-like pattern with a tail in the region at V_Z > 20 km s⁻¹. As expected from dynamical considerations (Binney & Tremaine 2008), and as previously reported in the literature using Gaia DR2 data (Laporte et al. 2019a), the Z-axis elongation increases with Galactic radius. Indeed, at large radii, dynamical timescales are longer due to a weaker self-gravity of the disc, inducing a less tightly wound phase-space spiral.

Furthermore, for the first time, we have the possibility to look for chemical signatures associated with the spiral phase-space structures discussed above. These are shown in the bottom panels of Fig. 20, where the Z versus V_Z distribution is colour coded according to ‘residual’ metallicity Δ[M/H], which is calculated as the difference between individual metallicities and the running median over bins of 100 pc in R. The phase space structures seen in the density and kinematic distributions appear to be linked to [M/H] variations, with more clear correlations in the inner and outer bins. In the inner Galaxy region, more metal-rich stars tend to be below the plane (−0.2 < Z < 0.0 kpc), where the spiral pattern produces a stellar density excess and lower V_ϕ values. In the external regions, stars with a metallicity excess with respect to the median [M/H] values can be found more than 0.5 kpc above the plane, coinciding again with a stellar density excess and lower rotation. Therefore, in both panels, the rotation velocity V_ϕ appears slightly anticorrelated with Δ[M/H], in the sense that stars with higher metallicity have slower rotation and vice versa. This is in agreement with the known V_ϕ − [M/H] correlation in the disc.

However, we note that an important warning has to be given with respect to the inner galaxy pattern, where a significant selection effect due to an inhomogeneous extinction could be biasing both the stellar density and the metallicity distribution. However, this is not the case in the external regions, where the metallicity excess correlation with the phase spiral seems more robust.

Finally, it is worth noting that Bland-Hawthorn et al. (2019) used chemical abundances to separate thin-disc (α-poor) and thick-disc (α-rich) populations, concluding that the phase spiral is confined to the thin disc. In this paper, we explore for the first time the chemical dependencies of the phase spiral within the thin disc at different radii, adding new and more focused constraints.

5.4. Chemo-kinematical ridges in V_ϕ versus R

Using Gaia DR2, diagonal ridge features were identified in the distribution of stars in the V_ϕ versus R space (Antoja et al. 2018; Kawata et al. 2018). Such features might reflect the ongoing phase mixing in the Galactic disc (Fux 2001; Minchev et al. 2009; Gómez et al. 2012), or, alternatively, be a consequence of resonant orbits induced by the bar and spiral arms (e.g. Dehnen 2000; Michtchenko et al. 2018; Fragkoudi et al. 2019; Hunt et al. 2019; Monari et al. 2019a). After their discovery, several properties of the V_ϕ versus R distribution were studied to shed light on their origin and dynamical nature (e.g. Ramos et al. 2018; Fragkoudi et al. 2019; Laporte et al. 2019b; Gaia Collaboration 2021c).

In this section, we map V_ϕ versus R ridges both in phase space and metallicity space using Gaia DR3. Figure 21 shows the diagonal ridges obtained using our sample of RGB stars (see Sect. 3). The results for the RGB stars at positive (negative) Galactic azimuths ϕ are shown in the right (left) column. Positive ϕ values are defined in the direction of Galactic rotation, and as a consequence, are closer to the major axis of the bar. The upper panels show the ridges in the overdensity , where N and are the bivariate kernel density estimators adopted in Appendix B1 of Poggio et al. (2021b) for the local and mean density, respectively. The adopted bandwidths for the local density are 100 pc in R and 2 km s⁻¹ in V_ϕ, while the corresponding values for the mean density are 700 pc and 14 km s⁻¹. The middle panels in Fig. 21 show the median metallicity [M/H], while the bottom panels indicate the residual metallicity, Δ[M/H], with respect to the median metallicity estimated as a function of Galactocentric distance R. Known ridges from Ramos et al. (2018) are overplotted. As we can see, there is a strong correlation between the ridges found in the overdensity and those observed in [M/H] and Δ[M/H]. Such a correlation was pointed out in previous studies (e.g. Khanna et al. 2019). While some ridges do not exhibit an obvious signature in metallicity (e.g. Arcturus), most of them are typically more metal rich (e.g. Hyades, Horn/Dehnen98). Moreover, some features are markedly different at positive and negative ϕ. For instance, the Horn/Dehnen98 ridge appears as a metal-rich feature, which is more prominent at ϕ > 0 than at ϕ < 0. Such variations in the metallicity of those ridges are detected here for the first time, taking advantage of the unprecedented coverage of Gaia DR3 chemical measurements. Furthermore, the residual [α/Fe] (Δ[α/Fe]) in the V_ϕ versus R space is shown in Fig. 22. The azimuthal dependency is not shown here in order to increase the statistics and better identify possible features. As we can see, diagonal ridges are apparent when stars are colour coded according to Δ[α/Fe]. In particular, we note that Hercules and Horn/Dehnen exhibit low [α/Fe].

For comparison, Fig. 23 shows a similar distribution for the sample of massive stars defined in Sect. 3. Ridges in the young populations were identified by Wang et al. (2019) and Gaia Collaboration (2021c) but here we now present them with an alternative and very robust selection of the young population. Due to the low number statistics, stars are plotted all together, without showing the azimuthal variations. While ridges are also present in this sample, we observe that the signal is somewhat noisier than the RGB sample, presumably due to the low number statistics.

The comparison between the ridges found in the young and old stellar populations potentially contains important clues as to their origin. Ridges in young stellar populations can be generated by either a recent perturbation (e.g. an interaction with a satellite galaxy) or an ongoing mechanism in the Galactic disc (e.g. resonances). Additionally, they could be related to star formation structures mixing into the disc, but in this case no match between the ridges in the young and old populations would be expected. Due to the low number of massive stars, it is difficult to establish a robust link between the young and old ridges in our data, although ridges such as Sirius and Arch1 seem to be present in both samples. An alternative scenario is that the global disc dynamics and the star formation processes are working so as to produce the observed similarities in the two samples.

The ridges shown in Figs. 21 and 23 are presumably related to the bumps found in the metallicity gradients as a function of guiding radius (Fig. 15; however, it should be noted that similar bumps can also be caused by the selection function when gradients are plotted as a function of observed radius R, as discussed in Sect. 4).

5.5. Rotation–metallicity correlations

Figure 24 shows a grid of [α/Fe] versus [M/H] chemical planes colour coded according to V_ϕ. The different panels show the distribution across the Milky Way for giants (log(g) < 3.5) of the High Quality sample, with the additional conditions on the astrometric parameters defined at the beginning of this section.

In the panels near the Sun and close to the Galactic plane with 6 < R < 10 kpc and |Z|< 1 kpc, red and orange colours represent fast rotating thin disc stars (V_ϕ ∼ 220 − 250 km s⁻¹) with low [α/Fe]. We notice that the azimuthal velocity of metal-rich thin-disc stars with [M/H]> 0 is slightly lower than for stars with [M/H]< 0. This effect is consistent with the well-known negative rotation–metallicity gradient of the thin-disc population, −20 ≲ ∂V_ϕ/∂[M/H] ≲ − 10 km s⁻¹, as reported by various studies (e.g. Lee et al. 2011; Allende Prieto et al. 2016; Kordopatis et al. 2017; Recio-Blanco et al. 2014) and is associated with blurring effects in the disc (i.e. stars on higher eccentricity visiting a given radius at their apocentre or pericentre (Lynden-Bell & Kalnajs 1972; Sellwood & Binney 2002).

The [α/Fe] abundance from GSP-Spec that mostly traces [Ca/Fe] does not show a clear gap between the thin- and thick-disc sequences as in other chemical elements, such as Mg, which is often adopted in the literature because it is mostly produced by high-mass stars, while Ca is also produced by SNe Ia on longer timescales. Therefore, we derived an approximate estimation of the rotation–metallicity gradient of the thin disc by selecting 763 845 stars close to the Galactic plane up to |Z|< 1 kpc, 8 < R < 10 kpc, −0.4< [M/H]<  + 0.6 dex, and V_ϕ > 100 km s⁻¹. We find a gradient of −13.9 ± 0.1 km s⁻¹ dex⁻¹, which is consistent with the values from the literature.

In the panels with 6 < R < 10 kpc and |Z|< 1 kpc, the low metallicity and α-enhanced stars belonging to the thick disc and halo are shown with green and blue colours that indicate a slower rotation, V_ϕ < 200 km s⁻¹. We notice that the fraction of thin-disc stars decreases significantly in the panels above the plane with 1 < |Z|< 3 kpc. Here, the thick disc is dominant and the colours provide clear evidence of the positive rotation–metallicity gradient of this population (Spagna et al. 2010; Kordopatis et al. 2011, 2013; Recio-Blanco et al. 2014; Lee et al. 2011). By selecting 1535 sources with 2 < |Z|< 3 kpc, 8 < R < 10 kpc, −1< [M/H]<  − 0.2 dex, and V_ϕ > 50 km s⁻¹ to reduce the contamination from thin-disc and halo stars, we measure a gradient of +33 ± 7 km s⁻¹ dex⁻¹ for the thick disc. This result is quite consistent with recent studies based on the chemical classification of thick-disc stars (e.g. Kordopatis et al. 2017; Re Fiorentin et al. 2019) and interpreted as a signature of inside-out formation and gas re-distribution in the primitive disc (Schönrich & McMillan 2017), or as a correlation resulting from the superposition of mono-abundance subpopulations with negative slopes (as in the thin disc, Minchev et al. 2019).

As discussed in Sects. 3 and 5.1, the thin-disc population is incomplete towards the inner disc because of interstellar extinction close to the Galactic plane. Instead, the thick disc can be observed and a positive rotation–metallicity gradient for the thick disc is present in the inner annulus, 4 < R < 6 kpc, as recently detected also by Re Fiorentin et al. (2019). Finally, because of the shorter radial scale length of the thick disc, the outer disc is mainly formed by thin-disc stars that show the typical negative rotation metallicity gradient until out to R = 12 kpc.

6. Chemo-dynamical analysis

Here we explore the relation between the orbital parameters and the stellar chemistry inferred by GSP-Spec using the General sample defined in Sect. 2.5. The orbital parameters of the full GSP-Spec sample were computed using the method described in Sect. 2.4. As the potential considered in the determination of the orbital parameters is stationary and axisymmetric (McMillan 2017), the vertical component of the angular momentum L_Z and the total energy E are integrals of motion; therefore, their values can be directly estimated from the initial conditions. For sake of simplicity, we adopt the convention L_Z > 0 for clockwise rotating systems and rescale the units of L_Z and E in terms of L_⊙ = R_⊙ V_⊙ and , respectively.

6.1. Merger debris in the Galactic halo

Figure 25 illustrates the distribution of the GSP-Spec sample in the E versus L_Z plane. In this diagram, positive energies correspond to stars not gravitationally bound to the Galaxy, while positive (negative) L_Z values correspond to co-rotating (counter-rotating) systems. As seen, this sample is clearly dominated by stars in the solar neighbourhood (L_Z = L_⊙, E ≈ −2.88) with almost circular orbits (E and L_Z close to the border of the forbidden area). By visual inspection, we can identify the overdensity of stars at low |L_Z| and −2.8 ≲ E ≲ −2.0 associated with the Gaia-Enceladus-Sausage (GES; Belokurov et al. 2018; Myeong et al. 2018; Helmi et al. 2018; Feuillet et al. 2020, 2021), as well as an extended structure of counter-rotating stars in the region L_Z ≲ −0.80L_⊙, E ≳ −2.8 likely related to the Sequoia stream (Myeong et al. 2019). Similarly, we infer the presence of the Thamnos structure (Koppelman et al. 2019; Helmi 2020) from the extension of the high-density region towards (L_Z, E) ≈ (−0.50L_⊙, −3.11).

The comparison of Fig. 25 with the E versus L_Z diagram of Fig. 2 in Koppelman et al. (2019) motivates the selection of one additional region associated with the Helmi stream (HStr; Helmi et al. 1999). In Fig. 26, we reproduce the E − L_Z density plot using the same selection criteria as those used by Koppelman et al. (2019) in order to focus on the stars with halo kinematics by imposing |V−V_circ| > 180 km s⁻¹, where V_circ is the circular velocity, and a minimum heliocentric distance of 3 kpc to exclude nearby sources. It is worthwhile noting that this comparison must be interpreted qualitatively due to the different assumptions made for R_⊙, V_⊙, and the values of the solar motion. Therefore, in our Galactic model, the HStr spans an elongated region between (L_Z, E)≈(0.40L_⊙, −2.84) and (0.86L_⊙, −2.30). Finally, we include a control sample of thin disc stars by selecting a circular area of radius 0.15 around the LSR in the E versus L_Z diagram (white circle in Fig. 25). We clean this sample by selecting only the giant stars (0.5 ≤log(g) < 1.8 dex, 4000 ≤T_eff< 8000 K) with good spectral fitting (logχ² < −3.5). For the sake of visualisation, we focus on a specific area that includes the five selected structures and analyse its chemical distribution.

Figure 27 illustrates the global distribution of metallicity (left panel) and [α/Fe] (right panel) of the Medium Quality GRVS sample in the E versus L_Z plane. This figure is complemented by the individual [α/Fe] versus [M/H] diagrams of the four selected structures and shown in the different panels of Fig. 28. As we can see, the Thamnos and Sequoia subsets are dominated by a metal-poor population (⟨[M/H]⟩∼–1.3 dex), with lower [α/Fe] values compared to the disc stars at the same metallicity (density plot in the background, corresponding to the white circle in Figs. 25 and 27). In contrast, the HStr subset shows a significant contamination by thin disc ([M/H] ≲ − 0.5 dex) and thick disc ([M/H]  ≳   − 0.5 dex, [α/Fe]  ≳  0.3 dex) stars due to its proximity to the disc region in the E–L_Z diagram (Fig. 25). However, we can infer a low-α population consistent with the stellar streams of accreted satellites in the metal-poor regime that traces a different trend in the [α/Fe] versus [M/H] diagram from that traced by the disc. Finally, the chemical pattern revealed by the Gaia-Sausage-Enceladus subset clearly differs from that of the disc despite its large dispersion. An in-depth analysis of accretion debris is beyond the scope of this paper. Nevertheless, the metallicity distributions of the low-[α/Fe] populations in the selected halo substructures show differences that could be attributed to the mass of the accreted satellites. In particular, Sequoia candidates have a metallicity distribution peaking around [M/H] ∼ − 1.4 dex for [α/Fe] ∼0, while Gaia-Enceladus-Sausage peaks at [M/H] ∼ − 1.2 dex at the same [α/Fe] values of ∼0. This difference of about 0.2 dex, suggests a more metal-poor knee in the [α/Fe] versus [M/H] distribution for Sequoia than for Gaia-Enceladus-Sausage. This is in agreement with literature results indicating a lower mass for Sequoia than for Gaia-Sausage-Enceladus using ground-based data (Myeong et al. 2019; Feuillet et al. 2020, 2021). In conclusion, Fig. 28 confirms that the Gaia chemical database possesses the required quality to provide chemical diagnostics of accretion through the abundance ratio of α-elements to iron.

6.2. Chemo-dynamics of globular clusters

The globular cluster (GC) system of the Milky Way has historically driven many Galactic and stellar physics studies. For this reason, we focus here on this particular population, illustrating the quality of the DR3 chemo-dynamical data available for these objects. Figure 29 shows the distribution of individual GC stars in the E versus L_Z diagram (coloured symbols). In the left panel, each cluster is denoted by a different colour and symbol, while in the right panel the colour code represents the [M/H] showing the metallicity coherence of the stars belonging to the same cluster. The density plot in the background corresponds to the General sample. Several clusters associated to accretion overdensities like Gaia-Enceladus (Helmi et al. 2018; Myeong et al. 2018) or Sequoia (Myeong et al. 2019) are included. GC members have been selected allowing a radial velocity dispersion of 20 km s⁻¹ with respect to the median value of each GC. It is worth noting that distance uncertainties spread the stars of the same cluster in both E and L_Z, as clearly seen for NGC 3201 inside the Sequoia orange ellipse.

GC stars are faint and are therefore generally in the low-S/N regime of the GSP-Spec sample, and require more deblending corrections of their spectra because of the crowding. For this reason, only a subsample of the clusters included in Fig. 29 have reliable [α/Fe] estimates, reducing the total sample to 14 clusters. Figure 30 shows the distribution of [α/Fe] abundances with respect to [M/H] for the subsample of GCs with the most reliable abundances. Error bars are based on the abundance dispersion within each cluster. We also note that GC stars with better quality spectra are often outside the applicability domain of the calibration polynomials in GSP-Spec. For this reason, we performed a specific calibration with respect to the literature metallicities; an [M/H] offset of 0.1 dex has been corrected). The colour code reflects the GC L_Z median values in solar units. In the background, a density plot with the RGB star sample distribution is presented as a visual reference. Although, the [α/Fe] dispersion is rather high in some cases, which is possibly due to the data quality and/or membership errors, the global [α/Fe] versus [M/H] trend follows the expected distribution, illustrating the potential of Gaia DR3 data for GC studies. Finally, it is important to note that a more thorough analysis focusing on each cluster individually will certainly improve the results. Moreover, in future Gaia data releases, GC will clearly benefit from the increase of the RVS spectra S/N, allowing a deeper study of this important Galactic population.

6.3. Chemo-dynamical signatures in radial actions

The actions (J_R, J_Z, L_Z) in static potentials are integrals of motion that characterise the orbit of the stars. In particular, the action J_R characterises the radial amplitude of the epicyclic orbits while the angular momentum L_Z sets the guiding radius (Binney & Tremaine 2008; Trick et al. 2021), a more robust estimate of the typical Galactic distance of the star than the present-day Galactocentric radius R.

The right panels of Fig. 31 illustrate the distribution of the General Sample in the J_R versus L_Z plane with colour coding reflecting stellar density (first row), the maximum orbital distance from the plane Z_max (second row), the metallicity (third row), and [α/Fe] (last row). We observe that the metallicity distribution follows the expected decreasing trend with Z_max, with the exception of the high angular momentum regime which shows relatively high metallicities. As this sample is biased by nearby stars, most of them have a vertical angular momentum L_Z ≈ L_⊙, the angular momentum of the LSR. These stars show almost circular orbits (J_R ∼ 0) between around 0.92 L_⊙ and 1.10 L_⊙, which implies a guiding radius of between 7.6 and 9.1 kpc in the assumed Galactic potential. However, the low-density region (≲10 counts) includes a wider variety of substructures.

First, the counter-rotating population (L_Z < 0), in which we identify the GC NGC 3201 in the elongated structure at L_Z ≲ −1.25L_⊙, J_R ≲ 0.5L_⊙ (Myeong et al. 2019). These retrograde stars are predominantly metal poor ([M/H]< −1.1 dex) as shown in the literature (Carollo et al. 2007; Morrison et al. 2009; Kordopatis et al. 2020).

Second, a region of almost zero angular momentum that extends towards large values of J_R (≳0.25L_⊙). This area is characterised by stars in highly eccentric orbits with large vertical displacements from the Galactic plane (second row in Fig. 31). The high inclination of these orbits and the low metallicity of the stars in this region (third row) are consistent with the halo population.

Third, an extension of the almost solar metallicity region [M/H] ≳ − 0.5 dex towards L_Z ≳ 1.2L_⊙, J_R ≳ 0.1L_⊙. This area is populated by stars from the outermost part of the Galactic disc (R ≳ 10 kpc) in relatively eccentric orbits ϵ ∼ 0.4. Although these properties are compatible with a radially migrated population from inner regions, further analyses are required to support this interpretation.

In addition to these features, a detailed exploration of the most populated region reveals the presence of ‘ridges’ in the distributions of stars (similar to those reported by Trick et al. 2019), Z_max , and metallicity (right panels in Fig. 31), while some of them have no counterpart in the [α/Fe] map. By visual inspection, we trace the configuration of these ridges (dashed lines) to illustrate the comparison of the maps. We find good agreement between the extension of the denser regions (upper right panel), the almost flat orbits (right panel in the second row), and the metal-rich stars (right panel in the third row). We can identify four prominent ridges (first, second, third, and fifth leftmost dashed lines) with negative slope in the J_R − L_Z plane, each of which can be seen clearly in both the metallicity and Z_max distributions. In particular, the fifth ridge delimits a transition in the metallicity map from [M/H] > –0.2 dex to [M/H]< –0.2 dex. Similarly, we discern a subtle ridge (fourth dashed line) with minor effects on the chemical pattern but significant imprint on Z_max. The sixth ridge corresponds to a relative minimum in metallicity and a change in slope in the iso-density contours. Finally, the seventh ridge indicates a weak increase in metallicity as well as another change of slope in the contour curves. In addition to the ridges, two high-metallicity and low-[α/Fe] areas centred at (L_Z, J_R)≈(0.8L_⊙, 0.003L_⊙) and (0.8L_⊙, 0.025L_⊙) can be identified, as well as a third region at (0.96L_⊙, 0.005L_⊙) with high metallicity but solar [α/Fe]. Although some of the ridges shown in Fig. 31 are similar to those reported by Trick et al. (2019), Trick (2022), who interpret them as signatures of the resonances, we cannot conclude our ridges are produced by the same mechanism, because: (1) we have not explored the influence of the bar and whether it could produce or alter the signatures we see; (2) we do not restrict our study to stars close to the Galactic plane, as in Trick (2022); and (3) some of our ridges have no counterpart in the predictions of Trick (2022). In any case, a complementary and detailed analysis of the orbital frequencies in the ridges would be required to confirm the resonances.

7. Individual abundances in the solar cylinder

We remind the reader that the GSP-Spec module was able to measure the individual abundances of 12 different chemical elements, including iron (Recio-Blanco et al. 2023). In this section, we turn our attention towards the solar cylinder, here defined by a distance to the Sun projected onto the Galactic plane smaller than 500 pc and no vertical cut. In this volume, we examine the distributions of nine abundance ratios, namely [N/Fe], [Mg/Fe], [Si/Fe], [S/Fe], [Ca/Fe], [Ti/Fe], [Cr/Fe], [Ni/Fe], and [Ce/Fe], as a function of mean metallicity [M/H] and their dependence on the maximum vertical distance to the plane (Z_max), the azimuthal velocity (V_ϕ), and the orbital eccentricity (e). [Zr/Fe] and [Nd/Fe] are not presented here, because the number of stars with measured abundances in the solar cylinder is too small.

We use the Individual abundance subsample (as defined in Sect. 2.5), to which we apply some additional filters: astometricparamssolved = 31 OR 95 (stars with respectively 6- or 5-parameter astrometric solutions), ruwe < 1.4, duplicatedsource IS false, and a projected distance to the Sun of less than 500 pc. The surface gravity range depends on the chemical element, as expected from the presence or not of the related spectral features. Figure G.1 shows the Kiel diagrams of the nine individual abundance samples. We also apply the abundance corrections as a function of surface gravity recommended in Recio-Blanco et al. (2023) and which are meant to set dwarf and giant stars on the same scale. Calibrations are provided for nitrogen, α-elements, and iron-peak elements, but not for heavy elements (see Recio-Blanco et al. 2023, for a discussion) and we therefore apply no correction to the cerium abundances. The samples were restricted to the range in surface gravity where the corrections are applicable.

Figures 32 and 33 show the distributions of the abundance ratios as a function of metallicity for [N/Fe], [Mg/Fe], [Si/Fe], [S/Fe], [Ca/Fe], and [Ti/Fe] and for [Cr/Fe], [Ni/Fe], and [Ce/Fe], respectively. Each row corresponds to one element. The colour bars encode: stellar density (col. 1), median maximum vertical distance to the plane (col. 2), median azimuthal velocity (col. 3), and median eccentricity (col. 4), per cell of 0.04 dex in [M/H] and per 0.02 dex in [X/Fe].

Nitrogen is produced through the CNO cycle. In Fig. 32, [N/Fe] shows a flat trend with [M/H]. This is in satisfactory agreement with the modest gradient measured by Ecuvillon et al. (2004), 0.10 ± 0.05, or Suárez-Andrés et al. (2016), 0.04 ± 0.07 and 0.13 ± 0.04 (respectively, for their control sample and for a sample of stars with detected planets). It should be noted that da Silva et al. (2015) measured a stronger gradient of 0.28 ± 0.06. The proportional enrichment in nitrogen and iron abundances supports the scenario in which nitrogen is produced in intermediate- and low-mass stars, rather than in rotating massive stars (see e.g. the discussion in Ecuvillon et al. 2004). In the solar cylinder sample, nitrogen is mainly measured in giants, including a small proportion of hot-massive giants. Their nitrogen surface abundances may have been enhanced during the dredge-up phases (Masseron & Gilmore 2015; Masseron et al. 2017; Shetrone et al. 2019).

We remind the reader that magnesium, silicon, sulfur, calcium, and titanium are α-elements. It can be noted that silicon, calcium, and titanium, which have several lines in the RVS wavelength range, have been measured in large sample of stars in the solar cylinder: from more than 80 000 for [Ti/Fe] to more than 140 000 for [Si/Fe]. As shown in Fig. 32, their abundance ratios decrease with metallicity consistently with previous studies in the solar neighbourhood (e.g. among others Adibekyan et al. 2012; Recio-Blanco et al. 2014; Bensby et al. 2003, 2014). This trend is produced by the progressive enrichment in iron of the interstellar medium by thermonuclear supernovae and is well reproduced by Galactic chemical evolution models (Snaith et al. 2022; Prantzos et al. 2018; Kobayashi et al. 2020; Matteucci 2021). As already discussed in Sect. 3, we again point out the continuous decrease in α-element abundances up to supersolar metallicities, as predicted by Galactic evolution models and already revealed by previous studies based on high-resolution spectra (see e.g. Santos-Peral et al. 2020; Perdigon et al. 2021). A sparsely populated feature is also visible in [Mg/Fe], [Ca/Fe], and [Ti/Fe] in the range [M/H] ∈[−0.5, 0] dex and at low abundance ratio values. It is made of a mix of dwarfs and giants in Mg and Ti and of giants in Ca (by construction of the sample). The Massive giants (as defined in Sect. 3) are contained within the feature, although they are not the unique contributors to it. This is due to the partial overlap between the Massive sample and the contiguous RGB stars in the Kiel diagram, as also seen in Fig. D.1. The relative contribution of the two populations varies from one element to another and this is why the median values of V_ϕ and eccentricity in Fig. 32 are slightly different from one element to another at that abundance locus. It is worth noting that the α-element enhancement and the metallicity compete to regulate the opacity of the stellar atmosphere. Both a higher [α/Fe] and a higher [M/H] produce a higher opacity and therefore a cooler effective temperature. As a consequence, low-[α/Fe] stars are hotter than high-[α/Fe] stars at the same metallicity, but they overlap in T_eff with lower metallicity stars. This is what happens in the Kiel diagram region at the limit of the Massive and RGB samples (see e.g. Fig. 22 in Recio-Blanco et al. 2023).

Chromium and nickel are iron-peak elements. As such, their abundance ratios are weakly dependent on metallicity, although some authors observe a slow rise in [Ni/Fe] at metallicities beyond solar (e.g. Adibekyan et al. 2012; Bensby et al. 2014). In Fig. 33, as expected, [Cr/Fe] abundances show no trend with metallicity. [Ni/Fe] shows a flat trend on the metal-rich side, while on the metal-poor side, the abundance ratios present a small positive offset.

Cerium is a heavy element mostly produced by the s-process (see Sect. 3). In Fig. 33, the bulk of solar cylinder stars show no significant trend with [M/H]. This is consistent with earlier works using high-resolution spectra (Battistini & Bensby 2016; Delgado Mena et al. 2017; Forsberg et al. 2019). In Fig. 33, the median cerium-to-iron-abundance ratio is around +0.2 dex. This is close to that measured by Forsberg et al. (2019) (for stars with [M/H] < 0 dex), but higher than Battistini & Bensby (2016) and Delgado Mena et al. (2017) whose average [Ce/Fe] level is around 0 dex or slightly below.

As discussed in Sect. 3, the possibility of estimating the abundance of a given chemical element depends on the presence of the lines of that element in the related spectra, and therefore it changes with stellar type and luminosity (see also Fig. G.1). As a consequence, the tracers of different chemical elements probe different ranges of maximum distance to the Galactic plane (Z_max). Eight out of the nine elements presented here (Na, Mg, Si, Ca, Ti, Cr, Ni and Ce) can be measured in giant stars. As shown in the second column of Figs. 32 and 33, on the metal-poor side, the median (Z_max) can reach up to 1 to 2 kpc. As already discussed for Fig. 9, sulfur is measured in hot F-G dwarf and turn-off stars from the RVS spectra (Recio-Blanco et al. 2023). It therefore probes a narrower range of Z_maxwith median (Z_max) being found at around a few hundred parsecs at all metallicities.

In the solar cylinder, the larger the range of Z_max probed, the more the thick disc is sampled. At low metallicity, the thick disc is clearly visible in the α-elements Mg, Si, Ca, and Ti as a lower median azimuthal velocity and higher median eccentricity structure (green/blue cells in cols. 3 and 4). At higher metallicity, the precision of the abundances does not allow us to see the gap between the high- and low-α sequences in stellar density space (col. 1). Nevertheless, in the azimuthal velocity space (col. 3), looking at Mg and Ti, the high-α sequence rotates a few tens of km s⁻¹ slower than the low-α sequence, as expected from previous studies (e.g. Adibekyan et al. 2013; Haywood et al. 2013; Recio-Blanco et al. 2014; Allende Prieto et al. 2016; Kordopatis et al. 2017; Re Fiorentin et al. 2019, and Sect. 5.5). This can be used to separate the thick and the thin disc in chemical space for these two elements. In the solar cylinder sample, sulfur is dominated by thin-disc stars where most hot dwarfs lie. It is interesting to note that the thin-disc negative azimuthal velocity gradient as a function of metallicity reported by Adibekyan et al. (2013), Haywood et al. (2013), Allende Prieto et al. (2016), Kordopatis et al. (2017), Re Fiorentin et al. (2019) for example and also measured in Sect. 5.5 is visible in the third column as a transition from red to yellow-green (from metal-poor to metal-rich).

In summary, the Gaia solar cylinder individual abundance distributions reveal the expected chemo-dynamical patterns of thin and thick disc stellar populations for a variety of chemical elements, highlighting the richness of the Gaia chemical space.

8. Chemo-kinematical cartography of the open clusters population

Open stellar clusters are groups of stars that formed together at the same time and from the same material. For this reason, their ages, distances, kinematics, and chemical patterns can be estimated easily and with very high precision. Given these properties, open clusters are widely considered excellent tracers of the Galactic disc, both in space and time (e.g. Friel 1995; Yong et al. 2012; Cunha et al. 2016; Donor et al. 2020; Casamiquela et al. 2021; Spina et al. 2021). In this section, we investigate the distribution of metals ([M/H]) traced by open clusters across the Milky Way disc.

To date, around 2600 open clusters spread throughout the Milky Way disc have been identified in the Gaia data set. Cluster ages, parallaxes, proper motions, and radial velocities have been derived by Cantat-Gaudin et al. (2020), Castro-Ginard et al. (2022), Tarricq et al. (2021) from stars with membership probability ≥0.7. These parameters have then been employed to derive the Galactocentric positions, velocities, and orbital parameters for 2162 open clusters through the same procedure described in Sect. 2.4. It was not possible to derive these quantities for the full sample of known open clusters because a few hundred of them have no radial velocity estimates.

As a result of the cluster membership analysis based on the Gaia astrometric solutions, Cantat-Gaudin et al. (2020), Castro-Ginard et al. (2022), Tarricq et al. (2021) also provide the list of stars that are possibly associated to these open clusters along with their membership probabilities. Of all these stars with membership probability ≥0.7, 6861 have mh_gspspec determinations. Among these latter, we select only those belonging to the High-Quality subsample (Sect. 2.5). From this sample, we further discard all sources with teff_gspspec ≥ 8000 and rv_expected_sig_to_noise ≥ 40. These two additional conditions further improve the quality of our final data set which is composed of 1613 stars. These objects are members of 597 open clusters for which we derive median metallicities [M/H] and relative uncertainties Δ [M/H]. This latter is defined for each open cluster as the standard deviation of the mh_gspspec values of its members. These open clusters have a number of members with GSP-Spec [M/H] ranging between 1 and 55, with an average of three members. When only one member is available, Δ[M/H] is assumed to be equal to the [M/H]_unc of the single star.

8.1. Galactic radial metallicity and [α/Fe] gradients

Here we use the data set described above to investigate the radial metallicity gradient traced by open clusters. The top panel of Fig. 34 shows the [M/H] abundance of the 503 open clusters older than 100 Myr and with R ≤ 12 kpc.

The cut in age was necessary because the spectroscopic analysis of very young stars is considerably complicated by high rotational velocities and by the effects of chromospheric activity which make stars appear more metal-poor than they are in reality (Yana Galarza et al. 2019; Baratella et al. 2020; Spina et al. 2020; Zhang et al. 2021). Instead, the cut in R allows us to include only the open clusters located before the break in the radial metallicity profile which is located at R ∼ 12−14 kpc (for a recent study, see e.g. Donor et al. 2020, and literature therein). As our sample does not include open clusters located at R > 14 kpc, the analysis of this section is focused on the metallicity distribution of the inner disc. Symbols in the top panel of Fig. 34 are colour coded as a function of age. The open cluster distribution is modelled with a Bayesian regression using the same procedure adopted by Spina et al. (2021). Namely, we employ a linear model y_i = α × x_i + β, where x_i and y_i are normal distributions centred on the R and [M/H] values of the ith cluster. The x_i normal distributions have standard deviations equal to . The ϵ parameter accounts for the intrinsic chemical scatter between clusters at fixed Galactic radius. The simulation is run with 10 000 samples, half of which are used for burn in, and a No-U-Turn Sampler (Hoffman & Gelman 2011). The script is written in Python using the pymc3 package (Salvatier et al. 2016). The 68% and 95% confidence intervals of the models resulting from the posteriors are represented in Fig. 34 with grey shaded areas. The results are also listed in Table 3, which reports the mean, standard deviation, and 95% confidence interval of each posterior distribution¹¹.

Interestingly, the radial metallicity gradient α = −0.054 ± 0.008 dex kpc⁻¹ derived from our sample of 503 open clusters is consistent with the value given by Casamiquela et al. (2019), that is −0.056 ± 0.011 dex kpc⁻¹ (18 open clusters), and that of Spina et al. (2022), namely −0.064 ± 0.007 dex kpc⁻¹ (175 open clusters). On the other hand, our slope estimation is higher than those provided by other recent works in the literature: for example −0.10 ± 0.02 dex kpc⁻¹ given by Jacobson et al. (2016) for 12 open clusters; −0.086 ± 0.009 by Netopil et al. (2016) for 88 open clusters); −0.077 ± 0.007 dex kpc⁻¹ by Carrera et al. (2019) for 90 open clusters; −0.068 ± 0.004 dex kpc⁻¹ by Donor et al. (2020) for 71 open clusters; −0.076 ± 0.009 dex kpc⁻¹ by Spina et al. (2021) for 134 open clusters; and −0.066 ± 0.005 dex kpc⁻¹ by Zhang et al. (2021) for 157 open clusters.

The great diversity of these estimations of the radial metallicity gradient highlights how studies based on open clusters are strongly affected by the samples employed in their analysis. In this regard, our study represents a significant step forward in the field as our data set is a few times larger than those employed in the past.

The full analysis is repeated for the distribution of open clusters in the [M/H]–r_guid plane, where r_guid is the guiding radius defined in Sect. 4 as (R_apo + R_peri)/2. The resulting posteriors are listed in Table 3 and are consistent within the uncertainties with those obtained from the [M/H]–R_Gal distribution. This is due to the fact that, with some rare exception, the vast majority of the open clusters in our sample are on almost circular orbits.

Finally, we analysed the radial distribution of [α/Fe] traced by open clusters. The resulting gradient is slightly positive, in agreement with recent studies (e.g. Donor et al. 2020; Spina et al. 2021). The results are plotted in the bottom panel of Fig. 34. The posteriors are listed in Table 3.

8.2. Temporal evolution of the Galactic radial metallicity gradients

The temporal evolution of Galactic metallicity gradients provides information about how the interstellar medium has been progressively enriched in metals, but also about how the metallicity tracers migrate across the Galactic disc. For instance, it is well known that the radial metallicity gradient traced by field stars gets flatter for older stars (e.g. Casagrande et al. 2011; Anders et al. 2017). This behaviour is due to radial migration which is expected to smooth the gradient with time (e.g. Schönrich & Binney 2009; Minchev et al. 2018). In contrast to what is observed for field stars, several works have shown that the radial metallicity gradient traced by open clusters flattens with decreasing age (e.g. Friel et al. 2002; Donor et al. 2020; Spina et al. 2021, 2022). However, the literature is not in complete consensus (see Salaris et al. 2004; Netopil et al. 2022).

We split the sample of open clusters into four groups depending on their ages: age < 1 Gyr, 1 ≤ age < 2, 2 ≤ age < 3, and age ≥3 Gyr. These bins are arbitrarily chosen to contain a balanced number of clusters. Figure 35 shows the metallicity gradients calculated with the same method outlined above. Based on our analysis, we confirm that both the [M/H]–R and [M/H]–r_guid gradients evolve with age in the same way, with the younger cluster gradients being flatter.

Interestingly, literature results for field stars seem to indicate an opposite temporal evolution of the gradient between open clusters and field stars (e.g. Casagrande et al. 2011; Anders et al. 2017; Minchev et al. 2018). This opposite behaviour could be due to a bias imposed by the Galaxy on the current open cluster population (Anders et al. 2017; Spina et al. 2021). In fact, while stars are free to migrate within the Galactic disc, the metal-rich clusters formed in the inner disc can survive only if they migrate outward where the Galactic potentials (e.g. spirals, bar, giant molecular clouds) are weaker and less destructive. In contrast, the metal-rich clusters migrating inward are quickly disrupted. The direct consequence of this Galactic bias is that the gradient traced by old clusters is steeper than the one seen for young clusters and field stars.

Nevertheless, our analysis of the young field stars in the Massive sample presented in Sect. 4 suggests a flattening of the gradient in recent times (a slope of −0.036 ± 0.002 dex kpc⁻¹ is obtained for the Massive sample, while it is of −0.055 ± 0.007 dex kpc⁻¹ for the entire field population considered in the Gradient sample). Therefore, from purely Gaia data, both open clusters and field stars seem to show a flattening of their radial metallicity gradient in the more recent epochs of Milky Way evolution. However, a careful analysis of stellar ages is necessary to put this result in robust context and to allow a significant comparison with the temporal evolution of the gradient traced by open clusters.

8.3. Kinematic features in the radial metallicity distribution

Similarly to what is observed for field stars, it is expected that the complex dynamical substructure of the Galactic disc discussed in Sect. 5.1 would produce visible signatures in the chemical distribution of elements traced by open clusters. In fact, based on a sample of 1382 open clusters, Tarricq et al. (2021) recently showed that they seem to follow the main diagonal ridges in the Galactic radial distribution of azimuthal velocities outlined by field stars.

The top panel of Fig. 36 shows the distribution of open clusters in the V_ϕ–R_Gal diagram (grey dots). Here we consider only the 322 open clusters with |ϕ − ϕ_⊙| < 10° older than 100 Gyr and with a metallicity estimation. We also plot the kernel density estimate of the open cluster distribution (in red) and the ridges (black lines) identified for field stars by Ramos et al. (2018). Similarly to what is observed by Tarricq et al. (2021), here we notice a steep decrease in the number of open clusters on the right of the ‘Sirius’ ridge (solid line). This could be due to a number of factors, such as the steep decrease in the number of open clusters beyond 9.5–10 kpc, and the fact that a significant fraction of our outermost clusters are older than a few gigayears in contrast to the rest of the sample. In any case, other than that, the density features outlined by open clusters do not perfectly match the locations of the other ridges. This may be related to the complex and unexplored selection function of the known open cluster sample which is presumably different from that of field stars. Nevertheless, what we want to highlight here is that the distribution of open clusters shows complicated patterns characterised by different overdensities and a V_ϕ range which is very restricted compared to that of field stars.

The complexity of this pattern is also shown in the middle panel of Fig. 36, where the solid orange line represents the locally weighted scatterplot smoothing (LOWESS) regression of the V_ϕ of clusters as a function of R_Gal. The shaded area shows the 90% confidence interval of the LOWESS functions obtained by sampling within the uncertainties of the two variables. As expected, the innermost and the outermost clusters have the highest and the lowest V_ϕ values, respectively. However, in between these extremes, the LOWESS function shows a noticeably wavy pattern which is particularly evident within 7 and 9.5 kpc. A similar wavy pattern was also noticed by Tarricq et al. (2021) who identified a significant dip at 7.5 kpc in the V_ϕ–R diagram with a departure of more than 2σ from the Galactic rotation curve derived analytically from a simple axisymmetric potential.

Interestingly, this wavy pattern seems to be linked to the metallicity distribution as shown in the bottom panel of Fig. 36. Here, we plot the LOWESS regression of [M/H]–[M/H]_grad as a function of R_Gal, where [M/H]_grad is the metallicity predicted by the linear model of Fig. 35 (top panel). As in the V_ϕ–R_Gal diagram, this LOWESS function also shows a wavy behaviour with the residual metallicity [M/H]–[M/H]_grad swinging around the zero value. We note that each crest in the [M/H]–[M/H]_grad versus R diagram is located near a trough in the V_ϕ–R_Gal plot.

This relation between the waves seen in the two LOWESS functions can be explained as the natural consequence of open cluster migration across the Milky Way disc. More specifically, if we consider only the type of migration that conserves angular momentum (i.e. blurring), the open clusters coming from the outer disc have higher V_ϕ and lower residual metallicity than adjacent clusters at a given R. On the other hand, open clusters coming from the inner disc have lower V_ϕ and higher residual metallicity than adjacent clusters.

However, we still need to explain the origin of the wavy V_ϕ–R distribution in the first place. It is likely that the V_ϕ oscillations are somehow linked to the orbital resonances and perturbations observed from field stars. Although the top panel of Fig. 36 does not provide strong evidence in that sense, the waves are likely the consequence of the same process that produced the diagonal ridges. A deeper analysis of the open cluster population, of their orbital properties, and knowledge of their selection function are necessary to provide additional clues on this topic.

9. Summary: The Gaia chemical cartography of the Milky Way

The previous sections reveal a number of chemo-dynamical features observed in the RVS GSP-Spec database. Some of these are the same underlying Galactic characteristics seen from different perspectives; this section summarises these characteristics, providing a more global view.

9.1. Chemical markers of Milky Way structure

The large radial and vertical coverage of the database, together with its substantial azimuthal coverage, allow us to explore the detail of the structure of the Milky Way within an unprecedentedly large volume. Sections 3–5 explore the Galaxy’s spatial symmetry through chemo-kinematical considerations. While radial and azimuthal substructures are detected (cf. Sect. 9.2), the Galaxy shows an important degree of vertical symmetry (agreement within 0.05 dex in [M/H], see for instance Figs. 13 and 14) in the region of radial distances of between 5 and 10 kpc from the centre, and within 2 kpc of the plane. Outside this region, selection functions and statistical fluctuations perturb the diagnostics. This vertical symmetry imposes an important constraint on simulations, in the context of the rich accretion history of the Milky Way. In particular, further detailed studies of chemo-dynamical correlations including halo stars would help us to explore the possible existence of a preferential axis of accretion, extending recent literature analyses hampered by low number statistics (e.g. Recio-Blanco et al. 2021).

Moreover, stars with typical chemo-dynamical characteristics of the disc population (fast rotation associated with metallicities higher than ∼–0.75 dex) are detected at large distances above and below the plane. The maximum height at which disc stars are detected grows with increasing radial distance. This behaviour is typical of a flared disc, as already detected in the Milky Way for both the stellar population and the gas (e.g. Grabelsky et al. 1987; Feast et al. 2014; López-Corredoira & Molgó 2014; Momany et al. 2006; Reylé et al. 2009; Li et al. 2019; Mosenkov et al. 2021). In the external regions, the disc population is detected up to 4 kpc from the plane (cf. Fig. E.1), suggesting a higher scale height of the flare than previously reported. We note however that an exact quantification of the disc shape is out of the scope of this paper. In addition, a certain contribution of the disc warp (e.g. Poggio et al. 2020; Reylé et al. 2009; López-Corredoira et al. 2020) in the external regions cannot be excluded.

Finally, the identification of massive stars (blue loop and other variable objects) in the GSP-Spec Kiel diagram has allowed us to detect the disc spiral structure. The spiral pattern traced by this population is in agreement with Poggio et al. (2021b), who used upper main sequence stars, open clusters, and classical Cepheids. A more detailed analysis of the spiral arms can be found in Gaia Collaboration (2023a). In the present paper, the metallicity distribution of the population of massive stars shows the expected radial gradient, although azimuthal variations also exist. In particular, at the solar Galactic radius, we find that stars closer in azimuth than the Sun to the major axis of the bar have higher metallicities (Sect. 5, Fig. 11). A detailed characterisation of these chemical inhomogeneities and their link to the spiral arms remains to be performed.

9.2. Chemical markers of disc kinematic disturbances

One of the most evident results of this analysis is that the known kinematic disturbances of the disc, including different manifestations of phase space correlations and kinematic substructure, are clearly associated with chemical patterns, acting as markers of the perturbations. From the observational point of view, it is important to highlight that chemistry allows a completely independent diagnostic from the kinematical one.

First, the V_Z versus V_ϕ distribution in the outer Galactic disc shows density substructures associated to vertical asymmetries and metallicity patterns (cf. Sect. 5.2). Second, the phase spiral in the Z versus V_Z plane (Antoja et al. 2018), already associated with V_R and V_ϕ patterns, appears to be linked to [M/H] variations (Fig. 20, bottom panels). In particular, in the external region, R > 9.25 kpc, where no significant bias from differential extinction exists, stars with metallicity excess with respect to the median [M/H] values clearly trace the observed phase-spiral pattern. Stars with metallicity excess are found more than 0.5 kpc above the plane, coinciding with stellar overdensity and lower rotation values (which seem to follow the expected metallicity trend with mean azimuthal velocities in the disc).

Third, the identified ridges in the V_ϕ versus R space are typically more metal-rich than their phase space surroundings and show azimuthal dependencies (Fig. 21). In addition, although the lower number statistics perturb the interpretation of results, substructures in the V_ϕ versus R plane are also seen in the young field stars tracing the spiral structure (Fig. 23) and in the open cluster population (Fig. 36), possibly associated to metallicity fluctuations.

Additional evidence of chemical fluctuations associated with the seismology of the Milky Way comes, for the first time, from the orbital space. The action plane defined by J_R versus L_Z (Fig. 31) shows clear ridges of higher stellar density characterised by stars in low Z_max orbits with metallicities higher than the surrounding median values.

In summary, these different diagnostics indicate that stars tracing disc kinematic disturbances have higher metallicities. This is true both for phase space correlation patterns and kinematic substructure. Although results from ground-based surveys (Khanna et al. 2019) already identified solar metallicity values in the local V_ϕ versus R ridges, this large-scale analysis of Gaia GSP-Spec data shows that stars in kinematic substructures are more metal-rich than the typical metallicity at a given radius, irrespective of the metallicity regime. The best chemical marker of disc kinematic disturbances is in fact the differential metallicity of the stars with respect to the median [M/H] value at each Galactic radius (as provided by the radial metallicity gradient).

To explain the higher metallicities tracing phase-spiral patterns and ridges, a younger age or an inner Galaxy origin could be invoked. Younger stars are indeed formed from a more chemically evolved (more metal-rich) interstellar medium. In addition, they are expected to be dynamically cooler, and therefore react more strongly than older dynamically hotter components when subject to either internal perturbations (Fragkoudi et al. 2019) or to satellite impacts (e.g. Chequers et al. 2018). A literature analysis using Gaia DR2 data and a heterogeneous compilation of stellar ages (Laporte et al. 2020) report slightly younger ages for stars in kinematical substructures. Extending previous results of Wang et al. (2019) and Gaia Collaboration (2021c), we show that objects formed only a few hundred million years ago (such as the massive population defined in Sect. 3 tracing the spiral arms or the open cluster population), also show kinematical substructure. On one hand, this is in agreement with N-body simulations of a recent (500–800 Myr, e.g. Laporte et al. 2019a) or a 1–2 Gyr old (e.g. Bland-Hawthorn & Tepper-García 2021) and repeated interaction of the Milky Way with a massive, Sagittarius-like satellite. On the other hand, internal mechanisms of disc perturbation such as the buckling of a stellar bar can also generate long-living phase-spirals (Khoperskov et al. 2020) affecting stars of different ages, in agreement with observations.

It is also worth mentioning that the observed azimuthal dependencies in kinematic substructure (cf. Fig. 21) already detected by previous studies for some ridges (e.g. Ramos et al. 2018; Monari et al. 2019b, for Hyades and Hercules) imply azimuthal metallicity variations due to the corresponding chemo-kinematical correlation. Moreover, azimuthal dependencies are also present in the metallicity distribution of disc RGB and young stars (cf. Sect. 3), with stars in azimuthal angles closer to the major axis of the bar being slightly more metal rich. Recent studies of the impact of Galactic disc asymmetries on azimuthal chemical variations (Spitoni et al. 2019) predict stronger azimuthal fluctuations near the spiral arm co-rotation radius. A detailed spatial analysis of chemical fluctuations is nevertheless out of the scope of this performance verification paper.

9.3. Chemical markers of satellite accretion

The large spatial coverage and large number statistics of the Gaia GSP-Spec database allow a chemo-dynamical study of halo substructure (cf. Sect. 6). From this analysis, which is not intended to be exhaustive, three main conclusions emerge.

First of all, halo substructure in the orbital space (cf. Figs. 26 and 31) is clearly observed, extending previous observations (cf. Koppelman et al. 2019) using Gaia DR2 and ground-based data. Second, the new Gaia chemical data are of sufficient quality to provide chemical diagnostics of accretion through the abundance ratio of α-elements to iron (cf. Fig. 28). This opens new horizons for in-depth analysis of accretion debris and stellar streams. In addition, the present study clearly detects stars with chemical patterns typical of the disc, but on halo orbits. The number of these puffed-up disc stars increases as we move towards the disc region in the energy versus L_Z plane. Finally, the chemo-dynamical properties of the Milky Way system of GCs can also be explored.

9.4. Chemo-dynamics of the most recent billion years

The Gaia chemical database includes 687 open clusters (cf. Sect. 8) and about 30 000 young stars tracing the spiral arms (cf. Sect. 3). These encode the evolution of the last billion years, anchoring the long Galactic history to the present.

On one hand, age estimates of Galactic open clusters allow us to assess the temporal evolution of the metallicity radial gradient of the cluster. Flatter gradients are progressively observed in the younger open cluster populations (cf. Fig. 35), in contrast to the opposite trend reported in the literature for field stars (Casagrande et al. 2011; Anders et al. 2017; Minchev et al. 2018). This could be due to differential effects of radial migration between the clusters and the field population (e.g. Anders et al. 2017; Spina et al. 2021).

On the other hand, interestingly, the young field population that is traced here by massive F- and G-type stars shows a chemical impoverishment with respect to older populations. This impoverishment is observed not only in the global metallicity, but also in the abundance ratio of several chemical species with respect to iron, including α-elements (cf. Fig. 11) and heavy elements like Cerium (cf. Fig. 12). Explaining this metallicity decrease requires a careful study with chemical evolution models. However, its link to the coeval open cluster population could be key to solving the chemo-dynamical puzzle of the most recent epochs of Galactic evolution.

9.5. The Sun in its Galactic environment

The results discussed above bring to the fore the question of the position of our Sun in the global chemo-dynamical evolution of the Milky Way. Spatially, the solar position in the local void outside the spiral arms is clearly visible in Fig. 11. From a chemical point of view, when the Sun is considered within the global stellar population in the solar cylinder (cf. Figs. 32 and 33), it has a chemical pattern and chemo-dynamical characteristics typical of a thin disc star. Interestingly, it has a higher metallicity than the corresponding median value at its Galactic position (assumed here to be R = 8.249 kpc) from the measured radial gradient (Δ[M/H] = 0.04 dex from the Gradient analysis sample and Δ[M/H] = 0.16 dex from the more homogeneous RGB sample). In addition, despite being born 4.6 Gyr ago, the Sun is also more chemically enriched than the young stars in the solar neighbourhood, even if they can be more than ∼4 Gyr younger. These two facts seem to invoke either a discontinuity in the chemical evolution or a change in the Sun’s orbital parameters (migration) or both phenomena at the same time.

From a kinematical point of view, the solar V_Z and Z values place the Sun near the central grid point of the phase spiral (medium panels Fig. 20). In addition, the Sun has a negative V_R velocity and a positive Δ[M/H], which is coherent with the typical values observed in its phase space surroundings. Moreover, as indicated by the solar position in the V_ϕ versus R plane (identified by a cross in Fig. 21), the Sun seems to lay near the Sirius ridge. In the action space, the solar position near (1,0) in Fig. 31 (in units of L_⊙), together with its low Z_max value of 0.12 kpc, again place the Sun near one of the identified ridges (the fifth one from the left in Fig. 31). The above information will be crucial in determining whether or not the Sun has been affected by the different perturbations at play.

10. Conclusions

This article explores the scientific quality of the Gaia chemical database produced by the GSP-Spec module from RVS spectra and published as part of the third Gaia data release (DR3). Combined with Gaia EDR3 astrometric data and DR3 radial velocities, the all-sky Gaia chemical cartography allows a powerful and precise chemo-dynamical view of the Milky Way with unprecedented spatial coverage.

In this framework, chemical abundances carry pivotal information on the origins of the stars. With the exception of some elements that suffer from internal mixing in certain stellar evolutionary phases, the chemical composition of a star shows a preserved pattern of element abundances encoding the conditions of the interstellar medium at the time and place of its formation. Therefore, chemical markers are crucial for validating kinematical and dynamical diagnostics, providing constraints on the properties of matter and their temporal evolution in theoretical modelling and analysis.

In addition, understanding the history of the Milky Way requires the exploration of a complex physico-chemical space that only huge data sets such as the Gaia catalogue can appropriately probe. Large and precise databases not only improve the confidence level of the derived general trends, but also allow researchers to reveal the crucial but rare non-standard populations in distribution tails. Adequately resolving the natural scatter in chemo-dynamical properties of Galactic stellar populations is essential for uncovering the imprint of physical processes of evolution. This scatter is indeed the consequence of star and gas dynamics, of discontinuous star formation histories, and of galaxy accretion, among other mechanisms.

The results of the present analysis illustrate how this chemo-dynamical distribution tail allows us (i) to unveil the pronouncedly flared structure of the disc, (ii) to robustly assess the Galaxy vertical symmetry, (iii) to chemically trace the large-scale motions and seismology of the disc, (iv) to uncover galaxy accretion debris and heated disc stars on halo orbits, and (v) to anchor the Galactic history to its present through the young but scarce and ephemeral massive stars and open clusters.

Gaia DR3 chemo-dynamical diagnostics open new horizons before the era of ground-based wide-field high-resolution spectroscopic surveys like WEAVE (Dalton et al. 2020) or 4MOST (de Jong et al. 2019) which will provide detailed abundances for millions of Milky Way stars. Moreover, the Gaia satellite allows all-sky, continuous data collection over timescales measurable in years, which is impossible to achieve from the ground. This results in an extraordinary homogeneity of data, and allows Gaia to include detailed bias modelling in its analysis and data-treatment procedures, which is at the core of its precision and observational success. A beautiful illustration of this are the Gaia DR3 heavy element abundances, which at first appeared to be out of reach for the RVS instrumental characteristics.

Finally, the results of the present work reinforce the vision of a heterogeneous Milky Way in permanent interaction with its environment, the outcome of an eventful history that continues to shape our Galaxy to the present day.

---

Acknowledgments

This work presents results from the European Space Agency (ESA) space mission Gaia (https://www.cosmos.esa.int/gaia). Gaia data are processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia archive website is https://archives.esac.esa.int/gaia. Acknowledgments from the financial institutions are given in Appendix H. Part of the calculations used for this article have been performed with the high-performance computing facility SIGAMM, hosted by the Observatoire de la Côte d’Azur.

References

Appendix A: The astrometric control sample

The astrometric control sample (ACS) is a subset of GSP-Spec whose distances are derived directly from the Gaia DR3 parallaxes, ϖ, after correcting for the individual parallax zero-point, Δ_ϖ, calculated according to Lindegren et al. (2021), that is, using the corrected parallaxes

ϖ′=ϖ − Δ_ϖ, and after including only sources with high-quality astrometric parameters, full three-band photometry, and radial velocities. In short, the ACS selection criteria are:

astrometric_params_solved = 31 AND
ruwes < 1.4 AND
parallax_over_error > 5 AND
duplicated_source IS false AND
radial_velocity IS NOT null  AND
phot_g_mean_mag IS NOT null AND
phot_bp_mean_mag IS NOT null AND
phot_rp_mean_mag IS NOT null

We note that astrometric₋params₋solved = 31 means that only five-parameter astrometric solutions are retained, and the criteria parallax₋over₋error > 5 immediately eliminates all negative parallaxes. This has no statistical consequence given the nature of the ACS sample and its use in this work. Trigonometric distances are then

(A.1)

and, as is often the case in the astronomical literature, the lower and upper distance errors are, respectively:

(A.2)

where α = σ_ϖ/ϖ is the parallax relative error. These distance errors correspond to the 1-σ parallax errors σ_ϖ with the asymmetry, , decreasing with α.

Table A.1 reports the basic statistics of the complete ACS and of the subsets ACS_MQ and ACS_HQ, that is, the lists with only Medium and High Quality samples, respectively (see Sect. 2.5). The full ACS includes 4 397 814 stars of the General sample, about half of which (53%) are astrometric primary flagged (astrometric₋primary₋flag = True). It reaches distances of up to ∼ 10 kpc from the Sun with a median of 1.3 kpc and 90% of its stars within ≃ 3.7 kpc. The parallax relative error α is below 7% for the 90% quantile. The ACS_MQ includes 3 362 012 stars up to 10 kpc with median distance r_ϖ = 1.5 kpc and mostly (90%) within ≃ 2.9 kpc. For this sample, α is below 5% for 90% of the sources. Finally, the ACS_HQ is formed by 1 815 072 stars up to 6 kpc with median distance r_ϖ = 1.5 kpc, mostly (90%) within ≃ 3.4 kpc, and α under 6% for 90% of the cases.

Table A.1 also includes the analogous statistics for the corresponding samples with Geo distances, r_Geo, from Bailer-Jones et al. (2021). These Geo samples contain sources with five- and six-parameter astrometric solutions, a quality index RUWE < 1.4, and no duplicated sources. The full Geo sample includes 4 654 150 stars, of which only 130 114 (3%) are from six-parameter solutions. The Medium and High Quality samples, Geo_MQ and Geo_HQ, list 3 469 143 and 1 894 809 sources, respectively. These Geo samples appear to be comparable to their corresponding astrometric control samples in terms of number of objects, spatial extension, and astrometric quality.

In order to compare trigonometric and Geo distances, we consider the fractional residuals, (r_Geo − r_ϖ)/r_ϖ, for the common objects. Median residuals over the whole distance range are 0.04% for both the ACS and ACS_MQ, while for the high-quality sample, ACS_HQ, the median residual with respect to Geo_HQ is 0.07%. The median absolute deviation is 0.2% for the three data sets.

The Galactic distribution of the mean fractional residuals is shown in Figure A.1 for the ACS and ACS_MQ on the left- and right-hand panel, respectively. The predominance of the colour green, with only traces of blue, indicates general good agreement (≲ 0.2%) between the two distance estimates away from the Galactic plane. The yellow-orange-red colour of the in-plane regions signifies that Bayesian and Geo distances are ≲ 4% larger than trigonometric parallax distances.

Figure A.2 shows the fractional residuals as a function of r_ϖ. The median is shown by the red line, while the yellow curves trace the 16% and 84% quantiles. The left panel shows the full ACS extending up to 10 kpc. At larger distances (> 6.5 kpc), trigonometric distances are larger than the Geo ones; conversely, below 6.5 kpc, r_ϖ is slightly smaller, that is by ∼1%. This systematic residual is produced by stars close to the Galactic plane (cf. Fig. A.1). Although better agreement is found between the cross-matched ACS_MQ and the Geo_MQ (right panels of Figs. A.1-A.2), similar systematic effects, such as a function of distance and Galactic latitude, are still present. These relatively small differences are consistent with the intrinsic biases expected in the trigonometric distances as a function of the parallax relative error (see e.g. Kovalevsky 1998), and the interplay with analogous effects in the Geo distances.

As for possible consequences of those differences, we recall that the chemo-kinematical distributions discussed in Sect. 5 were scrutinised using both r_ϖ and r_Geo distances and were found to be statistically indistinguishable. Therefore, the present comparison allows us to conclude that the above-mentioned biases do not affect the different analysis of this paper nor its conclusions.

Appendix B: Table queries of data samples

The following ADQL queries to the Astrophysical_parameters table can be used to retrieve the data samples employed in this work. We first present the definition of the different working samples and then provide the ADQL queries at the end of this Appendix.

General sample (5,527,090 stars, Sects. 3 and 6):

From the complete \gspspec\ sample, we have removed:
(i) very cool dwarfs defined by (\T < 3800K \& \g>3.5)
(ii) the {\it KM-extended flag}-stars defined by
(\T<4150K \& 2.4 < \g < 3.6).

From this query, objects without r_Geo distances were also removed. Medium quality sample (4,140,759 stars, Sects. 3, 5, and 6):

General sample & Teff>3500 & 0<log(g)<5 &
Teff_unc < 375 & log(g)_unc < 0.5 & [M/H]_unc < 0.25 &
vbroadT < 2 & vbroadG < 2 & vbroadM < 2 &
vradT < 2 &vradG < 2 &vradM < 2 &
fluxNoise < 4 &extrapol < 3 &KMtypestars < 2

Gradient analysis sample (2,762,809 stars, Sect. 4):

vbroadT < 2 &vbroadG < 2 & vbroadM < 2 &
vradT < 2 & vradG < 2 & vradM<2 &
fluxNoise<2 & extrapol=0 & KMtypestars=0 &
NegFlux=0 & nanFlux=0 & emission=0 &
nullFluxErr=0 & RUWE<=1.4 & Teff>4000

High Quality sample (2,218,573 stars, Sect. 5):

Teff>3500 & 0<log(g)<5 AND
vbroadT=0 & vbroadG=0 & vbroadM=0 AND
vradT=0 & vradG=0 & vradM=0 AND
fluxNoise=0 & extrapol<3 & KMtypestars<2 AND
NegFlux=0 & nanFlux=0 & emission=0 & nullFluxErr=0

Individual abundances samples (Sects. 3 and 7):

High Quality & Teff<7000 AND
[X/Fe]_unc<0.2 & XUpLim<2 & XUncer=0 AND
IF [Ca/Fe]:
cafe_gspspec_nlines>=3 AND
cafe_gspspec_linescatter<0.1 AND
(logg_gspspec < 3.5) AND
((teff_gspspec != 4250) | (logg_gspspec != 1.5))
IF [N/Fe]:
nfe_gspspec_nlines>=2 AND
nfe_gspspec_linescatter<0.1 AND
[N/Fe]_unc<0.15
IF [Cr/Fe]: teff_gspspec > 4100 & logg_gspspec < 2.5
IF [Ni/Fe]: nife_gspspec_unc < 0.1 & logg_gspspec < 3
IF [Ce/Fe]: CeUpLim<=2 AND
logg_gspspec >= 1.21 & logg_gspspec <= 2.46 AND
((teff_gspspec != 4250 | logg_gspspec != 1.5))

ADQL queries for the working samples from the Gaia Archive

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE ((teff_gspspec>=3800) OR (logg_gspspec>=3.5))
    AND ((teff_gspspec>=4150) OR (logg_gspspec>=3.6)
    OR (logg_gspspec<=2.4))

Listing 1. ADQL query for the General sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500) AND (logg_gspspec>0) AND
    (logg_gspspec<5) AND
    ((teff_gspspec_upper-teff_gspspec_lower)<750)
    AND ((logg_gspspec_upper-logg_gspspec_lower)<1.)
    AND ((mh_gspspec_upper-mh_gspspec_lower)<.5) AND
    (teff_gspspec>=3800 OR logg_gspspec<=3.5) AND
    (teff_gspspec>=4150 OR logg_gspspec<=2.4 OR
    logg_gspspec>=3.6 ) AND ((flags_gspspec LIKE
    "____________0%") OR (flags_gspspec LIKE
    "____________1%")) AND ((flags_gspspec LIKE
    "0%") OR (flags_gspspec LIKE "1%")) AND
    ((flags_gspspec LIKE "_0%") OR (flags_gspspec
    LIKE "_1%")) AND ((flags_gspspec LIKE "__0%") OR
    (flags_gspspec LIKE "__1%")) AND ((flags_gspspec
    LIKE "___0%") OR (flags_gspspec LIKE "___1%"))
    AND ((flags_gspspec LIKE "____0%") OR
    (flags_gspspec LIKE "____1%")) AND
    ((flags_gspspec LIKE "_____0%") OR
    (flags_gspspec LIKE "_____1%")) AND
    ((flags_gspspec LIKE "______0%") OR
    (flags_gspspec LIKE "______1%") OR
    (flags_gspspec LIKE "______2%") OR
    (flags_gspspec LIKE "______3%")) AND
    ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%"))

Listing 2. ADQL query for the Medium quality sample.

SELECT gaia.source_id
FROM gaiadr3.gaia_source as gaia INNER JOIN (
SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE ((flags_gspspec LIKE "0%") OR (flags_gspspec
    LIKE "1%")) AND ((flags_gspspec LIKE "_0%") OR
    (flags_gspspec LIKE "_1%")) AND ((flags_gspspec
    LIKE "__0%") OR (flags_gspspec LIKE "__1%")) AND
    ((flags_gspspec LIKE "___0%") OR (flags_gspspec
    LIKE "___1%")) AND ((flags_gspspec LIKE
    "____0%") OR (flags_gspspec LIKE "____1%")) AND
    ((flags_gspspec LIKE "_____0%") OR
    (flags_gspspec LIKE "_____1%")) AND
    ((flags_gspspec LIKE "______0%") OR
    (flags_gspspec LIKE "______1%")) AND
    (flags_gspspec LIKE "_______0%") AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    (flags_gspspec LIKE "____________0%") AND
    (teff_gspspec>4000) ) as aux USING (source_id)
    WHERE(gaia.ruwe < =1.4)

Listing 3. ADQL query for the Gradient analysis sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500) AND (logg_gspspec>0 AND
    logg_gspspec<5) AND (flags_gspspec LIKE "0%")
    AND (flags_gspspec LIKE "_0%") AND
    (flags_gspspec LIKE "__0%") AND (flags_gspspec
    LIKE "___0%") AND (flags_gspspec LIKE "____0%")
    AND (flags_gspspec LIKE "_____0%") AND
    (flags_gspspec LIKE "______0%") AND
    ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%")

Listing 4. ADQL query for the High quality sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500 AND teff_gspspec<7000) AND
    (logg_gspspec>0 AND logg_gspspec<5) AND
    (flags_gspspec LIKE "0%") AND (flags_gspspec
    LIKE "_0%") AND (flags_gspspec LIKE "__0%") AND
    (flags_gspspec LIKE "___0%") AND (flags_gspspec
    LIKE "____0%") AND (flags_gspspec LIKE
    "_____0%") AND (flags_gspspec LIKE "______0%")
    AND ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    ((cafe_gspspec_upper-cafe_gspspec_lower)<0.4)
    AND ((flags_gspspec LIKE
    "_____________________0%") OR (flags_gspspec
    LIKE "_____________________1%")) AND
    (flags_gspspec LIKE "______________________0%")
    AND (cafe_gspspec_nlines>=3) AND
    (cafe_gspspec_linescatter<0.1) AND
    (logg_gspspec<3.5) AND ((teff_gspspec!=4250) OR
    (logg_gspspec!=1.5))

Listing 5. ADQL query for the Calcium abundance sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500 AND teff_gspspec<7000) AND
    (logg_gspspec>0 AND logg_gspspec<5) AND
    (flags_gspspec LIKE "0%") AND (flags_gspspec
    LIKE "_0%") AND (flags_gspspec LIKE "__0%") AND
    (flags_gspspec LIKE "___0%") AND (flags_gspspec
    LIKE "____0%") AND (flags_gspspec LIKE
    "_____0%") AND (flags_gspspec LIKE "______0%")
    AND ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    ((nfe_gspspec_upper-nfe_gspspec_lower)<0.30) AND
    ((flags_gspspec LIKE "_____________0%") OR
    (flags_gspspec LIKE "_____________1%")) AND
    (flags_gspspec LIKE "______________0%") AND
    (nfe_gspspec_nlines>=2) AND
    (nfe_gspspec_linescatter<0.1)

Listing 6. ADQL query for the Nitrogen abundance sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500 AND teff_gspspec<7000) AND
    (logg_gspspec>0 AND logg_gspspec<5) AND
    (flags_gspspec LIKE "0%") AND (flags_gspspec
    LIKE "_0%") AND (flags_gspspec LIKE "__0%") AND
    (flags_gspspec LIKE "___0%") AND (flags_gspspec
    LIKE "____0%") AND (flags_gspspec LIKE
    "_____0%") AND (flags_gspspec LIKE "______0%")
    AND ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    ((crfe_gspspec_upper-crfe_gspspec_lower)<0.4)
    AND ((flags_gspspec LIKE
    "_________________________0%") OR (flags_gspspec
    LIKE "_________________________1%")) AND
    (flags_gspspec LIKE
    "_________________________0%") AND
    ((teff_gspspec>4100) AND (logg_gspspec<2.5))
    

Listing 7. ADQL query for the Chromium abundance sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500 AND teff_gspspec<7000) AND
    (logg_gspspec>0 AND logg_gspspec<5) AND
    (flags_gspspec LIKE "0%") AND (flags_gspspec
    LIKE "_0%") AND (flags_gspspec LIKE "__0%") AND
    (flags_gspspec LIKE "___0%") AND (flags_gspspec
    LIKE "____0%") AND (flags_gspspec LIKE
    "_____0%") AND (flags_gspspec LIKE "______0%")
    AND ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    ((nife_gspspec_upper-nife_gspspec_lower)<0.2)
    AND ((flags_gspspec LIKE "%0_________") OR
    (flags_gspspec LIKE "%1_________")) AND
    (flags_gspspec LIKE "%0________") AND
    (logg_gspspec<3)

Listing 8. ADQL query for the Nickel abundance sample.

SELECT source_id
FROM gaiadr3.astrophysical_parameters
WHERE (teff_gspspec>3500 AND teff_gspspec<7000) AND
    (logg_gspspec>0 AND logg_gspspec<5) AND
    (flags_gspspec LIKE "0%") AND (flags_gspspec
    LIKE "_0%") AND (flags_gspspec LIKE "__0%") AND
    (flags_gspspec LIKE "___0%") AND (flags_gspspec
    LIKE "____0%") AND (flags_gspspec LIKE
    "_____0%") AND (flags_gspspec LIKE "______0%")
    AND ((flags_gspspec LIKE "_______0%") OR
    (flags_gspspec LIKE "_______1%") OR
    (flags_gspspec LIKE "_______2%") ) AND
    ((flags_gspspec LIKE "____________0%") OR
    (flags_gspspec LIKE "____________1%") ) AND
    (flags_gspspec LIKE "________0%") AND
    (flags_gspspec LIKE "_________0%") AND
    (flags_gspspec LIKE "__________0%") AND
    (flags_gspspec LIKE "___________0%") AND
    ((cefe_gspspec_upper-cefe_gspspec_lower)<0.4)
    AND ((flags_gspspec LIKE "%0_____") OR
    (flags_gspspec LIKE "%1_____") OR (flags_gspspec
    LIKE "%2_____")) AND (flags_gspspec LIKE
    "%0____") AND ((logg_gspspec>=1.21) AND
    (logg_gspspec<=2.46)) AND ((teff_gspspec!=4250)
    OR (logg_gspspec!=1.5))

Listing 9. ADQL query for the Cerium abundance sample.

Appendix C: The Milky Way in [α/Fe] for different distance intervals

Here, we provide a complementary figure to Fig. 4 in Sect. 3, showing the [α/Fe] content of the Milky Way, cut in various distance intervals.

Appendix D: Impact of data calibrations on [α/Fe] versus [M/H] trends

Discussing the different abundance patterns observed in the investigated Galactic stellar populations as presented in Sect.3.2 requires knowledge of the used abundance calibrations and their relative impact on the studied samples. This is presented in Fig. D.1 for the observed distributions in the [α/Fe] versus [M/H] plane for (i) the raw uncalibrated data (upper panel), (ii) an [α/Fe] calibration based on log(g) (left middle panel) or T_eff(left lower panel), and (iii) these two [α/Fe] calibration flavours combined with a [M/H] calibration based on log(g) (middle and lower right panels). The calibration polynomial are taken from Tables 3 and 4 in Recio-Blanco et al. (2023). In addition, the figure allows to compare the three samples defined in 3.2 (HotTO in orange, RGB in grey and Massive in turquoise) and how their relative positions in the [α/Fe] versus [M/H] plane depend on the calibrations. To help this comparison, the locus of the thin-disc sequence in the RGB sample is highlighted, distinguishing an outer (green points) and an inner (pink) thin disc, as selected through the stellar orbit apocentres (< 7.5 kpc and > 12.5 kpc for the inner and outer disc respectively) and Zmax values (< 0.5 kpc).

Figure D.1 shows that the Massive sample is, in all the combinations of calibration and no calibration flavours, impoverished in [M/H] with respect to the thin disc sequence. In addition, the maximum metallicity reached by the Massive sample is also lower than that of the RGB and HotTO samples, independently of the chosen calibration procedure.

Appendix E: Disc flaring

Figure E.1 presents a kinematical evidence of the disc flaring discussed in Sect. 5.1.

Appendix F: Chemo-kinematical asymmetries in the outer disc

Figure F.1 shows the V_Z versus V_ϕ distributions from R = 10 kpc to 14.5 kpc and |Z|< 1 kpc for the Medium Quality stars with good astrometric parameters, i.e. not duplicated sources with RUWE < 1.4. These plots are consistent with the asymmetric velocity distributions found in Gaia Collaboration (2021c) towards the Galactic anticenter.

Appendix G: Kiel diagrams of individual abundance samples in the solar cylinder

For a better understanding of the abundance plots in the solar cylinder presented in Section 7, Figure G.1 shows the Kiel diagrams of the analysed individual abundance samples. The differences in the underlying selection function are clearly illustrated with this figure. They are imposed by the dectectability of the spectral lines of each element in the RVS spectra of the stars, which depends on the stellar parameters, the element abundance and the S/N.

Appendix H: Full acknowledgements

The Gaia mission and data processing have financially been supported by, in alphabetical order by country:

– the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique of Bouzareah Observatory;

– the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha Firnberg Programme through grants T359, P20046, and P23737;

– the BELgian federal Science Policy Office (BELSPO) through various PROgramme de Développement d’Expériences scientifiques (PRODEX) grants, the Research Foundation Flanders (Fonds Wetenschappelijk Onderzoek) through grant VS.091.16N, the Fonds de la Recherche Scientifique (FNRS), and the Research Council of Katholieke Universiteit (KU) Leuven through grant C16/18/005 (Pushing AsteRoseismology to the next level with TESS, GaiA, and the Sloan DIgital Sky SurvEy – PARADISE);

– the Brazil-France exchange programmes Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível Superior (CAPES) - Comité Français d’Evaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB);

– the Chilean Agencia Nacional de Investigación y Desarrollo (ANID) through Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) Regular Project 1210992 (L. Chemin);

– the National Natural Science Foundation of China (NSFC) through grants 11573054, 11703065, and 12173069, the China Scholarship Council through grant 201806040200, and the Natural Science Foundation of Shanghai through grant 21ZR1474100;

– the Tenure Track Pilot Programme of the Croatian Science Foundation and the École Polytechnique Fédérale de Lausanne and the project TTP-2018-07-1171 ‘Mining the Variable Sky’, with the funds of the Croatian-Swiss Research Programme;

– the Czech-Republic Ministry of Education, Youth, and Sports through grant LG 15010 and INTER-EXCELLENCE grant LTAUSA18093, and the Czech Space Office through ESA PECS contract 98058;

– the Danish Ministry of Science;

– the Estonian Ministry of Education and Research through grant IUT40-1;

– the European Commission’s Sixth Framework Programme through the European Leadership in Space Astrometry (ELSA) Marie Curie Research Training Network (MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European Commission’s Seventh Framework Programme through grant FP7-606740 (FP7-SPACE-2013-1) for the Gaia European Network for Improved data User Services (GENIUS) and through grant 264895 for the Gaia Research for European Astronomy Training (GREAT-ITN) network;

– the European Cooperation in Science and Technology (COST) through COST Action CA18104 ‘Revealing the Milky Way with Gaia (MW-Gaia)’;

– the European Research Council (ERC) through grants 320360, 647208, and 834148 and through the European Union’s Horizon 2020 research and innovation and excellent science programmes through Marie Skłodowska-Curie grant 745617 (Our Galaxy at full HD – Gal-HD) and 895174 (The build-up and fate of self-gravitating systems in the Universe) as well as grants 687378 (Small Bodies: Near and Far), 682115 (Using the Magellanic Clouds to Understand the Interaction of Galaxies), 695099 (A sub-percent distance scale from binaries and Cepheids – CepBin), 716155 (Structured ACCREtion Disks – SACCRED), 951549 (Sub-percent calibration of the extragalactic distance scale in the era of big surveys – UniverScale), and 101004214 (Innovative Scientific Data Exploration and Exploitation Applications for Space Sciences – EXPLORE);

– the European Science Foundation (ESF), in the framework of the Gaia Research for European Astronomy Training Research Network Programme (GREAT-ESF);

– the European Space Agency (ESA) in the framework of the Gaia project, through the Plan for European Cooperating States (PECS) programme through contracts C98090 and 4000106398/12/NL/KML for Hungary, through contract 4000115263/15/NL/IB for Germany, and through PROgramme de Développement d’Expériences scientifiques (PRODEX) grant 4000127986 for Slovenia;

– the Academy of Finland through grants 299543, 307157, 325805, 328654, 336546, and 345115 and the Magnus Ehrnrooth Foundation;

– the French Centre National d’Études Spatiales (CNES), the Agence Nationale de la Recherche (ANR) through grant ANR-10-IDEX-0001-02 for the ‘Investissements d’avenir’ programme, through grant ANR-15-CE31-0007 for project ‘Modelling the Milky Way in the Gaia era’ (MOD4Gaia), through grant ANR-14-CE33-0014-01 for project ‘The Milky Way disc formation in the Gaia era’ (ARCHEOGAL), through grant ANR-15-CE31-0012-01 for project ‘Unlocking the potential of Cepheids as primary distance calibrators’ (UnlockCepheids), through grant ANR-19-CE31-0017 for project ‘Secular evolution of galaxies’ (SEGAL), and through grant ANR-18-CE31-0006 for project ‘Galactic Dark Matter’ (GaDaMa), the Centre National de la Recherche Scientifique (CNRS) and its SNO Gaia of the Institut des Sciences de l’Univers (INSU), its Programmes Nationaux: Cosmologie et Galaxies (PNCG), Gravitation Références Astronomie Métrologie (PNGRAM), Planétologie (PNP), Physique et Chimie du Milieu Interstellaire (PCMI), and Physique Stellaire (PNPS), the ‘Action Fédératrice Gaia ’ of the Observatoire de Paris, the Région de Franche-Comté, the Institut National Polytechnique (INP) and the Institut National de Physique nucléaire et de Physique des Particules (IN2P3) co-funded by CNES;

– the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901, 50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, 50QG1404, 50QG1904, 50QG2101, 50QG2102, and 50QG2202, and the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität Dresden for generous allocations of computer time;

– the Hungarian Academy of Sciences through the Lendület Programme grants LP2014-17 and LP2018-7 and the Hungarian National Research, Development, and Innovation Office (NKFIH) through grant KKP-137523 (‘SeismoLab’);

– the Science Foundation Ireland (SFI) through a Royal Society - SFI University Research Fellowship (M. Fraser);

– the Israel Ministry of Science and Technology through grant 3-18143 and the Tel Aviv University Center for Artificial Intelligence and Data Science (TAD) through a grant;

– the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0, 2014-025-R.0, 2014-025-R.1.2015, and 2018-24-HH.0 to the Italian Istituto Nazionale di Astrofisica (INAF), contract 2014-049-R.0/1/2 to INAF for the Space Science Data Centre (SSDC, formerly known as the ASI Science Data Center, ASDC), contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), INAF, and the Italian Ministry of Education, University, and Research (Ministero dell’Istruzione, dell’Università e della Ricerca) through the Premiale project ‘MIning The Cosmos Big Data and Innovative Italian Technology for Frontier Astrophysics and Cosmology’ (MITiC);

– the Netherlands Organisation for Scientific Research (NWO) through grant NWO-M-614.061.414, through a VICI grant (A. Helmi), and through a Spinoza prize (A. Helmi), and the Netherlands Research School for Astronomy (NOVA);

– the Polish National Science Centre through HARMONIA grant 2018/30/M/ST9/00311 and DAINA grant 2017/27/L/ST9/03221 and the Ministry of Science and Higher Education (MNiSW) through grant DIR/WK/2018/12;

– the Portuguese Fundação para a Ciência e a Tecnologia (FCT) through national funds, grants SFRH/BD/128840/2017 and PTDC/FIS-AST/30389/2017, and work contract DL 57/2016/CP1364/CT0006, the Fundo Europeu de Desenvolvimento Regional (FEDER) through grant POCI-01-0145-FEDER-030389 and its Programa Operacional Competitividade e Internacionalização (COMPETE2020) through grants UIDB/04434/2020 and UIDP/04434/2020, and the Strategic Programme UIDB/00099/2020 for the Centro de Astrofísica e Gravitação (CENTRA);

– the Slovenian Research Agency through grant P1-0188;

– the Spanish Ministry of Economy (MINECO/FEDER, UE), the Spanish Ministry of Science and Innovation (MICIN), the Spanish Ministry of Education, Culture, and Sports, and the Spanish Government through grants BES-2016-078499, BES-2017-083126, BES-C-2017-0085, ESP2016-80079-C2-1-R, ESP2016-80079-C2-2-R, FPU16/03827, PDC2021-121059-C22, RTI2018-095076-B-C22, and TIN2015-65316-P (‘Computación de Altas Prestaciones VII’), the Juan de la Cierva Incorporación Programme (FJCI-2015-2671 and IJC2019-04862-I for F. Anders), the Severo Ochoa Centre of Excellence Programme (SEV2015-0493), and MICIN/AEI/10.13039/501100011033 (and the European Union through European Regional Development Fund ‘A way of making Europe’) through grant RTI2018-095076-B-C21, the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M, the University of Barcelona’s official doctoral programme for the development of an R+D+i project through an Ajuts de Personal Investigador en Formació (APIF) grant, the Spanish Virtual Observatory through project AyA2017-84089, the Galician Regional Government, Xunta de Galicia, through grants ED431B-2021/36, ED481A-2019/155, and ED481A-2021/296, the Centro de Investigación en Tecnologías de la Información y las Comunicaciones (CITIC), funded by the Xunta de Galicia and the European Union (European Regional Development Fund – Galicia 2014-2020 Programme), through grant ED431G-2019/01, the Red Española de Supercomputación (RES) computer resources at MareNostrum, the Barcelona Supercomputing Centre - Centro Nacional de Supercomputación (BSC-CNS) through activities AECT-2017-2-0002, AECT-2017-3-0006, AECT-2018-1-0017, AECT-2018-2-0013, AECT-2018-3-0011, AECT-2019-1-0010, AECT-2019-2-0014, AECT-2019-3-0003, AECT-2020-1-0004, and DATA-2020-1-0010, the Departament d’Innovació, Universitats i Empresa de la Generalitat de Catalunya through grant 2014-SGR-1051 for project ‘Models de Programació i Entorns d’Execució Parallels’ (MPEXPAR), and Ramon y Cajal Fellowship RYC2018-025968-I funded by MICIN/AEI/10.13039/501100011033 and the European Science Foundation (‘Investing in your future’);

– the Swedish National Space Agency (SNSA/Rymdstyrelsen);

– the Swiss State Secretariat for Education, Research, and Innovation through the Swiss Activités Nationales Complémentaires and the Swiss National Science Foundation through an Eccellenza Professorial Fellowship (award PCEFP2_194638 for R. Anderson);

– the United Kingdom Particle Physics and Astronomy Research Council (PPARC), the United Kingdom Science and Technology Facilities Council (STFC), and the United Kingdom Space Agency (UKSA) through the following grants to the University of Bristol, the University of Cambridge, the University of Edinburgh, the University of Leicester, the Mullard Space Sciences Laboratory of University College London, and the United Kingdom Rutherford Appleton Laboratory (RAL): PP/D006511/1, PP/D006546/1, PP/D006570/1, ST/I000852/1, ST/J005045/1, ST/K00056X/1, ST/K000209/1, ST/K000756/1, ST/L006561/1, ST/N000595/1, ST/N000641/1, ST/N000978/1, ST/N001117/1, ST/S000089/1, ST/S000976/1, ST/S000984/1, ST/S001123/1, ST/S001948/1, ST/S001980/1, ST/S002103/1, ST/V000969/1, ST/W002469/1, ST/W002493/1, ST/W002671/1, ST/W002809/1, and EP/V520342/1.

The GBOT programme uses observations collected at (i) the European Organisation for Astronomical Research in the Southern Hemisphere (ESO) with the VLT Survey Telescope (VST), under ESO programmes 092.B-0165, 093.B-0236, 094.B-0181, 095.B-0046, 096.B-0162, 097.B-0304, 098.B-0030, 099.B-0034, 0100.B-0131, 0101.B-0156, 0102.B-0174, and 0103.B-0165; and (ii) the Liverpool Telescope, which is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias with financial support from the United Kingdom Science and Technology Facilities Council, and (iii) telescopes of the Las Cumbres Observatory Global Telescope Network.

All Tables

All Figures