Convective core entrainment in 1D main sequence stellar models

3D hydrodynamics models of deep stellar convection exhibit turbulent entrainment at the convective-radiative boundary which follows the entrainment law, varying with boundary penetrability. We implement the entrainment law in the 1D Geneva stellar evolution code. We then calculate models between 1.5 and 60 M$_{\odot}$ at solar metallicity ($Z=0.014$) and compare them to previous generations of models and observations on the main sequence. The boundary penetrability, quantified by the bulk Richardson number, $Ri_{\mathrm{B}}$, varies with mass and to a smaller extent with time. The variation of $Ri_{\mathrm{B}}$ with mass is due to the mass dependence of typical convective velocities in the core and hence the luminosity of the star. The chemical gradient above the convective core dominates the variation of $Ri_{\mathrm{B}}$ with time. An entrainment law method can therefore explain the apparent mass dependence of convective boundary mixing through $Ri_{\mathrm{B}}$. New models including entrainment can better reproduce the mass dependence of the main sequence width using entrainment law parameters $A \sim 2 \times 10^{-4}$ and $n=1$. We compare these empirically constrained values to the results of 3D hydrodynamics simulations and discuss implications.


INTRODUCTION
It has long been known that convective boundary mixing (CBM) must be included into stellar models in order to reproduce observations. The main sequence (MS) width of clusters is one of the best-known examples of such observations; other examples include large samples of wide binaries and asteroseismic measurements (e. g. Claret & Torres 2019;Deheuvels et al. 2016). As a result, stellar models' CBM schemes are calibrated to give results consistent with the observed reality. Castro et al. (2014) showed that current generations of models have MS widths on the Hertzsprung-Russell diagram (HRD) which are too narrow for high mass stars. The discrepancy in width grows larger with mass.
Currently, CBM is usually implemented in 1D stellar evolution codes in one of two ways. The first of these is step overshoot, which is an extension of the convective core by E-mail: l.j.a.scott@keele.ac.uk some fraction of a pressure scale height (see e. g. Ekström et al. 2012). Depending on the code, mixing in the overshoot region could be instantaneous or diffusive. The second is exponentially decaying diffusion, where mixing is governed by a diffusion coefficient which decays exponentially from a value near the Schwarzschild boundary (Freytag et al. 1996;Herwig 2000). The parameters in both can be calibrated in order to match observations such as post-MS spin down ( Fig. 1 in Brott et al. 2011) and asteroseismic frequencies (Aerts et al. 2018).
Both 3D hydrodynamics simulations and observations can be compared to 1D models incorporating CBM and other mixing processes, such as waves (Meakin & Arnett 2007;Jones et al. 2017;Edelmann et al. 2019;Müller 2020;Pratt et al. 2020). Incorporating 3D hydrodynamics results, such as convective regions which grow as a result of entrainment, into 1D models also allows them to be studied on evolutionary timescales. It has been shown that the rate of entrainment of material at convective borders is dependent on the bulk Richardson number, Ri B , a dimensionless mea-sure of the penetrability of a boundary by convection. For example, Cristini et al. (2019) showed that both the upper and lower boundaries in a convective shell followed the same entrainment law, suggesting that CBM is controlled by the global properties of the convective region. Despite these results from simulation, the entrainment law is not widely used in 1D stellar evolution codes.
Prior to this study, only Staritsin (2013) has published 1D entrainment law stellar models. Staritsin's models of 16 and 24 M main sequence stars with entrainment were calibrated using asteroseismology values for the extent of mixing. In these models, the extent of extra mixing beyond the formally convective region (in units of pressure scale heights) decreased as the models evolved. This contrasts traditional CBM which typically stays constant.
In this paper, we investigate entrainment in 1D main sequence models from 1.5 to 60 M using the Geneva stellar evolution code. We then compare our new models with entrainment to models including standard overshoot and constrain our entrainment parameters using these models. We further constrain our entrainment models with comparison to observations, in particular the MS width. We contrast entrainment parameters constrained by observations with those obtained in 3D simulations and discuss implications. Section 2 explains the definition and calculation of Ri B and the entrainment algorithm, along with the parameters of the model grid. Section 3 discusses the properties of the models, focussing on Ri B and the entrainment parameters. We present our conclusions in Section 4 and discuss the the plausibility and implications of our entrainment algorithm.

Calculation of Bulk Richardson number
The bulk Richardson number is defined as where l c is a length-scale for turbulent motions in the convective region, v c is the typical speed of convective flows and ∆b is the buoyancy jump. For l c we use 0.5H P,b , where H P,b is the pressure scale height at r = r b and r b is the radius of the convective core. If there is no CBM included in the model, r b is equivalent to the Schwarzschild boundary. Otherwise, it is the radius to which the CBM extends. l c represents the length-scale of the largest fluid elements in the turbulent region. The motivation for our choice of l c = 0.5H P.b comes from the results of Meakin & Arnett (2007), who found that the horizontal correlation lengthscale for velocity in their simulation of convection was approximately half a pressure scale height. We aim to be consistent with Cristini et al. (2019) by using this estimate of the horizontal correlation length-scale as a proxy for l c .
The buoyancy jump is an integral of the squared buoyancy frequency N 2 with respect to radius r, given by where r 1 and r 2 encompass the boundary of the convective core, centred at r = r b . The upper limit r 2 is equal to r b plus some fraction of a pressure scale height; in this study we use r 2 = r b + 0.25H P,b to be consistent with previous work (Cristini et al. 2019). Conversely, r 1 is the larger of either r b − 0.25H P,b or the Schwarzschild boundary. Using this maximum prevents negative N 2 regions, which do not contribute to buoyancy braking, from being included in the buoyancy jump integration. In our prescription, the size of the integration region encompassed by r 1 and r 2 is between 0.25 and 0.5H P,b , depending on the size of the entrainment region. This is supposed to encompass the part of the boundary in which fluid elements are decelerated and turned back towards the convective region by buoyancy. Cristini et al. (2019) Table 2 gives examples of their simulation boundary widths, which are all a fraction (0.1 to 0.6) of a pressure scale height. We cannot use a boundary width of ∼ 2v c /N as in Staritsin (2013), since at our boundary we have N = 0, so we use the approach of Cristini et al. (2019). However, we cannot be sure if the boundary width does not vary with mass (other than what is already contained in the mass dependence of H P,b ), and the integration region size must still be considered a free parameter. The buoyancy jump and its dependence on these parameters is discussed in more detail in Section 3.1 and Fig. 2.
For the buoyancy frequency, we use where g is the gravitational acceleration, H P the pressure scale height, δ the density gradient with respect to temperature (− ∂ ln ρ ∂ ln T ), ∇ ad the adiabatic temperature gradient, ∇ the actual temperature gradient, and ∇ µ the mean molecular weight gradient.
For v c , we use a mass-weighted root mean square of the mixing length theory (MLT) velocity, v MLT , in the core: where i represents the model mesh point and ∆m is the mass contained between the midpoints of shells i + 1 and i, and i − 1 and i. This sum over i is taken from the centre of the core to the Schwarzschild boundary. This average value is less subject to fluctuation due to numerical factors such as zoning than a single value chosen at some distance from the boundary. v c estimates the typical speed of the flow in the convective region which is responsible for entrainment.

Entrainment law algorithm
Ri B is a good measure of how difficult it is for convective flows to entrain material from the stable region. Indeed, the numerator, l c ∆b, measures the stability or stiffness of the convective boundary region via the buoyancy frequency N 2 . The denominator, v 2 c (∝ specific kinetic energy), measures the vigour of the convective flows approaching the boundary. A higher Ri B value thus means that it is harder for convection to entrain material from the stable region above.
We then use Ri B in the entrainment law 1 (e. g. Fernando 1991) to calculate an entrainment rate. The entrainment law where v e is the convective boundary progression speed and A and n are parameters controlling the entrainment rate. Note that if n = 1 (as is the case for most of our models), any uncertainty in l c in Eq. 1 would inversely scale A. However, we are not targeting exact values for these parameters, and we can be fairly certain given the results of Meakin & Arnett (2007) that l c ∼ H P,b as we have assumed. A mass entrainment rate,Ṁ ent , can be derived from this to giveṀ with ρ b being the density at r = r b . The mass contained within the entrained region, M ent , is then where j denotes the model time step with length ∆t. This region is then considered part of the convective core. This means that the region is then instantaneously mixed and the temperature gradient is set to ∇ ad (further discussed in Section 4). In our implementation of the entrainment law, the entrained mass accumulates over the lifetime of the core with each time step, according to Eq. 7. Since the value of Ri B con-trolsṀ ent, j rather than M ent directly, any previous history of entrainment in the models is unaffected by the instantaneous value of Ri B . This contrasts the previous implementation of Staritsin (2013), in which the entrained distance at any time step is equal to v e ∆t. Thus, our prescription can be viewed as cumulative entrainment and Staritsin's as instantaneous (meaning that it depends only on the stellar structure at the current time step). 3D hydrodynamic simulations of stellar convection which exhibit entrainment show that the convective region continuously accumulates material. This is the motivation for a cumulative entrainment method, as once material is entrained, it stays well-mixed. However, it is not known whether this holds true on evolutionary time-scales so it is not clear at this point which approach is more appropriate. Our method allows us to investigate the consequences of cumulative entrainment which is controlled by the changing value of the bulk Richardson number and we compare our results to Staritsin (2013) in Section 4.

Geneva code model grid
We use the Geneva stellar evolution code (GENEC, Eggenberger et al. 2008) to compute a grid of non-rotating MS models with solar metallicity (Z = 0.014). The masses included are 1. 5, 2.5, 8, 15, 25, 32, 40 and 60 M . For each mass we compute at least one standard CBM model and one entrainment model.
The standard CBM prescription in GENEC is step overshoot, where the convective core is extended by some distance α ov H P,b . In GENEC, the default value for α ov is 0.1 for law. However, since this is the currently accepted terminology in other fields such as geophysics, we will continue to use it. models with initial mass M ini ≥ 1.7 M , 0.05 for 1.7 M > M ini ≥ 1.25 M and 0 for M ini < 1.25 M . These default values were calibrated using the MS width of low mass stars (for details see Ekström et al. 2012). As in the core, the CBM (a.k.a. overshoot) region is mixed instantaneously (for both chemical species and entropy) and uses the adiabatic temperature gradient. Table 1 lists the models computed and their key properties. The first four columns of Table 1 define the initial parameters of the model. These are the initial mass M ini and the CBM parameters (either α ov for step overshoot models or a combination of A and n for entrainment models). The τ MS column is the main sequence lifetime. This is defined as the age of the model when the central hydrogen mass fraction has reached 10 −4 . The next column, T eff,min , is the minimum effective temperature reached by the model during the MS. Next is the mean of the bulk Richardson number, Ri B , taken over the duration of τ MS , along with the means of its components, v c and l c ∆b . The final three columns pertain to the model attributes at the end of the MS. These include the final mass, M fin , the mass of the helium core, M He , and the total mass entrained, M ent,tot . M He is defined as the mass of the convective core at a central hydrogen mass fraction of one per cent.
Both a default step overshoot model and an entrainment model with A = 10 −4 and n = 1 were calculated for each mass. This value of A was chosen to reproduce the MS lifetime of the 2.5 M standard overshoot model, as this mass is within the mass range originally used to calibrate the step overshoot. A = 2 × 10 −4 was also used for some masses to explore the widening of the MS in the high-mass range.
Previous simulations of convection have found that n ∼ 1, which guided our choice to keep n = 1 for the majority of our grid. However, the A values used for our 1D MS models (A ∼ 10 −4 ) are substantially lower than those derived from 3D simulations. A values derived from 3D simulations include A = 1.06 (Meakin & Arnett 2007, oxygen burning), A ≈ 0.1 (Müller et al. 2016, oxygen burning) and A = 0.05 (Cristini et al. 2019, carbon burning). The difference could simply be a matter of evolutionary phase, since these 3D simulations are all of later stages than the MS. One potential confounding factor is radiative diffusion. Since the burning stages from carbon onward are neutrino-cooled, the effect of radiative diffusion on the mixing process is minimal, in contrast to the MS. Another point is partial degeneracy, which plays a part in later-stage stellar evolution but not in MS convective cores. Finally, the entrainment law may not keep the same slope for all Ri B values. Our 1D models have Ri B in the range of ∼10 4 to ∼10 7 , which is substantially higher than the upper limit of Ri B ∼ 1000 in the 3D simulations and may represent a different entrainment law regime. Alternatively, there may be other important physics which is not encompassed by the entrainment law in its current form.
See Section 3.3 for more details on the chosen entrainment parameter values. Appendix A contains details on model resolution.

Time dependence of boundary penetrability and mass entrainment rate
The time dependence of the bulk Richardson number, Ri B , for two 15 M models (step overshoot with α ov = 0.1, entrainment with A = 10 −4 and n = 1) is presented in Fig. 1.
Over the MS, the variations in Ri B are modest, within one order of magnitude. Nevertheless, we can see in Fig. 1 that Ri B initially increases and later on decreases.
The increase in Ri B can be understood by considering the evolution of the buoyancy jump ∆b (the length-scale for turbulent motions, l c , which is set to half of a pressure scale height, is roughly constant during the MS), which is an inte-gration of the buoyancy frequency N 2 over the boundary region. In a massive star such as the 15 M model plotted, the convective core continuously recedes in mass over the MS. As the convective core recedes, it leaves behind a chemical gradient which contributes to an increase in N 2 and hence ∆b. This leads to an increase in l c ∆b (the numerator in Ri B shown in the second row of Fig. 1), which is strongest at the very beginning of the main sequence since there is no chemical composition gradient to start with. After some time (age ∼ 6.5 Myr), the core recedes far enough that the outermost limit of the buoyancy jump integration is lower than the original extent of the convective core. From this point onward, ∆b remains roughly constant since the full extent of its integration region is already occupied by the chemical gra- dient left by the convective core. Note that this saturation would likely occur earlier in the evolution if the size of the integration region was smaller; see the text below Eq. 2 in Section 2.1. The transient spikes in Ri B also come from spikes in ∆b. These originate from the finite differencing used in the code, since the boundary lies between two grid points. Fortunately, they have no impact on the results since they cause a temporary decrease in the entrainment mass rate (bottom row in Fig. 1). The spikes can be further explained by considering the integration of the squared buoyancy frequency. The buoyancy frequency depends on the gradient of the mean molecular weight, ∇ µ (see Eq. 3), which becomes the dominant part of N 2 at the upper edge of the CBM region. In the absence of mixing above this edge, ∇ µ can experience large local spikes. This is reflected in the mean molecular weight µ as step-like features rather than a smooth profile, and can cause transient increases in Ri B . These perturbations in Ri B do not cause pathological changes in the core mass, which evolves smoothly (see Fig. 4). We intentionally did not include any shear mixing be-yond the entrained region to study the effects of entrainment without any additional extension the MS lifetime from other factors. Any additional mixing processes such as shear would make it difficult to determine how much the entrainment itself affects the MS width and lifetime. However, we know from 3D simulations that there is a shear layer, which will smooth composition and structure profiles and probably prevent these spikes in the models (Arnett & Moravveji 2017;Jones et al. 2017). This shear layer could be modelled using an exponentially decaying diffusion coefficient (exp-D hereinafter, Freytag et al. 1996;Herwig 2000) at the edge of the entrained region. Preliminary results suggest that a combination of entrainment and exp-D improves the transient spikes in Ri B seen in pure entrainment models. Since exp-D provides an extra source of CBM, smaller values of A might be needed in these combination models to produce the required MS widths. Figure 2 shows the buoyancy jump integration region at three stages of the evolution of the 15 M entrainment model with A = 2 × 10 −4 . The dashed lines represent the position of the edge of the entrained region at r = r b . The dotted lines represent the upper and lower limits of the integration. Both N 2 (blue) and ∆b (green) are plotted, with ∆b being the value obtained when integrating from the convective boundary to the corresponding radius on the x-axis. Hence, the value for ∆b at (r − r b )/H P,b = 0.25 is the value plotted in Fig. 1. Values for ∆b at higher radius would be obtained if the upper limit of the integration was larger.
Since the temperature gradient in the entrained region is adiabatic, N 2 is only positive in the stable region, which is the only region to contribute to the buoyancy jump in our current models. Fig. 2 also shows that the main contribution to the buoyancy jump is from the region close to r = r b . Outside our chosen integration region, ∆b remains at a similar order of magnitude. If the integration region is in fact significantly smaller than our chosen value, e.g. 0.05H P , then the buoyancy jump would also be significantly smaller.
Finally, the modest decrease in Ri B towards the end of the MS is due to the gradual increase in convective velocities (third row in Fig. 1). The increase in convective velocities is due to the luminosity of the star gradually increasing over the MS. Since velocity and luminosity are related by v 3 c ∝ L (Biermann 1932), convective velocities also increase over the MS. Compared to ∆b, however, the variation in v c is small, which explains why ∆b has the greatest effect on the overall changes of Ri B during the MS.
Over the MS, Ri B varies between a few tens of thousands and a few hundreds of thousands (excluding short spikes explained above). Using the entrainment law (Eq. 5) with A = 10 −4 and n = 1, this leads to mass entrainment rates between 10 −6.3 and 10 −7.5 M yr −1 in a 15 M model. The mass entrainment rate, which is inversely proportional to Ri B , first decreases during the first part of the MS and later on increases slightly. The mass entrainment rate in this model leads to a total entrained mass of 0.960 M (see last column of Table 1).

Mass dependence of boundary penetrability
Current observations seem to suggest that convective boundary mixing is mass dependent.  2014) performed a large study on Milky Way stars and found significant broadening of the MS at higher masses; we compare to this work in particular in Section 3.4. In this section, we explore if this dependence can be explained by the mass dependence of stellar structure and properties. It is well known that the luminosity has a strong mass dependence. For low-mass stars, the dependence is steep with L ∝ M 3 . For massive stars, it flattens and approaches a linear dependence with mass above about 20 M (see Fig. 6 in Yusof et al. 2013). A higher luminosity leads to higher convective velocities (v 3 c ∝ L, Biermann 1932). Since the bulk Richardson number, Ri B , contains a velocity term, it would also be expected to show mass dependence.
The left panel of Figure 3 shows the logarithm of the time average of two values: Ri B and v 2 c (sign reversed, since it is the denominator of Ri, and scaled by a constant value to fit on the same axis). This panel demonstrates that Ri B is mass dependent and its dependence is dominated by the velocity term. The right panel shows the velocity term compared to total luminosity (again scaled by a constant), demonstrating that the mass dependence of velocity is also very similar to that of luminosity, as expected from the mass luminosity relation. Conversely, the buoyancy jump term also plotted in the right-hand panel does not demonstrate mass dependence since its logarithm varies by less than 0.5 dex. Despite this, the buoyancy jump term does dominate the variation of mass entrainment rate with time ( Fig. 1) and so cannot be ignored when considering entrainment at the convective boundary. Note also that this only holds if our assumptions on the buoyancy jump integration region (see text below Eq. 2 in Section 2.1) are correct.
In this section, we showed that convective boundary properties have a clear mass dependence, which can be measured via Ri B . Next, we want to explore whether the entrainment law, which uses Ri B can provide the mass dependence of the convective boundary mixing needed to reproduce the observed MS width. We can already note that Ri B decreasing with initial mass will lead to higher entrainment rates for more massive stars, which goes in the right direction.

Dependence of entrainment on the entrainment law parameters
Both 3D simulations and theoretical studies determined various values for the entrainment law parameters A and n. From a theoretical energy balance argument, n should be 1 (Stevens & Lenschow 2001 In this study, we start by taking n = 1 and use published 1D GENEC evolution models with step overshoot and α ov = 0.1 to determine a value of A that would reproduce the published models. The value of α ov = 0.1 in GENEC models is constrained using the main sequence width for low/intermediate-mass stars (Ekström et al. 2012). The same value of α ov is then applied to all higher masses (at all metallicities) in the published grids of GENEC models. Therefore, 2.5 M models were used to constrain an A value in entrainment models that matches the general properties of the 2.5 M GENEC model with step overshoot and α ov = 0.1: MS width in the Hertzsprung-Russell (HR) diagram, core masses and MS lifetime. Table 1 and Fig. 4 show the comparison between en-trainment models with different values of A and the default overshoot model. They confirm that the minimum effective temperatures reached by the models with α ov = 0.1 and A = 10 −4 , n = 1 are very similar. Table 1 also indicates that the MS lifetimes are similar. Figure 4 shows the evolution of the convective core mass in 2.5 M entrainment models. The left-hand panel shows how the entrainment depends on the value of A with values of A ranging from zero (no CBM) to 3×10 −4 (all models with n = 1). As expected, a larger value of A leads to more entrainment and thus larger convective core masses and longer lifetimes. One point to note is that since entrainment rate is reduced if Ri B increases, the potential problem of the convective region quickly encompassing the whole star can be avoided. Indeed, as the entrainment extends further, the jump in composition and entropy at the boundary increases and makes the boundary stiffer, which makes it harder for additional entrainment. The use of the entrainment law thus provides an important feedback. This is best seen for the A = 3 × 10 −4 model (brown curve), where entrainment leads to core growth only during the first part of the MS. After a while, the entrained mass plateaus since the entrainment rate drops significantly and the convective regions shrinks in mass due to the Schwarzschild boundary receding as in the step overshoot models. Note that much larger values of A may still lead to the entire model becoming convective. Much larger values of A are not needed or supported by observations anyway as discussed in Section 3.4.
Keeping n = 1, the value A = 10 −4 provides the closest match to the default step overshoot model in terms of MS lifetime. We see, however, that the time evolution of the convective core is very different in entrainment models compared to the step overshoot model, as shown in the righthand panel of Fig. 4. The step overshoot model assumes that mixing is an instantaneous process (compared to the MS lifetime) and thus the convective core is significantly larger on the ZAMS in these models. On the other hand, entrainment is a time-dependent process (in that the size of the entrained region is dependent on the earlier entrainment history) and builds up over the entire MS, as shown in Eq. 7. The dashedred line indicates the Schwarzschild boundary in the A = 10 −4 model and shows how the entrained region (region between the solid and dashed red lines) grows in mass with time. This means that for a given lifetime, the entrainment models start with smaller and end with larger convective core masses compared to step overshoot models (see Table 1 and Section 3.5).
We also tested the dependence of entrainment on the value of n with various 15 M models with values of n = 0.9, 1.2 and 1.5 (keeping A = 10 −4 , see Table 1). The dependence on n is strong. Indeed, values of n slightly larger than 1 (1.2 or 1.5) strongly reduce the total entrained mass (by a factor of 10 or more, see last column of Table 1) and values of n slightly smaller than 1 (0.9) lead to significantly more entrainment (by a factor of more than 3). While the dependence on n and A are not independent, our models tend to show that n cannot be too far from 1. We will compare the values determined in this study to observational constraints and hydrodynamic simulations in the discussion.

Impact of entrainment on main sequence width
One of the main observational constraints on stellar models is the MS width. Castro et al. (2014) represents one of the most comprehensive study of MS width at solar metallicity. One of their key findings is that models using a massindependent value of step overshoot (Ekström et al. 2012;Brott et al. 2011) do not reproduce the observed MS width.
Instead, it appears that CBM must increase with initial mass. While the sample used in Castro et al. (2014) is far from complete, it is worth comparing our new entrainment  Table 1). The mixing schemes used are denoted by different coloured markers as in the legend. The one-off α ov = 0.7 model with 25 M is shown with the black clover symbol (see Section 4). models with models with various amounts of step overshoot and the MS width deducted by Castro et al. (2014). Castro et al. (2014) find that the MS generally extends to a lg (T eff ) ∼ 4.3 over a range of luminosities, which corresponds to stars in the mass range ∼ 10−20 M . Above 20 M , the MS does not seem to have a well-defined cool end and instead appears to extend down to very cool temperature. Figure 5 shows the minimum effective temperature, T eff,min reached on the MS by each model in the grid. T eff,min for the default step overshoot models is shown with blue disks and we can see that they indeed predict an MS cool edge that deviates from the observed lg (T eff ) ∼ 4.3 further as the mass of the model increases from 10 M upwards. We also see that these models do not predict the observed widening of the MS above 20 M . As discussed in the previous section, entrainment models with A = 10 −4 and n = 1 (red pluses) reproduce the main features of the default (α ov = 0.1) step overshoot models (MS lifetime and HRD tracks). Thus as expected, the T eff,min of the entrainment models is also hotter than the observed one for stars between 10 and 20 M . One difference appears for stars above 20 M with the entrainment models predicting a cooler edge for the 32 M model and a very cool edge for the 40 M . The 60 M models do not follow this trend because they experience strong mass loss towards the end of the MS, which keeps the models on the hot side of the HRD.
Increasing the value of A from 10 −4 to 2 × 10 −4 (purple crosses) provides a reasonable match to the observed MS edge at T eff,min ∼ 4.3. Indeed, 4.34 ≤ T eff,min ≤ 4.24 in the models between 8 and 25 M . Furthermore, the 32 M model now extends to very cool T eff . While the observational constraints are not very tight, the MS width for lower masses is slightly wider than observations so a larger value of A would not be favoured. The broader MS width can be reproduced with an increased value of α ov (e.g. with α ov = 0.5, green squares) but in this case, the MS width for lower mass stars would be too wide. The reason why the entrainment models have  Table 1 for the polynomial coefficients of the three TAMS lines used in this figure. The dots, pluses and crosses have been placed where the model reaches 90%, 95% and 99% of the MS lifetime respectively.
broader MS width for more massive stars is due to the mass dependence of Ri B discussed in Section 3.2, which is used in the entrainment law. This means that the entrainment law provides a partial physical explanation for the apparent mass dependence of the overshooting parameters and a way of providing a much better fit to the observations with a single value of the parameters A and n, which is harder for other CBM such as step overshoot or exponentially decaying diffusion coefficients. Whilst the fit to the TAMS edge could be improved, for example by varying n in addition to A, the usefulness of this approach would be limited. Other factors such as rotation, metallicity variations and different mass loss prescriptions could provide additional mass-dependent factors which are not included in these models. Castro et al. (2014) gathered observations of galactic stars and placed them on a spectroscopic HRD (sHRD) (Langer & Kudritzki 2014), in which they show the density of observed stars in each region of the HRD. Since the sample is incomplete and possibly biased (Vink et al. 2010;McEvoy et al. 2015), it is difficult to compare densities of stars across the HRD. Nevertheless, it is still interesting to determine what fraction of the MS lifetime models spend in a given location in the HRD. This is indicated in Fig. 6 (dots at 90% of the MS lifetime, pluses at at 95% and crosses at 99%). Figure 6 shows our models on an sHRD so that we can compare directly to these observations. We focus on the 8 to 32 M range, which encompasses the region of the sHRD in which there is a clear observed TAMS boundary (marked on Fig. 6 with the dash-dotted line). Also shown are the TAMS boundaries obtained by Castro et al. (2014) Castro et al. (2014) TAMS in the middle of the mass range, deviating at both the high-mass and low-mass extremes. This suggests that when using step overshoot to determine CBM, a massdependent α ov is needed in models to reproduce observations over a large mass range.
The entrainment law naturally accounts for the mass dependence of CBM through the mass dependence of Ri B (see Section 3.2). In Fig. 6, the entrainment models have a markedly different TAMS boundary shape (approximated by the positions of the cross markers) with an increased widening of the MS with increasing mass. In particular, the

Impact of entrainment on helium core masses and lifetimes
The type and degree of CBM also affects the mass of the helium core at the end of the MS. The size of this core, whilst not directly observable, has very important implications for post-MS evolution. The compactness and explodability of pre-supernova models is dependent on the post-MS structure, in which the helium core plays an important role driven by the conditions in the core, CBM parameters that produce large cores can mimic the results of more massive models with less CBM.
In Table 1, the helium core mass at the end of the MS, M He , is given in the penultimate column. We define M He as the mass of the convective core (including CBM) when the central hydrogen mass fraction drops to one per cent. This definition gives similar results to taking the mass coordinate at which the hydrogen mass fraction drops to one per cent at the last time step of the MS. Figure 7 shows both M He (left) and M He divided by its value in the default step overshoot model (right). As expected, the left panel shows that larger amounts of CBM produce larger core masses at the end of the MS. In absolute terms, this increase in core mass is greatest in the more massive stars. In the right-hand panel, the majority of CBM choices show the opposite trend, with a greater effect of CBM on relative core mass for the lower-mass models. This is particularly true for α ov = 0.5. In contrast, the A = 10 −4 models increase M He by ∼30 per cent across the mass range of the grid, except for the 1.5 M model which displays a milder change in core mass.
The value of entrainment that best produces the Castro et al. (2014) MS width in Fig. 6 is A = 2 × 10 −4 . As can be seen in the right-hand panel of Fig. 7, this value of A creates helium cores which are a factor of 1.6 to 1.8 larger than using default step overshoot models. The 20 M point (taken from Ekström et al. 2012) in the left-hand panel of Fig. 7 illustrates the implications of this; a 15 M model with A = 2 × 10 −4 has a similar helium core mass to an α ov = 0.1 model, which is 5 M more massive initially. This shift of at least 5 M has wide-ranging implications for massive star evolution and their fate. Examples include the upper mass limit of observed supernova progenitors (Smartt 2009)  and the mass range of black hole production (e. g. Chieffi & Limongi 2020). Figure 8 shows the MS lifetime, τ MS (column 5 in Table 1), relative to the default step overshoot case for various CBM parameters across the mass range of the grid. This figure shows similar trends to the right-hand panel of Fig. 7, but with more CBM producing longer lifetimes rather than larger cores. When comparing step overshoot models only, it is clear that the relative increase in lifetime is smaller for higher-mass stars. The entrainment models show more complicated non-monotonic behaviour. However, it is important to note that the relative effect of increasing CBM on lifetime is milder compared to the effect on helium core masses; the maximum relative increase in lifetime in Fig. 8 is nearly 15 per cent for the 2.5 M model, whereas M He is increased by a ∼80 per cent for the same model in Fig. 7.
The strong effect on core masses and more modest effect on lifetimes can be understood from the difference between step overshoot and entrainment discussed in Section 3.3 and highlighted in Fig. 4. Whilst the mass contained within the CBM region in the step overshoot models decreases with time, entrainment is a cumulative process, which builds up over the main sequence and thus leads to much larger final core masses.
Another important difference for the later evolution is the initial sizes of convective cores. The step overshoot model starts with a much larger core. This will leave an imprint on the structure of that part of the star, which will affect the behaviour of the intermediate convective zone (Kaiser et al. 2020).

DISCUSSION AND CONCLUSIONS
We have calculated a grid of 1D stellar models using the Geneva stellar evolution code with masses between 1.5 and 60 M and at solar metallicity (Z = 0.014). We have shown that the boundary penetrability by convective flows, quantified by the bulk Richardson number Ri B , decreases monotonically with increasing mass. This decrease is dominated by the increase in typical convective velocities due to the steep mass-luminosity relation for stars in the 1 to 20 M range. The boundary stiffness, l c ∆b, is nearly invariant with mass in the range studied.
Due to the decrease in Ri B with mass, models with entrainment experience a mass-dependent increase in mixing. This is reflected in a corresponding mass-dependent MS widening in the HRD. For our models, we find the value of A which best reproduces the observed MS widths of massive stars is A = 2 × 10 −4 , with n = 1. Note however that more extended samples are desired to place a very tight constraint on A and that the effects of rotation were not considered in this work (Martinet et al. in prep.).
The choice of temperature gradient in the entrained region is also an important factor in the implementation of entrainment. As explained in Section 2.2, we use ∇ = ∇ ad in the entrained region, since 3D simulations show that entropy is well mixed as the convective region grows (Cristini et al. 2017). 3D simulations also show a narrow boundary above the entrained region with a smooth chemical gradient rather than a switch from one µ to another; it is likely that the mixing of entropy is similarly slowed compared to the entrained region in this boundary. Indeed, asteroseismic observations support MS convective cores with a smooth ∇ profile in the CBM region (Arnett & Moravveji 2017).
In standard models, the global evolutionary effect of a slight change in ∇ in the CBM region is subtle, especially if the CBM region is small. In entrainment models, however, the size of the CBM region towards the end of the MS can be significant, especially with larger values of A. The choice of ∇ may also have a more important role in entrainment models due to its effect on the buoyancy jump, ∆b. In our current implementation, the CBM region has no contribution to ∆b whatsoever, since it is fully mixed (∇ µ = 0) and ∇ = ∇ ad . This means that the buoyancy frequency in the entrained region is 0. If the temperature gradient were to instead transition smoothly from ∇ ad to ∇ rad within the entrained region (as explored in Michielsen et al. 2019), there would be a contribution to ∆b from the entrained region. This contribution would grow larger as the entrained region grows in size, therefore providing more feedback slowing the entrainment for larger values of A. Consequently, these models would require larger values of A than models with ∇ = ∇ ad to produce the same MS width.
Since we have demonstrated in Section 3.2 that the mass dependence of Ri B is dominated by the change in typical convective velocities with mass, it is interesting to test whether a scaling based on v c could provide a simpler alternative to entrainment. This seems reasonable since the dependence of Ri B with mass is almost entirely controlled by v c , with l c ∆b staying nearly constant with mass. To constrain this scaling, we take the value of α ov = 0.5 for 15 M , as this most closely matches the Castro et al. (2014) observational lg (T eff,min ) ∼ 4.3. The scaled overshoot parameter for each mass, M ini , is then given by where α 15 M = 0.5 and v m (M ini ) is the average of the convective velocity to the power m over the MS of the model of initial mass M ini . Various scenarios support different values for m. According to Eq. 6, the mass entrainment rateṀ ent is proportional to v 3 c . If α ov in the step overshoot case most closely corresponds toṀ ent in entrainment models, then m = 3 is appropriate. However, m = 2 would be supported if α ov corresponds best to the penetrability of the boundary (Ri B is inversely proportional to v 2 c if n = 1). The m = 1 case of α ov ∝ v c would be similar to the findings of Denissenkov et al. (2019), who reported that the exp-D f parameter scales linearly with the cube root of the convective driving luminosity, or equivalently f ∝ v c . Figure 9 shows the predicted values of α ov according to Eq. 8 with m = 1, 2 and 3. The m = 2 and m = 3 cases quickly reach very high values of α ov above 15 M ; thus the y-axis scale is zoomed onto the lower α ov range. The values of α ov = 0.05 and α ov = 0.1 for 1.5 M and 2.5 M respectively have already been calibrated (Ekström et al. 2012), but are also underestimated by the m = 2 and m = 3 cases. Only m = 1 matches the already-known values for the lower mass range and does not produce extremely high values in the higher mass range.
Since the m = 1 case seems the most reasonable, we have provided a polynomial fit to this scaling, described in the caption of Fig. 9. We emphasise that this scaling should only be considered a temporary fix to the problem of massdependent CBM and behaviour of the mass range above 60 M is unknown. Whilst the m = 1 scaling supports previous findings (Denissenkov et al. 2019), the step overshoot values at M ini ≥ 30 M are already much larger than the value of α ov = 0.5 favoured by Higgins & Vink (2019). Eq. 8 also does not take the stiffness of the boundary into account. This may be less of a problem for MS cores, but convective shells which have two boundaries are known to have different stiffnesses for each and different entrainment rates according to the entrainment law (Cristini et al. 2019). In Step overshoot parameter α ov scaled using Eq.8. The polynomial fit to the m = 1 points uses the equation α ov (M ini ) = −0.00037867M 2 ini +0.03885918M ini −0.01237484. Previous studies such as Claret & Torres (2017) and Moravveji et al. (2016) show that similar results are obtained using an exp-D f parameter which is roughly a factor of 10 to 15 smaller than the equivalent step overshoot α ov . Therefore a fit for f (M ini ) would be roughly 1/10 to 1/15 of α ov (M ini ).
addition, the possible mass dependence of H P,b (which is used to determine the total overshooting distance, d ov = α ov H P,b ) should not be discounted, as it also contains information on the stellar structure near the boundary.
Nevertheless, we have calculated an additional model at 25 M with α ov = 0.7, which is approximately the value suggested by the m = 1 case of Eq. 8. This model can be found in Table 1 and in Figures 5, 7 and 8 represented by a black clover symbol. Fig. 5 in particular shows that this model produces a very broad MS with a minimum lg T eff ∼ 4, which is consistent with Castro et al. (2014).
The focus of this study is on entrainment at the convective core boundary during the MS, but many 3D simulations which resulted in entrainment were concerned with later evolutionary phases. The effects of entrainment in post-MS 1D models are unknown, but may be similar to that of other CBM with phenomena such as increased likelihood of convective shell mergers. In convective envelopes, the lengthscales and pressure stratification can be significantly different to convective cores. The relatively high importance of thermal diffusion may mean that entrainment is not a suitable CBM prescription in convective envelopes (Viallet et al. 2015).
Since our entrainment implementation is cumulative, it is interesting to compare our results to those of Staritsin (2013), who used instantaneous entrainment. Staritsin's values for A were also much smaller than the results of 3D simulations, with A = 4.425 × 10 −4 for the 16 M model and A = 4.054 × 10 −4 for the 24 M . This is not dissimilar to our value of 2 × 10 −4 , perhaps due to the similarity in calibration: Staritsin required that the entrained distance at the ZAMS was 0.1 H P , guided by asteroseismic results for the star HD 46202 (Briquet et al. 2011). We also based our initial estimate of A = 10 −4 on the MS lifetime of models with α ov = 0.1, as explained in Section 3.3.
However, there are also significant differences between our models and the models of Staritsin (2013). The key result of Staritsin (2013) was an entrainment region which decreased with time as the model evolved; we instead see the opposite, since the mass of our entrained region can only ever increase (by construction). As such, Staritsin's entrainment models produced less Helium overall than standard α ov = 0.1 models, whereas our estimations for Helium core sizes were much greater in the entrainment models (see Table 1). In addition, the buoyancy jump continuously increases in Staritsin (2013), whereas we see a plateau in the buoyancy jump near the middle of the MS (as explained in Section 3.1. This difference could be due to the buoyancy jump integration distance used by Staritsin, h ∼ 2v c /N. Since v c grows with time during the MS (in Staritsin's models as well as ours), the integration length h would similarly increase with time, potentially leading to the increase in ∆b.
To conclude, the entrainment law, through its dependence on the bulk Richardson number, produces models with a wider MS for high mass stars than standard models. In addition, the extension of the MS increases with mass, as required by observation. However, the value of the entrainment law A parameter required to produce the correct MS width for the lower mass stars in our grid is orders of magnitude smaller than the value derived from 3D simulations of convection in the later stages of stellar evolution. This value may change further if more aspects of convective boundary physics are included, such as shear. Although these models are not complete, they are an important step in the right direction since they show that convective boundary penetrability is a key part of the physics behind the mass dependence of CBM.
tions grant ST/K003267/1, and Durham University. DiRAC is part of the National E Infrastructure.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.

A2 Time Resolution
The time step in GENEC is controlled using the energy generation rate in the centre. It is generally set so that the MS is split in several thousand time steps. Figure A2 shows the effect on MS lifetime of changing time step length by a factor λ t in a spatially resolved (λ s = 1) 15 M model. The total number of time steps for each of these models is given in Table A2. As with the spatial resolution, we chose a temporal resolution based on the lifetimes of the step overshoot models shown in Fig. A2 before calculating our grid. However, the entrainment model lifetimes shown for comparison display a similar behaviour to the step overshoot models. The temporal resolution corresponding to λ t = 1 results in ∼ 10 4 time steps in most cases, although the 1.5 M models generally have around half the number of steps compared to the other masses.
This paper has been typeset from a T E X/L A T E X file prepared by the author.