Lithium depletion and angular momentum transport in F-type and G-type stars in Galactic open clusters

Open clusters provide clues to understand the evolution of Li7 at the surface of low-mass stars and its possible correlation with stellar rotation, which is a challenge for both stellar hydrodynamics and Galactic chemical evolution. We aim to quantify the efficiency of the transport processes for both angular momentum and chemicals that are required to explain simultaneously the observed behaviour of surface Li7 and rotation as well as the internal rotation profiles inferred from helio- and asteroseismology in F- and G-type main sequence stars. We apply the model for the transport of angular momentum and chemicals that we tailored in a previous work for solar-type stars to an extended range of initial masses and metallicities corresponding to F- an G-type stars in a sample of 20 Galactic open clusters. We evaluate its ability to explain the Li7, Be9, and rotation periods observations. Over the entire range of masses, metallicities, and ages explored, we reproduce the evolution of the surface rotation rates and predict, for the first time, the observed anti-correlation between the surface rotation rate and Li7 depletion as a consequence of the penetrative convection prescription. However, the ability of the model to reproduce the so-called Li7 dip centred around 6600K strongly depends on the adopted prescriptions for shear turbulence. It also requires a stellar mass dependence for the viscosity adopted for the transport of angular momentum, similar to the behaviour predicted for the generation and luminosity of internal gravity waves generated by stellar convective envelopes. We provide an efficient way to model G-type stars of different ages and metallicities successfully. However, the Li7 and Be9 dip constraints call for further hydrodynamical studies to better model turbulence in stars.


Introduction
Understanding the evolution of 7 Li (hereafter Li) in lowmass stars is one of the main challenges for stellar and Galactic astrophysics. Despite Li being a very scarce element, it is a tracer of Galactic chemical evolution (e.g. Spite & Spite 1982;Matteucci et al. 1995;Romano et al. 2001;Travaglio et al. 2001;Prantzos 2012) and of transport processes that occur in stellar interiors (e.g. Charbonneau & Michaud 1990; Montalban & Schatzman 1996;Montalbán & Schatzman 2000;Charbonnel & Talon 2005;Talon & Charbonnel 2010;Castro et al. 2016;Dumont et al. 2021, hereafter Paper I, and references therein). An overall picture of the different possibly involved processes is described in the reviews by Mathis (2013) and Aerts et al. (2019).
Observations of open clusters have provided numerous Li abundance data for stars of different ages over a large range of masses and metallicities (e.g. Boesgaard 1976;Duncan & Jones 1983;Cayrel et al. 1984;Balachandran et al. 1988Balachandran et al. , 2011 ⋆ e-mail: thibaut.dumont@unige.ch Soderblom et al. 1990Soderblom et al. , 1993bBoesgaard 1991;Thorburn et al. 1993;Garcia Lopez et al. 1994;Swenson et al. 1994;Jones et al. 1999;Boesgaard et al. 2003a,b;Sestito & Randich 2005;Anthony-Twarog et al. 2009Cummings et al. 2012Cummings et al. , 2017Deliyannis et al. 2019;Randich et al. 2020). It has been clearly evidenced that Li depletion increases with time and is linked to stellar mass (e.g. Deliyannis et al. 2000, for a review). At a given age, the Li behaviour as a function of the stellar effective temperature (T eff ) shows two specific patterns. On the one hand, photospheric Li abundances of G-type stars decrease with decreasing effective temperature (decreasing mass). On the other hand, a group of F-type stars with T eff centred around 6'600 K fall into the so-called Li dip which appears in open clusters older than ∼ 200 Myrs (e.g. Wallerstein et al. 1965;Boesgaard & Tripicco 1986;Hobbs & Pilachowski 1986;Soderblom et al. 1993a;Balachandran 1995;Steinhauer & Deliyannis 2004;Boesgaard et al. 2016).
Classical stellar evolution (accounting only for convection as a mixing process in stellar interiors) predicts noticeable Li A&A proofs: manuscript no. main depletion during the pre-main sequence (PMS) when the base of the convective envelope is deep enough to reach the Li-burning temperature, which happens only for less massive stars (i.e. below ∼ 0.9 -1.0 M ⊙ at solar metallicity; e.g. Bodenheimer 1965;Pinsonneault 1997). However, it does not predict any further surface Li variation until the first dredge-up occurs after the stars leave the main sequence (MS), in striking contrast with the observational evidence of the steepening along the main sequence of the Li-T eff trend of the cool side of the Li dip and of the Li dip itself.
The key role of rotation was pointed out as the cool edge of the Li dip coinciding with the so-called Kraft rotation break, as observed for instance in the Hyades (Boesgaard 1987) and in NGC 752 (Hobbs & Pilachowski 1986). The Kraft break (Kraft 1967) corresponds to the transition stellar mass (≃ 1.2 M ⊙ at solar metallicity) where important structural changes occur in main sequence stars. In particular, stars more massive than the transition value have an extremely thin convective envelope, implying weaker or ineffective magnetic braking compared to cooler, less massive stars with a thick convective envelope that can sustain efficient stellar wind magnetic braking (e.g . Schatzman 1962;Weber & Davis 1967;Matt et al. 2015;Kawaler 1988;Cummings et al. 2017;Deliyannis et al. 2019). The cool side of the dip also corresponds to the transition region where the convective envelope can efficiently generate internal gravity waves that transport angular momentum , 2005, and to an internal structure where the stellar core is radiative while hotter MS stars host a convective core.
Importantly, Li appears to be only one piece of a bigger puzzle that should also include the constraints from 9 Be (hereafter Be) that burns at a slightly higher temperature than Li (∼ 3.5 MK and 2.5 MK, respectively; e.g. Pinsonneault 1997). Indeed, while the Be dip coincides with the Li dip (e.g. Deliyannis et al. 1998;Boesgaard et al. 2001Boesgaard et al. , 2004bBoesgaard et al. , 2020, Be is hardly depleted in the Sun, solar-like stars, and G-type stars (Balachandran & Bell 1998;Boesgaard et al. 2003aBoesgaard et al. ,b, 2004aBoesgaard et al. , 2016Boesgaard et al. , 2020. The Be behaviour thus brings additional constraints to the possible origin of the observed behaviour of Li in F-and G-type stars and to the depth of the required mixing process. Last but not least, the difficulty is to find the actual connection between the transport of chemicals and the transport of angular momentum at play in stars of different spectral types along their evolution, to account simulta-neously for the internal rotation profiles that can be deduced for the Sun and for some other stars thanks to asteroseismology (Kosovichev 1988;Elsworth et al. 1995;Mathur et al. 2008;Eff-Darwich et al. 2008;Marques et al. 2013;Benomar et al. 2015;Eggenberger et al. 2019a;García & Ballot 2019).
In this work, we explore the chemical and rotational evolution of low-mass stars on the PMS and the MS for different stellar masses and metallicities that cover the range of Galactic open clusters with ages between 5 Myrs and 4 Gyrs. We study the effects and the relevance for these stars of different transport processes that we already explored and validated for the specific case of the solar-type stars (Paper I), and that depend on both mass and metallicity. In Sect. 2, we present the observational data that we used to constrain model predictions. In Sect. 3 we describe the input physics of the stellar models and the different processes implemented in the evolution code STAREVOL used for this work. In Sect. 4, we compare the observations for the Hyades and Praesepe open clusters with the predictions for Li, Be, and surface rotation of our non-standard model including meridional circulation, shear-induced turbulence, atomic diffusion, overshoot, and parametric viscosity and turbulence. In Sect. 5, we compare the model predictions with Li and rotational periods (P rot ) data in a sample of open clusters over a broad range in age and metallicity. In Sect. 6, we discuss model predictions for the internal rotation and compare to asteroseismic measurements. We summarise our results and conclude in Sect. 7.

Observational data
We use observational data for a sample of Galactic open clusters. Their names, ages, metallicities, and distances to the Galactic centre are given in Table 1 along with bibliographical references from which the Li and Be surface abundances and rotation periods of individual stars were taken. We only take into account non-binary stars with confirmed membership as mentioned or flagged in the reference papers cited in Table 1.

Lithium and Beryllium abundances
In this work, we consider Li spectroscopic data for a sample of 14 open clusters (see Tab (Asplund et al. 2009) as the original abundance of lithium. All the original papers give 1D local thermodynamic equilibrium (LTE) Li abundances, except Jeffries et al. (2002), where they give non-local thermodynamic equilibrium (NLTE) Li abundances for NGC 6633 stars. As lithium abundances are known to be sensitive to non-LTE effects, we computed the 3D NLTE corrections (∆ NLTE ) to all the 1D LTE lithium abundances, using the code Breidablik 2 by Wang et al. (2021) 3 . To do so, for each star, we used the [Fe/H], T eff , and logg values given in the original papers also providing the lithium abundances (see 1 A(X) = log 10 (N X /N H ) + 12, where N X is the number density of element X. 2 https://github.com/ellawang44/Breidablik 3 In the case of the NGC 6633 stars from Jeffries et al. (2002), we first reversed their NLTE correction using the code of Carlsson et al. (1994) to obtain the 1D LTE abundances and then computed the same 3D NLTE corrections as for the other data.  (Boesgaard et al. 2003a) and M 50 (Douglas et al. 2016). Distances to the Galactic centre are from Cantat-Gaudin et al. (2020), except for UMa (Ujjwal et al. 2020) and they are based on the Gaia DR2 measurements. Tab. 1). In the specific case of the Pleiades, Bouvier et al. (2018) do not give the surface gravity. We thus adopted log g = 4.4 for its stars, which appears to describe the typical surface gravity of F and G dwarfs 4 well. The absolute values of NLTE corrections rarely exceed 0.1 dex; we consequently do not expect a significant impact on the results. We also used the Be spectroscopic data as an additional constraint for the Hyades and Praesepe from Boesgaard et al. (2004a) and Boesgaard et al. (2016). Observations were directly extracted from Boesgaard et al. (2004aBoesgaard et al. ( , 2016 without any modification. We adopted the meteoritic abundance A( 9 Be) = 1.41 (Asplund et al. 2009) as the original abundance of beryllium.

Age
There is no self-consistent determination in the literature of the ages of all the clusters we consider here, and this task is out of the scope of this work. To be consistent with Paper I, we used the ages from Bossini et al. (2019, when available) for the clusters with Li abundance measurements. For the clusters with rotation period measurements, we used the ages from Godoy-Rivera et al. (2021, when available). For the five clusters with Li data that were not studied in those two works, we took the ages quoted in the respective observation papers (see Table 2). Except for the Pleiades, the differences in age determinations from the literature weakly affect our conclusions (see discussion in Sect. 5.1).

Effective temperature T eff
The effective temperatures used in this study were taken from the same original sources as the lithium (or beryllium) abundances for consistency (see Tab. 1). In rare cases, spectroscopic T eff are available (this is the case of M67 and NGC 2420). In most cases, T eff were determined using colour-temperature calibration relations from (B-V), (V-I c ), or (V-K) colour indices.
For the 14 open clusters presented in Table 1, at least six different relations and methods were used. These relations differ by their metallicity dependence, the zero-point of their temperature scale, the colours used, etc. The relative consistency between the different calibration relations was tested in several previous stud-A&A proofs: manuscript no. main ies (Huang et al. 2015;Casagrande et al. 2010;Meléndez et al. 2010). These studies emphasise that maximum differences of about 100 K can be found when applying the different temperature scales to dwarf stars in the metallicity and photometric domain considered here. We thus consider that the general agreement between the different temperature scales adopted in the original papers that we use in this work is satisfactory.

Surface and internal rotation
We compared the surface rotation predicted by our models with the observational data set gathered by , Gallet & Bouvier (2015), and Godoy-Rivera et al.  Table 1). We constrained our model predictions for internal rotation using the asteroseismic measurements obtained by Benomar et al. (2015) for a sample of MS field stars observed by Kepler. We selected a sub-sample of four stars with metallicities close to that of the Hyades.

General assumptions
We followed the conclusions of Paper I on different transport processes of chemicals and angular momentum (namely: rotation, penetrative convection, parametric turbulence, parametric viscosity, and atomic diffusion, as described below), which were tested with respect to Li abundances and surface and internal rotation constraints for solar-type stars. We used the nomenclature described in Paper I. For instance, model M ν R1 T 6.42 A is a model computed with an R1 prescription for rotation-induced turbulence and parametric viscosity (ν=ν add ), penetrative convection (A), parametric turbulence down to log T 0 = 6.42, at median initial rotation velocity (M), and with atomic diffusion. All the details about the corresponding prescriptions for the different processes are given below. To a broader range of masses and metallicities, we applied the prescriptions that were identified to be the most relevant for solar-type stars. Models were computed in the mass range between 0.8 M ⊙ and 1.5 M ⊙ (δM = 0.1M ⊙ ) for seven values of [Fe/H] (-0.4, -0.2, -0.1, -0.05, 0, +0.10, +0.15) that cover the metallicity range of the Galactic open clusters described in Table 1. Computations started on the Hayashi track at the beginning of the deuterium burning phase on the PMS that we consider as the time zero of the evolution.

Input physics
We used the same updated version of the stellar evolution code STAREVOL as in Paper I (for general information and previous versions see Siess et al. 2000;Palacios et al. 2006;Decressin et al. 2009;Lagarde et al. 2012;Amard et al. 2019), which we refer to for details and references. The models presented in this work were computed with the same inputs physics (equation of state, opacities, nuclear reactions, model atmosphere, and mass loss). We used the same values for the mixing length parameter (α MLT , assuming the Schwarzschild criteria for convective stability) and the initial abundances that resulted from model calibrations on the Sun that were carried out for classical models (no transport of chemicals beyond convection), and for models including both rotation and atomic diffusion (respectively models C and R1 from Paper I).

Atomic diffusion, penetrative convection, parametric turbulence
Atomic diffusion was implemented according to Paquette et al. (1986) and Thoul et al. (1994). We did not take radiative accelerations into account. Their impact starts, however, to be non-negligible for stars with effective temperature higher than ∼ 6'800 K, which corresponds to ∼ 1.4 M ⊙ at solar metallicity Richard et al. 2002;Deal et al. 2018Deal et al. , 2020. This is discussed in Sect. 7.
Penetrative convection was treated as an overshoot and computed following the formalism of Augustson & Mathis (2019) with the following expression for the diffusion coefficient: where D 0 = (υ conv × H p × α MLT )/3 is the convective turbulent diffusivity, with υ conv being the mean velocity of the convective elements obtained from MLT, α MLT being the mixing length parameter, and H p being the pressure scale height; r is the local radius; r bcz is the radius at the base of the convective zone; (v/v 0 ) is the ratio of the velocity of the convective elements when taking rotation into account for the non-rotating inviscid value; and d ov = 0.0325 (as calibrated by Paper I to reproduce the Li abundance in solar-type stars in very young open clusters). The coefficients λ = 6 × 10 −3 and µ = 5 × 10 −3 are as prescribed by Baraffe et al. (2017) based on the simulations of Pratt et al. (2017).
Parametric turbulence is defined according to Richer et al. (2000) and Richard et al. (2005) with the following prescription for the diffusion coefficient: where T 0 is a free parameter that sets the depth of the maximum efficiency of the mixing depending on the value of the atomic diffusion coefficient He (D He ). Initially introduced by Richer et al. (2000) to counteract an impact of atomic diffusion in AmFm stars that was too strong, we introduced it in our models with rotation to increase the surface Li depletion predicted on the MS and reproduce the observations for solar-type stars and for the Sun. For all the reference models, we adopted log T 0 = 6.42 as calibrated by Paper I. We also present a set of models for the Hyades metallicity ([Fe/H] = + 0.15 dex), where we increased this depth to log T 0 = 6.5.

Angular momentum evolution and rotation-induced mixing
Stellar rotation was implemented in STAREVOL as described by Amard et al. (2016Amard et al. ( , 2019 and Paper I. We adopted the shellular rotation hypothesis developed by Zahn (1992), Maeder & Zahn (1998), and Mathis & Zahn (2004) to describe the transport of angular momentum and chemicals by meridional circulation, treated as an advective process for angular momentum, and turbulent shear (vertical and horizontal), treated as diffusive processes. We followed Eggenberger et al. (2012bEggenberger et al. ( , 2019b and Spada et al. (2016) who introduced an additional parametric vertical viscosity ν add in order to flatten the internal rotation profile as evidenced by helio-and asterosismology of low-mass stars. The transport of angular momentum hence obeys the following advection-diffusion equation: where ρ, r, Ω, U 2 , and ν v are the density, the radius, the angular velocity, the meridional circulation velocity, and the vertical shellular component of the turbulent viscosity, respectively. For the reference models, we adopted ν add = 3.5 × 10 4 cm 2 s −1 in the entire radiative region. This value was calibrated in Paper I on the solar angular velocity profile provided by helioseismology. We also present a set of models for the Hyades metallicity ([Fe/H] = + 0.15 dex) where we modified this value depending on the initial stellar mass (ν ′ add = 3.5 × 10 5 , 1.0 × 10 5 , 2.5 × 10 5 , 6.5 × 10 5 , and 8.5 × 10 5 cm 2 s −1 for the 1.5, 1.4, 1.35, 1.3, and 1.2 M ⊙ models, respectively). These values were kept constant during the entire evolution.
We explored two combinations of prescriptions for turbulence shear referred to as R1 and R2 as in Paper I; R1 includes prescriptions from Mathis et al. (2018) and Zahn (1992) for the horizontal diffusivity D h and the vertical diffusivity D v , respectively, and R2 includes prescriptions from Zahn (1992) and Talon & Zahn (1997) for D h and D v , respectively. The detailed expressions of the four turbulent diffusion coefficients can be found in Appendix B of Paper I. We used the same parameters as in Paper I for the extraction of angular momentum at the stellar surface due to magnetised winds (m=0.22, p=2.1, χ = 14, K = 7.5 × 10 30 erg unless otherwise indicated) according to the Matt et al. (2015) formalism. Models were computed for three values of the initial rotation period on the PMS: 1.6, 4.5, and 9.0 days, which are referred to as fast ( F R), median ( M R), and slow ( S R) rotating models, respectively. The disc coupling timescale was set at τ disc = 2.5 Myrs for the fast rotators and at τ disc = 5 Myrs for the median and the slow rotators, in agreement with Gallet & Bouvier (2015) and Amard et al. (2019).
To summarise, the best models for solar-type stars that we developed in Paper I include the self-consistent treatment of physical processes as well as parametrised additional transports for angular momentum and chemicals as summarised in Tab. 2. The parameters are strongly constrained by observations (see col. 3 Tab. 2) and cannot be varied freely. In the present study, we explore a possible variation in the efficiency of these additional processes with metallicity and initial mass.

Model predictions for Li, Be, and surface rotation -The Hyades and Praesepe test case
In this section, we explore the predictions of model ν R1 T 6.42−6.5 A over a range of masses and metallicities for three different rotation rates. We recall that this model was developed for the specific case of solar-type stars in Paper I. We discuss the general impact of mass and metallicity on Li depletion, and then we compare the predictions for the behaviour of surface Li, Be, and rotation rates to the observations in the Hyades and Praesepe, and to the predictions of classical models and of model M ν R2 T 6.42 A computed with different prescriptions for rotation-induced turbulence (see Tab. 3).
4.1. Impact of metallicity and mass on Li evolution Figure 1 shows the evolutionary tracks in the Hertzsprung-Russell diagram (PMS and MS) and the predicted surface Li abundance as a function of T eff for a selection of models (all computed with median rotation). We see the well-known impact of varying the mass and the metallicity on the evolution tracks, and the global consequence on Li depletion. When metallicity decreases for a given mass, or when the mass increases for a given metallicity, the model is hotter and brighter. As a consequence, its convective envelope retracts more rapidly at the beginning of the Hayashi track on the PMS, and it is thinner on the MS. Hence, its base is more distant from the layers where Li is burning, which overall leads to lower Li depletion, although the quantitative details depend on the efficiency of the transport process(es) that connect(s) the base of the convective envelope to the Li-free region (see below).
In all the models shown in Fig. 1, Li decreases along two successive episodes, first on the PMS, and later along the MS. The PMS depletion episode is first due to Li burning in the convective envelope at the beginning of the Hayashi track and later to overshoots that are both dependent on the size of the convective envelope and its timescale for retraction along the PMS. As a consequence, PMS depletion is minute for stars with ZAMS 5 T eff higher than ∼ 6'500 K, and it increases with decreasing stellar mass and increasing metallicity. The MS depletion episode is due to the combined effects of rotation-induced mixing and parametric turbulence. The first process is more efficient in cooler stars (i.e. less massive at a given metallicity, or more metal-rich at a given mass), which undergo more significant extraction of angular momentum by their magnetised winds (see e.g. Fig. 6 in Amard et al. 2019). Similar mass and metallicity dependencies exist in the case of the second process because in all of our models we assume, due to Eq. (2), the same value for the free parameter T 0 which is closer to the temperature of the base of the convective envelope in less massive or more metal-rich stars. In summary, overall Li depletion is stronger in less massive, more metal-rich stars.

The Hyades and Praesepe
Here we focus on the well-studied Hyades and Praesepe for which Li, Be, and surface rotation data are available. These two open clusters are close in age (0.72 Gyr and 0.75 Gyr, respectively) and metallicity ([Fe/H] = +0.13±0.05 and +0.16±0.08 dex, respectively) according to the references quoted in Table 1. Figure 2 shows the observed Li abundances and rotation velocities (V sin i) as a function of the effective temperature from Cummings et al. (2017). Both clusters exhibit similar patterns, namely the so-called Li dip between ∼ 6'400 and 6'800 K, the Li-T eff decrease in G-type stars on the cool side of the dip, and the well-know break in rotation velocity (≈ 6'400-6'500 K) for dwarf stars later than the F4 spectral type (Schatzman 1959(Schatzman , 1962Kraft 1967;Boesgaard 1987). In the same figure, we show the predictions at 0.75 Gyr of the classical models, and of the complete models ν R1 T 6.42  Notes. Transport processes (Column 1) and values of the free parameters (when relevant) (Column 3) with the corresponding observational constraints adjusted (Column 4) for the best model for solar-type stars (Paper I) and adopted in this study. The flags in Column 2 indicate the transported quantity (C for chemicals and AM for angular momentum). As has been long known in the literature, and as shown in Fig. 2, classical models (model C) that account solely for convection predict Li depletion on the PMS only, when the size of the convective envelope is large enough to reach the Li-burning temperatures along the Hayashi track (see Sect. 4.1). While the predicted mass and T eff dependence of the Li abundance after the PMS depletion is similar to the observed pattern for G-type stars, the predicted surface Li abundances are too high compared to the observations in the Hyades and Praesepe, as well as in the other clusters of different ages considered in this work (see Sect. 5).
Models ν R1 T 6.42 A and ν R1 T 6.5 A predict rotation velocities that account for the observed V sin i trends well (including the position in effective temperature of the Kraft break) in both open clusters, thanks to the extraction of angular momentum due to magnetised winds (see also Fig. 5). However, on the cool side of the Kraft rotation break (see insert in Fig. 2), models reproduce the lower observational envelope well, but not the entire observed spread. It is mainly the case for the lowest stellar masses where the parametric viscosity is more efficient.
Models ν R1 T 6.42 A and ν R1 T 6.5 A also provide a very good fit to the Li data in G-type stars on the cool side of the Li dip where angular momentum extraction is maximal, that is, for stellar masses lower than 1.3M ⊙ for the metallicity of these clusters. At the age of the Hyades and Praesepe, the increased (with respect to the classical model C) Li depletion is mainly due to penetrative convection (Eq. 1) along the PMS, and only slightly due to turbulent mixing at the beginning of the MS. As discussed in Paper I, the Augustson & Mathis (2019) expression for penetrative convection leads to larger Li depletion for slow rotators than for fast rotators because of the influence of the rotation rate on the depth of the overshoot via the ratio v/v 0 in Eq. (1). This anticorrelation with the rotation rate, which has long been observed (e.g. Bouvier 2008;Bouvier et al. 2018;Arancibia-Silva et al. 2020), is obtained here for the first time for the entire mass range considered at the Hyades and Praesepe age. 6 In addition, the 6 This anti-correlation was initially obtained by Somers & Pinsonneault (2015) who invoked a radius inflation de-predicted Li spread induced by the different initial rotation rates assumed here (Sect. 3.4) is amplified in lower mass stars where the base of the convective envelope is deeper and closer to the Li burning layers, which also leads to more efficient parametric turbulence on the lower mass end. For this later process, we present model predictions for two values of the parameter T 0 from Eq. (2), that is, log T 0 = 6.42 as calibrated in the Sun and log T 0 = 6.5 which better fits the data spread in cool stars of the Hyades and Praesepe. Overall, our complete model ν R1 T 6.42−6.5 A reproduces fairly well the Li-T eff -V sin i trend and the observed Li dispersion on the cool side of the Li dip (i.e., G-type and late F-type stars). In this T eff domain, it also accounts for the Be data in the Hyades and Praesepe, as shown in Figs. 3 and 4.
However, models ν R1 T 6.42−6.5 A predict little Li depletion (both on the PMS and the MS at the ages of the Hyades and Praesepe) for stars more massive than 1.3M ⊙ that have small convection zones that are vanishing on the MS and which undergo a negligible angular momentum braking. As a consequence, it does not reproduce the depth of the Li dip, nor the Li rise on the cool edge of the dip. Around 6'600 K, where many stars only have Li upper limits, the maximum Li depletion predicted reaches only ∼ 1 dex in the case where we assume slightly deeper parametric turbulence (i.e. log T 0 = 6.5). This low depletion mainly results from the choice of prescriptions for the horizontal and vertical diffusivities D h and D v (taken here from Mathis et al. 2018 andZahn 1992, respectively;see Sect. 3.4). This is in contrast with previous studies on the hot edge of the Li dip (Palacios et al. 2003), where the use of the prescriptions by Zahn (1992) and Talon & Zahn (1997) for D h and D v , respectively, led to a good match between model predictions and observations. As can be seen in Fig. 3, much stronger Li and Be depletion is indeed obpendent on the rotation velocity to explain the anti-correlation at the Pleiades age. However, even though it has been confirmed in the Pleiades (  where we used the same prescriptions as in that earlier paper. As already evidenced in Paper I for the solar case, this combination however leads to Li depletion that is too strong in F-and G-type stars for the value of ν add adopted in that paper and here (3.5 × 10 4 cm 2 .s −1 ). Additionally, reducing the torque parameter K for the wind from Matt et al. (2015), from 7.5 × 10 30 erg to 1.5 × 10 30 erg, improves the comparison on the hot edge of the Li dip and of the Be dip, but not on the cool edge as can be seen in Figs. 3 and 4 for model M ν R2 T 6.42 A.K ′ . Talon & Charbonnel (1998) reached similar conclusions as Palacios et al. (2003), using the same turbulence prescriptions (D h and D v from Zahn 1992 andZahn 1997, respectively). These authors interpreted the rise in Li on the cool edge of the dip as the signature of the appearance of a process that transports angular momentum more efficiently than meridional circulation and turbulent shear. , 2005 then showed that the generation of internal gravity waves (hereafter IGW) by the stellar convective envelope becomes efficient inside the Li dip and increases on its cool edge. In their model, the maximum efficiency of the IGWs in terms of angular momentum transport is expected at T eff around 5'800-5'900 K, before it decreases in cooler stars (see Fig. 6 of Talon & Charbonnel 2003 and Fig. 4 of Talon & Charbonnel 2004). Last but not least, their model also accounts for the internal solar rotation (Charbonnel & Talon 2005). It thus has the proper mass and T eff dependence to be the required transport mechanism involved, contrary to other processes that could also potentially transport angular momentum (e.g. the Taylor-Spruit dynamo Spruit 2002;Eggenberger et al. 2005Eggenberger et al. , 2010. Following this theoretical trend for IGWs, we computed the M ν ′ R2 T 6.42 A of various masses assuming a T eff dependence of the parametric viscosity ν ′ add (increasing value with decreasing mass on the cool edge of the Li dip) in order to mimic a transport of angular momentum as predicted for low-mass stars. Although the use of a uniform and constant parametric viscosity within the radiative layers leads to a very different angular momentum profile than that shaped by IGW (compare Fig.7  efficient internal transport of angular momentum compared to that driven by meridional circulation and turbulent shear (see Fig. 7 which is discussed in detail in Sect. 7), hence reducing their efficiency for the transport of chemicals. The ν ′ add values required to reproduce the shape of the Li dip are 3.5 × 10 5 cm 2 .s −1 and 1.0 × 10 5 cm 2 .s −1 for the 1.5 and 1.4 M ⊙ models, respectively, on the hot edge of the Li dip, and 2.5 × 10 5 , 6.5 × 10 5 , and 8.5 × 10 5 cm 2 .s −1 for the 1.35 M ⊙ , 1.3 M ⊙ , and 1.2 M ⊙ models, respectively, on the cool edge of the dip. On the other hand, the Li abundance in the 1.0 M ⊙ model is well reproduced for ν ′ add = 2.5 × 10 5 cm 2 .s −1 . However for this mass, it leads to a reduction that is too strong of the surface velocity (see Figs. 5 and 7), which makes it hard to reconcile with P rot measurements in the Hyades and Praesepe. As shown in Figs. 3 and 4, these assumptions explain the Be plateau in cool stars, within the observational uncertainties. However, the M ν ′ R2 T 6.42 A models do not account for the large Be depletion in the Be dip nor in its cool edge. This conclusion holds for all the initial rotation rates explored here.
It thus appears very difficult to reconcile the Li, Be, and surface rotation rate constraints in the test case of the Hyades and Praesepe. While the ν R1 T 6.42 A models account for the surface rotation well over the entire mass domain explored, and for the Li and Be abundances in G-type stars, they fail to explain the Li and Be dips in these two clusters. On the other hand, M ν ′ R2 T 6.42 A models with the parametric viscosity for the transport of angular momentum, dependent on the stellar mass of the star and adjusted to fit the Li abundances along the entire mass range, predict a Be dip that is not as deep as the one observed, and they lead the 1 M ⊙ model to rotate too slowly at the age of the two clusters (0.72-0.75 Gyr). In Sect. 5 we evaluate the respective compatibility of the models ν R1 T 6.42 A and M ν ′ R2 T 6.42 A with the surface Li and rotation rates in a sample of open clusters, before discussing in Sect. 6 the impact of the different assumptions for the transport of angular momentum on the internal rotation profiles of the models in the light of asteroseismic constraints. In Sect. 7, we provide a summary of the successes and difficulties of the different assumptions, and we discuss their meaning in terms of the uncertainties of the different processes involved.

Model predictions for Li and surface rotation -Comparison to open clusters of different ages and metallicities
In this section, we compare the predictions of models ν R1 T 6.42 A and M ν ′ R2 T 6.42 A to rotation rates and surface Li abundances for a sample of open clusters with different ages and metallicities (see Table 1). Be cannot be used here, as there are not enough available data for clusters other than the Hyades and Praesepe. In the ν R1 T 6.42 A case, we computed models for the actual [Fe/H] of the individual clusters and for three initial rotation rates (slow, median, and fast), while in the M ν ′ R2 T 6.42 A case, we used the models at the metallicity of the Hyades and Praesepe ([Fe/H] = +0.15 dex), and for median rotation discussed previously. The values of A(Li), T eff , and P rot predicted by the M ν R1 T 6.42 A model are given in Appendix A for the different masses and cluster ages.

P rot versus mass
In Fig. 5 we compare the predicted surface rotation periods of the ν R1 T 6.42  Table 1 for the P rot . Godoy-Rivera et al. (2021) have studied the difference in mass estimates that can be derived from using different families of isochrones and they show that it is quite small (∆M ≈ 0.05M ⊙ ).
We confirm the conclusions reached in previous studies (e.g. Krishnamurthi et al. 1997;Allain 1998;Amard et al. 2016Amard et al. , 2019: the evolution of the rotation period can be divided in several phases that our models succeed to predict. Firstly, the large dispersion at very young ages for all masses (NGC 2362 to IC 2391, first row of Fig. 5) results in our models from the assumption of constant surface angular velocity as long as disc-locking is efficient. Secondly, the progressive slow down of first, slow and median, and finally fast rotators of decreasing mass in clusters of increasing age (from α Per to M 37, second row of Fig. 5) is very well reproduced by our models as the result of the secular evolution of the radius along the PMS and the efficient wind braking on the early MS. Finally, for clusters older than Praesepe (≈ 750 Myrs), we observe the convergence of the three families of rotators into one single sequence recovering the observational law by Skumanich (1972). In addition, the sequence presents a change in slope around M ≈ 1.2 M ⊙ in the P rot versus mass diagram (last row of Fig. 5), resulting from the magnetic braking by the stellar winds. This change in slope corresponds to the Kraft break and to the mass domain where the convective envelope becomes very thin. In the Hyades and Praesepe, we also show the predictions of the model M ν ′ R2 T 6.42 A introduced to fit the cool edge of the Li dip (see Sect. 4.2). The predicted rotation periods for models more mas-Article number, page 9 of 18 A&A proofs: manuscript no. main  Table 1. sive than 1.2 M ⊙ are fully compatible with observations, while the 1 M ⊙ median rotator model is far too slow, as already discussed in Sect. 4.2. As in Amard et al. (2019), who use the same prescriptions for horizontal and vertical turbulence and for the braking law as in model ν R1 T 6.42 A (with a parameter K=7.0×10 30 erg instead of 7.5×10 30 erg) but who do not include the parametric viscosity for the internal transport of angular momentum, our model predictions for the older clusters present a larger negative slope on the lower mass side, with the lower mass models spinning slower than the observed stars. For the more massive stars in NGC 6819, the extraction of angular momentum, which is directly linked to the depth of the convective envelope, becomes inefficient in the models, which spin faster than the observed stars. Finally, for the Pleiades we obtain a poor agreement with observations for the lowest stellar masses (0.8 and 0.9 M ⊙ ) compared to what we get for all the other clusters and to what was obtained in Amard et al. (2019). This is essentially due to the difference in age adopted for this cluster in both studies (125 Myrs from Barrado y Navascués et al. (2004) in Amard et al. (2019) instead of 87 Myrs as assumed here), as shown in Fig. 5 where the orange diamonds are the predictions of our models at 125 Myrs, which better account for the data.
To summarise, the comparison between the predictions of our ν R1 T 6.42 A and M ν ′ R2 T 6.42 A , and the observational data for the rotation period is rather satisfactory, although some peculiar features remain difficult to reproduce without further adjustments. This work confirms the global theoretical behaviour obtained by Amard et al. (2016Amard et al. ( , 2019 despite the fact that we included an additional diffusive source of transport of angular momentum in the present study.

Lithium -T eff
In Fig. 6, we compare the predictions of models ν R1 T 6.42 A for the surface Li abundances to data within the cluster sample (see Sect. 2 and Table 1), and we report the predictions of models M ν ′ R2 T 6.42 A discussed in Sect. 4 for the most metal-rich clusters (αPer, Hyades, Praesepe, and NGC 6819). Importantly, over the  Table 1 for references; grey squares and triangles are abundance determinations and upper limits, respectively). The typical observation errorbars from the original papers are indicated by a cross in each panel. Coloured diamonds are the predictions of the models at the age and metallicity of the corresponding clusters. The ν R1 T 6.42 A models are shown for three initial velocities (slow, median, and fast are the cyan, blue, and violet diamonds, respectively), and models M ν ′ R2 T 6.42 A for median rotation only (magenta diamonds). Models from left to right in each panel correspond to masses between 1.5 M ⊙ (warmer) and 0.8 M ⊙ (cooler). entire mass and metallicity range explored, the faster the initial rotation velocity, the lower the predicted Li depletion. As already discussed in Sect. 4.2, this is due to penetrative convection acting on the PMS, which we simulated with the Augustson & Mathis (2019) prescription.
The Li-T eff relation for G-type stars is well reproduced by models ν R1 T 6.42 A at all ages from the youngest (IC 2602 and IC 2391) to the oldest cluster (M67) with solar or super-solar metallicities. For the young clusters with subsolar metallicity (M35, NGC 6633, and UMa), the cooler models, corresponding to lower mass stars (see Table in Appendix A), however predict larger A(Li) values than observed. In these models, the convective envelope is thinner due to the lower metallicity and the penetrative convection, which is the main process responsible for the Li depletion on the PMS and it is less efficient than at higher metallicity. This can be directly seen when comparing the panels showing NGC 6633 and UMa clusters as well as the Hyades and Praesepe clusters, which have similar ages but different metallicities. This suggests a stronger metallicity dependence of the Li depletion in cool stars compared to that predicted. Data for the cooler stars in the most metal-poor and oldest cluster (NGC 2243) would be required to confirm this trend at older ages.
We see in Fig. 6 that the Li dip is not present in the youngest clusters, and that its depth increases with time along the MS (see also Wallerstein et al. 1965;Boesgaard & Tripicco 1986;Hobbs & Pilachowski 1986;Soderblom et al. 1993a;Balachandran 1995;Steinhauer & Deliyannis 2004;Anthony-Twarog et al. 2009;Cummings et al. 2012;Boesgaard et al. 2016 higher T eff compared to the younger higher metallicity clusters (Cummings et al. 2012;François et al. 2013). Since models ν R1 T 6.42 A lie below the observational points on the cold edge of the Li dip, we expect that model M ν ′ R2 T 6.42 A , assuming a similar dependence between ν ′ add and T eff as assumed to fit the Hyades and Praesepe, would predict an even stronger Li depletion and would be incompatible with the data of NGC 2243.
We thus reach similar conclusions as in Sect. 4.2, for a relatively large range of metallicities and ages. Model ν R1 T 6.42 A initially developed to reproduce surface and internal constraints for solar-type stars nicely explains the general trend of the data for both the rotation and Li behaviours over the mass and the metallicity ranges probed in this study, except for the Li dip. On the other hand, model M ν ′ R2 T 6.42 A with turbulent viscosity for the transport of angular momentum parametrised to fit the dip in the Hyades and the Pleiades predicts rotation rates that are too slow for the lower mass models. Importantly, and for the first time, we predict that over the entire mass and metallicity range explored here, penetrative convection that is efficient at the early stages of the PMS leads to larger Li depletion in the slower rotators. This anti-correlation between the rotation rate and the Li abundance, which builds early on the PMS, remains constant over time.

Model predictions for internal rotation and comparison to asteroseismic constraints
In this section, we present the predictions of models M ν R1 T 6.42  In Fig. 7, we present the internal rotation profiles for the 1.0, 1.2, 1.3, and 1.4 M ⊙ models at three different ages on the MS. For G-type stars, the M ν R1 T 6.42 A model that reproduces all surface observational constraints (see previous sections) predicts a small radial differential rotation that decreases with time and is compatible with helioseismic data as already shown in Paper I. For the more massive models corresponding to the Li-dip stars in the Hyades, Praesepe, and NGC 6819, this model predicts the same behaviour. The M ν R2 T 6.42 A model predicts larger radial differential rotation for all masses and ages than the M ν R1 T 6.42 A model, which is in agreement with the associated strong Li depletion discussed in previous sections. This results from the much higher efficiency of the vertical turbulent shear D v when using the Talon & Zahn (1997) prescription (see also Paper I). For the M ν ′ R2 T 6.42 A model, which introduces an enhanced parametric turbulent transport of angular momentum to fit the cold edge of the Li dip, the radial differential rotation is strongly reduced at all ages, especially in the 1.2 M ⊙ model with the larger value of ν ′ add , which presents an almost flat rotation profile at all ages displayed on the figure. Concerning the case of the 1.4 M ⊙ model, although the angular momentum contrast between the surface and the core is of the same order in models M ν R1 T 6.42 A and M ν ′ R2 T 6.42 A , the angular velocity gradient near the base of the convective envelope is much larger in the later model (see change in slope in magenta profiles around r = 0.9 R ⋆ in Fig. 7), thus leading to an enhanced shear and associated turbulent transport in the exact region where the connection between the convective envelope and Li burning zone is made. This leads to the strong Li depletion previously discussed in this model (see Sect. 4.2). On the contrary, the shear is near to null at the base of the convective envelope in model  Fig. 8 shows the mean angular velocity in the radiative interior (i.e. excluding the convective core if present 7 , which is consistent with asteroseismic data that do not take it into account) as a function of the central hydrogen mass fraction in the 1.2, 1.3, and 1.4 M ⊙ models at [Fe/H]=+0.15 dex. We compare this to a selected sub-sample of four field stars with metallicities close to that of the Hyades for which the average rotation in the radiative interior was estimated from asteroseismology by Benomar et al. (2015). Our predictions are of the same order of magnitude as the results from Benomar et al. (2015), but the predicted rotation is faster than the values they inferred for the selected targets. Models M ν ′ R2 T 6.42 A lead to a better agreement with the values derived from asteroseismic data than models M ν R2 T 6.42

A
. A more precise comparison with the results in Benomar et al. (2015) is prevented by the use of MESA stellar models in their study including very different input physics and settings. In particular, core hydrogen mass fractions X c cannot be compared.
Additional data for different stellar masses would be required to better discriminate between the different prescriptions for the transport of angular momentum, especially for cool stars. Nevertheless, this comparison supports the need for strong transport of angular momentum on the MS for stars in the Li-dip region, as obtained when increasing the strength of ν add as assumed in models M ν ′ R2 T 6.42 A .

Summary and discussion
The need for additional transport processes beyond atomic diffusion and so-called Type-I rotation-induced processes (turbulent shear and meridional circulation) for the transport of chemicals and angular momentum has long been reported (see references in Sect. 1). This work follows Paper I analysing the impact of several processes, using state-of-the-art prescriptions for specific mechanisms (rotation-induced turbulence and penetrative convection in particular), as well as parametric prescriptions for others (parametric turbulence for chemicals and parametric viscosity for angular momentum) that were calibrated to reproduce the evolution of the surface Li abundance and rotation rates as well as the internal angular velocity profiles in the Sun and solar-type stars. The aim of the present paper is to study the impact of these processes (using the solar-type calibrations) on Fand G-type MS stars at different metallicities and to test them against data in Galactic open clusters over a large range in ages. Although the adopted ages come from different references and are not fully consistent, our conclusions are robust given the uncertainties on their estimate. Our predictions support, however, the Pleiades to be older than assumed in this work, which is in agreement with other age determinations for this cluster. The so-called ν R1 T 6.42−6.5 Li, Be, and surface rotation rates available in a sample of open clusters of different ages and metallicities. We also compared the predicted internal rotation profiles with asteroseismic constraints in field MS stars with [Fe/H] close to that of the Hyades. Both ν R1 T 6.42−4.5 A and M ν ′ R2 T 6.42 A explain the main general trends observed between Li depletion and stellar mass, age, rotation, and metallicity covered in this study. Thanks to the prescription we used for penetrative convection (Augustson & Mathis 2019), an anti-correlation between Li depletion and surface rotation build up early on the PMS, and this remains throughout the evolution on the MS. Theoretical surface Li abundance and T eff were correlated as observed in cluster stars, with cooler and lower mass stars being more Li depleted than hotter and more massive ones (except in the Li dip when this feature is present). Our models also recover the anti-correlation between metallicity and Li depletion efficiency that is observed in open clusters stars. Metal-poor stars are less Li-depleted than their metal-rich counterparts, as expected from the metallicity dependence of the location of the base of the convective envelope that affects the PMS Li depletion. However, Li depletion in our metal-poor models is more modest than in open cluster stars.
Model ν R1 T 6.42−6.5 A succeeds to reproduce the Li, Be, and rotation rates observations in G-type stars at all ages, and in F-type stars in the young clusters where the Li dip has not started to form yet. This is achieved thanks to the combined effects of penetrative convection, rotation-induced processes, and parametric turbulence, whose efficiencies vary with both the stellar mass and the metallicity. The impact of these mechanisms on Li and Be depletion becomes minute in the more massive stars that have a very thin convective envelope, more distant from the Li and Be burning layers. As a result, model ν R1 T 6.42−6.5 A is not able to reproduce the Li and the Be dips centred around ∼6'600 K that appear later on the MS, although it reproduces the rotation rates over the entire mass range covered in this study. On the other hand, model M ν R2 T 6.42 A computed with the same value of the parametric viscosity, but with different prescriptions for horizontal and vertical turbulence (D h and D v ) that lead to more differential rotation (hence more efficient rotation-induced mixing) in the stellar interior, predicts too much Li depletion on the cool side of the dip, which is as expected from previous studies.
Building on the assumption that an additional transport process such as internal gravity waves would be the main driver for the transport of angular momentum in stellar interiors on the cold edge of the Li dip, we introduced a mass-dependent viscosity ν ′ add with a maximum efficiency at ≈ 5'900 K as predicted in  and Talon & Charbonnel (2004) for the IGW excitation by the convective envelope and their luminosity. Although the assumption we make of a parametric viscosity ν ( ′ ) add being uniform within the radiative interior and constant with time is a very crude parametrisation of a transport of angular momentum within stars, the M ν ′ R2 T 6.42 A model predicts a good agreement with Li abundances on the cold edge of the dip in open clusters of different ages at the metallicity of the Hyades (this was not tested for other metallicities). This is due to the flattening of the angular velocity gradient within the stellar interior, and in particular between the base of the convective envelope and the Li burning depth, implied by higher values of the parametric viscosity. A further improvement to our model now requires the full treatment of the transport of angular momentum by IGWs as in Charbonnel & Talon (2005), also taking into account recent developments on the generation and the behaviour of IGWs (e.g. Pinçon et al. 2016;Augustson et al. 2020;Ratnasingam et al. 2020). The Tayler-Spruit dynamo process that was updated in recent studies (e.g. Fuller et al. 2019; A&A proofs: manuscript no. main Eggenberger et al. 2019a,b,c) is also a candidate for additional transport of angular momentum and should also be tested regarding the results of the present work. The Be depletion obtained with model M ν ′ R2 T 6.42 A is on the other hand too modest to explain the depth of the Be dip observed in a couple of clusters. The acquisition of more Be observational data in open clusters of different ages and metallicities, in addition to those in field stars (e.g. Delgado Mena et al. 2012), would be an important key to better constrain models.
Importantly, the transport of angular momentum simulated with the parametric viscosity values adopted in our different models is much more efficient than that driven by meridional circulation and turbulence. Both the ν R1 T 6.42 A and the M ν ′ R2 T 6.42 A models predict mean angular rotation velocity in the radiative interior that are compatible with the values inferred from asteroseismology in field MS with the same metallicity as the Hyades and the Pleiades. Additional asteroseismic data to probe the internal rotation of young MS stars over a large metallicity range would be valuable.
Reflecting on the modelling of turbulent shear, we have to conclude that none of the combinations for the prescriptions of the horizontal and vertical turbulence (D h and D v ) we used is able to reconcile all the surface constraints used in this study. Model ν R1 T 6.42−6.5 A , tailored for solar-type stars, includes the Mathis et al. (2018) prescription for the horizontal shear diffusivity D h , which is the only prescription including the contributions of both the horizontal and vertical shear on the horizontal turbulent transport, and the Zahn (1992) prescription for the vertical shear diffusivity D v , which does not take into account the effects of thermal and molecular diffusivity on the vertical turbulent transport, but it is the only one that has been validated by direct numerical simulations (Prat & Lignières 2013;Garaud et al. 2017). This model thus includes the apparently most robust available prescriptions, but it cannot explain the Li dip even when combined with additional sources for turbulent transport such as D T 0 and ν add .
On the other hand, the use of the Zahn (1992) prescription for the horizontal shear diffusivity (D h ) and that of Talon & Zahn (1997), which includes the effect of thermal and molecular diffusion for the vertical shear diffusivity (D v ) in model M ν ′ R2 T 6.42 A , requires higher values for the parametric viscosity for the internal transport of angular momentum, with a mass-dependent efficiency as expected from IGW. The choice of these prescriptions thus impacts our ability to constrain the other transport processes (see Meynet et al. 2013;Amard et al. 2016) as additionally illustrated by the predictions obtained in this work and for instance by Semenova et al. (2020) for the specific case of NGC 2420. Insight and validation of the different available prescriptions for the vertical turbulent shear transport from hydrodynamicists and multi-dimensional numerical simulations, similar to the ongoing work on the horizontal turbulent shear (Park et al. 2020(Park et al. , 2021, are now mandatory in order to overcome this impasse. Finally, two additional promising leads to understand the formation of the Li and Be dips, in particular, should be explored, which were neglected here. First, we do not consider radiative accelerations on heavy elements. Their effects are, however, nonnegligible for stars with an effective temperature higher than ≈ 6'800 K, which corresponds to ≈ 1.4 M ⊙ at solar metallicity, close to the hot edge of the Li dip (e.g. Richer & Michaud 1993;Deal et al. 2018Deal et al. , 2020. Their impact on the stellar opacity could partly modify some of our conclusions concerning the most massive models. Second, the process sustaining the parametric turbulence used in this work and others in the literature has yet to be identified. In particular, even if the tachocline mixing (Spiegel & Zahn 1992;Brun et al. 1999;Garaud 2020) was shown to not be adapted in solar-type stars (Paper I), this process could play a role in the building of the Li and Be dips because of its dependence at the latitudinal differential rotation that is predicted to scale inversely with rotation for F-type stars (Augustson et al. 2012) compared to G-type stars.