Issue 
A&A
Volume 539, March 2012



Article Number  A102  
Number of page(s)  15  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201117868  
Published online  01 March 2012 
Limb darkening laws for two exoplanet host stars derived from 3D stellar model atmospheres
Comparison with 1D models and HST light curve observations^{⋆}
^{1}
Astrophysics Group, School of Physics, University of Exeter,
Stocker Road,
Exeter
EX4 4QL,
UK
email: hayek@astro.ex.ac.uk; sing@astro.ex.ac.uk; fpont@astro.ex.ac.uk
^{2}
Research School of Astronomy & Astrophysics, Cotter Road,
Weston Creek
2611,
Australia
email: martin@mso.anu.edu.au
^{3}
Max Planck Institut für Astrophysik, KarlSchwarzschildStr. 1, 85741
Garching,
Germany
Received: 10 August 2011
Accepted: 14 January 2012
We compare limb darkening laws derived from 3D hydrodynamical model atmospheres and 1D hydrostatic MARCS models for the host stars of two wellstudied transiting exoplanet systems, the latetype 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 sampling^{⋆}. We test our theoretical predictions using leastsquares fits of model light curves to wavelengthintegrated 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 2900 Å and 5700 Å 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 2900 Å and 5700 Å, 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 5350 Å and 10 500 Å 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 latetype 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 2900–5700 Å region are not sufficient to distinguish more clearly between the 3D and 1D limb darkening predictions.
Key words: stars: atmospheres / stars: individual: HD 209458 / stars: individual: HD 189733 / planets and satellites: individual: HD 209458b / planets and satellites: individual: HD 189733b
Full theoretical spectra for both stars are available in electronic form at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/539/A102, as well as at www.astro.ex.ac.uk/people/sing.
© ESO, 2012
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 flatbottom 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 signaltonoise 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 loworder 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 wavelengthdependent degeneracy of limb darkening with transit depth and thus with the important planettostar 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 2900 Å and 5700 Å covered by the HST STIS G430L grating, where limb darkening is more sensitive to the photospheric temperature gradient, using leastsquares fits of a nonlinear law to the numerical calculations. We additionally investigate the region between 5350 Å and 10 500 Å 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
Fig. 1
Temperature distributions (gray shades) on surfaces with constant optical depth at 5000 Å 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 log τ_{5000} = 0 and log τ_{5000} = −2, from which the dominant contribution to the continuum surface brightness between the disk center at μ = 1.0 and the limb near μ = 0.01 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 latetype 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 timedependent radiationhydrodynamical 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 timeindependent radiative transfer are solved simultaneously in a planeparallel box with a resolution of 240 × 240 × 240 grid points, covering ~10 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, intensitymean 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 12 Mm × 12 Mm (HD 209458) and 4 Mm × 4 Mm (HD 189733) on equidistant axes with periodic boundaries. On the vertical axis, which has finer resolution around the continuum optical surface of the models, 5.9 Mm (HD 209458) and 2.2 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 10^{6} ≲ τ_{5000} ≲ 10^{7} at 5000 Å 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 signaltonoise HARPS spectra^{1}. The observed equivalent widths, derived from leastsquares 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 χ_{ex}, where low excitation lines exhibit the strongest sensitivity to temperature due to the excitationionization 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.
Abundances of neutral and ionized iron lines, derived from 3D hydrodynamical and 1D hydrostatic model atmospheres.
Literature values compared to our adopted effective temperature, surface gravity, and metallicity.
A large number of parameter determinations for the Gdwarf 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 ±150 K) produce some scatter with values ranging between approximately 5990 K and 6120 K. We obtain a flat distribution of Fe i abundances with increasing excitation level for a 3D model with timeaveraged ⟨ T_{eff} ⟩ = 6095 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 T_{eff} = 6000 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 5000 Å (green line) with the corresponding Tτ 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 log τ_{5000} = 0. 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 log g = 4.3 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 ±0.1 in log g. Southworth (2009) found log g = 4.368 from a light curve analysis, while Casagrande et al. (2011) derived log g = 4.33 from a Hipparcos parallax. Asplund et al. (2005) recommend a solar iron abundance of log ϵ(Fe ii)_{⊙} = 7.45, which is sufficiently close to our 3D and 1D Fe ii results for HD 209458, given the uncertainties in log g 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 ξ = 1.4 km s^{1} for the 1D case from the distribution of linetoline abundances as a function of equivalent width; a MARCS model atmosphere with slightly lower microturbulence (ξ = 1.0 km s^{1}) 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 Kdwarf HD 189733 exist in the literature (see lower panel in Table 2), which deviate by up to about 110 K and are therefore in reasonable agreement considering the systematic effects of the different methods. We determine ⟨ T_{eff} ⟩ = 5050 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 T_{eff} = 5000 K.
We find good agreement for our 3D ionization equilibrium (see lower panel of Table 1) for the surface gravity log g = 4.53 of Bouchy et al. (2005), which is also compatible with the light curve fit of Southworth (2010) and the parallaxbased result of Casagrande et al. (2011). Similar agreement is found for the 1D model with its effective temperature of 5000 K. The accuracy of our log g 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 ξ = 0.9 km s^{1}; we therefore choose again a grid model with ξ = 1.0 km s^{1} for the 1D calculations.
Fig. 2
Left: theoretical average surface brightness distribution as a function of projection factor μ integrated over the spectrum between 3400 Å and 3500 Å, computed using our “standard” numerical resolution for an arbitrary snapshot of the 3D model of HD 189733 (crosses), and leastsquares fits of the Claret (2000) nonlinear law (solid line) and a linear law (dotted line). The grayshaded 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 N_{φ} = 8 direction angles (solid line) and with the full horizontal resolution of N_{x} = N_{y} = 240 grid points (dotted line) from the interpolated “standard” setting. 
3. LTE spectrum formation and limb darkening laws
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 Trampedach (priv. comm.), as well as sampled line opacities from 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 ξ = 1.0 km s^{1} to approximate its contribution to spectral line broadening. The inaccuracies of this method should be small as we only investigate broadbandintegrated 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 planeparallel stratification inherently limits the realism of surface brightness calculations near the limb, where a sharp dropoff 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 μ_{d} 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 μ_{d} can be obtained by assuming a transparent atmosphere of one pressure scale height, which leads to the expression (1)where R is the stellar radius and H is the pressure scale height. Adopting R_{HD 209458} ≈ 1.23 R_{⊙} and R_{HD 189733} ≈ 0.77 R_{⊙} from van Belle & von Braun (2009), as well as choosing H_{HD 209458} ≈ 200 km and H_{HD 189733} ≈ 100 km at the surface as predicted by our 3D simulations, we obtain a dropoff at μ_{d} ≈ 0.02 for both stars. This result is in good agreement with Claret & Hauschildt (2003), who find μ_{d} ≈ 0.05 for a star with T_{eff} = 5000 K and log g = 4.0. 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 planeparallel stratification should therefore predict surface brightness distributions with sufficient realism for light curve analyses, except for some extreme cases of grazing transits.
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 N_{μ}.
3.1. Analytical limb darkening laws
SCATE yields monochromatic surface intensities I_{λ}(x,y,μ,φ,t) at wavelength λ for a range of discrete positions on the stellar disk, assuming planeparallel geometry as the underlying 3D model atmospheres, with the projection factor μ = cosθ for angle θ in radial direction off the disk center. The resulting specific intensities are averaged over direction angle φ, horizontal dimensions x and y at each angle, and time t: (2)where N_{t}, N_{φ}, N_{x}, and N_{y} 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 I_{k}(μ)/I_{k}(1) in a given wavelength band k are then derived by numerically integrating the averaged surface intensities ⟨I_{λ}⟩(μ) over the band: (3)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, (4)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 planeparallel calculations, (5)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 percoefficient 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 lowerorder 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 N_{μ} = 17 nonequidistant points in the range 0.01 ≤ μ ≤ 1.0, adopting the same set of angles that is used for the Kurucz grid (see the left panel of Fig. 2), and N_{φ} = 4 angles per point for all μ < 1.0. We reduce the horizontal resolution N_{x} = N_{y} = 240 of the model atmospheres by half to accelerate the computation. The resulting intensities are averaged over time sequences with N_{t} = 10 snapshots that span ≈6 periods of the fundamental pressure oscillation mode in each model. This is equivalent to ≈70 min of stellar time in the case of HD 209458 and ≈30 min for HD 189733.
We test the accuracy of our choice of numerical resolution by comparing cubic splineinterpolated theoretical limb darkening predictions in the spectral region between 3400 Å and 3500 Å, 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 N_{t} = 1, and with our “standard” resolution N_{x} = N_{y} = 120, N_{φ} = 4 and N_{μ} = 17 is plotted as crosses in the left panel of Fig. 2. The solid line shows a leastsquares 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 grayshaded area above μ = 0.86 approximately indicates the inner region of the stellar disk that is never reached by the transiting planet, assuming an impact parameter b = 0.663 and a planettostar radius ratio R_{Planet}/R_{Star} = 0.158 (see Table 5).
The right panel of Fig. 2 compares splineinterpolated computations with N_{φ} = 8 (solid line) and with the full horizontal resolution of the 3D models (N_{x} = N_{y} = 240, dotted line), keeping the other parameters constant. The relative deviation from the standard resolution is smaller than 0.2% in either case, indicating that limb darkening is wellrepresented 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 N_{μ} = 17 points (solid line) and N_{μ} = 51 (dotted line) from N_{μ} = 101. The node distribution of the N_{μ} = 51 and N_{μ} = 101 cases follows GaussLegendre quadrature (plus disk center at μ = 1.0). The N_{μ} = 17 set reaches 1% relative deviation near the limb at μ = 0.01 as the interpolated curve is visibly less constrained between the nodes (marked by grey vertical lines in the plot). The differences above μ ≳ 0.1 are practically negligible, as are those between the N_{μ} = 51 and the N_{μ} = 101 cases (dotted line). The nonlinear law cannot represent such smallscale variations in the limb darkening curve; increasing N_{μ} leads to slightly larger deviation of the computed intensity samples from the bestfit curve (right panel of Fig. 3), while the curve itself changes weakly. Table 3 compares the bestfit coefficients of the nonlinear law (Eq. (5)) and the reduced χ^{2} that results from the different choices of N_{μ}. The degeneracy of the four coefficients does not allow a direct comparison; they change by up to ≈30% (c_{4}) between N_{μ} = 17 and N_{μ} = 101, which is much larger than the actual pointtopoint difference in the curves. The inability of the nonlinear law to take smallscale variations into account is reflected by an increasing for increasing N_{μ}.
Fig. 3
Left: relative deviation of the splineinterpolated limb darkening predictions for N_{μ} = 17 (solid line) and N_{μ} = 51 (dotted line) from the case with N_{μ} = 101 as a function of μ. Grey vertical lines mark the position of the 17 μ angles used in the “standard” resolution. Right: relative deviation of the leastsquaresfitted limb darkening laws from the computed data points for N_{μ} = 17 (crosses), N_{μ} = 51 (diamonds) and N_{μ} = 101 (triangles) as a function of μ. All calculations are based on the 3D model of HD 189733. 
Table 3 also compares the surface flux predicted by the μintegral of the nonlinear law, (6)with the direct numerical integral F_{num} of the 3D computation, which uses the trapezoid rule for N_{μ} = 17 and GaussLegendre quadrature for N_{μ} = 51 and N_{μ} = 101 (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 10^{4} for N_{μ} = 51 and N_{μ} = 101. The trapezoid rule with N_{μ} = 17 leads to less accurate integration, as the sampling of the surface brightness distribution is less optimal compared to GaussLegendre quadrature. The integral of the fitted nonlinear law, however, deviates much less from the computations with N_{μ} = 51 and N_{μ} = 101, 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 ≈108 000 wavelengths over a spectral range between 910 Å and 20 μm with a constant sampling rate of λ/Δλ = 20 000. 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 Δλ_{i} = λ_{i}/200, where λ_{i} is the center wavelength in bin i, by reducing the number of opacity samples per bin by factors of 3. Errors reach the 5%/ level at λ = 4000 Å with Δλ = 20 Å, they increase towards the UV and decrease towards longer wavelengths.
The wavelength sampling rate affects the strength of bandintegrated 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 3400 Å and 3500 Å by reducing the number of sampling points by 2 and 4, resulting in sampling rates of R = 10 000 and R = 5000, respectively. The relative deviations of splineinterpolated 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 R = 10 000 and dotted lines show the four possible choices for R = 5000. The deviation can reach up to 4.5% 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.
Fig. 4
Left: relative deviation of the splineinterpolated limb darkening predictions computed for a single snapshot of the 3D model for HD 189733 with reduced wavelength sampling of R = 10 000 (solid lines) and R = 5000 (dotted lines) from the full sampling with R = 20 000. Each set of resampled wavelengths is shown in either case, resulting in two or four choices, respectively. Right: relative deviation of the splineinterpolated limb darkening predictions computed for the first 10 snapshots of a time sequence with N_{t} = 20 (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. 
We finally investigate the influence of timedependent variation in the atmospheric stratification. Similar to real stars, 3D timedependent hydrodynamical models exhibit pressuremode oscillations that lead to periodic steepening of the mean Tτ 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 N_{t} = 10 snapshots in our standard setting in order to obtain a robust representation, typically covering ~6 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 1.5% 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 μ ≲ 0.05, a region where the assumption of planeparallel 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 leastsquares 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: timevariation of the atmospheric temperature gradient may be a significant issue for very weak transit signals which occur in systems with smaller planettostar radius ratios, such as earthsized planets with solartype 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
Fig. 5
Surface intensities I_{λ}(μ)/I_{λ}(1) 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. Timeaveraged 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 5700 Å, 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 λ/Δλ = 250 for clarity. 
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 timedependent 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 wavelengthdependent approximate relation (7)for the steepness u_{λ} of a linear limb darkening law (cf. Eq. (4)). B_{λ} is the Planck function and T(τ_{λ}) 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 τ_{λ} = 1. The wavelengthdependence of limb darkening stems from the relative temperature sensitivity of the Planck function, expressed through the term (dB_{λ}/dT)/B_{λ}, 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 2900 Å and 5700 Å 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 wavelengthdependence of limb darkening is clearly visible in the decreasing relative intensity towards shorter wavelengths for given μ, resulting in a stronger effect. The opacitydependence of limb darkening manifests itself in the reverse Balmer jump at 3646 Å: on the blue side, the additional H i boundfree 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 3D1D 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 solid 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 timedependent 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 5700 Å.
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 log τ_{5000} ≲ −2 remove or even reverse the 3D1D difference; the effect is visible in the features around 3800 Å in the right panel of Fig. 5.
Wavelengthintegrated limb darkening curves for HD 209458 are shown in the left panel of Fig. 7, the monochromatic intensity distributions between 2900 Å and 5700 Å 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.
Fig. 7
Left: wavelengthintegrated 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 2900 Å and 5700 Å 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. 
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 μ = 0.1 compared to the region around, e.g., μ = 0.5. 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 4000 Å.
Spectral line absorption is overall stronger in HD 189733 due to its cooler atmospheric temperatures. The broadbandintegrated 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 μ ≲ 0.1.
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 2900 Å and 5700 Å 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 A.1). Using large spectral bandwidth in the comparison maximizes the signaltonoise 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 3D1D comparison. To complement the sparse transit coverage of our HST STIS data for HD 189733, we also include light curve fits of observations between 5350 Å and 10 500 Å 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 leastsquares fits with the LevenbergMarquardt 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 planettostar radius ratio, stellar baseline flux and several polynomial coefficients that correct for the 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 planettostar 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 outoftransit 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 planettostar 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
Fig. 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 2900 Å and 5700 Å 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 rereduced the HD 209458b dataset in the same manner as in Sing et al. (2011) for HD 189733b, producing a wavelengthintegrated 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 3D1D effect on the limb darkening curve (cf. Figs. 5 and 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 χ^{2} 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 3D1D 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 planettostar 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 R_{Planet}/R_{Star} 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 wavelengthdependent. 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
Fig. 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 2900 Å and 5700 Å 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 5350 Å and 10 500 Å. 
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 5350 Å and 10 500 Å (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 3D1D difference and due to stronger obstruction through spectral lines (Figs. 6 and 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 χ^{2} 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 temperaturesensitivity of the Planck function. The 3D1D difference is smaller for the same reason (blue line).
Contrary to HD 209458b, the light curves with 3D limb darkening yield slightly smaller planettostar 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 χ^{2}, 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 timedependent 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 0.1 ≤ μ ≤ 0.9, 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 planettostar radii, such as superearth planets with solartype 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 Ftype to late Ktype, and metallicities between solar and very metalpoor; see Collet et al. (2011) and Magic et al. (in prep.) for details. The models will soon be used for an investigation of limb darkening across a wider range of stellar parameters.
Acknowledgments
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 coadded 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/20072013 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.
References
 Asplund, M. 2005, ARA&A, 43, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Nordlund, Å., Trampedach, R., & Stein, R. F. 2000, A&A, 359, 743 [NASA ADS] [Google Scholar]
 Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes III, & F. N. Bash, ASP Conf. Ser., 336, 25 [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ballester, G. E., Sing, D. K., & Herbert, F. 2007, Nature, 445, 511 [Google Scholar]
 Blackwell, D. E., LynasGray, A. E., & Smith, G. 1995, A&A, 296, 217 [NASA ADS] [Google Scholar]
 Bouchy, F., Udry, S., Mayor, M., et al. 2005, A&A, 444, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, T. M., Charbonneau, D., Gilliland, R. L., Noyes, R. W., & Burrows, A. 2001, ApJ, 552, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Claret, A. 2000, A&A, 363, 1081 [NASA ADS] [Google Scholar]
 Claret, A. 2009, A&A, 506, 1335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., & Hauschildt, P. H. 2003, A&A, 412, 241 [Google Scholar]
 Collet, R., Magic, Z., & Asplund, M. 2011, J. Phys. Conf. Ser., 328, 012003 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, G., Laws, C., Tyagi, S., & Reddy, B. E. 2001, AJ, 121, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 [NASA ADS] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W., Asplund, M., Collet, R., & Nordlund, Å. 2011, A&A, 529, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heiter, U., & Luck, R. E. 2003, AJ, 126, 2015 [NASA ADS] [CrossRef] [Google Scholar]
 Heyrovský, D. 2007, ApJ, 656, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Holweger, H., Kock, M., & Bard, A. 1995, A&A, 296, 233 [NASA ADS] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Noyes, R. W., Brown, T. M., & Gilliland, R. L. 2007, ApJ, 655, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 1996, in Stellar Surface Structure, ed. K. G. Strassmeier, & J. L. Linsky, IAU Symp., 176, 523 [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Mashonkina, L., & Gehren, T. 2001, A&A, 376, 232 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mazeh, T., Naef, D., Torres, G., et al. 2000, ApJ, 532, L55 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Meléndez, J., & Barbuy, B. 2009, A&A, 497, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mihalas, D. 1978, Stellar atmospheres, 2nd edition, ed. J. Hevelius (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Neckel, H., & Labs, D. 1994, Sol. Phys., 153, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A. 1982, A&A, 107, 1 [NASA ADS] [Google Scholar]
 Nordlund, Å., & Galsgaard, K. 1995, A 3D MHD Code for Parallel Computers, Tech. rep., Astronomical Observatory, Copenhagen University [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 Pereira, T. M. D., Asplund, M., & Kiselman, D. 2009a, A&A, 508, 1403 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pereira, T. M. D., Kiselman, D., & Asplund, M. 2009b, A&A, 507, 417 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plez, B. 2008, Phys. Scr. 133, 014003 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Gilliland, R. L., Moutou, C., et al. 2007, A&A, 476, 1347 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Pont, F., Knutson, H., Gilliland, R. L., Moutou, C., & Charbonneau, D. 2008, MNRAS, 385, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Ramírez, I., & Meléndez, J. 2004, ApJ, 609, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Sadakane, K., Ohkubo, M., Takeda, Y., et al. 2002, PASJ, 54, 911 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, N. C., Israelian, G., & Mayor, M. 2004, A&A, 415, 1153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sing, D. K. 2010, A&A, 510, A21 [Google Scholar]
 Sing, D. K., VidalMadjar, A., Désert, J., Lecavelier des Etangs, A., & Ballester, G. 2008, ApJ, 686, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Sing, D. K., Pont, F., Aigrain, S., et al. 2011, MNRAS, 1159 [Google Scholar]
 Skartlien, R. 2000, ApJ, 536, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Southworth, J. 2008, MNRAS, 386, 1644 [Google Scholar]
 Southworth, J. 2009, MNRAS, 394, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Southworth, J. 2010, MNRAS, 408, 1689 [NASA ADS] [CrossRef] [Google Scholar]
 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 van Belle, G. T., & von Braun, K. 2009, ApJ, 694, 1085 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Limb darkening coefficients for standard broadband filters and different instruments
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 leastsquares fits to the wavelengthintegrated limb darkening law, (A.1)where S_{λ} is the filter or instrument response curve, respectively, and the factor λ converts the incoming intensity I_{λ} into a photon flux; see also Sect. 3. The effective wavelength of each filter was calculated using the integral (A.2)Table A.1 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 irtfweb.ifa.hawaii.edu/~nsfcam/filters.html. The STIS and ACS sensitivity curves were downloaded from the STScI calibration database. The Kepler sensitivity curve was downloaded from the Kepler website (keplergo.arc.nasa.gov/kepler_response_hires1.txt). The CoRoT sensitivity curve was extracted from Auvergne et al. (2009). The Spitzer sensitivity curves were downloaded from irsa.ipac.caltech.edu/data/SPITZER/docs/irac/calibrationfiles/spectralresponse/.
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 EddingtonBarbier approximation, the depthdependence of the Planck function is simplified by the linear expression (B.1)where τ_{λ} is the vertical monochromatic optical depth with τ_{λ} = 0 outside the star, and τ_{λ} = 1 marks the optical surface at wavelength λ. a_{λ} and b_{λ} are constants, which can be expressed as (B.2)and (B.3)where B_{λ}(τ_{λ} = 1) is the Planck function at the optical surface in the stellar disk center. In planeparallel geometry, the source function (Eq. (B.1)) produces the surface radiation field (B.4)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 a_{λ} and b_{λ}, one obtains (B.5)which is equivalent to a limb darkening law (B.6)Using the chain rule, one can express the strength u_{λ} of of the limb darkening effect through (B.7)where dT/dlog τ_{λ} is the vertical atmospheric temperature gradient at the monochromatic optical surface τ_{λ} = 1; 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 (dB_{λ}/dT)/B_{λ}, 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, (B.8)leads to the relation (B.9)between the two optical depths. Using Eqs. (B.7) and (B.9), the limb darkening strength u′ is given by (B.10)For k > 1, the optical surface moves to larger atmospheric height at τ_{λ} = 1/k on the original scale where thermal emission B_{λ} is smaller if gas temperature decreases monotonously outward. At the same time, the source function gradient decreases by a factor of k with respect to the original gradient dB_{λ}/dτ_{λ}. Using Eq. (B.1), (B.2), (B.3) and (B.7), the source function is expressed in terms of u_{λ}: (B.11)Inserting this result into Eq. (B.10), one obtains the simple relation (B.12)It is easy to verify that for 0 ≤ u_{λ} ≤ 1 and k > 1. only holds in the limits of an isothermal atmosphere (u_{λ} = 0) or with maximal limb darkening and zero emission at the limb (u_{λ} = 1). Between these extremes, limb darkening is always weaker in line cores than in the continuum, if temperature decreases monotonously with increasing atmospheric height.
All Tables
Abundances of neutral and ionized iron lines, derived from 3D hydrodynamical and 1D hydrostatic model atmospheres.
Literature values compared to our adopted effective temperature, surface gravity, and metallicity.
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 N_{μ}.
Theoretical limb darkening coefficients for a range of standard broadband filters and instruments.
All Figures
Fig. 1
Temperature distributions (gray shades) on surfaces with constant optical depth at 5000 Å 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 log τ_{5000} = 0 and log τ_{5000} = −2, from which the dominant contribution to the continuum surface brightness between the disk center at μ = 1.0 and the limb near μ = 0.01 is emitted. 

In the text 
Fig. 2
Left: theoretical average surface brightness distribution as a function of projection factor μ integrated over the spectrum between 3400 Å and 3500 Å, computed using our “standard” numerical resolution for an arbitrary snapshot of the 3D model of HD 189733 (crosses), and leastsquares fits of the Claret (2000) nonlinear law (solid line) and a linear law (dotted line). The grayshaded 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 N_{φ} = 8 direction angles (solid line) and with the full horizontal resolution of N_{x} = N_{y} = 240 grid points (dotted line) from the interpolated “standard” setting. 

In the text 
Fig. 3
Left: relative deviation of the splineinterpolated limb darkening predictions for N_{μ} = 17 (solid line) and N_{μ} = 51 (dotted line) from the case with N_{μ} = 101 as a function of μ. Grey vertical lines mark the position of the 17 μ angles used in the “standard” resolution. Right: relative deviation of the leastsquaresfitted limb darkening laws from the computed data points for N_{μ} = 17 (crosses), N_{μ} = 51 (diamonds) and N_{μ} = 101 (triangles) as a function of μ. All calculations are based on the 3D model of HD 189733. 

In the text 
Fig. 4
Left: relative deviation of the splineinterpolated limb darkening predictions computed for a single snapshot of the 3D model for HD 189733 with reduced wavelength sampling of R = 10 000 (solid lines) and R = 5000 (dotted lines) from the full sampling with R = 20 000. Each set of resampled wavelengths is shown in either case, resulting in two or four choices, respectively. Right: relative deviation of the splineinterpolated limb darkening predictions computed for the first 10 snapshots of a time sequence with N_{t} = 20 (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. 

In the text 
Fig. 5
Surface intensities I_{λ}(μ)/I_{λ}(1) 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. Timeaveraged 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 5700 Å, 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 λ/Δλ = 250 for clarity. 

In the text 
Fig. 6
Same as Fig. 5, computed for HD 189733. 

In the text 
Fig. 7
Left: wavelengthintegrated 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 2900 Å and 5700 Å 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. 

In the text 
Fig. 8
Same as Fig. 7 for the case of HD 189733. 

In the text 
Fig. 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 2900 Å and 5700 Å 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. 

In the text 
Fig. 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 2900 Å and 5700 Å 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 5350 Å and 10 500 Å. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.