Issue 
A&A
Volume 574, February 2015



Article Number  A35  
Number of page(s)  22  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201323127  
Published online  20 January 2015 
A nongrey analytical model for irradiated atmospheres^{⋆,}^{⋆⋆}
II. Analytical vs. numerical solutions
^{1} Laboratoire J.L. Lagrange, Université de NiceSophia Antipolis, CNRS, Observatoire de la Côte d’Azur, BP 4229, 06304 Nice, France
email: vivien.parmentier@oca.eu
^{2} Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA
^{3} NASA Ames Research Center, MS2453, Mofett Field, CA 94035, USA
Received: 25 November 2013
Accepted: 29 October 2014
Context. The recent discovery and characterization of the diversity of the atmospheres of exoplanets and brown dwarfs calls for the development of fast and accurate analytical models.
Aims. We wish to assess the goodness of the different approximations used to solve the radiative transfer problem in irradiated atmospheres analytically, and we aim to provide a useful tool for a fast computation of analytical temperature profiles that remains correct over a wide range of atmospheric characteristics.
Methods. We quantify the accuracy of the analytical solution derived in paper I for an irradiated, nongrey atmosphere by comparing it to a stateoftheart radiative transfer model. Then, using a grid of numerical models, we calibrate the different coefficients of our analytical model for irradiated solarcomposition atmospheres of giant exoplanets and brown dwarfs.
Results. We show that the socalled Eddington approximation used to solve the angular dependency of the radiation field leads to relative errors of up to ~5% on the temperature profile. For grey or semigrey atmospheres (i.e., when the visible and thermal opacities, respectively, can be considered independent of wavelength), we show that the presence of a convective zone has a limited effect on the radiative atmosphere above it and leads to modifications of the radiative temperature profile of approximately ~2%. However, for realistic nongrey planetary atmospheres, the presence of a convective zone that extends to optical depths smaller than unity can lead to changes in the radiative temperature profile on the order of 20% or more. When the convective zone is located at deeper levels (such as for strongly irradiated hot Jupiters), its effect on the radiative atmosphere is again on the same order (~2%) as in the semigrey case. We show that the temperature inversion induced by a strong absorber in the optical, such as TiO or VO is mainly due to nongrey thermal effects reducing the ability of the upper atmosphere to cool down rather than an enhanced absorption of the stellar light as previously thought. Finally, we provide a functional form for the coefficients of our analytical model for solarcomposition giant exoplanets and brown dwarfs. This leads to fully analytical pressure–temperature profiles for irradiated atmospheres with a relative accuracy better than 10% for gravities between 2.5 m s^{2} and 250 m s^{2} and effective temperatures between 100 K and 3000 K. This is a great improvement over the commonly used Eddington boundary condition.
Key words: radiative transfer / planets and satellites: atmospheres / stars: atmospheres / planetstar interactions
A FORTRAN implementation of the analytical model is only available 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/574/A35 or at http://www.oca.eu/parmentier/nongrey.
Appendix A is available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
The large diversity of exoplanets in terms of irradiation temperature, gravity, and chemical composition discovered around stars with different properties calls for the development of fast, accurate, and versatile atmospheric models.
In Paper I (Parmentier & Guillot 2014), we derived a new analytical model for irradiated atmospheres. Unlike previous models, our model takes into account nongrey opacities both in the visible and in the thermal frequency ranges. Using two different opacity bands in the thermal frequency range, we highlighted the dual role of thermal nongrey opacities in shaping the thermal structure of the atmosphere. Opacities dominated by lines (i.e., opacities where the lowest of the two values is dominant) enable the upper atmosphere to cool down significantly compared to a grey atmosphere whereas opacities dominated by bands (i.e., opacities where the highest of the two values is dominant) lead to a significant cooling of the upper atmosphere and to a significant heating of the deep atmosphere.
The pressure and temperature dependent linebyline opacities that are used in numerical models to compute accurate temperature profiles are represented in analytical models by only a handful of parameters. Thus, to compute accurate temperature structure from our analytical model for specific planetary atmospheres, we need to know how those parameters vary with the physical properties of the planet.
In this study, we apply our model to irradiated, solarcomposition, semiinfinite atmospheres e.g., brown dwarfs, giant planets or planets with a surface situated in the optically thick region of the atmosphere. Based on the results from a stateoftheart numerical model, we assess the goodness of the different approximations inherent in analytical solutions of the radiative transfer equation. Then, using a grid of numerical models, we calibrate the different coefficients of our analytical model and provide a useful tool for a fast computation of analytical temperature profiles for planet atmospheres that remains correct over a wide range of gravity and irradiation temperatures.
As a first step, in Sect. 3 we quantify the accuracy of models derived with the Eddington approximation, a common simplification of the radiative transfer equation in analytical model atmospheres. Then in Sect. 4 we build a simple radiative/convective model where the radiative solution of Paper I is replaced by a convective solution whenever the Schwarzschild criterion is verified. We then discuss and quantify the intrinsic error of such a simple model of convective adjustment. Then, guided by a stateofthe art numerical integration of the radiative transfer equation, we constrain the parameters of the analytical solution of Paper I to develop a fully analytical solution for the atmospheric temperature/pressure profiles of irradiated giant planets. The solution presented in Sect. 5 reproduces with a 10% accuracy the numerical solutions over a wide range of gravity and irradiations. Finally, in Sect. 6 we highlight the important role of nongrey thermal effects in the influence of titanium oxide on the temperature profiles of irradiated planets.
2. Models
2.1. Setting
We consider the case of a planet with a thick atmosphere (i.e., a planet with no surface or with a surface at very high optical depth) orbiting at a distance a from its host star of radius R_{∗} and effective temperature T_{∗}. The total flux received by the planet is where the equilibrium temperature for zero albedo is defined as (1)At a given point in the planet, the incoming stellar flux is characterized by a temperature (2)This flux hits the atmosphere with an angle θ_{∗} with the vertical direction. We define T_{μ∗ 0} as the projected flux hitting the top of the atmosphere: (3)Where μ_{∗} = cosθ_{∗}. A part 1 − A_{μ∗} of this flux is reflected back to space, where A_{μ∗} is the angledependent reflectivity. We characterize the flux that penetrates the atmosphere by the temperature T_{μ∗}: (4)where .
While these definitions are useful to calculate the temperature profile at a given location in the atmosphere, they need to be averaged over μ_{∗} to calculate the mean state of the atmosphere. As shown by Guillot (2010), the temperature profile obtained for a mean stellar angle of and an incoming flux equal to a fraction of the total incoming flux is a reasonable approximation of the exact mean temperature profile of the planet. When considering average profiles we use the socalled isotropic approximation (e.g., Guillot 2010): (5)where A_{B} is the Bond albedo of the planet. The relationship between A_{B} and A_{μ∗} is not straightforward and will be discussed in more details in Sect. 5.3. f is a parameter smaller than one. For f = 0.25 we have T_{μ∗} = T_{eq} and the thermal profile is close to the planet average profile. For f = 0.5, the thermal profile corresponds to the dayside average profile. The average profiles are obtained by setting and in Eq. (76) of Paper I^{1}.
The atmosphere is also heated by the planet interior and receives a flux from below. The total energy budget of the atmosphere is characterized by its effective temperature (6)which is valid for both the averaged and the nonaveraged cases. Finally, we define an effective temperature for zero albedo: (7)All the quantities defined for zero albedo, denoted by the subscript 0 are independent of the properties of the atmosphere and can be calculated a priori. All other quantities are characteristic of a given atmosphere and need to be constrained either by the observations or determined by numerical calculations.
2.2. Opacities
The interaction between photons and atmospheric gas is described by opacities which are a function of the wavelength of the radiation and of the temperature, pressure and composition of the gas. Although the variety of mixtures and cases to be considered is infinite, we choose to limit the present study to one set of opacities because of its very extensive use both in the context of giant exoplanets and brown dwarfs, i.e., the solarcomposition opacities provided by Freedman et al. (2008). These opacities have been calculated for a solarcomposition mixture in chemical equilibrium. They do not account for the presence of clouds, and any chemical species that condenses at a given temperature and pressure is taken out of the mixture. Although clouds are thought to exist in planetary atmospheres (see Marley et al. 2013, for a review) and should affect the thermal structure of their atmosphere (e.g., Heng et al. 2012), we do not take into account scattering by cloud particles in this study. The first order effect of clouds is to reflect part of the incoming stellar light to the space, what can be taken into account by specifying the relevant Bond Albedo when calculating the effective temperature of the model. Scattering by the gas, however, is taken into account.
While tens of millions of lines have been used for the calculation of these opacities, we choose to show them in Fig. 1 in the same form as they are used by the numerical code described hereafter in Sect. 2.4: in the socalled correlatedk method, the opacities values are sorted from the lowest to the highest values within a limited number of spectral bins (in our case 196), this is similar to the opacity distribution function (ODF) method in the stellar atmospheres modeling (Strom & Kurucz 1966). As long as the spectral bins are small compared to the width of the local Planck function, the error made on the wavelength corresponding to a given opacity is expected to be small and the consequences for the computed temperature profile limited (see Goody & Yung 1989). Figure 1 shows the opacities for different pressure and temperature points taken along a selected planetary temperature/pressure profile corresponding approximately to a solar composition 1Jupiter mass and radius planet at 0.05 AU from a Sunlike star. The wavelength range in which the Planck function has 90% and 99% of the total energy is shown by the thick and thin horizontal bars, respectively, for the different temperatures considered. The contribution of the spectral lines to the opacities shapes the cumulative distribution function inside each bin. As one moves progressively from the top to the bottom of the atmosphere, pressure (always) and temperature (generally) increase which broadens the spectral line profiles. This results in a flattening of the cumulative opacity distribution function within each bin. Opacities in the upper atmosphere are characterized by very strong variations with wavelength and a comblike structure. Deeperdown, the wavelength dependence is mostly due to the presence of molecular bands and takes place on scales significantly larger than our bin size.
Fig. 1 Linebyline opacities as a function of wavelength for five different conditions corresponding to different points in the PT profile of a giant planet with g = 25 m/s^{2}, , T_{int} = 100 K and T_{μ∗} = 1231 K, corresponding to the dayside average profile of a planet orbiting at 0.053 AU from a sunlike star. Inside each bin of frequency, we plot the cumulative distribution function of the opacities instead of the linebyline opacity function. The purple bar represents the wavelength range where 90% (thick line) and 99% (thin line) of the stellar energy is emitted. The other horizontal bars represents the wavelength range of the thermal emission of the planet at different locations along the PT profile. 

Open with DEXTER 
An important feature of these opacities is that most of the variations of the opacity with wavelength take place on scales shorter than the characteristic wavelength range of the Planck function. This is certainly the case at lowpressures when the opacity varies extremely quickly with wavelength, but it remains true (to some extent) at high pressures in the band regime. Another feature of irradiated atmospheres is that the temperature variations remain limited so that there is always a significant overlap between the Planck function from the low to the high optical depth levels. These two features justify the use of the picketfence approximation, and hence of the analytical model of Paper I.
2.3. Analytical model
Although analytical models of irradiated atmospheres can only be obtained for very restrictive approximations on the opacities, they provide nonetheless a useful tool to understand the physics of the radiative transfer and to compute with a low computational cost temperature profiles for a large variety of atmospheric properties. In the particular model derived in Paper I the linebyline opacities are modeled by two different homogeneous set of lines, the full opacity function being described by six independent parameters.
The first set of lines, described by three parameters, represents the thermal part of the opacities, i.e., the part of the opacity function in the frequency range covered by the local Planck function of the atmospheric thermal emission. The Rosseland mean opacity κ_{R}(P,T) is the only one of those parameters that can vary with depth in the atmosphere. In particular, it is the relevant opacity to describe accurately the energy transport in the optically thick part of the atmosphere (Mihalas 1978; Mihalas & Mihalas 1984). The other two parameters describe the nongreyness of the opacities, i.e., their variation in frequency. The first one, γ_{p} is the ratio of the Planck mean opacity to the Rosseland mean opacity, where the Planck mean opacity is dominated by the highest values of the opacities whereas the Rosseland mean is dominated by he lowest values of the opacities. Thus, grey opacities have γ_{p} = 1 and any departure from the grey model increases γ_{p}. The second parameter, β, is the relative width of the opacity lines. Values of β lower than 0.1 represent opacities dominated by atomic lines whereas values of β between 0.1 and 0.9 correspond to opacities dominated by molecular bands. In the following sections, the parameter γ_{p} will sometimes be replaced by an equivalent parameter: κ_{1}/κ_{2}, where κ_{1} is the highest of the two opacities and κ_{2} the lowest. Value of κ_{1}/κ_{2} between 10^{4} − 10^{5} in the upper atmosphere and between 10 − 100 in the deep atmosphere can be estimated from Fig. 1 for a typical hotJupiter. The simple relationship between κ_{1}/κ_{2} and γ_{p} is described by Eq. (87) of Paper I.
The second set of lines represent the visible parts of the opacities, i.e., the part of the opacity function in the frequency range covered by the Planck function of the stellar irradiation. Since the planet’s atmosphere is usually cooler than the stellar photosphere, the two set of opacity lines can be considered independent of each other. Whereas the model cannot take into account more than two thermal opacity bands, it can model as many visible bands as needed. Here we choose to represent the visible opacities with three different bands of adjustable strength represented by γ_{v1}, γ_{v2} and γ_{v3} the ratio of the highest, medium and lowest opacity to the thermal Rosseland mean opacity. Each band is supposed to have the same spectral width described by β_{v1} = β_{v2} = β_{v3} = 1/3.
Although we differentiate the visible and the thermal opacities throughout the paper, the difference is not based on a spectral difference but rather on a geometrical one^{2}. As shown in Paper I the stellar flux is a collimated beam propagating downward in the atmosphere and can be traced down until it is absorbed by the atmosphere. At each level, it’s absorption is proportional to the remaining flux times the local opacities. The mean opacities relevant to understand the absorption of the stellar flux at a given atmospheric level are therefore a combination of the spectrally dependent remaining flux and opacities at this level. The other part of the radiation, sometimes called the diffuse component has a more complex geometry that is often approximated via the Eddington approximation (see Sect. 3). The thermal opacities characterize how the diffuse radiation is absorbed and emitted. Both visible and thermal radiation can be related as they both depend on the physical and chemical properties of the atmosphere but do not have to be equal, even when the stellar flux and the planetary emission overlap in wavelength.
In our analytical model, the Rosseland mean opacity can vary with pressure and temperature. Thus, physical processes producing an overall increase in the opacities, such as the increasing importance of the collision induced absorption with pressure can be accurately taken into account. Our model is the first analytical model to take into account nongrey thermal opacities in irradiated atmosphere. However, the variation of the opacity with frequency cannot change through the atmosphere. Therefore, all the other coefficients must remain constant in the whole atmosphere and a physical phenomenon such as the variation of the pressure or thermal broadening of the lines through the atmosphere cannot be taken into account.
2.4. Numerical model
Whereas analytical models are confined to model atmospheres with very simplified opacities, the radiative transfer equation can be solved by numerical integration using the full, linebyline, frequency, pressureand temperaturedependent opacities described in Sect. 2.2. Moreover, numerical models can integrate the radiative transfer equation by taking into account an arbitrary high number of angular directions, with no need to invoke the Eddington approximation.
Here, we use the extrasolar giant planet (EGP) code initially developed by McKay et al. (1989) for the study of Titan’s atmosphere. Since then, it has been extensively modified and adapted for the study of giant planets (Marley & McKay 1999), brown dwarfs (Marley et al. 1996, 2002; Burrows et al. 1997), and hot Jupiters (e.g., Fortney et al. 2005, 2008; Showman et al. 2009). The version of the code we employ solves the radiative transfer equation using the deltadiscrete ordinates method of Toon et al. (1989) for the incident stellar radiation and the twostream source function method, also of Toon et al. (1989), for the thermal radiative transfer. In some cases incident stellar and emitted thermal radiation bands may overlap, but the radiative transfer is solved separately for each radiation source. Opacities are treated using the correlatedk method (e.g., Goody & Yung 1989). We consider 196 frequency bins ranging from 0.26 to 300 μm; within each bin, the information of typically 10 000 to 100 000 frequency points is compressed inside a single cumulative distribution function that is then interpolated using 8 kcoefficients. The angular dependency is computed using the Gauss quadrature formula for the fluxes. This formula allows us to transform an integral over μ = cosθ into a simple sum over angles (8)with the ω_{i} and the μ_{i} being tabulated in Abramowitz & Stegun (1965). Here we use 5 Gauss points. The EGP model calculates a selfconsistent radiative/convective solution, deriving the adiabatic gradient using the equation of state of Saumon et al. (1995) but can also look for a fully radiative solution.
Although numerical models were built in order to incorporate the full complexity of the opacity function, it can nonetheless solve the radiative transfer equation with the same simplifications than the ones used in the analytical models. In particular, the kcoefficient method can be used to easily implement the simplified opacities of Parmentier & Guillot (2014) by setting a given number of kcoefficients at κ_{1} and the other ones at κ_{2} in each frequency bin. Moreover, the opacities used to compute the absorption of the stellar flux can be independent from the opacities used to compute the thermal fluxes. Simplified opacities as in the analytical case can therefore be used.
Fig. 2 Radiative numerical temperature profile in units of effective temperature as a function of optical depth compared to the analytical solution from Chandrasekhar (1960) using the discrete ordinate in the first, second, third, fourth and fifth approximation. The left panel shows the profiles whereas the right panel shows their relative difference (T_{a}/T_{n} − 1 where T_{a} is the analytical solution and T_{n} is the numerical solution). 

Open with DEXTER 
2.5. Comparison to an asymptotically exact solution
In order to test the validity of the radiative solution found by the numerical model, we compare it to the analytical solution obtained by the method of discrete ordinates in the grey case (Chandrasekhar 1960). This method solves the moment equations with grey opacities for a nonirradiated atmosphere by replacing the integrals over angle by a Gaussian sum. By increasing the number of terms in the sum (i.e., the order of the calculation), it converges towards the exact solution. The first order solution being equivalent to the Eddington approximation.
In Fig. 2, we compare the numerical model for the grey, nonirradiated case to these analytical solutions up to the 5th order. The first order analytical solution deviates from the others and from the numerical result by about 2% with the maximum deviation occurring near optical depth unity. We can therefore expect the analytical models based on the Eddington approximation to differ from the exact results by about this value at least – we will come back to that in Sect. 3. The higher order analytical solutions appear to smoothly converge towards the exact solution, but the numerical solution is found to be about ~0.5% warmer at low optical depths. This discrepancy arises from a different use of the Gaussian quadrature formula in the two approaches. Whereas the analytical solution uses the Gaussian quadrature to compute the integral , the numerical code uses the quadrature formula to compute the flux integral . Therefore, the 5th order analytical solution is formally not the same as the five Gauss points numerical model and does not converge toward the same solution. We tested that using eight Gauss points in the numerical model leads to a solution that is correct to 0.1% when compared to the 8th order analytical solution.
Because a 0.5% error is significantly smaller than the other sources of uncertainties in the model (the first one being due to the use of the Eddington approximation) we chose to use the five Gauss points numerical model. We note that this kind of test is unfortunately not possible in the irradiated case (even in the grey approximation) for which no exact analytical solution is known.
3. Consequences of the Eddington approximation
We have seen in Sect. 2.5 that an asymptotically exact solution of the radiative transfer problem can be found in the grey, nonirradiated case. Unfortunately, such a solution does not exist when accounting for external irradiation. The angle dependency of the radiative transfer problem therefore has to be approximated. Analytical models such as the one in Paper I use a closure relation between two moments of the intensity field I_{ν}(μ) (with ν the frequency of the radiation): (9)This approximation is exact in two specific cases: when the radiation field is isotropic (I_{ν}(μ) = cte ∀μ) and when radiation field is semiisotropic (I_{ν}(μ) = I^{+}∀μ> 0 and I_{ν}(μ) = I^{−}∀μ< 0). In the deep atmosphere, the radiation is quasiisotropic and this approximation holds. Toward the top of the atmosphere, most of the thermal radiation comes from the deep layers and is therefore close to be semiisotropic. In between, the solution is only an approximation. In addition, a boundary condition relating two other moments of the intensity field must be adopted: (10)These two conditions form what is called the Eddington approximation.
In the grey, nonirradiated case, those two approximations are linked and . However using Eq. (9) and imposing leads to the exact solution at the top of the atmosphere, even though it lacks of selfconsistency. In the irradiated case and in the nongrey case the two approximations are independent and is usually set to either 1/2 or , following the grey, nonirradiated case (see Paper I, for a complete discussion).
As discussed in Sect. 2.5, the relative uncertainty on the temperature profile resulting from the Eddington approximation is ~2% in the grey, nonirradiated case. In order to estimate its magnitude in the grey and nongrey irradiated cases, we must rely on comparison with numerical models. We hereafter adopt the EGP numerical model with 5 Gauss points.
Now we compare the radiative solutions from our numerical model and different analytical models using the simplified opacities described in Sect. 2.3. Thus the solution can be expressed as a function of the Rosseland optical depth τ only and is independent of the Rosseland mean opacity or of the gravity. Once normalized by the effective temperature, the temperature as a function of the optical depth in each model only depends on the values of γ_{v},γ_{p} (or κ_{2}/κ_{1}), β and the ratio T_{irr}/T_{int}.
3.1. Irradiated semigrey solutions
In the semigrey case, several analytical models have been developed (e.g., Hansen 2008; Guillot 2010; Robinson & Catling 2012). As reviewed in Paper I, those models differ mainly by their choice of and their choice of the upper boundary condition. For simplicity, we will compare only three of them: the two different versions of Eq. (27) of Guillot (2010; with f_{H} = 1/2 or ) and the semigrey limit of the model derived in Paper I with which uses as upper boundary condition a mix between the model of Guillot (2010) and the one of Hansen (2008). We compare those models for two different values of the main parameter of semigrey models: the ratio of the visible to the thermal opacities, γ_{v}.
Fig. 3 Comparison between the radiative numerical solution (black line), our work (red line) and Guillot (2010) model for two different values of f_{H} (blue and green lines) for a fully radiative semigrey atmosphere with γ_{v} = 0.25 (plain lines) or γ_{v} = 10 (dashed lines). We set , T_{irr} = 1288 K and T_{int} = 500 K. 

Open with DEXTER 
Figure 3 compares these models for a typical irradiated Jupitermass exoplanet close to a solartype star and shows the magnitude of the error which is due to the Eddington approximation – both the closure relation defined by Eq. (9)and the adopted value of – and the chosen upper boundary condition, different between Guillot (2010) and Paper I. The left panel shows the temperature profiles as a function of optical depth which mainly depends on the magnitude of the greenhouse effect: when γ_{v} is small, most of the incoming irradiation is absorbed deep in the atmosphere, the temperature increases monotonously with increasing depth, and the solution behaves like the nonirradiated solution with the same effective temperature (see the 1st order case of Fig. 2). When γ_{v} is large, most of the incoming stellar light is absorbed high up, creating a temperature inversion around visible optical depth unity (and thus thermal optical depth τ = 1 /γ_{v}).
The right panel of Fig. 3 shows that the magnitude of the difference between the numerical solution and the analytical ones strongly depends on the choice of and of the top boundary condition, but remains of the same orderofmagnitude as for the nonirradiated grey case of Sect. 2.5. Specifically, the Eddington approximation is found to lead to a ~4% uncertainty on the temperature profile and always converges towards zero at large optical depths. Except for the solution, all other analytical solutions, including the one from Paper I, systematically underestimate the temperature at a given depth.
It is obvious from Fig. 3 that, unlike the nonirradiated case, no choice of can yield an exact skin temperature T(τ = 0) (see related discussion in Paper I).
3.2. Irradiated nongrey solutions
We now test the analytical model of Paper I in the nongrey case. In order to do so, we compare the analytical model to the numerical model for different values of the ratio of the thermal opacities κ_{2}/κ_{1} and a single visible channel to the numerical model with the same thermal and visible opacities. We adopt β = 0.86 and γ_{v} = 0.25, typical values needed to reproduce detailed models of hot Jupiters (see Sect. 5) and the same irradiation and internal temperature as in the previous section.
Fig. 4 Comparison between the analytical model (plain lines) and the radiative numerical model (dashed lines) for different values of κ_{2}/κ_{1} (left panel). The right panel shows the relative difference between the analytical and the numerical solution for each case. We used , T_{irr} = 1288K, T_{int} = 500K, γ_{v} = 0.25 and β = 0.86. 

Open with DEXTER 
Figure 4 shows the resulting temperatureoptical depth profiles and the relative difference between the numerical and analytical solutions. As κ_{2}/κ_{1} increases, the temperature profile gets cooler in the upper atmosphere and warmer in the deep atmosphere, an effect described in details in Paper I.
The red curve (κ_{2}/κ_{1} = 1) corresponds to the semigrey solution already seen in Sect. 3.1 and Fig. 3. As shown in the right panel, the discrepancy between the analytical and numerical models increases with the nongreyness of the opacities. The maximum error (in absolute terms) increases from about ~2% to a little bit less than ~5% when κ_{2}/κ_{1} is increased from 1 to 10^{5}. Moreover, the optical depth range for which the discrepancy is larger than 1% increases at the same time from [0.1 ~ 10] to [10^{5} ~ 10^{3}].
This increase in the extent of the region in which the temperature profile departs from the numerical solution is a direct consequence of the Eddington approximation in the two thermal channels with opacities κ_{1} and κ_{2}, respectively. At high optical depth, in the diffusion limit, the radiation field is isotropic in each thermal channel and the Eddington approximation is valid. At very low optical depth radiation comes mostly from the levels where the first and the second thermal channels become optically thin, much deeper in the atmosphere. Therefore radiation in the optically thin layers is close to be semiisotropic which validates the choice of the Eddington approximation. Inbetween, the difference between the analytical and the numerical solutions exhibits two maxima. Those maxima correspond to the levels where the first and the second thermal bands become optically thin. As the ratio κ_{1}/κ_{2} increases, the first channel becomes optically thin at higher Rosseland optical depth and the second channel becomes optically thin at lower Rosseland optical depth, creating the twopeak feature of Fig. 4.
We see however that the error induced by the Eddington approximation remains lower than 5%, with the deep temperatures being colder in the analytical model than in the numerical model. Compared to other sources of uncertainty (in particular our assumptions that β and κ_{2}/κ_{1} are uniform in the atmosphere), this is an acceptable level of uncertainty.
4. Consequences of convection on the overlaying radiative solution
At highenough optical depth, the deep atmospheres of giant planets and brown dwarfs become convective (e.g., Guillot 2005), a consequence of the increase in the opacity with pressure (see Rauscher & Menou 2012). This increase in the opacity in substellar atmospheres is due both to collisioninduced absorption by hydrogen molecules increasing with density (above roughly 10^{3} g cm^{3}) and eventually to new opacity sources linked to a larger abundance of electrons at temperatures ~2000 K and above. Generally, exoplanets and brown dwarfs with lowirradiation levels (i.e., such that T_{irr}_{~}^{<}T_{int}) have a convective zone extending all the way from the deep interior to the τ ~ 1 optical depths. This is for example the case of Jupiter, whose atmosphere becomes convective at pressures of order P ~ 0.3 bar – but with considerable heterogeneity depending on the latitude and longitude on the planet (e.g., Magalhaes et al. 2002; West et al. 2004). However, in very closein exoplanets and brown dwarfs, the high stellar irradiation maintains the atmosphere in a very hot state and pushes the radiative/convective transition down to very high optical depths (see Guillot et al. 1996; Guillot 2005).
Numerical models naturally account for these convective zones by imposing a temperature gradient set by convection when a condition such as the Schwarzschild or Ledoux criterion is met. The temperature profile in the radiative part(s) of the atmosphere is then recalculated taking into account the presence of convective zone(s). The procedure is applied iteratively until a full radiative/convective equilibrium is reached. While it is easy to implement the first condition (imposing a convective gradient when necessary) in analytical atmospheric models, it is generally not possible to implement the second one and modify the radiative solution because of the presence of a convective region. In the specific case of the grey and semigrey models, Robinson & Catling (2012) recently derived a radiativeconvective model that satisfies these two conditions, although it necessitates a small numerical integration. For nongrey thermal opacities, no analytical model solves selfconsistently for the convective and the radiative parts of the atmosphere. In the specific case of the model of Paper I, the boundary condition of the radiative atmosphere lays in the optically thick layers and a the solution cannot be modified to account for a change in the temperature gradient at deep levels.
We want to estimate the error made when the presence of a convective zone is neglected when calculating the temperature profile in the radiative parts of the atmosphere. We build our analytical radiative/convective model by switching from our radiative solution to the adiabatic solution whenever the convective gradient becomes lower than the radiative one. We compare the resulting analytical solution to the numerical solution in which both the depth of the radiative/convective boundary and the atmospheric temperature profile are converged iteratively. As the presence and depth of a convective zone depends on the exact value of the opacities, we need to specify the Rosseland mean opacity in our model. In this section, in order to facilitate the comparison, we fix the Rosseland mean opacity as a function of pressure to its value in our fiducial model, described in Fig. 1. However, all our results will be relative to the depth of the convective zone and thus independent from the exact Rosseland mean opacity function used.
We only consider the case for which the atmosphere transitions from being radiative at high altitudes to being convective at depth (i.e., we do not include the possibility of alternating radiative and convective zones). In the convective zone, we assume that the temperature gradient is exactly adiabatic (i.e., we do not account for the superadiabatic gradient required to transport the heat flux – e.g., Guillot (2005)).
4.1. Nonirradiated grey case
We first compare the solutions obtained in the nonirradiated grey case. In order to see how the location of the radiative/convective zone influences the solutions, we artificially modify the adiabatic gradient by a factor that varies from 1/4th to 4.
Fig. 5 Comparison of our numerical and analytical radiativeconvective models for different adiabatic gradients in the nonirradiated, grey case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{int} = 500 K and g = 25 m/s^{2}. Note that the cases ∇_{ad} × 2 (green) and ∇_{ad} × 4 (red) are superimposed. 

Open with DEXTER 
When the radiative/convective transition occurs below optical depth unity (red, green, blue and purple curves in Fig. 5), the difference between the analytical and numerical solutions is unchanged (the corresponding curves are indistinguishable on the right panel) and entirely due to the Eddington approximation as discussed in the previous section. This error is frozen at the radiative/convective boundary and propagates in the convective zone leading to an estimate of the deep temperature profile that is at most 2% percent off. For a convective zone that crosses the τ ≈ 1 limit (orange and black curves of Fig. 5), the lower boundary condition used in the analytical radiative model – that the deep atmosphere reaches the diffusion limit – is no more valid. The error becomes dependent on the location of the radiative/convective transition (and value of the adiabatic gradient). It however remains of the same order as the one due to the Eddington approximation. This validates models calculating the radiative/convective boundary of the deep convective zone without recalculating the upper radiative profile. However, the presence of detached convective zones cannot be modeled correctly with this method, and an approach similar to Robinson & Catling (2012) is needed.
4.2. Nonirradiated nongrey case
We now turn, with Fig. 6, to the nonirradiated nongrey case, using the fiducial values κ_{2}/κ_{1} = 10^{2} and β = 0.83. As in the grey case, the errors are dominated by the Eddington approximation as long as the radiative/convective boundary occurs at optical depths larger than unity in the two thermal channels that are considered. The error at low optical depths is larger because the error due to the Eddington approximation is greater in the nongrey case. However, as soon as the convective zone extends to levels of optical depth unity or smaller, the discrepancy between the analytical and numerical solutions increases significantly: the upper atmosphere warms up and our analytical solution is no more a good representation of the radiative atmosphere. This is clearly due to a nongrey effect. In a given spectral interval, the thermal flux present at a given pressure is set by the integrated thermal emission of all the atmospheric layers below it in this specific spectral interval. At large optical depth, the emission is thermalized and the thermal flux per wavelength emitted in both spectral channels is the same regardless of the temperature gradient. At optical depth close to unity, the thermal flux in each channel depends on the actual temperature gradient. The analytical solution assumes that the temperature gradient is set by radiation transport everywhere and thus calculates inaccurately the flux emitted in the two spectral channels if convection extends to optical depths smaller than unity. The resulting temperature profile can differ by tens of percentage points from the numerical one. In addition, because the relative error is frozen at the one obtained at the radiative/convective transition, it does not tend towards zero at large optical depths as was the case with the purely radiative solutions.
Considerable caution should therefore be exerted when switching from radiative to convective gradient without recalculating the radiative solution in the general (nongrey) case. Specifically, when the atmosphere becomes convective at optical depths smaller than unity, the resulting temperature profile may be inaccurate by several tens of percentage points.
Fig. 6 Comparison of our numerical and analytical radiativeconvective solutions for different adiabatic gradients in the nonirradiated, nongrey case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{int} = 500 K, g = 25 m/s^{2}, κ_{2}/κ_{1} = 10^{2}, and β = 0.83. The cases ∇_{ad} × 2 (green) and ∇_{ad} × 4 (red) are superimposed. 

Open with DEXTER 
4.3. Irradiated nongrey case
We now consider the effect of irradiation with our fiducial hot Jupiter atmosphere. As already discussed, the strong irradiation tends to push the radiative/convective zone towards deep levels (see Guillot 2005). This is seen in the profiles of Fig. 6 for which the radiative/convective boundary always occur at optical depths ~100 or deeper, with only a small dependence on the value of the chosen adiabatic gradient. As expected, this suppresses the changes in the temperature profile in the purely radiative atmosphere. The errors are almost independent of the assumed adiabatic gradient and mostly due to the Eddington approximation. For hot Jupiters, and generally for strongly irradiated atmospheres, the presence of a deep convective zone may be accounted for by adopting a purely radiative solution and switching to the convective one when the Schwarzschild criterion is verified.
Of course, for a smaller irradiation level and/or larger values of the κ_{2}/κ_{1} ratio, the presence of a convective zone reaching optical depths closer to unity (in one of the thermal channels at least) will lead to an increase on the error of the calculated temperature profile. We expect this error to be approximately bounded by that of the nongrey, nonirradiated case.
Fig. 7 Comparison of our numerical and analytical radiative/convective solutions for different adiabatic gradients in the nongrey, irradiated case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{irr} = 1288 K, T_{int} = 500 K, γ_{v} = 0.25, , κ_{2}/κ_{1} = 10^{2} and β = 0.83. 

Open with DEXTER 
5. Modeling the nongrey effects
Analytical model atmospheres are useful to understand the key physical processes at stake in planetary atmospheres. Unfortunately they cannot take into account the complex variations of the opacities with frequency, temperature and pressure. However, when modeling a specific planet atmosphere with a given chemical composition, the knowledge of the linebyline opacities should guide one in a proper choice of parameters when using the analytical models. In this section we wish to understand what characteristics of the opacities shape the temperature/pressure profile of a planet atmosphere and find a method to derive the simplified opacities of our analytical model from the linebyline opacities. Ideally, the resulting analytical temperature/pressure profile should be a good approximation of the numerical solution computed with the full frequency, temperature and pressure dependent opacities.
A first approach to determine our coefficients is an a posteriori determination i.e., to choose the coefficients such that the analytical and the numerical profiles match correctly. Although this should give the best results in terms of goodness of the fit, the retrieved coefficients might not be physically realistic and it could be difficult to relate them to the real atmospheric opacities. Another approach is to find a priori values, directly from the opacities. This requires a deep understanding of the opacities and how they shape the temperature profile. A last possibility is to combine the two approaches: using an a priori determination when possible and adjusting the remaining coefficients a posteriori to fit the numerical profile.
5.1. A priori determination of the coefficients
5.1.1. Visible coefficients
The visible coefficients control at which depth the stellar flux is absorbed in the atmosphere. When the visible absorption is strong, the stellar flux is absorbed in the upper part of the atmosphere and radiated back to space. At the opposite, when the visible absorption is weak, the incoming irradiation is deposited at depth where the thermal optical depth is large and the deep atmosphere warms up. This is the wellknown greenhouse effect.
When taking into account only one visible band (i.e., β_{v} = 1, as in Guillot 2010), a natural choice for the parameter γ_{v1} is the ratio of the mean Rosseland visible opacity (using the Planck function to weight the linebyline opacities) to the mean Rosseland thermal opacity (using the local Planck function at the stellar effective temperature to weight the linebyline opacities). Unfortunately, this ratio can vary significantly with height. We find that choosing the ratio at τ_{v} = 2/3 (where τ_{v} is the Rosseland visible optical depth) leads to a correct representation of the absorbed stellar flux at depth and could be used, together with a proper modeling of thermal nongrey effects, to get a first guess of the deep temperature. The reason is that the part of the stellar flux that heats up the deep atmosphere is the one that propagates down to the τ> 1 level. Thus, the opacities that determine the relevant strength of the visible absorption are the lowest visible opacities. The Rosseland mean is a good estimate of the weakest opacities over a given frequency range and is thus a suitable estimate.
However, when a significant portion of the stellar radiation is absorbed in the upper atmosphere of the planet, for exemple when strong visible absorbers such as titanium oxide or sodium are present in the atmosphere, the stellar flux that reaches the τ> 1 level depends strongly on the amount of absorption in each spectral channels in the upper atmosphere. The knowledge of γ_{v} at a given level is not sufficient anymore for a correct estimate of the deep temperature. A more sophisticated model of the visible absorption is then needed.
As shown in Paper I, the moment equations are linear with respect to the absorbed stellar flux. Thus our model can take into account as many spectral bands in the visible as needed with the condition that the different values γ_{vi} = κ_{vi}/κ_{R} in each visible bands are constant through the atmosphere. If well chosen, constant nongrey visible opacities can relatively well approximate the absorbed stellar flux at all atmospheric levels. We therefore adopt the following method: using the linebyline opacities from Freedman et al. (2008) and the actual numerical PT profile, we calculate the total absorbed flux at each layer of the atmosphere. We then adjust the relative contributions of the different visible opacity bands in order to correctly match the absorbed visible flux from the numerical simulation. The stellar flux absorbed by n spectral bands of width β_{vi} is (11)where the visible bands are homogeneously distributed in frequency (similar to the thermal bands), F_{0} is the total incident stellar flux and the β_{vi} must verify: ∑ _{i}β_{vi} = 1. We apply this method using one to four opacity bands. As seen in Fig. 8, the absorbed flux is poorly represented when considering only one opacity band. When using two, three or four bands, the absorbed flux is described with a 4%, 1%, and 0.5% accuracy, respectively. In the following we limit ourself to three opacity bands of constant width β_{vi} = 1/3 for the visible opacity.
Figure 9 compares the numerical model in black, taking into account all the linebyline opacities and the semigrey model (blue line) where the visible opacities are adjusted in order to have the same absorbed flux as in the numerical model but where the thermal opacities remain grey. The semigrey model, even though it models correctly the absorbed flux as a function of depth, lays far from the numerical solution. Clearly, nongrey thermal opacities are needed.
Fig. 8 Absorbed stellar flux from the numerical model (dots) and from the analytical model (lines) considering respectively 1, 2, 3 or 4 absorption bands in the visible for our fiducial model atmosphere (see Fig. 1). 

Open with DEXTER 
Fig. 9 Pressuretemperature profiles calculated using the numerical model and the full set of opacities (black), the semigrey (blue) and the nongrey (red) analytical models. As un Fig. 1, g = 25 m/s^{2}, , T_{int} = 100 K and T_{μ∗} = 1253 K. The coefficients used for the analytical models are taken from Table 1 with γ_{p} = 1 for the semigrey case. The nongrey model is a much better match to the numerical profile than the semigrey one. 

Open with DEXTER 
5.1.2. Thermal coefficients
Fig. 10 Cumulative distribution function of the opacities at P = 0.85 bar and T = 1464 K, corresponding to the τ = 2/3 level of an atmosphere with T_{μ∗} = 1253 K and a gravity g = 25 m/s^{2}. The Yaxis represents the fraction of frequency where the monochromatic opacities are lower than the corresponding κ_{0} of the Xaxis. The red line shows the value of the Rosseland mean opacity and the black line the Planck mean opacity. 25% of the frequency range have monochromatic opacities smaller than the Rosseland mean opacity whereas 90% have monochromatic opacities smaller than the Planck mean opacity. 

Open with DEXTER 
Fig. 11 Top panel: coefficients γ_{v3}, γ_{v2}, γ_{v1}, β, and γ_{p} obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets of solar composition with different gravities and an internal temperature of 100 K. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4} bar <P< 10^{2} bar, the second one for 10^{2} bar <P < 10^{0} bar and the third one for 10^{0}bar <P< 10^{2} bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/g) ≲ τ_{R} ≲ 10 (25 ms^{2}/g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/g). 

Open with DEXTER 
The thermal coefficients describe how well the atmosphere is able to retain its energy. As explained qualitatively by Pierrehumbert (2010) and quantitatively in Paper I, the presence of nongrey thermal opacities can strongly affect the temperature profile of the planet. Because it is tied to the emission and absorption of the thermal flux, only the opacity variations that have an extent smaller or comparable to the local Planck function can contribute to the nongrey effects. The cumulative distribution function of the opacities in the frequency range covered by the local Planck function should thus contain enough information to constrain the nongrey effects. As a grey atmosphere cools down principally by emission from the τ = 2/3 level, the thermal opacities at this level should determine the strength of the nongrey effects.
We plot in Fig. 10 the cumulative distribution function of the opacities at the τ = 2/3 level for our fiducial model atmosphere. It represents the relative spectral width over which the opacities are lower than a given opacity κ_{0} as a function of κ_{0}. The opacities cover a wide range of value (6 orders of magnitude in the specific example shown in Fig. 10). Our analytical model can describe the nongrey thermal opacities with only two parameters: the ratio of the Planck mean opacity to the Rosseland mean opacity, γ_{p} and the relative size of the two bands, β. Unlike in the visible case, the thermal effects are local effects that do not depend on the behavior of the rest of the atmosphere. The value of γ_{p} can hence be calculated as a function of pressure and temperature from tables available in the community (e.g., Freedman et al. 2008).
The parameter β describes the relative amount of the opacities which are in the first band compared to the second band. The Rosseland mean opacity is determined by the smallest values of the opacities, which is the second band opacity in our model. We decide to use as β the fraction of the opacities in the spectral range covered by the local Planck function that are higher than the Rosseland mean opacity. This can be derived directly from the cumulative distribution function of the opacities plotted in Fig. 10. In the specific example of Fig. 10, 25% of the opacities lay below the Rosseland mean opacity hence β = 0.75.
5.2. Application/different models
In order to test the goodness of our analytical model and derive reasonable estimates of the coefficients, we use the EGP numerical code to build a grid of atmospheric radiative/convective models for giant planets with a solar composition atmosphere, three different gravity (2.5, 25, and 250m/s^{2}), and an internal temperature of T_{int} = 100 K. We consider the case of a planet orbiting a sunlike star at various distances corresponding to irradiation temperatures from 100 K to 3000 K. All the profiles were calculated using and therefore represent global or dayside averaged temperature profiles (see Sect. 2.1 for more details). Figure 11 shows different models obtained for different estimates of our coefficients (top panel) and a comparison between the numerical profiles and the resulting analytical profiles (bottom panel). In all models but model D, we use as Rosseland mean opacity the one calculated by the numerical model directly from the linebyline opacities. In model D, we use the functional fit of the Rosseland mean opacities of Freedman et al. (2008) provided by Valencia et al. (2013). Similarly, model D uses a functional form for the convective gradient when convective adjustment is necessary (see Sect. 5.4). Model D is therefore a fully analytical model that can be downloaded and implemented by the community. We now describe the different models.
Model A: in this model, we adjust all our coefficients a posteriori in order to have the best match to the numerical profiles. It leads to temperature profiles in agreement within 5% with the numerical ones. It shows that our analytical model can represent a large variety of atmospheric temperature profiles and goes beyond the limitation of previous semigrey models (Parmentier & Guillot 2014). However, the spread of the retrieved value of the coefficients makes it difficult to derive a trustable functional form and a better approach is needed in order to get a fully analytical model.
Model B: here we use the methods of Sect. 5.1 to determine a priori the various coefficients. The visible coefficients, γ_{v1}, γ_{v2}, and γ_{v3} are determined respectively by the first, the second, and the third thirds of the incoming stellar energy that are absorbed by the atmosphere as the stellar irradiation propagates downward. Thus γ_{v1} is determined by the highest values of the visible opacities, γ_{v2} by the median ones and γ_{v3} by the lowest ones. All of them are also proportional to the inverse of the thermal Rosseland mean opacities. They exhibit four different types of behavior as a function of the effective temperature. Those behaviors reflect changes in chemical composition with the irradiation temperature (a plot of the linebyline opacities for the four different regimes is shown in Appendix A):

For T_{eff}< 250 K, the visible opacities are dominated by Rayleigh scatteringand exhibit a slight dependence with gravity. As the effectivetemperature increases, more absorbers are present in gaseousform and the visible opacities increase.

For 250 K <T_{eff}< 600 K, the visible opacities are dominated by the sodium lines at great depth, where the profile is warm enough to have sodium in gaseous state, whereas it is dominated by much weaker lines in the upper atmosphere. As T_{eff} increases, the lines broaden, increasing the lowest values of the opacities and decreasing the highest values of the opacities. As a consequence, γ_{v2} and γ_{v3} increase with T_{eff} whereas γ_{v1} decreases.

For 600 K <T_{eff}< 1700 K the sodium and potassium become the main gaseous absorbers in the upper atmosphere, leading to a strong visible absorption and an increase in the visible opacities with T_{eff}.

For T_{eff}> 1700 K titanium and vanadium oxides become the main gaseous absorbers in the upper atmosphere, creating a substantial increase in the visible opacities. This is reflected by the sudden increase in γ_{v1} and γ_{v2}.
The thermal coefficients have a smoother behavior with T_{eff}. β is rather constant and equal to ≈0.8. This high value of β can be interpreted as a predominance of the molecular bands (i.e., the water and methane bands) to the atomic lines in the thermal opacities. Around T_{eff} = 300 K, β reaches values even closer to 1. At these temperatures, the Planck function at the atmospheric levels of τ ≈ 2/3 overlaps with the 5 μm window in the opacities which is consistent with large values of β. The variations of γ_{p} with T_{eff} have a saddlelike shape with two maxima at 300 K and 1500 K. At low temperatures, the Planck function of the atmosphere shifts towards large wavelengths (>10 μm) for which the opacities are almost constant, leading to small values of γ_{p}. At very high T_{eff}, the Planck function of the atmosphere shifts toward smaller wavelengths (<1 μm) for which the TiO broadband absorption significantly flattens the opacities, leading to small values of γ_{p}. In between, when the atmospheric Planck function is between 1 and 10 μm, the opacities are dominated by the water and methane bands, which raises the value of γ_{p} up to ≈100.
Fig. 12 Comparison between the numerical solutions (dashed lines) and the analytical solutions of model D (using the functional form of the coefficients given in Table 1) over a wide range of irradiation temperatures for a giant planet of solar composition orbiting a sunlike star. Here we used g = 25 m/s^{2}, T_{int} = 100 K, and . Thick and thin lines represent convective and radiative zones, respectively. 

Open with DEXTER 
Although this model gives a correct estimate of the profile at high pressure, it leads to errors of ≈40% at medium and low pressures. Given that the coefficients were all determined a priori, reaching a 40% accuracy can be a fair, first guess of the temperature profile. This method could be extended to planets with very different opacities without going through the whole numerical integration of the radiative transfer equation. However, as proven by model A, a much better accuracy can be obtained by the analytical profile and a mixed method with some coefficients derived a priori and others a posteriori can be a good compromise.
Model C: in this model we use a mixed method to derive the coefficients of the analytical model, with some of them being derived a priori and some of them a posteriori. The method to determine the visible coefficients seems robust, as it can give the correct absorbed flux as a function of optical depth in the atmosphere with a 1% accuracy. The method to determine the thermal coefficients is more subject to caution as it is unclear whether the value of γ_{p} in our analytical model should correspond to the value of γ_{p} derived from the real opacities. Moreover, our criteria to choose β (the fraction of the opacities that are higher than the Rosseland mean opacity) is ad hoc and does not rely on strong physical arguments. At last, there is no strong argument to choose the depth at which those coefficients are calculated. We thus decided to obtain the visible coefficients from the a priori solution and to fit the thermal ones by adjusting the analytical profile to the numerical profile. The resulting analytical solutions lead to an estimate of the temperature profile that always differs by less than 10% from the numerical solution.
Compared to model B, only the thermal coefficients are changed in model C. Below 200 K, γ_{p} is small and β is not well defined. For 200 K <T_{eff}< 2000 K, β is roughly constant with T_{eff} with values of approximately 0.8, in agreement with the a priori determination in model B. For T_{eff}> 2000 K, β decreases slightly with T_{eff}, showing that nongrey effects become more important in the upper atmosphere than in the deep atmosphere.
γ_{p}keeps the same dependency with T_{eff} than in model B but is one order of magnitude smaller. The high values of γ_{p} derived in model B seem relevant to understand the deep atmospheric structure but are too high to represent the upper atmosphere. Large values of γ_{p} lead to a stronger cooling of the upper atmosphere and in temperatures cooler than expected from the numerical model.
Model D: in order to produce a fully analytical model we fit a functional form to the coefficients derived in model C as a function of T_{eff}. Following the different regimes that we described in model B, we fit different affine functions to the visible coefficients. Although the transitions at ≈250 and ≈1800 K are discontinuous, we decided to provide a continuous fit by introducing two different transition zones for 200 K <T_{eff}< 300 K and 1400 K <T_{eff}< 2000 K. Ensuring continuity in the fit is important to increase the numerical convergence of numerical models where this fit can be used as a boundary condition (e.g., internal structure and evolution models of giant planets). The thermal coefficients having a much smoother variation with the irradiation temperature, we use a constant value with an affine function for T_{eff}> 2000 K to represent the parameter β and we use a 2nd order polynomial to represent γ_{p} over the whole temperature range. The functional form of the coefficients are presented in Table 1.
Functional form of the planet planeparallel albedo log _{10}(A_{μ∗}) = a + bX where X = log _{10}(T_{eff0}) and g is the gravity of the planet in m/s^{2}.
The resulting model matches the numerical profiles over a wide range of irradiation temperature and planet gravity with an accuracy always better than 10% for pressures ranging from 100 bar to 10^{2} bar and generally better than 10% for pressures ranging from 10^{2} bar to 10^{4} bar (see Fig. 12 and Col. D of Fig. 11). The fit is slightly worst at low T_{eff} where we do not model the gravity dependance of the visible coefficients and at the transition zone 1400 K <T_{eff}< 2000 K where we model as continuous a discontinuous transition. Whereas for previous models the Rosseland mean opacities used in the analytical solution were calculated by the numerical model directly from the linebyline opacities, Model D uses the fit of the Freedman et al. (2008) Rosseland mean opacities provided by Valencia et al. (2013). Moreover, model D includes a functional form of the convective gradient (see Sect. 5.4 for more details). This makes model D a selfconsistent fully analytical model.
Model E: here, the importance of nongrey effects are tested. This model has grey thermal opacities (i.e., γ_{p} = 1) but uses the functional form derived in model D for the visible coefficients. Therefore, model E is a good representation of the absorption of the stellar irradiation by the atmosphere but lacks the nongrey effects. At medium and large pressures, model E gives a reasonable estimate of the temperature profile, although worst than model D. At low pressures model E provides a very bad estimate of the temperature, with discrepancies of up to a factor of two with the numerical model. The nongrey absorption of the stellar irradiation therefore cannot by itself explain the temperature structure in planetary atmospheres. Nongrey thermal effects, such as the ones considered in model D, are necessary.
Model F: in this model, a comparison with the previous estimate of Guillot (2010) is done. We use grey thermal opacities and the visible coefficients provided by Guillot (2010). Those coefficients were derived in order to match the deep temperature of highly irradiated planets, which it does well. However, at low pressures model F always fails to represent the numerical temperature profile and, at temperatures smaller than 2000 K the discrepancy can reach 40% at all atmospheric pressures.
5.3. Albedo
Although our analytical model has been derived for a purely absorbing atmosphere, we indirectly take into account scattering via the value of the Bond albedo. Our numerical model calculates the reflectivity A_{μ∗} of a planeparallel atmosphere irradiated with an angle . Exact calculation of the Bond albedo of a planet involves an integration of the radiative transfer equation for numerous angles of the incident stellar beam and is beyond the scope of this paper (e.g., Cahoy et al. 2010). We hereafter approximate the Bond albedo of the planet by the planeparallel albedo calculated for . This is similar to the socalled isotropic approximation for the mean temperature profile (e.g., Guillot 2010). Although the approximation is a strong one, we expect A_{B} and A_{μ∗} to have similar values and to follow the same trends.
Fig. 13 Planeparallel albedos for solarcomposition, clearsky atmospheres with different values of gravity and internal temperature, numerical model (dots) and analytical fit described in Table 2 (lines). In the inset, the albedo plotted in linear scale. Here we approximate the Bond albedo of the planet by the planeparallel albedo. 

Open with DEXTER 
As shown in Fig. 13, hotter planets have lower albedos. We consider cloudless atmospheres, neglecting all types of clouds that might be present in this temperature range. The derived albedos are therefore lower bounds as any additional diffusion by clouds should increase the planet albedo. Here, the albedo is set by the competition between Rayleigh scattering and atomic/molecular absorption. Hotter planets have larger abundances of gaseous absorbers in their atmospheres and have therefore a lower albedo. At high effective temperatures (T_{eff} ≳ 1250 K), titanium and vanadium oxides broadband absorption together with collisioninduced absorption dominate the optical opacities, leading to Bond albedos lower than 0.1. As the effective temperature decreases, TiO and VO become less abundant and the albedo increases. From 750 K ≲ T_{eff} ≲ 1250 K, the visible opacities are dominated by the absorption lines of sodium and potassium, leading to a plateau in the value of the Bond albedo. At lower effective temperatures, the sodium and the potassium slowly disappear from the atmosphere and, at T_{eff} ≲ 250 K, the albedo reaches its maximum value.
We also see in Fig. 13 that planets with higher gravities have lower albedos. This trend is particularly important around T_{eff} ≈ 1000 K. At these effective temperatures, absorption in the atmosphere is dominated by the pressurebroaden sodium and potassium lines. The pressure of the layer where the photons at a given wavelength are deposited is proportional to the gravity of the planet. For higher gravity planets the stellar photons reach deeper layers where the alkali lines are broader and the absorption higher. Since Rayleigh scattering is independent of pressure, planets with higher gravities absorb more efficiently the stellar irradiation and have lower albedos.
Interestingly, the albedo is very insensitive to the internal temperature of the planet. For a given effective temperature, changing the internal temperature does not change the temperature profile enough to modify the chemical equilibrium of the atmosphere and change its composition.
Our albedos agree well with the calculations of Marley & McKay (1999) for T_{eff} = 200 K. For T_{eff} = 1000 K our albedos are much smaller than the ones of Marley & McKay (1999). As shown in Sudarsky et al. (2000) alkali metals, not considered in the previous study increase significantly the absorption at those effective temperatures and explain the discrepancy.
In Table 2 we provide a fit of the albedos for different effective temperature and gravities. This fit can be used together with Model D of Sect. 5.2 to obtain a selfconsistent model of irradiated planet atmospheres.
5.4. Radiative/convective model
As discussed in Sect. 4, the deep atmosphere of substellar objects is convective and the temperaturepressure profile can be determined by integrating the adiabatic gradient: (12)This gradient is an intrinsic property of the fluid and can be derived from the equation of state. Based on a fit of the Saumon et al. (1995) equation of state at high pressures, we use the functional form: (13)Our radiative/convective model consist of a convective zone from the bottom of the model up to the radiative/convective boundary at a pressure P_{R/C}. For P<P_{R/C} the analytical model with the coefficients of model D is used. To choose P_{R/C} we compare the radiative gradient and the convective gradient given by Eq. (13)from the bottom to the top of the model. As long as the convective gradient is smaller than the radiative gradient the atmosphere is convective. Whenever the radiative gradient becomes smaller than the convective one the radiative/convective boundary is reached.
Although this method should provide a good estimate of P_{R/C}, it fails at low effective temperatures when used with model D. At low effective temperatures, a radiative zone squeezed between two convective zones develops around P ≈ 100 bar. In the numerical model, the two convective zones merge and the radiative zone disappears for T_{eff}< 200 K. In the analytical model, however, the radiative zone remains and our estimate of P_{R/C} is biased. This discrepancy is due to the accuracy of the fit of the Rosseland optical depth provided by Valencia et al. (2013). At P ≈ 100 bar and T ≈ 1000 K the fit systematically underestimate the Rosseland mean opacities by ≈ 30 − 50% (see Fig. 8 of Freedman et al. 2014). This underestimation of the Rosseland mean opacity directly translates to an underestimation of the radiative gradient and to the appearance of a spurious deep radiative zone and a bias in our estimate of P_{R/C}. In order to recover the behavior of the numerical model, we consider that a radiative zone that is in between two convectively unstable zone must have a radiative gradient smaller than 0.7 times the adiabatic gradient, where the 0.7 factor corresponds to the ≈30% error in the fit of the Rosseland mean opacity. This reflects the sensitivity of the deep atmospheric structure to the Rosseland mean opacities calculations. As shown by Freedman et al. (2014), a variation of ≈10–50% is common between different calculations of the Rosseland mean opacities.
Several trends concerning the depth of the radiative/convective boundary are apparent in Fig. 14. For a given effective temperature, the radiative/convective boundary is much deeper for irradiated planets than for nonirradiated planet (as shown by e.g., Guillot & Showman 2002). Indeed the pressure of the radiative/convective boundary is almost constant while T_{eff}>T_{int} and decreases suddenly when T_{eff} ≈ T_{int} and thus T_{μ∗} ≲ T_{int}. The value of P_{R/C} also increases with gravity, which is directly linked to the decrease in the radiative gradient with gravity. Finally, we see that our analytical radiative/convective model properly predicts the depth of the radiative/convective boundary and its variation with gravity, internal temperature and effective temperature.
Fig. 14 Pressure of the radiative/convective boundary as a function of the effective temperature of the planet for different gravities and internal flux obtained with the numerical model (dashed lines) and with the analytical one (plain lines). 

Open with DEXTER 
5.5. Low irradiation planets and brown dwarfs
Gravitational contraction and deuterium burning can be a significant source of internal luminosity in young giant planets and brown dwarfs, respectively. This luminosity can overtake the stellar irradiation as the dominant heating source in the atmosphere. We calculated temperature/pressure profiles for planets with an internal temperature of 300 K and 1000 K and derived the same analytical models as in the case with T_{int} = 100 K presented in Sect. 5.2. Figs. 15 and 16 show that at pressures of 1 − 100 bar model D of Sect. 5.2 – derived considering an internal temperature of 100 K – correctly matches the numerical temperature/pressure profile for higher internal temperatures and can therefore be used as a boundary condition for internal structure model. In the T_{int} = 1000 K low gravity case we observe a discrepancy of ≈20% between the numerical and the analytical model. This discrepancy is due to a the ionization of hydrogen not taken into account in our fit of the convective gradient. The upper atmosphere is well represented as long as T_{int} ≲ T_{μ∗} (i.e., T_{eff} ≈ T_{μ∗}). When the internal luminosity becomes the dominant energy source of the atmosphere, our analytical model becomes systematically hotter than the numerical model at low pressures with a discrepancy up to 40% between the two models (see also Fig. A.3).
Fig. 15 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for T_{int} = 300 K. For a full description see Fig. 12. 

Open with DEXTER 
Fig. 16 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for T_{int} = 1000 K. For a full description see Fig. 12. 

Open with DEXTER 
5.6. Application to the atmospheres of solarsystem giant planets
In Fig. 17 we use model D of Sect. 5.2 without further adjustments to calculate the temperaturepressure profiles of the four giant planets of our solarsystem. Our model pertains to cloudless, solar composition atmospheres, whereas Jupiter, Saturn, Uranus and Neptune have clouds and are all found to be significantly enriched in heavy elements. In spite of this, the analytical temperature profiles match the observations well (within 5 to 10%) for pressures higher than about 0.05 bar for Jupiter and Saturn and pressures higher than 0.1 bar for Uranus and Neptune. Specifically, the temperature at 1 bar is found to be of 163 K, 131 K, 85 K and 77 K for the four planets. The corresponding temperatures inferred from observations are 165 ± 5 K, 135 ± 5 K, 76 ± 2 K, and 72 ± 2 K, respectively (e.g., Lindal 1992; Guillot 2005). Our model can therefore be used as a boundary condition for internal structure models of the solarsystem giant planets. We predict that the radiative zone of Uranus extends down to P ≈ 100 bar, what might affect its cooling history and outward thermal flux. Since our model cannot take into account the presence of a detached convective zone, it predicts an unphysical, higherthanadiabatic gradient between 0.5 bar and 7 bar. For an internal flux equal to its observed 1 − σ upper value, the convective zone extends continuously from the deep layers to the ≈0.5 bar level, leading to a more realistic temperature profile in Uranus’ deep interior. A similar behavior is recovered when the opacities are increased by 50%.
At low pressures, our model is systematically cooler than the observations, revealing the limits of our approach. In those cold planets, the photon deposition layer spans several orders of magnitudes in pressure and is poorly modeled by considering only three visible opacities. Moreover, our coefficients are constant with gravity whereas it is apparent from Fig. 11 that the dependance with gravity is stronger at low effective temperature.
Fig. 17 Pressure temperature profiles of the four giant planets of the solarsystem as derived from observations (dotted line, Lindal 1992) and from model D of Sect. 5.2 (plain lines). The thin lines are radiative regions whereas the thick lines are convective regions. For Uranus we plot two profiles: one with the measured value of the internal temperature T_{int} = 29.4 K (the hotter profile) and one with the 1σ upper limit T_{int} = 35.5 K (the cooler profile). The fiducial model of Uranus have a deep radiative zone, it illustrates the difficulties of the method when a detached convective zone appears (here between 0.5 bar and 7 bar), since it is not taken into account by the analytical solution. 

Open with DEXTER 
5.7. Recommended model
When modeling gas giant planets of solar composition, we recommend the use of model D of Sect. 5.2. This model uses the solution of the radiative transfer equation given by Paper I where the first five parameters describing the opacities are expressed as a function of the effective temperature (see Table 1) whereas the analytical Rosseland mean opacities are given by Valencia et al. (2013). Model D is fully analytical, yet achieves an overall accuracy of 10% in temperature (at a given pressure) for irradiated giant planet atmospheres of solar composition with gravities in the range 2.5 − 250 m/s^{2}, internal temperatures of 100 to 1000 K and effective temperatures from 100 to 3000 K. When the internal flux dominates over the external flux (i.e., T_{int}>T_{μ∗}), model D becomes less accurate in the medium and upper atmosphere with an error that can reach ≈20% and ≈40%, respectively. In the deep atmospheres (P ≈ 1 − 100 bar) it remains accurate and can be used as boundary condition for internal structure models.
This accuracy is to be compared to that of simpler models. For example, models where the temperature is set to the effective temperature at τ = 2/3 and the profile is assumed to follow the diffusion approximation below (i.e., the socalled Eddington boundary condition). We calculated that this commonly used prescription (e.g., Bodenheimer et al. 2003; Batygin et al. 2011,among many others ) lead to an error in the temperature profile below the τ = 2/3 level on the order of ≈30% except fortuitously for 800 K <T_{eff}< 1200 K where the error is lower than 10%. Such an error on the boundary condition of interior models can strongly affect internal structure and planetary evolution calculations. Semigrey models (e.g., Hansen 2008; Guillot 2010) cannot reach an accuracy better than 20%, even with adjusted variable opacity coefficients.
A FORTRAN implementation of model D is available for download on the internet^{3}.
6. The role of TiO and VO
The most irradiated planets have dayside atmospheric temperatures high enough such that, for a solar composition atmosphere some metal oxides such as titanium and vanadium oxides (TiO and VO, respectively), are chemically stable in gas phase (Lodders 2002). To this date, there is no firm direct detection of TiO at the atmospheric limb of an exoplanet (see Désert et al. 2008; Huitson et al. 2013; Sing et al. 2013). Condensation in a vertical cold trap (Spiegel et al. 2009), in an horizontal cold trap (Parmentier et al. 2013), dissociation by stellar radiation (Knutson et al. 2010) or the presence of a high C/O ratio (Madhusudhan 2012) have been proposed as mechanisms to deplete TiO and VO from the upper atmosphere of irradiated planets (see also the review in Parmentier et al. 2014).
Several studies (e.g., Hubeny et al. 2003; Fortney et al. 2008) show that, if present in solar abundance in the upper atmosphere of irradiated planets, titanium and vanadium oxides should change significantly the temperature structure of those atmospheres, creating a strong thermal inversion (also called stratosphere) at low pressures and reducing the temperature in the deep atmosphere. The signature of such a stratosphere have been searched for in the secondary eclipse spectra of a dozen of planets. Although no unambiguous evidence for a stratosphere has yet been found, their presence have not been ruled out either (Hansen et al. 2014). Here, we examine how TiO and VO opacities shape the atmospheric temperature profile of irradiated planets.
Fig. 18 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for an atmosphere without TiO/VO. For a full description see Fig. 12. 

Open with DEXTER 
6.1. Retrieved coefficients
We calculated a grid of pressure/temperature profiles and derived the same analytical models as in the previous section in the case where TiO and VO have been removed from the whole atmosphere by any of the aforementioned processes. The resulting analytical model provides a good fit of the numerical profile, as seen in Fig. 18 and quantified in Fig A.2. The corresponding coefficients are presented in Table A.1 and compared with the solar composition case in Fig. 19.
As expected, the absence of TiO and VO changes significantly how the atmosphere absorbs the stellar irradiation: γ_{v1} and γ_{v2} are 10 times smaller than in the solar composition case, showing that the first two thirds of the stellar energy are deposited at deeper levels than when TiO/VO are present. Conversely, γ_{v3} remains almost constant: because the spectral width of the TiO band is smaller than the spectral width of the incoming stellar flux (see Fig A.1), part of the stellar flux remains unaffected by the presence of TiO and is absorbed at the same depth in the two cases.
Usually not considered as important, the thermal coefficients are also affected by the presence of TiO and VO. β is ≈20% smaller and γ_{p} is ≈10 times higher than in the solarcomposition case. Although TiO affects principally the optical region of the spectrum, hot Jupiters are hot enough that the local Planck function of their atmospheres can overlap with the TiO band. The broadband opacities of TiO have two main effects: first, it fills the gap in the opacities seen between the Rayleigh scattering dominated opacities at small wavelength and the water bands beginning at ≈1 μm, leading to a grayer atmosphere and thus a smaller value of γ_{p}. Second, it changes the shape of the opacities: in the case without TiO, the opacities in the optical spectral range are dominated by sodium and potassium lines, leading to smaller values of β whereas in the solar composition case the optical opacities are dominated by the wide TiO band, leading to larger values of β.
Fig. 19 Coefficients from our analytical model corresponding to the solar composition case (dashed lines, from Table 1) and the case where TiO and VO have been removed from the atmosphere (plain lines, Table A.1). 

Open with DEXTER 
6.2. Effects of TiO and VO on the temperature profile
We show that both the visible and the thermal coefficients are affected by the presence of TiO in the atmosphere. Thus both the greenhouse effect, the upper atmospheric cooling and the blanketing effect should contribute to the difference in the temperature profile between a solarcomposition atmosphere and an atmosphere without TiO/VO. As shown by the blue and the red curves of Fig. 20, the disappearance of TiO and VO from the atmosphere reduces the temperature in the upper atmosphere (removing the temperature inversion) and increases the temperature in the deep atmosphere.
We now wish to disentangle the contribution of the visible and the thermal effects in shaping the thermal profile between these two cases. We calculate an intermediate temperature profile in which the stellar flux is absorbed at the same Rosseland optical depth than in the solarcomposition case but where TiO and VO have been removed from the thermal opacities. Thus, the grey model in Fig. 20 have the same greenhouse effect than the solarcomposition profile but have the same thermal effects than in the case without TiO/VO. This grey profile lay approximately at the midpoint between the profile with TiO/VO and the profile without TiO/VO. Nongrey thermal effects are therefore as important as nongrey absorption of the stellar light to set the difference between the cases with and without TiO/VO.
In the upper atmosphere, our intermediate profile does not show any temperature inversion although it absorbs the stellar irradiation with the same strength than in the solarcomposition case. The visible effects of TiO/VO alone cannot explain the presence of a temperature inversion in hotJupiters atmospheres. The inversion mainly happens because the broadband absorption of TiO increases the grayness of the opacities, reducing the ability of the upper atmosphere to cool down efficiently.
In the deep atmosphere, the absence of TiO and VO yields both a stronger greenhouse effect and a stronger blanketing effect. Both effects contribute equally to increase the deep temperature profile, leading to higher interior temperatures for the same global effective temperature and therefore to a slower evolution (Parmentier & Guillot 2011; Budaj et al. 2012; Spiegel & Burrows 2013). In the absence of TiO and VO, hot Jupiters should be on average larger than if these gases are present.
7. Conclusion
Analytical solutions of the radiative transfer equation, although derived using very restrictive (but necessary) approximations, offer a deep insight in the physical processes shaping the temperature profile of planetary atmospheres and can provide fast and roughly accurate solutions to be incorporated in more complex planetary models.
In this study we used handinhand the analytical model derived in Paper I, that includes nongrey effects due to the visible and thermal opacities, and a stateoftheart numerical model that solves the radiative transfer equation considering their full frequency and angular dependency.
We first quantified the validity of the Eddington approximation. We showed that this approximation leads to errors in the temperature profile of at most 2% in the grey case and 4% in the nongrey case.
Planets with a thick atmosphere usually become convective below a certain depth. Thus, a common way to produce a radiative/convective temperature profile is to switch from a radiative solution to a convective solution whenever the Schwarzchild criterion is met, considering that the radiative solution remains unaffected by the presence of a convective zone below it. We showed that this approach is always valid in the grey case – the error due to the Eddington approximation being frozen at the radiative/convective boundary and propagated along the convective zone. However, for nongrey atmospheric opacities, we showed that this method is valid only as long as the radiative/convective boundary remains in the optically thick layer of the atmosphere. When the radiative/convective boundary is in the optically thin region of the atmosphere, the radiative solution is very sensitive to the precise location of the radiative/convective boundary and this common approach can lead to relative errors of tens of percentage points when estimating the upper, radiative, atmospheric temperatures.
We showed that semigrey effects (greenhouse and antigreenhouse effects) are not sufficient to explain the atmospheric temperature profiles calculated with the full frequency dependent opacities and that nongrey thermal effects need to be taken into account. We provided a reliable method to obtain the visible coefficients of our analytical model directly from the opacities and explored how the thermal coefficients could also be directly derived from the knowledge of the linebyline atmospheric opacities.
Fig. 20 Temperaturepressure profile for T_{eff} = 2268 K, g = 25 m/s^{2} and T_{int} = 100 K for a solarcomposition atmosphere (red), an atmosphere where TiO and VO have been removed from the atmosphere (blue) and and atmosphere where the stellar irradiation is absorbed at the same optical depth as in the case with TiO/VO but where TiO/VO have been removed from the thermal opacities (grey). Our analytical model is shown as plain lines whereas the numerical solution is shown as dashed lines. 

Open with DEXTER 
In particular, we showed that the presence of TiO can warm up the upper atmosphere and cool down the deep atmosphere not only because it absorbs a significant amount of stellar irradiation in the upper atmosphere, but also because its broad band opacity reduces the nongrey thermal blanketing effect. The presence of a thermal inversion in hotJupiters induced by an extra absorber in the optical opacities cannot be explained solely by an increased absorption of the stellar light and depends mostly on the ability of the atmosphere to cool down via the nongrey thermal effects.
Finally, using an a priori determination of the visible coefficients and an a posteriori determination of the thermal coefficients, we provide a fully analytical model for solar composition optically thick atmospheres. This model agrees with the numerical calculations within 10% over a wide range of gravities and effective temperatures. Our model leads to a much better estimate of the deep temperature profile than the previous analytical estimates. Therefore, when modeling the atmospheric structure of giant planets, we recommend the use of Model D described in Sect. 5.7 that uses the analytical expressions derived in Paper I with the five parameters given in Table 1 for a solarcomposition atmosphere and by Table A.1 in the case where TiO and VO have been removed from the atmosphere. The Rosseland mean opacities of model D are given by the functional form of Valencia et al. (2013) and the convective gradient follows Eq. (13). For convenience, we provide an implementation in FORTRAN of our model in the CDS and at the address http://www.oca.eu/parmentier/nongrey.
A factor was missing in the expression of T_{irr} in Sect. 3.7 of Paper I.
Acknowledgments
We acknowledge the anonymous referee for his general comments that increased the quality of this manuscript. This work was performed in part thanks to a joint Fulbright Fellowship to V.P. and T.G. The whole project would not have been possible without the help and support of Douglas Lin. We also acknowledge University of California Santa Cruz for hosting us while this work was carried out. During the last part of this work, V.P. was under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.
References
 Abramowitz, M., & Stegun, I. A. 1965, Handbook of mathematical functions with formulas, graphs, and mathematical tables, eds. M. Abramowitz, & I. A. Stegun (New York: Dover) [Google Scholar]
 Batygin, K., Stevenson, D. J., & Bodenheimer, P. H. 2011, ApJ, 738, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P., Laughlin, G., & Lin, D. N. C. 2003, ApJ, 592, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Budaj, J., Hubeny, I., & Burrows, A. 2012, A&A, 537, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856 [NASA ADS] [CrossRef] [Google Scholar]
 Cahoy, K. L., Marley, M. S., & Fortney, J. J. 2010, ApJ, 724, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative transfer (New York: Dover) [Google Scholar]
 Désert, J.M., VidalMadjar, A., Lecavelier Des Etangs, A., et al. 2008, A&A, 492, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., LustigYaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Goody, R. M., & Yung, Y. L. 1989, Atmospheric radiation: theoretical basis [Google Scholar]
 Guillot, T. 2005, Ann. Rev. Earth Planet. Sci., 33, 493 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Guillot, T. 2010, A&A, 520, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillot, T., & Showman, A. P. 2002, A&A, 385, 156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillot, T., Burrows, A., Hubbard, W. B., Lunine, J. I., & Saumon, D. 1996, ApJ, 459, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, B. M. S. 2008, ApJS, 179, 484 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, C. J., Schwartz, J. C., & Cowan, N. B. 2014, MNRAS, 444, 3632 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, MNRAS, 420, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011 [NASA ADS] [CrossRef] [Google Scholar]
 Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252 [NASA ADS] [CrossRef] [Google Scholar]
 Knutson, H. A., Howard, A. W., & Isaacson, H. 2010, ApJ, 720, 1569 [NASA ADS] [CrossRef] [Google Scholar]
 Lindal, G. F. 1992, AJ, 103, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Lodders, K. 2002, ApJ, 577, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N. 2012, ApJ, 758, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Magalhaes, J. A., Seiff, A., & Young, R. E. 2002, Icarus, 158, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., Ackerman, A. S., Cuzzi, J. N., & Kitzmann, D. 2013, Comparative Climatology of Terrestrial Planets, eds. S. J. Mackwell, A. A. SimonMiller, J. W. Harder, & M. A. Bullock (Tucson: University Arizona Press), 610, 367 [Google Scholar]
 McKay, C. P., Pollack, J. B., & Courtin, R. 1989, Icarus, 80, 23 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Mihalas, D. 1978, Stellar atmospheres. 2nd edn. (San Francisco: W. H. Freeman and Co) [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Parmentier, V., & Guillot, T. 2011, in EPSCDPS Joint Meeting 2011, 1367 [Google Scholar]
 Parmentier, V., & Guillot, T. 2014, A&A, 562, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parmentier, V., Showman, A. P., & de Wit, J. 2014, Exp. Astron., in press [arXiv:1401.3673] [Google Scholar]
 Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press) [Google Scholar]
 Rauscher, E., & Menou, K. 2012, ApJ, 750, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Robinson, T. D., & Catling, D. C. 2012, ApJ, 757, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Sing, D. K., Lecavelier des Etangs, A., Fortney, J. J., et al. 2013, MNRAS [Google Scholar]
 Spiegel, D. S., & Burrows, A. 2013, ApJ, 772, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Strom, S. E., & Kurucz, R. L. 1966, J. Quant. Spectr. Rad. Transf., 6, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Sudarsky, D., Burrows, A., & Pinto, P. 2000, ApJ, 538, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Toon, O. B., McKay, C. P., Ackerman, T. P., & Santhanam, K. 1989, J. Geophys. Res., 94, 16287 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10 [NASA ADS] [CrossRef] [Google Scholar]
 West, R. A., Baines, K. H., Friedson, A. J., et al. 2004, Jovian clouds and haze, eds. F. Bagenal, T. E. Dowling, & W. B. McKinnon (Cambridge: Cambridge University Press), 79 [Google Scholar]
Online material
Appendix A: Additional material
The opacities in the form of kcoefficients used in the numerical model and discussed in Sect. 5.2 are presented in Fig. A.1.
The analytical model adjusted to match the temperature/pressure profile of an atmosphere without TiO, discussed in Sect. 6 is presented in Fig. A.2 and Table A.1.
The effect of a strong internal luminosity on the analytical model, discussed in Sect. 5.5 is presented if Fig. A.3.
Fig. A.1 Opacities from Freedman et al. (2008) organized as kcoefficient inside each bin of wavelength atmospheres with different irradiation. The first five panels are for solarcomposition atmospheres whereas the bottom right panel is for an atmospheres depleted in TiO/VO. The different colors are for different temperature and pressure taken along the corresponding numerical PT profile. The thick bars on top represents the wavelength range where 90% of the thermal flux is emitted, the thin bars where 99% of the thermal flux is emitted. We used T_{int} = 100K, and g = 25m/s^{2}. 

Open with DEXTER 
Fig. A.2 Top panel: coefficients obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets with different gravities and an internal temperature of 100 K. TiO and VO have been artificially removed from the atmosphere. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4}bar <P< 10^{2}bar, the second line for 10^{2}bar <P< 10^{0}bar and the third line for 10^{0}bar <P< 10^{2}bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/ g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10 (25 ms^{2}/ g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/ g). 

Open with DEXTER 
Fig. A.3 Top panel: coefficients obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets with a solar composition atmosphere with different gravities and an internal temperature of 100 K, 300 K and 1000 K. The model of Col. D uses the functional form of the coefficients derived in the case T_{int} = 100 K only. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4}bar <P< 10^{2}bar, the second line for 10^{2}bar <P< 10^{0}bar and the third line for 10^{0}bar <P< 10^{2}bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/ g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10 (25 ms^{2}/ g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/ g) When T_{eff} ≈ T_{int} the low pressures are not properly represented by our model (see Sect. 5.5 for more details). 

Open with DEXTER 
All Tables
Functional form of the planet planeparallel albedo log _{10}(A_{μ∗}) = a + bX where X = log _{10}(T_{eff0}) and g is the gravity of the planet in m/s^{2}.
All Figures
Fig. 1 Linebyline opacities as a function of wavelength for five different conditions corresponding to different points in the PT profile of a giant planet with g = 25 m/s^{2}, , T_{int} = 100 K and T_{μ∗} = 1231 K, corresponding to the dayside average profile of a planet orbiting at 0.053 AU from a sunlike star. Inside each bin of frequency, we plot the cumulative distribution function of the opacities instead of the linebyline opacity function. The purple bar represents the wavelength range where 90% (thick line) and 99% (thin line) of the stellar energy is emitted. The other horizontal bars represents the wavelength range of the thermal emission of the planet at different locations along the PT profile. 

Open with DEXTER  
In the text 
Fig. 2 Radiative numerical temperature profile in units of effective temperature as a function of optical depth compared to the analytical solution from Chandrasekhar (1960) using the discrete ordinate in the first, second, third, fourth and fifth approximation. The left panel shows the profiles whereas the right panel shows their relative difference (T_{a}/T_{n} − 1 where T_{a} is the analytical solution and T_{n} is the numerical solution). 

Open with DEXTER  
In the text 
Fig. 3 Comparison between the radiative numerical solution (black line), our work (red line) and Guillot (2010) model for two different values of f_{H} (blue and green lines) for a fully radiative semigrey atmosphere with γ_{v} = 0.25 (plain lines) or γ_{v} = 10 (dashed lines). We set , T_{irr} = 1288 K and T_{int} = 500 K. 

Open with DEXTER  
In the text 
Fig. 4 Comparison between the analytical model (plain lines) and the radiative numerical model (dashed lines) for different values of κ_{2}/κ_{1} (left panel). The right panel shows the relative difference between the analytical and the numerical solution for each case. We used , T_{irr} = 1288K, T_{int} = 500K, γ_{v} = 0.25 and β = 0.86. 

Open with DEXTER  
In the text 
Fig. 5 Comparison of our numerical and analytical radiativeconvective models for different adiabatic gradients in the nonirradiated, grey case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{int} = 500 K and g = 25 m/s^{2}. Note that the cases ∇_{ad} × 2 (green) and ∇_{ad} × 4 (red) are superimposed. 

Open with DEXTER  
In the text 
Fig. 6 Comparison of our numerical and analytical radiativeconvective solutions for different adiabatic gradients in the nonirradiated, nongrey case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{int} = 500 K, g = 25 m/s^{2}, κ_{2}/κ_{1} = 10^{2}, and β = 0.83. The cases ∇_{ad} × 2 (green) and ∇_{ad} × 4 (red) are superimposed. 

Open with DEXTER  
In the text 
Fig. 7 Comparison of our numerical and analytical radiative/convective solutions for different adiabatic gradients in the nongrey, irradiated case. The thin line is the radiative zone and the thick one represents the convective zone. We used T_{irr} = 1288 K, T_{int} = 500 K, γ_{v} = 0.25, , κ_{2}/κ_{1} = 10^{2} and β = 0.83. 

Open with DEXTER  
In the text 
Fig. 8 Absorbed stellar flux from the numerical model (dots) and from the analytical model (lines) considering respectively 1, 2, 3 or 4 absorption bands in the visible for our fiducial model atmosphere (see Fig. 1). 

Open with DEXTER  
In the text 
Fig. 9 Pressuretemperature profiles calculated using the numerical model and the full set of opacities (black), the semigrey (blue) and the nongrey (red) analytical models. As un Fig. 1, g = 25 m/s^{2}, , T_{int} = 100 K and T_{μ∗} = 1253 K. The coefficients used for the analytical models are taken from Table 1 with γ_{p} = 1 for the semigrey case. The nongrey model is a much better match to the numerical profile than the semigrey one. 

Open with DEXTER  
In the text 
Fig. 10 Cumulative distribution function of the opacities at P = 0.85 bar and T = 1464 K, corresponding to the τ = 2/3 level of an atmosphere with T_{μ∗} = 1253 K and a gravity g = 25 m/s^{2}. The Yaxis represents the fraction of frequency where the monochromatic opacities are lower than the corresponding κ_{0} of the Xaxis. The red line shows the value of the Rosseland mean opacity and the black line the Planck mean opacity. 25% of the frequency range have monochromatic opacities smaller than the Rosseland mean opacity whereas 90% have monochromatic opacities smaller than the Planck mean opacity. 

Open with DEXTER  
In the text 
Fig. 11 Top panel: coefficients γ_{v3}, γ_{v2}, γ_{v1}, β, and γ_{p} obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets of solar composition with different gravities and an internal temperature of 100 K. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4} bar <P< 10^{2} bar, the second one for 10^{2} bar <P < 10^{0} bar and the third one for 10^{0}bar <P< 10^{2} bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/g) ≲ τ_{R} ≲ 10 (25 ms^{2}/g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/g). 

Open with DEXTER  
In the text 
Fig. 12 Comparison between the numerical solutions (dashed lines) and the analytical solutions of model D (using the functional form of the coefficients given in Table 1) over a wide range of irradiation temperatures for a giant planet of solar composition orbiting a sunlike star. Here we used g = 25 m/s^{2}, T_{int} = 100 K, and . Thick and thin lines represent convective and radiative zones, respectively. 

Open with DEXTER  
In the text 
Fig. 13 Planeparallel albedos for solarcomposition, clearsky atmospheres with different values of gravity and internal temperature, numerical model (dots) and analytical fit described in Table 2 (lines). In the inset, the albedo plotted in linear scale. Here we approximate the Bond albedo of the planet by the planeparallel albedo. 

Open with DEXTER  
In the text 
Fig. 14 Pressure of the radiative/convective boundary as a function of the effective temperature of the planet for different gravities and internal flux obtained with the numerical model (dashed lines) and with the analytical one (plain lines). 

Open with DEXTER  
In the text 
Fig. 15 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for T_{int} = 300 K. For a full description see Fig. 12. 

Open with DEXTER  
In the text 
Fig. 16 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for T_{int} = 1000 K. For a full description see Fig. 12. 

Open with DEXTER  
In the text 
Fig. 17 Pressure temperature profiles of the four giant planets of the solarsystem as derived from observations (dotted line, Lindal 1992) and from model D of Sect. 5.2 (plain lines). The thin lines are radiative regions whereas the thick lines are convective regions. For Uranus we plot two profiles: one with the measured value of the internal temperature T_{int} = 29.4 K (the hotter profile) and one with the 1σ upper limit T_{int} = 35.5 K (the cooler profile). The fiducial model of Uranus have a deep radiative zone, it illustrates the difficulties of the method when a detached convective zone appears (here between 0.5 bar and 7 bar), since it is not taken into account by the analytical solution. 

Open with DEXTER  
In the text 
Fig. 18 Comparison of our numerical (dashed lines) and analytical (plain lines) solutions for an atmosphere without TiO/VO. For a full description see Fig. 12. 

Open with DEXTER  
In the text 
Fig. 19 Coefficients from our analytical model corresponding to the solar composition case (dashed lines, from Table 1) and the case where TiO and VO have been removed from the atmosphere (plain lines, Table A.1). 

Open with DEXTER  
In the text 
Fig. 20 Temperaturepressure profile for T_{eff} = 2268 K, g = 25 m/s^{2} and T_{int} = 100 K for a solarcomposition atmosphere (red), an atmosphere where TiO and VO have been removed from the atmosphere (blue) and and atmosphere where the stellar irradiation is absorbed at the same optical depth as in the case with TiO/VO but where TiO/VO have been removed from the thermal opacities (grey). Our analytical model is shown as plain lines whereas the numerical solution is shown as dashed lines. 

Open with DEXTER  
In the text 
Fig. A.1 Opacities from Freedman et al. (2008) organized as kcoefficient inside each bin of wavelength atmospheres with different irradiation. The first five panels are for solarcomposition atmospheres whereas the bottom right panel is for an atmospheres depleted in TiO/VO. The different colors are for different temperature and pressure taken along the corresponding numerical PT profile. The thick bars on top represents the wavelength range where 90% of the thermal flux is emitted, the thin bars where 99% of the thermal flux is emitted. We used T_{int} = 100K, and g = 25m/s^{2}. 

Open with DEXTER  
In the text 
Fig. A.2 Top panel: coefficients obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets with different gravities and an internal temperature of 100 K. TiO and VO have been artificially removed from the atmosphere. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4}bar <P< 10^{2}bar, the second line for 10^{2}bar <P< 10^{0}bar and the third line for 10^{0}bar <P< 10^{2}bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/ g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10 (25 ms^{2}/ g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/ g). 

Open with DEXTER  
In the text 
Fig. A.3 Top panel: coefficients obtained for the six different models described in Sect. 5.2 as a function of the irradiation temperature for planets with a solar composition atmosphere with different gravities and an internal temperature of 100 K, 300 K and 1000 K. The model of Col. D uses the functional form of the coefficients derived in the case T_{int} = 100 K only. Bottom panel: mean relative difference between the numerical and the analytical model for the six different models described in Sect. 5.2. The first line is the mean difference for 10^{4}bar <P< 10^{2}bar, the second line for 10^{2}bar <P< 10^{0}bar and the third line for 10^{0}bar <P< 10^{2}bar. In terms of Rosseland optical depth, the low pressure zone corresponds to the optically thin part of the atmosphere with 10^{8} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{2} (25 ms^{2}/ g). The medium pressure zone corresponds to the transition from optically thin to optically thick with 10^{4} (25 ms^{2}/ g) ≲ τ_{R} ≲ 10 (25 ms^{2}/ g). The high pressure zone corresponds to the optically thick part of the atmosphere with (25 ms^{2}/ g) ≲ τ_{R} ≲ 10^{4} (25 ms^{2}/ g) When T_{eff} ≈ T_{int} the low pressures are not properly represented by our model (see Sect. 5.5 for more details). 

Open with DEXTER  
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.