Context. Recent evidence from spectroscopic surveys points towards the presence of a metalpoor young stellar population in the low  chemically thin disk. In this context the investigation of the spatial distribution and time evolution of precise unbiased abundances is fundamental to disentangle the scenarios of formation and evolution of the Galaxy. Aims. We study the evolution of abundance gradients in the Milky Way by taking advantage of a large sample of open star clusters which are among the best tracers for this purpose. In particular we used data from the last release of the GaiaESO survey. Methods. We performed a careful selection of open cluster member stars excluding those members that may be affected by biases in spectral analysis. We compared the cleaned open cluster sample with detailed chemical evolution models for the Milky Way using welltested stellar yields and prescription for radial migration. We tested different scenarios of Galaxy evolution to explain the data namely the twoinfall and the threeinfall frameworks which suggest the chemical thin disk is formed by one or two subsequent gas accretion episodes respectively. Results. With the performed selection in cluster member stars we still find a metallicity decrease between intermediateage 1 lt AgeGyr lt 3 and young Age lt 1 Gyr open clusters. This decrease cannot be explained in the context of the twoinfall scenario even by accounting for the effect of migration and yield prescriptions. The threeinfall framework with its late gas accretion in the last 3 Gyr is able to explain the low metallic content in young clusters. However we have invoked a milder metal dilution for this gas infall episode relative to previous findings. Conclusions. To explain the observed low metallic content in young clusters we propose that a late gas accretion episode triggering a metal dilution would have taken place extending the framework of the threeinfall model for the first time to the entire Galactic disk.
['stars: abundances', 'Galaxy: abundances', 'Galaxy: disk', 'Galaxy: evolution', 'open clusters and associations: general', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
Mapping radial abundance gradients with GaiaESO open clusters Evidence of recent gas accretion in the Milky Way disk
{'Evidence of recent gas accretion in the Milky Way disk': "M. Palla 1 , 2 , L. Magrini 3 , E. Spitoni 4 , F. Matteucci 4 , 5 , 6 , C. Viscasillas Vázquez 7 , M. Franchini 4 , M. Molero 4 , 8 , and S. Randich 3  \n- 1 Dipartimento di Fisica e Astronomia 'Augusto Righi', Alma Mater Studiorum, Università di Bologna, Via Gobetti 93 / 2, 40129 Bologna, Italy\n- 2 INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93 / 3, 40129 Bologna, Italy e-mail: marco.palla@inaf.it\n- 3 INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy\n- 4 INAF - Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy\n- 5 Dipartimento di Fisica, Sezione di Astronomia, Università degli studi di Trieste, Via G.B. Tiepolo 11, I-34143 Trieste, Italy\n- 6 INFN - Sezione di Trieste, via A. Valerio 2, I-34100, Trieste, Italy\n- 7\n- Institute of Theoretical Physics and Astronomy, Vilnius University, Sauletekio av. 3, 10257 Vilnius, Lithuania\n- 8 Institut für Kernphysik, Technische Universität Darmstadt, Schlossgartenstr. 2, Darmstadt 64289, Germany  \nReceived XXX; accepted XXX", 'ABSTRACT': 'Context. Recent evidences from spectroscopic surveys point towards the presence of a metal-poor, young stellar population in the lowα / chemical thin disk. In this context, the investigation of the spatial distribution and time evolution of precise, unbiased abundances is fundamental to disentangle the scenarios of formation and evolution of the Galaxy.  \nAims. We study the evolution of abundance gradients in the Milky Way by taking advantage of a large sample of open star clusters, which are among the best tracers for this purpose. In particular, we use data from the last release of the Gaia -ESO survey.  \nMethods. We perform careful selection of open cluster member stars excluding those members that may be a ff ected by biases in spectral analysis. The cleaned open clusters sample is compared with detailed chemical evolution models for the Milky Way, using well tested stellar yields and prescription for radial migration. Di ff erent scenarios of Galaxy evolution are tested to explain the data, i.e. the two-infall and the three-infall frameworks, suggesting that the chemical thin disk is formed by one or two subsequent gas accretion episodes, respectively.  \nResults. With the performed selection in cluster member stars, we still find a metallicity decrease between intermediate age (1 < Age / Gyr < 3) and young (Age < 1 Gyr) open clusters. This decrease cannot be explained in the context of the two-infall scenario, even by accounting for the e ff ect of migration and yield prescriptions. The three-infall framework, with its late gas accretion in the last 3 Gyr, can explain the low metallic content in young clusters. However, we invoke a milder metal dilution for this gas infall episode relative to previous findings.  \nConclusions. To explain the observed low metallic content in young clusters, we propose that a late gas accretion episode triggering a metal dilution should have taken place, extending the framework of the three-infall model for the first time to the entire Galactic disk.  \nKey words. Galaxy: disk -- Galaxy: abundances - Galaxy: evolution - stars: abundances - open clusters and associations: general', '1. Introduction': 'A fundamental constraint to study the formation and chemical evolution of the Galaxy are abundance gradients along the Galactic disk.  \nDi ff erent stellar and nebular Galactic tracers that correspond to di ff erent epochs in the evolution of our Galaxy have been used to probe radial abundance gradients. These are Open Clusters (OCs, e.g. Randich et al. 2003, 2022; Magrini et al. 2010; Yong et al. 2012), HII regions (e.g. Balser et al. 2011; Esteban et al. 2017; Méndez-Delgado et al. 2022), young massive O and B stars (e.g. Daflon & Cunha 2004; Bragança et al. 2019), Classical Cepheids (CCs, e.g. Lemasle et al. 2007, 2008; Luck & Lambert 2011; Genovali et al. 2015; Kovtyukh et al. 2022), planetary nebulae (PNe, e.g. Maciel et al. 2003; Henry et al. 2010; Stanghellini & Haywood 2010, 2018), and also field stars with  \nprecise stellar ages (e.g. Anders et al. 2017; Santos-Peral et al. 2021).  \nThe plethora of information from all these tracers has been carefully interpreted by means of models, allowing us to understand fundamental properties for the Milky Way (MW) disk formation. Among these, the inside-out mechanism (e.g. Matteucci & Francois 1989; Chiappini et al. 2001; Schönrich & McMillan 2017), variable star formation e ffi ciency (SFE), i.e. higher in the inner regions than in the outer ones (e.g. Colavitti et al. 2009; Grisoni et al. 2018; Palla et al. 2020b) and radial gas flows (e.g. Portinari & Chiosi 2000; Spitoni & Matteucci 2011; Bilitewski & Schönrich 2012; Cavichia et al. 2014). Moreover, studies on tracers of the "old" stellar gradients have given important information on the impact of the process of stellar radial migration (e.g. Minchev et al. 2018; Willett et al. 2023).  \nOn the other side, the analysis of spectroscopic data from ground based surveys, such as the APOGEE (e.g. Hayden et al. 2015; Queiroz et al. 2020), the Gaia -ESO (GES, e.g. RecioBlanco et al. 2014; Kordopatis et al. 2015; Rojas-Arriagada et al. 2016), and the GALAH (Buder et al. 2019, 2021) ones, the accurate asteroseismic stellar ages (e.g. Pinsonneault et al. 2014, 2018; Miglio et al. 2021) and the kinematics and dynamical properties provided by the Gaia mission (Gaia Collaboration et al. 2016, 2018, 2021, 2023b) have pointed towards the existence of two sequences of stars in the [ α / Fe] 1 versus [Fe / H] abundance pattern in the MW disk: the so-called highα and lowα sequences.  \nTo explain the wealth of these data, Palla et al. (2020b, 2022) (see also Spitoni et al. 2019, 2021) suggested the that presence of this feature (also known as α -bimodality) may be connected to to a delayed ( ≳ 3 Gyr) accretion of gas. The latter forges the lowα sequence of stars, while highα stars are formed promptly in a gas infall episode occurring in the first phases of Galactic formation. This scenario is confirmed by several models and simulations of the evolution of galactic disks (e.g. Noguchi 2018; Grand et al. 2018; Mackereth et al. 2018; Buck 2020) which suggested that the bimodality may be strictly connected to a delayed accretion of gas of primordial / metal-poor chemical composition.  \nRecently, to reproduce the chemical abundances from Gaia DR3 Radial Velocity Spectrometer (RVS) spectra in the solar vicinity (Gaia Collaboration et al. 2023a; Recio-Blanco et al. 2023) and in particular a population of massive stars with evidences of a recent chemical impoverishment, Spitoni et al. (2023) suggested a novel scenario of chemical evolution in which the lowα population of stars is generated by two distinct gas accretion episodes, with the latter infall happening at very recent times ( < 3 Gyr of age). This scenario is constrained by the star formation histories for disk stars inferred from Gaia DR1 and DR2 color-magnitude diagrams (CMDs, Bernard 2017; Ruiz-Lara et al. 2020), which show evidence for short episodes of enhanced star formation (hereafter, SF) in recent times. In this model, the enhanced SF is triggered by the gas infall, which in turn causes the visible chemical impoverishment in the young stellar population as observed by Gaia RVS.  \nIn the light of these recent development for the MW disk formation and evolution, we exploit the OCs from the last data release of the Gaia -ESO survey (Randich et al. 2022) to investigate the late evolution of radial chemical gradients in the Galaxy. OCs are in fact considered excellent tracers of the chemical properties of the disk stellar populations of our Galaxy, including the spatial distribution of elemental abundances, especially when observed with high-resolution spectroscopy (Spina et al. 2022). To this regard, the Gaia -ESO survey (Gilmore et al. 2012, 2022; Randich et al. 2013, 2022), is the only survey performed on a 8 m-class telescope, which put specific focus on the population of Galactic OCs. Gaia -ESO targeted OCs over a wide range of ages, distances, masses, and metallicities, observing large unbiased samples of cluster candidates, with a well-defined selection function (Bragaglia et al. 2022; Randich et al. 2022).  \nOf this sample, Magrini et al. (2023) selected 62 OCs, with extended radial (up to R > 15 kpc) and age (up to 7 Gyr) ranges, analysing the shape of radial gradients in chemical elements spanning di ff erent nucleosynthetic origin (from Oxygen to Europium) and their time evolution. They found that the gradients of most of chemical elements, including the metallicity [Fe / H],  \ncan be better approximated with a two-slope shape, steeper in the inner regions and rather flat in the outer ones. Shallower gradient slopes in outer regions were also observed in other studies using di ff erent tracers (OCs, e.g. Carbajo-Hijarrubia et al. 2024; CCs, e.g. da Silva et al. 2023), suggesting a flat SFE law for large Galactocentric distances at variance with previously theorised (e.g. Grisoni et al. 2018; Palla et al. 2020b). In addition, Magrini et al. (2023) found that the youngest clusters in the sample (age < 1 Gyr) have lower metallicity than their older counterpart, even though the e ff ect could be mitigated by avoid considering stars with low surface gravity.  \nIn this work, we compare the above mentioned data sample with detailed chemical evolution models, which also account for the e ff ect of stellar radial migration. We start from well tested models under the revised two-infall scenario (Palla et al. 2020b; da Silva et al. 2023), which successfully reproduce data from high-resolution surveys in di ff erent Galactocentric regions (see also Spitoni et al. 2021) and were already adopted in the context of radial abundance gradients (Palla et al. 2020b; da Silva et al. 2023). To better investigate the observed behaviour in Gaia -ESO OCs, we then extend the comparison to the newly proposed scenario of the three-infall model (Spitoni et al. 2023), expanding this framework for the first time to the whole MW disk and discussing its feasibility in the context of the adopted dataset.  \nThe paper is organised as follows. In Section 2 we describe the Gaia -ESO OCs sample and the additional dataset adopted in this work. In Section 3, we present the model framework used, i.e. from the model scenarios to the nucleosynthesis and radial migration prescriptions. In Section 4 we present the comparison between the di ff erent model predictions and the observations, also discussing the implications of the obtained results. Finally, in Section 5 we draw some conclusions.', '2.1. Gaia -ESO open clusters': 'Gaia -ESO is a large public spectroscopic survey that observed the major components of our Galaxy with FLAMES@VLT from 2011-2018 (see Randich et al. 2022; Gilmore et al. 2022, for a full description of the survey). The final release dr 5.1 is public and available at the ESO website from June 2023. Gaia -ESO made use of FLAMES with both the high-resolution spectrograph UVES (operating at a resolving power, R = 47 000) and the medium-resolution spectrograph GIRAFFE (R ∼ 20 000). It observed open star clusters for about 30% of its 340 nights. The observed clusters cover a wide range in age, distance, mass and metallicity (see Randich et al. 2022), with an unbiased selection of cluster candidates. Each cluster was observed with both UVES and GIRAFFE. In particular, UVES spectra cover a wide spectral range from 480.0 nm to 680.0 nm (U580) or from 420.0 nm to 620.0 nm (U520). The large spectral interval, combined with the high signal-to-noise ratios (S / N) and with the high-resolution, have enabled an unprecedented characterisation of a large sample of open clusters, observed and analysed in a homogeneous way (Bragaglia et al. 2022; Hourihane et al. 2023). For the about 80 observed clusters, it was possible to obtain precise stellar parameters and abundances of more than 30 di ff erent ions, including those of elements belonging to all the main nucleosynthesis channels, from the lightest ones, such as Li to the heaviest one, such as Eu. The detailed chemistry of the open cluster sample in Gaia -ESO, combined with uniform ages  \nand precise distances from Gaia (e.g. Cantat-Gaudin et al. 2020) has been used, e.g., to calibrated age-sensitive abundance ratios, the so-called chemical clocks (Casali et al. 2019, 2020b; Viscasillas Vázquez et al. 2022), to investigate the nucleosynthesis of neutron-capture elements (Magrini et al. 2018, 2022; Van der Swaelmen et al. 2023; Molero et al. 2023), to study the evolution of Li abundance (Randich et al. 2020; Romano et al. 2021; Magrini et al. 2021b,a), and to study the shape and the evolution of the radial abundance gradients (Jacobson et al. 2016; Overbeek et al. 2017; Magrini et al. 2017; Spina et al. 2017, 2022; Magrini et al. 2023; da Silva et al. 2023).  \nIn the present work, we make us of the sample of open clusters used in Magrini et al. (2023, hereafter Ma23), in which abundances for 25 chemical elements were provided for the 62 Gaia -ESO OCs older than 100 Myr. The adopted cluster ages are the same as in Viscasillas Vázquez et al. (2022) and are derived from the homogeneous analysis of Gaia photometric and astrometric data by Cantat-Gaudin et al. (2020) by means of an artificial neural network trained on a set of objects with well-determined parameters in the literature. Here, we consider the guiding radius (R guide ) as tracer of the OCs Galactocentic distance. These are computed as defined by Halle et al. (2015, 2018), i.e. as the average between the minimum and maximum radius of the orbits. For the latter, we calculate them by using the G alpy code with the axisymmetric potential MWP otential 2014 (Bovy 2015). For further information on the sample cluster parameters, as well as on their distribution in di ff erent quantities, we refer to Viscasillas Vázquez et al. (2022). For each cluster, the membership analysis was done using both radial velocities from Gaia -ESO and proper motions and parallaxes from Gaia edr 3 (Gaia Collaboration et al. 2021), as described in Magrini et al. (2021b), Viscasillas Vázquez et al. (2022) and Jackson et al. (2022).', '2.1.1. A restricted sample to avoid observational biases': 'From a purely observational point of view, the cluster sample in Ma23 showed an unexpected behaviour regarding the evolution with time of the metallicity gradient: the youngest clusters (age < 1 Gyr) in the inner disk have lower metallicity than their older counterparts and they outline a flatter gradient below that of the older population  \nTo distinguish the real evolution of the gradient from possible spectral analysis e ff ects, Ma23 restricted the sample of member stars per cluster. In fact, an investigation of the internal abundances of each cluster showed that there are trends of [Fe / H] versus stellar parameters, in particular gravity and microturbulence. These trends are not specific to Gaia -ESO but are present in all the considered spectroscopic surveys (see Fig. 10 and Appendix in Ma23). This might be due to two e ff ects: the former is related to problems in modelling the atmospheres of giant stars, which a ff ect the spectral analysis of low-gravity giant stars, as noticed already in Casali et al. (2020a) and Spina et al. (2022), and the latter to the e ff ects of magnetic activity in young massive giant stars. Hence, the spectral analysis of giants (log g < 2.5) likely underestimates their [Fe / H] of about 0.1-0.2 dex due to the combination of these e ff ects.  \nTherefore, in a conservative approach, in order to preserve the actual mean abundance of clusters, we consider a restricted sample of member stars, where only stars with log g > 2.5 and ξ < 1.8 km -1 are taken into account to compute the mean cluster abundances. The average abundances used in the present work, restricted in stellar parameters, are reported in Tab. C.1. In the Table we also report the di ff erent cluster parameters (i.e. age,  \ndistances, orbital parameters), which are the same as the ones introduced earlier in 2.1 for the original Gaia -ESO OC sample.', '2.2.1. Gaia -ESO field stars': 'We select a sample of about 3800 field stars in Gaia -ESO dr 5.1 following the criteria of Viscasillas Vázquez et al. (2022), to which we refer for a complete description. Here are the basic steps of our selection and the computation of their ages. Our set of field stars is composed by both field stars (GES\\_FLD keywords GES\\_MW for general MW fields, GES\\_MW\\_BL for fields in the direction of the Galactic bulge, GES\\_K2 for stars observed in Kepler2 (K2) fields, GES\\_CR for stars observed in CoRoT fields) and non members of the 62 open clusters considered in this work (age > 100 Myr). We applied some further quality selections: SNR > 20; eT e ff < 150 K, e log g < 0 . 25, e [Fe / H] < 0 . 20 and e ξ < 0.20 km s -1 . We also apply a further cut in abundance errors, considering only those values that have an eA ( El ) < 0 . 1.  \nThe selection function adopted in the Gaia -ESO survey for UVES observations (see Stonkut˙e et al. 2016) favours the main sequence turn-o ff (MSTO), which constitute the majority of the sample. By construction, the sample of field stars has a limited extent in Galactocentric distances, and is, therefore, used in our analysis for comparison with the model in the [X / Fe] vs [Fe / H] plane.', '2.2.2. Classical Cepheids in Da Silva et al. (2023)': "da Silva et al. (2023, hereafter DS23) provided the largest (1118 spectra, 356 objects) and most homogeneous spectroscopic sample for Galactic CCs with measured metallicity from optical high-resolution, high-S / N spectra. The sample is distributed across the four Galactic quadrants and it ranges from the inner (R ∼ 5 kpc) to the outer (R ∼ 25 kpc) disk. For the distances, measurements were based either on trigonometric parallaxes from Gaia DR3 or on near-infrared period-luminosity relations (Ripepi et al. 2022). Due to the steadily variation in target's physical properties due to their natural radial oscillations, special care was dedicated to the estimate of the di ff erent atmospheric parameters ( Tef f , log g , ξ ), which were verified using di ff erent approaches and / or diagnostics (see also da Silva et al. 2022).  \nIn this work, we take advantage of the CCs sample from DS23 to have an additional observational probe for present-day metallicity gradients, in addition to the young OCs from Ma23.", '3. Chemical evolution of the Milky Way disk': 'In this Section, we present the main assumptions and features of the multi-zone chemical evolution models adopted in this work. In particular, in 3.1 we provide the details of the revised two-infall model proposed by Palla et al. (2020b) (see also DS23), whereas in 3.2 we describe the details of the threeinfall framework from Spitoni et al. (2023), which we expand throughout this work.  \nFor both of the chemical evolution models listed above, the basic equations that describe the chemical evolution of a given element i are:  \n˙ Gi (R , t ) = -ψ (R , t ) Xi (R , t ) + Ri (R , t ) + ˙ Gi , in f (R , t ) + ˙ Gi , Rf , (1)  \nTable 1. Summary of the main parameters of the two-infall model (2INF) adopted in this study.Notes. All the above parameters are the same to the ones adopted in DS23. The negative sign on the radial flow speed indicates inward flows.  \nwhere Gi (R , t ) = Xi (R , t ) G (R , t ) is the fraction of the gas mass in the form of an element i and G (R , t ) is the fractional gas mass. Xi (R , t ) represents the abundance fraction in mass of a given element i , with the summation over all elements in the gas mixture being equal to unity.  \nThe first term on the right-hand side of Eq. (1) corresponds to the rate at which an element i is removed from the ISM due to star formation. The star formation rate (hereafter, SFR) is parametrised according to the Schmidt-Kennicutt law (Kennicutt 1998):  \nψ (R , t ) = ν Σ gas (R , t ) k , (2)  \nwhere Σ gas is the surface gas density, k = 1 . 5 is the law index and ν is the star formation e ffi ciency (SFE).  \nRi (R , t) (see Palla et al. 2020a for the complete expression) takes into account the nucleosynthesis from low-intermediate mass stars (LIMS, m < 8M ⊙ ), core collapse (CC) SNe (Type II and Ib / c, m > 8M ⊙ ) and Type Ia SNe. For these latter, we assume the single-degenerate scenario and in particular the delay-timedistribution (DTD) by Matteucci & Recchi (2001). This choice can be considered a good compromise to describe the delayed pollution from the entire Type Ia SN population as it enables us to obtain abundance patterns that are similar to those obtained with other DTDs (see Palla 2021 for details). Ri (R , t ) output is also weighted by the initial mass function (IMF). Here, we adopt the IMF by Kroupa et al. (1993), which is preferred to reproduce the characteristics of the MW disk (Romano et al. 2005).  \nThe last term of Eq. (1) refers to radial inflows of gas, which here are implemented following Portinari & Chiosi (2000) (see also Palla et al. 2020b for a detailed description of the implementation). In our models, we use a constant speed pattern, i.e. with vflow = 1 kms -1 across all radii, as suggested by Palla et al. (2020b). Low speeds, i.e. small radial inflow motion, are also suggested by previous chemical evolution studies (e.g. Bilitewski & Schönrich 2012; Mott et al. 2013; Vincenzo & Kobayashi 2020) as well as observations of external galaxies (e.g. Wong et al. 2004; Di Teodoro & Peek 2021).  \nIn general, we ignore the e ff ect of Galactic winds on chemical evolution of the MW disk. In fact, by studying the Galactic fountains originated by the explosions of Type II SNe in OB associations, Melioli et al. (2008, 2009) and Spitoni et al. (2009) found that metals fall back to approximately the same Galactocentric region from where they were ejected. Moreover, Spitoni et al. (2009) computed the typical timescale of the fallback of this material, finding a value of 0.1 Gyr. These results were also later supported by cosmological simulation of galactic discs of virial mass > 10 11 M ⊙ , i.e. encompassing the MW, which showed that the majority of the mass ejected by the disc is reaccreted on short timescales and close to the ejection location (e.g. Hopkins et al. 2023 and references therein). Therefore, galactic winds are likely to produce a negligible e ff ect on the global chemical evolution of the Galaxy.', '3.1. The two-infall model': 'In the two-infall model formalism, the model assumes two consecutive gas accretion episodes feeding the MW disk, forming the so-called highα and lowα sequences observed in the Galactic disk. Therefore, the third term on the right-hand side of Eq. (1) can be expressed in this way:  \n˙ Gi , in f (R , t ) = A (R) Xi , 1 in f (R) e -t τ 1 + + θ ( t -tmax , 1) B (R) Xi , 2 in f (R) e -t -t max , 1 τ 2 , (3)  \nwhere Gi , in f (R , t ) is the infalling material in the form of element i and Xi , Jin f is the abundance of the same element for the J -th infall. τ 1 and τ 2 are the timescales of the two infall episodes, while tmax , 1 indicates the time of maximum infall, which is also the delay between the first and the second infall episodes. θ is the Heavyside step function, while A (R) and B (R) coe ffi cients are set to reproduce the present-day total surface mass density of the highα and lowα disks at di ff erent Galactocentric radii. These latter are assumed to have exponential profiles, in this way:  \nΣ J (R) M ⊙ pc -2 = Σ 0 , J e -R / R d , J , (4)  \nwhere the disk scale length R d , J is 3.5 kpc for the lowα disk ( J = 2) and 2.3 kpc for the highα one ( J = 1, see Palla et al. 2020b and references therein). The quantities Σ 0 , J are tuned to reproduce the total surface mass density in the solar neighbourhood as provided by McKee et al. (2015) of 47 . 1 ± 3 . 4 M ⊙ pc -2 .  \nIn this work, we take advantage of the model prescriptions adopted in Palla et al. (2020b) and later revised in DS23. The model is thus a revised version of the two-infall paradigm (see also Spitoni et al. 2019) in which two consecutive gas accretion episodes are separated by a delay of tmax , 1 = 3 . 25 Gyr. The fairly large delay relative to the "classical" two-infall paradigm (of 1 Gyr, see Chiappini et al. 1997; Romano et al. 2010), allow us to reproduce [ α / Fe] vs. [Fe / H] abundance diagrams (e.g. Palla et al. 2020b; Spitoni et al. 2021) throughout the MW disk as well as stellar ages trends (e.g. Spitoni et al. 2019; Nissen et al. 2020) in the solar vicinity.  \nGoing into more detail, in the first infall (forging the highα sequence) the timescale of gas accretion is fixed at τ 1 = 1 Gyr at all radii, with also a fixed SFE of ν = 2 Gyr -1 . For the second gas-infall episode (forming the lowα sequence), the timescale for gas accretion τ 2 increases linearly with radius according to the inside-out scenario (following Romano et al. 2000; Chiappini et al. 2001 law). In order to reproduce the slope observed in radial abundance gradients of CCs in DS23, as well as the gradients in other physical quantities such as SFR and gas density (e.g. Stahler & Palla 2005; Nakanishi & Sofue 2003, 2006), the SFE of the second infall episode is variable depending on the galactocentric radius, with values between ν = 5 Gyr -1 (at R = 4  \nkpc) and ν = 0 . 4 Gyr -1 (from R > 12 kpc). In addition, inward radial gas flows with a constant velocity of vflow = 1 km s -1 are also needed to reproduce the gradients, as discussed above.  \nA summary of all the main parameters adopted in this work for the two-infall model are listed in Table 1.', '3.2. The three-infall model': 'In their recent work, Spitoni et al. (2023, hereafter Sp23) proposed a new chemical evolution framework for the Galactic disk components, constrained by the star formation histories inferred from CMD analyses in Gaia DR1 and DR2 (Bernard 2017; RuizLara et al. 2020). These works revealed enhanced SF activity within the last 2-3 Gyrs. This is mimicked in the chemical evolution model by including a recent (age < 3 Gyr) gas infall episode, which triggers this enhanced SF at late times. In this way, Sp23 were also able to reproduce the new abundance ratios provided by the General Stellar Parametriser-spectroscopy module for the Gaia DR3 (Gaia Collaboration et al. 2023a; RecioBlanco et al. 2023), which show a chemical impoverishment in the young population of stars in the solar neighbourhood.  \nIn this work, we extend the approach presented in Sp23 to the whole MW disk, in order to investigate abundance gradients.  \nIn this model, the functional form of the gas infall rate is:  \n˙ Gi , in f (R , t ) = A (R) Xi , 1 in f (R) e -t τ 1 + + θ ( t -tmax , 1) B (R) Xi , 2 in f (R) e -t -t max , 1 τ 2 + + θ ( t -tmax , 2) C (R) Xi , 3 in f (R) e -t -t max , 2 τ 3 ,  \n(5)  \nwhere τ 3 is the timescales of the third gas accretion episode, tmax , 2 is the Galactic time associated to the start of the third infall, and C (R) is the coe ffi cient to reproduce the present-day total surface density of the third accretion episode Σ 3. Here, the sum between the latter and the total surface density of the second infall ( Σ 2) is equal to the density profile as described in 3.1 for the lowα disk in the two-infall model. All the other variables are as in Eq. (3). In fact, the new model uses the framework presented in 3.1, but split the lowα sequence into two distinct gas accretion episodes in order to mimic the recent enhanced SF activity. Therefore, we will leave unchanged all the parameters adopted for the first two infall episodes, such as the infall timescales and the SFEs at di ff erent Galactocentric radii, as well as the IMF (from Kroupa et al. 1993).  \nFor what concern the third additional gas accretion, instead, we set its starting time at tmax , 2 = 11 Gyr (as in Sp23) for all the Galactocentric radii. For the other parameters, i.e. the SFE ν 3, the timescale of infall τ 3 and the total surface mass density accreted by the third gas infall, we test di ff erent parametrisations, which are listed in Tab. 2:  \n- 1. First, we adopt a setup very similar to the one proposed by Sp23 (3INF-1), which is shown in Tab. 2 upper row. Here, the parameters of the third infall, i.e. tmax , 2, τ 3, Σ 2 /Σ 3 are the same of the latter paper. Concerning the SFE ν 3, rather than fixing the value adopted in Sp23, we fix the proportion between the SFE in the second and third infall episode to be similar to the one used in that paper;\n- 2. to test further the viability of the three-infall scenario in the context of radial gradients, we allow to vary the third infall parameters relative to the values adopted in the model 3INF1. While leaving constant the starting time of the third infall (to be consistent with SF peaks as show by, e.g., Ruiz-Lara  \nTable 2. 3rd infall parameters adopted in the models in this paper.  \nNotes. For 3INF-1, τ 3 and Σ 2 /Σ 3 parameters are from Spitoni et al. (2023). For ν 3, we consider a similar proportion to the one used by Spitoni et al. (2023) for the solar vicinity, i.e. R = 8 kpc.  \net al. 2020) we act on all the other physical parameters, i.e. the infall timescale τ 3, the SFE ν 3 and for the ratio between the baryonic mass accreted by the second and third infall Σ 2 /Σ 3. The parameter for this setup (3INF-2) are shown in Tab. 2 bottom row.', '3.3. Nucleosynthesis prescriptions': 'The nucleosynthesis prescriptions and the implementation of the stellar yields are fundamental ingredients for chemical evolution models. LIMS, massive stars and Type Ia SNe play a fundamental role in shaping the [X / Fe] vs. [Fe / H] abundance patterns as well as radial abundance gradients of the elements of the Periodic Table.  \nIn this work, we mainly adopt the prescriptions listed below:  \n- -for LIMS, we use the yield set from Ventura et al. (2013, 2018, 2020), also comprising the domains of super-AGB stars (6 < m / M ⊙ < 8 -9) and supersolar metallicities.\n- -for massive stars we use the stellar yields from Limongi & Chie ffi (2018). In particular, we adopt the "mixed vrot 2 set", i.e. the one used in the best model (MWG-12) in Romano et al. (2019) for the MW disk.\n- -for Type Ia SNe we adopt the stellar yields from Iwamoto et al. (1999) (W7 model) which are extensively used in chemical evolution literature (e.g. Romano et al. 2010; Prantzos et al. 2018; Palla et al. 2020b among others).  \nHowever, underlying uncertainties in stellar evolution and nucleosynthesis theory may limit our chances of getting a firm grasp on the evolutionary scenario for the Galaxy.  \nFor this reason, we also test other stellar yields for massive stars and Type Ia SNe. In particular, we run additional models taking advantage of the stellar yields from Kobayashi et al. (2006, 2011) for massive stars. Concerning Type Ia SNe, instead, we also adopt either yields from Leung & Nomoto (2018) (benchmark model) or Leung & Nomoto (2020) (bubble detonation pattern model), which represent some of the most recent Type Ia SN models for di ff erent progenitor classes, i.e. nearChandrasekhar mass white dwarfs (nearMCh , Leung & Nomoto 2018) and sub-Chandrasekhar mass white dwarfs (subMCh , Leung & Nomoto 2020, see Kobayashi et al. 2020; Palla 2021 for more details).', '3.4. Accounting for migration and observational uncertainties': 'The overall picture on MW shows also evident signatures of stellar migration, both on the theoretical and observational sides  \n(e.g. Schönrich & Binney 2009; Minchev et al. 2011, 2018; Kordopatis et al. 2015), even when specifically focusing on OCs (e.g. Anders et al. 2017; Spina et al. 2021; Myers et al. 2022). Therefore, we include in the models stellar radial migration prescriptions from the literature to account for this phenomenon. In particular, we implement migration in the chemical evolution model by adopting the approach already tested in Palla et al. (2022) for MW disk stars. It is worth noting that OCs are more massive than single stars, and therefore the e ff ect of the interactions with perturbing structures should be in principle less pronounced than for field stars (e.g. Zhang et al. 2021; Ma23), especially for young-intermediate ages (see Viscasillas Vázquez et al. 2023). However, a model for radial migration valid for field stars is appropriate within our work, as it gives at least a robust upper limit on the e ff ect of radial migration on OCs, which is a necessary mechanism to properly explain their radii, abundances and ages trends in the Galactic disk.  \nIn the following, we provide some details on the implementation, which is extensively described in Frankel et al. (2018). Here, migration is seen as a result of a di ff usion process, which is the e ff ect of repeated and transient torques on stars by features such as spiral arms or a bar, and it is treated in a parametrical way. Following Sanders & Binney (2015) and adapting their parameterization to a Galactocentric radius coordinate, the probability for a star to be currently at a Galactocentric radius R f , given that it was born at a radius R0 and at a certain age can be written as:  \nln p (R f | R0 , Age ) = ln( c 3) -(R f -R0) 2 2 σ RM ( Age / 10 Gyr) , (6)  \nwhere σ RM is the radial migration strength and c 3 a normalization constant ensuring that stars do not migrate to negative radii. For σ RM we adopt a value of 3, as found by Frankel et al. (2020) as a result of their Bayesian fitting procedure of APOGEE red clump stars in the lowα disk. The value above mentioned refers to the churning 3 strength, as we already account for blurring effect by considering OCs guiding radius (R guide ) rather than their present-day Galactocentric radius R GC .  \nIt is worth noting that the migration framework adopted is in 1D, i.e. it is not considering azimuthal variations in migration strength. In this way, it can be fully integrated within the chemical evolution models described previously in this Section. The 1D assumption is robust, as various works showed rather small azimuthal abundance variations in galaxies ISM (e.g. Kreckel et al. 2019). Another possible limitation of the adopted migration model is the assumption of no radial or temporal dependence on the migration strength. However, including these dependencies means adding further and uncertain assumptions on the inventory of speeds and strengths of spiral and bar patterns during Galactic evolution (see Frankel et al. 2018). In any case, the calibration of the model parameters on a sample of lowα disc stars extending up to radii ≃ 15 kpc (Frankel et al. 2020) allows for a radial and temporally averaged estimate of migration in a radial range that comprises all the OCs described in Section 2, except one (Br29).  \nWe also account for the e ff ect of observational uncertainties on the chemical abundances of the predicted stellar populations in our model.  \nFig. 1. Time evolution of the radial [Fe / H] gradient as predicted by the two-infall model. Filled circles with errorbars are the OC sample by Magrini et al. (2023), which are divided in three age bins: young ( Age < 1 Gyr, blue points), intermediate (1 < Age / Gyr < 3, green points) and old ( Age > 3 Gyr, red points). Solid lines are the results for the [Fe / H] gradient as predicted by the two-infall model at 0.5 Gyr (blue lines), 2 Gyr (green lines) and 5 Gyr (red lines). In this plot, and in the following Figures, we use the guiding radius of the orbit, computed as the average between the Apogalacticon and Perigalacticon radii (see 2.1), as an indication of the location of each cluster in the disk.  \n<!-- image -->  \nIn particular, we add at each Galactic time t a random error to the abundances of the stars formed at t (see also Spitoni et al. 2019; Palla et al. 2022). In this way, we have for each chemical element a "new abundance", defined as:  \n[X / H] new ( t ) = [X / H]( t ) + N ([X / H] , σ [X / H]) , (7)  \nwhere N is a random function with Normal distribution. In order to have fair comparison with the Ma23 data set adopted in this study, the standard deviation σ [X / H] corresponds to average spread observed within the OCs in the age intervals of interest, i.e. 0 . 1 < Age / Gyr < 1, 1 < Age / Gyr < 3, 3 < Age / Gyr < 7.', '4. Results': 'In this Section, we show the results of the comparison between the models presented in Section 3 and the Ma23 OC sample. In particular, in 4.1 we look at the prediction of the two-infall model. In 4.2 we instead show the gradients obtained by means of the three-infall scenario in the light of the data at disposal, discussing such a scenario in 4.3.', '4.1. Comparing observed gradients with the standard two-infall scenario': 'In Fig. 1, we show the time evolution of the [Fe / H] as predicted by the two-infall model and compared to the OC data presented in Ma23. From this sample, we remove the Blanco 1 cluster due to its extremely high internal spread observed in metallicity ( σ > 0 . 5 dex), which may hide problems in abundance derivation of their members. Coming back to Fig. 1, we see that the model captures the general trend of the data: the gradient slope clearly decreases going towards larger radii, in agreement with the trend shown by OCs. This result confirms the conclusion by DS23 of a flattening of the gradient at R ≳ 12 kpc, which requires a flat behavior of the SFE at large Galactocentric radii.  \nHowever, if we focus on di ff erent age bins, we note that OCs with Age < 1 Gyr (blue points) are clearly below the prediction from the model. Moreover, we note that the metallicity of these OCs is lower than that of older clusters, i.e. the ones with 1 < Age / Gyr < 3 (green points). Therefore, such a decrease in metallicity with decreasing age cannot be reproduced by genuine chemical evolution model predictions in a scenario of continuous SF, as the one of the two-infall model in the age ranges investigated in this study: subsequent stellar generations are progressively enriching the ISM in metals, as demonstrated by model lines in Fig. 1.  \nFor this reason, we try to explore if the inclusion of e ff ect of stellar migration and abundance uncertainties / spread within OCs (see 3.4 for details) may reconcile predictions and observations. The results are shown in Fig. 2. Here, the density plots represent the probability of finding a star with an abundance [X / H] at a Galactocentric radius R in a certain age bin, i.e. Age > 3 Gyr (left panel) 1 < Age / Gyr < 3 (central panel) and Age < 1 Gyr (right panel), according to the predictions by the two-infall model. In order to highlight the spread caused by the inclusion of the di ff erent e ff ects, the density plot in Fig. 2 (as well as subsequent Figures) is represented in log scale. By looking at the left and central panels of the Figure, we can see that the model predictions generally well capture the observed slope and spread within Ma23 OC sample, with just the exception of few objects with 10 ≲ R / kpc ≲ 12 in the intermediate age bin. On the other hand, the right panel clearly shows that young clusters are clearly overestimated by model predictions at di ff erent radii. This is also demonstrated by the plotted solid and dashed contour lines, showing the limits in which are contained the 68 and 95% of the predicted stars, respectively. Almost all the observed clusters are in fact outside the 2 σ of the distribution.  \nHowever, as noted by Ma23, abundances in stars with lower surface gravity log g and higher mictroturbolence parameter ξ are susceptible to artifacts in stellar spectral analysis (see 2.1.1). In turn, this may a ff ect the reliability of the obtained OC abundances especially for young clusters, which are more prone to contain such stars. Therefore, we decide to remove stars with log g < 2 . 5 and ξ > 1 . 8 km s -1 , as suggested by Ma23, to minimise the possible bias introduced in stellar abundance determination within clusters. Despite such a choice costs a not negligible number of OCs for the data-model comparison, this is the only way we can assure robust estimations for di ff erent chemical abundances.  \nIn Fig. 3, we show the comparison between the models already shown in Fig. 2 and the Ma23 OC sample with imposed cuts on individual stellar parameters (hereafter, restricted sample). As mentioned above, the old and intermediate age bins are barely a ff ected by the cut in stellar parameters and therefore the data-model comparison the left and central panel is very similar to what already shown in Fig. 2. For the right panel, we see instead that the restricted sample shows in general larger metallicities than those observed in the full OC sample. Nonetheless, the observed [Fe / H] are still clearly overestimated by the prediction of two-infall model, with most of the clusters still outside the 2 σ of the predicted distribution.  \nThis is also highlighted in Tab. 3, where we show intercepts and slopes of gradient linear fits for the full and restricted OC samples and the two-infall model predictions in di ff erent age bins. The fits are obtained in the radial range below R < 12 kpc: this choice is done to avoid the influence of the gradient flattening at large radii, which is found to start around this Galactocen-  \ntric distance (see, e.g. Carraro et al. 2007; Magrini et al. 2017; Donor et al. 2020; DS23). For the young age bin, the restricted sample shows a significant increase in the fit intercept relative to the full sample, but still this is around 0.2 dex lower than what predicted by the two-infall model. Therefore, even by considering stellar migration and abundance uncertainties e ff ects, these are not su ffi cient to explain the low [Fe / H] imprinted in the data sample.  \nWe also checked whether uncertainties on cluster ages can be a source of bias for such a result. We consider either the cases of (i) "internal" uncertainties for the age derivation method adopted in this work and (ii) adoption of di ff erent isochrones grids and / or di ff erent photometric datasets to derive ages. For (i), we perform 1000 Monte Carlo resamplings of the age of the clusters according to their uncertainties. We consider an age uncertainty of log( Age ) = 0 . 2 dex, as found in Cantat-Gaudin et al. (2020) for clusters older than log( Age / yr) > 8 . 5 in their validation sample 4 . Even by randomly perturbing the cluster ages according to the uncertainties, we still find a decrease > 0 . 15 dex (precisely, from 0.70 to 0.52 dex) between the intercepts in the intermediate and young age bins relative to the ones shown in Table 3 for the restricted sample. For (ii), instead, we refer to Je ff ries et al. (2023), who compared di ff erent sets of age determinations for OCs in Gaia -ESO, including a large fraction of the OCs used in this work. We find that the di ff erence between the mean literature ages (from Je ff ries et al. 2023) and the ages adopted in this work is always lower than the uncertainty of log( Age ) = 0 . 2 dex adopted for our resampling test: therefore, the distribution of gradient in di ff erent age bins will remain similar even when changing the ingredients in the age derivation.  \nTo further inspect the low [Fe / H] abundances by OCs at young ages, we also consider the metallicities derived for the CCs sample in DS23 (see Section 2). It is worth reminding that the latter sample provided homogeneous derivation of Fe abundances relative to the OC sample adopted here, preventing any additional observational bias. Both the datasets are shown in Fig. 4, together with the prediction for the young age bin by the two-infall model. It is worth noting that here and in the following Figures, we avoid showing the contours of the predicted stellar distribution to avoid overcrowding the plots. Despite the much larger abundance spread observed in DS23 sample, the Figure shows that the bulk of Cepheids observed in DS23 up to a radius about ∼ 10 kpc are also underabundant relative to the model predictions. For larger radii instead, the very large spread in [Fe / H] and the relevant abundance uncertainties for distant Cepheids (we do not show them here for sake of readability, but they arrive up to 0.25 dex), do not allow us to draw strong conclusions, despite the hint of a flattening of the radial metallicity gradient at large Galactocentric distances remains (see DS23 and references therein). In any case, the superposition between the two datasets, which we remind are obtained using di ff erent gradient tracers, strengthen the findings described in this Section.  \nAs described in 3.3, we also run additional simulations using di ff erent yields for massive stars and Type Ia SNe. In this way, we test the dependence of the observed present-day gradient overestimation on nucleosynthetic calculations.  \nFig. 2. Time evolution of the radial [Fe / H] gradient as predicted by the two-infall model at Age > 3 Gyr (left panel), 1 < Age / Gyr < 3 (central panel), and Age < 1 Gyr (right panel) for the two-infall model, including stellar migration and abundance uncertainties (see 3.4). Density plots show the normalised density (in log scale) of stars as predicted by the model in given Galactocentric bin of 0.2 kpc width. White contour lines show the limits within are contained the 68% (solid) and the 95% (dashed) of the predicted stellar distribution in a given radial bin. Solid lines show the results for the [Fe / H] gradient as predicted by the genuine chemical tracks, i.e. the ones shown in Fig. 1, at Ages = 0.5 Gyr (cyan line), 2 Gyr (green line) and 5 Gyr (magenta line). Data are the same as in Fig. 1  \n<!-- image -->  \nFig. 3. Same as Fig. 2, but with data restricted to OCs members with log g > 2.5 and ξ < 1 . 8 km s -1 (restricted sample).  \n<!-- image -->  \nTable 3. Comparison between OC sample by Magrini et al. (2023) and two-infall model [Fe / H] gradient slopes and intercepts for R < 12 kpc at di ff erent ages. Data fits are shown for both the full OC sample and the one restricted to OCs members with log g > 2.5 and ξ < 1 . 8 km s -1 .  \nResults for these runs are shown in Appendix A. They highlight that the theoretical predictions with di ff erent stellar yields are very similar, with an overestimation of the present-day gradient in the two-infall model and a similar evolution of the predicted gradient through cosmic time. This denotes that the conclusions of the evolutionary scenarios on the [Fe / H] gradient are marginally a ff ected by the nucleosynthesis prescriptions. Therefore, in the rest of paper, we will proceed by using the sets of yields adopted throughout this Section.', '4.2. A late time metal dilution: the three-infall model': 'To explain the unexpected decrease in the [Fe / H] gradient at late times, we explore the scenario proposed by Sp23 for the solar vicinity, where a late-time burst of SF (suggested by Gaia CMD analysis, e.g. Ruiz-Lara et al. 2020) is fueled by gas  \naccretion, which in turn causes a metal dilution in the ISM gas. In particular, we extend such a scenario to the whole disk than the solar neighbourhood.  \nWe start by adopting a setup for the third gas accretion very similar to the one by Spitoni et al. (2023), which is shown in Table 2 upper row. Here, the parameters of the third infall are the same of the latter paper, except for the SFE ν 3, for which we fix the proportion between the SFE in the second and third infall episode to be similar to the one used in Spitoni et al. (2023). In this way, we are able to preserve the gradient slope observed in the data at di ff erent ages: if we apply a flat SFE with radius for the third infall, this will significantly decrease the gradient slope, which is not really seen in the observations (see Table 3). Finally, at variance with Sp23, we adopt a primordial chemical composition for the third gas infall. However, we also perform an  \nFig. 4. Same as Fig. 3, but showing only the results for the young age bin ( < 1 Gyr), with the addition of Classical Cepheids data sample from da Silva et al. (2023) (pink filled circles).  \n<!-- image -->  \nadditional run with an infall enrichment as in Sp23 (1 / 5 enriched with an abundance pattern as the one predicted for the highα phase at [Fe / H] = -0 . 75 dex), finding negligible di ff erences.  \nIn Fig. 5, we show the predicted evolution of the [Fe / H] gradient by the three-infall model with parameters as for the 3INF1 setup. The results are shown for the last 3 Gyr, i.e. the ages at which this late gas accretion episode is actually acting and changing the gradient evolution. Fig. 5 shows that the e ff ect of metal dilution by the gas accretion is too strong. In particular, the upper panel highlights that we miss to reproduce the metal-rich stars in the intermediate age bin. This is due to the very strong gas dilution happening 2.7 Gyr ago, which leads the bulk of stellar production at subsolar metallicity, even at small Galactocentric radii. This is instead not seen in OC observations which, despite of the spread seen at individual radial distances, only marginally cover the region where the model expects most of the data.  \nThis is also reflected by looking at Fig. 6, where we show the Age-[Fe / H] relation as predicted by the model 3INF-1 in di ff erent radial bins, i.e. inner (R < 7 kpc, left panel), solar (7 < R / kpc < 9, central panel) and outer (R > 9 kpc, right panel). Even accounting for migration and abundance uncertainties e ff ects (see density bins in the Figure), the prominent metal dilution prevents a good agreement with the data. This is especially evident at small Galactocentric radii (left panel), where the model fails to reproduce the age-metallicity trend within the last 3 Gyr. The situation is less dramatic at larger radii, where a large fraction of the OCs falls within the range of values allowed by the predictions. However, in the central and right panels we also observe that the genuine chemical evolution predictions by the model at 8 and 10 kpc (solid lines in central and left panels, respectively), slightly underestimate the metal content observed in young clusters. The slight underestimation of the present-day gradient by the model 3INF-1 is in fact observed in Fig. 5 lower panel.  \nTo test further the viability of the three-infall scenario in the context of radial gradients, we allow to vary the third infall parameters relative to the values proposed in Sp23. While leaving constant the starting time of the third infall (to be consistent with SF peaks as show by, e.g., Ruiz-Lara et al. 2020) we act on all the other physical parameters, i.e. the infall timescale τ 3, the SFE  \n<!-- image -->  \nFig. 5. Time evolution of the radial [Fe / H] gradients at 1 < Age / Gyr < 3 (upper panel) and Age < 1 Gyr (lower panel) for the model 3INF-1, including stellar migration and average OC spread (see 3.4). Lines show the results for the [Fe / H] gradient as predicted by the model at Ages = 0.5 Gyr (cyan line), 1.8 Gyr (green dashed line) and 2.8 Gyr (green solid line). Data are the same as in Fig. 3.  \n<!-- image -->  \nν 3 and for the ratio between the baryonic mass accreted by the second and third infall Σ 2 /Σ 3, as already shown in Tab. 2 bottom row.  \nIn addition, we allow a mild chemical enrichment for the infalling gas during the third gas accretion, with half of the gas enriched at a level of [Fe / H] = -0.75 dex with abundance pattern as the one predicted for the highα phase at that metallicity. This assumption is justified in the light of two possible invoked physical mechanisms behind the observed recent peak in SF (e.g. Isern 2019; Mor et al. 2019; Ruiz-Lara et al. 2020). On one hand, it is suggested a tight connection between the SF peak and the last pericentric passage of Sagittarius dSph (Ruiz-Lara et al. 2020; Roca-Fàbrega et al. 2021): the gas retained by Sagittarius after its first encounter with the Galaxy may have been definitely stripped in its second pericentric passage, contributing to the gas accretion together with steady cooling flow of gas from the hot corona. On the other side, Nepal et al. (2024) proposed that the peak in MW disk SF was triggered by intense MW bar activity, which was shown to trigger enhanced SF in galaxies both from a theoretical (e.g. Baba & Kawata 2020) and observational (with Integral Field Spectroscopy on local galaxies, e.g. Lin et al. 2020) points of view. Without discriminating between the two  \nFig. 6. Age-[Fe / H] relation for the model 3INF-1 in di ff erent radial ranges, i.e. R < 7 kpc (left panel), 7 < R / kpc < 9 (central panel) and R > 9 kpc (right panel), including stellar migration and OC spread (see 3.4). The density plot show the normalised density of stars (in log scale) as predicted by the model in a given age bin of 0.25 Gyr width. Lines are genuine chemical tracks at 6, 8 and 10 kpc (solid) and 4 and 12 kpc (dashed). Filled circles represent the restricted sample within OCs in Magrini et al. (2023).  \n<!-- image -->  \nscenarios, it is therefore likely that a mild chemical enrichment have to be present in this late stage of Galactic evolution.  \nIn Fig. 7 we show the predicted evolution of the [Fe / H] gradient by the three-infall model with parameters as for the 3INF-2 setup. In the upper panel, showing the metallicity gradients in the age bin 1-3 Gyr, we note that the dilution e ff ect is milder relative to what seen in Fig. 5. This is highlighted by the green solid and dashed lines, showing that the chemical evolution track before and almost 1 Gyr after the onset of the third infall, which have a metallicity di ff erence of around 0.2 dex. The smaller dilution allows us to capture the trend and spread observed in OCs in this age range: here, the location of data points broadly correspond to the regions with predicted higher probability density of stars, at variance with what happens in Fig. 5. The good agreement between OC data and the 3INF-2 model is also seen in Fig. 7 lower panel, with all the OCs in the youngest age bin falling within the range of values allowed by the model.  \nAs done for model 3INF-1, we also compare the observed age-metallicity relations at di ff erent radii with the prediction of the model 3INF-2. These are shown in Fig. 8. As already noted in Fig. 7, here the dilution by the third gas accretion episode is much less prominent relative to the one seen in Fig. 6 for the model 3INF-1. This di ff erence can be explained by the larger infall timescale for model 3INF-2: here, the slower gas accumulation allows dying stars to progressively pollute with metals the ISM, leaving a more prolonged but less marked decrease in metallicty. This behaviour better reproduces the observed age-metallicity trends, as no sharp dilution e ff ects are seen in our restricted sample. Moreover, the slightly larger SFE ν 3 in model 3INF-2 allows a better agreement between the chemical evolution tracks at di ff erent radii (solid and dashed lines in Fig. 8) and the metal content in young clusters, which is instead underestimated by model 3INF-1.  \nTo probe even further the three-infall scenario, we looked also at other chemical elements than Fe. However, in doing this we are more prone to the intrinsic uncertainties in models related to stellar nucleosynthesis, which are less significant when probing the Fe gradient (see Appendix A). To limit as much as possible this problem, we focus on two α -elements whose  \n<!-- image -->  \nFig. 7. Same as Fig. 5, but for the model 3INF-2.  \n<!-- image -->  \n[ α / Fe] vs. [Fe / H] abundance patterns are generally well reproduced throughout the whole MW metallicity range with the main sets of yield adopted in this work (see 3.3), i.e. O and Si (see  \nFig. 8. Same as Fig. 6, but for the model 3INF-2.  \n<!-- image -->  \nPrantzos et al. 2018; Romano et al. 2019). Other α -elements, such as Mg, are severely underestimated when adopting massive stars Limongi & Chie ffi (2018) prescriptions (Prantzos et al. 2018; Palla et al. 2022), while most of Fe-peak elements also suffer of the additional uncertainty on the progenitor nature of Type Ia SNe, which severely alters the abundance trends (Kobayashi et al. 2020; Palla 2021).  \nIn Fig. 9 we show the evolution of [O / H] (upper panels) and [Si / H] (lower panels) as predicted by the model 3INF-2, compared with our restricted sample. The Figure highlights a global agreement between the model scenario and the data. [Si / H] gradient resembles the trends already described in this Section for Fe, sharing similar gradients slopes with the [Fe / H] gradient in di ff erent age ranges. This reflects in a rather flat [Si / Fe] gradient throghout the last Gyr of galactic evolution. On the other hand, the predicted and observed [O / H] gradient is much shallower, resulting in a positive [O / Fe] gradient. For the O gradient, we also note that in the 1-3 Gyr age bin (Fig. 9 upper left panel) all OC abundances lay in the upper end of the range of values allowed by the model, not showing any sign of particular spread as opposed to Fe and Si. However, it is worth noting that for O in general we have to rely on a much lower number of abundance data, which limit us in understanding the real motivation of this discrepancy. Nonetheless, the [O / H] gradient in the youngest age bin (Fig. 9 upper right panel) is reproduced by the model predictions, well in agreement with our proposed scenario.  \nIn Fig. 10 we also report the evolution of the OCs and the three-infall scenario in the [X / Fe] vs. [Fe / H] diagram in di ff erent radial regions at di ff erent ages. In particular we show the results for Si, as we can rely on a much larger sampling of OCs, for a more thorough comparison relative to O. In Fig. 10 upper panels we show the comparison for ages larger than 3 Gyr, i.e. before the onset of the third infall. Here, the model predictions including the e ff ects of migration and abundance uncertainties are shown in each radial region with colored areas. Moreover, to guide the eye, the genuine chemical tracks at 6, 8 and 10 kpc (in the left, central and right panel, respectively) are highlighted with color between 3 and 7 Gyr, i.e. the age range considered in these panels. Also, in the left and right panel tracks for 4 and 12 kpc in the age range considered are shown. Despite of the small data sampling in this age range, we can say that data are well within the predicted [Si / Fe] vs. [Fe / H] ranges considering  \nmigration and abundance uncertainties e ff ects. In Fig. 10 lower panel we show instead the predictions for ages below 3 Gyr, the ones interested by the late gas accretion. As for the upper panel, we highlight the parts of the tracks with ages < 3 Gyr with color. Also in this case, the agreement between predictions and observations is remarkable, with all the clusters falling within the range of values allowed by model predictions.', '4.3. Discussion': 'Throughout this work we investigated the behaviour in radial abundance gradients as traced by Ma23 Gaia -ESO OCs, which show a decrease in their metal content towards the youngest population. Within the framework of our models of chemical evolution, that also include the e ff ects of stellar radial migration (following Frankel et al. 2018, 2020), we observe that the observed metal impoverishment of young OCs should be mainly caused by a chemical dilution from a recent gas infall episode, which in turn is triggering a late enhanced SF activity. These results confirm the conclusion by Sp23, who found the necessity of very late ( ≲ 3 Gyr of age) gas accretion episode to explain the abundance ratios in the solar neighbourhood observed in Gaia DR3 RVS spectra, and extend it to a much larger range of Galactocentric radii.  \nHowever, our preferred scenario generally requires a milder metal dilution than in Sp23 in this late gas infall. In particular, our best scenario results to be more similar to the "weak dilution" scenario in Sp23, with smaller contribution of the late infall to the total mass budget and larger SFE relative to their best model. The di ff erence in the results can be explained by the di ff erent data set adopted in this work to calibrate our best model. In fact, in Sp23 the model was tuned to reproduce the abundances from young massive stars from Gaia RVS spectra (Gaia Collaboration et al. 2023a; Recio-Blanco et al. 2023) and cross matched targets in APOGEE DR17 (Abdurro\'uf et al. 2022, see their Fig. 14). Ma23 showed that such a class of objects may su ff er of biases in the determination of stellar abundances, with an evident decrease in their metallic content, and this in turn critically influence our view on [X / Fe] vs. [Fe / H] diagrams and abundance gradients. This is also shown in this work, where we see how massive giants with low log g can influence our view on the gradients of young stellar populations (see 4.1). As shown in Magrini et al.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 9. Radial and [O / H] (upper panels) and [Si / H] (lower panels) gradients at 1 < Age / Gyr < 3 (left panels) and Age < 1 Gyr (right panels) for the model 3INF-2, including stellar migration and average OC spread (see 3.4). Lines show the results for the [Fe / H] gradient as predicted by the model at Ages = 0.5 Gyr (cyan lines), 1.8 Gyr (green dashed lines) and 2.8 Gyr (green solid lines). Filled circles represent the restricted sample within OCs in Magrini et al. (2023).  \n<!-- image -->  \n(2023), this problem concerns spectral analysis in general and is present in various surveys (e.g. APOGEE, GALAH) regardless of the analysis method and spectral range. An e ff ort is needed to improve the model atmospheres of low gravity giant stars, and to include the e ff ect of the magnetic field to advance the analysis of these stars. Spectral analysis of these stars will have even more implications in the future because in the ELT era the brightest giant stars will give us detailed information about chemical abundances in distant galaxies (Roederer et al. 2024).  \nIn the current situation, the choice of imposing motivated cuts in the stellar parameters of OC members allows us to be confident of being less prone to abundance systematics. Moreover, we compare the results of our model with all the di ff erent diagnostics available, i.e. abundance gradients in di ff erent chemical abundances as well as age-metallicity relations and [X / Fe] vs. [Fe / H] abundance patterns in di ff erent Galactocentric regions, all of them showing a good agreement between data and models. We also verify that the proposed scenario agrees with the observed values at solar radius and gradients of di ff erent physical quantities, such as SFR and gas surface densities (see Appendix B). In this way, we ensure that our proposed scenario of chemical evolution is robust in the context of the MW disk formation.  \nNonetheless, the search for unbiased samples with abundances from high-resolution spectroscopy comes at expenses of  \nthe sample size. In fact, our results are still limited by the moderate sample size of OCs within our restricted sample ( ≃ 50). More data are certainly needed to better probe the three-infall scenario of chemical evolution, imposing more stringent limits and maybe revising the values of the physical parameters that have been adopted in this work. In particular, further sampling in age and in Galactocentric radii will be fundamental to pursue this goal. As for example, it is worth noting that for moderate to young ages we clearly lack of data at R ≳ 12 kpc. This crucially limits our ability to draw firm conclusions on the evolution of Galactic outer regions, whose evolutionary trends are only suggested following a framework similar to that found for the innermost disk. However, reliable tracers of gradient evolution as OCs are clearly missing from the outer galaxy (see, e.g. Cantat-Gaudin et al. 2020) The CCs have started to be observed with high statistical significance also in the outermost regions of the Galaxy (e.g. da Silva et al. 2023; Trentin et al. 2024), but these objects only inform us on the present-day situation. Therefore, accurate abundances and ages of the stellar populations of the outermost Galaxy are needed, and planned (e.g. 4MOST, de Jong et al. 2012; WEAVE, Jin et al. 2024; PLATO, Rauer et al. 2022 ) and proposed (e.g. WST, Mainieri et al. 2024) facilities will undoubtedly help us in this search.  \nFig. 10. [Si / Fe] vs. [Fe / H] evolution for the model 3INF-2 in di ff erent radial ranges, i.e. R < 7 kpc (left panels), 7 < R / kpc < 9 (central panels) and R > 9 kpc (right panels). The upper panels show results for ages > 3 Gyr, while the lower panels for ages < 3 Gyr. The shaded cyan, light green and magenta areas are the model prediction in a certain radial range taking into account the e ff ect of stellar migration and OC spread (see 3.4). The solid lines represent genuine chemical evolution tracks at 6, 8 and 10 kpc and are colored in the age range considered in the respective panel. Colored dashed lines are the same as colored solid lines but for the radii of 4 kpc (left panels) and 12 kpc (right panel). Colored filled circles represent the restricted sample within OCs in Magrini et al. (2023). Grey points are selected field stars from the Gaia -ESO survey.  \n<!-- image -->', '5. Summary and conclusions': 'In this paper we studied the evolution of radial abundance gradients in the Milky Way (MW) by taking advantage of the sample of open clusters (OCs) from the last data release of the Gaia -ESO survey (Randich et al. 2022; Magrini et al. 2023). OCs, in fact, are among the best tracers of the shape and time evolution of the radial metallicity gradient due to their ages and distances, which can be properly measured by isochrone fitting of the complete sequence, and to high-resolution spectroscopic observations that provide precise abundances. Gaia -ESO dedicated about 30% of its observing time to provide the largest sample of precise and homogeneous stellar parameters and abundances of member stars in Galactic OCs.  \nFrom a theoretical perspective, we started from the well tested revised two-infall model of chemical evolution (e.g. Palla et al. 2020b, see also Spitoni et al. 2019), which successfully reproduces data from high-resolution surveys, such as APOGEE (Ahumada et al. 2020), and extended the comparison to the newly proposed scenario of the three-infall model (Spitoni et al. 2023), proposed in the light of constraints given by Gaia star formation history (e.g. Ruiz-Lara et al. 2020) and abundance ratios (Gaia Collaboration et al. 2023a). For all the probed  \nscenarios, our models take into account the e ff ects of stellar radial migration, by including well tested prescription from the literature (Frankel et al. 2018, 2020), allowing us to directly see its impact on the late evolution of radial abundance gradients.  \nOur main considerations and conclusions are thus summarised as follows:  \n- 1. conservatively, we excluded stars whose spectral analysis can generate metallicity bias. This choice reduces the trends between stellar parameters and metallicity within the same cluster. This solution does not solve the problem of analysing massive and / or low-gravity giants, but allows us to circumvent it until progress is made on atmosphere models for these stars and including magnetic activity;\n- 2. despite of restricting the OC sample to stars which should not hide biases in spectroscopic analysis, we find a metallicity decrease between intermediate age (1-3 Gyr) and young ( < 1 Gyr) OCs. We show that the radial metallicity gradient as traced by young OCs is overestimated by the predictions of the two-infall scenario, even by accounting for the e ff ect of stellar migration in the model. We also check whether the adoption of di ff erent nucleosynthetic yields for massive stars and Type Ia SNe, i.e. the main contributors to Fe enrichment,  \nmay a ff ect this result. We find negligible di ff erences between the di ff erent runs, confirming the conclusion;  \n- 3. to explain the observed low metallic content in young clusters, we propose that a late gas accretion episode triggering a metal dilution should have taken place. This is in agreement with the proposed scenario of the three-infall model by Spitoni et al. (2023) for the solar vicinity, which explains the recent star formation history and abundance ratios as derived by Gaia satellite, as a consequence of recent gas infall episode triggering enhanced SF at recent ages. It is worth noting that in this work, for the first time we extend the threeinfall scenario to the whole MW disk;\n- 4. at variance with the best model presented in Spitoni et al. (2023), we invoke a milder metal dilution for this late gas infall episode. In particular, our best scenario requires smaller contribution of the late infall to the mass budget forging the lowα disk (factor ∼ 1.5-4 lower), while larger infall timescale ( τ 3 ≃ 1 Gyr instead of τ 3 ≃ 0 . 1 Gyr) and star formation e ffi ciency ( ν 3 ≃ (2 / 3) ν 2 instead of ν 3 ≃ (1 / 2) ν 2). The di ff erence in this results can be explained by the di ff erent data sample adopted in this work. In fact, our model are thought to reproduce a sample of OCs cleaned of stars subject to biases in chemical abundance determination (Magrini et al. 2023, see also point 1 of this Section), whose class is instead considered in Spitoni et al. (2023) (see Gaia Collaboration et al. 2023a; Recio-Blanco et al. 2023).  \nFurther data are definitely needed to probe the new threeinfall scenario of chemical evolution on the whole MW disk, imposing more stringent limits and maybe revising the values of the physical parameters that have been adopted in this work.  \nHowever, the constraints coming from high-resolution, unbiased chemical abundances, precise ages and star formation histories from multiple tracers (e.g. Isern 2019; Mor et al. 2019; Ruiz-Lara et al. 2020) allow us to consider the proposed model as a robust and viable scenario for the MW disk formation.  \nAcknowledgements. The authors thank the referee for the careful reading of the manuscript and the useful comments improving the paper content. MP acknowledges financial support from the project "LEGO - Reconstructing the building blocks of the Galaxy by chemical tagging" granted by the Italian MUR through contract PRIN2022LLP8TK\\_001. LM, ES, MF, and SR thank INAF for the support (Large Grant EPOCH) and MP, LM, and CVV for the MiniGrant Checs. LM, MF, and SR acknowledge financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union - NextGenerationEU- Project \'Cosmic POT\' Grant Assignment Decree No. 2022X4TM3H by the Italian Ministry of Ministry of University and Research (MUR). MM thanks the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID 279384907 - SFB 1245, the State of Hessen within the Research Cluster ELEMENTS (Project ID 500 / 10.006) for financial support. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany\'s Excellence Strategy EXC-2094 - 390783311. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programmes 188.B3002, 193.B-0936, and 197.B-1074. These data products have been processed by the Cambridge Astronomy Survey Unit (CASU) at the Institute of Astronomy, University of Cambridge, and by the FLAMES / UVES reduction team at INAF / Osservatorio Astrofisico di Arcetri. These data have been obtained from the Gaia -ESO Survey Data Archive, prepared and hosted by the Wide Field Astronomy Unit, Institute for Astronomy, University of Edinburgh, which is funded by the UK Science and Technology Facilities Council. This work makes use of results from the European Space Agency (ESA) space mission Gaia . Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https: // www.cosmos.esa.int / gaia. The Gaia archive website is https: // archives.esac.esa.int / gaia.', 'References': "Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35 Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3 Anders, F., Chiappini, C., Minchev, I., et al. 2017, A&A, 600, A70 Baba, J. & Kawata, D. 2020, MNRAS, 492, 4500 Balser, D. S., Rood, R. T., Bania, T. M., & Anderson, L. D. 2011, ApJ, 738, 27  \n- Bernard, E. J. 2017, in SF2A-2017: Proceedings of the Annual meeting of the  \nFrench Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo,  \nF. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di  \nBilitewski, T. & Schönrich, R. 2012, MNRAS, 426, 2266  \n- Bovy, J. 2015, ApJS, 216, 29\n- Bragaglia, A., Alfaro, E. J., Flaccomio, E., et al. 2022, A&A, 659, A200\n- Bragança, G. A., Daflon, S., Lanz, T., et al. 2019, A&A, 625, A120\n- Buck, T. 2020, MNRAS, 491, 5435\n- Buder, S., Lind, K., Ness, M. K., et al. 2019, A&A, 624, A19\n- Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150\n- Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1 Carbajo-Hijarrubia, J., Casamiquela, L., Carrera, R., et al. 2024, arXiv e-prints,\n- arXiv:2405.00110\n- Carraro, G., Geisler, D., Villanova, S., Frinchaboy, P. M., & Majewski, S. R. 2007, A&A, 476, 217\n- Casali, G., Magrini, L., Frasca, A., et al. 2020a, A&A, 643, A12\n- Casali, G., Magrini, L., Tognelli, E., et al. 2019, A&A, 629, A62\n- Casali, G., Spina, L., Magrini, L., et al. 2020b, A&A, 639, A127\n- Cavichia, O., Mollá, M., Costa, R. D. D., & Maciel, W. J. 2014, MNRAS, 437, 3688\n- Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765\n- Chiappini, C., Matteucci, F., & Romano, D. 2001, ApJ, 554, 1044\n- Colavitti, E., Cescutti, G., Matteucci, F., & Murante, G. 2009, A&A, 496, 429\n- da Silva, R., Crestani, J., Bono, G., et al. 2022, A&A, 661, A104  \nda Silva, R., D'Orazi, V., Palla, M., et al. 2023, A&A, 678, A195  \nDaflon, S. & Cunha, K. 2004, ApJ, 617, 1115  \n- Dame, T. M. 1993, in American Institute of Physics Conference Series, Vol. 278, Back to the Galaxy, ed. S. S. Holt & F. Verter, 267-278\n- de Jong, R. S., Bellido-Tirado, O., Chiappini, C., et al. 2012, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84460T\n- Di Teodoro, E. M. & Peek, J. E. G. 2021, ApJ, 923, 220  \nDonor, J., Frinchaboy, P. M., Cunha, K., et al. 2020, AJ, 159, 199  \n- Esteban, C., Fang, X., García-Rojas, J., & Toribio San Cipriano, L. 2017, MNRAS, 471, 987\n- Frankel, N., Rix, H.-W., Ting, Y.-S., Ness, M., & Hogg, D. W. 2018, ApJ, 865, 96\n- Frankel, N., Sanders, J., Ting, Y.-S., & Rix, H.-W. 2020, ApJ, 896, 15\n- Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1\n- Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1\n- Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2\n- Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., et al. 2023a, A&A, 674, A38  \nGaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023b, A&A, 674, A1 Genovali, K., Lemasle, B., da Silva, R., et al. 2015, A&A, 580, A17 Gilmore, G., Randich, S., Asplund, M., et al. 2012, The Messenger, 147, 25 Gilmore, G., Randich, S., Worley, C. C., et al. 2022, A&A, 666, A120 Green, D. A. 2014, in Supernova Environmental Impacts, ed. A. Ray & R. A.  \n- Grand, R. J. J., Bustamante, S., Gómez, F. A., et al. 2018, MNRAS, 474, 3629 McCray, Vol. 296, 188-196  \nGrisoni, V., Spitoni, E., & Matteucci, F. 2018, MNRAS, 481, 2570 Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2015, A&A, 578, A58 Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 616, A86 Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132 Henry, R. B. C., Kwitter, K. B., Jaskot, A. E., et al. 2010, ApJ, 724, 748 Hopkins, P. F., Gurvich, A. B., Shen, X., et al. 2023, MNRAS, 525, 2241 Hourihane, A., François, P., Worley, C. C., et al. 2023, A&A, 676, A129  \nIsern, J. 2019, ApJ, 878, L11 Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439 Jackson, R. J., Je ff ries, R. D., Wright, N. J., et al. 2022, MNRAS, 509, 1664 Jacobson, H. R., Friel, E. D., Jílková, L., et al. 2016, A&A, 591, A37 Je ff ries, R. D., Jackson, R. J., Wright, N. J., et al. 2023, MNRAS, 523, 802 Jin, S., Trager, S. C., Dalton, G. B., et al. 2024, MNRAS, 530, 2688 Kennicutt, Robert C., J. 1998, ApJ, 498, 541 Kobayashi, C., Karakas, A. I., & Umeda, H. 2011, MNRAS, 414, 3231 Kobayashi, C., Leung, S.-C., & Nomoto, K. 2020, ApJ, 895, 138 Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145 Kordopatis, G., Binney, J., Gilmore, G., et al. 2015, MNRAS, 447, 3526 Kovtyukh, V., Lemasle, B., Bono, G., et al. 2022, MNRAS, 510, 1894 Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2019, ApJ, 887, 80 Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545  \n- Lemasle, B., François, P., Bono, G., et al. 2007, A&A, 467, 283\n- Lemasle, B., François, P., Piersimoni, A., et al. 2008, A&A, 490, 613\n- Leung, S.-C. & Nomoto, K. 2018, ApJ, 861, 143\n- Leung, S.-C. & Nomoto, K. 2020, ApJ, 888, 80\n- Limongi, M. & Chie ffi , A. 2018, ApJS, 237, 13\n- Lin, L., Li, C., Du, C., et al. 2020, MNRAS, 499, 1406\n- Luck, R. E. & Lambert, D. L. 2011, AJ, 142, 136\n- Maciel, W. J., Costa, R. D. D., & Uchida, M. M. M. 2003, A&A, 397, 667\n- Mackereth, J. T., Crain, R. A., Schiavon, R. P., et al. 2018, MNRAS, 477, 5072\n- Magrini, L., Randich, S., Kordopatis, G., et al. 2017, A&A, 603, A2\n- Magrini, L., Randich, S., Zoccali, M., et al. 2010, A&A, 523, A11\n- Magrini, L., Smiljanic, R., Franciosini, E., et al. 2021a, A&A, 655, A23\n- Magrini, L., Spina, L., Randich, S., et al. 2018, A&A, 617, A106\n- Magrini, L., Vázquez, C. V., Casali, G., et al. 2022, Universe, 8, 64\n- Magrini, L., Vescovi, D., Casali, G., et al. 2021b, A&A, 646, L2\n- Magrini, L., Viscasillas Vázquez, C., Spina, L., et al. 2023, A&A, 669, A119\n- Mainieri, V., Anderson, R. I., Brinchmann, J., et al. 2024, arXiv e-prints, arXiv:2403.05398\n- Matteucci, F. & Francois, P. 1989, MNRAS, 239, 885\n- Matteucci, F. & Recchi, S. 2001, ApJ, 558, 351\n- McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13\n- Melioli, C., Brighenti, F., D'Ercole, A., & de Gouveia Dal Pino, E. M. 2008, MNRAS, 388, 573\n- Melioli, C., Brighenti, F., D'Ercole, A., & de Gouveia Dal Pino, E. M. 2009, MNRAS, 399, 1089\n- Méndez-Delgado, J. E., Amayo, A., Arellano-Córdova, K. Z., et al. 2022, MNRAS, 510, 4436\n- Miglio, A., Chiappini, C., Mackereth, J. T., et al. 2021, A&A, 645, A85\n- Minchev, I., Anders, F., Recio-Blanco, A., et al. 2018, MNRAS, 481, 1645\n- Minchev, I., Famaey, B., Combes, F., et al. 2011, A&A, 527, A147\n- Molero, M., Magrini, L., Matteucci, F., et al. 2023, MNRAS, 523, 2974\n- Mor, R., Robin, A. C., Figueras, F., Roca-Fàbrega, S., & Luri, X. 2019, A&A, 624, L1\n- Mott, A., Spitoni, E., & Matteucci, F. 2013, MNRAS, 435, 2918\n- Myers, N., Donor, J., Spoo, T., et al. 2022, AJ, 164, 85\n- Nakanishi, H. & Sofue, Y. 2003, PASJ, 55, 191\n- Nakanishi, H. & Sofue, Y. 2006, PASJ, 58, 847\n- Nepal, S., Chiappini, C., Guiglion, G., et al. 2024, A&A, 681, L8\n- Nissen, P. E., Christensen-Dalsgaard, J., Mosumgaard, J. R., et al. 2020, A&A, 640, A81  \nNoguchi, M. 2018, Nature, 559, 585  \nOverbeek, J. C., Friel, E. D., Donati, P., et al. 2017, A&A, 598, A68  \n- Palla, M. 2021, MNRAS, 503, 3216\n- Palla, M., Matteucci, F., Calura, F., & Longo, F. 2020a, ApJ, 889, 4\n- Palla, M., Matteucci, F., Spitoni, E., Vincenzo, F., & Grisoni, V. 2020b, MNRAS, 498, 1710\n- Palla, M., Santos-Peral, P., Recio-Blanco, A., & Matteucci, F. 2022, A&A, 663, A125  \nPinsonneault, M. H., Elsworth, Y., Epstein, C., et al. 2014, ApJS, 215, 19  \nPinsonneault, M. H., Elsworth, Y. P., Tayar, J., et al. 2018, ApJS, 239, 32  \nPortinari, L. & Chiosi, C. 2000, A&A, 355, 929  \n- Prantzos, N., Abia, C., Limongi, M., Chie ffi , A., & Cristallo, S. 2018, MNRAS, 476, 3432\n- Queiroz, A. B. A., Anders, F., Chiappini, C., et al. 2020, A&A, 638, A76 Rana, N. C. 1991, ARA&A, 29, 129\n- Randich, S., Gilmore, G., & Gaia-ESO Consortium. 2013, The Messenger, 154, 47\n- Randich, S., Gilmore, G., Magrini, L., et al. 2022, A&A, 666, A121\n- Randich, S., Pasquini, L., Franciosini, E., et al. 2020, A&A, 640, L1\n- Randich, S., Sestito, P., & Pallavicini, R. 2003, A&A, 399, 133\n- Rauer, H., Aerts, C., Deleuil, M., et al. 2022, in European Planetary Science Congress, EPSC2022-453\n- Recio-Blanco, A., de Laverny, P., Kordopatis, G., et al. 2014, A&A, 567, A5 Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2023, A&A, 674, A29 Ripepi, V., Chemin, L., Molinaro, R., et al. 2022, MNRAS, 512, 563\n- Roca-Fàbrega, S., Llorente de Andrés, F., Chavero, C., Cifuentes, C., & de la Reza, R. 2021, A&A, 656, A64\n- Roederer, I. U., Alvarado-Gómez, J. D., Allende Prieto, C., et al. 2024, Experimental Astronomy, 57, 17\n- Rojas-Arriagada, A., Recio-Blanco, A., de Laverny, P., et al. 2016, A&A, 586, A39\n- Romano, D., Chiappini, C., Matteucci, F., & Tosi, M. 2005, A&A, 430, 491\n- Romano, D., Karakas, A. I., Tosi, M., & Matteucci, F. 2010, A&A, 522, A32\n- Romano, D., Magrini, L., Randich, S., et al. 2021, A&A, 653, A72\n- Romano, D., Matteucci, F., Salucci, P., & Chiappini, C. 2000, ApJ, 539, 235\n- Romano, D., Matteucci, F., Zhang, Z.-Y., Ivison, R. J., & Ventura, P. 2019, MNRAS, 490, 2838\n- Ruiz-Lara, T., Gallart, C., Bernard, E. J., & Cassisi, S. 2020, Nature Astronomy, 4, 965\n- Sanders, J. L. & Binney, J. 2015, MNRAS, 449, 3479\n- Santos-Peral, P., Recio-Blanco, A., Kordopatis, G., Fernández-Alvar, E., & de Laverny, P. 2021, A&A, arXiv:2106.13677\n- Schönrich, R. & Binney, J. 2009, MNRAS, 399, 1145\n- Schönrich, R. & McMillan, P. J. 2017, MNRAS, 467, 1154\n- Spina, L., Magrini, L., & Cunha, K. 2022, Universe, 8, 87\n- Spina, L., Randich, S., Magrini, L., et al. 2017, A&A, 601, A70\n- Spina, L., Ting, Y. S., De Silva, G. M., et al. 2021, MNRAS, 503, 3279\n- Spitoni, E. & Matteucci, F. 2011, A&A, 531, A72\n- Spitoni, E., Matteucci, F., Recchi, S., Cescutti, G., & Pipino, A. 2009, A&A, 504, 87\n- Spitoni, E., Recio-Blanco, A., de Laverny, P., et al. 2023, A&A, 670, A109\n- Spitoni, E., Silva Aguirre, V., Matteucci, F., Calura, F., & Grisoni, V. 2019, A&A, 623, A60\n- Spitoni, E., Verma, K., Silva Aguirre, V., et al. 2021, A&A, 647, A73\n- Stahler, S. W. & Palla, F. 2005, The Formation of Stars\n- Stanghellini, L. & Haywood, M. 2010, ApJ, 714, 1096\n- Stanghellini, L. & Haywood, M. 2018, ApJ, 862, 45\n- Stonkut˙e, E., Koposov, S. E., Howes, L. M., et al. 2016, MNRAS, 460, 1131\n- Trentin, E., Catanzaro, G., Ripepi, V., et al. 2024, arXiv e-prints, arXiv:2404.17299\n- Van der Swaelmen, M., Viscasillas Vázquez, C., Cescutti, G., et al. 2023, A&A, 670, A129\n- Ventura, P., Dell'Agli, F., Lugaro, M., et al. 2020, A&A, 641, A103\n- Ventura, P., Di Criscienzo, M., Carini, R., & D'Antona, F. 2013, MNRAS, 431, 3642\n- Ventura, P., Karakas, A., Dell'Agli, F., García-Hernández, D. A., & GuzmanRamirez, L. 2018, MNRAS, 475, 2282\n- Vincenzo, F. & Kobayashi, C. 2020, MNRAS, 496, 80\n- Viscasillas Vázquez, C., Magrini, L., Casali, G., et al. 2022, A&A, 660, A135\n- Viscasillas Vázquez, C., Magrini, L., Spina, L., et al. 2023, A&A, 679, A122\n- Willett, E., Miglio, A., Mackereth, J. T., et al. 2023, MNRAS, 526, 2141\n- Wong, T., Blitz, L., & Bosma, A. 2004, ApJ, 605, 183\n- Yong, D., Carney, B. W., & Friel, E. D. 2012, AJ, 144, 95\n- Zhang, M., Xiang, M., Zhang, H.-W., et al. 2021, ApJ, 922, 145", 'Appendix A: Probing present-day gradients overestimation with different yield sets': 'As mentioned in 3.3 and 4.1, we perform additional model runs testing di ff erent massive stars and Type Ia SN yields than the reference ones used in this paper (Limongi & Chie ffi 2018 and Iwamoto et al. 1999, respectively).  \nIn the following Figures we show the results of this experiment, by reporting the gradient evolution (left panels) and the gradient in the young age bin considering the e ff ects of stellar migration and abundance uncertainties in the model (right panels) In particular, Fig. A.1 shows the predicted radial [Fe / H] gradient obtained by adopting the massive star yields from Kobayashi et al. (2006) instead of those from Limongi & Chie ffi (2018). Fig. A.2 instead reports the predicted radial [Fe / H] gradient obtained by adopting the Type Ia SN yields from an equal mixture of Leung & Nomoto (2018) (nearMCh progenitors) and Leung & Nomoto (2020) (subMCh progenitors) instead of those Iwamoto et al. (1999). Regarding the latter Figure, it is worth mentioning that no significant changes are obtained by employing di ff erent mixtures in Type Ia yields.', 'Appendix B: Gradients in physical quantities from the three-infall model': 'As mentioned in 4.3, we checked whether our three-infall scenario reproduces the constraints in other gradients than those in chemical abundances, i.e. the present-day gradients in the physical quantities.  \nResults of the comparison between our best model 3INF-2 and literature constraints are shown in Fig. B.1 and B.2. In Fig. B.1, we show the prediction of the model for the SFR surface density. In the left panel, we show the time evolution of this quantity at di ff erent radii, with the present-day SFR observed in the solar vicinity (see Prantzos et al. 2018) shown as reference. In the right panel, we show instead the comparison between the predicted present-day gradient and measurements from the literature (Rana 1991; Stahler & Palla 2005; Green 2014). In Fig. B.2, we display the same scheme as Fig. B.1, but for the gas surface density, where the reference value for the solar vicinity and the observed gradients are taken from Dame (1993); Nakanishi & Sofue (2003, 2006). For a detailed discussion on the adopted data sets, we refer to Palla et al. (2020b).', 'Appendix C: The restricted OCs sample': 'As described in 2.1.1, to avoid observational biases in the computation of the mean abundance of OCs, we build a restricted sample of member stars, where only stars with log g > 2.5 and ξ < 1.8 km -1 are considered to compute the mean cluster abundances.  \nIn Tab. C.1, we provide the average abundances obtained with this membership selection, as well as all the obtained cluster parameters (see also Viscasillas Vázquez et al. 2022).  \n<!-- image -->  \nFig. A.1. Predicted radial [Fe / H] gradient from the model 2INF using Kobayashi et al. (2006) instead of Limongi & Chie ffi (2018) yields for massive stars. Left panel: same as Fig. 1, but showing the data for the restricted OC sample. Right panel: same of Fig. 3, but showing only the results for the young age bin ( < 1 Gyr).  \n<!-- image -->  \n<!-- image -->  \nFig. A.2. Predicted radial [Fe / H] gradient from the model 2INF using Leung & Nomoto (2018, 2020) instead of Iwamoto et al. (1999) yields for Type Ia SNe. Left panel: same as Fig. 1, but showing the data for the restricted OC sample. Right panel: same of Fig. 3, but showing only the results for the young age bin ( < 1 Gyr).  \n<!-- image -->  \n<!-- image -->  \nFig. B.1. Left panel: SFR surface density time evolution at 4, 8, 12 and 16 kpc. The present-day observed value for the solar vicinity is taken from Prantzos et al. (2018). Right panel: present-day radial SFR surface density gradient. Data are from Rana (1991) (black points with errorbar), Stahler & Palla (2005) (black squares with errorbar), Green (2014) (black crosses).  \n<!-- image -->  \n<!-- image -->  \nFig. B.2. Left panel: gas surface density time evolution at 4, 8, 12 and 16 kpc. The present-day observed value for the solar vicinity is taken from Dame (1993); Nakanishi & Sofue (2003, 2006). Right panel: present-day radial gas surface density gradient. The dashed curve is the average between the Dame (1993) and Nakanishi & Sofue (2003, 2006) data sets. The grey shaded region represents the typical uncertainty at each radius (see Palla et al. 2020b for more details).  \n<!-- image -->  \nTable C.1. Excerpt from the list mean cluster abundances and cluster parameters from our adopted restricted sample (see 2.1.1).  \nNotes. For the orbital parameters, we use the G alpy code, with the axis-symmetric potential MWP otential 2014 (Bovy 2015). The complete Table will be available at the CDS.'}