The Gaia-ESO Survey: A new approach to chemically characterising young open clusters II. Abundances of the neutron-capture elements Cu, Sr, Y, Zr, Ba, La, and Ce

Young open clusters (t<200 Myr) have been observed to exhibit several peculiarities in their chemical compositions, from a slightly sub-solar iron content, super-solar abundances of some atomic species (e.g. ionised chromium), and atypical enhancements of [Ba/Fe], with values up to +0.7 dex. Regarding the behaviour of the other $s$-process elements like yttrium, zirconium, lanthanum, and cerium, there is general disagreement in the literature. In this work we expand upon our previous analysis of a sample of five young open clusters (IC2391, IC2602, IC4665, NGC2516, and NGC2547) and one star-forming region (NGC2264), with the aim of determining abundances of different neutron-capture elements, mainly CuI, SrI, SrII, YII, ZrII, BaII, LaII, and CeII. We analysed high-resolution, high signal-to-noise spectra of 23 solar-type stars observed within the \textit{Gaia}-ESO survey. We find that our clusters have solar [Cu/Fe] within the uncertainties, while we confirm the super-solar [Ba/Fe] values (from +0.22 to +0.64 dex). Our analysis also points to mildly enhanced [Y/Fe] values (from 0 and +0.3 dex). For the other $s$-process elements we find that [X/Fe] ratios are solar at all ages. It is not possible to reconcile the anomalous behaviour of Ba and Y at young ages with standard stellar yields and Galactic chemical evolution model predictions. Thus, we explore different possible scenarios related to the behaviour of spectral lines, from the sensitivity to the presence of magnetic fields to the first ionisation potential effect. We also investigate the possibility that they may arise from alterations of the structure of the stellar photosphere due to higher levels of activity in such young stars. We are still unable to explain these enhancements, but we suggest that other elements (i.e. La) might be more reliable tracer of the $s$-process at young ages and encourage further observations.


Introduction
Open clusters (OCs) are excellent tracers of the chemical properties of the Galactic disc and their time evolution. Thanks to dedicated spectroscopic surveys (e.g. the APO Galactic Evolution Experiment, APOGEE, Cunha et al. 2016, Donor et al. 2018; the Open Clusters Chemical Abundances from Spanish Observatories, OCCASO, Casamiquela et al. 2019; GALactic Archaelogy with HERMES, GALAH, Spina et al. 2021) we can analyse these systems with a large amount of data. In particular, within the Gaia-ESO Survey (Gilmore et al. 2012;Randich et al. 2013) almost 80 OCs with Based on observations collected with the FLAMES instrument at VLT/UT2 telescope (Paranal Observatory, ESO, Chile), for the  up to 100 members, spanning ages between a few million to several billion years, have been homogeneously analysed. However, different studies over the last 15 years seem to indicate that young stars within 500 pc share a slightly sub-solar metal content (with [Fe/H] between −0.05 and −0.10 dex), both in OCs, moving groups and associations (e.g. James et al. 2006;Santos et al. 2008;Biazzo et al. 2011;Spina et al. 2014aSpina et al. ,b, 2017. This is in contrast with what is expected from the standard Galactic chemical evolution (GCE) models that predict an enrichment of the interstellar medium of 0.10-0.15 dex over the last 4-5 Gyr (e.g. Minchev et al. 2013).
Another intriguing aspect of young OCs (YOCs, i.e. OCs with ages 200 Myr) is the behaviour of the elements mainly produced via the slow neutron-capture process (hereafter sprocess elements, Käppeler et al. 2011, and references therein). Early analytical models found that the solar system abundances of the whole s-process elements could be explained by the contribution of the weak, the main and the strong components. The weak component accounts for the formation of elements up to the atomic mass A∼90 (from Fe to Sr) and it takes place mostly in massive stars during convective He core and C shell burning phases (e.g. The et al. 2007, Pignatari et al. 2010Sukhbold et al. 2016;Limongi & Chieffi 2018a). Most of the copper, gallium, and germanium in the solar system is made by the weak sprocess in massive stars (Pignatari et al. 2010). In particular, copper was thought to be mostly made by thermonuclear supernovae since the s-process contribution was limited (Matteucci et al. 1993). However, thanks to a new generation of neutron-capture reaction rates the s-process production of copper in massive stars was revised (Heil et al. 2008). Therefore, present s-process calculations in massive stars means that the missing copper and the solar abundances can be explained (Bisterzo et al. 2005;Romano & Matteucci 2007;Pignatari et al. 2010).
Starting from the pioneering work of D' Orazi et al. (2009), it has been confirmed that the observed [Ba/Fe] ratios dramatically increase at decreasing ages, reaching values up to +0.60 dex in very young clusters like IC 2391 and IC 2602 (ages of ∼ 30 − 50 Myr). Conversely, older clusters with ages 1 Gyr exhibit solarscaled abundances. These extraordinary enhancements cannot be explained, neither with non-local thermodynamic equilibrium (NLTE) effects nor with stellar nucleosynthesis and GCE models (Travaglio et al. 1999;Busso et al. 2001). As discussed in D' Orazi et al. (2009), increasing the stellar yields by a factor ∼6 for AGB stars with masses of 1 − 1.5 M , a GCE model is able to reproduce the observed abundances up to 500 − 600 Myr, but not the measured massive overproduction of Ba in the last 50−100 Myr. The Ba overabundance has subsequently been confirmed by other studies (e.g. Yong et al. 2012;Jacobson & Friel 2013, among the others).
Interest has also moved toward the behaviour of other sprocess elements (Y, Zr, La, and Ce). At present the abundance evolution of these elements with respect to age is matter of debate. Maiorca et al. (2011) measured abundance ratios for these elements in a sample of 19 OCs, with ages from 0.7 to 8.4 Gyr, and they found a steep growth at younger ages. Similar conclusions have been reached by Magrini et al. (2018), who analysed a sample of 22 OCs with ages spanning from 0.1 to 7 Gyr. Recently, Frasca et al. (2019) studied the young cluster ASCC 123 (age ∼ 150 Myr) and found an overabundance of Sr, Y, and Zr, with values between +0.3 and +0.5 dex. On the other hand, other studies have confirmed that young clusters and local moving groups display first and second peak elements with different behaviour to Ba, in all cases showing solar values of Y, Zr, La, and Ce (e.g. D'Orazi et al. 2012;Yong et al. 2012;Jacobson & Friel 2013;D'Orazi et al. 2017). Reddy & Lambert (2015) analysed stars belonging to five local associations  and they found a large spread in [Ba/Fe] ratios, from +0.07 to +0.32 dex. Mishenina et al. (2015) again confirmed the trend of increasing Ba at decreasing ages from the analysis of giant stars in five OCs, together with solar-like abundances of La. From the stellar nucleosynthesis point of view the most puzzling signature to explain is not the intrinsic enrichment of Ba, but the production of Ba disentangled from La. For pure nuclear physics reasons, this cannot be done in s-process conditions (e.g. Käppeler et al. 2011). At the same time, observations of old metal-poor rprocess-rich stars have confirmed that the r-process co-produces Ba and La in similar amounts (e.g. Sneden et al. 2008, and references therein). In the light of these considerations, Mishenina et al. (2015) proposed the intermediate (i-) process as a possible explanation of the Ba enrichment in OCs.
The i-process was first introduced by Cowan & Rose (1977), and it is characterised by neutron densities that are intermediate between the s-process and the r-process, of the order of 10 14−16 neutrons cm −3 . Under these conditions Ba production is disentangled from La, and it is indeed possible to reproduce the high [Ba/La] ratios seen in stars hosted by YOCs (Bertolli et al. 2013;Denissenkov et al. 2021). Different types of stars have been proposed as possible stellar hosts of the i-process: post-AGB stars (Herwig et al. 2011) and low-mass AGB stars (e.g. Lugaro et al. 2015;Cristallo et al. 2016;Choplin et al. 2021), super-AGB stars (Jones et al. 2016), rapidly-accreting white dwarfs (e.g. Denissenkov et al. 2017;Côté et al. 2018;Denissenkov et al. 2019) and massive stars (Roederer et al. 2016;Clarkson et al. 2018;Banerjee et al. 2018). However, in the context of YOCs, which stellar site where the i-process has become so relevant only in the last ∼Gyr is still a mystery. From the analysis of solar-twin stars, Reddy & Lambert (2017) found a mild increase in La, Ce, Nd, and Sm with decreasing ages, while the trend for [Ba/Fe] is more evident and confirms all the previous findings. They also provided an important piece of evidence in trying to solve the so-called barium puzzle. These authors detected a positive correlation between the activity index of the stars and their [Ba/Fe] ratios. A similar trend between [Ba/H] and chromospheric and accretion diagnostics were also found in the Lupus star-forming region (SFR) by Biazzo et al. (2017).
In Baratella et al. (2020b) (hereafter Paper I) we demonstrated that the higher levels of stellar activity could affect the formation of spectral lines forming in the upper layers of the photosphere. The difference of the equivalent width (EW) of the strong Fe lines between a 30 Myr solar analogue and the Sun increases at decreasing optical depth (i.e. moving up in the photosphere). The direct consequence is that the microturbulence velocity (ξ) parameter 1 should be increased when derived by imposing that strong and weak Fe lines provide the same abundance. Values of the order of 2.0-2.5 km s −1 have been found for ξ in young solar-type, main-sequence stars. The net effect is an underestimation of [Fe/H], with the various abundance ratios [X/Fe] rescaling accordingly. The same result was confirmed by Yana Galarza et al. (2019) and Spina et al. (2020), who proposed that the magnetic intensification could be responsible for the observed patterns. In both studies intermediate-age stars (∼ 400 Myr) were analysed; in our case we are dealing with much younger stars and the effects of activity could be so intense that the solution of magnetic intensification is not sufficient.  In Paper I we used the Gaia-ESO available reduced spectra of solar-type stars belonging to five YOCs and one SFR, and developed a new spectroscopic approach to overcome the above-mentioned issues. Here, we expand the analysis of these stars to derive the abundances of the heavy elements Cu, Sr, Y, Zr, Ba, La, and Ce. The main goal of this new investigation is to shed light on the behaviour of the s-process domi-nated elements. To our knowledge, for the cluster NGC 2547 and the SFR NGC 2264 no previous studies focusing on the heavy element abundances have been published to date. The exceptions are IC 2391 and IC 2602, which will be used as calibrators; NGC 2516, for which only few elements have been investigated in the past; and IC 4665, recently analysed by Spina et al. (2021) within the GALAH DR3 (Buder et al. 2021).
The paper is organised as follows. In Sect. 2 we present the stellar sample analysed, adopting the procedure described in Sect. 3: line-list selection, computation of the optical depth of line formation, and measurement of the element abundance. Our results, along with a comparison with literature estimates, are reported in Sect. 4. We discuss our findings and their scientific implications in Sect. 5, while Sect. 6 contains our conclusions.

Data
In this work we analysed high-resolution, high signal-to-noise ratio (S/N) spectra of 23 solar-type dwarf stars, with spectral types from F9 to K1. The selected targets are five YOCs (IC 2391, IC 2602, IC 4665, NGC 2516, and NGC 2547 and the SFR NGC 2264. The stellar spectra were acquired with the 580setup (spectral coverage 4800 − 6800 Å) of the FLAMES-UVES spectrograph (nominal resolution R = 47 000; Pasquini et al. 2002). The data reduction was performed by the Gaia-ESO consortium (see Sacco et al. 2014). In Paper I we selected only targets with rotational velocities v sini < 20 km s −1 to avoid significant line blending, and with S/N per pixel > 50. We also analysed spectra of the Sun and the four (old and slow-rotating) Gaia benchmark stars (hereafter GBS), namely α Cen A, τ Cet, β Hyi, and 18 Sco, exploiting UVES spectra taken from Blanco-Cuaresma et al. (2014b). Out of 34 GBS, our selection was restricted only to those targets with atmospheric parameters similar to our stars (Jofré et al. 2015b(Jofré et al. , 2018. To analyse the largest set of spectral lines possible, we included in our analysis lines found in the bluer region of the spectra. It is well known that the majority of strong, clean, and isolated atomic lines for heavy-element abundance determination are located in the wavelength range 3800 − 4800 Å, which is not accessible to our spectral setup. For this reason we searched through the ESO archive for further observational datasets of the cluster stars. We find that only star 08440521−5253171 (IC 2391) and star 10440681−6359351 (IC 2602) have been observed with the UVES, FEROS (R ∼ 48 000, λ = 3500 − 9200 Å - Kaufer et al. 1999), or HARPS (R ∼ 115 000, λ = 3830 − 6900 Å - Mayor et al. 2003) spectrographs. In addition, we also re-analysed three stars of IC 2391 that had been previously published in D' Orazi et al. (2017), namely PMM 1142, PMM 665, and PMM 4362. These stars were observed with UVES (blue setup λλ=3900 Å) in the framework of the programme ID 082.C-0218 (PI Melo).

Analysis
In the following, we describe extensively the selected lines used in our analysis (Sect. 3.1), the procedure to compute the optical depths (Sect. 3.2), and the details of the abundance measurements (Sect. 3.3).

Selection of the spectral lines
We carried out a careful selection of spectral lines in the redder part of the spectra, searching in the official Gaia-ESO survey master line list (Heiter et al. 2020) for lines of Cu i, Sr i and Sr ii, Y ii, Zr i and Zr ii, Ba ii, La ii, and Ce ii. From this list we selected only lines with highly accurate measurements of the atomic data (g f _ f lag=Y) and those that were mostly unblended (syn f lag=Y or U).
The Gaia-ESO line list covers the wavelength range 4750 − 6850 Å for the region of the UVES-580 setting. For the bluer part of the spectrum, we included lines that have been extensively used and that have been proven to be reliable (see e.g. D' Orazi et al. 2017). Despite the large number of lines available for each atomic species, we selected only those that are moderately strong and not blended in the solar spectrum. Since our stellar sample includes stars with v sini up to 20 km s −1 and some of the spectra are noisy, most of the pre-selected lines are too broad and not measurable. We show two examples for lines 5274.23 Å of Ce ii and 6390.48 Å of La ii in the left and right panels of Fig. 1, respectively. In this figure the observed spectra of the Sun, τ Cet, and star 08440521−5253171 of IC 2391 (v sini =16.7km s −1 and S/N=260) are displayed. As can be seen, both lines are already relatively weak, though still usable, in the Sun, but they disappear at higher values of v sini.
In the following we report the details of the lines used in the analysis (all the line lists are available upon request).
Copper: In the solar spectrum only lines of neutral Cu were identified. For this element we only relied on the line at 5105.54 Å since the one at 5700.24 Å is heavily blended in the Sun, whereas the line at 5782.13 Å falls in the UVES wavelength gap. For this element we considered the isotopic solar mixture of 69% of 63 Cu and 31% of 65 Cu . Copper is affected by isotopic broadening and hyperfine structure (HFS), for which we adopted values from Kurucz (2011). According to Shi et al. (2014), the NLTE corrections for line 5105.54 Å are small in the Sun, being of the order of +0.02 dex.
Strontium: We measured lines 4607.33 Å of Sr i and 4215.52 Å of Sr ii. According to Bergemann et al. (2012) the Sr i line has a NLTE correction of +0.10 dex in dwarf stars with solar metallicity, while line 4215.52 Å has negligible NLTE corrections. We note that the line of Sr ii is very strong and is almost saturated. It is also blended with a nearby Fe i line at 4215.42 Å and with the CN molecular lines. Both features have been accounted for in the spectral synthesis.
Yttrium: For Y ii we selected the lines at 4398.01 Å, 4883.69 Å, 4900.12 Å, 5087.42 Å, 5289.82 Å, and 5728.89 Å. The line at 4900.12 Å is blended with a nearby Ti i line at 4899.91Å, which becomes significant for v sini > 4 − 5 km s −1 . The lines at 5289.82 Å and 5728.89 Å are instead very weak in the Sun, both having EW∼ 4 mÅ, and thus we are not able to measure them in our targets. The HFS for Y can be neglected, as previously discussed in several papers (e.g. D' Orazi et al. 2017).
Zirconium: Zirconium is present in the form of neutral and ionised species in the solar photosphere. However, the available reliable lines of Zr i (at 6127.44 Å, 6134.55 Å, 6140.46 Å, 6143.2 Å, and 6445.74 Å ) are too weak to be measured in our sample stars. In the bluer region we used lines 4050.32 Å and 4208.98 Å of Zr ii. The line at 5112.27 Å of Zr ii is also weak, and we were able to measure it only in the Sun, in the GBS sample, and in one target. According to Velichko et al. (2010), Zr ii lines form under LTE conditions in solar-type stars.
Barium: For ionised barium, instead, we used only 5853.7 Å, which is not blended and does not experience severe HFS or isotopic shifts. To our knowledge this line is the best diagnostic to measure the Ba abundance. There are other lines of Ba ii in our spectral range. However, the resonance Ba ii line at 4554.03 Å is almost saturated; the line at 6141.7 Å is known to be blended Article number, page 4 of 26 M. Baratella et al.: The Gaia-ESO Survey: A new approach to chemically characterising young open clusters with a strong Fe i line; the line at 6496.9 Å is also blended with an iron line, and it is affected by strong NLTE effects. Reddy & Lambert (2015) explored the possible detection of the line of neutral Ba at 5535 Å in a sample of F-G dwarfs. Even so, this line is blended with a strong Fe i line that dominates the profile and its abundance shows a significant correlation with T eff (see their Fig. 7), most likely caused by large NLTE effects. Therefore, as already pointed out by the authors, in the absence of NLTE corrections it is not suitable to derive accurate abundances and to solve the Ba puzzle. Nevertheless, to obtain more accurate abundances, we also considered the HFS data from McWilliam (1998) and we adopted the isotopic solar mixture of 81% for ( 134 Ba + 136 Ba + 138 Ba) and 19% for ( 135 Ba + 137 Ba) (see Grevesse et al. 2015 for further details). According to Korotin et al. (2015), the NLTE corrections are small for stars in the parameter space covered by our sample. Gallagher et al. (2020) derived NLTE corrections for the Sun and found that the ∆ 1DNLTE =−0.11 dex and ∆ 3DNLTE =0.03 dex. However, there are no available tables of NLTE corrections for stars with parameters similar to our sample. We note in this context that the NLTE corrections are not sufficient to solve the Ba puzzle. For these reasons, we report the LTE Ba abundances.
Lanthanum: For La ii we selected lines at 4804.04 Å, 4920.98 Å, 5122.99 Å, and 6390.48 Å. Unfortunately, none of them is strong enough to be measured in our stars in the range 4800 − 6800 Å. Instead, in the bluer part we relied on the measurements of lines 3988.51 Å and 4086.71 Å. Lanthanum has one single isotope 139 La that accounts for 99.9% of the total La abundance in the solar material, and it is strongly affected by HFS. We followed the prescriptions by Lawler et al. (2001).
In Table 1 only the lines for which we obtained more than one measurement in our stellar sample are indicated; the element (Column 1), the corresponding wavelength (Column 2), the excitation potential energy (E.P., Column 3), the oscillator strength log g f (Column 4), references for the log g f values (Column 5), and the solar log(X) (Column 6) of each individual line are given (the average solar abundances are in Table 3). For each line we computed the Landé factor g L following Landi Degl'Innocenti (1982) (Column 7). Instead, the first ionisation potential (FIP) values are taken from

Computation of the optical depths of line formation
Our working hypothesis is that lines forming in the upper layers of the photosphere are more influenced by the higher levels of the activity present in young stars. Therefore, these lines are stronger than those observed in the spectra of old and quiet stars, affecting the derivation of the stellar parameters and, finally, the abundances (Yana Galarza et al. 2019;Baratella et al. 2020b;Spina et al. 2020). We computed the optical depth of line formation log τ 5000 of all the selected lines in a consistent way following the prescriptions by Gurtovenko & Sheminova (2015).
Calculations of the average formation depth of the absorption line are based on the contribution function (CF), which describes the contribution of the various layers of the stellar atmosphere to the absorption line (or line depression). Gurtovenko et al. (1974) suggested to use the CF as the integrand of the emergent line depression in the solar disc centre, computed as where R=1 − I l /I c is the line l depression and η is the ratio of the coefficient of the selective absorption to the coefficient of continuum c absorption. In the same formula, g is the Unsold weighting function for LTE (Unsöld 1932), multiplied by the emergent intensity in the continuum I c (τ c = 0). This weighting function is expressed as where B(τ c ) is the Plank function.
When we are dealing with the interpretation of observed line profiles or line depth in its centre or equivalent width, we use the average depth of the layers contributing to the absorption line. The average depth at a given wavelength position of the line profile ∆λ and at given position on the stellar disc µ = cos θ is calculated by the following formula: ( 3) If we consider the integrated line profile (i.e. the EW), its average optical depth is calculated as Here R(∆λ, µ) is the line depression at ∆λ and µ, while λ 1 and λ 2 are the initial and final wavelength positions of the line profile, respectively. To obtain the average depth of formation of the line depression observed in the stellar spectra at a given ∆λ, we use the following formula: Instead, to have the average formation depth of the whole line profile using EWs, we use where R is the line depression in the spectra of stellar flux. In this work we computed the average depth of line formation of our lines both in the core (log τ 5000 core ) and in the whole profile (log τ 5000 full ), reported in Column 7 and 8 of Table 2, respectively. For this calculation we assumed the LTE approximation and we considered the damping constant associated with the van der Waals force between the absorbing and perturbing atoms to be equal to γ 6 according to the classical Unsold approximation. Our assumptions are acceptable to estimate the average Table 2: Optical depths of line formation log τ 5000 core of the line core and log τ 5000 full of the full line profile. depth of the formation of weak, moderate, and moderately strong lines, like those analysed here. We measured the EW obs and R (Column 5 and 6, respectively) of each line using the ARESv.2 software (Sousa et al. 2015) and performed the synthesis in the solar spectrum with the SPANSAT code of Gadun & Sheminova (1988). The log τ 5000 value was then derived from the abundance obtained when the EW of the synthetic line matches the EW obs . We adopted the MARCS solar atmosphere model with the chemical composition taken from Lodders (2019), and with the stellar parameters reported in Table A.1. We also considered the macroturbulent velocity equal to 2 km/s and v sini = 1.84 km s −1 (Sheminova 2019). We considered the HFS indirectly in our computations for those lines labelled y in Column 9 of Table 2. Lines affected by strong HFS are split into multiple components, resulting in larger EWs. Thus, completely neglecting the HFS when deriving the abundances from EWs will result in overestimated values Jofré et al. 2017), and all lines form in higher layers of the photosphere since the overestimated abundances (corresponding to the larger EW obs ) are used in the computation of the optical depth. At the same time, the abundance calculated from the fitting of central depths (R ) of the lines can be underestimated and the depths of the formation of the lines will be large (i.e. the line will form in deeper layers). Therefore, these values are indicative and should be considered with some caution. We believe that the possible errors due to the adopted approximations in the above-mentioned computations are small and that they can be neglected within the limits of specific calculations of the depths of line formation.

Abundance measurements
The abundances of the s-process elements were derived using the technique of synthetic spectrum line profile fitting through the driver synth in MOOG (version 2017(version , Sneden 1973Sobeck et al. 2011). We used 1D-LTE plane-parallel MARCS model atmospheres (Gustafsson et al. 2008), fixing the atmospheric parameters to the values we found in Paper I, both for the cluster stars and the GBS. We report all the stellar parameters in Table A.1 for completeness. All abundances were computed in a strictly differential way (i.e. line-by-line) with respect to the Sun as [X/H] =log(X) − log(X) (using the individual abundances log(X) in column 6 of Table 1 (2019), corresponding to the Gaia G BP pass-band (their Table 6). The rotational broadening profiles were calculated using the v sini measured by the Osservatorio Astrofisico di Catania (OACT) Node of the Gaia-ESO consortium, measured using the routine ROTFIT (see e.g. Frasca et al. 2015, and references therein, for more details); the values are reported in Column 4 of Table 5.

Error budget
There are two sources of internal uncertainties affecting the [X/Fe] ratios derived via spectral synthesis. The first kind of error, σ 1 , is related to the best fit procedure, and spans values from ±0.06 to ±0.3 dex depending on the quality of the spectra, mainly the S/N, which affects the continuum placement, and on the individual spectral features under consideration.
The second kind of error, σ 2 , is related to uncertainties in the stellar parameters (Table A. 1). We calculated these uncertainties in a conservative way by varying each quantity separately, leaving the others unchanged, and evaluating the abundance sensitivity to those changes as We report both errors in Tables 1, 4, 5, and 6.

The Sun and the Gaia benchmarks
In Table 3 we report the average solar abundance values and we compare them with the photospheric abundances from Grevesse et al. (2015), the meteoritic abundances from Lodders (2019), and the results reported in the Gaia-ESO internal Data Release 4 (for La and Sr) and 5 (iDR4 and iDR5, respectively). For our values of Y, Zr, La, and Ce in Table 3 we take the simple mean of abundances derived in the bluer and in the redder regions; the uncertainties are computed as the errors on the mean. For Sr the value in the Table 3 is the average between the Sr i (corrected for NLTE) and the Sr ii results. For Ba and Cu, for which we analysed only one line, we report the individual abundance; the uncertainty is the error on the fitting procedure (σ 1 ). As can be seen from the comparison with the literature values, our mean solar abundances agree well with the photospheric values from Grevesse et al. (2015), and with the meteoritic abundances from Lodders (2019). They are also in fair agreement with the Gaia-ESO iDR5 and iDR4 results. For the GBS we list in Table 4 the abundance ratios [X/Fe] obtained for each line separately in the blue and red spectral ranges, while the average values and the comparison with the literature values are provided in Table 7 (the GBS atmospheric parameters can be found in the Appendix; see Table A.1). For Cu and Ba, for which we have only one line, the uncertainties in the table represent the total uncertainties (computed as the square root of the sum of the squares of σ 1 and σ 2 ). Instead, for the other elements the error is the standard deviation of the abundances obtained from the different lines. As can be seen, for the GBS that are old and for quiet stars we obtain solar-scaled abundances for all the heavy elements. Our estimates for α Cen A, τ Cet, and 18 Sco are in fair agreement with Luck (2018), who analysed high-resolution (R= 115 000), high S/N HARPS spectra of a sample of 907 F-G-K dwarf stars in the solar neighbourhood. Our estimates are also similar to the results by Casali et al. (2020), who performed a detailed spectral analysis of HARPS spectra of 560 solar-type stars. The marginal discrepancies of the individual abundances could be related to the different line lists (inclusion of HFS and isotopic splitting) and methods used (i.e. EW versus spectral synthesis). We also find that our values are in good agreement with Casamiquela et al. (2020), who analysed a sample of GBS. They used the same spectra and Gaia-ESO line list as we did, but they analysed different lines (e.g. the Cu i lines at 5218.20 Å and at 5220.07 Å). Thus, the small discrepancies that can be seen in Table 7 (for example for [Cu/Fe] of α Cen A or for [La/Fe] of τ Cet) can be related to the different lines used. β Hyi is a slightly evolved GBS, which is in the sub-giant branch phase of its evolution. In the literature we find few measurements of the s-process elements: our [Ba/Fe] estimate is in fair agreement with Bensby et al. (2014) and Jofré et al. (2015a). To our knowledge these are the only results of the heavy elements abundances for β Hyi.

The young clusters
As already mentioned in Sect.2, in addition to the stars observed within the GES, we analysed three stars of IC 2391 taken from D' Orazi et al. (2017). Firstly, we re-determined the atmospheric parameters by applying our new spectroscopic approach using titanium lines, as described in Paper I. Our final values can be found in Table 8. The largest differences are seen for the ξ parameter, for which we obtained a difference ∆ξ(our-DS13)=−0.59 km s −1 for star PMM 665, −0.43 km s −1 for PMM 4362, and −0.55 km s −1 for star PMM 1142. Instead, for star PMM 1142 our T eff value is 300 K lower than the value found by DS13, mostly due to the different line lists used. However, our spectroscopic T eff is corroborated by the estimates derived using 2MASS photometry (Cutri et al. 2003) and using the relations by Casagrande et al. (2010), which are equal to T (J − K) = 5378 ± 105 K, T (V − J) = 5386 ± 135 K, T (V − H) = 5408 ± 153 K, and T (V − K) = 5380 ± 155 K. An increase of ∼ 0.6 km s −1 in the ξ values results in a decrease of ∼ 0.2 dex in [Fe/H] and of ∼ 0.3 dex in the Y abundance, which is derived from moderately strong lines. Instead, the same variation produces a negligible change in La abundance, for which only weak lines have been used. Thus, the sensitivity to variations in the ξ parameter depends mainly on the strength of the line.
The final abundances of the individual stars can be found in Table 5 for the red spectral setup and Table 6 for the blue range. The uncertainties σ 1 and σ 2 indicated in the two tables are calculated as described in Sect. 3.4. The mean cluster abundances for each elements are given in Table 7; in this case, the uncertainties were computed as the error on the mean. We note that these mean values come from the averaged abundances of redder and bluer spectral ranges. For NGC 2264 we analysed only one star; therefore, we assumed as a conservative error value the uncertainty in the fitting procedure. For completeness, we also report the v sini measured by us along with values derived by the OACT node of the Gaia-ESO consortium (v sini lit. ), as well as the S/N. The mean difference (and error on the mean) between our v sini and the OACT values is equal to 0.6 ± 0.2 km s −1 , with Notes. The first source of uncertainty is due to the best fit procedure, while the second is related to uncertainties in the stellar parameters (details on the computations can be found in Sect. 3.4).
a standard deviation σ=0.8 km s −1 . Overall, our measurements are in good agreement with the OACT node results. As can be seen in Table 7 derived a log(Ba) equal to 2.22 dex, which is 0.09 dex lower than our adopted value. For IC 2391 we find a value of [Ba/Fe] that is 0.09 dex lower than that found by De Silva et al. (2013), and that is in fair agreement within the observational uncertainties.
When focusing on the other s-process elements, we detect a mild enhancement for [Y/Fe], at ∼ 0.3 dex level, higher with respect to the values reported by D'Orazi et al. (2017) Maiorca et al. (2011) and Magrini et al. (2018) reported +0.18 ± 0.02 and +0.38 ± 0.10, respectively. However, we measured Ce only in one star, 07544342−6024437. The discrepancy between our results and the literature values could be due to the different techniques used.
The cluster IC 4665 was recently used by Spina et al. (2021) in his study of Galactic OCs observed within the GALAH survey (Buder et al. 2021) (for which we considered only the results obtained with the dwarf stars). As can be seen in Table 7, there is a large discrepancy between our [Cu/Fe] value and that of Spina et al. (2021), possibly due to the different line lists, techniques, stars per clusters analysed, and the lower resolution of the spectra analysed, all of which can affect the measurements in crowded regions. The measured [Y/Fe] ratios in the two studies agree very well, and we also found fairly good agreement for [Ba/Fe]. In the table we also report the value of [La/Fe], for which caution should be used, however, as stated in Buder et al. (2021). These measurements could be affected by heavy blends of La lines used in the analysis in the GALAH survey.
Finally, for NGC 2547 and the SFR NGC 2264, this is the first time the abundances of the s-process elements have been determined. We

Trends with stellar parameters
In Figs. 2, 3, and 4, we plot the [X/Fe] ratios as a function of the stellar parameters T eff , log g, and v sini, respectively. In these plots the points are colour-coded according to age. There are no  Values reported in both the iDR5 and iDR6, no error provided because it is an upper limit (VSINI FLAG=1).
Article number, page 9 of 26 A&A proofs: manuscript no. baratella   Notes. The uncertainties for the GBS stars is the standard deviation of different lines (for Cu and Ba, for which we only analysed one line, the error is the quadratic sum of σ 1 and σ 2 ). The uncertainties of the mean abundances computed for our clusters are the errors on the mean, representing the star-to-star variation. *: Since we analysed only one star in this SFR, we report the individual abundance values and the uncertainties on the fitting procedure. **: Measurement based on only one target. ***: Values affected by strong blends (Buder et al. 2021 For both elements (Y and Ba) a sharp separation between the blue dots (the youngest stars, with ages of less than 50 Myr) and the red dots (the oldest stars in our sample, with ages of ∼150 Myr) can be seen. For Cu, instead, the results are homogeneously distributed with the age.

Discussion
As already mentioned in Sect. 1, it is not possible to reconcile the super-solar abundances of Ba together with solar La and Ce abundances with the predictions of the s-process and r-process nucleosynthesis models (e.g. Kobayashi et al. 2020 and references therein), without invoking other processes. On the other hand, the enrichment of Y with respect to Sr and Zr is less clear and needs to be studied in more detail. As discussed in this work, the Y enhancement could be caused by some observational is- PMM 1142 5400 ± 100 4.28 ± 0.07 0.95 ± 0.10 0.00 ± 0.02 ± 0.07 +0.11 ± 0.03 ± 0.06 +0.04 ± 0.02 ± 0.10 +0.06 ± 0.03 ± 0.05 PMM 665 5425 ± 100 4.47 ± 0.05 1.01 ± 0.15 +0.08 ± 0.02 ± 0.07 +0.14 ± 0.05 ± 0.06 +0.10 ± 0.02 ± 0.11 +0.11 ± 0.05 ± 0.05 PMM 4362 5550 ± 100 4.35 ± 0.05 0.97 ± 0.15 +0.13 ± 0.02 ± 0.08 +0.14 ± 0.04 ± 0.06 +0.09 ± 0.03 ± 0.11 +0.12 ± 0.02 ± 0.05 Notes. Atmospheric parameters and abundances of Fe and Ti derived with the same methodology as in Paper i of the stars from D' Orazi et al. (2017). The uncertainties are due to the dispersion among different atomic lines (first value) and to uncertainties on the stellar parameters (second value). sue. Moreover, a large variety of processes can contribute to the production of elements in the Sr-Y-Zr region, which could have caused the variations observed in YOCs. Together with the nucleosynthesis processes mentioned in the previous sections (the s-process, the i-process, and the r-process), elements in this mass region may also be made in different neutrino-wind components from CCSNe (e.g. Farouqi et al. 2009;Roberts et al. 2010; Arcones & Montes 2011) and in electron-capture supernovae (e.g. Wanajo et al. 2011;Jones et al. 2019). From the observational point of view, the difficulty may be related to the spectral analysis and mechanisms that are at work on the photosphere of a young star to magnify the Ba abundance or to the nucleosynthesis processes that produce the elements heavier than Fe in the Galaxy. We discuss both possibilities below; in Sect. 5.1 the problems that may be related to the spectra lines and in Sec. 5.2 the potential issues with GCE models.

Behaviour of spectral lines
Since we cannot reconcile the observed Ba overabundance and the slightly super-solar values of Y with nucleosynthesis models predictions, we believe that the key to understanding the Ba puzzle might rely on the age of the stars. The younger the star, the higher the levels of its activity, at the level of the chromosphere or of photospheric magnetic fields. This in turn can result in an alteration of the photospheric structure, and consequently in an alteration of the profile (i.e. strengths) of the spectral lines. Thus, in young stars it is important to know where the line forms in the photosphere and how it is affected by magnetic activity, as already demonstrated by Spina et al. (2020).
Looking at Table 2, where the optical depths of line formation are listed, we infer that the line at 5105.54 Å of Cu and line at 5853.7 Å of Ba form at similar depths, with log τ 5000 core of −3.4 and −3.2, respectively. Thus, we would expect to observe the same effects in the two elemental abundances derived through these lines. Nonetheless, we obtain solar values of [Cu/Fe], while [Ba/Fe] is enhanced (between +0.25 and +0.70 dex). This is also confirmed by Fig. 5, where we compare the spectrum of the Sun (light pink line) with the solar analogue 10442256−6415301 (black line) that belongs to IC 2602 (age ∼ 35 Myr). In these plots the solar spectrum was convolved with a rotational profile with v sini =11 km s −1 , to match the v sini of the young star, and we added Gaussian noise to obtain S/N=110 (we use the iSpec tool by Blanco-Cuaresma et al. It has been observed that strong ionised lines of Fe, Ti, and Cr yield large abundances in young and cool (T eff <5400 K) stars (see e.g. D'Orazi et al. 2009;Schuler et al. 2010;Tsantaki et al. 2019;Baratella et al. 2020a). This is referred to as the ionisation balance problem. Aleo et al. (2017) explain this effect as being due to the presence of undetected blends in the ionised lines (in particular of Fe) that become more severe in the cool regime. We selected a set of seven strong Fe ii lines that were initially excluded from our analysis in Paper i. For each line in Fig. A.1 we compare the spectra of the Sun and star 10442256−46415301, as in Fig. 5. Most of the Fe ii lines are deeper in the young star than in the Sun, as observed for the Ba ii line. This is also corroborated by the measured EWs and abundances obtained from the Fe ii lines (computed by adopting the stellar parameters of Paper i), as can be seen in Table 9.
However, we note that even lines of neutral species show a behaviour similar to that of the ionised lines. In Fig. A.1 it is evident that Fe i (in the first and second panels of the figure) is stronger in the young star than in the Sun. In the bottom panel of Fig. 5 the blend Ca i+Ni i at 5857.5 Å behaves similarly to Fe i. We then compare the two spectra in small windows around two Mg i lines in Fig. A.3, and eight Ca i lines in Fig. A Ca we note that the weak lines (panels on the left) have similar depths, while the strong lines (panels on the right) are deeper in the young star (black line) than in the Sun (light pink line). Barium and strontium belong to the same group in the periodic table, so their outermost electron shells have similar configurations. Furthermore, they share similar nucleosynthetic origins. Hence, we should witness some effects also on Sr abundances, for which we exploit lines at 4607.33 Å of Sr i (∆ NLTE =+0.1 dex) and 4215.52 Å of Sr ii. As can be seen in Table 6, we find good agreement between the Sr i and Sr ii abundances, at least in stars 0844052−5253171 (IC 2391) and 10440681−6359351 (IC 2602). In general, we find that [Sr/Fe] is solar in both cases. Looking at the depth of formation (Table  2), the Sr i line forms deeper in the photopshere than the Ba line. This might in principle explain why we obtain solar-scaled abundances. Conversely, the Sr ii line forms at log τ 5000 =−5.2, in the upper layers. However, when deriving the abundance of Sr ii, we obtain again a solar composition, as can be seen from Fig. A.2 where we plot the best fit models of Sr ii lines in the Sun, star 08440521-5253171 (IC 2391), and τ Cet.
Regarding the Y lines at 4883.69 Å and at 5087.42 Å, we found that they form at a similar depth, with log τ 5000 core equal to −2.6 and −2.1, respectively. From the comparison of the solar and the young solar analogue spectra in Fig. 5 (central panel), it can be seen that the Y ii line at 4883.69 Å is stronger in the young star than in the Sun. This is in agreement with our derived abundance estimates. Interestingly, the La ii lines analysed in this work should exhibit a behaviour similar to the Y ii features (they share formation depths, ionisation stages, and nucleosynthesis channels). Nevertheless, the La abundances are solar, whereas [Y/Fe] values are enhanced.
Finally, in Fig. A.5 we compare Sc i and Sc ii lines in the same way as in Fig. 5. Scandium has the same electronic configuration as Y, and is similar to La, so we expect these elements to show similar behaviours. Instead, from Fig. A.5 it is evident how the profile of both neutral and ionised Sc lines are the same in the Sun and in the young star. In this case there are no differences between the two spectra, for any line, with the exception of 5658.361 Å. The small difference we see in this line may be due to a blend with a nearby Fe i line.
We conclude that line formation depth and ionisation stages of the elements are not able to fully explain the very peculiar pattern of s-process elements in young open clusters. On the other hand, we cannot exclude that different (conspiring) mechanisms could be simultaneously at work. Figure 6 shows that there seems to be a correlation between the larger abundances of Y and Ba with decreasing log R HK (i.e. higher levels of stellar activity). We considered the log R HK values already computed in Paper I, which were derived from the X-ray luminosities found in the literature and using the conversion relation by Mamajek & Hillenbrand (2008). We note that the X-ray luminosities (and hence the log R HK values) are not synchronous to our spectra, and consequently, to the derived abundances. Nevertheless, globally, we can conclude that indeed there is an indication of a correlation between the [Y/Fe] and [Ba/Fe] and log R HK (not so evident for [Cu/Fe]). This plot surely deserves an in-depth investigation, and further observations are needed to study the behaviour at log R HK > −4.0 and log R HK < −4.4. Spina et al. (2020) proposed magnetic intensification as a possible explanation of the anomalous abundances of Ba (and of other elements). In Biazzo et al. (2017), this possibility is also explored. Given the young age of our targets, this seems to be a promising solution. The presence of magnetic fields causes the atomic levels to be split into different components according to the Zeeman effect. This results in a broadening of the spectral line, with increased EWs and reduced line depths, which however is not seen in our lines (see Fig.5). The amount of splitting is directly proportional to g L , the square of the wavelength and the magnitude B of the magnetic field. To evaluate the sensitivity to magnetic fields of each line, we compute the Landé factor g L (Column 7 of Table 1), as described in Sect. 3.1. Overall, our lines have g L < 1.3, which is relatively low. Looking at the Ba and Cu lines, we note that they have very similar g L values, 1.07 and 1.10, respectively. Nevertheless, we obtained super-solar Ba abundance and solar Cu. The La line at 3988.5 Å has g L = 1.33, which is the highest value, but we obtain solar abundance. Thus, we conclude that the magnetic intensification cannot explain the high values we obtain.
Another possibility is the first ionisation potential (FIP) effect; our lines, in particular Ba, seem to be suitable candidates to show this. It has been shown that the coronal abundances derived from lines with FIP values below 10 eV in the Sun are enhanced with respect to the photospheric values (see e.g. the review by Laming 2015). Sheminova & Solanki (1999) explored the idea that the gas exhibiting the FIP effect in the corona is connected to the photosphere through magnetic flux tubes, generated from magnetic elements or sunspots present on the surface. In principle, the same enhancement observed in the corona could also be found for photospheric abundances. We retrieve the FIP values (Column 8 of Table 1) from Gray (1992). As can be seen, Ba and La lines have similar FIP values; therefore, this does not explain their discrepant abundances. However, the higher levels of activity due to the very young ages of our stars could be completely different than what is observed in other active, older stars.
In summary, all the possible effects described above may play a role; however, there is no convincing evidence that any of them provides a definitive solution, yet.

The Galactic chemical evolution of s-process elements at young ages
In Fig. 7 we plot the Ba, Y, La, and Cu abundance ratios as a function of the age of the Galactic OCs taken from different sources in the literature. For all the clusters in the different samples we considered the ages from Cantat-Gaudin et al. (2020), who report that the uncertainties in log t for young clus-ters ranges from 0.15 to 0.25, while for old clusters from 0.1 to 0.2.
The [Ba/Fe] time evolution, with increasing values at decreasing ages, is confirmed by the observational data and we can now also detect a significant scatter at young ages. In particular, for the SFR NGC 2264 (age ≈ 5 Myr), the Lupus region (age ∼3 Myr), and the Orion subgroup Ic (age ∼3 Myr) the values are [Ba/Fe]≈+0.4 dex (this study), [Ba/Fe]≈0.7 dex (Biazzo et al. 2017), and [Ba/Fe]≈0.1 dex (Reddy & Lambert 2015), respectively. This large scatter between the different [Ba/Fe] at log t ∼ 6.5 cannot be fully explained with the adopted microturbulence values, and it certainly reflects fundamental issues in the analysis of such young stars. The solution of magnetic intensification proposed by Spina et al. (2020) can only partly explain the problem: our analysis suggests that we are witnessing an additional effect. A similar rising trend, although of much smaller extent, emerges when Y is considered. Like Ba, La belongs to the second peak of the s-process elements. As can be seen in the bottom left panel of Fig. 7, our measurements confirm solar values even at ages where [Ba/Fe] is extremely enhanced, as indicated in previous works (e.g. D'Orazi et al. 2012;Reddy & Lambert 2015;Mishenina et al. 2015). Finally, we find in the literature only a few studies where the Cu abundance was derived for OCs, namely Frasca et al. (2019), Casamiquela et al. (2020), Donor et al. (2020) (near-infrared measurements), and Spina et al. (2021). From the bottom right panel of Fig. 7, it is evident that [Cu/Fe] is solar (within the uncertainties) at all ages; we note the large scatter of the measurements, especially at younger ages.
In Fig. 8 Magrini et al. (2021). The recent production of the first-peak, Y, and second peak, La and Ba, s-process elements is mainly driven by the evolution of AGB stars, with lower percentages coming from massive stars in the early Galactic epochs (see e.g. Cescutti & Chiappini 2014 for a summary of the variety of possible scenarios). In low-mass AGB stars, the neutrons necessary for the production of s-process elements are mainly provided by the so-called 13 C pocket, which forms at the bottom of the convective envelope after each third dredge-up episode (Cristallo et al. 2009). The extension of the 13 C pocket plays a major role in the final production of neutron capture elements, and it can be parametrised in different ways. The GCE models adopted here      Fig.7. The cyan and magenta lines are the GCE models described in Magrini et al. (2021) with the MAGN stellar yields (Vescovi et al. 2020; continuous curves) and FRUITY (Cristallo et al. 2009;dot-dashed curves). In the [Cu/Mg] vs log t the models from Romano et al. (2010) with different stellar yields are considered: model 1 (green curve) with the Woosley & Weaver (1995) yields, model 5 with yields from Kobayashi et al. (2006) with ( consider the s-process yields from the FRUITY models, calculated by applying a simple exponentially decreasing profile of the convective velocities at the inner border of the convective envelope (Cristallo et al. 2009), and the MAGN models, a recent revision of the FRUITY yields which include the mixing triggered by magnetic fields (Vescovi et al. 2020), which could explain the behaviour of [Y/Mg] in open clusters at different Galactocentric distances (Magrini et al. 2021).
As can be seen, the GCE models considered here cannot reconcile the time evolution of Ba with that of La. In Fig. 8 we show the [Ba/La] time evolution; as expected, the production of Ba and La in the models is the same, thus neither of them can predict the (apparent) massive Ba production and the observed [Ba/La] rise in the last 100 Myr. As discussed in Mishenina et al. (2015), according to any s-process predictions high Ba yields should always be accompanied by high La and Ce yields, due to the presence of the magic number of neutrons (82) corresponding to these elements (e.g. Busso et al. 2001). This is at odds with what is observed in OCs. Mishenina et al. (2015) proposed that the intermediate neutron-capture (i-) process, which proceeds along a different path of neutron captures than the sprocess, is an additional source of Ba. According to their analysis, a combination of s-, r-, and i-processes may be able to reproduce the [Ba/La]> +0.20 dex observed in OCs for [Eu/La] ranging from −0.4 and +0.4 dex (see Fig.5 and Fig.6 in Mishenina et al. 2015). We cannot confirm this analysis in more detail because we cannot measure Eu in our stellar sample since the available Eu lines are too weak to be detected in stars with mild rotations as those in our sample. There are still large uncertainties concerning what stellar sources can host the i-process, and their efficiency in producing heavy elements. In particular, to explain the Ba excess in the YOCs a site that has become relevant only in the last 200 Myr would be needed, and that was not yet effective in contributing to the Ba abundance in the solar system. The discovery that the Ba enhancement could be explained by some observational issues (i.e. alteration of the spectral line) would help to solve the Ba puzzle without affecting our present understanding of the nucleosynthesis of Ba and La in the Galactic disc. From our discussion presented above we cannot make this conclusion either, and therefore the Ba puzzle remain unsolved.
Since our analysis has provided no clear answer for the excess of Ba, in order to trace the s-process element abundances at young ages, in particular for ages less than 500 Myr, we suggest looking for other elements and investigating further La and Ce. Theoretical models for Galactic chemical enrichment of heavy elements and theoretical GCE models including only the s-process should consider Ba and Y with extreme caution. The anomalous abundances of Y can also have an impact on the use of some chemical clocks, [Y/Mg] for example, as age indicators for young stars. As different studies have demonstrated (see e.g. Nissen 2015; Spina et al. 2018;Casamiquela et al. 2021 for solar twin stars), the [Y/Mg] ratio properly traces the age up to 500-700 Myr. Unfortunately, the two studies do not consider ages as young as the clusters we have analysed in this work. Our analysis suggests that the effects of alteration of the spectral lines could affect the relation [Y/Mg] versus age, with larger impact below 100 Myr, but probably also between 100 and 500 Myr. This is particularly evident in the bottom right panel in Fig. 8, where it is clear that the adopted GCE models cannot reproduce the increased abundance at young ages. Thus, we suggest that caution should be applied when using Y as a tracer of the s-process and as an age indicator below 500 Myr (e.g. [Y/Mg] and all the other ratios based on the Y abundances).

The time evolution of Cu
To follow the evolution of Cu, we adopt the GCE models presented by Romano & Matteucci (2007) and Romano et al. (2010Romano et al. ( , 2019. We have seen that most Cu production on Galactic scales is due to the weak s-process acting in massive stars. This mechanism depends on the initial metallicity of the stars. The neutrons originate mainly from the reaction 22 Ne(α,n) 25 Mg; the large abundance of 22 Ne during He-burning cores derives from the original CNO nuclei transmuted into 14 N in the H-burning ashes, followed by double α-capture on 14 N. Copper produced in this way is thus a secondary element. A small primary yield of Cu, 5 to 10% of its solar abundance derives from explosive nucleosynthesis in the inner regions of core-collapse supernovae (see e.g. Woosley & Weaver 1995, Rauscher et al. 2002, Pignatari et al. 2010. SN Ia models predict a negligible production of Cu during thermonuclear explosions (e.g. Iwamoto et al. 1999;Travaglio et al. 2005). Low-and intermediate-mass stars produce minor quantities of Cu as well. Thus, the models adopted in this paper assume that almost all Cu comes from massive stars, with a minor contribution from SNe Ia. After a short early phase in which the primary contribution from explosive nucleosynthesis in core-collapse SNe dominates, the evolution of Cu is regulated by the weak s-process. The models shown in Fig.8 include massive star yields from (i) Woosley & Weaver (1995) (green curve), (ii) Kobayashi et al. (2006) with (orange curve), or (iii) without (blue curve) the contribution from hypernovae (cf. models 4 and 5 of Romano et al. 2010), and (iv) the yields from rotating (for [Fe/H]< −1) and non-rotating (for [Fe/H]> −1) corecollapse SN progenitor provided by Limongi & Chieffi (2018b) (red line, corresponding to model MWG-11 of Romano et al. 2019 to which we refer for more details).
We note that some data in the range log t ∼ 7.5 − 8 (from this work and from Spina et al. 2021) fall below the predictions of the GCE models by about 10-60% (corresponding to 0.05-0.2 dex). Considering the uncertainties on the atomic physics affecting the measured abundances from Cu or Mg (or both), all the abundance ratios could be increased by 20-30%. Both Mg and Cu are mostly made in CCSNe, but they are produced through different nucleosynthesis processes. For instance, the present uncertainties in the 22 Ne(α,n) 25 Mg reaction (responsible for the Cu production) could justify a reduction of the [Cu/Mg] curve (see Fig.18 in Talwar et al. 2016). Additional variations of the order of the discrepancy observed could derive from the limitations of the GCE simulations adopted here, which predict average trends and cannot deal with local inhomogeinities that could play a role at this level of variation. Overall, taking into account the observational uncertainties and the limitations of the models, the red line (model MWG-11 adopting the yields by Limongi & Chieffi 2018b) is in fair agreement with the full dataset.

Concluding remarks
In this work we expanded our investigation of the chemical composition of young OCs, started in Baratella et al. (2020b), to shed light on the behaviour of the s-process elements at very young ages (t 200 Myr). In particular, we derived abundances of Cu i, Sr i and ii, Y ii, Zr ii, Ba ii, La ii, and Ce ii. For all the clusters we reported the very first determinations of [Cu/Fe]. Regarding IC 2391 and IC 2602, which are the most studied clusters in our sample, we presented the first determinations of [Sr/Fe]. On the other hand, we presented for the first time heavy-element abundances for NGC 2264 and NGC 2547.
Our measurements confirm the super-solar (0.25-0.70 dex) [Ba/Fe] abundances in the youngest population, a mild enhancement of [Y/Fe] (between 0 and 0.30 dex), and a solar-scaled abundance pattern for all the other s-process elements. We investigated several aspects in order to envisage possible solutions to the anomalous behaviour of the s-process element Ba.
From the comparison of spectral lines in the Sun and a solar analogue of ∼ 30 Myr, we note that the lines of some elements, for example Fe ii, Ca i, Ba ii, and Y ii, are stronger in the young star than in the Sun. On the other hand, La, Sr, and relatively weak lines of other elements, like Sc i or Sc ii, are almost identical. Sc lines can form at different depths to the other elements for which we observe a remarkable difference. Looking at the results obtained for the GBS, which are old and quiet stars, we do not reveal any anomalous trend. It is not clear what is altering the structure of the photosphere in very young stars and how this can modify the profiles of the spectral lines. From our analysis, the situation appears rather complex: the magnetic intensification is not sufficient to fully explain the large abundances. Both ionised and neutral lines are altered in the same way, but this alteration may vary with the optical depth, in the sense that lines at smaller log τ 5000 (i.e. in the upper layers) are more affected than lines forming deeper in the photosphere.
The solution proposed by Spina et al. (2020) of magnetic intensification and a pure photospheric, 1D-LTE treatment involving microturbulence does not seem to be sufficient to account for the observed pattern since in our case we are dealing with much younger stars (t 200 Myr). Recently,Şenavcı et al. (2021) analysed the young (∼ 30 Myr), active, and relatively fast rotator (∼ 17km s −1 ), the solar analogue EK Draconis. They derived precise atmospheric parameters and chemical abundances, and they studied the spot distribution on the stellar surface.  (2017), they concluded that the Ba overabundance is likely due to the assumption of depth-independent microturbulence velocity, but our analysis suggests that the explanation may be far more complex.
Overall, the anomalous behaviour of the s-process elements, in particular of Ba and Y with respect to La, cannot be reconciled with only nucleosynthesis models and Galactic chemical evolution predictions. NLTE effects cannot be similarly invoked since the corrections are not sufficiently large. Thus, we suggest that Ba should not be used as a tracer of the s-process elements for young stars along with Y. We instead promote the use of La or Ce as the most reliable tracer for the investigation of time evolution of the s-process elements (especially at recent Galactic ages). Possible solutions of the Ba puzzle, both from the spectra and the nucleosynthesis prospective, are still under investigation. Masseron et al. (2020b) reported anomalous Ba enhancements (ten times higher than the other s-process elements) in stars that are also anomalous in the lighter elements, specifically those rich in P (Masseron et al. 2020a). As discussed above, this pattern cannot be reproduced by an s-process model because of the nuclear properties of the isotopes involved. It remains to be seen whether the Ba anomaly in these P-rich stars has any connection with the Ba anomaly in YOCs; while the P-rich stars have [Fe/H] values of roughly −1, and are therefore not young, their overall puzzling nucleosynthetic pattern may represent a clue to the site of the i-process.
Finally, we note that Zr lines might provide us with reliable diagnostics for the first-peak elements because its lines form deeper in the photosphere than the two Y lines we have used here. However, we should also increase the number of spectro-scopic observations of very young objects in the 3000-5000 Å, where the best La, Zr, and Ce lines are. Future multi-object, highresolution spectrographs in the blue wavelength domain will give fundamental contributions to this framework.