Gaia Data Release 3: Chemical cartography of the Milky Way

Gaia DR3 opens a new era of all-sky spectral analysis of stellar populations thanks to the nearly 5.6 million stars observed by the RVS and parametrised by the GSP-spec module. The all-sky Gaia chemical cartography allows a powerful and precise chemo-dynamical view of the Milky Way with unprecedented spatial coverage and statistical robustness. First, it reveals the strong vertical symmetry of the Galaxy and the flared structure of the disc. Second, the observed kinematic disturbances of the disc -- seen as phase space correlations -- and kinematic or orbital substructures are associated with chemical patterns that favour stars with enhanced metallicities and lower [alpha/Fe] abundance ratios compared to the median values in the radial distributions. This is detected both for young objects that trace the spiral arms and older populations. Several alpha, iron-peak elements and at least one heavy element trace the thin and thick disc properties in the solar cylinder. Third, young disc stars show a recent chemical impoverishment in several elements. Fourth, the largest chemo-dynamical sample of open clusters analysed so far shows a steepening of the radial metallicity gradient with age, which is also observed in the young field population. Finally, the Gaia chemical data have the required coverage and precision to unveil galaxy accretion debris and heated disc stars on halo orbits through their [alpha/Fe] ratio, and to allow the study of the chemo-dynamical properties of globular clusters. Gaia DR3 chemo-dynamical diagnostics open new horizons before the era of ground-based wide-field spectroscopic surveys. They unveil a complex Milky Way that is the outcome of an eventful evolution, shaping it to the present day (abridged).


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, Brown, A. G. A. et al. 2021).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. 2018a;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. 2019b;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 Article number, page 2 of 49 Gaia Collaboration et al.: Gaia Data Release 3: Chemical cartography of the Milky Way responsible for the rapidly precessing disc warp (Poggio et al. 2020(Poggio et al. , 2021b)).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 1925b,a), 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 hydrogen1 ; 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 conserved2 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. 2021), 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, Vallenari et al. 2022) 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. 2022) and parametrised by the General Stellar Parametriser -spectroscopy (GSP-Spec; Recio-Blanco et al. 2022).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.63 , 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 diagnostics4 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. (2022).
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 chemokinematical correlations, while in Sect.6 we explore the relation between the orbital parameters and stellar chemistry.Sub-sequently, Sect.7 is dedicated to chemo-dynamical relations in solar cylinder populations using individual element abundances, and in Section 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 chemodynamical 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.

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 astrophysi-cal_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 (mh_gspphot field).Although GSP-Phot metallicities are suitable for different scientific purposes, their application to large-scale Galactic chemodynamical 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) 5 , global metallicity [M/H]6 , and the global abundance of α-elements7 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 Ma-tisseGauguin GSP-Spec workflow and is described in detail in the GSP-Spec processing paper (Recio-Blanco et al. 2022).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 teff_gspspec field; log(g) comes from the logg_gspspec field; [M/H] is taken from mh_gspspec; and [α/Fe] corresponds to alphafe_gspspec with a calibration8 imposed that requires a zero average value for [α/Fe] in the solar neighbourhood for any gravity (see Recio-Blanco et al. 2022, Table 16).
In addition, we make use of the GSP-Spec quality flags reported in the flags_gspspec string chain (which consists of 41 characters) defined in Recio-Blanco et al. (2022).For example, we make use of the first three characters in this chain (that is, flags_gspspec[0], flags_gspspec[1], and flags_gspspec[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. 2022).
Finally, the uncertainty on any derived parameter (Param_unc) or abundance ([X/Fe]_unc) is defined as half of the difference between its upper and lower confidence levels (e.g.Teff_unc = [teff_gspspec_upper − teff_gspspec_lower]/2).

Radial velocities
The complete GSP-Spec sample contains 5 594 205 stars in total (based on the query flags_gspspec is not null).DR3 provides radial velocities (radial_velocity, V Rad ) for 33 834 834 stars (Katz et al. 2022) 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. 2022).This was found to be very similar to rv_expected_sig_to_noise9 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. 2022).The sensitivity of GSP-Spec parametrisation to spectra that are not perfectly corrected for their radial velocity shift is flagged in the GSP-Spec flags_gspspec string.In particular, characters 3 to 5 (flags_gspspec [3], flags_gspspec[4] and flags_gspspec [5]), called vradT, vradG, and vradM respectively, report on the degree of parameter biases from V Rad uncertainties (see Recio-Blanco et al. 2022, for more details on these flags).

Astrometric data and distances
High-precision astrometric parameters (α, δ, µ α * , µ δ , ) from Gaia EDR3 (Gaia Collaboration, Brown, A. G. A. et al. 2021) 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 crossmatching 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.

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 et al. 2021;Bennett & Bovy 2019) and Galactocentric cylindrical velocities 9.5, 250.7, 8.56) km s −1 (Reid & Brunthaler 2020;Gravity Collaboration et al. 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 −1 (Gravity Collaboration et al. 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-ofsight velocities, and a broken Gaussian distribution for distance (based on the lower and upper confidence limits of Bailer-Jones et al. 2021).

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.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 f luxNoise = 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.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.

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, Prusti, T. et al. 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, 2, and 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  Figures 4 and C.1 in the Appendix 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 Article number, page 6 of 49 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, Luri, X. et al. 2021), 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 lowmetallicity thick disc and halo stars, while for R>12 kpc, flared disc stars of higher metallicity and lower [α/Fe] values dominate.

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, Article number, page 7 of 49 in this section we explore its impact so as to help DR3 users and to support the chemo-dynamical analysis of this paper.

Stellar parameter distribution across the Galaxy
To dissect the covered spatial volume, we divided the (R, Z) plane into bins of 2x2 kpc 2 .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 de-creasing trend of [α/Fe] with metallicity can be seen across the Galaxy.

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 Figures 9,10,and 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 Article number, page 8 of 49 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 (c.f.Sect.2.5), the number of reported abundances is huge (from ∼ 10 4 to a couple of ∼ 10 6 , 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. ( 2022) 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 agree-ment with the spatial maps derived by Poggio et al. (2021a) 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, Schultheis et al. 2022).

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 Article number, page 10 of 49

[Ca/Fe] distributions
We now look at the [Ca/Fe] abundances presented in the upper right panels of Figures 9,10,and 11 (individual abundances samples).Concerning HotTO stars, Figure 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.0dex at [M/H]∼+0.25dex.
Finally, Fig. 11 illustrates the [Ca/Fe] abundance for about 7 400 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 Sect.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 Figures 9,10,and 11.For HotTO stars, Figure 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.06dex).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 4 000 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.
Article number, page 12 of 49 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 7 000 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.

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. 2022), 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. 2021).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   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 −1 .This is slightly different from but still consistent with the metallicity gradients found using Cepheids (e.g.−0.060 ± 0.002 dex kpc  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. 2020).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 −1 ).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 7kpc, we find a value of +0.007 ± 0.  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 dominates10 (> 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]     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 −1 ) 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.

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.

Toomre diagram and Galactic populations
Although the mean disc rotation varies as a function of R and Z, here we adopt a constant value, V 0 = 238.5 km s −1 , 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. 2018b;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 Article number, page 15 of 49 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.

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, Antoja, T. et al. (2021).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, Antoja, T. et al. (2021) 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, Antoja, T. et al. (2021) 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 −1 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 −1 , is more significant.This is clear evidence of the bimodal distribution found by Gaia Collaboration, Antoja, T. et al. (2021) 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 Figure E.1 which shows the V Z versus V φ distribution for a set of radial bins from R = 10 kpc to 14 kpc.

Vertical velocity oscillations
We focus on the Z versus V Z distribution, where vertical wavelike motions have been detected in the form of a snail-like pat-tern 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 −1 , 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  On the right panel, that is 9.25 < R < 11.25 kpc, stars with V φ ∼ 230 km s −1 follow a spiral-like pattern with a tail in the region at V Z > 20 km s −1 .As expected from dynamical considerations (Binney & Tremaine 2008), and as previously reported in the literature using Gaia DR2 data (Laporte et al. 2019b), the Z-axis elongation increases with Galactic radius.Indeed, at large radii, dynamical timescales are longer due to a weaker selfgravity 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 phasespace 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 metalrich 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.

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, T. A. et al. 2018;Fragkoudi et al. 2019;Hunt et al. 2019;Monari et al. 2019b).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 Article number, page 18 of 49 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.  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 Article number, page 20 of 49  2021) 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 Fig. 21 and Fig. 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 Section 4).High Quality sample, with the additional conditions on the astrometric parameters defined at the beginning of this section.

Rotation-metallicity correlations
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 −1 ) with low [α/Fe].We notice that the azimuthal velocity of metalrich 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 −1 , 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 thickdisc sequences as in other chemical elements, such as Mg, which Article number, page 21 of 49 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 −1 .We find a gradient of −13.9 ± 0.1 km s −1 dex −1 , 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 −1 .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. 2011Kordopatis et al. , 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 −1 to reduce the contamination from thin-disc and halo stars, we measure a gradient of +33 ± 7 km s −1 dex −1 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 Sect. 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.

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 V 2 , respectively.

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 (counterrotating) systems.As seen, this sample is clearly dominated by stars in the solar neighbourhood (L Z = L , E ≈ −2.88V 2 ) 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.8V 2 E −2.0V 2 associated with the Gaia-Enceladus-Sausage (GES; Belokurov et al. 2018a;Myeong et al. 2018;Helmi et al. 2018;Feuillet et al.Article number, page 22 of 49 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 −1 , 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.84V 2 ) and (0.86L , −2.30V 2 ).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χ 2 < −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 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. 2020Feuillet et al. , 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.
Article number, page 23 of 49  (2019).The colour code is saturated in the high-density regions to emphasise the details of the selected areas.No solar neighbourhood area is selected in this figure because it would contain no stars due to the applied selection criteria.

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.

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.2dex to [M/H]<-0.2dex.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.

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. 2022).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 We use the Individual abundance subsample (as defined in Sect.2.5), to which we apply some additional filters: astometric_params_solved = 31 OR 95 (stars with respectively 6-or 5-parameter astrometric solutions), ruwe < 1.4, duplicated_source 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.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 dredgeup 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. 2003Bensby et al. , 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. 2021;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. 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.
Article number, page 26 of 49  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 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  tic 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. 2022).It therefore probes a narrower range of Z max with median (Z max ) being found at around a few hundred parsecs at all metallicities.
Article number, page 28 of 49 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 Article number, page 29 of 49 median azimuthal velocity and higher median eccentricity structure (green/blue cells in columns 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 −1 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.  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.

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  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.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 i th cluster.The x i normal distributions have standard deviations equal to ∆[M/H] 2 + 2 .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 11 .
Interestingly, the radial metallicity gradient α=-0.054±0.008dex kpc −1 derived from our sample of 503 open clusters is consistent with the value given by Casamiquela et al. (2019), that is −0.056±0.011dex kpc −1 (18 open clusters), and that of Spina et al. (2022), namely −0.064±0.007dex kpc −1 (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.02dex kpc −1 given by Jacobson et al. (2016) (Donor et al. 2020;Spina et al. 2021, e.g.).The results are plotted in the bottom panel of Fig. 34.The posteriors are listed in Table 3.
Article number, page 31 of 49 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. 2021Spina et al. , 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.002dex kpc −1 is obtained for the Massive sample, while it is of -0.055±0.007dex kpc −1 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.

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 Section 5. 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 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.

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.

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, 4, and 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 chemodynamical 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.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;  3 as a function of R. As in the middle panel, the solid line shows the LOWESS function and the shaded area shows its 90% confidence interval. 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. (2021a), 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, Drimmel R. & et al. (2022).In the present paper, the metallicity distribution of the population Article number, page 33 of 49 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.

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, Antoja, T. et al. (2021), 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. 2019b) 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. 2019a, for Hyades and Hercules) imply azimuthal metallicity variations due to the corresponding chemokinematical 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.

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 chemodynamical 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 (c.f 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 migra-Article number, page 34 of 49 tion 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.

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.249kpc) 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.

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 forma-tion.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 chemodynamical 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 largescale 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.

Fig. 1 .
Fig. 1.Milky Way stars in the GSP-Spec database (General sample).This HEALPix map in Galactic coordinates has a resolution of 0.46 • .The colour code corresponds to the median of the stellar distance from the Sun.This colour code enhances thin and thick disc populations in the Galactic plane, far from the solar neighbourhood.The Galactic bulge stars are also visible in the central region.Finally, the higher distance of the Magellanic Clouds and of several globular clusters reveal their presence in the regions far from the Galactic plane otherwise dominated by foreground stars of the solar neighbourhood.
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. 2022).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.

Fig. 2 .
Fig.2.Same as Fig.1but colour coded with the median of the stellar metallicity [M/H] in each pixel.In the Galactic plane, the higher metallicity values of the thin disc are visible.In the central Galactic regions, a more metal-poor mix of bulge and thick disc populations is present.Far from the disc plane, highmetallicity thin disc stars at low distances from the Sun and more distant metal-poor halo stars are present.

Fig. 3 .
Fig. 3. Same as Fig.2 but colour-coded with the median of [α/Fe] in each pixel.The large structures close to the ecliptic poles are artifacts caused by the Gaia scanning law (see text).Thin disc stars are visible in the plane thanks to their low [α/Fe] values.
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).

Fig. 4 .
Fig. 4. Same as Fig. 2 but for different distance intervals shown from the closest (top panel of left column) to the most distant (bottom panel of right column).The distance ranges adopted in each panel are indicated in their lower right corner together with the number of stars in their bottom left corner.The colour code corresponds to the median of the metallicity and is continuous between the first four subpanels showing rather metal-rich stars, and the two last panels dominated by metal-poor populations.

Fig. 5 .
Fig. 5. Galactic maps for the General sample stars colour coded according to stellar count, the median of [M/H], and the median of [α/Fe] (from left to right, respectively).The top row shows the maps in Cartesian coordinates (X,Y), with the background image being the Milky Way sketch from Churchwell et al. (2009).In the bottom row, R and Z are the distances from the Galactic centre and plane, respectively.

Fig. 6 .
Fig. 6.Kiel diagrams across the Milky Way for the Medium Quality sample, colour coded according to the median of the mean metallicity.The number of stars in each subpanel is indicated in its upper left corner.The global (R, Z) distribution of Fig. 5 (bottom row) is shown in the background with grey levels.Moving away from the Sun, the sample becomes dominated by intrinsically bright giants.

Fig. 8 .
Fig. 8. Kiel diagram of all the Medium Quality stars located within 1 kpc of the Galactic plane.The location of the HotTO, RGB, and Massive subsamples is illustrated.Their spatial distribution and chemical properties are detailed in Figs. 9, 10, and 11, respectively.

Fig. 9 .
Fig. 9. HotTO stars of the Medium Quality sample located within 1 kpc of the Galactic plane.Their XY spatial distribution is shown within ±4 kpc of the Sun, colour coded according to stellar count and the median of the mean metallicity (left and central panels, respectively).The solar position is indicated by a filled star coloured according to solar metallicity ([M/H]=0.0dex).The right panel shows the evolution of calcium (top) and sulfur (bottom) abundances with respect to iron abundances versus the mean metallicity.

Fig. 10 .
Fig. 10.Same as Fig. 9 but for RGB stars within 1 kpc of the Galactic plane and over a more extended XY domain (±6 kpc from the Sun).Nickel abundances are shown in the bottom-right panel.

Fig. 11 .
Fig. 11.Same as Fig. 9 but for Massive stars located within 400 pc of the Galactic plane, thus focusing attention on the thin disc.The spiral arm structure is visible in the left and central panels.Calcium and nitrogen abundance distributions are shown in the right top and bottom panels, respectively.For comparison, the distribution of the RGB stars in terms of calcium abundance is shown in grey in the top panel .

Fig. 12 .
Fig. 12. Cerium abundances for the coolest Massive and RGB stars of the Medium Quality sample located close to the Galactic plane and colour coded according to stellar density and the median of [Ca/Fe] (top and bottom panel, respectively).

GaiaFig. 13 .
Fig. 13.Radial gradients for metallicity (top) and [α/Fe] (bottom) for different distances from the Galactic midplane (in kpc; see the legend).The trends are computed as running medians in bins of 0.5 kpc, with a 40 percent overlap, provided that at least 50 stars are available to compute the median.The shaded areas represent the Poisson uncertainty on the trends.The colours, which are associated to Z distances from the plane, are similar if the |Z| range is the same.

Fig. 14 .
Fig. 14.Vertical gradients for metallicity (top) and [α/Fe] (bottom), for different Galactocentric radial ranges.The trends are computed as running medians in bins of 0.5 kpc, with a 40 percent overlap, provided that at least 50 stars are available to compute the median.Shaded areas represent the Poisson uncertainties on the plotted values.
Figure13shows 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. 2021).Figure14is 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 fromHaywood 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 Article number, page 13 of 49 −1 Genovali et al. 2014, and −0.045 ± 0.07 dex kpc −1 , Lemasle et al. 2018 ), open clusters (e.g.−0.076 ± 0.09 dex kpc −1 Spina et al. 2021, and −0.068 ± 0.001 dex kpc −1 , Donor et al. 2020), or thin disc (i.e.low α−sequence) field stars (e.g.∼ −0.053 to ∼ −0.068 dex kpc −1 Bergemann et al. 2014; Recio-Blanco et al. 2014; 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. (

Fig. 15 .
Fig.15.Same as Fig.13, selecting this time the stars as a function of Z max and plotting the guiding radius (defined as (R apo + R peri )/2) instead of the observed radius.

Figure 17
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 V 2 R + V 2 Z 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:

Fig. 17 .
Fig. 17.Galactic map of the Toomre diagrams colour coded according to [M/H].The circular dashed line delimits the regions where the thin-and thick-disc stars are dominant.The vertical dashed line, V φ = 0, separates prograde and retrograde rotating stars.Halo stars seem preferentially located at positive V φ in many panels.

Fig. 20 )
follows a curled, spiral-shaped distribution particularly pronounced up to |V Z | 40 km s −1 .In the left panel, that is, 5.25 < R < 7.25 kpc, we still clearly see a snail-like structure Article number, page 17 of 49

Fig. 19 .
Fig. 19.Velocity distribution V Z vs. V φ of the Medium quality sample for |Z| < 1 and 11 < R < 11.5 kpc (top), 11.5 < R < 12 kpc (middle), 12 < R < 12.5 kpc (bottom).The colour coding for stars above and below the plane represents actual counts (first and second columns) and median [M/H] (last two columns), respectively.The bimodality in the velocity distribution observed below the plane and revealed by Gaia Collaboration, Antoja, T. et al. (2021) is accompanied by a metallicity bimodality as seen in the rightmost panels.

Fig. 20 .Fig. 21 .
Fig. 20.Galactic map of the velocity distribution V Z vs. Z of the Medium Quality sample for 5.25 < R < 11.25 kpc.From left to right, the radial range is shown in three rings at R = 6.25, 8.25, and 10.25 kpc, each of 2 kpc in width.From top to bottom: Normalised density distribution, Z − V Z plane coloured as a function of median V R , V φ , and ∆[M/H].Article number, page 19 of 49 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 ∆N = (N − N)/ N, where N and N are the bivariate kernel density estimators adopted in Appendix B1 of Poggio et al. (2021a) for the local and mean density, respectively.The adopted bandwidths for the local density are 100 pc in R and 2 km s −1 in V φ , while the corresponding values for the mean density are 700 pc and 14 km s −1 .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. (

Figure 24
Figure24shows 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

Fig. 23 .
Fig. 23.Velocity distribution V φ vs. R for massive giants within |Z| < 1 kpc.Top panel: Density distribution of the sample.Middle panel: Distribution colour coded according to median metallicity.Bottom panel: Residual metallicity, ∆[M/H] (see text).Black lines shows the known ridges from Ramos et al. (2018).

Fig. 24 .
Fig. 24.Galactic map of the [α/Fe] vs. [M/H] distribution of the High Quality sample colour coded according to median V φ .

Fig. 25 .
Fig. 25.Density map in the (L Z , E) diagram for the GSP-Spec sample.Dashed ellipses represent the contours of regions discussed in Sect.6.The shaded area corresponds to pairs (L Z , E) with no physical meaning, in which the upper boundary (solid black line) implies perfect circular orbits.
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.3dex), with lower [α/Fe] values compared to the disc stars at the same metallicity (density plot in the background, corresponding to the white circle in Fig. 25 and Fig. 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.

Fig. 26 .
Fig. 26.Zoomed-in version of the area enclosed by grey rectangles in Fig. 25 imposing the selection criteria described in Koppelman et al. (2019).The colour code is saturated in the high-density regions to emphasise the details of the selected areas.No solar neighbourhood area is selected in this figure because it would contain no stars due to the applied selection criteria.

Fig. 27 .
Fig. 27.Distribution of median metallicities (left panel) and α−element enrichment with respect to iron (right panel) in the energy-angular momentum (E, L Z ) plane for the General sample stars without those sources with log(g)<0.5.The coloured ellipses illustrate the selected areas associated with Gaia-Sausage-Enceladus (black), Thamnos (magenta), the Helmi stream (red), Sequoia (orange), and the solar neighbourhood sample (white).
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.
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. (2022) 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. 2022, 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.Article number, page 25 of 49

Fig. 28 .
Fig. 28.[α/Fe] versus [M/H] diagram for all the stars contained in the regions shown in Fig. 25 that satisfy the Medium Quality selection criteria.The colour code represents the vertical component of the angular momentum L Z .The density plot in the background corresponds to the Medium Quality subsample of solar neighbourhood stars.
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. 2022).

Fig. 29 .
Fig. 29.Distribution of globular cluster stars in the E vs. 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 metallicity.The density plot in the background corresponds to the General sample.

Fig. 30 .
Fig. 30.Distribution of [α/Fe] abundances with respect to [M/H] for the subsample of GCs with more reliable abundances.Error bars are based on the abundance dispersion within the cluster.A metallicity offset of 0.1 dex with respect to the literature is observed and corrected.Globular cluster stars are generally in the low-S/N regime of the GSP-Spec sample and need more deblending corrections of their spectra due to the crowding.The colour code reflects the GC L Z median values in solar units.In the background, a density plot with the RGB stars sample distribution is presented.
Figure G.1).As a consequence, the tracers of different chemical elements probe different ranges of maximum distance to the Galac-Article number, page 27 of 49

Fig. 31 .
Fig. 31.Distribution of the General sample in the J R − L Z diagram colour coded according to density (upper panels), median maximum distance above the Galactic plane Z max (second panels), median metallicity (third panels), and median [α/Fe] (lower panels).The grey boxes in the left panels indicate the areas for which a zoomed-in version is shown in the right column.The contour lines of the zoomed density plot are included in the right column as a visual reference, while the tentative positions of the ridges are denoted by the dashed lines.Bins with less than ten stars are omitted in the zoomed density plot to enhance the gradient in the colour code.The position of the LSR at (J R , L Z )=(0, L ) is denoted by the star symbol.
(2013); Haywood et al. (2013); Allende Prieto et al. (2016); Kordopatis et al. (2017); Re Fiorentin et al. ( 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. (2021); 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 Section 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. (2021); Tarricq et al. ( 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 Article number, page 30 of 49 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.

Fig. 34 .
Fig. 34.Galactic radial metallicity and [α/Fe] gradients traced by open clusters.Top: Cluster [M/H] values as a function of their Galactocentric distance R. Clusters are marked by circles colour coded according to their age.The grey shaded areas represent the 68% and 95% confidence intervals of the linear models resulting from the Bayesian regression, while the black dashed line traces the most probable model.Bottom: Cluster [α/Fe] values as a function of their Galactocentric distance R. Symbol colours and shaded areas are the same as those in the top panel.

Fig. 35 .
Fig. 35.Age dependence of the Galactic metallicity gradient traced by open clusters in the [M/H]-R (blue dots) and [M/H]-r guid (red dots) diagrams.
Figure 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.Article number, page 32 of 49 Gaia Collaboration et al.: Gaia Data Release 3: Chemical cartography of the Milky Way

Fig. 36 .
Fig. 36.Kinematic features in the radial metallicity distribution as traced by open clusters.Top: Distribution of open clusters in the V φ -R diagram.The 332 open clusters with [M/H] estimation are shown as grey dots.The red area shows the kernel density estimate of that sample.The Galactic ridges identified byRamos et al. (2018) are plotted as in Fig.21.Middle: Locally weighted scatterplot smoothing (LOWESS) regression of cluster V φ as a function of R (solid line).The shaded area shows the 90% confidence interval of the LOWESS functions obtained by sampling within the uncertainties of the two variables.Bottom: Residuals between the cluster [M/H] values and the metallicities predicted by the radial metallicity gradient as in Table3as a function of R. As in the middle panel, the solid line shows the LOWESS function and the shaded area shows its 90% confidence interval.

Table 1 .
Radial metallicity and [α/Fe]gradients for R > 7 kpc and all of the stellar types, at different distances from the plane.The (a, b) columns refer to the linear coefficients of the fit and the (e a , e b ) ones to their associated errors.

Table 2 .
Radial metallicity and [α/Fe] gradients for R > 7 kpc, at different distances from the plane when only giants are considered.