© The Authors 2023

1 Introduction

Stars show a decrease in specific intensity from their centre to their limb. This is because opacity causes the sightlines of the observer to penetrate less far when they enter the photosphere at an angle. This effect, known as limb-darkening, occurs both for real stars and for theoretical stellar atmosphere models, and is important in any situation where a star is spatially resolved, such as the study of eclipsing binaries, transiting planetary systems, stellar interferometry, microlensing, and spectral line profiles.

The specific intensity profiles generated by model atmospheres with spherical symmetry differ from their plane-parallel equivalents mainly in the region near the limb, where the former shows a sharp drop while the latter presents a finite specific intensity. This important difference is due to the fact that in the plane-parallel case, the curvature of the atmosphere is neglected. Neglecting curvature effects is justified in stars with low atmospheric scale heights, but not for cool giants or supergiants. In plane-parallel model calculations, the medium is infinite, there is no natural boundary condition, and the specific intensity is not zero at the limb but depends on the local temperature in this region. However, in the case of models with spherical symmetry, the temperature tends to zero at small optical depths (see Eq. (7.182) in Mihalas 1978). A concise explanation for this crucial difference can be found in Mihalas (1978, p. 246–247); see also Larson (1969). In the theoretical-observational context of flux and polarisation calculations, Kostogryz et al. (2017) report the superiority of spherical models over plane-parallel models in some configurations, such as systems with grazing eclipses, transits with Earth-size planets, or for hotter planet host stars (effective temperatures T_eff > 6000 K).

In Claret & Southworth 2022 (hereafter Paper I) we showed the superiority of the power-2 law (Hestroffer 1997) in terms of the quality of fits of limb-darkening coefficients (LDCs) to theoretical predictions over the other two-parameter laws in the case of plane-parallel model atmospheres. However there remain few sources of LDCs calculated for the power-2 law, especially using spherical model atmospheres. The exceptions are LDCs for DA, DB, and DBA white dwarfs using plane-parallel models (Claret et al. 2020), solar-type stars using 3D models (Maxted 2018), extensive tables using the ATLAS plane-parallel models (Paper I), and LDCs for the CHEOPS space mission, where spherical models are adopted (Claret 2021).

The aim of this work is to provide users with LDCs for the power-2 law based on spherical model atmospheres, covering a wide range of T_eff values, surface gravities, and passbands. The passbands adopted here are Strömgren uυby, Johnson–Cousins UBVRIJHK, and SDSS ugriz photometric systems, plus those for the space missions Gaia, Kepler, TESS, and CHEOPS. The results of the present study supersede those of Claret (2021) for the CHEOPS mission.

The structure of the paper is as follows. Section 2 is dedicated to the description of the spherical models (PHOENIX-COND) and the numerical details of the calculations of the LDCs. A comparison between the power-2 and four-parameter LDCs is presented in Sect. 3. Section 4 gives a comparison between the plane-parallel and spherical LDCs, and Sect. 4 summarises our findings.

2 The spherical atmosphere models and numerical methods

In this paper, we use specific intensity profiles from the PHOENIX-COND spherical model atmospheres. The main characteristics of these models are described in Husser et al. (2013); see also Claret et al. (2012). We adopt a Levenberg-Marquardt least-squares minimisation method to compute the LDCs. The specific intensities for the uυby, UBVRIJHK, ugriz, Gaia, Kepler, TESS, and CHEOPS photometric systems were integrated using the following equation:

(1)

where h is Planck’s constant, c is the speed of light in vacuum, λ is the wavelength in Å, and µ is given by cos(γ), where γ is the angle between the line of sight and the outward surface normal. I_a(µ) is the specific intensity for a given passband a, I(λ,µ) is the monochromatic specific intensity, and S_a(λ) is the response function. For the uυby, UBVRIJHK, and ugriz photometric systems, the response function includes the transmission of one airmass of the Earth’s atmosphere. The passbands used were obtained from Spanish Virtual Observatory Filter Profile Service¹. For the specific case of Strömgren uυby, the response function was obtained from the Observatorio de Sierra Nevada, Granada, Spain (C. Cárdenas, priv. comm.). For JHK, the corresponding S(λ) were obtained from the Observatorio del Teide-IAC, Spain (Alonso et al. 1994).

The specific intensities generated by the spherical model atmospheres show a much more pronounced curvature near the limb than the plane-parallel models, making it much more difficult to obtain a good fit to them. For example, the four-parameter law is able to fit the full profiles well but only for some filters (see e.g. Fig. 1 in Claret et al. 2012). Some alternative methods have been proposed to better describe the specific intensity profiles for models with spherical symmetry, such as Claret & Hauschildt (2003), where the concept of quasi-spherical models was introduced.

Later, Wittkowski et al. (2004) introduced a more elaborate method: instead of truncating the models at a certain value of µ, the truncation was defined by searching for the maximum of the derivative of the specific intensity with respect to r, where . This point corresponds to τ_R ≈ 1.0, where τ_R is the Rosseland mean optical depth. Fortunately, the PHOENIX models were computed with sufficient points in the drop-off region to enable us to accurately determine the corresponding derivatives. Wittkowski et al. (2004) gave two ways to determine this critical point. The first method (M1) uses the point where the derivative is maximum, which is known as the point µ_cri, while the second method (M2) uses the midpoint between the point µ_cri and the point located at µ_cri–1 (G. Morello, priv. comm.). The specific intensity using M2 at this average point is larger than that obtained by the M1 method. Therefore, it gives rescaled profiles that are more similar to the plane-parallel models than those provided by M1. However, M2 does not accurately represent the specific intensity at the critical point as defined by Wittkowski et al. (2004). The corresponding critical points for both methods are characterised by µ_cri,1 and µ_cri,2. In the current work, we present calculations for both M1 and M2 for completeness.

The LDC calculations were performed adopting 823 PHOENIX stellar model atmospheres for solar metallicity, surface gravities log g from 0.0 to 6.0, T_eff values from 2300 K to 12 000 K, and microturbulent velocities V_ξ = 2.0 km s⁻¹. We adopted 900 equally spaced-points in µ, rather than the 78 original ones adopted in the PHOENIX models. In our previous papers on LDCs from spherical models, we used only 100 µ points, and so we have updated our calculations using 900 points for the photometric systems studied in this work.

The power-2 LD law introduced by Hestroffer (1997) is

(2)

where g and h are the corresponding LDCs. Given that the profiles generated by PHOENIX are much more complicated and difficult to adjust than in the plane-parallel approximation, we also calculated LDCs using the four-parameter law (Claret 2000) which is

(3)

where a_k are the four LDCs. The greater flexibility of the four-parameter law allows it to reproduce the specific intensity profiles to a higher precision than limb-darkening laws using only two LDCs.

We define a merit function, which measures the quality of the fit to a given set of I(µ) values as

(4)

where y_i is the model intensity at point i, Y_i is the fitted function at the same point, and N is the number of µ points.

Fig. 1 Angular distribution of the specific intensity for a model with Teff = 4500 K, log g = 4.5, [M/H] = 0.0, and Vξ = 2 km s¯1 for the TESS passband. The red line represents the specific intensity distribution and the red cross indicates µcri,1. The green line denotes the fitting adopting the power-2 law, while the green cross indicates the fitting at µcri,1. The blue line denotes the four-parameter law approach and the blue cross indicates the fitting at µcri,1. This plot is for case M1, and the specific intensity profiles have not been re-scaled.

Fig. 2 Same as in Fig. 1, except for case M2 and µcri,2 instead of case M1 and µcri,1.

3 Results

An extensive comparison of the different limb-darkening laws was carried out in Paper I. It was found that the power-2 law produced better fits to the specific intensities from plane-parallel model atmospheres than the other two-parameter laws (quadratic, logarithmic, and square-root). Given this result and the increased complexity of the specific intensities from spherical model atmospheres, it was not thought necessary to repeat that analysis here. However, we find it useful to investigate the relative success of the power-2 and four-parameter laws, as this information is helpful to users of the LDCs.

The resulting χ² values are generally small and are dominated by the difficulty of fitting the µ_cri,1 and µ_cri,2 points, which is most conspicuous for the power-2 laws in the M1 case. Figures 1–4 show the specific intensity profiles for the M1 and M2 cases using both the power-2 (Eq. (2)) and four-parameter (Eq. (3)) law. For the sake of clarity, only the regions near the limb are shown, because the fits for µ > µ_cri are almost perfect. For main-sequence stars (Figs. 1 and 2), and in both the M1 and M2 cases, the four-parameter law (Eq. (3)) provides better fits than those provided by the power-2 law (Eq. (2)). The situation is the same but even clearer for the giant stars (Figs. 3 and 4). It is important to emphasise that one of the main differences between the plane-parallel models and the spherical ones is in the dropoff region, and so it is necessary to characterise this region well through the LDCs. The best test of the spherical models would be achieved using Eq. (3) and case M1.

The relative quality of the fits can be computed using the following expression:

(5)

Figures 5–10 show the quantity log α for the M1 case, a range of T_eff and log g values, and for the Gaia G_BP, G, G_RP, Kepler, TESS, and CHEOPS passbands. It can be seen that the χ² provided by Eq. (3) is typically a factor of 8 ± 2 smaller than the χ² given by Eq. (2), depending on the effective wavelength of the passband. We also see that there is a common maximum in all passbands centred on log T_eff ≈ 3.65 (≈4500K), where the four-parameter law performs much better than the power-2 law. A comparison of the effects for the SDSS u and z passbands (not shown) indicates that the ratio α is larger at shorter wavelengths, which is in agreement with the results for the Gaia G_BP pass-band (Fig. 5). It is also apparent that the power-2 law has a better relative performance at higher T_eff values, but this is still poorer than the four-parameter law.

In Figs. 5–10, we trace the relative performance for two surface gravities: log g = 6.0 for the most compact stars and log g = 3.5 to represent subgiants. This shows how the relative quality of the fits depends on the surface gravity of a star. In general, the power-2 law is better at higher surface gravities, but is still not as good as the four-parameter law. This is the same in all passbands except for the interval 3.4 ⪆ log T ⪆ 3.55 in the G_BP band. For models in radiative equilibrium, there is also a tendency for the relative quality of the fits to not depend strongly on the surface gravity.

Fig. 3 Same as in Fig. 1, except for log g = 2.5 instead of log g = 4.5.

Fig. 4 Same as in Fig. 2, except for log g = 2.5 instead of log g = 4.5.

Fig. 5 Comparison of the quality of fits provided for the Gaia GBP passband. The quantity α is the ratio between the χ2 for the power-2 law and the χ2 for the four-parameter law (both for case M1), for all 823 models. The red line traces the results for the most compact models (log g = 6.0) whilst the green line indicates the results for subgiant ones (log g = 3.5).

Fig. 6 Same as in Fig. 5, but for the Gaia G passband.

Fig. 7 Same as in Fig. 5, but for the Gaia GRP passband.

Fig. 8 Same as in Fig. 5, but for the Kepler passband.

Fig. 9 Same as in Fig. 5, but for the TESS passband.

Fig. 10 Same as in Fig. 5, but for the CHEOPS passband.

Fig. 11 Comparison between the power-2 g (upper panel) and h (lower panel) LDCs in the case of plane-parallel model atmospheres (solid lines) and spherical model atmospheres for cases M1 (dotted lines) and M2 (dashed lines).

4 Comparison between plane-parallel and spherical LDCs

A comparison between the power-2 LDCs predicted by plane-parallel and spherical model atmospheres is useful. For the plane-parallel model atmospheres, we adopted the results from Paper I, which were calculated using the ATLAS code. For the spherical model atmospheres, we used the results from the current work for both cases M1 and M2.

Figure 11 shows the variation of the g and h LDCs as a function of T_eff, for log g = 4.0, solar metallicity, and V_ξ = 2.0 km s⁻¹. Cases M1 and M2 show much closer agreement than the ATLAS LDCs, in particular at lower T_eff values where spherical model atmospheres are more reliable. The anti-correlation between g and h is also clear, particularly at low T_eff values for M1 and M2, and at 6000 < T_eff < 8000 K for the ATLAS LDCs.

An important use of the LDCs presented in the current work is in calculating the light curves of planetary transits (see Fig. 12). We explored this by using the JKTEBOP code (Southworth 2013) to calculate the transit light curve of a system similar to HAT-P-7 (Pál et al. 2008), using the transit parameters from Southworth (2011). This system was chosen because the small ratio of the radii of the planet and star gives a good spatial sampling of the limb-darkening over the disc of the star. To avoid interpolating LDCs to specific atmospheric parameters, we adopted fixed values of T_eff = 6000 K, log g = 4.0, solar metallicity, and V_ξ = 2.0 km s⁻¹. The differences between the ATLAS and PHOENIX LDCs cause a change in the transit shape that reaches approximately 70 ppm around second and third contact. In principle, this change is large enough to be measured from space-based missions such as Kepler and TESS, especially when many transits can be averaged together, but correlations between the measured LDCs and other parameters of the system will reduce the significance of the signal. The M1 and M2 LDCs give much more similar transit light curves – the greatest deviation is only 4 ppm –, which suggests that the choice of which to use is unimportant.

Fig. 12 Top panel: comparison between the transit shapes using the power-2 LDCs for the plane-parallel (ATLAS) and spherical (cases M1 and M2) model atmospheres (labelled). Lower panels: differences between the transit shapes in each case (labelled).

5 Summary and final remarks

We computed LDCs for the power-2 law using specific intensities from spherical model atmospheres, as such results were only previously available for a small range of stellar parameters. To this end, we used the PHOENIX-COND models for T_eff values from 2300K to 12000K, logy values from 0.0 to 6.0, solar metallicity, V_ξ = 2.0 km s⁻¹, and with 900 µ points. We also computed LDCs for the four-parameter law, which supersede our previous calculations based on only 100 µ points. For both laws, we computed coefficients for the uυby, UBVRIJHK, and SDSS ugriz passbands, and for the Gaia, Kepler, TESS, and CHEOPS photometric systems. Two approaches were used to define the limb of the star based on the drop-off of intensity in the spherical model atmosphere predictions. We emphasise that the µ values have been rescaled to remove those beyond the limb of the star, and so this need not be done by users of the LDCs.

We performed comparisons between the merit function χ² for the power-2 and the four-parameter laws, finding that the four-parameter law is superior in all cases to the power-2 law. The performance of the power-2 law is better at higher temperatures and at redder wavelengths. Additional calculations of LDCs for other photometric systems not covered in this paper are available upon request to the first author.

We recommend that the four-parameter law be used whenever high precision is needed for specific intensity profiles. If a less complex law is preferred, we recommend the power-2 law as it reproduces the theoretical intensity profiles to a higher precision than the other two-parameter laws (see Paper I); we also recommend not using the linear and quadratic laws because they provide the least good fit to theoretical specific intensity profiles. For modelling eclipsing binary star and transiting planetary systems, the choice of which limb-darkening law to use has a relatively small effect on the best-fitting parameter values (Southworth 2023).

The current work, together with Paper I, provides three sets of LDCs for any given set of stellar parameters and passbands, within those covered in this work. These are LDCs from planeparallel ATLAS model atmospheres (Paper I) and LDCs from spherical model atmospheres (this work) calculated using methods M1 and M2. Those from methods M1 and M2 are very similar, whereas the LDCs from ATLAS are moderately different. We show that these LDCs could be tested using high-precision photometry of transits or eclipses. In our experience, planetary transits are the best approach for this because the small size of the planet relative to the star provides a finer spatial sampling of the star’s specific intensity profile, but the size of the differences between the transit light curves from different LDCs is not large.

Acknowledgements

We thank the anonymous referee for his/her helpful comments that have improved the manuscript. The Spanish MEC (AYA2015-71718-R, ESP2017-87676-C5-2-R, PID2019-107061GB-C64, and PID2019-109522GB-C52) is gratefully acknowledged for its support during the development of this work. A.C. acknowledges financial support from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033. This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France, of NASA’s Astrophysics Data System Abstract Service and of SVO Filter Profile supported from the Spanish MINECO through grant AYA2017-84089.

Appendix A Description of Tables 1–12 (available at the CDS)

References

All Tables

All Figures

Fig. 1 Angular distribution of the specific intensity for a model with Teff = 4500 K, log g = 4.5, [M/H] = 0.0, and Vξ = 2 km s¯1 for the TESS passband. The red line represents the specific intensity distribution and the red cross indicates µcri,1. The green line denotes the fitting adopting the power-2 law, while the green cross indicates the fitting at µcri,1. The blue line denotes the four-parameter law approach and the blue cross indicates the fitting at µcri,1. This plot is for case M1, and the specific intensity profiles have not been re-scaled.

In the text

	Fig. 2 Same as in Fig. 1, except for case M2 and µcri,2 instead of case M1 and µcri,1.
In the text

	Fig. 3 Same as in Fig. 1, except for log g = 2.5 instead of log g = 4.5.
In the text

	Fig. 4 Same as in Fig. 2, except for log g = 2.5 instead of log g = 4.5.
In the text

	Fig. 5 Comparison of the quality of fits provided for the Gaia GBP passband. The quantity α is the ratio between the χ2 for the power-2 law and the χ2 for the four-parameter law (both for case M1), for all 823 models. The red line traces the results for the most compact models (log g = 6.0) whilst the green line indicates the results for subgiant ones (log g = 3.5).
In the text

	Fig. 6 Same as in Fig. 5, but for the Gaia G passband.
In the text

	Fig. 7 Same as in Fig. 5, but for the Gaia GRP passband.
In the text

	Fig. 8 Same as in Fig. 5, but for the Kepler passband.
In the text

	Fig. 9 Same as in Fig. 5, but for the TESS passband.
In the text

	Fig. 10 Same as in Fig. 5, but for the CHEOPS passband.
In the text

	Fig. 11 Comparison between the power-2 g (upper panel) and h (lower panel) LDCs in the case of plane-parallel model atmospheres (solid lines) and spherical model atmospheres for cases M1 (dotted lines) and M2 (dashed lines).
In the text

	Fig. 12 Top panel: comparison between the transit shapes using the power-2 LDCs for the plane-parallel (ATLAS) and spherical (cases M1 and M2) model atmospheres (labelled). Lower panels: differences between the transit shapes in each case (labelled).
In the text