© ESO, 2015

1. Introduction

Low-mass, main sequence stars, with their deep convective envelopes, are astrophysical laboratories that allow us to investigate weak mixing processes under stellar conditions. The small distance between the convection zone proper and the Li-burning layer in these stars makes the surface abundance of Li a sensitive indicator of any mixing process that may be at work just below the convection zone. The surface abundance of Be, which is burnt at a somewhat higher temperature, provides an additional constraint on the extent of this mixing.

The abundance of Li in the Sun is more than two orders of magnitude lower than in meteorites (Greenstein & Richardson 1951; Anders & Grevesse 1989; Asplund et al. 2009). Lithium abundances observed in open clusters (Herbig 1965; Danziger 1969; Zappala 1972) and solar analogues and twins (Baumann et al. 2010; Meléndez et al. 2010, 2014; Monroe et al. 2013) of different ages show that Li depletion takes place systematically on the main sequence, although sufficiently cool stars also deplete significant amounts of Li in the pre-main-sequence phase (Magazzu & Rebolo 1989; Martin et al. 1994; Sestito & Randich 2005). The situation is further complicated by the observed spread in the Li abundances in cool stars of the same age, composition, and effective temperature (Duncan & Jones 1983; Soderblom et al. 1993a,b; Garcia Lopez et al. 1994; Jones et al. 1999). At the same effective temperature, fast rotators are more Li-rich than the slow ones in the Pleiades and in α Per (Soderblom et al. 1993a; Balachandran et al. 2011). Lithium distributions in open clusters also contain a prominent feature in the temperature range between 6200 K and 7000 K, the so-called Li dip, in which the Li abundance drops to very low levels (Boesgaard & Tripicco 1986).

The abundance of Be is notoriously hard to measure. Its value in the Sun is equal to the meteoritic value within 0.1–0.2 dex (Balachandran & Bell 1998; Asplund 2004). Beryllium studies in open clusters have mostly been focused on stars hotter than 6000 K with the general conclusion that the Li dip stars are also Be-deficient (Boesgaard 1976; Stephens et al. 1997) and that there is a Be-Li correlation on the cool side of the dip (Boesgaard et al. 2001, 2004a,b). Cool (T_eff ≲ 6000 K), main-sequence field stars show significant Be depletion (Santos et al. 2004; Delgado Mena et al. 2012), although much smaller than that of Li.

The standard solar model predicts negligible Li and Be depletion on the main sequence (Schwarzschild et al. 1957; Weymann & Sears 1965; Pinsonneault 1997; Schlattl & Weiss 1999). To explain the observations quoted above, we have to include an additional, slow mixing process into the standard models of low-mass stars. It should mix the material from the convection zone proper down to the Li- and Be-burning layers and the mixing rate should quickly decrease with depth.

Several candidate processes have been proposed (see Chaboyer1998, for a more comprehensive review). The most popular is rotation-induced mixing (Zahn 1992; Charbonnel et al. 1994; Meynet & Maeder 1997; Maeder 1997), especially for its ability to naturally explain the observed Li spread in stars of the same mass and age by ascribing it to the wide range of rotation periods stars are born with. One would expect the depletion rate due to rotation-induced mixing to increase with decreasing rotation period (the opposite to what is observed, see above), but it is possible to reverse this trend by assuming that the mixing-length parameter α_MLT decreases with decreasing rotation period (Somers & Pinsonneault 2015).

Alternatively, the mixing might be caused by some form of convective overshooting, usually considered in the diffusion approximation suggested by Freytag et al. (1996). If this is the case, the overshooting parameter f_ov must take a different value on the pre-main sequence than on the main sequence (Schlattl & Weiss 1999). Mixing induced by internal gravity waves has also been employed to explain the Li- and Be-depletion patterns in low-mass stars, including the Li dip (Press 1981; Garcia Lopez & Spruit 1991; Montalban 1994; Schatzman 1996). In this paper we focus on the process we call “convective settling” proposed by Spruit (1997) and elaborated in Andrássy & Spruit (2013, hereafter Paper I).

Envelope convection zones are dominated by large-scale downflows spanning the whole convection zone (see Nordlund & Stein 1997; Trampedach et al. 2014). They are generated by the strong radiative cooling in the photosphere, which makes their initial entropy much lower than that of the nearly isentropic upflow. As they sink they merge, entrain mass from the hot upflow, and are heated by radiative diffusion on small scales. If a small fraction (≈10^-7, see Paper I) of the photospheric downflows retain a substantial entropy deficit until they have arrived at the base of the convection zone, they will continue sinking until each of them has settled on its level of neutral buoyancy, hence the name “convective settling”. The entropy the downflows start with is low enough for them to reach the Li-burning layer. In general, there will be a broad distribution of entropy contrasts at the base of the convection zone, spanning from the low values predicted by the mixing-length theory (MLT) up to the highest values the downflows start with in the photosphere. They will settle at a range of depths, and mass conservation will enforce an upflow carrying the Li- and Be-depleted material back to the convection zone, reducing the surface abundances.

The envelope convection problem is too difficult even for the state-of-the-art numerical simulations. Hence, the distribution mentioned above can only be parametrised. In Paper I, we show that a power-law parametrisation leads to a model that can explain the slow, main-sequence Li depletion in the Sun without changing the thermal stratification enough to come into conflict with the results of helioseismology. The simplified model presented in Paper I, however, is based on an approximate, non-evolving model of the Sun and it does not utilise the Be constraint. In this paper, we construct a model of the convective settling process that takes stellar evolution into account and applies both the Li and Be constraints. We also extend the computation to a range of stellar masses (from 0.8 M_⊙ to 1.2 M_⊙) in order to compare predictions of the model with abundance measurements in open clusters, solar twins, and field stars. We do not attempt to reproduce the Li/Be dip, which is probably due to a separate process (Garcia Lopez & Spruit 1991).

2. Model

2.1. Overview

We construct a two-component, kinematic model of the convective settling process. The first component is an ensemble of downflows leaving the convection zone and sinking towards the deeper layers. They are created by the rapid cooling in the photosphere, and the entrainment and heating processes change their distribution on their way through the convection zone. We do not model these processes. Instead, we parametrise their output by a distribution of a mass flow rate¹ in the downflows with respect to their entropy contrast and model how they settle below the convection zone. Each downflow settles at the point where its entropy equals that of its surroundings, i.e. when it becomes neutrally buoyant. The other component of the model is an upflow due to mass conservation, the strength of which at any given depth depends on the total settling rate below that depth. The upflow advects the Li- and Be-depleted gas back to the convection zone.

The distribution of the mass flow rate is a central element in our model. The downflows at the base of an envelope convection zone are characterised by a wide range of entropy contrasts δs with respect to the isentropic upflow. This range can be split into two parts. One part represents the well-mixed material and the fluctuations created locally by turbulent convection. Their typical amplitude (δs)_typ is very small owing to the high efficiency of convection deep in the stellar interior. Therefore, they settle in a very thin layer below the convection zone. Downflows with an entropy contrast significantly greater than (δs)_typ correspond to the incompletely mixed remnants of the photospheric downflows. They span several orders of magnitude in δs, up to (δs)_max, which is the maximum entropy contrast reached by the downflows just below the photosphere. The coolest of these downflows settle deep below the convection zone. We ignore the well-mixed component in our model and parametrise the distribution of the mass flow rate carried by the incompletely mixed downflows. We do not model the entrainment and heating processes explicitly; the mass flow rate as a function of the contrast value as parametrised in the model refers to the point where the downflow ultimately settles. We call this distribution the settling rate distribution.

We assume that the settling rate distribution is a power law defined over the range of downflow entropy contrasts (δs)_min ≤ δs ≤ (δs)_max. The minimal entropy contrast considered, (δs)_min, corresponds to the transition point to the range of the well-mixed downflows, thus we set (δs)_min = (δs)_typ. We estimate (δs)_typ using the MLT and extract (δs)_max from radiation-hydrodynamic simulations of stellar photospheres, which are readily available today (see Sect. 2.3).

In the settling paradigm, the slope β of the settling rate distribution should be a characteristic value resulting from the physics of the entrainment and heating processes in the convection zone and in the settling region. Therefore, we keep this value – albeit unknown – constant for all stars and all their evolutionary stages, and investigate what influence β has on the results.

The total mass flow rate Ṁ of the distribution is given by the mass flow rate that is leaving the photosphere, but it is considerably modified by the entrainment and heating processes in the convection zone and in the settling region. One could argue that the entropy of a downflow depends on the radiative heat exchange on small scales between the core of the downflow and the entrained material, which takes place on the downflow’s way from the photosphere to its settling point. The relative importance of radiation in the convection zone is parametrised through the value of (δs)_typ. The stronger the radiative heating of the downflows, the higher the value of (δs)_typ. On the other hand, the stronger the radiative heating, the smaller the fraction of the cold downflows that reach the settling region. To take this into account we scale Ṁ both in proportion to the mass downflow rate in the photosphere and in inverse proportion to (δs)_typ. We use this scaling so as to include the qualitative effect of the heating process, although we realise that more parameters are likely to play a role and the dependence is much more complicated in reality. The constant of proportionality in the scaling is adjusted until the model reproduces the observed Li depletion in the Sun for a given value of β.

The convective settling process, depending on its strength, can change the thermal stratification below the convection zone. In Paper I, we quantified this effect and showed it to be negligible in a solar model calibrated to reproduce the Li depletion observed in the solar photosphere if the input settling rate distribution is not too steep (β ≲ 2.5). In this work, we neglect the influence of convective settling on the thermal stratification and compute Li and Be depletion using a few pre-computed stellar-evolution models.

2.2. Stellar models

We consider the main-sequence evolution of five solar-composition stars with masses of 0.8 M_⊙, 0.9 M_⊙, 1.0 M_⊙, 1.1 M_⊙, and 1.2 M_⊙. The models have been computed with the stellar-evolution code GARSTEC (Weiss & Schlattl 2008), neglecting the processes of convective overshooting and gravitational settling (also called sedimentation or diffusion). Since the code does not output the stratification of the specific entropy s in the star, we compute an approximation to it from the stratification of the temperature T and pressure p using the ideal-gas expression $\mathit{s}\mathrm{=}R\ln \left(\frac{{\mathit{T}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\mathit{p}}\right)\mathit{,}$(1)where R= k_B/ (μm_u) is the gas constant including the mean molecular weight μ, k_B the Boltzmann constant, and m_u the atomic mass unit. The mean molecular weight is constant in the region we are interested in. Equation (1) assumes a constant level of ionisation, which is a good assumption at the temperatures and densities prevailing below the convection zones of the stars considered.

2.3. Mathematical formulation

2.3.1. Parametric model

The settling rate ṁ in the downflows is distributed as $\mathrm{d}\mathit{ṁ}\mathrm{=}\mathit{Ṁ}\mathit{f}\mathrm{\left(}\mathit{\delta s}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\delta s,}$(2)where Ṁ ≥ 0 is the total settling rate, and δs = s₀−s_d is the entropy contrast with s₀ the entropy of the stratification at the bottom of the convection zone and s_d the entropy of the downflow. As discussed in Sect. 2.1, the distribution function f(δs) is assumed to be a power law, $\mathit{f}\mathrm{\left(}\mathit{\delta s}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{for}\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\min }\mathrm{\le }\mathit{\delta s}\mathrm{\le }\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\max }\mathit{,}\\ & \mathrm{otherwise}\mathit{,}\end{array}$(3)where N is a normalisation factor, (δs)_min> 0 and (δs)_max> 0 are the bounds, and β> 0 is the slope of the distribution. We require ${{}^{\mathrm{\int }}}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathit{f}\mathrm{\left(}\mathit{\delta s}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\delta s}\mathrm{=}\mathrm{1}$, so that $\mathit{N}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{for}\mathit{\beta }\ne \mathrm{1}\mathit{,}\\ & \mathrm{for}\mathit{\beta }\mathrm{=}\mathrm{1.}\end{array}$(4)A downflow of entropy s_d = s₀−δs reaches neutral buoyancy and settles down at the point where the entropy of the surrounding stratification s equals s_d. Mass conservation requires the upward mass flow rate at this point to be ṀF(s), where F(s) is the fraction of downflows that settle below² this point and is given by the cumulative distribution function $\mathit{F}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathit{s}}{\mathrm{\int }}}\mathit{f}\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\mathit{.}$(5)

2.3.2. Parameter scaling

As described above, there are four parameters to be specified: (δs)_min, (δs)_max, Ṁ, and β. In Paper I, we used fixed values because we only modelled one star, the evolution of which was also neglected. Now, we intend to model a range of different stars and follow their evolution, hence we have to adapt the parameter values to the changing physical conditions in the convection zone. The only exception is β, which is held constant (see Sect. 2.1).

The entropy contrast of the coldest downflow in the ensemble, (δs)_max, is given by the maximum entropy contrast that the cooling process in the photosphere can create. This value reaches a well-defined maximum (in a time-averaged sense) just below the photosphere and can be extracted from radiation-hydrodynamic simulations of stellar photospheres. We use the grid of models computed by Magic et al. (2013) and approximate the dependence of (δs)_max on the effective temperature T_eff of the star and on its surface gravity log g by the fitting function $\log \hspace{0.17em}\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\max }\mathrm{=}{\mathit{a}}_{\mathrm{0}}\mathrm{+}{\mathit{a}}_{\mathrm{1}}\mathit{x}\mathrm{+}{\mathit{a}}_{\mathrm{2}}\mathit{y,}$(6)where x = (T_eff−5777)/1000, y = log g−4.44, a₀ = 8.164, a₁ = 0.491, and a₂ = −0.461 with cgs units assumed throughout.

We set the lowest entropy contrast considered, (δs)_min, equal to a typical entropy contrast (δs)_typ predicted by the MLT (see Sect. 2.1). We use the MLT formulation of Kippenhahn et al. (2012) with α_MLT = 1.65 to estimate the super-adiabatic temperature gradient Δ∇ = ∇−∇_ad at the point where the pressure p = p₀ e^{− 1/2}, where p₀ is the pressure at the bottom of the convection zone. We then estimate the entropy contrast an adiabatic convective element would reach after having overcome approximately one pressure scale height in such an environment, $\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\min }\mathrm{\equiv }\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\mathrm{typ}}\mathrm{=}\frac{R}{{\mathrm{\nabla }}_{\mathrm{ad}}}\mathrm{\Delta \nabla }\mathit{.}$(7)The total mass settling rate Ṁ is scaled as (see Sect. 2.1) $\mathit{Ṁ}\mathrm{=}\mathit{Ṁ}\mathrm{0}\hspace{0.17em}\frac{\mathit{Ṁ}\mathrm{phot}}{\mathrm{2.04}\mathrm{\times }{\mathrm{10}}^{\mathrm{21}}\hspace{0.17em}\mathrm{g}\hspace{0.17em}{\mathrm{s}}^{-1}}{\left(\frac{\mathrm{(}\mathit{\delta s}{\mathrm{)}}_{\mathrm{typ}}}{\mathrm{7.96}\mathrm{\times }{\mathrm{10}}^{\mathrm{1}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{g}}^{-1}\hspace{0.17em}{\mathrm{K}}^{-1}}\right)}^{-1}\mathit{,}$(8)where Ṁ_phot is the mass downflow rate at the point (close to the photosphere) where the downflows reach the maximum entropy contrast (δs)_max, and the constant Ṁ₀ is adjusted until the solar model reproduces the observed Li abundance in the Sun. The numbers in the denominators in Eq. (8) correspond to the current solar values. The downflow mass flux in the photosphere, ℱ_phot, is computed in the same way as (δs)_max (Eq. (6)), $\log \hspace{0.17em}{\mathrm{ℱ}}_{\mathrm{phot}}\mathrm{=}{\mathit{b}}_{\mathrm{0}}\mathrm{+}{\mathit{b}}_{\mathrm{1}}\mathit{x}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathit{y,}$(9)where b₀ = −1.475, b₁ = −0.239, and b₂ = 0.511 in cgs units. No extrapolation is needed when using Eqs. (6) or (9) with the stellar models considered in this work. The mass downflow rate Ṁ_phot is then $\mathit{Ṁ}\mathrm{phot}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}\hspace{0.17em}{\mathrm{ℱ}}_{\mathrm{phot}}\mathit{,}$(10)where R_∗ is the radius of the star.

2.3.3. Computational approach

We neglect the thermodynamic response of the star to the convective settling process. We do, however, estimate the convective flux that would be caused by this process to further constrain our model and check the plausibility of our assumptions (Sect. 3.4). All downflows considered pass through the bottom of the convection zone (as defined by the Schwarzschild criterion) and their entropy contrast with respect to their surroundings is the highest at that point, hence the convective flux reaches a maximum there. Its value relative to the total flux of energy is estimated to be $\mathrm{ℱ̂}\mathrm{conv}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{p}}\mathit{T}\overline{\frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}}\mathit{Ṁ}}{\mathit{L}}\mathit{,}$(11)where c_p is the heat capacity at constant pressure, $\overline{\mathrm{\Delta }\mathit{T}\mathit{/}\mathit{T}}$ the mean temperature contrast in the distribution (weighted by the mass flow rate), and L the luminosity of the star, all evaluated at the bottom of the convection zone. As explained in Sect. 2.1, the mass flow rate for every downflow in our distribution corresponds to the point where the downflow settles, because we do not model mass entrainment explicitly. Most of the convective flux, however, is carried by the downflows that settle very close to the bottom of the convection zone, i.e. to the place where we compute . Hence, we regard Eq. (11) as a reasonable order-of-magnitude estimate.

To compute the Li and Be burning, we map the GARSTEC models (Sect. 2.2) on a grid equidistant in the mass fraction q³. We only include the outermost 10 –15% of the stellar mass, which are relevant for the convective settling process. Interpolation of the models in time is done via the nearest-neighbour algorithm. We model the burning and transport of Li and Be using the set of equations $\frac{\mathrm{d}{\mathit{A}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{R}}_{\mathrm{b}\mathit{,i}}\mathrm{+}{\mathit{R}}_{\mathrm{m}\mathit{,i}}\mathrm{+}{\mathit{R}}_{\mathrm{s}\mathit{,i}}\mathrm{+}{\mathit{R}}_{\mathrm{a}\mathit{,i}}\mathit{,}$(12)where A_i = N_i/N_H is the abundance of Li or Be in the ith grid cell, N_i the number of Li or Be nuclei in the ith grid cell, N_H the number of hydrogen nuclei per grid cell⁴, i increases with radius, t is the time, R_b,i the burning rate, R_m,i a mixing rate, R_s,i the settling rate, and R_a,i the advection rate due to the upflow. The burning rate ${\mathit{R}}_{\mathrm{b}\mathit{,i}}\mathrm{=}\mathrm{-}\frac{{\mathit{A}}_{\mathit{i}}}{{\mathit{\tau }}_{\mathrm{b}\mathit{,i}}}$(13)is related to the nuclear-burning time scale τ_b,i, which we compute using the standard expressions for low-energy nuclear reaction rates in an ideal gas that can be found in e.g. Hansen & Kawaler (1994). We consider the burning of ⁷Li by the reaction ⁷Li(p,α)α and the burning of ⁹Be by the reactions ⁹Be(p,α)⁶Li and ⁹Be(p,d)⁸Be.⁵ The products of the last two reactions are quickly transformed to ³He and ⁴He nuclei, respectively, and are of no interest for this work. The low-energy astrophysical S-factors of the three reactions are taken from the NACRE-II database (Xu et al. 2013). Electron screening is neglected. The burning time scale below the settling layer is set equal to that at the bottom of the layer for numerical reasons. The rate at which the downflows settle in the ith grid cell is $\mathit{ṁ}\mathit{i}\mathrm{=}\mathit{Ṁ}\left[\mathit{F}\mathrm{\left(}{\mathit{s}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{F}\mathrm{\left(}{\mathit{s}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right)}\right]\mathit{,}$(14)where the cumulative distribution function F(s) is given by Eq. (5), s is the local entropy of the stratification, and the index i + ¹/ 2 refers to the top and the index i−¹/ 2 to the bottom of the ith grid cell. Every single cell of our equidistant grid contains a mass of Δm, so we can define a “recycling” time scale${\mathit{\tau }}_{\mathrm{r}\mathit{,i}}\mathrm{=}\frac{\mathrm{\Delta }\mathit{m}}{\mathit{ṁ}\mathit{i}}\mathit{,}$(15)which is the time it takes the convective settling process to completely replace the content of the ith grid cell by “fresh” material from the convection zone. Most of the downflows settle just below the convection zone and τ_r,i can become very short there. To avoid severe time-step restrictions, we define a well-mixed zone, which is composed of the convection zone and the region where τ_r< 10⁵ yr and τ_b> 10⁸ yr. The fast recycling together with the slow burning prevent the formation of any significant gradients in the abundances of Li and Be in the well-mixed zone. We homogenise this zone using the artificial mixing term ${\mathit{R}}_{\mathrm{m}\mathit{,i}}\mathrm{=}\frac{{\overline{\mathit{A}}}_{\mathrm{wmz}}\mathrm{-}{\mathit{A}}_{\mathit{i}}}{\mathrm{\Delta }\mathit{t}}\mathit{,}$(16)where ${\overline{\mathit{A}}}_{\mathrm{wmz}}$ is the average abundance of Li or Be in the well-mixed zone and Δt is the length of the current time step. We set R_m,i = 0 outside the well-mixed zone. The downflows in our model bring the Li- and Be-rich material from the convection zone and deposit it in the settling layer. The settling rate of Li and Be nuclei is $\mathit{Ṅ}\mathrm{s}\mathit{,i}\mathrm{=}{\mathit{A}}_{\mathrm{cz}}\frac{\mathit{X}\mathit{ṁ}\mathit{i}}{{\mathit{m}}_{\mathrm{p}}}\mathit{,}$(17)where A_cz is the abundance of Li or Be in the convection zone, X the hydrogen mass fraction, and m_p the proton mass. The rate of change of the abundance due to settling is then $\mathit{Ṙ}\mathrm{s}\mathit{,i}\mathrm{=}\frac{\mathit{Ṅ}\mathrm{s}\mathit{,i}}{{\mathit{N}}_{\mathrm{H}}}\mathrm{=}{\mathit{A}}_{\mathrm{cz}}\frac{\mathit{ṁ}\mathit{i}}{\mathrm{\Delta }\mathit{m}}\mathrm{=}\frac{{\mathit{A}}_{\mathrm{cz}}}{{\mathit{\tau }}_{\mathrm{r}\mathit{,i}}}\mathrm{·}$(18)Finally, there is an upflow due to mass conservation, which we model by the advection rate of Li or Be nuclei $\mathit{Ṅ}\mathrm{a}\mathit{,i}\mathrm{=}{\mathit{A}}_{\mathit{i}\mathrm{-}\mathrm{1}}\frac{\mathit{X}{\mathit{\sigma }}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{m}}_{\mathrm{p}}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\frac{\mathit{X}{\mathit{\sigma }}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{m}}_{\mathrm{p}}}\mathit{,}$(19)where ${\mathit{\sigma }}_{\mathit{i}\mathrm{-}{}^{\mathrm{1}}\mathit{/}\mathrm{2}}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathit{i}\mathrm{-}\mathrm{1}}\mathit{ṁ}{}_{\mathit{k}}$ is the mass inflow rate to the ith grid cell at its bottom boundary and ${\mathit{\sigma }}_{\mathit{i}\mathrm{+}{}^{\mathrm{1}}\mathit{/}\mathrm{2}}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathit{i}}\mathit{ṁ}{}_{\mathit{k}}$ is the mass outflow rate from the ith grid cell at its top boundary. The rate of change of the abundance due to advection is then ${\mathit{R}}_{\mathrm{a}\mathit{,i}}\mathrm{=}{\mathit{A}}_{\mathit{i}\mathrm{-}\mathrm{1}}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathit{i}\mathrm{-}\mathrm{1}}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{r}\mathit{,}\mathrm{i}}}\mathrm{-}{\mathit{A}}_{\mathit{i}}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathit{i}}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{r}\mathit{,}\mathrm{i}}}\mathrm{·}$(20)Equation (12) is independent of the absolute abundance scale (see also Eqs. (13), (16), (18), and (20)). Therefore, we start all our calculations with A_i = 1 and integrate Eq. (12) using the standard, fourth-order Runge-Kutta method. The initial condition that we want to impose (see Sect. 3.1) is then taken care of simply by rescaling the results accordingly. We use the usual astrophysical notation log ϵ = 12 + log (A_cz) in the rest of the text.

3. Results

3.1. Preliminaries

Cool stars deplete significant amounts of their Li during their pre-main-sequence (PMS) evolution. Changes in the structure of a PMS star, however, are rather dramatic and show very large scatter even in a single cluster (Kenyon et al. 1990; Kenyon & Hartmann 1995; Dunham et al. 2008; Baraffe et al. 2009; Evans et al. 2009). Applying our simple scaling relations to such a wide range of conditions would be questionable. Instead, we start all our computations at the age of the Pleiades, for which we adopt a value of 100 Myr. We assume that the observed Li distribution in this cluster is a reasonable approximation to the PMS Li depletion. Ignoring the scatter at T_eff = const., we model the mean trend by a smooth function (see Fig. 1) and use it as an initial condition for the Li abundance. The depletion of Be is much lower and difficult to measure at low T_eff, so we use the meteoritic value log ϵ_Be = 1.30 (Asplund et al. 2009) as an initial condition for the Be abundance. We assume that both Li and Be are homogeneously distributed in the interior of the star at the start of the computation.

Fig. 1 Lithium distribution in the Pleiades (Sestito & Randich 2005). The mean trend is approximated by a smooth function, which provides an initial condition for our model.

Fig. 2 Age dependence of (δs)max in the stellar models used in this work. The discontinuities are caused by the nearest-neighbour interpolation that we use.

The dependences of (δs)_max, (δs)_min, and Ṁ on stellar mass and age, as defined by Eqs. (6)–(8), are shown in Figs. 2–4. The values of Ṁ_phot in Eq. (8) turn out to lie within ~30% of one another for all of the stars considered, although the spread in the values of ℱ_phot in Eq. (10) is an order of magnitude larger. The scaling of Ṁ in Eq. (8) is thus dominated by the (δs)_typ factor.

Fig. 3 As in Fig. 2, but (δs)min is plotted.

Fig. 4 As in Fig. 3, but Ṁ is plotted. This figure assumes β = 2.0. The curves would be shifted downwards by 2.8 dex at β = 1.5 and upwards by 2.9 dex at β = 2.5.

3.2. Li and Be depletion in the Sun

The influence of stellar evolution on the Li depletion in the Sun is illustrated in Fig. 5. The depletion rate becomes quasi-stationary after ~200 Myr and slowly decreases as the Sun ages. The depletion rate is hardly influenced by the assumed value of β. The observational data over-plotted in Fig. 5 suggest a somewhat more pronounced slowdown in the depletion rate, although the error bars are quite large. We also show a non-evolving model, in which the stratification is given by the solar-structure model at an age of 4.6 Gyr and is not allowed to evolve during the computation. In this case, the depletion rate becomes constant after the initial transition, as could be expected.

The abundance of Be predicted by the evolving, 1.0 M_⊙ model at an age of 4.6 Gyr ranges from 1.15 at β = 1.5 to 1.17 at β = 2.5, which deviates from the observed value of 1.38 ± 0.09 (Asplund et al. 2009) by −2.5σ. The meteoritic value is only 1.30 ± 0.03 (Asplund et al. 2009).

Fig. 5 Age dependence of the Li depletion in the Sun computed using our parametric model of convective settling compared with the Li abundances in four solar twins and in solar-type stars of three open clusters (the Sun is a calibration point and the Pleiades set the initial condition for Li in our model; see Sects. 2.3.2 and 3.1). The thermal stratification in the model plotted by the dotted line is not allowed to evolve in time. The abundances in the open clusters are from Sestito & Randich (2005) and correspond to the solar effective temperature at the age of the cluster; the error bars show the typical scatter in the data. The measurements in 18 Sco and HIP 102152 are from Monroe et al. (2013), and those in HIP 56948 and HIP 114328 from Meléndez et al. (2012, 2014).

Fig. 6 Dependence on the effective temperature of the Li depletion computed using our parametric model of convective settling compared with the Li abundances in three 600-Myr-old open clusters as determined by Sestito & Randich (2005). Almost no dependence on β can be seen at this age.

Fig. 7 As in Fig. 6, but the comparison is made with two 2-Gyr-old open clusters; the measurements are from Sestito & Randich (2005). A slight dependence on β becomes visible at this age.

3.3. Mass dependence of Li and Be depletion

Figures 6 and 7 show the dependence of the Li depletion on the effective temperature of the star (hence on its mass) at a fixed age as compared with Li abundances observed in open clusters. The metallicity of the stars may influence the extent of Li depletion, because the higher the metallicity, the higher the opacity and the deeper the convection zone. The metallicities of the clusters used in this work are summarised in Table 1. Our model underestimates the Li depletion at an age of 600 Myr (Fig. 6) independently of the value of β. This may stem from our assumption of a homogeneous Li distribution in the stellar interior at the start of the computation. The fit to the Hyades data (see Figs. 5 and 6) would also improve if the age of this cluster were not 600 Myr as we assume, but 950 Myr as Brandt & Huang (2015) recently suggested. The data are scarce at an age of 2 Gyr (Fig. 7), but the overall trend fits better. A weak dependence on β can be seen at high effective temperatures.

The depletion of Be is much smaller and much more difficult to measure than the depletion of Li. Therefore, we resort to a comparison with field stars in Fig. 8. We show models with β = 2.0 only, because their dependence on β is rather weak. The trend in the Be depletion is well reproduced assuming that the stars with T_eff ≲ 5500 K are older than ~5 Gyr, which is a reasonable assumption for cool field stars.

Fig. 8 Dependence on the effective temperature of the Be depletion computed using our parametric model of convective settling in a range of ages compared with the Be abundances measured in a sample of field stars from Santos et al. (2004). The fastest evolving 1.2 M⊙ star is computed up to 5 Gyr only.

3.4. Heat flux due to convective settling

The convective flux at the bottom of the convection zone, (Eq. (11)), unlike the depletion of Li and Be, is very sensitive to the assumed value of β, see Fig. 9. This comes about because the settling rate distribution spans several orders of magnitude in the downflow entropy contrast δs and the slope β of the distribution is constant over the whole range (see Sect. 2.1). The calibration of the 1 M_⊙ model to the observed solar Li depletion sets the total amount of material that has to settle in the Li-burning layer. Thus, the calibration fixes the tail of the distribution where the entropy contrast is high, of the order of (δs)_max. If we increase the value of β and recalibrate the model, the settling rate integrated over the Li-burning layer will not change much, but that just below the convection zone, where δs ≈ (δs)_min ≪ (δs)_max, will increase considerably. The main cause of the mass dependence of seen in Fig. 9 is the variation of (δs)_min between different stars. This parameter changes both directly, by changing the normalisation factor N of the settling rate distribution (Eq. (4)), and indirectly, via the scaling relation for Ṁ, see Eq. (8) and the related discussion in Sect. 2.3. Since the values of (δs)_max span a much narrower range than those of (δs)_min, the influence of this parameter on is correspondingly weaker. The differences in the thermal stratification between the stars considered are significantly larger than those due to their main-sequence evolution. The dependence of on the age of the star is therefore much weaker than on its mass (see Fig. 9). The models with approaching or even exceeding unity are in conflict with our assumption of no thermal feedback of convective settling on the star (see the last paragraph of Sect. 2.1), so their predictions should be taken with a grain of salt.

Fig. 9 Dependence on the stellar mass of the convective flux due to the convective settling process at the bottom of the convection zone. Three sets of models with β ∈ {1.5,2.0,2.5} are shown at ages of 1 Gyr (thin lines) and 5 Gyr (thick lines).

4. Summary and discussion

The process of convective settling is based on the idea that the envelope convection zones of low-mass, main sequence stars are dominated by large-scale downflows spanning the whole convection zone. If a tiny fraction of the low-entropy, photospheric downflows crosses the whole convection zone without having experienced much heating, their strongly negative buoyancy will make them sink and settle as deep as in the Li- and Be-burning layers. Mass conservation implies an upflow of the Li- and Be-depleted material back into the convection zone, reducing the photospheric abundances.

Building on the results of Paper I, we have explored the dependence of the Li and Be depletion on stellar mass in the range from 0.8 M_⊙ to 1.2 M_⊙. We assume that the mass settling rate is distributed with respect to the entropy contrast of the downflow as a power law. It spans a wide range of entropy contrasts from (δs)_min, given by the MLT at the bottom of the convection zone, to (δs)_max, given by the maximum entropy contrast the downflows reach just below the photosphere. The slope β of the distribution, assumed to be a constant, parametrises the physics of entrainment and heating processes acting on the downflows on their way to their settling points. In the absence of information on the dependence of the total mass settling rate Ṁ on the structure of the star, we scale this parameter in proportion to the mass downflow rate in the photosphere and in inverse proportion to (δs)_min, so that we qualitatively capture the effect of radiative diffusion on the heating of downflows. The mass downflow rate in the photosphere turns out to be essentially constant for the stars considered. We calibrate the scaling of Ṁ so as to reproduce the solar Li depletion. We allow for stellar evolution, but neglect the thermal feedback of convective settling on the star because we have shown in Paper I that the feedback is negligible provided that β ≤ 2.5. The computation is started at 100 Myr, using the observed Li distribution in the Pleiades and the meteoritic Be abundance as initial conditions for Li and Be burning, respectively.

Changes in the solar structure cause a slowdown in the Li-depletion rate as the Sun ages. This slowdown seems to be somewhat milder in the model compared to the observed Li evolution in open clusters and solar twins. The main discrepancy occurs in the first ~1 Gyr when real stars deplete Li faster than those in our model. This may be a consequence of our assumption that both Li and Be are homogeneously distributed in the star when the computation is started. Nevertheless, it is encouraging that the model can nicely reproduce the observed dependence of the Be depletion on the effective temperature in a sample of field stars. The current abundance of Be in the Sun predicted by the model is also compatible with the observed value. These conclusions are essentially independent of the assumed value of β, which is likely caused by the similarity of the internal structures of the stars considered.

We show that the convective flux at the bottom of the convection zone is very sensitive to β. The mass-flow-rate distribution in our model does not include the downflows with an entropy contrast lower than the MLT-based estimate (δs)_min. Such downflows carry a significant portion of the total flux even close to the bottom of the convection zone. Thus, one would expect the convective flux in our model to fall roughly into the interval to be compatible with our conceptual picture of convective settling. This corresponds to 2.0 ≲ β ≲ 2.5 in the solar model, in agreement with the conclusions of Paper I.

The main caveat of our present analysis is the qualitative nature of the scaling of the total mass settling rate with the properties of the star. It determines the sensitivity of the predicted Li- and Be-depletion rates and of the convective flux on the stellar mass and age. We would need to model the details of the downflows’ mass entrainment and heating to shed light on this issue, and to answer the question of whether the mass-settling-rate distribution is a power law in the first place.

Acknowledgments

We thank Achim Weiss for the GARSTEC models used in this work, Zazralt Magic for the grid of model photospheres, and the anonymous referee for critical comments that improved the introduction section and the overall presentation of the text.

References

All Tables

All Figures

	Fig. 1 Lithium distribution in the Pleiades (Sestito & Randich 2005). The mean trend is approximated by a smooth function, which provides an initial condition for our model.
In the text

	Fig. 2 Age dependence of (δs)max in the stellar models used in this work. The discontinuities are caused by the nearest-neighbour interpolation that we use.
In the text

	Fig. 3 As in Fig. 2, but (δs)min is plotted.
In the text

	Fig. 4 As in Fig. 3, but Ṁ is plotted. This figure assumes β = 2.0. The curves would be shifted downwards by 2.8 dex at β = 1.5 and upwards by 2.9 dex at β = 2.5.
In the text

Fig. 5 Age dependence of the Li depletion in the Sun computed using our parametric model of convective settling compared with the Li abundances in four solar twins and in solar-type stars of three open clusters (the Sun is a calibration point and the Pleiades set the initial condition for Li in our model; see Sects. 2.3.2 and 3.1). The thermal stratification in the model plotted by the dotted line is not allowed to evolve in time. The abundances in the open clusters are from Sestito & Randich (2005) and correspond to the solar effective temperature at the age of the cluster; the error bars show the typical scatter in the data. The measurements in 18 Sco and HIP 102152 are from Monroe et al. (2013), and those in HIP 56948 and HIP 114328 from Meléndez et al. (2012, 2014).

In the text

	Fig. 6 Dependence on the effective temperature of the Li depletion computed using our parametric model of convective settling compared with the Li abundances in three 600-Myr-old open clusters as determined by Sestito & Randich (2005). Almost no dependence on β can be seen at this age.
In the text

	Fig. 7 As in Fig. 6, but the comparison is made with two 2-Gyr-old open clusters; the measurements are from Sestito & Randich (2005). A slight dependence on β becomes visible at this age.
In the text

	Fig. 8 Dependence on the effective temperature of the Be depletion computed using our parametric model of convective settling in a range of ages compared with the Be abundances measured in a sample of field stars from Santos et al. (2004). The fastest evolving 1.2 M⊙ star is computed up to 5 Gyr only.
In the text

	Fig. 9 Dependence on the stellar mass of the convective flux due to the convective settling process at the bottom of the convection zone. Three sets of models with β ∈ {1.5,2.0,2.5} are shown at ages of 1 Gyr (thin lines) and 5 Gyr (thick lines).
In the text