Limb darkening laws for two exoplanet host stars derived from 3D stellar model atmospheres

Limb darkening laws for two exoplanet host stars derived from 3D stellar model atmospheres

Comparison with 1D models and HST light curve observations
W. Hayek Astrophysics Group, School of Physics, University of Exeter, Stocker Road, Exeter, EX4 4QL
[hayek,sing,fpont] Max Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany
   D. Sing Astrophysics Group, School of Physics, University of Exeter, Stocker Road, Exeter, EX4 4QL
   F. Pont Astrophysics Group, School of Physics, University of Exeter, Stocker Road, Exeter, EX4 4QL
   M. Asplund Research School of Astronomy & Astrophysics, Cotter Road, Weston Creek 2611, Australia Max Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany
July 4, 2019
Key Words.:
Stars: atmospheres – Stars: individual: HD 209458, HD 189733 – Planets and satellites: individual: HD 209458b, HD 189733b

We compare limb darkening laws derived from 3D hydrodynamical model atmospheres and 1D hydrostatic MARCS models for the host stars of two well-studied transiting exoplanet systems, the late-type dwarfs HD~209458 and HD~189733. The surface brightness distribution of the stellar disks is calculated for a wide spectral range using 3D LTE spectrum formation and opacity samplingthanks: Full theoretical spectra for both stars are available in electronic form at the CDS via anonymous ftp to ( or via, as well as at We test our theoretical predictions using least-squares fits of model light curves to wavelength-integrated primary eclipses that were observed with the Hubble Space Telescope (HST). The limb darkening law derived from the 3D model of HD 209458 in the spectral region between  Å and  Å produces significantly better fits to the HST data, removing systematic residuals that were previously observed for model light curves based on 1D limb darkening predictions. This difference arises mainly from the shallower mean temperature structure of the 3D model, which is a consequence of the explicit simulation of stellar surface granulation where 1D models need to rely on simplified recipes. In the case of HD 189733, the model atmospheres produce practically equivalent limb darkening curves between  Å and  Å, partly due to obstruction by spectral lines, and the data are not sufficient to distinguish between the light curves. We also analyze HST observations between  Å and  Å for this star; the 3D model leads to a better fit compared to 1D limb darkening predictions. The significant improvement of fit quality for the HD 209458 system demonstrates the higher degree of realism of 3D hydrodynamical models and the importance of surface granulation for the formation of the atmospheric radiation field of late-type stars. This result agrees well with recent investigations of limb darkening in the solar continuum and other observational tests of the 3D models. The case of HD 189733 is no contradiction as the model light curves are less sensitive to the temperature stratification of the stellar atmosphere and the observed data in the  Å -  Å region are not sufficient to distinguish more clearly between the 3D and 1D limb darkening predictions.

1 Introduction

The surface brightness of a star varies across its visible disc as the optical path length along the line of sight depends on the viewing angle onto the atmosphere. Closer to the limb, radiative flux emerges from higher and cooler regions in the stellar photosphere, therefore causing “limb darkening”. This angular dependence is commonly described with an analytical, linear or nonlinear law, which forms an essential ingredient for analyzing the light curves of transiting exoplanets and eclipsing binary stars. The occulting body indirectly samples radiative intensities in the observer’s direction as it travels across the stellar disk, leaving a clear signature in the shape of the light curve. A uniform intensity distribution would result in an almost rectangular shape with a flat-bottom light curve; the limb darkening effect causes shallower slopes during the ingress and egress phases and a curved signal in the center of the transit.

Light curve analyses rely on theoretical predictions for limb darkening in general. With few exceptions, these are still mostly derived from classical 1D hydrostatic stellar model atmospheres, for which extensive grids are available (see, e.g., Claret 2000; Sing 2010). Very high signal-to-noise light curve observations were obtained with the Hubble Space Telescope (HST) for HD 209458 (Brown et al. 2001; Knutson et al. 2007) and for HD 189733 (Pont et al. 2008; Sing et al. 2011). The data quality enables a direct fit of at least a low-order limb darkening law (see, e.g., the discussion in Southworth 2008) along with the other transit parameters, and therefore an assessment of the accuracy of model predictions for limb darkening. It turns out that using 1D models leads to substantially higher residuals in the light curve fits compared to directly fitted limb darkening laws. Moreover, these residuals exhibit deviations in the ingress and egress phases of the transit that are characteristic for incorrect limb darkening (Knutson et al. 2007). Accurate independent model predictions for limb darkening are nevertheless generally preferable to fitting a law to the data in order to reduce the number of free parameters in the fit and avoid the lower accuracy of the linear or quadratic laws that are typically used in the method. This would also eliminate the wavelength-dependent degeneracy of limb darkening with transit depth and thus with the important planet-to-star radius ratio on which transit spectroscopy of planetary atmospheres is based.

As the limb darkening curve is mainly determined by the photospheric temperature gradient, the systematic residuals in the light curve fits of Knutson et al. (2007) for HD 209458b point towards shortcomings in the structure of 1D stellar model atmospheres. Claret (2009) further investigated the case, comparing limb darkening coefficients derived from various 1D model atmospheres and different analytical laws with coefficients that resulted from empirical fits by Southworth (2008), and concluded that current 1D atmosphere models of HD 209458 cannot produce consistent results.

The Sun has traditionally been the most important bench mark for testing the realism of model atmospheres as, e.g., the surface brightness distribution can be readily measured in great detail (Neckel & Labs 1994). Comparisons with different solar model atmospheres by Asplund et al. (2009) revealed that the latest generation of 3D hydrodynamical models is capable of satisfying this important test with very good accuracy (see their Fig. 2), while the limb darkening predicted by 1D hydrostatic MARCS models is too strong; Sing et al. (2008) found a similar result for a solar 1D ATLAS model. 3D models take the effect of convective motions in the surface granulation explicitly into account and are thus able to reproduce the solar atmosphere with a higher degree of realism; see Nordlund et al. (2009) for a detailed description of the physics and methods.

We compute detailed theoretical surface intensity distributions from 3D hydrodynamical model atmospheres of HD 209458 and HD 189733 across a large wavelength range and derive limb darkening coefficients for the spectral band between  Å and  Å covered by the HST STIS G430L grating, where limb darkening is more sensitive to the photospheric temperature gradient, using least-squares fits of a nonlinear law to the numerical calculations. We additionally investigate the region between  Å and  Å using ACS HRC G800L data in the case of HD 189733, as our STIS observations include only a small fraction of the transit. The models and the choice of stellar parameters are described in Sect. 2, followed by a discussion of the 3D radiative transfer computations needed for obtaining theoretical spectra in Sect. 3. We compare our predictions of 3D limb darkening with the results from 1D MARCS models in Sect. 4 and test their accuracy using fits to HST transiting light curves in Sect. 5. The paper concludes with a summary of our results in Sect. 6. Limb darkening coefficients for a selection of various standard broadband filters and instrument sensitivities are presented in Appendix A.

2 Stellar model atmospheres

Figure 1: Temperature distributions (gray shades) on surfaces with constant optical depth at  Å for arbitrary snapshots of the 3D models for HD 209458 (left) and HD 189733 (right). The spatial and temporal mean temperature stratification of the 3D models (green solid lines) is compared to 1D MARCS models (red solid lines) with the same stellar parameters. The dashed lines bracket the approximate height region between and , from which the dominant contribution to the continuum surface brightness between the disk center at and the limb near is emitted.

Theoretical models of stellar atmospheres have traditionally been based on the approximation of horizontal homogeneity and plane parallel or spherically symmetric geometry. The problem is thereby reduced to one dimension, which enables very fast computation on modern computers and a very detailed treatment of radiative energy transfer. Large grids of stellar atmospheres were created, such as the Kurucz grid (Kurucz 1996), the MARCS grid (Gustafsson et al. 2008) and the NextGen grid (Hauschildt et al. 1999).

However, the atmospheres of late-type stars are known to develop a convection zone in the envelope that reaches the stellar surface, producing a granulation pattern that causes strong horizontal inhomogeneities in the atmospheric structure, which strongly influences the outgoing radiation field (see the discussion in, e.g., Nordlund et al. 2009, for the important case of the Sun). 3D time-dependent radiation-hydrodynamical model atmospheres take these convective motions explicitly into account.

We construct two 3D models to represent the exoplanet host stars HD 209458 and HD 189733 using the Stagger Code (Nordlund & Galsgaard 1995). The coupled equations of compressible hydrodynamics and time-independent radiative transfer are solved simultaneously in a plane-parallel box with a resolution of grid points, covering granules at a time to represent stellar surface convection. The dependence of continuous and spectral line opacities on wavelength and atmospheric height is approximated using the opacity binning method (Nordlund 1982). Opacities were taken from the same opacity package that was used for the spectrum calculations (see Sect. 3) and sorted into 12 groups as a function of wavelength and Rosseland optical depth of their monochromatic optical surfaces. Based on the spatial and temporal mean stratification of the 3D model, intensity-mean and Rosseland mean opacities are computed for each group. They are subsequently combined into a single opacity using exponential bridging as a function of optical depth; see Skartlien (2000) for a detailed description of the method and Hayek et al. (2010) for a summary of opacity sources. The chemical composition is based on the solar abundances of Asplund et al. (2005). The simulations reach a horizontal extent of (HD 209458) and (HD 189733) on equidistant axes with periodic boundaries. On the vertical axis, which has finer resolution around the continuum optical surface of the models,  Mm (HD 209458) and  Mm (HD 189733) are covered with open boundaries to allow free inflow and outflow of stellar gas. This corresponds to an optical depth range of about at  Å for both models, thus including all atmospheric layers that are relevant for spectrum formation above the optical surface and for surface granulation below. The simulations assume magnetically quiet convection, which ignores the effects of magnetic fields on plasma dynamics. This simplification is valid for stars like the Sun or HD 209458, where this activity should only weakly influence the outgoing radiation field in most parts of the spectrum. The case of HD 189733, where transit light curves show the existence of numerous star spots (see Sing et al. 2011), may require more scrutiny in the future with larger models that allow the inclusion of star spots.

The 1D atmosphere models are based on interpolations in the MARCS grid; see Gustafsson et al. (2008) for a description of the physics, input data and methods.

In order to test the consistency of stellar parameters, we perform LTE (local thermodynamic equilibrium) abundance analyses with lines of neutral and ionized iron on high resolution, high signal-to-noise HARPS spectra111Based on data obtained from the ESO Science Archive Facility.. The observed equivalent widths, derived from least-squares fits of Voigt profiles to the data, are compared to synthetic LTE line profiles that we compute with the SCATE line formation code (for a description see Hayek et al. 2011). SCATE uses the same opacity package on which the MARCS grid and the 3D model atmospheres are based. The code was also used to compute the intensity spectra needed for deriving the limb darkening laws (see Sect. 3), ensuring complete consistency in our investigation.

We base our abundance analysis on the Fe I line list of Asplund et al. (2000), with line data taken from Blackwell et al. (1995) and Holweger et al. (1995), as well as on the compilation of Fe II laboratory data of Meléndez & Barbuy (2009). The effective temperature is constrained by removing trends of the Fe I abundance with excitation potential , where low excitation lines exhibit the strongest sensitivity to temperature due to the excitation-ionization balance. The correct surface gravity is found by requiring matching abundances of Fe I and Fe II through LTE ionization equilibrium, where Fe II lines exhibit the strongest sensitivity. It is well known that departures from LTE affect absorber populations in stellar atmospheres, in particular in the case of neutral iron. These effects should be small for ionized iron, which thus yields a more reliable abundance measurement (see the review of Asplund 2005). For our goal of deriving theoretical limb darkening coefficients from LTE spectrum formation, achieving LTE ionization equilibrium for Fe I and Fe II still leads to consistent results: the observed spectra of HD 209458 and HD 189733 constrain the atmospheric stratification of the models through their stellar parameters, while the same radiative transfer methods and opacity data are then used in turn to predict the stellar surface brightness distribution that enters the model light curves.

HD 209458
Data set
Fe I B 15
Fe I H 23
Fe II 16
HD 189733
Data set
Fe I B 11
Fe I H 16
Fe II 4

(1) Blackwell et al. (1995); (2) Holweger et al. (1995); (3) Meléndez & Barbuy (2009)

Table 1: Abundances of neutral and ionized iron lines, derived from 3D hydrodynamical and 1D hydrostatic model atmospheres.
HD 209458
[K] [cm s] [Fe/H] Reference Method
Mazeh et al. (2000) spectroscopic
Gonzalez et al. (2001) spectroscopic
Mashonkina & Gehren (2001) spectroscopic
Sadakane et al. (2002) spectroscopic
Heiter & Luck (2003) spectroscopic
Santos et al. (2004) spectroscopic
Ramírez & Meléndez (2004) photometric
Valenti & Fischer (2005) spectroscopic
Sousa et al. (2008) spectroscopic
Southworth (2009) light curve
Casagrande et al. (2011) photometric
adopted (3D model) spectroscopic
adopted (1D model) spectroscopic
HD 189733
[K] [cm s] [Fe/H] Reference Method
Bouchy et al. (2005) spectroscopic
van Belle & von Braun (2009) photometric
Southworth (2010) light curve
Casagrande et al. (2011) photometric
adopted (3D model) spectroscopic
adopted (1D model) spectroscopic
222 Quoted metallicities could not be corrected for differences in the assumed solar composition, as the latter was not always specified in the cited publications, leading to artificial offsets in [Fe/H].
Table 2: Literature values compared to our adopted effective temperature, surface gravity, and metallicity.

A large number of parameter determinations for the G-dwarf HD 209458 are available in the literature (see upper panel in Table 2). Most of the effective temperature measurements are based on spectroscopic methods similar to our determination. Inherent inaccuracies of the method (uncertainties typically reach the order of  K) produce some scatter with values ranging between approximately  K and  K. We obtain a flat distribution of Fe I abundances with increasing excitation level for a 3D model with time-averaged  K, which agrees well with the majority of spectroscopic results and with the recent photometric measurement of Casagrande et al. (2011). The 1D model yields a smaller effective temperature of  K. This deviation between the 3D and 1D MARCS results stems from their different temperature gradients and the inhomogeneities of the 3D model that arise from the convective motions. The left panel in Fig. 1 compares the spatial and temporal average temperature stratification of the 3D model as a function of optical depth at  Å (green line) with the corresponding - relation of a 1D MARCS model (red line) with the same parameters as the 3D model, which exhibits a steeper temperature gradient around the optical surface at . The gray shades in the plot show the temperature distribution of an arbitrary snapshot of the 3D simulation. The dashed lines mark the approximate region that is relevant for the surface brightness distribution across the stellar disk in the continuum (see Sect. 4).

Choosing a surface gravity of leads to reasonably consistent ionization equilibria for both the 3D and 1D models (upper panel of Table 1) and is in good agreement with the majority of literature results, considering a measurement uncertainty of at least in . Southworth (2009) found from a light curve analysis, while Casagrande et al. (2011) derived from a Hipparcos parallax. Asplund et al. (2005) recommend a solar iron abundance of , which is sufficiently close to our 3D and 1D Fe II results for HD 209458, given the uncertainties in and the rather small number of Fe II lines; we therefore adopt their solar composition for the model atmospheres and spectrum computations.

We obtain a microturbulence parameter  km s for the 1D case from the distribution of line-to-line abundances as a function of equivalent width; a MARCS model atmosphere with slightly lower microturbulence ( km s) is used for all 1D spectrum computations to ensure consistency with the opacity sampling tables (see the discussion in Sect. 3). 3D line formation calculations do not require microturbulence as its contribution to profile broadening is naturally produced by 3D model atmospheres through Doppler shifts in the convective plasma flow.

To our knowledge, three independent spectroscopic and photometric measurements of the effective temperature of the K-dwarf HD 189733 exist in the literature (see lower panel in Table 2), which deviate by up to about  K and are therefore in reasonable agreement considering the systematic effects of the different methods. We determine  K for our 3D model, which is coincidentally equal to the spectroscopic value of Bouchy et al. (2005) that was derived with a 1D ATLAS model. The 1D MARCS model produces a slightly lower effective temperature of  K.

We find good agreement for our 3D ionization equilibrium (see lower panel of Table 1) for the surface gravity of Bouchy et al. (2005), which is also compatible with the light curve fit of Southworth (2010) and the parallax-based result of Casagrande et al. (2011). Similar agreement is found for the 1D model with its effective temperature of  K. The accuracy of our determinations is limited due to the very low number of usable Fe II lines that we measured in the spectrum, but the stellar surface gravity is only of secondary importance for the atmospheric temperature gradient of the model. The solar iron abundance of Asplund et al. (2005) agrees again reasonably well with our measurement. We determine a microturbulence parameter of  km s; we therefore choose again a grid model with  km s for the 1D calculations.

3 LTE spectrum formation and limb darkening laws

Figure 2: Left: Theoretical average surface brightness distribution as a function of projection factor integrated over the spectrum between  Å and  Å, computed using our “standard” numerical resolution for an arbitrary snapshot of the 3D model of HD 189733 (crosses), and least-squares fits of the Claret (2000) nonlinear law (solid line) and a linear law (dotted line). The gray-shaded area indicates the central part of the stellar disk that is not reached during the transit. Right: Relative deviation of spline interpolations of theoretical limb darkening computed with direction angles (solid line) and with the full horizontal resolution of grid points (dotted line) from the interpolated “standard” setting.
Figure 3: Left: Relative deviation of the spline-interpolated limb darkening predictions for (solid line) and (dotted line) from the case with as a function of . Grey vertical lines mark the position of the 17 angles used in the “standard” resolution. Right: Relative deviation of the least-squares-fitted limb darkening laws from the computed data points for (crosses), (diamonds) and (triangles) as a function of . All calculations are based on the 3D model of HD 189733.
Figure 4: Left: Relative deviation of the spline-interpolated limb darkening predictions computed for a single snapshot of the 3D model for HD 189733 with reduced wavelength sampling of (solid lines) and (dotted lines) from the full sampling with . Each set of resampled wavelengths is shown in either case, resulting in two or four choices, respectively. Right: Relative deviation of the spline-interpolated limb darkening predictions computed for the first 10 snapshots of a time sequence with (solid line) and deviation of the predictions computed for the first and last 10 snapshots of the same sequence (dotted line). All calculations are based on the 3D model of HD 189733.
117 0.5999 1.7021 0.64738 0.65040
151 0.6479 1.8951 0.64747 0.64751
101 0.6548 1.9230 0.64748 0.64751
333 All calculations are based on the 3D model of HD 189733 and integrated over the wavelength band between  Å and  Å. Radiative fluxes are divided by to obtain the same normalization as the limb darkening laws.
Table 3: Comparison of 4th order limb darkening coefficients, goodness of fit, radiative fluxes from the nonlinear law and numerical quadrature, and their relative deviation for various choices of .

We use the SCATE line formation code to compute surface intensities for all models at different angles with respect to the vertical axis, corresponding to varying positions on the stellar disk. The calculations are based on continuous opacities from Gustafsson et al. (1975) with updates by R. Trampedach (priv. comm.), as well as sampled line opacities from B. Plez (priv. comm.); see also Gustafsson et al. (2008). The chemical composition assumes the solar abundances of Asplund et al. (2005). In order to reduce the computational expense of full spectrum calculations in 3D geometry, all opacities were precomputed and tabulated as a function of wavelength, pressure and temperature using the approximation of local thermal equilibrium (LTE) and neglecting the effects of scattering. This approach precludes a detailed treatment of Doppler shifts, contrary to our 3D calculations of iron lines. We use opacity sampling tables with a microturbulence parameter of  km s to approximate its contribution to spectral line broadening. The inaccuracies of this method should be small as we only investigate broadband-integrated spectra. Moreover, our wavelength resolution is only sufficient for a statistical representation of the spectra and does not resolve the profile shapes of most lines.

A representation of a stellar atmosphere with a plane-parallel stratification inherently limits the realism of surface brightness calculations near the limb, where a sharp drop-off should appear as soon as the line of sight no longer reaches below the atmosphere. Claret & Hauschildt (2003) found such an intensity drop using 1D models with spherical symmetry. The position angle of the drop depends on the stellar parameters, which influence the pressure scale height of the atmosphere, as well as on wavelength through the gas opacity. A simple estimate for can be obtained by assuming a transparent atmosphere of one pressure scale height, which leads to the expression


where is the stellar radius and is the pressure scale height. Adopting and from van Belle & von Braun (2009), as well as choosing  km and  km at the surface as predicted by our 3D simulations, we obtain a drop-off at for both stars. This result is in good agreement with Claret & Hauschildt (2003), who find for a star with  K and . The lower surface gravity should yield a larger atmospheric pressure scale height compared to HD 189733, which has a very similar effective temperature. Our 3D models with plane-parallel stratification should therefore predict surface brightness distributions with sufficient realism for light curve analyses, except for some extreme cases of grazing transits.

3.1 Analytical limb darkening laws

SCATE yields monochromatic surface intensities at wavelength for a range of discrete positions on the stellar disk, assuming plane-parallel geometry as the underlying 3D model atmospheres, with the projection factor for angle in radial direction off the disk center. The resulting specific intensities are averaged over direction angle , horizontal dimensions and at each angle, and time :


where , , , and denote the number of snapshots in the 3D time series, the number of azimuth angles, and the number of points on the horizontal axes. Limb darkening laws in a given wavelength band are then derived by numerically integrating the averaged surface intensities over the band:


The result is a smooth function of . Interpolation of the limb darkening samples in using, e.g., cubic splines yields the best representation of the numerical results, but this method is more difficult to handle in the computation of model light curves. Sufficient accuracy can still be achieved with simple analytical fit formulae for stellar types that are relevant to planetary transit analyses, even though this approach may be more problematic for other cases (see, e.g., the results for supergiant stars in Chiavassa et al. 2009). The simple linear law,


is generally not a good approximation as it cannot capture the inherent nonlinearity of stellar limb darkening (see Fig. 2), which may lead to unphysical results in the light curve analysis (Southworth 2008). We use a widely applied formula proposed by Claret (2000) for plane-parallel calculations,


which is fitted to the intensity points with equal weighting. The linear combination of 4 power laws exhibits some degeneracy in the coefficients, which can impede a per-coefficient comparison of limb darkening laws predicted by different models. However, the law is sufficiently versatile to fit a large range of brightness distributions with very small residuals and very good flux conservation (see Table 3), which avoids the inaccuracies and ambiguities of using lower-order laws; see, e.g., the discussion in Heyrovský (2007).

3.2 Numerical resolution of the 3D calculations

The monochromatic surface brightness distribution on the stellar disk is resolved with non-equidistant points in the range , adopting the same set of angles that is used for the Kurucz grid (see the left panel of Fig. 2), and angles per point for all . We reduce the horizontal resolution of the model atmospheres by half to accelerate the computation. The resulting intensities are averaged over time sequences with snapshots that span periods of the fundamental pressure oscillation mode in each model. This is equivalent to  min of stellar time in the case of HD 209458 and  min for HD 189733.

We test the accuracy of our choice of numerical resolution by comparing cubic spline-interpolated theoretical limb darkening predictions in the spectral region between  Å and  Å, where limb darkening is strong and deviates significantly from a linear shape, using an arbitrary snapshot of the 3D model for HD 189733.

The numerical result for the computation with , and with our “standard” resolution , and is plotted as crosses in the left panel of Fig. 2. The solid line shows a least-squares fit using the nonlinear law (Eq. 5), which yields a good representation of the model limb darkening. The dotted line shows a linear law (Eq. 4), which clearly does not provide sufficient accuracy. The gray-shaded area above approximately indicates the inner region of the stellar disk that is never reached by the transiting planet, assuming an impact parameter and a planet-to-star radius ratio (see Table 5).

The right panel of Fig. 2 compares spline-interpolated computations with (solid line) and with the full horizontal resolution of the 3D models (, dotted line), keeping the other parameters constant. The relative deviation from the standard resolution is smaller than  % in either case, indicating that limb darkening is well-represented by our choices.

Larger sensitivity of the interpolated limb darkening law is observed when the number of angles is increased, refining the radial sampling of the surface brightness distribution. The left panel in Fig. 3 shows the relative deviation of interpolated computations with points (solid line) and (dotted line) from . The node distribution of the and cases follows Gauss-Legendre quadrature (plus disk center at ). The set reaches 1 % relative deviation near the limb at as the interpolated curve is visibly less constrained between the nodes (marked by grey vertical lines in the plot). The differences above are practically negligible, as are those between the and the cases (dotted line). The nonlinear law cannot represent such small-scale variations in the limb darkening curve; increasing leads to slightly larger deviation of the computed intensity samples from the best-fit curve (right panel of Fig. 3), while the curve itself changes weakly. Table 3 compares the best-fit coefficients of the nonlinear law (Eq. 5) and the reduced that results from the different choices of . The degeneracy of the four coefficients does not allow a direct comparison; they change by up to  % () between and , which is much larger than the actual point-to-point difference in the curves. The inability of the nonlinear law to take small-scale variations into account is reflected by an increasing for increasing .

Table 3 also compares the surface flux predicted by the -integral of the nonlinear law,


with the direct numerical integral of the 3D computation, which uses the trapezoid rule for and Gauss-Legendre quadrature for and (the disk center has zero weight in the latter two cases). We obtain excellent flux conservation for the nonlinear law with a relative deviation of less than for and . The trapezoid rule with leads to less accurate integration, as the sampling of the surface brightness distribution is less optimal compared to Gauss-Legendre quadrature. The integral of the fitted nonlinear law, however, deviates much less from the computations with and , as the fit appears less sensitive to the sampling of the surface brightness distribution, which emphasizes again that the number of angles in our standard resolution is sufficient.

Our full spectrum computations contain surface intensities for a set of wavelengths over a spectral range between  Å and with a constant sampling rate of . This statistical treatment of the stellar spectrum leads to inaccuracies in the spectral energy distribution (SED) when narrow wavelength ranges are investigated due to missing information on spectral line profiles between the opacity samples. Plez (2008) estimated the uncertainty of the solar radiative flux derived from a MARCS model for constant wavelength bin sizes of , where is the center wavelength in bin , by reducing the number of opacity samples per bin by factors of 3. Errors reach the level at  Å with  Å, they increase towards the UV and decrease towards longer wavelengths.

The wavelength sampling rate affects the strength of band-integrated limb darkening, which is weaker in strong lines than at continuum wavelengths (see the discussion in Appendix B). More frequent and stronger lines in the opacity sampling thus effectively lead to weaker limb darkening. We test the importance of wavelength resolution between  Å and  Å by reducing the number of sampling points by 2 and 4, resulting in sampling rates of and , respectively. The relative deviations of spline-interpolated limb darkening predictions from full wavelength resolution are compared in the left panel of Fig. 4, where solid lines show the two possible choices of resampling for and dotted lines show the four possible choices for . The deviation can reach up to  % near the limb at the lowest sampling rate. It is important to keep in mind that this result strongly depends on the width of the wavelength band for which limb darkening is calculated and on the wavelength region itself through the presence or absence of spectral lines.

We finally investigate the influence of time-dependent variation in the atmospheric stratification. Similar to real stars, 3D time-dependent hydrodynamical models exhibit pressure-mode oscillations that lead to periodic steepening of the mean - gradient and, as a consequence, to slight variations in the strength of limb darkening (see also Sect. 4). Surface intensities are averaged over a sample of snapshots in our standard setting in order to obtain a robust representation, typically covering periods of the fundamental oscillation mode. We compare the resulting limb darkening predictions with a time sequence that is twice as long (solid line in the right panel of Fig. 4) and find a relative deviation of up to  % near the limb, declining quickly towards the disk center. The difference between the two sets is slightly larger (dotted line).

3.3 Summary

Most uncertainties of numerical resolution affect the limb darkening predictions only very close to the limb beneath , a region where the assumption of plane-parallel geometry starts to yield generally unrealistic radiative intensities as it was discussed above. This affects the earliest ingress and latest egress phases during the transit of HD 209458 and HD 189733. Characteristic residuals in the model light curve that result from inaccurate predictions near the limb are not observed (see Sect. 5), although the available data only covers parts of the critical regions and may not be sufficiently accurate to clearly highlight discrepancies. There is still some influence of the calculations at small on the entire limb darkening law due to the finer sampling, giving it more relative weight in the least-squares fit. As HST light curve analyses based on our 3D model for HD 189733 were shown to yield excellent results that rival the quality of direct fits of limb darkening (see Table 2 in Sing et al. 2011), we trust that our calculations provide a sufficiently accurate and realistic description of the effect.

The assumptions and simplifications used in this work may fail in other cases: time-variation of the atmospheric temperature gradient may be a significant issue for very weak transit signals which occur in systems with smaller planet-to-star radius ratios, such as earth-sized planets with solar-type host stars. The effects of spherical geometry will certainly become more important for transiting planets that reach only the outermost parts of the stellar disk. We also stress that an investigation of limb darkening for increasingly narrow wavelength bands requires more detailed spectrum synthesis with better wavelength resolution, as the wavelength sampling effects will become increasingly important. This issue affects observations with future generations of telescopes that will allow narrowband transit spectroscopy of exoplanetary atmospheres.

4 Limb darkening predictions of the 3D models and 1D models

Figure 5: Surface intensities for HD 209458 at different angles, including only continuous absorption (left panel) as well as continuous and spectral line absorption (right panel) in the radiative transfer computation. Time-averaged 3D spectra (black solid lines) are compared to spectra derived from the mean 3D model (green solid lines) and the 1D MARCS model (red solid lines). Dotted black lines show spectra derived from 3D snapshots that exhibit intensity minima and maxima at  Å, bracketing the time variation of the 3D surface brightness due to oscillations in the model atmosphere. Note that spectra in the right panel were convolved with a Gaussian kernel with for clarity.
Figure 6: Same as Fig. 5, computed for HD 189733.
Figure 7: Left: Wavelength-integrated limb darkening curves for the 3D model of HD 209458 (black solid line), the mean 3D model (green solid line) and the 1D MARCS model (red solid line). Surface intensities between  Å and  Å were weighted with the sensitivity function of the HST STIS instrument and wavelength (see Appendix A). Right: Relative deviation between the limb darkening predictions of the 1D MARCS model with respect to the 3D model.
Figure 8: Same as Fig. 7 for the case of HD 189733.

Classical 1D models of the solar atmosphere are known to overestimate the effect of limb darkening compared to solar continuum observations of Neckel & Labs (1994); Sing et al. (2008) examine an ATLAS model and Asplund et al. (2009) a MARCS model. The latest generation of 3D hydrodynamical model atmospheres, however, exhibits very good agreement with the solar continuum (see Fig. 2 in Asplund et al. 2009) and other observational tests (e.g. Pereira et al. 2009a, b).

3D models take surface granulation explicitly into account rather than parameterizing convection with a simplified recipe, producing strong horizontal temperature fluctuations and time-dependent pressure oscillations that affect the temperature structure. The spatial and temporal mean temperature gradient thus deviates from a 1D model with the same stellar parameters. While this often accounts for the leading order effect, the temperature fluctuations with respect to the mean stratification can also have an impact on the emitted radiation.

Quantitative limb darkening predictions from numerical models rely on a complex nonlinear coupling of radiative transfer, plasma dynamics and microphysics. The main mechanisms can still be understood using basic results from 1D stellar atmosphere theory, which yields the wavelength-dependent approximate relation


for the steepness of a linear limb darkening law (cf. Eq. (4)). is the Planck function and is the temperature stratification of the model atmosphere as a function of monochromatic vertical optical depth ; see Appendix B for a derivation. In this linear model, the strength of limb darkening at a given wavelength scales directly with the temperature gradient at the optical surface . The wavelength-dependence of limb darkening stems from the relative temperature sensitivity of the Planck function, expressed through the term , which generally leads to an increase towards shorter wavelengths, as well as from opacity variations that change the monochromatic optical depth. The latter effect can play an important role in spectral features, such as strong absorption lines, and leads to weaker limb darkening if the temperature decreases monotonously outward; see Appendix B for a discussion. The simple model behind Eq. (7) assumes linearity of the Planck function with optical depth and thus cannot include any variation of the gradients found in real stellar atmospheres; its applicability is therefore limited to a more or less narrow region around the optical surface.

4.1 Hd 209458

We compute surface intensities in the spectral range between  Å and  Å for the 3D model of HD 209458, its spatially and temporally averaged 1D stratification (mean 3D model), and the 1D MARCS model which has a lower effective temperature (see Table 2). The left panel in Fig. 5 compares continuum intensities as a function of wavelength and for the different models.

The wavelength-dependence of limb darkening is clearly visible in the decreasing relative intensity towards shorter wavelengths for given , resulting in a stronger effect. The opacity-dependence of limb darkening manifests itself in the reverse Balmer jump at  Å: on the blue side, the additional H I bound-free opacity leads to weaker limb darkening compared to the red side where the atmosphere is more transparent.

The 3D model (black solid lines in Fig. 5, the spectrum of each individual 3D snapshot is computed individually and then averaged over time) predicts significantly weaker limb darkening than the 1D MARCS model (red solid lines) in the range shown in the plot. The deviations decrease towards the limb as the 3D and 1D temperature gradients in the higher atmosphere become more similar. Using a 1D MARCS model with the same effective temperature as the 3D model (not shown in the plot) only leads to a small deviation from the cooler 1D model. The result of the 3D-1D comparison is similar to the solar case (see Fig. 2 in Asplund et al. 2009), as expected from the similarity of stellar parameters.

The mean 3D model (green sold lines in Fig. 5) yields almost the same continuous spectrum as the full 3D model well down to the Balmer jump, indicating that the mean temperature stratification dominates the 3D effect in this part of the spectrum. At shorter wavelengths, temperature fluctuations play a greater role for continuum formation, leading to larger deviations between the 3D prediction and the mean 3D result. The time-dependent pressure oscillations and temperature fluctuations of the 3D model lead to slight variation in the normalized intensities; the black dotted lines indicate the spectra with minimum and maximum intensity at  Å.

Strong absorption features appear inverted in the right panel of Fig. 5 due to weaker limb darkening for larger opacity. Broadband limb darkening curves therefore integrate a variety of limb darkening strengths if spectral lines are present, depending on number and strength of the features. It is thus essential to use the same line data and chemical composition for a model comparison as it is performed in this paper. The deviation between 3D spectra and 1D spectra decreases in the presence of line absorption: a larger atmospheric height range influences the outgoing radiation field, with temperature gradients of the 3D model and 1D MARCS model becoming more similar higher up in the atmosphere. The slightly hotter temperatures of the 1D MARCS model beneath remove or even reverse the 3D-1D difference; the effect is visible in the features around  Å in the right panel of Fig. 5.

Wavelength-integrated limb darkening curves for HD 209458 are shown in the left panel of Fig. 7, the monochromatic intensity distributions between  Å and  Å were weighted with the sensitivity function of the STIS instrument and with wavelength (see Appendix A). The shape of the 3D curve is nonlinear across the stellar disk, with a slightly steeper drop near the limb. A comparison of the 3D and 1D MARCS results immediately reveals a clear distinction between the two models, despite the influence of spectral lines. The mean 3D model produces almost the same limb darkening curve as the full 3D model. The right panel of Fig. 7 shows the relative deviation between the 3D and 1D results, which decreases close to the limb as the temperature stratification of the 1D MARCS model becomes shallower compared to the 3D model in the higher atmosphere.

4.2 Hd 189733

The same computation is repeated for HD 189733; the left panel of Fig. 6 shows continuum intensities, the right panel continuum and line intensities. The 3D model (black solid lines) again predicts overall weaker limb darkening in the continuum, but the difference between the 3D and 1D MARCS (red solid lines) results is visibly smaller compared to HD 209458, owing to the smaller difference in the atmospheric temperature stratifications of the models deeper in the photosphere (see Fig. 1). The mean 3D model (green solid lines) represents the continuous spectrum well also below the Balmer jump, indicating that the mean temperature gradient of the 3D model has the dominant influence. The time variation of the 3D curves, indicated by black dotted lines in Fig. 6, is small.

The 1D MARCS model is again warmer than the 3D model in the higher atmosphere, but the difference is much more pronounced for HD 189733. The shallow 1D temperature gradient leads to significantly smaller deviations between the 3D and 1D continuum intensities at compared to the region around, e.g., . The consequences of the hotter 1D temperatures become even more apparent when spectral lines are included in the computation (right panel of Fig. 6). The cores of the strongest lines that form in the higher photosphere exhibit weaker limb darkening in 1D than in 3D, visible in the features beneath  Å.

Spectral line absorption is overall stronger in HD 189733 due to its cooler atmospheric temperatures. The broadband-integrated limb darkening curves therefore include a larger variety of limb darkening strengths and are more heavily influenced by temperature gradients at different atmospheric heights. The deviation between 3D and 1D limb darkening curves (black and red solid line in the left panel of Fig. 8) is significantly smaller compared to HD 209458; the differences in the atmospheric temperature stratification are partially hidden through the influence of spectral lines. The mean 3D model (green solid line) predicts practically the same limb darkening curve as the full 3D model, with the exception of the region closest to the limb. The right panel of Fig. 8 shows again the relative deviation between the 3D and 1D limb darkening curves. The hotter atmospheric temperatures of the 1D MARCS model lead to visibly weaker broadband limb darkening compared to the 3D model beneath .

5 Transit light curve fits

We test the theoretical limb darkening predictions of the 3D models and the 1D MARCS models using HST light curves for HD~209458b and HD~189733b. Observations from the Space Telescope Imaging Spectrograph (STIS) with the G430L grating were integrated over the wavelength range between  Å and  Å where limb darkening is strong; theoretical intensity spectra were weighted with the sensitivity function of the instrument (see Appendix A; limb darkening coefficients for other instruments and a range of standard broadband filters are listed in Table 6). Using large spectral bandwidth in the comparison maximizes the signal-to-noise ratio of the observations and reduces the uncertainties of the opacity sampling technique (Sect. 3.2). The differences between the atmospheric stratifications of 3D models and 1D models can be obscured by such broadband integration if spectral lines are present, as discussed in Sect. 4, but light curve fits in smaller wavelength bands yielded essentially the same results for the 3D-1D comparison. To complement the sparse transit coverage of our HST STIS data for HD 189733, we also include light curve fits of observations between  Å and  Å that were obtained with the HST Advanced Camera for Surveys (ACS) and the HRC G800L grating, although the limb darkening effect is weaker towards longer wavelengths.

Observed light curves are analyzed with the models of Mandel & Agol (2002) and using least-squares fits with the Levenberg-Marquardt algorithm. The ratio of orbital to stellar radius, the orbital inclination and the stellar density are well known for both systems and therefore fixed to their literature values. The central transit time needs to be fitted for HD 209458b and the ACS data for HD 189733b. The remaining free parameters are the planet-to-star radius ratio, stellar baseline flux and several polynomial coefficients that correct for the above mentioned systematic effects of the instrument (see Sing et al. 2011).

Limb darkening affects the entire transit of the planet; its generally growing strength towards shorter wavelengths leads to an increasingly round shape of the light curve (see Fig. 3 in Knutson et al. 2007). This causes some degeneracy between limb darkening and the planet-to-star radius ratio. The latter parameter thus varies slightly between a direct fit of limb darkening and using model predictions. In addition to that, the removal of instrument systematics from the observed stellar flux with a polynomial function is partly degenerate with the shape of the transit light curve, allowing for compensation of inaccurate limb darkening predictions in a simultaneous fit. We therefore first remove the systematics using the out-of-transit points in each visit and then fit the resulting light curves with limb darkening predictions from 3D and 1D models, keeping only the baseline flux and the planet-to-star radius ratio as free parameters in each case. While this split in the process does not necessarily produce the best fit to the data, it allows for a clear comparison of limb darkening predictions without any obliteration by other free parameters.

5.1 HD 209458b

HD 209458
3D 1D
Parameter Visit 1 Visit 2 Visit 1 Visit 2
Incl. [deg]
 [g cm]
DOFDegrees of freedom.
444 Fitted parameters have error bars. Parameters kept constant for different fits appear repeatedly to improve readability of the table.
Table 4: Parameters for the HST transit light curves of HD 209458b, shown in Fig. 9.
Figure 9: Transit light curve fits for two visits of HD 209458b (dots in the lower panels of each plot), integrated over the wavelength range between  Å and  Å with limb darkening coefficients derived from the 3D model (green lines) and the 1D MARCS model (red lines); see Table 4 for the fit parameters. The residuals of the 3D fits and the 1D fits are shown in the upper panels and center panels in each plot, the blue lines show the deviation between the 3D and 1D model light curves.

We use the HST data of Knutson et al. (2007) for our investigation, taken as part of GO 9447 (PI Charbonneau). Results from the HD 209458b STIS dataset have also been published in Ballester et al. (2007) and Sing et al. (2008). Notably, we re-reduced the HD 209458b dataset in the same manner as in Sing et al. (2011) for HD 189733b, producing a wavelength-integrated light curve which also has the detector position systematic trends taken into account. We find that this addition noticeably improves the quality. Transits were observed in two visits, the results of our individual fits are shown in Fig. 9. The two visits together cover almost the entire light curve. Combined with the large 3D-1D effect on the limb darkening curve (cf. Fig. 5 and Fig. 7), the HD 209458 system allows for a clear distinction between the stellar model atmospheres. The improvement of the 3D model is visible in the comparison of residuals of the 3D fit and the 1D fit (black dots in upper and center panel of Fig. 9), which is accompanied by a significant reduction of the parameter from 3D to 1D (see Table 4). The differences between the two light curves are visualized by the blue lines in Fig. 9, which closely follow the characteristic shape of the 1D fit residuals (see also Fig. 3 of Knutson et al. (2007). The 3D-1D effect is clearly strongest in the ingress and egress phases, but also leads to a noticeable deviation around the central transit. Our result strongly indicates that the 3D model provides a more realistic description of the atmospheric temperature structure of HD 209458, which agrees well with the case of the Sun (Asplund et al. 2009).

The shallower model light curve that results from overall weaker 3D limb darkening requires a slightly larger planet-to-star radius ratio to fit the observations. Averaged over the two visits in either case, the parameter increases by 0.2 %. Even though the overall magnitude of the 3D correction on is small, the effect can still be significant for an accurate characterization of the planetary atmosphere through transmission spectroscopy, in particular if the correction is wavelength-dependent. It would therefore be interesting to investigate possible consequences for the determination of the atmospheric structure and composition of HD 209458b, but such analyses lie beyond the scope of this paper.

We note that a comparison of predicted linear and quadratic limb darkening coefficients with empirical results as in Southworth (2008) and Claret (2009) leads to similarly strong disagreement between our 3D model and the observational data, despite the clear outcome of the transit light curve fits. We suspect that inaccuracies of the linear law and also, to some degree, of the quadratic law in fitting the actual limb darkening curve and the interaction between the different free parameters of a direct fit of limb darkening are responsible for this discrepancy. Southworth (2008) indeed warns that the two coefficients of the quadratic law are correlated and therefore cannot always be uniquely determined, which inhibits a clearer comparison.

5.2 HD 189733b

HD 189733
3D 1D
Incl. [deg]
 [g cm]
555 Fitted parameters have error bars. Parameters kept constant for different fits appear repeatedly to improve readability of the table.
Table 5: Parameters for the HST transit light curves of HD 189733b, shown in Fig. 10.
Figure 10: Left: Transit light curve fits for one HST STIS visit of HD 189733 (dots in the lower panel, data from Sing et al. 2011), integrated over the wavelength range between  Å and  Å with limb darkening coefficients derived from the 3D model (green line) and the 1D MARCS model (red line); see Table 5 for the fit parameters. The residuals of the 3D fit and the 1D fit are shown in the upper panel and center panel, the blue line shows the deviation between the 3D and 1D model light curves. Right: Transit light curve fits for one HST ACS visit of HD 189733 (data from Pont et al. 2007), integrated over the wavelength range between  Å and  Å.

We use transit observations of HD 189733b from program GO 11740 (PI FP), with the data reduction method as detailed in Sing et al. (2011). The HST STIS light curves are strongly distorted by star spots (see Fig. 3 in Sing et al. 2011). While the second visit is affected in its entirety, we can use the first visit by ignoring the peak near the transit center. The observation (black dots in the lower panel on the left side of Fig. 10) thus only includes a part of the egress phase which results in weaker constraints on limb darkening compared to the HD 209458 system. We also analyze the first visit of the HST ACS data of Pont et al. (2007), which is not distorted by star spots and covers parts of the ingress phase between  Å and  Å (right side of Fig. 10). The model light curves based on 3D predictions are represented by green lines in the figure, red lines show the 1D MARCS case.

Both models reach equivalent fit quality for the STIS data (see Table 5). The fit residuals are shown in the upper and center panels on the left side of Fig. 10. The deviation between the two model light curves is represented by the blue solid line in the center panel, exhibiting again the characteristic shape of a too strong 1D limb darkening prediction. The magnitude of the effect is smaller compared to HD 209458b due to the smaller 3D-1D difference and due to stronger obstruction through spectral lines (Fig. 6 and Fig. 8). This test is not decisive enough to clearly distinguish between the temperature structures of the 3D and 1D models, as both cases deliver equivalent results with respect to the comparatively sparse observational data.

Fits of the ACS data yield a smaller for 3D limb darkening (see Table 5), the residuals are shown on the right side of Fig. 10. The light curve at longer wavelengths has a flatter bottom, as limb darkening is weaker due to the lower temperature-sensitivity of the Planck function. The 3D-1D difference is smaller for the same reason (blue line).

Contrary to HD 209458b, the light curves with 3D limb darkening yield slightly smaller planet-to-star radius ratios; the differences reach again 0.2 % for both the STIS data and ACS data.

Sing et al. (2011) compared light curves based on direct fits of limb darkening with a linear law and our 3D model predictions with the nonlinear law for various narrower wavelength bands within the spectral region of the STIS G430L grating. They showed that 3D limb darkening yields model light curves with similar fit quality in terms of the achieved , even though one less free parameter is needed (see their Table 2 for details). We therefore conclude that the 3D model delivers realistic limb darkening predictions, but the 1D MARCS model reaches similar fit quality for the STIS data due to the low sensitivity of the test.

6 Conclusions

We compute theoretical surface intensity distributions for the two exoplanet host stars HD 209458 and HD 189733 using 3D time-dependent hydrodynamical and 1D hydrostatic MARCS model atmospheres. The stellar parameters are determined using Fe I and Fe II abundances based on HARPS spectra, leading to slightly lower effective temperatures of the 1D models due to their different temperature stratification.

A comparison between predicted continuum intensities as a function of wavelength and radial position on the stellar disk between , where is the cosine of the polar angle with respect to the vertical axis in the atmosphere, shows that the 3D models exhibit generally weaker limb darkening in the visual region than the 1D models, similar to an analysis of the solar atmosphere by Asplund et al. (2009). The overall shallower temperature structure of the temporal and spatial mean 3D models around the continuum optical surface in the disk center is to large parts responsible for this result: 3D models take convective motions in the atmosphere explicitly into account, rather than using a simplified prescription of convective heat flux that is inherent to 1D hydrostatic models. The 3D models exhibit slightly cooler temperatures higher up in the atmosphere compared to 1D models, which can result in a stronger 3D darkening effect close to the limb.

The presence of spectral lines partly obstructs a comparison of models if theoretical spectra are integrated over broad wavelength bands. These features form over a wide atmospheric height range and lead to a blend of limb darkening strengths that partly hides deviations in the temperature structure. It is therefore important to only compare theoretical spectra that were computed with the same set of opacities and chemical composition if the realism of stellar model atmospheres is investigated.

We analyze the results of model light curve fits to HST transit observations of the exoplanets HD 209458b and HD 189733b based on 3D and 1D limb darkening predictions. The 3D model produces a significantly better fit to the STIS light curve of HD 209458b and is capable of removing systematic residuals that appear when 1D limb darkening is used, which strongly indicates that the 3D model provides a more realistic description of the stellar temperature stratification. The improvement with 3D limb darkening is weaker for the ACS data of HD 189733b, and both models deliver equivalent fit quality for the STIS light curves. While the differences between 3D and 1D limb darkening predictions are generally smaller for this star, they are further reduced by spectral line absorption at shorter wavelengths or by the weaker limb darkening effect at longer wavelengths, respectively. A previous investigation of the same STIS data by Sing et al. (2011) showed that our 3D predictions rival the quality of a direct fit of limb darkening in different wavelength bands even though one free parameter less is used, which underlines the realism of our calculations.

Full intensity spectra for both stars will be made available online. However, it is important to emphasize that the opacity sampling technique used for our computations is not sufficiently accurate to represent spectral line absorption in narrow wavelength bands. Obtaining limb darkening coefficients for these cases requires a more detailed computation of stellar surface intensities with much better wavelength resolution, a detailed line list and the inclusion of Doppler shifts by the convective motions in the atmosphere. Such computations are available from the authors on request.

The weak transit signal of exoplanets with small planet-to-star radii, such as super-earth planets with solar-type host stars, will require better accuracy of the limb darkening predictions of stellar model atmospheres than what 1D models are able to provide. Already existing and future planet finding missions, such as Kepler or Plato, and the new James Webb Space Telescope (JWST), underline the importance of such efforts. A large grid of 3D hydrodynamical models is currently under construction, with surface gravities that cover the range from dwarf stars to red giant stars, effective temperatures from F-type to late K-type, and metallicities between solar and very metal-poor; see Collet et al. (2011) and Magic et al. (2012, in preparation) for details. The models will soon be used for an investigation of limb darkening across a wider range of stellar parameters.

The authors would like to thank I. Baraffe for inspiring the project, the Rechenzentrum Garching (RZG) for providing the large computational resources that were essential for this study and P. Baumann for providing co-added HARPS spectra. We would also like to thank the referee for helpful comments. W.H. acknowledges support by the European Research Council under the European Community’s 7th Framework Programme (FP7/2007- 2013 Grant Agreement no. 247060). This work is based on observations with the NASA/ESA Hubble Space Telescope, obtained at the STScI operated by AURA, Inc.


  • Asplund (2005) Asplund, M. 2005, ARA&A, 43, 481
  • Asplund et al. (2005) Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes III & F. N. Bash, 25
  • Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  • Asplund et al. (2000) Asplund, M., Nordlund, Å., Trampedach, R., & Stein, R. F. 2000, A&A, 359, 743
  • Auvergne et al. (2009) Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411
  • Ballester et al. (2007) Ballester, G. E., Sing, D. K., & Herbert, F. 2007, Nature, 445, 511
  • Blackwell et al. (1995) Blackwell, D. E., Lynas-Gray, A. E., & Smith, G. 1995, A&A, 296, 217
  • Bouchy et al. (2005) Bouchy, F., Udry, S., Mayor, M., et al. 2005, A&A, 444, L15
  • Brown et al. (2001) Brown, T. M., Charbonneau, D., Gilliland, R. L., Noyes, R. W., & Burrows, A. 2001, ApJ, 552, 699
  • Casagrande et al. (2011) Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138
  • Chiavassa et al. (2009) Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351
  • Claret (2000) Claret, A. 2000, A&A, 363, 1081
  • Claret (2009) Claret, A. 2009, A&A, 506, 1335
  • Claret & Hauschildt (2003) Claret, A. & Hauschildt, P. H. 2003, A&A, 412, 241
  • Collet et al. (2011) Collet, R., Magic, Z., & Asplund, M. 2011, Journal of Physics Conference Series, 328, 012003
  • Gonzalez et al. (2001) Gonzalez, G., Laws, C., Tyagi, S., & Reddy, B. E. 2001, AJ, 121, 432
  • Gustafsson et al. (1975) Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407
  • Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951
  • Hauschildt et al. (1999) Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377
  • Hayek et al. (2010) Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49
  • Hayek et al. (2011) Hayek, W., Asplund, M., Collet, R., & Nordlund, Å. 2011, A&A, 529, A158
  • Heiter & Luck (2003) Heiter, U. & Luck, R. E. 2003, AJ, 126, 2015
  • Heyrovský (2007) Heyrovský, D. 2007, ApJ, 656, 483
  • Holweger et al. (1995) Holweger, H., Kock, M., & Bard, A. 1995, A&A, 296, 233
  • Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Noyes, R. W., Brown, T. M., & Gilliland, R. L. 2007, ApJ, 655, 564
  • Kurucz (1996) Kurucz, R. L. 1996, in IAU Symposium, Vol. 176, Stellar Surface Structure, ed. K. G. Strassmeier & J. L. Linsky, 523
  • Mandel & Agol (2002) Mandel, K. & Agol, E. 2002, ApJ, 580, L171
  • Mashonkina & Gehren (2001) Mashonkina, L. & Gehren, T. 2001, A&A, 376, 232
  • Mazeh et al. (2000) Mazeh, T., Naef, D., Torres, G., et al. 2000, ApJ, 532, L55
  • Meléndez & Barbuy (2009) Meléndez, J. & Barbuy, B. 2009, A&A, 497, 611
  • Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres, 2nd edition, ed. J. Hevelius (San Francisco: W. H. Freeman and Co.)
  • Neckel & Labs (1994) Neckel, H. & Labs, D. 1994, Sol. Phys., 153, 91
  • Nordlund (1982) Nordlund, A. 1982, A&A, 107, 1
  • Nordlund & Galsgaard (1995) Nordlund, Å. & Galsgaard, K. 1995, A 3D MHD Code for Parallel Computers, Tech. rep., Astronomical Observatory, Copenhagen University
  • Nordlund et al. (2009) Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Living Reviews in Solar Physics, 6, 2
  • Pereira et al. (2009a) Pereira, T. M. D., Asplund, M., & Kiselman, D. 2009a, A&A, 508, 1403
  • Pereira et al. (2009b) Pereira, T. M. D., Kiselman, D., & Asplund, M. 2009b, A&A, 507, 417
  • Plez (2008) Plez, B. 2008, Physica Scripta Volume T, 133, 014003
  • Pont et al. (2007) Pont, F., Gilliland, R. L., Moutou, C., et al. 2007, A&A, 476, 1347
  • Pont et al. (2008) Pont, F., Knutson, H., Gilliland, R. L., Moutou, C., & Charbonneau, D. 2008, MNRAS, 385, 109
  • Ramírez & Meléndez (2004) Ramírez, I. & Meléndez, J. 2004, ApJ, 609, 417
  • Sadakane et al. (2002) Sadakane, K., Ohkubo, M., Takeda, Y., et al. 2002, PASJ, 54, 911
  • Santos et al. (2004) Santos, N. C., Israelian, G., & Mayor, M. 2004, A&A, 415, 1153
  • Sing (2010) Sing, D. K. 2010, A&A, 510, A21
  • Sing et al. (2011) Sing, D. K., Pont, F., Aigrain, S., et al. 2011, MNRAS, 1159
  • Sing et al. (2008) Sing, D. K., Vidal-Madjar, A., Désert, J., Lecavelier des Etangs, A., & Ballester, G. 2008, ApJ, 686, 658
  • Skartlien (2000) Skartlien, R. 2000, ApJ, 536, 465
  • Sousa et al. (2008) Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373
  • Southworth (2008) Southworth, J. 2008, MNRAS, 386, 1644
  • Southworth (2009) Southworth, J. 2009, MNRAS, 394, 272
  • Southworth (2010) Southworth, J. 2010, MNRAS, 408, 1689
  • Valenti & Fischer (2005) Valenti, J. A. & Fischer, D. A. 2005, ApJS, 159, 141
  • van Belle & von Braun (2009) van Belle, G. T. & von Braun, K. 2009, ApJ, 694, 1085

Appendix A Limb darkening coefficients for standard broadband filters and different instruments

HD 209458 HD 189733
Johnson U  Å
Johnson B  Å
Johnson V  Å
Johnson I  Å
Strömgren u  Å
Strömgren v  Å
Strömgren b  Å
Strömgren y  Å
SDSS u’  Å
SDSS g’  Å
SDSS r’  Å
SDSS i’  Å
SDSS z’  Å
Bessell J m
Bessell H m
Bessell K m
MK J m
MK H m
MK K’ m
MK Ks m
MK K m
MK L’ m
MK M’ m
STIS G430LTheoretical intensity spectra were integrated over the entire accessible wavelength range of the grating, leading to small deviations from the coefficients in Table 4 and Table 5.  Å
ACS HRC G800LTheoretical intensity spectra were integrated over the entire accessible wavelength range of the grating, leading to small deviations from the coefficients in Table 4 and Table 5.  Å
Kepler  Å
Corot  Å
Spitzer IRAC 1 m
Spitzer IRAC 2 m
Spitzer IRAC 3 m
Spitzer IRAC 4 m
666All computations are based on the 3D models of HD 209458 and HD 189733.
Table 6: Theoretical limb darkening coefficients for a range of standard broadband filters and instruments.

Limb darkening coefficients for the Claret (2000) law were computed for a range of standard broadband filters and instruments, based on the 3D models of HD 209458 and HD 189733. The coefficients are derived from least-squares fits to the wavelength-integrated limb darkening law,


where is the filter or instrument response curve, respectively, and the factor converts the incoming intensity into a photon flux; see also Sect. 3. The effective wavelength of each filter was calculated using the integral


Table 6 gives a summary of the results. Standard filter curves were taken from the SYNPHOT package (SYNPHOT is a product of the Space Telescope Science Institute, which is operated by AURA for NASA). The New Mauna Kea (MK) filter set was downloaded from The STIS and ACS sensitivity curves were downloaded from the STScI calibration database. The Kepler sensitivity curve was downloaded from the Kepler website ( The Corot sensitivity curve was extracted from Auvergne et al. (2009). The Spitzer sensitivity curves were downloaded from

Appendix B Dependence of limb darkening on the atmospheric temperature gradient

The strength of limb darkening is closely related to the vertical atmospheric temperature gradient near the optical surface and depends on wavelength, which will be established in the following using basic stellar atmosphere theory (see, e.g., Mihalas 1978). In the Eddington-Barbier approximation, the depth-dependence of the Planck function is simplified by the linear expression


where is the vertical monochromatic optical depth with outside the star, and marks the optical surface at wavelength . and are constants, which can be expressed as




where is the Planck function at the optical surface in the stellar disk center. In plane-parallel geometry, the source function (Eq. (10)) produces the surface radiation field


which has linear anisotropy and consequently leads to a linear limb darkening relation. The second equality means that the surface intensity observed at angle is emitted from vertical optical depth in the atmosphere. Inserting the constants and , one obtains


which is equivalent to a limb darkening law


Using the chain rule, one can express the strength of of the limb darkening effect through


where is the vertical atmospheric temperature gradient at the monochromatic optical surface ; the logarithmic optical depth is commonly used as a depth variable and facilitates comparison with Fig. 1. A steeper gradient at the optical surface will therefore produce a stronger limb darkening effect. Likewise, if the opacity (and thus the atmospheric height of the optical surface) varies only weakly with wavelength, the increasing relative temperature sensitivity of the Planck function, expressed through the factor , increases the limb darkening effect monotonously towards shorter wavelengths.

The dependence of limb darkening for a given atmospheric stratification on monochromatic opacity, which can vary strongly between a spectral line core and the continuum, can be understood with similar arguments. Assuming a second opacity at nearly the same wavelength that scales linearly with respect to the first one, independent of depth,


leads to the relation


between the two optical depths. Using Eq. (16) and Eq. (18), the limb darkening strength is given by


For , the optical surface moves to larger atmospheric height at on the original scale where thermal emission is smaller if gas temperature decreases monotonously outward. At the same time, the source function gradient decreases by a factor of with respect to the original gradient . Using Eq. (10) through Eq. (12) and Eq. (16), the source function is expressed in terms of :


Inserting this result into Eq. (19), one obtains the simple relation


It is easy to verify that for and . only holds in the limits of an isothermal atmosphere () or with maximal limb darkening and zero emission at the limb (). Between these extremes, limb darkening is always weaker in line cores than in the continuum, if temperature decreases monotonously with increasing atmospheric height.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description