© ESO 2022

1. Introduction

The transport of chemical elements inside stars plays a fundamental role in determining their structure and evolution. The distribution of chemical elements in stellar interiors is the result of complex interactions between macroscopic and microscopic processes, which are still not fully understood and therefore poorly modelled.

The proper modelling of the transport of chemical elements is becoming crucial for characterising stars, taking advantage of the high-quality asteroseismic data obtained by space missions such as Kepler (Borucki et al. 2010) or the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) mission, as well as for future such missions. As an example, the PLAnetary Transits and Oscillations of stars (PLATO) mission, which is an ESA space mission set to launch in 2026, aims to detect Earth-like planets in the habitable zone of Sun-like stars. In order to characterise the detected exoplanets, it is crucial to determine the properties of the host star (namely its mass, radius, and age) with high accuracy. The goal of this mission is to reach an accuracy of about or lower than 3% in radius and 10% in mass and age for the planets. This translates into the need of reaching accuracies of about or lower than 2% in radius, 15% in mass, and 10% in age for host stars similar to the Sun (Rauer et al. 2014). However, this error budget can easily be exceeded when some transport processes (e.g. atomic diffusion) are neglected, especially for the age (e.g. Nsamba et al. 2018; Deal et al. 2018; Aerts et al. 2018; Deal et al. 2020).

Atomic diffusion selectively affects chemical elements, and is a direct consequence of internal gradients in stars, such as pressure and temperature. Helioseismology has shown in the past decades that atomic diffusion must be included in solar models (Christensen-Dalsgaard et al. 1993). Other examples of cases in which atomic diffusion is important are in G-, F-, and A-type main-sequence stars, where not taking it into account has a significant effect on their structure and abundance profiles (e.g. Théado et al. 2009; Richer et al. 2000; Richard et al. 2001; Michaud et al. 2011; Deal et al. 2016, 2018). When correctly modelled, atomic diffusion can be used to infer the efficiency of competing macroscopic transport processes (Richard et al. 2005; Michaud et al. 2011; Verma & Silva Aguirre 2019; Deal et al. 2020; Semenova et al. 2020).

The actual diffusion velocity of an element mainly depends on the competition between two main forces (or accelerations): gravity, and radiative accelerations (g_rad). The latter is a transfer of momentum between photons and ions and counteracts gravitational settling. It is different for each element. Thus, diffusion has a direct effect on the abundance profiles in the star, which in turn leads to modifications in the Rosseland opacities.

Atomic diffusion is efficient in the outer layers of stars because the diffusion timescale is approximately proportional to the density of protons (Deal et al. 2018). Therefore, accurate predictions of surface abundances, which are necessary for the inference of stellar properties, also require taking into account atomic diffusion, including radiative accelerations (e.g. Richard et al. 2001; Michaud et al. 2011). Surface metallicity is another fundamental stellar parameter that is affected by atomic diffusion, and it is important for the study of exoplanetary systems because the probability of a star to have a detectable planetary system is linked to its metal content (e.g. Alecian & Michaud 2005).

The number of stars for which atomic diffusion (including radiative accelerations) has an important effect ranges from 33% up to 59% of the PLATO core program star sample (Deal et al. 2018). This means that for a large number of the stars observed by PLATO, an accurate modelling of atomic diffusion will be necessary to achieve the requirements of the mission.

Despite their importance, only a few evolution codes currently incorporate consistent computations of stellar models that include the complete treatment of atomic diffusion. The code Modules for Experiments in Stellar Astrophysics (hereafter MESA; Paxton et al. 2018, 2019) is one of them. The aim of this paper is to compare MESA models with models computed using the Montreal/Montpellier evolution code (Turcotte & Richer 1998; Richer et al. 2000) and the Code d’Evolution Stellaire Adaptatif et Modulaire (hereafter CESTAM, the “T” stands for transport; Morel & Lebreton 2008; Marques et al. 2013; Deal et al. 2018), all models including atomic diffusion with the contribution of the radiative accelerations.

The first objective is to show how MESA compares to the two other codes at a given input physics (these two other codes were already compared, giving very similar results, in Deal et al. 2018). The second objective is to provide a comprehensive assessment of the implications of using the different MESA atomic diffusion options dedicated to speed up these computations with some approximations, and pointing out those that may affect the resulting models most.

This paper is organised as follows: in Sects. 2 and 3, we present a brief introduction to atomic diffusion and the way in which MESA, CESTAM, and the Montreal/Montpellier code calculate it, as well as a direct comparison between them. Then, in Sect. 4, we present the comparison between the three codes. Section 5 addresses the impact of MESA atomic diffusion options that may have a significant effect on the resulting abundance profiles, which is then continued in Appendix A. Finally, Sect. 6 is dedicated to the conclusion.

2. Atomic diffusion

Atomic diffusion is calculated from first principles of physics, being a direct consequence of the internal gradients of pressure, temperature, and concentration in stars. This can be seen in the diffusion equation,

$$\begin{array}{c}\hfill \rho \frac{\partial X_i}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\rho D_{\mathrm{turb}}\frac{\partial X_i}{\partial r}\right]-\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\rho v_i\right]+A_i{m}_p{\displaystyle [\sum _j^{}(r_{\mathit{ji}}-r_{\mathit{ij}})],}\end{array}$$(1)

where X_i is the mass fraction of element i, v_i is its atomic diffusion velocity, ρ is the density in the considered layer, D_turb is the turbulent diffusion coefficient (i.e. the combined effect of the macroscopic transport processes, which cannot be simply treated as their sum due to their possible couplings and interactions), A_i is its atomic mass, m_p is the mass of a proton, and r_ij is the reaction rate of the reaction that transforms element i into j.

Additionally, the diffusion velocity of a trace element i is given by

$$\begin{array}{c}\hfill v_i=D_{\mathit{ip}}[-\frac{\partial X_i}{\partial r}+\frac{A_i{m}_{\mathrm{p}}}{\mathit{kT}}(g_{\mathrm{rad},i}-g)+\frac{(\overline{Z_i}+1)m_{p}g}{2kT}+{\kappa }_T\frac{\partial lnT}{\partial r}],\end{array}$$(2)

where D_ip is the diffusion coefficient of element i in a proton dominated plasma. The variable g is the local gravity, $\overline{Z_i}$ is the average charge (in proton charge units) of the element i (roughly equal to the charge of the dominant ion), k is the Boltzmann constant, T is the temperature, κ_T is the thermal diffusivity, and g_rad, i is the radiative acceleration on element i (Richer et al. 1998),

$$\begin{array}{c}\hfill g_{\mathrm{rad},i}=\frac{1}{4\pi r^2}\frac{L_r^{\mathrm{rad}}}{c}\frac{{\kappa }_R}{X_i}{\displaystyle {\int }_0^{\infty }\frac{{\kappa }_{u,i}}{{\kappa }_{u,\mathrm{total}}}\mathcal{P}(u)\mathrm{d}u,}\end{array}$$(3)

where $L_r^{\mathrm{rad}}/(4\pi r^{2}c)$ is the total radiative momentum flux at radius r, u ≡hν/(kT) is the dimensionless frequency variable, 𝒫(u) is the normalised blackbody flux, κ_u, i is the monochromatic opacity of element i, κ_u, total is the total monochromatic opacity, and κ_R is the Rosseland opacity, which can be obtained by integrating over the spectrum,

$$\begin{array}{c}\hfill \frac{1}{{\kappa }_R}={\displaystyle \int }\frac{1}{{\kappa }_{\nu }}\frac{\mathrm{d}{B}_{\nu }}{\mathrm{d}T}\mathrm{d}\nu /{\displaystyle \int }\frac{\mathrm{d}{B}_{\nu }}{\mathrm{d}T}\mathrm{d}\nu .\end{array}$$(4)

Therefore, the diffusion velocity of an element depends mainly on two opposite forces (or accelerations): gravity, which causes the element to migrate towards the centre of the star (gravitational settling), and radiative accelerations. Radiative accelerations result from a selective absorption by the various chemical species present in the star of part of the net outgoing momentum flux carried by photons. Usually, this absorption produces a net outward force on each absorbing species, which can exceed its weight. If radiative accelerations are present, both over- and underabundances (compared to the initial abundance) appear in the abundance profiles of the star, mostly as a result of the imbalance between radiative forces and gravity. The importance of radiative accelerations comes from the fact that they dominate atomic diffusion in many stars, especially in the hottest white dwarfs and in main-sequence stars with T_eff > 6000 K (Richer et al. 1998). Nevertheless, competing hydrodynamic processes (e.g. rotationally induced mixing; Deal et al. 2020) may moderate the expected effects of atomic diffusion.

3. Stellar evolution codes

There are different sets of equations to solve the diffusion equations and different approaches for calculating g_rad. Some give priority to either accuracy or computation time, while others try to find a compromise between them (for an overview, see Alecian 2018). In this section we present the methods used by three codes, MESA (r11701), Montreal/Montpellier and CESTAM, to treat atomic diffusion, including g_rad, before summarising the differences between them.

3.1. MESA

Atomic diffusion is calculated in MESA by solving the Burgers equations (Burgers 1969) using the Thoul et al. (1994) formalism, while radiative accelerations are calculated using the OPCDv3.3 package (Seaton 2005). It follows an opacity sampling method with opacity tables at fixed frequency grids. In order to compute the integral of Eq. (3), MESA includes the modifications of the original OPCD routines described in Hu et al. (2011). Radiative accelerations are computed for H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ca, Fe, and Ni (see Sect. 4.3 for the case of nickel). This method is less accurate than the direct use of atomic data, but it has the advantage of being faster to compute. MESA also offers a vast number of customizable options, including options for atomic diffusion. In Appendix A we focus on some of them. A more detailed description of the treatment of atomic diffusion in MESA can be found in Paxton et al. (2011, 2018).

3.2. Montreal/Montpellier code

In the Montreal/Montpellier, atomic diffusion is calculated by solving the full set of Burgers equations (Burgers 1969) for H, He, C, N, O, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Ti, Cr, Mn, Fe, and Ni. Radiative accelerations are also computed with an opacity sampling method using the OPAL monochromatic opacity table (Iglesias & Rogers 1996, not publicly available), with a fixed frequency grid. A more detailed description of the treatment of atomic diffusion in the Montreal/Montpellier code can be found in Turcotte & Richer (1998) and Richer et al. (1998).

3.3. CESTAM

In the CESTAM code, atomic diffusion is computed following the formalism presented in Michaud et al. (1993), which is an approximation of the Burgers equations, for H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ca, and Fe. Atomic diffusion can also be computed solving the full set of Burgers equations (Burgers 1969). The models used in this study are computed with the first option. Radiative accelerations are computed using the single-valued parameter (SVP) derived in LeBlanc & Alecian (2004), which is based on a simplified form of Eq. (3) obtained by separating the terms involving the atomic quantities from those describing the local plasma. The use of the resulting parametric equations allows for faster computations when compared to opacity sampling because it involves very short pre-calculated tables. The SVP method is only valid for stellar interiors. It has to be noted that the CESTAM models presented in this study take the latest version of the SVP method into account (Alecian & LeBlanc 2020). A more detailed description of the treatment of atomic diffusion in CESTAM can be found in Morel & Lebreton (2008) and Deal et al. (2018).

3.4. Differences between the three codes

The treatment of atomic diffusion is similar in the three codes, being based on the Burgers equations, with slight differences: MESA uses the Thoul et al. (1994) formalism and CESTAM the Michaud et al. (1993) formalism. It has already been shown by code comparisons within the CoRoT/ESTA program that these formalisms give very similar results (Lebreton et al. 2008). To compute radiative accelerations, MESA and the Montreal/Montpellier code both use the opacity sampling method, while CESTAM uses the SVP approximation. The SVP method is numerically faster, but introduces a maximum uncertainty of approximately 30% on the g_rad (Alecian & LeBlanc 2020). We show in the next sections that this uncertainty is generally lower.

In addition to the different computation methods, the codes also use different opacity data. MESA and the SVP approximation implemented in CESTAM use the Opacity Project (OP) data, while the Montreal/Montpellier code uses the OPAL data. Nevertheless, the difference between the OPAL and OP data is close in terms of Rosseland mean opacities (e.g. 5% at the bottom of the convective zone of the Sun, Seaton & Badnell 2004). Therefore, we expect the differences in the radiative acceleration profiles to rather come from the difference in the monochromatic opacity of each element (integral of Eq. (3)) than from the differences in the Rosseland mean opacities between OP and OPAL. We also expect the abundance profiles of the MESA models to be closer to those of CESTAM because radiative accelerations are computed from the same opacity data (OP).

4. Comparison of MESA models with other codes

In this section, we first present the physical input parameters of the models computed with the three evolution codes (MESA, Montreal/Montpellier, and CESTAM). Then, we compare their internal abundance profiles and surface abundances, before proceeding to analyse the particular case of nickel in more detail.

4.1. Model input physics

In order to perform our comparison, we selected models of two different masses (1.1 and 1.4 M_⊙). The mass of 1.1 M_⊙ was chosen because it develops a deeper surface convective zone (SCZ) during the main sequence, mainly to compare the effect of gravitational settling, while the 1.4 M_⊙ model allows us to make a direct comparison with the models in Deal et al. (2018). The surface convective zone of these models is small enough to make the effects of radiative accelerations noticeable, but not too small to prevent unrealistic abundance variations appearing in models that neglect known additional transport processes, such as rotation-induced mixing.

All the models used for this comparison, as well as their selected physical input parameters, are presented in Table 1:

Montreal/Montpellier models. The reference models for the comparisons are computed with the Montreal/Montpellier evolution code. The 1.4 M_⊙ model (B0) is the same as was used in the comparison with CESTAM model in Deal et al. (2018).

MESA models. MESA models are computed with the same initial chemical composition, metal mixture (Grevesse & Noels 1993), and atmospheres (Eddington) as the Montreal/Montpellier models. The main differences are in the monochromatic opacity tables, the equation of states, and the nuclear reaction rates. We also calibrated the value of α_MLT in order to obtain the position of the bottom of the surface convective zone as close as possible to the Montreal/Montpellier models. This step is necessary in order to compare the abundance profiles and distinguish the differences coming from the structure and coming from the opacity data used to compute radiative accelerations. The monochromatic opacities were considered for a temperature range [10⁴,10^6.3] K, using options 3, 4, 5, and 6 of Table A.1. The default model (AD) includes the default value for atomic diffusion options, which speeds up the computations, while the other models (A0 and A1) include an optimised setup for some of these options. These specific options are discussed in Sect. 5 and are detailed in Appendix A. The MESA model A0¹ (see Table 1) was used as a reference to test the impact of the different atomic diffusion options available in MESA.

CESTAM models. The CESTAM models are computed with a more recent version of the code than was used in the Deal et al. (2018) comparison. This version of the code includes the improved method described in Hui-Bon-Hoa (2021) to calculate the Rosseland mean opacity with the OPCD monochromatic opacity tables (Seaton 2005). It also includes the latest version of the SVP tables (Alecian & LeBlanc 2020). Similarly to the MESA models, α_MLT was adjusted, and the same temperature range was defined for the use of monochromatic opacities to compute the Rosseland mean opacity.

4.2. Comparison of abundance profiles

The abundance profiles of ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²³Na, ²⁴Mg, ²⁷Al, ²⁸Si, ³²S, ⁴⁰Ca, and ⁵⁶Fe are presented in Fig. 1 for the 1.4 M_⊙ models. As expected, the agreement between A0 and B0/C0 is very satisfactory. The depletion in ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²³Na, ²⁴Mg, and ⁴⁰Ca, as well as the accumulations of ²⁷Al and ⁵⁶Fe, are well reproduced by the A0 model. However, the profiles of ²⁸Si and ³²S in the A0 model are not consistent with the B0 model, while they are consistent with the C0 model (with larger differences in the case of ²⁷Al). As we explain below, this is mainly due to the differences in the monochromatic opacity data used by the codes (Table 1). The AD and A0 models show identical abundance profiles. This shows that when they are properly set up, the options allowing us to speed up atomic diffusion computations have no effect on the abundance profiles (see Appendix A for more details about these options).

Fig. 1. Abundance profiles according to log(ΔM/M*) (with ΔM being the mass above the radius r) obtained with MESA (Ad model, dot-dashed orange lines; A0 model, solid red lines), the Montreal/Montpellier code (B0 model, dashed blue lines) and CESTAM (C0 model, dotted green lines) at 420 Myr. The dot-dashed grey lines represent the abundance profiles for the AD model without radiative accelerations. The vertical dotted grey line represents the average position of the surface convective zone of the three models. The surface is towards the left part of the plots.

When comparing all the curves with the dot-dashed grey lines (model without radiative accelerations), we clearly see that the effect of radiative accelerations is different from one element to the next. For He, N, O, and Ne, in this specific case, radiative accelerations are negligible.

Figure 2 shows the difference in surface abundances between the Montreal/Montpellier model (B0) and the MESA (A0) and CESTAM (C0) models². The differences in temperature, density, gravity, and Rosseland mean opacity around the bottom of the surface convective zones are smaller than 3%, 4%, 1% and 1%, respectively. The differences in the gravity profiles (in direct competition with radiative accelerations) cannot induce large differences in the abundance profiles in these conditions (see Fig. C.1, for which the gravity profile of the three models overlap). It leads to surface abundance differences lower than 0.01 dex for the elements not affected by radiative accelerations (e.g. He, C, and O). For almost all of the elements in the MESA model, the differences are smaller, except for ²⁸Si and ³²S.

Fig. 2. Difference of surface abundances2 between the model B0 and the 2 other models (A0: red symbols, and C0: green symbols) shown in Fig. 1. The error bars of the C0 model are estimated using the 30% uncertainties on the radiative accelerations values obtained with the SVP method. The dotted grey lines show the ±0.01 dex interval.

For ²⁸Si and ³²S, this is the indication that these differences come from the monochromatic opacity data, affecting the integral in Eq. (3). The radiative acceleration profiles for the A0 and C0 models are close but slightly larger than gravity at the bottom of the SCZ, while this is the opposite for the B0 model (see Fig. C.1). This explains why these two elements are accumulated in the A0 and C0 models, while that is not the case with the B0 model. The comparison with the CESTAM model (C0) confirms the origin of the behaviour of ²⁸Si and ³²S, as both codes use the same opacity data to compute radiative accelerations and show the same behaviour. In the case of ³²S, this difference can be seen in Fig. 7 of Seaton & Badnell (2004), where its contribution to the Rosseland mean opacity at the temperature of the BSCZ is slightly higher with OP than OPAL opacities, which leads to slightly higher values of the radiative acceleration for this element.

For ²⁷Al, the surface abundance of the C0 model is different from that of A0 (by ≈0.03 dex) despite the fact A0 and B0 are close. In this case, the difference cannot come from the monochromatic opacity data. On the other hand, the radiative acceleration profiles shown in Fig. C.1 indicate that the C0 model is closer to the B0 models than to that of A0. Nevertheless, the ²⁷Al abundance of the C0 model is consistent with the A0 and B0 abundances when considering the 30% uncertainty of the SVP method. This behaviour of the ²⁷Al abundance for A0 and C0 may be a combination of the effect of evolution and the radiative acceleration computation. Moreover, most of the time, surface abundances predicted with the SVP approximation are very close to those of MESA. This indicates that the 30% uncertainties of this method (see Sect. 3.4) are often smaller, and it confirms the robustness of the method.

The abundance profiles are even closer for the 1.1 M_⊙ models (see Fig. B.2). In this case, gravitational settling dominates the transport (which is less subject to uncertainties than radiative accelerations), hence the smaller differences. Their average differences are smaller than 0.005 dex (Fig. B.1). This also confirms the origin of the behaviour of ²⁸Si and ³²S in the A0 and C0 models compared to that of B0.

4.3. Case of nickel

The nickel monochromatic opacity data available in the OPCD package are not determined similarly to the other elements, but are derived from those of iron. Figure 3 shows the comparison of the abundance profile of nickel between the A0 and B0 models. Even though the abundance profile pattern is similar, the surface abundances of the two models are significantly different (10% of the surface mass fraction). Similarly to the previous section, Fig. 2 shows that the differences in the predicted Ni surface abundance (≈0.07 dex) cannot be related to the structure differences, are more than two times larger than the effect of the difference in the opacity data we see for ²⁸Si and ³²S, and are larger than the difference of ²⁷Al in the C0 model. For this reason, the surface abundances of nickel predicted by MESA models, and more generally by the OPCD package, should be taken with caution. This is also the reason why nickel is currently not available in C0 and CESTAM models (treated as an element at equilibrium) using the SVP method (Alecian & LeBlanc 2020).

Fig. 3. Same legend as Fig. 1 for nickel.

5. Atomic diffusion options in MESA

In this section we briefly address the impact of atomic diffusion options dedicated to computational optimisation available in MESA, which may have a significant effect on the resulting abundance profiles unless they are carefully set up. Quantitative results, along with a more detailed explanation of the options and their effects on surface abundances, are presented in Appendix A.

Two of the most important options allow consideration of a region at the surface as a single cell for atomic diffusion calculations (options 1 and 2 of Table A.1). The size of this region is defined either by a mass fraction of the total mass or by a specific temperature. While this helps minimise the computation time, these values must be chosen thoughtfully. If this value is higher than the fraction of mass of the SCZ, errors may occur. In this case, the mixed region is artificially forced to be deeper, introducing an artificial and efficient full mixing in regions where it should not be present. This option is further analysed in Appendix A, as well as other options whose effect on the surface abundances can be significant.

6. Conclusion

We computed stellar models with the MESA evolution code taking atomic diffusion into account (including radiative accelerations). We compared these models to others computed with the Montreal/Montpellier code and CESTAM.

The models computed with the three codes compared in this study (MESA, Montreal/Montpellier, and CESTAM) show a very high level of agreement in terms of abundance profiles and surface abundances (Sect. 4.2). However, due to differences in the opacity data used by the Montreal/Montpellier code, its abundance profiles for ²⁸Si and ³²S show slight differences from those of MESA and CESTAM.

We showed that some options dedicated to speed up atomic diffusion computation in MESA need to be set carefully. The most important are options 1 and 2 of Table A.1. These two options act on the treatment of atomic diffusion in the surface layers, making it possible to treat a certain region at the surface as a single cell for atomic diffusion calculations. This is very useful to speed up the computation, but leads to wrong abundance profiles when this region is larger than the natural surface convective zone. It mimics the addition of very efficient mixing at the bottom of the SCZ, modifying the predicted surface abundances. This also causes structural modifications, as the new chemical mixture changes the opacity in this zone, altering its size and the effective temperature.

The code comparison carried out in this paper showed that the MESA stellar evolution code compares well with the Montreal/Montpellier code and CESTAM. This kind of comparison and benchmarking of the codes is fundamental to assess the robustness of stellar evolution computation using different numerical methods, and atomic data for opacities and atomic diffusion calculation. It is also important in order to analyse present and future observations that missions such as Kepler, TESS, or PLATO provide or will provide.

Acknowledgments

This work was supported by FCT/MCTES through the research grants UIDB/04434/2020, UIDP/04434/2020 and PTDC/FIS-AST/30389/2017, and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalização (grant: POCI-01-0145-FEDER-030389). M.D. and D.B. are supported by national funds through FCT in the form of a work contract. We thank an anonymous referee for valuable comments which helped to improve the paper.

References

Appendix A: Analysis of the impact of MESA atomic diffusion options on surface abundance profiles

In this appendix we address the impact of the options mentioned in Sect. 5 in more detail. We primarily focus on the physical significance of the code options, before presenting quantitative results applied to the models. The list of these options, as well as their their default values, is presented in Table A.1. The models are defined in Table A.2, as well as a summary of the test results. Since some of these options have high computational costs, it may not be always feasible to choose values for these options that provide the maximum accuracy. Therefore, this appendix can be used as an indication, that should be adapted for each case, in order to reach the desired balance between the computational cost and the accuracy.

A.1. Surface mass and temperature cell control parameters

This appendix is dedicated to options 1 and 2 of Table A.1, introduced in Sect. 5. In the limiting case where the option 1 is equal to the fraction of mass corresponding to the surface convective zone (model A11 in Fig.3), the whole zone will be treated as a single cell for the diffusion computation. Since convection completely homogenizes the chemical composition, the treatment of atomic diffusion is then done correctly. However, if this value is larger than the fraction of mass of the SCZ, errors may occur. In this case, the mixed region is artificially forced to be deeper, introducing an artificial and efficient full mixing in regions where it should not be present (see models A12 and A13 in Fig.3) and incorrectly uniformizing the composition of the star.

This option can be particularly critical in more massive stars, where the SCZ are very thin. This is why the default value is suitable as far as the size of the surface convective zone does not become smaller than this value. Analogously, one must not set a value for option 2 that is higher than the temperature at the bottom of the surface convective zone (see Fig. A.2).

Indeed, the abundance differences caused by this artificial full mixing are already significant for the A12 and A22 models (Figs. A.1 and A.2), where the values of the options 1 and 2 are twice the mass of the SCZ and the temperature at its bottom, respectively.

Fig. A.1. Relative difference between the abundance profiles resulting from option 1 (diffusion_min_dq_at_surface) tests at 420 Myr for 1.4M⊙ models and the A0 model. The option 1 values for each model are given in Table A.2. The vertical dotted line represents the bottom of the surface convection zone of the A0 model.

Fig. A.2. Relative difference between the abundance profiles resulting from option 2 (diffusion_min_T_at_surface) tests at 420 Myr for 1.4M⊙ models and the A0 model. The option 2 values for each model are given in Table A.2. The vertical dotted line represents the bottom of the surface convection zone of the A0 model.

The average surface abundance differences in terms of mass fraction between these models and A0 are about 6% for A12, and 16% for A22 (Table A.2). For ¹²C, ¹⁴N and ⁴⁰Ca, these differences are almost doubled. This translates to a surface abundance difference up to 0.06 dex for the A12 model, 0.14 dex for the A22 model and 0.16 dex for models A13 and A23 (see Fig. D.1).

It should also be noticed that setting options 1 or 2 to M_SCZ or T_BSCZ, respectively, allows to save up to 15% of the computation time, compared to the default settings, without affecting significantly the abundance profiles. Moreover, the deeper the SCZ is, the larger the value of options 1 and 2 can be, hence a faster computation. This is particularly interesting for the less massive models.

A.2. Opacity upper temperature limit

The computation of the Rosseland mean opacity using OPCD3 monochromatic tables, in the A0 model, was restricted to the temperature interval [10⁴, 10^6.3] K. Outside this interval, opacity tables with a fixed metal mixture are used. This range is defined by options 3 to 6 in the following manner: if log(T) is in the interval [O₆, O₄], where O_i represents the value of the ith option, "OP_mono" opacities will be used. Standard opacity tables are used if log(T) > O₃ or log(T) < O₅, and in the remaining intervals both are partially used. In the models of this paper, given that O₃ = O₄ and O₅ = O₆, the "OP_mono" opacities are either fully used or not at all. The tests carried out in this appendix involved changing the upper limit (options 3 and 4, simultaneously).

This possibility of adapting the temperature range where monochromatic opacities are used can be a way to save computational time. In the regions of the model where the mixture of metal is not strongly affected by atomic diffusion, it is a reasonable approximation (see Hui-Bon-Hoa 2021 for more massive stars). In the following, we quantify its effect on the surface abundances.

Figure A.3 shows the relative differences between the abundance profiles of models including larger temperature intervals for the use of monochromatic opacities and the A0 model. The main impact of altering the temperature range is to slightly modify the stellar structure in the regions where the opacity tables are different (T > 10^6.3 K). This has a limited effect on the surface abundances, which indicates that for a reasonable accuracy/computation time compromise, an upper temperature of T = 10^6.3 K is adequate. It is important to note that changing this upper temperature value has an impact on the mixing length parameter obtained with a solar calibration. Hence, these differences in surface abundance may be even smaller.

Fig. A.3. Relative difference between the abundance profiles resulting from varying option 3 (high_logT_op_mono_full_off) at 410 Myr for 1.4M⊙ models and the A0 model. The models are defined in Table A.2. The vertical line represents the bottom of the SCZ of the A0 model.

Indeed, the average surface abundance difference between the A0 and A3 models is only 0.9%, while it is 4.1% between the A0 and A32. This translates into 0.008 and 0.04 dex, respectively (See Fig. D.1). Extending the upper limit to 10^7.3 K doubles the A0 model computation time. Furthermore, the A3 model, which only uses standard opacity tables, is about four times faster than the A0 one.

Lastly, while we have discussed the effects of altering the upper limit for opacity calculations using monochromatic opacities (options 3 and 4), there is also a lower temperature limit (10⁴K, options 5 and 6) implemented in the A0 model. We decided not to change it since, at lower temperatures, molecules would start to have a significant impact on the opacity. Given that their contribution is not taken into account in the current version of OP monochromatic opacities, it is more reasonable to use standard opacity tables in this temperature range.

A.3. Extra options

In this appendix we discuss the effects of two additional options that, while not directly related to atomic diffusion, are nevertheless worth mentioning, given that they can significantly alter the surface abundances of the models.

A.3.1. Radiation turbulence coefficient

Turbulent mixing counteracts the depletion of helium and metals in the outer envelope of stars. Radiative viscosity was thought to induce turbulent mixing of chemical elements (Morel & Thévenin 2002). However, the description of the effect of radiative viscosity on chemical element transport, as it was described in that paper, has been shown to be erroneous (Alecian & Michaud 2005).

It is possible to include this effect with an option available in MESA, i.e. option 7 (radiative_turbulence_coeff) of Table A.1. In order to test its effects we adopted two values: 0 (the default one) and 1. The average impact on the surface abundances is 4.2% (up to 9% for ¹²C and ⁴⁰Ca). Nevertheless, since this process is no longer physically justified (Alecian & Michaud 2005), its usage is not recommended.

A.3.2. Diffusion classes

In MESA, the treatment atomic diffusion is done with element classes (options 8, 9 and 10 of Table A.1), which give the possibility to group elements and follow their evolution as a single one. They define the number of representative classes of species for diffusion calculations, as well as their atomic mass. In the A0 models, all considered elements were associated with their own class. This is the only way to have atomic diffusion treated properly for these elements. These options are convenient when some elements are not needed for specific analysis, because they save computational time. In order to avoid any inconsistency in the opacity profiles of the models, the usage of one class per element is recommended.

Appendix B: Comparison between MESA, the Montreal/Montpellier code and CESTAM modes at 1.1 M_⊙

This appendix is dedicated to the comparison between the three codes for models A01, B01 and C01. The comparison of the abundance profiles is presented in Fig. B.2. The relative difference between the Montreal/Montpellier model and the two other codes is presented in Fig. B.1.

Fig. B.1. Same legend as Fig. 2 for models A01 and C01.

Fig. B.2. Same legend as Fig. 1 for models A01, B01 and C01.

Appendix C: Radiative acceleration profiles for A0, B0 and C0 models

In this appendix we present a comparison of the radiative acceleration profiles for models A0, B0 and C0 (Fig. C.1).

Fig. C.1. Radiative acceleration profiles according to log(ΔM/M*) (with ΔM being the mass above the radius r) obtained with MESA (A0 model, red solid lines), the Montreal/Montpellier code (B0 model, blue dashed lines) and CESTAM (C0 model, green dotted lines). The black solid line represent the gravity. The vertical grey dotted line represents the average position of the surface convective zone of the three models. The surface is towards the left part of the plots.

Appendix D: Surface abundance differences between A0 and, A1X, A2X and A3X models

In this appendix we present the surface abundances differences (in dex) between A0 and all models in Table A.1. In Figure D.1, [A/H]${}_{\mathrm{A}0}={log}_{10}{\left(\frac{X(A)}{X(H)}\right)}_{\mathrm{AX}}-{log}_{10}{\left(\frac{X(A)}{X(H)}\right)}_{\mathrm{A}0}$, where X(A) and X(H) are the mass fractions of the element A and hydrogen, respectively.

Fig. D.1. Difference of surface abundances between the model A0 and all the models presented in Table A.1. The grey dotted lines show the ±0.01 dex interval.

All Tables

All Figures

Fig. 1. Abundance profiles according to log(ΔM/M*) (with ΔM being the mass above the radius r) obtained with MESA (Ad model, dot-dashed orange lines; A0 model, solid red lines), the Montreal/Montpellier code (B0 model, dashed blue lines) and CESTAM (C0 model, dotted green lines) at 420 Myr. The dot-dashed grey lines represent the abundance profiles for the AD model without radiative accelerations. The vertical dotted grey line represents the average position of the surface convective zone of the three models. The surface is towards the left part of the plots.

In the text

	Fig. 2. Difference of surface abundances2 between the model B0 and the 2 other models (A0: red symbols, and C0: green symbols) shown in Fig. 1. The error bars of the C0 model are estimated using the 30% uncertainties on the radiative accelerations values obtained with the SVP method. The dotted grey lines show the ±0.01 dex interval.
In the text

	Fig. 3. Same legend as Fig. 1 for nickel.
In the text

	Fig. A.1. Relative difference between the abundance profiles resulting from option 1 (diffusion_min_dq_at_surface) tests at 420 Myr for 1.4M⊙ models and the A0 model. The option 1 values for each model are given in Table A.2. The vertical dotted line represents the bottom of the surface convection zone of the A0 model.
In the text

	Fig. A.2. Relative difference between the abundance profiles resulting from option 2 (diffusion_min_T_at_surface) tests at 420 Myr for 1.4M⊙ models and the A0 model. The option 2 values for each model are given in Table A.2. The vertical dotted line represents the bottom of the surface convection zone of the A0 model.
In the text

	Fig. A.3. Relative difference between the abundance profiles resulting from varying option 3 (high_logT_op_mono_full_off) at 410 Myr for 1.4M⊙ models and the A0 model. The models are defined in Table A.2. The vertical line represents the bottom of the SCZ of the A0 model.
In the text

	Fig. B.1. Same legend as Fig. 2 for models A01 and C01.
In the text

	Fig. B.2. Same legend as Fig. 1 for models A01, B01 and C01.
In the text

Fig. C.1. Radiative acceleration profiles according to log(ΔM/M*) (with ΔM being the mass above the radius r) obtained with MESA (A0 model, red solid lines), the Montreal/Montpellier code (B0 model, blue dashed lines) and CESTAM (C0 model, green dotted lines). The black solid line represent the gravity. The vertical grey dotted line represents the average position of the surface convective zone of the three models. The surface is towards the left part of the plots.

In the text

	Fig. D.1. Difference of surface abundances between the model A0 and all the models presented in Table A.1. The grey dotted lines show the ±0.01 dex interval.
In the text