© 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 non-grey 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 non-grey 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 line-by-line 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, solar-composition, semi-infinite 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 state-of-the-art 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 state-of-the 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 non-grey 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 $\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{e}{\mathrm{q}}_{\mathrm{0}}}^{\mathrm{4}}$ where the equilibrium temperature for zero albedo is defined as ${\mathit{T}}_{\mathrm{e}{\mathrm{q}}_{\mathrm{0}}}^{\mathrm{4}}\mathrm{\equiv }\frac{{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}}{\mathrm{4}}{\left(\frac{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{a}}\right)}^{\mathrm{2}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(1)At a given point in the planet, the incoming stellar flux is characterized by a temperature ${\mathit{T}}_{\mathrm{ir}{\mathrm{r}}_{\mathrm{0}}}^{\mathrm{4}}\mathrm{=}\mathrm{4}{\mathit{T}}_{\mathrm{e}{\mathrm{q}}_{\mathrm{0}}}^{\mathrm{4}}\mathit{.}$(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: ${\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }\mathrm{0}}}^{\mathrm{4}}\mathrm{=}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{ir}{\mathrm{r}}_{\mathrm{0}}}^{\mathrm{4}}\mathit{.}$(3)Where μ_∗ = cosθ_∗. A part 1 − A_μ∗ of this flux is reflected back to space, where A_μ∗ is the angle-dependent reflectivity. We characterize the flux that penetrates the atmosphere by the temperature T_μ∗: ${\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{=}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathit{.}$(4)where ${\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{{\mathit{\mu }}_{\mathrm{\ast }}}\mathrm{)}{\mathit{T}}_{\mathrm{ir}{\mathrm{r}}_{\mathrm{0}}}^{\mathrm{4}}$.

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 ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ 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 so-called isotropic approximation (e.g., Guillot 2010): $\{\begin{array}{c}\\ {\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}\\ {\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{B}}\mathrm{)}\mathrm{4}\mathit{f}{\mathit{T}}_{\mathrm{e}{\mathrm{q}}_{\mathrm{0}}}^{\mathrm{4}}\end{array}\mathit{,}$(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 ${\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}$ and ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ in Eq. (76) of Paper I¹.

The atmosphere is also heated by the planet interior and receives a flux $\mathit{\sigma }{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}$ from below. The total energy budget of the atmosphere is characterized by its effective temperature ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{,}$(6)which is valid for both the averaged and the non-averaged cases. Finally, we define an effective temperature for zero albedo: ${\mathit{T}}_{\mathrm{e}{\mathrm{ff}}_{\mathrm{0}}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{{\mathit{\mu }}_{\mathrm{\ast }\mathrm{0}}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{.}$(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 solar-composition opacities provided by Freedman et al. (2008). These opacities have been calculated for a solar-composition 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 so-called correlated-k 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 1-Jupiter mass and radius planet at 0.05 AU from a Sun-like 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 comb-like structure. Deeper-down, the wavelength dependence is mostly due to the presence of molecular bands and takes place on scales significantly larger than our bin size.

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 low-pressures 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 picket-fence 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 line-by-line 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 non-greyness 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: κ₁/κ₂, where κ₁ is the highest of the two opacities and κ₂ the lowest. Value of κ₁/κ₂ between 10⁴ − 10⁵ in the upper atmosphere and between 10 − 100 in the deep atmosphere can be estimated from Fig. 1 for a typical hot-Jupiter. The simple relationship between κ₁/κ₂ 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². 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 non-grey 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, line-by-line, frequency, pressure-and temperature-dependent 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 delta-discrete ordinates method of Toon et al. (1989) for the incident stellar radiation and the two-stream 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 correlated-k 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 k-coefficients. 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 ${\mathrm{\int }}_{-1}^{\mathrm{1}}\mathit{\mu }{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{\omega }}_{\mathit{i}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathrm{\right)}$(8)with the ω_i and the μ_i being tabulated in Abramowitz & Stegun (1965). Here we use 5 Gauss points. The EGP model calculates a self-consistent 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 k-coefficient method can be used to easily implement the simplified opacities of Parmentier & Guillot (2014) by setting a given number of k-coefficients at κ₁ and the other ones at κ₂ 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.

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 non-irradiated 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, non-irradiated 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 ${{}^{\mathrm{\int }}}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }$, the numerical code uses the quadrature formula to compute the flux integral ${{}^{\mathrm{\int }}}_{-1}^{\mathrm{1}}\mathit{\mu }{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }$. 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, non-irradiated 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): ${\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}{\mathit{\mu }}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{3}}{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathit{.}$(9)This approximation is exact in two specific cases: when the radiation field is isotropic (I_ν(μ) = cte ∀μ) and when radiation field is semi-isotropic (I_ν(μ) = I⁺∀μ> 0 and I_ν(μ) = I⁻∀μ< 0). In the deep atmosphere, the radiation is quasi-isotropic 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 semi-isotropic. In between, the solution is only an approximation. In addition, a boundary condition relating two other moments of the intensity field must be adopted: ${\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathit{\mu }\hspace{0.17em}\mathrm{d}\mathit{\mu }\hspace{0.17em}{|}_{\mathrm{top}}\mathrm{\approx }{\mathit{f}}_{\mathrm{H}}{\mathrm{\int }}_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\hspace{0.17em}{|}_{\mathrm{top}}\mathit{.}$(10)These two conditions form what is called the Eddington approximation.

In the grey, non-irradiated case, those two approximations are linked and $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$. However using Eq. (9) and imposing $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ leads to the exact solution at the top of the atmosphere, even though it lacks of self-consistency. In the irradiated case and in the non-grey case the two approximations are independent and $\stackrel{\mathrm{h}}{.}$ is usually set to either 1/2 or $\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$, following the grey, non-irradiated 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, non-irradiated case. In order to estimate its magnitude in the grey and non-grey 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 κ₂/κ₁), β and the ratio T_irr/T_int.

3.1. Irradiated semi-grey solutions

In the semi-grey 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 $\stackrel{\mathrm{h}}{.}$ 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 ${\mathit{f}}_{\mathrm{H}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$) and the semi-grey limit of the model derived in Paper I with $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$ 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 semi-grey models: the ratio of the visible to the thermal opacities, γ_v.

Figure 3 compares these models for a typical irradiated Jupiter-mass exoplanet close to a solar-type 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 $\stackrel{\mathrm{h}}{.}$ – 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 non-irradiated 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 $\stackrel{\mathrm{h}}{.}$ and of the top boundary condition, but remains of the same order-of-magnitude as for the non-irradiated 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 ${\mathit{f}}_{\mathrm{H}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ 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 non-irradiated case, no choice of $\stackrel{\mathrm{h}}{.}$ can yield an exact skin temperature T(τ = 0) (see related discussion in Paper I).

3.2. Irradiated non-grey solutions

We now test the analytical model of Paper I in the non-grey 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 κ₂/κ₁ 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.

Figure 4 shows the resulting temperature-optical depth profiles and the relative difference between the numerical and analytical solutions. As κ₂/κ₁ 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 (κ₂/κ₁ = 1) corresponds to the semi-grey 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 non-greyness of the opacities. The maximum error (in absolute terms) increases from about ~2% to a little bit less than ~5% when κ₂/κ₁ is increased from 1 to 10⁵. 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³].

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 κ₁ and κ₂, 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 semi-isotropic 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 κ₁/κ₂ 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 two-peak 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 κ₂/κ₁ are uniform in the atmosphere), this is an acceptable level of uncertainty.

4. Consequences of convection on the overlaying radiative solution

At high-enough 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 collision-induced 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 low-irradiation 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 close-in 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 semi-grey models, Robinson & Catling (2012) recently derived a radiative-convective model that satisfies these two conditions, although it necessitates a small numerical integration. For non-grey thermal opacities, no analytical model solves self-consistently 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. Non-irradiated grey case

We first compare the solutions obtained in the non-irradiated 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.

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 re-calculating 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. Non-irradiated non-grey case

We now turn, with Fig. 6, to the non-irradiated non-grey case, using the fiducial values κ₂/κ₁ = 10² 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 non-grey 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 non-grey 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 (non-grey) 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.

4.3. Irradiated non-grey 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 κ₂/κ₁ 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 non-grey, non-irradiated case.

5. Modeling the non-grey 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 line-by-line 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 line-by-line 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 well-known 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 line-by-line opacities) to the mean Rosseland thermal opacity (using the local Planck function at the stellar effective temperature to weight the line-by-line 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 non-grey 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 non-grey visible opacities can relatively well approximate the absorbed stellar flux at all atmospheric levels. We therefore adopt the following method: using the line-by-line 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 $\mathit{F}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{0}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{\beta }}_{\mathrm{v}\mathit{i}}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}\mathit{i}}\mathit{\tau }\mathit{/}{\mathit{\mu }}_{\mathrm{\ast }}}\mathit{,}$(11)where the visible bands are homogeneously distributed in frequency (similar to the thermal bands), F₀ 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 line-by-line opacities and the semi-grey 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 semi-grey model, even though it models correctly the absorbed flux as a function of depth, lays far from the numerical solution. Clearly, non-grey thermal opacities are needed.

5.1.2. Thermal coefficients

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 non-grey 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 non-grey 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 non-grey 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 non-grey 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 κ₀ as a function of κ₀. 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 non-grey 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²), and an internal temperature of T_int = 100 K. We consider the case of a planet orbiting a sun-like star at various distances corresponding to irradiation temperatures from 100 K to 3000 K. All the profiles were calculated using ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$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 line-by-line 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 semi-grey 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 line-by-line 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 saddle-like 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 broad-band 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.

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 non-grey effects become more important in the upper atmosphere than in the deep atmosphere.

γ_pkeeps 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.

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 line-by-line 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 self-consistent fully analytical model.

Model E: here, the importance of non-grey 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 non-grey 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 non-grey absorption of the stellar irradiation therefore cannot by itself explain the temperature structure in planetary atmospheres. Non-grey 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 plane-parallel atmosphere irradiated with an angle ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$. 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 plane-parallel albedo calculated for ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$. This is similar to the so-called 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.

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 collision-induced 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 pressure-broaden 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 self-consistent 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 temperature-pressure profile can be determined by integrating the adiabatic gradient: ${\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{=}{\left(\frac{\mathit{\partial }\log \mathit{T}}{\mathit{\partial }\log \mathit{P}}\right)}_{\mathit{S}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(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: ${\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{\approx }\mathrm{0.32}\mathrm{-}\mathrm{0.1}\left(\frac{\mathit{T}}{\mathrm{3000}\hspace{0.17em}\mathrm{K}}\right)\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(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 non-irradiated 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.

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).

5.6. Application to the atmospheres of solar-system giant planets

In Fig. 17 we use model D of Sect. 5.2 without further adjustments to calculate the temperature-pressure profiles of the four giant planets of our solar-system. 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 solar-system 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, higher-than-adiabatic 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.

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², 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 so-called 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. Semi-grey 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³.

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.

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 solar-composition 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 broad-band 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 β.

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 solar-composition 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 solar-composition 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 solar-composition 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. Non-grey thermal effects are therefore as important as non-grey 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 solar-composition case. The visible effects of TiO/VO alone cannot explain the presence of a temperature inversion in hot-Jupiters 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 hand-in-hand the analytical model derived in Paper I, that includes non-grey effects due to the visible and thermal opacities, and a state-of-the-art 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 non-grey 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 non-grey 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 semi-grey effects (greenhouse and anti-greenhouse effects) are not sufficient to explain the atmospheric temperature profiles calculated with the full frequency dependent opacities and that non-grey 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 line-by-line atmospheric opacities.

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 non-grey thermal blanketing effect. The presence of a thermal inversion in hot-Jupiters 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 non-grey 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 solar-composition 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.

Online material

Appendix A: Additional material

The opacities in the form of k-coefficients 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.

---

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

All Tables

All Figures