© N. M. Kostogryz et al. 2022

1 Introduction

The discovery of the first exoplanet orbiting a solar-type star (Mayor & Queloz 1995) has led to the rapid development of various new methods for planet detection and characterization. By now more than 3000 exoplanets have been discovered and confirmed (see the Exoplanet Orbit Database¹; EOD). One of the most efficient ways of detecting exoplanets is offered by the transit method (Charbonneau et al. 2000), which became particularly useful together with the high-precision photometry made possible by space telescopes such as CoRoT (Baglin & Fridlund 2006), Kepler (Borucki et al. 2010), TESS (Ricker et al. 2015), and CHEOPS (Cessa et al. 2017). The transit photometry method aims to detect periodic drops in stellar flux due to the partial occultation of the stellar disk by the planet. The amplitude of the flux drop can be used to determine the planet-to-star radius ratio, which is one of the most important parameters for characterizing an exoplanet. The transit method is also used to study atmospheres of exoplanets, using transit spectroscopy (i.e., by simultaneously observing transits at several wavelengths). If the exoplanet has an atmosphere, the apparent size of the planet appears larger at wavelengths where stellar radiation is effectively absorbed or scattered by the planetary atmosphere (see, e.g., Seager & Sasselov 2000; Brown 2001). Consequently, studying the dependence of the apparent planetary radius on wavelength allows the chemical composition of the planet’s atmosphere to be constrained. The determination of the planetary radius from transit measurements is one of the most important tasks in transit photometry and spectroscopy. Such a determination requires accurate knowledge of the center-to-limb variation of stellar intensity (hereafter, limb darkening) since it defines the shape of the transit light curve and affects its depth.

High-precision photometric measurements have shown that the present treatment of limb darkening is not sufficiently accurate and leads to systematic errors in the derived parameters of the exoplanets (e.g., Espinoza & Jordán 2016; Morello et al. 2017a; Maxted 2018). Typically, limb darkening is represented by rather simple laws, such as a linear law (Schwarzschild 1906), a quadratic law (Kopal 1950), a square-root law (Diaz-Cordoves & Gimenez 1992), a power-2 law (Hestroffer 1997), or a four-coefficients law (Claret 2000), so that when modeling the transit light curves, limb darkening is parameterized by some set of coefficients. Ideally, these coefficients should be constrained by the stellar modeling. Consequently, many libraries of limb-darkening coefficients covering a wide range of effective temperatures (T_eff), surface gravity (log g), and metal-licities (M/H) have been produced (e.g., Claret 2000; Sing 2010; Claret & Bloemen 2011; Magic et al. 2015) using various radiative transfer codes, such as ATLAS (Kurucz 1993), NextGen (Hauschildt et al. 1999), PHOENIX (Husser et al. 2013), MARCS (Gustafsson et al. 2008), and STAGGER (Magic et al. 2013). However, the limb-darkening parameters diverge between different libraries and often lead to an inadequate quality of fits to the observed transit profiles (Csizmadia et al. 2013). First, this can be due to errors in limb-darkening coefficients introduced by interpolation from the grid of stellar parameters used in these libraries to the actual stellar fundamental parameters. Second, available stellar calculations may not treat mechanisms that affect limb darkening with sufficient accuracy, for instance, convection (Pereira et al. 2013; Chiavassa et al. 2017) or magnetic activity (Csizmadia et al. 2013). Therefore, limb-darkening coefficients are often left as free parameters in a least-squares fit to observed light curves (e.g., Southworth 2008; Claret 2009; Cabrera et al. 2010; Gillon et al. 2010; Csizmadia et al. 2013; Maxted 2018). While this method usually leads to a good quality fit to observed transit profiles, it introduces additional free parameters, resulting in possible biases and degeneracies in the returned planetary radii (Espinoza & Jordán 2015; Morello et al. 2017b). The way to reduce these biases and reliably determine planetary radii is to improve theoretical computations of stellar limb darkening.

With this work, we begin our effort to produce a new library of stellar limb darkening. We employ the recently developed 1D Merged Parallelized Simplified ATLAS (MPS-ATLAS, Witzke et al. 2021) code to create a grid of model atmospheres and the corresponding limb darkening for a wide range of stellar fundamental parameters, utilizing the most accurate information on stellar abundances and mixing lengths. The code has been intensively tested against solar and stellar observations, and incorporates improved molecular and atomic data. It is substantially faster in comparison to other available codes, allowing us to calculate limb darkening on a very fine grid of stellar parameters and reduce potential errors in interpolation.

This work is organized as follows. In Sect. 2 we describe the computations of stellar model atmospheres on a grid of fundamental stellar parameters. We note that while the main goal of this paper is to achieve accurate calculations of the limb darkening, we also provide the grids of stellar atmospheric models, which can be used for a large range of applications. The synthesis of stellar spectra and their convolution with different passbands is presented in Sect. 3. Section 4 describes the derivation of coefficients for several limb-darkening laws. We compare our calculations to solar (Neckel & Labs 1994) and stellar (Maxted 2018) observations, as well as to various available libraries of limb darkening in Sect. 5. Finally, in Sect. 6 we summarize our results and briefly describe the outlook.

2 Model Atmospheres

To compute the limb darkening in different broadband passbands, we first computed 1D plane-parallel model atmospheres on a grid of fundamental stellar parameters, namely effective temperature (T_eff), metallicity (M/H), and surface gravity (log g), that was a lot finer than previous model atmosphere grids. In Table 1 we summarize the properties of our grid of stellar models. It covers the stellar parameters of most of the known exoplanet host stars (see EOD), which are mainly dwarfs and subgiants of spectral types M–F. In Fig. 1 we compare our grid of stellar models on metallicity and effective temperature to other two grids that are widely used by the community: PHOENIX (Husser et al. 2013) and MARCS (Gustafsson et al. 2008). Clearly, our grid covers the largest range of metallici-ties with a finer sampling than the other grids. In particular, for the range of metallicities of all known stars with exoplanets (M/H = [−1.0 : 0.5], see EOD) it has a particularly fine sampling of 0.05. Foreffective temperatures ourgrid covers a similar range and has the same temperature step as the PHOENIX models, but with 100 K it has only half the temperature step of the MARCS models, which have a larger stepping (200 K) for models with T_eff > 4000 K. For surface gravity, we also have a better sampling than the other grids, especially in the range of log g of most known stars (log g = [4.0–5.0]). Although the surface gravity has a small effect on stellar limb darkening, such a fine grid of model atmospheres might be useful for other purposes, such as for calculations of center-to-limb variations of polarization (see, e.g., Kostogryz & Berdyugina 2015).

For the grid of stellar parameters, we computed stellar models using the MPS-ATLAS code. This code is based on the ATLAS9 (Kurucz 1993) and DFSYNTHE (Castelli 2005; Kurucz 2005) codes, which were merged into MPS-ATLAS. MPS-ATLAS is also an extension of these codes, for instance, in its ability to treat optimum opacity distribution functions (ODFs, Cernetic et al. 2019). The parallelization of the code allows computations to be performed much faster. Moreover, it was made more flexible and user-friendly. MPS-ATLAS includes treatment of the convection described by the mixing-length theory of turbulent transport (Böhm-Vitense 1958), and parameterized by the mixing-length parameter α. This parameter characterizes the efficiency of energy transport in the atmosphere. It also allows the convective material to overshoot when computing the model atmosphere.

The line blanketing, which is very important in the construction of model atmospheres, was taken into account by using the ODF approach. We computed ODFs using more than 100 million of atomic and molecular transitions from Kurucz (2005). The line shape was approximated by the Voigt profile using a micro-turbulent velocity of 2 km s⁻¹ (see Kurucz 2005; Witzke et al. 2021). The continuum opacity was calculated taking into account contributions from both absorption and scattering processes. For the continuum absorption, we considered free-free and bound-free transitions in H⁻, H , He, He⁻, C, N, O, Ne, Mg, A1, Si, Ca, Fe, the molecules CH, OH, and NH, and their ions. For the scattering contribution, we included electron scattering and Rayleigh scattering on HI, HeI, and H₂. A more detailed description of all the physical processes and numerical schemes applied in MPS-ATLAS is given in Witzke et al. (2021).

In addition to the grid of fundamental parameters, there were two further important inputs to our simulations, namely chemical composition and mixing-length parameters. The previous grid of ATLAS9 model atmospheres by Castelli & Kurucz (2003) was computed utilizing Grevesse & Sauval (1998) abundances and prescribing a mixing-length parameter of α = 1.25 for all computed atmospheres. As MPS-ATLAS is based on the ATLAS9 code, we first presented the set of models with the same abundances and mixing-length parameter. Hereafter, we call this set of models as “set 1”. Our second set of models (“set 2”) was based on more recent chemical abundances (Asplund et al. 2009). For set 2 we also accounted for the dependence of the mixing-length parameter on stellar fundamental parameters (T_eff, M/H, and log g) using an approximation on α/α_⊙ derived by Viani et al. (2018) from asteroseismic measurements of about 450 stars in the range of stellar parameters of 4500 K ≤ T_eff ≤ 7000 K, −0.7 ≤ M/H ≤ 0.5, and 3.3 ≤ log g ≤ 4.5. We scaled these dependencies to α_⊙ = 1.25 (as proposed by Kurucz 1993). In Fig. 2 we present an example of α dependence on the stellar effective temperature for three different values of metallicity and the same surface gravity (log g = 4.5). We did not extrapolate α outside the ranges of stellar parameters for which the relationship was derived, but took the boundary value as a constant α for all models outside the range (the dashed lines in Fig. 2). In Table 2 we summarize the description of different sets of stellar model grid computations.

To compute a model atmosphere, the MPS-ATLAS code takes an initial starting model and iterates it until radiative equilibrium (RE) is satisfied. We note that in the upper convection region, MPS-ATLAS allows the model to deviate from RE due to convection and overshooting (see details in Witzke et al. 2021). The model is converged if the maximum relative temperature deviation at each point in the atmosphere is less than 10⁻⁵ or until it reaches the maximum number of iterations, which is set to 1000. The latter case happens rarely and mostly for hot metal-rich stars. For comparison, we note that ATLAS9 models by Kurucz (1993) and Castelli (2005) are computed considering only 15 iterations. For set 1, we took starting models with the closest fundamental parameters from the ATLAS9 grid. For set 2, the starting model was taken from the set 1 model with the corresponding fundamental parameters. We note, however, that the change of the starting model does not affect the final model atmosphere structure (but might strongly affect the number of iterations needed to reach the convergence).

In Fig. 3 we present an example of the temperature structures computed for the model atmospheres of set 1. To illustrate how the temperature stratification depends on different stellar parameters, namely the effective temperature, surface gravity, and metallicity, we fixed two of them and varied the third. Figure 3 shows that increasing the T_eff shifts the entire structure to higher temperatures, but also produces some changes in the details of the stratification. The effect of metallicity on temperature structures is different. For a fixed effective temperature, metal-poor stars have smaller opacities, and therefore require a smaller temperature gradient in order to transport the same energy through the atmosphere (middle panel in Fig. 3). Vice versa, metal-rich stars have larger opacities and, thus, steeper temperature gradients. The surface gravity does not affect the model structure computation, except for small effects in the region where convection and overshooting play a role. We note that we allow the convective material to overshoot up to one pressure scale height for almost all models except for some models with low effective temperatures and low metallicities, for which we suppress overshoot to avoid divergence in iterations of model atmosphere structures. The energy influx brought by the overshoot causes a decrease in the temperature gradient. It is visible as a bump at τ_Ross ≈ 1 (where τ_Ross is the Rosseland optical depth) in the atmospheric temperature structures of stars with T_eff ≳ 6000 K. For the same stars at τ_Ross ≈ 10 another decrease in the temperature gradient caused by hydrogen ionization happens and becomes prominent as a second bump in the temperature structure.

The two sets of stellar model atmospheres have been uploaded as online material of this paper at the Centre de donnees astronomiques de Strasbourg (CDS). In addition to model structures, the uploaded files also contain information about chemical composition, convection, and overshoot parameters, as well as maximum temperature alteration at the last iteration (which is 10⁻⁵ in most cases but can be larger for a couple of models if 1000 iterations were not enough for convergence). In Appendix D, we describe an example of a model structure with all additional parameters we used.

Fig. 1 Representation of metallicity and effective temperature grids. Black dots show the grid from this paper. Blue dots depict the MARCS grid from Gustafsson et al. (2008) and red dots represent the PHOENIX grid from Husser et al. (2013).

Fig. 2 Mixing-length parameter dependence on the stellar effective temperature for three values of metallicities (M/H = [−0.5, 0.0, 0.5]) and a surface gravity of log g = 4.5. The solar mixing-length parameter (α⊙ = 1.25) is taken from Kurucz (1993). Different colors correspond to different metallicities, as described in the legend. Horizontal dashed lines show the α values that are used in computations outside the region of effective temperatures for which the mixing-length behavior was derived by Viani et al. (2018).

Fig. 3 Temperature structure of several models from set 1 plotted as a function of Rosseland optical depth. The three panels show models when one of the three stellar parameters is varied and the other two are fixed. Fixed parameters are presented in the title to each panel and the varying parameter in the legend. We plot every fifth model from the effective temperature and metallicity grids.

3 Spectral Synthesis

The spectral synthesis computations for the two sets of model atmospheres were performed using the same radiative transfer solver, turbulent velocity, and continuum opacities used for computing atmospheric models in Sect. 2. We computed the limb-darkening profiles for the passbands of instruments that are widely used to detect and characterize transiting exoplanets, namely, Kepler (Borucki et al. 2010), CHEOPS (Cessa et al. 2017), and TESS (Ricker et al. 2015). We also provide preliminary (not official) limb darkening profiles for the upcoming PLATO mission (Rauer et al. 2014). The detectors of the telescopes are charged-coupled devices (CCDs), which count photons and do not measure energy directly. Thus, to compute the specific intensity, I(λ,µ), emitted from the stellar atmosphere at different wavelengths (λ) and view angles (µ = cos(θ), where θ is the angle between the normal to the stellar surface and the direction to the observer), we need to convert units from the energy flux to photon number flux:

$$I_{\text{pb}}\mu ={\displaystyle {\int }_0^{\infty }{R}_{pb}\lambda I\lambda ,\mu /hc/\lambda \text{d}\lambda },$$(1)

where R_pb is the response function of the instrument with a CCD detector, “pb” denotes the passband of different space missions, h is the Planck constant, and c is the speed of light. This transformation strongly affects the spectral profile of emergent intensity in broadband passbands.

In Fig. 4 we present the transmission curves of the (broad) passbands considered in this paper, together with the example of limb darkening curves in these filters for a solar analog with M/H = 0.0, T_eff = 5800 K, and log g = 4.5. We note that for the PLATO passband, we used the preliminary transmission curve published by Marchiori et al. (2019), which may change in the future. Figure 4 shows that the limb darkening in Kepler, CHEOPS, and PLATO are very close to each other, with some small deviation at the limb in the PLATO passband, while limb darkening in the TESS passband deviates substantially for the same target. The deviations in the TESS passband are because it is shifted to longer wavelengths and transmits stellar radiation in the near-infrared region, cutting all the radiation below 600 nm (see left panel in Fig. 4).

We utilized the ODFs to calculate the specific intensity, Ι(λ, µ), in Eq. (1). While the ODF approach has been routinely used in calculating grids of limb darkening (e.g., Claret & Bloemen 2011), there have been concerns that its performance might deteriorate for cool metal-rich stars (Ekberg et al. 1986). In Appendix A we introduce a modification of the ODF setup for cool metal-rich stars, which leads to a more accurate spectral synthesis. We tested limb darkening profiles calculated with the ODF approach against those calculated utilizing spectral calculations with high-spectral resolution (R = 500 000) and show that our setup allows very accurate calculations of limb darkening.

In the online tables (CDS), we provide the absolute value of intensity at the disk center and the center-to-limb intensity normalized to the disk-center intensity for 24 disk positions in different passbands for two sets of our grid. We present an example of the table in Appendix E.

Fig. 4 Normalized broadband transmission curves of the exoplanetary space missions, Kepler, TESS, CHEOPS, and PLATO, considered in this study (left panel) and the corresponding limb darkening for a star with M/H = 0.0, Teff = 5800 K, and log g = 4.5 (right panel). Different colors indicate different filters described in the legend.

4 Parameterization of Stellar Limb Darkening

Intensities, normalized to the disk center, I_pb(µ)/I_pb(1.0) are often approximated by different limb-darkening laws that are linear or nonlinear polynomial expansions in µ. Different authors have presented a variety of limb darkening-laws, for example, linear (Milne 1921), quadratic (Kopal 1950), square root (Diaz-Cordoves & Gimenez 1992), logarithmic (Klinglesmith & Sobieski 1970), power-2 (Hestroffer 1997; Maxted 2018), and four-parameter nonlinear (Claret 2000). The four-parameters law describes the detailed shape of the limb-darkening profile for a given filter in the most accurate way (e.g., Morello et al. 2020). It is written as follows:

$$\frac{I_{\text{pb}}\mu }{I_{\text{pb}}1.0}=1-{\displaystyle \sum _{j=1}^4{a}_{j}1-{\mu }^{j/2},}$$(2)

where a_j are four coefficients of the limb-darkening law. Among limb-darkening laws with two coefficients, Morello et al. (2017a) found that a power-2 law outperforms other laws, in particular, for cool stars. This law is presented as

$$\frac{I_{\text{pb}}\mu }{I_{\text{pb}}1.0}=1-c1-{\mu }^{\alpha },$$(3)

where c and α are the two coefficients that are a function of stellar parameters. This law accurately matches the limb-darkening profile due to the exponent of µ.

In this paper, the center-to-limb variations of intensity convolved with the transmission profiles of different space mission instruments were fitted by both limb-darkening laws, given by Eq. (2) and by Eq. (3). We used a nonlinear least squares fitting procedure from the python scipy library² to fit the limb darkening curves of our grids by these two laws. The derived coefficients of the power-2 (Eq. (3)) and the four-coefficients (Eq. (2)) laws are provided in online tables for the two grid sets in all passbands. An example of the table is given in Appendix F. We note, however, that since a limb darkening curve cannot be exactly represented by parameterization, the fitting procedure might introduce biases (see Appendix C for an example). We provide a list of all online tables with titles in Table 3.

An important feature of our calculations is that they are available on a very fine grid of fundamental parameters (in particular in metallicity, see Fig. 1). This allows us to minimize the error associated with the interpolation of the limb darkening coefficients from the grid points to specific stellar parameters (see Sect. 1). As an example in Fig. 5, we show how the power-2 law coefficients, namely α and c from set 1 in the Kepler passband, vary with stellar metallicity and effective temperature. It is important to note that the variation of the coefficients on stellar surface gravity is tiny, so we do not present it in Fig. 5. One can see that the dependence of the coefficients on stellar parameters is not linear. This nonlinearity causes the errors when using linear interpolation, especially when the step between the grid points is not fine enough.

To give an example of possible values of the interpolation errors, we selected our limb-darkening coefficients from “set 1” on the metallicity grids of MARCS and PHOENIX models and linearly interpolated them to our metallicity grid (see top panels of Fig. 6). In the bottom panels of Fig. 6, we present the ratio of the limb darkening curves directly derived from our “set 1” coefficients to the limb darkening derived from the interpolated coefficients for metallicity in the range from −2.0 to 0.5 for the MARCS models and −2.0 to 1.0 for the PHOENIX models. Interestingly, Fig. 6 shows that for stars with solar effective temperature, the errors associated with the interpolation from MARCS and PHOENIX grids can reach 0.4 and 0.6%, respectively. The interpolation errors increase toward later spectral types and decrease toward earlier spectral types (compare Fig. B.1 to Figs. B.2 and B.3).

Fig. 5 Limb-darkening coefficients of the power-2 law on the stellar metallicity (top panel) and effective temperature (bottom panel) grids in the Kepler passband. The different colors depict the coefficients (α and c). Fixed parameters are labeled on each panel.

5 Validation of Limb Darkening Computations

In this section, we present the validation of our computed limb darkening for the Sun and other stars by comparing the computations with solar (Neckel & Labs 1994) and stellar (Maxted 2018) measurements, respectively. We also compared the computed limb-darkening curves with other computations (e.g., Sing 2010; Claret & Bloemen 2011; Claret 2017).

5.1 Solar Computation-to-measurements Comparison

Here we compare our calculations for a star with solar fundamental parameters to available solar measurements. Various authors (e.g., Pierce et al. 1977; Neckel & Labs 1994) measured the solar continuum spectrum in a broad wavelength range at different view angles. In this work, we make use of the observations by Neckel & Labs (1994) as they cover the entire wavelength range of the passbands considered in our study. The limb-darkening computations are done for the solar parameters of M/H = 0.0, T_eff = 5777 K, log g = 4.4377 and solar abundances from Asplund et al. (2009), meaning they correspond to set 2.

In the upper panel of Fig. 7, we show the solar spectrum computed with the MPS-ATLAS code, marking wavelengths at which the Neckel & Labs (1994) measurements were performed. In the lower panel of Fig. 7, we show the computed and measured values of intensity at several view angles, normalized to the disk-center intensity. In addition to our final product (i.e., computations with convection and overshoot), we also plot limb darkening computed assuming RE (i.e., no convection and no overshoot) and limb darkening computed taking convection into account but neglecting overshoot. One can see that MPS-ATLAS solar limb-darkening computations with overshoot accurately agree with the observations at all wavelengths and view angles. The convection without overshoot heats the atmosphere in very deep layers, which contribute only marginally to the emergent radiation. Shortward of ~400 nm, the continuum opacity increases due to multiple photoionization processes and line haze, so that the contribution from the layers affected by convection can be neglected. Longward of ~550 nm, the temperature sensitivity of the Planck function decreases so that the temperature change due to convection does not noticeably affect limb darkening. Therefore, we only see the deviations between pure RE and convection models, between ~400 and ~550 nm. In comparison to the convection, the overshoot affects the temperature structure of higher layers so that the agreement with measurements significantly deteriorates if overshoot is neglected. This result is in line with Pereira et al. (2013), who showed that limb darkening computed assuming RE does not agree well with observations.

5.2 Computation-to-computation Comparison for a Solar–analog Star

After testing our computations against the solar measurements, we performed a validation test for a solar-analog star (M/H = 0.0, T_eff = 5800 K, log g = 4.5) against other computations in the Kepler and TESS passbands. First, we took the coefficients for the four-parameters limb-darkening law (Eq. (2)) from the most widely used libraries in the Kepler (Claret & Bloemen 2011) and in TESS (Claret 2017) passbands for the closest stellar parameters to our solar analog, namely for M/H = 0.0, T_eff = 5750 K, and log g = 4.5. While Claret & Bloemen (2011) also give coefficients based on the PHOENIX library, we selected the coefficients computed based on the ATLAS9 library. Using these coefficients, we computed limb darkening curves and we refer to them as “CB11“ (available for Kepler passband) and “C17“ (available for TESS passband). Second, we took the limb darkening coefficients derived by fitting the same four-parameters law in each passband from our sets and computed limb darkening using these coefficients. We refer to these computations as “LD1” and “LD2” for the limb darkening computed from set 1 and set 2 libraries, respectively.

In Fig. 8 we compare LD1 and LD2 computed with and without overshoot to CB11 and C17 limb darkening in the Kepler and TESS passbands, respectively. The left panels of Fig. 8 show that in the Kepler passband, the CB11 limb darkening is closer to LD1 and LD2 computed with overshoot. The maximum difference between CB11 and our calculations with overshoot is slightly above, 1% close to the stellar limb. At the same time, our computation without overshoot deviates up to ~3% from CB11. The right panels of Fig. 8 show that in the TESS passband, the maximum difference between our computations and C17 is about 2%. Interestingly, the difference between LD1 and LD2, arising mainly due to the difference between Grevesse & Sauval (1998) and Asplund et al. (2009) abundances, is significantly more prominent for the models with overshoot.

The discrepancy of limb darkening in both passbands between different computations are quite significant and can lead to bias in exoplanet radius determination. In order to understand how accurate our limb darkening (LD1 and LD2) are, we compared them to the data that mimic solar observations in the Kepler passband. Since the Sun was never observed in the Kepler passband, we emulated the observations by applying the following method. We computed the difference between the solar observations (Neckel & Labs 1994) and our computations with overshoot (from Sect. 5.1) at each wavelength between 400 nm and 900 nm where observations were available and for each view angle. The difference monotonously changed with wavelength at each µ. This monotonicity prompted us to interpolate the difference to the MPS-ATLAS wavelength grid in the Kepler passband, obtaining corrections to our computations as a function of wavelength and view angle. After that, we applied these corrections to our calculations, convolved them with the Kepler transmission passband, and applied the fitting procedure for the four-parameters limb-darkening law. We refer to the derived limb darkening as “NL_Kepler”. In Fig. 9 we present the ratio of LD2 and CB11 to NL_Kepler. In addition, we present the limb darkening of a solar-analog star in the Kepler filter computed by Sing (2010) using ATLAS9 models. We show that the deviations between our limb darkening and NL_Kepler reach up to 0.5%, while for other previous calculations, larger deviations occur, especially closer to the limb.

Fig. 6 Comparison of the limb darkening at our, MARCS, and PHOENIX grids of stellar parameters in the Kepler passband. Top panels: limb darkening coefficients of the power-2 law as a function of stellar metallicity. The vertical dashed lines show the metallicity grid points of the MARCS and PHOENIX models (left and right panels, respectively). Bottom panels: ratios between the limb darkening curves produced from coefficients directly calculated on our fine metallicity grid and those produced from the same coefficients but selected on the MARCS and PHOENIX metallicity grid points, and linearly interpolated to our fine metallicity grid (left and right panels, respectively). Ratios corresponding to different metallicity values are indicated by different colors. Coefficients and the limb darkening ratios are presented for the fixed Teff = 5800 K and log g = 4.5.

5.3 Stellar Computation-to-measurements Comparison

Recently, Maxted (2018) presented limb darkening coefficients

$$h_1=\frac{I_{\text{pb}}0.5}{I_{\text{pb}}1.0}=1-c1-2^{-\alpha }$$(4)

and

$$h_2=\frac{I_{\text{pb}}0.5}{I_{\text{pb}}1.0}-\frac{I_{\text{pb}}0.0}{I_{\text{pb}}1.0}-c{2}^{-\alpha }$$(5)

of the power-2 law derived from light curves for transiting exo-planet systems observed with Kepler. He selected the Kepler targets for which the stellar parameters, such as M/H, T_eff, and log g were determined using high-precision spectroscopy. He compared the obtained coefficients to STAGGER calculations by Magic et al. (2015) interpolated to the fundamental parameters of the Kepler targets (see, Fig. 4 in Maxted 2018).

In order to compare our models to observations, we chose to calculate the h₁ and h₂ coefficients for the same set of stars as in Maxted (2018). Since 1D MPS-ATLAS computations are very fast, we performed model atmospheres and limb darkening computations directly for the considered Kepler targets, avoiding any interpolation. We directly inferred the h₁ and h₂ coefficients using Eqs. (4)–(5) from our calculations. In Fig. 10 we present the difference between observations and the two computations (STAGGER and MPS-ATLAS) for h₁ and h₂. For most cases, our computations are closer to observations. However, for several Kepler targets, the 3D STAGGER limb-darkening coefficients show better agreement with observations. We note that, in contrast to the MPS-ATLAS computations, the difference between STAGGER computations and observations also contains uncertainties brought about by the interpolation. Thus, our result does not imply that 1D MPS-ATLAS computations are more accurate than 3D computations, but rather indicates that the combined error of 3D computations and interpolation (which is unavoidable for grids of computationally expensive 3D models) might be even larger than the error of 1D computations.

Interestingly, Fig. 10 shows that there are systematic differences between the h₁ and h₂ parameters deduced from observations and from modeling. Namely, h₁ values from observations are systematically larger than those from modeling, while h₂ values from observations are systematically lower than those from modeling. It should be noted that neither the error bars of measurements (see Table A.1. in Maxted 2018) nor the computation uncertainties coming from nonideal stellar parameter determination (priv. comm. with H.-G. Ludwig) can explain the systematic bias of the coefficients. In a follow-up paper, we will show that deviations in h₁ can be explained by stellar surface magnetic fields (which are not accounted for in MPS-ATLAS or in STAGGER computations). The interpretation of the deviation in the h₂ values is not so clear because the value obtained by fitting the transit light curve of a hot Jupiter may be biased by the extraction procedure (see, Appendix C). This conclusion is based on using a power-2 transit model to fit simulated transit light curves based on a realistic solar limb-darkening profile. The bias in the h₁ values obtained from transit fitting is much lower (<0.001), provided that the impact parameter for the transit is b ⪅ 0.8. More work is needed to better understand the best choice of parameters for testing limb-darkening profiles from stellar models against observations of transit light curves for real stars.

Fig. 7 Solar spectrum (top panel) and a comparison of solar limb-darkening computations with solar measurements by Neckel & Labs (1994) (bottom panel). Stars in both panels depict the solar measurements, and lines represent computations. Measurements and computations of solar spectrum in the top panel are shown at the disk center. Spectra at selected view angles, namely the µ = 0.8, 0.6, 0.4, and 0.2, normalized to the solar spectrum at the disk center, are shown in the bottom panel. The different line types marked in the legend to the figure correspond to the computation with different approximations: RE, with convection (α = 1.25), and with convection and overshoots. Intensity is measured in [erg cm−3 s−1 sr−1].

6 Summary and Outlook

In this paper we have employed our recently developed MPS-ATLAS code to compute a new library of stellar model structures and to synthesize stellar limb darkening in different broad passbands (Kepler, TESS, CHEOPS, and PLATO) for an extensive and fine grid of stellar parameters. We fit the “power-2” and “4-coefficients” laws to the computed limb darkening and derived the limb-darkening coefficients. We created two grids of stellar atmospheres and limb darkening parameters with different solar abundances (Grevesse & Sauval 1998; Asplund et al. 2009) and different mixing-length parameter considerations (constant or varying for different stellar parameters). We have made the full set of computations (model atmosphere structure, limb darkening, and their coefficients) for our two sets available at the CDS.

Our computations successfully passed different validation tests, such as comparisons with solar measurements, with other computations of stellar limb darkening, and with stellar limb-darkening coefficients derived from Kepler observations. We found a very good agreement between our computed solar limb darkening and measurements at different wavelengths and viewing angles. Our computations agree with Kepler measurements at the same level as 3D STAGGER calculations, and for many of the stars even slightly better than STAGGER. However, the sets of coefficients derived from STAGGER and our computations show a systematic offset relative to the observational data. This offset can come from ignoring the effect of magnetic activity in limb darkening computations. In a forthcoming paper, we will quantify the effects of magnetic activity on stellar limb darkening using a 3D radiative-magnetohydrodynamic (MHD) simulation with the MURaM code (Vögler et al. 2005). Subsequently, we will extend the library presented here to the case of magnetized stars.

Fig. 8 Comparison of the MPS-ATLAS limb-darkening computations to that by CB11 (for the Kepler passband, left panels) and by C17 (for the TESS passband, right panels) for a star with solar fundamental parameters. Top panels: describe the limb darkening curves and the bottom panels present their ratios. Calculations with MPS-ATLAS (LD1 and LD2) are given, taking the overshoot into account (blue) as well as ignoring it (red). See text for the explanation of the abbreviations.

Fig. 9 Ratios between limb darkening from LD2, CB11, and Sing (2010) (Sing) to the NL_Kepler limb darkening. Different colors correspond to different ratios (see legends).

Fig. 10 Comparison of the limb darkening coefficients. The stars symbols show the difference between observations (Maxted 2018) and computations done with STAGGER (blue) and MPS-ATLAS (red) codes in Kepler passbands. The horizontal dashed line represents the zero level.

Acknowledgements

We thank the referee, Giuseppe Morello, for constructive and useful comments on the manuscript. We are also thankful to H.-G. Ludwig for our fruitful discussion and to Jesper Schou for useful comments on our manuscript. The preliminary PLATO passband transmission curve (Marchiori et al. 2019) was kindly made available to us by Martin Pertenais. N.M.K. and L.G. acknowledge support from the German Aerospace Center (DLR FKZ 50OP1902). V.W. and A.I.S. were funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant no. 715947). L.G. acknowledges funding from ERC Synergy Grant WHOLE SUN 810218 and from NYUAD Institute Grant G1502.

Appendix A High-resolution Calculations Versus the ODF Approach

In this section we test the ODF approach, used to produce limb-darkening dependencies presented in our study, against direct (i.e., without utilizing ODF) calculations on a fine spectral grid (hereafter, high-resolution calculations).

For the high-resolution calculations, we computed I(λ, µ) in Eq. 1 with a spectral resolution of R=500000. For the ODF calculations presented in this study, we used the standard “little” wavelength binning grid (Castelli 2005). In each bin, we first sorted the opacities in ascending order, and then we split the bin into further sub-bins. In each of the sub-bins, we then calculated the mean value of the opacity and, instead of solving the radiative transfer equation on the high-resolution spectral grid, we solved it for each of the sub-bins. Summing up the weighted values of intensity from each of the sub-bins within the bin (with weights given by the ratios of sub-bin and bin sizes), we then obtained approximate values of the intensity in the bins (see the detailed description of the ODF procedure and its performance in Cernetic et al. 2019). For most of the calculations (see below for the explanation of the exceptions), we utilized the ODF setup proposed by Castelli & Kurucz (2003) splitting each bin into 12 sub-bins (with a 1/10 bin width for the first nine sub-bins with the lowest opacity values, and then with 1/20, 1/30, and 1/60 bin widths for the remaining three sub-bins). Such a nonuniform splitting of bins into sub-bins allowed us to accurately account for the cores of strong spectral lines (see Cernetic et al. 2019).

Fig. A.1 Ratios of intensity values in the Kepler passband computed using ODF to that using high-resolution opacity (resolving power R=500000). Calculations were performed for several pairs of metallicities and effective temperatures, but the same value of the surface gravity (log g = 4.5).

Figure A.1 shows that the ODF approach using the setup by Castelli & Kurucz (2003) gives very accurate values of intensity (with errors below 0.2%) except for cold and/or metal-rich stars. The deterioration of the ODF performance for cold stars is not surprising. The error of the ODF approach has two main components. The first component is due to the replacement of the intensity averaged over the sub-bin by intensity calculated using opacity averaged over the sub-bin. The second component is brought about by the failure of the ODF approximation. According to this approximation, the sorting of opacities is the same over the entire atmospheric region, contributing to the emerging intensity (see detailed discussion in Ekberg et al. 1986). The decrease in the sub-bin widths leads to the decrease in the first component of the error (since smaller sub-bins correspond to averaging of the closer values of opacity) but leads to the increase in the second component (since the allocation of frequencies into sub-bins becomes more susceptible to the failure of the ODF approximation). For the cold and metal-rich stars, the opacity is brought about by the mixture of lines from different species. These species (in particular, molecules which are the most important opacity source in the atmospheres of cold stars) have different temperature sensitivities and, thus, form in different parts of the stellar atmosphere. This leads to a significant change in the opacity profile within the atmosphere, violating the ODF approximation. Consequently, the second component of the ODF error dominates over the first one. Cernetic et al. (2019) showed that the best strategy in such a case is to split the bin into sub-bins of equal width. Here we followed up on this result and, in addition to the sub-binning proposed by Castelli & Kurucz (2003), also utilized an alternative sub-binning, splitting the bins into four sub-bins of equal width (hereafter, optimized sub-binning).

Figure A.2 shows a comparison of the performance of the Castelli & Kurucz (2003) and the optimized sub-binning for stars with T_eff = 3500 K. One can see that the optimized sub-binning performs better for the disk-center calculations (top panels of Fig. A.2) but the Castelli & Kurucz (2003) sub-binning performs better for the near-limb calculations (bottom panels of Fig. A.2). This is because photons emitted close to the limb come from a relatively narrow region of stellar atmospheres (in other words, the contribution function for the µ = 0.2 intensity is narrower than that for the µ = 1 intensity). As a result, the ODF approximation is more precise for the near-limb calculations and smaller sub-bin widths lead to better results by reducing the first component of the ODF error.

Figure A.3 shows how the RMS error of ODF calculations (averaged over the 400–550 nm spectral domain) depends on the view angles. In line with the discussion above, while the error of the ODF calculations with the optimized sub-binning only barely depends on the disk position, the error of the calculations with the Castelli & Kurucz (2003) sub-binning substantially decreases from the disk center to the limb.

All in all, we opted to replace the calculations based on the Castelli & Kurucz (2003) sub-binning with calculations based on the optimal sub-binning for µ ≥ 0.5 for the following values of M/H and T_eff (independently of log g values):

$$\begin{array}{c}0.3\le \text{M/H}\le 0.45:{\text{T}}_{\text{eff}}\le 4200\text{\hspace{0.17em}K};\\ 0.5\le \text{M/H}\le 0.9:{\text{T}}_{\text{eff}}\le 4500\text{\hspace{0.17em}K};\\ 1.0\le \text{M/H}\le 1.5:{\text{T}}_{\text{eff}}\le 5000\text{\hspace{0.17em}K}.\end{array}$$(A.1)

Fig. A.2 Ratio of emergent intensities at µ = 1.0 (top panels), µ = 0.5 (middle panels), and µ = 0.2 (bottom panels), calculated using the ODF approach to intensities calculated on a fine spectral grid and then averaged over bins used in the ODF calculations. The ratios are shown for the models with M/H = 0.5 (left column) and M/H = 1.0 (right column), Teff = 3500 K, and log g = 4.5. ODF calculations between 400 and 550 nm were performed with two different setups: using sub-binning proposed by Castelli & Kurucz (2003) (red) and using sub-binning optimized for cold stars (i.e., four sub-bins of equal size, blue). ODF calculations longward 550 nm were performed only with Castelli & Kurucz (2003) sub-binning.

Fig. A.3 RMS error of the ODF calculations in the 400 – –550 nm spectral domain. The errors are shown of the calculations with the Castelli & Kurucz (2003) sub-binning (blue) and optimized sub-binning (orange), computed relative to the high-resolution calculations.

Appendix B Interpolation Errors of Stellar Limb Darkening

In Sect. 4 we discussed the errors in limb-darkening curves associated with the interpolation of the power-2 law coefficients from the MARCS and PHOENIX metallicity grids at solar effective temperature (see Fig. 6). In Figs.B.1 and B.2, we show the corresponding errors for hotter (9000 K) and cooler (3500 K) stars, respectively. In addition, we present interpolation errors for several selected metallicity values in Fig. B.3.

We note that α and c are not independent (Maxted 2018). As a result, the dependence of the limb-darkening coefficients on metallicity contains wiggles (see Fig. 6 in the main text and Fig. B.2) arising from the fitting procedure. This, however, does not imply any abrupt changes in the limb darkening curves since the wiggles in one of the parameters are compensated by the wiggles in another.

Fig. B.1 Same as for Fig. 6, but for a star with effective temperature of 9000 K.

Fig. B.2 Same as for Fig. 6, but for a star with effective temperature of 3500 K.

Fig. B.3 Same as for the lower panels in Fig. B.2, but for selected values of metallicity.

Appendix C Biases of Limb-darkening Coefficients

The fitting of the limb darkening curve by any parameterized limb-darkening law is bound to introduce biases. Fig. C.1 gives an example of the biases arising from the representation of the limb-darkening curve by the power-2 law (see Eq. (3)). The biases are given in terms of the h₁ and h₂ coefficients (see Sect. 5.3). The h₁ coefficient is defined by Eq. (4). Since our limb-darkening curves stop at µ = 0.01, instead of using the h₂ coefficient given by Eq. (5), we introduce the $h_2^{*}$ coefficient:

$$h_2^{*}=\frac{I_{\text{pb}}0.5}{I_{\text{pb}}1.0}-\frac{I_{\text{pb}}0.1}{I_{\text{pb}}1.0}-c{2}^{-\alpha }-{100}^{-\alpha }.$$(C.1)

We compare the coefficients derived directly from the actual intensity profiles with the coefficients derived from the c and α coefficients extracted from the power-2 law fit. Figure C.1 shows that the h₁ coefficient is stable and can be obtained from the fitting procedure (see left panels of Fig. C.1). However, the fitting procedure cannot return accurate values of $h_2^{*}$ (see right panels of Fig. C.1).

Fig. C.1 Dependence of the h1 (left panels) and $h_2^{*}$ (right panels, see text for the definition) coefficients on metallicity (top panels) and effective temperature (bottom panels). Dashed lines indicate the solar value of the metallicity (M/H = 0, top panels) and the effective temperature of 5800 K (bottom panels). The h1 and $h_2^{*}$ values are shown, obtained from the actual intensity profiles (red curves) and the power-2 law fit (blue curves).

Appendix D Model Atmosphere Structure

An example of the model atmosphere structure from set 2 for a star with M/H = 0.0, T_eff = 5800 K, and log g = 4.5 is presented in Table D.1. All model structures from both set 1 and set 2 are available at the CDS. For each model we also provide the stellar parameters for which the model was computed, the values of the mixing-length parameter and overshoot, abundances, and the convergence criteria (maximum temperature deviation).

Appendix E Stellar Limb Darkening

Appendix F Parameterized Stellar Limb Darkening

References

All Tables

All Figures

	Fig. 1 Representation of metallicity and effective temperature grids. Black dots show the grid from this paper. Blue dots depict the MARCS grid from Gustafsson et al. (2008) and red dots represent the PHOENIX grid from Husser et al. (2013).
In the text

Fig. 2 Mixing-length parameter dependence on the stellar effective temperature for three values of metallicities (M/H = [−0.5, 0.0, 0.5]) and a surface gravity of log g = 4.5. The solar mixing-length parameter (α⊙ = 1.25) is taken from Kurucz (1993). Different colors correspond to different metallicities, as described in the legend. Horizontal dashed lines show the α values that are used in computations outside the region of effective temperatures for which the mixing-length behavior was derived by Viani et al. (2018).

In the text

	Fig. 3 Temperature structure of several models from set 1 plotted as a function of Rosseland optical depth. The three panels show models when one of the three stellar parameters is varied and the other two are fixed. Fixed parameters are presented in the title to each panel and the varying parameter in the legend. We plot every fifth model from the effective temperature and metallicity grids.
In the text

	Fig. 4 Normalized broadband transmission curves of the exoplanetary space missions, Kepler, TESS, CHEOPS, and PLATO, considered in this study (left panel) and the corresponding limb darkening for a star with M/H = 0.0, Teff = 5800 K, and log g = 4.5 (right panel). Different colors indicate different filters described in the legend.
In the text

	Fig. 5 Limb-darkening coefficients of the power-2 law on the stellar metallicity (top panel) and effective temperature (bottom panel) grids in the Kepler passband. The different colors depict the coefficients (α and c). Fixed parameters are labeled on each panel.
In the text

Fig. 6 Comparison of the limb darkening at our, MARCS, and PHOENIX grids of stellar parameters in the Kepler passband. Top panels: limb darkening coefficients of the power-2 law as a function of stellar metallicity. The vertical dashed lines show the metallicity grid points of the MARCS and PHOENIX models (left and right panels, respectively). Bottom panels: ratios between the limb darkening curves produced from coefficients directly calculated on our fine metallicity grid and those produced from the same coefficients but selected on the MARCS and PHOENIX metallicity grid points, and linearly interpolated to our fine metallicity grid (left and right panels, respectively). Ratios corresponding to different metallicity values are indicated by different colors. Coefficients and the limb darkening ratios are presented for the fixed Teff = 5800 K and log g = 4.5.

In the text

Fig. 7 Solar spectrum (top panel) and a comparison of solar limb-darkening computations with solar measurements by Neckel & Labs (1994) (bottom panel). Stars in both panels depict the solar measurements, and lines represent computations. Measurements and computations of solar spectrum in the top panel are shown at the disk center. Spectra at selected view angles, namely the µ = 0.8, 0.6, 0.4, and 0.2, normalized to the solar spectrum at the disk center, are shown in the bottom panel. The different line types marked in the legend to the figure correspond to the computation with different approximations: RE, with convection (α = 1.25), and with convection and overshoots. Intensity is measured in [erg cm−3 s−1 sr−1].

In the text

Fig. 8 Comparison of the MPS-ATLAS limb-darkening computations to that by CB11 (for the Kepler passband, left panels) and by C17 (for the TESS passband, right panels) for a star with solar fundamental parameters. Top panels: describe the limb darkening curves and the bottom panels present their ratios. Calculations with MPS-ATLAS (LD1 and LD2) are given, taking the overshoot into account (blue) as well as ignoring it (red). See text for the explanation of the abbreviations.

In the text

	Fig. 9 Ratios between limb darkening from LD2, CB11, and Sing (2010) (Sing) to the NL_Kepler limb darkening. Different colors correspond to different ratios (see legends).
In the text

	Fig. 10 Comparison of the limb darkening coefficients. The stars symbols show the difference between observations (Maxted 2018) and computations done with STAGGER (blue) and MPS-ATLAS (red) codes in Kepler passbands. The horizontal dashed line represents the zero level.
In the text

	Fig. A.1 Ratios of intensity values in the Kepler passband computed using ODF to that using high-resolution opacity (resolving power R=500000). Calculations were performed for several pairs of metallicities and effective temperatures, but the same value of the surface gravity (log g = 4.5).
In the text

Fig. A.2 Ratio of emergent intensities at µ = 1.0 (top panels), µ = 0.5 (middle panels), and µ = 0.2 (bottom panels), calculated using the ODF approach to intensities calculated on a fine spectral grid and then averaged over bins used in the ODF calculations. The ratios are shown for the models with M/H = 0.5 (left column) and M/H = 1.0 (right column), Teff = 3500 K, and log g = 4.5. ODF calculations between 400 and 550 nm were performed with two different setups: using sub-binning proposed by Castelli & Kurucz (2003) (red) and using sub-binning optimized for cold stars (i.e., four sub-bins of equal size, blue). ODF calculations longward 550 nm were performed only with Castelli & Kurucz (2003) sub-binning.

In the text

	Fig. A.3 RMS error of the ODF calculations in the 400 – –550 nm spectral domain. The errors are shown of the calculations with the Castelli & Kurucz (2003) sub-binning (blue) and optimized sub-binning (orange), computed relative to the high-resolution calculations.
In the text

	Fig. B.1 Same as for Fig. 6, but for a star with effective temperature of 9000 K.
In the text

	Fig. B.2 Same as for Fig. 6, but for a star with effective temperature of 3500 K.
In the text

	Fig. B.3 Same as for the lower panels in Fig. B.2, but for selected values of metallicity.
In the text

Fig. C.1 Dependence of the h1 (left panels) and $h_2^{*}$ (right panels, see text for the definition) coefficients on metallicity (top panels) and effective temperature (bottom panels). Dashed lines indicate the solar value of the metallicity (M/H = 0, top panels) and the effective temperature of 5800 K (bottom panels). The h1 and $h_2^{*}$ values are shown, obtained from the actual intensity profiles (red curves) and the power-2 law fit (blue curves).

In the text