© ESO, 2014

1. Introduction

The discovery of numerous star-planet systems and the possibility of characterizing the planets’ atmospheric properties has led to a great many publications using radiative transfer calculations, often taken “off-the-shelf” from numerical models. Given the infinite amount of possible compositions for mostly unknown exoplanetary atmospheres, it is highly valuable to be able to perform very fast calculations and also to understand what determines the thermal structure of an irradiated atmosphere.

Analytical radiative transfer solutions for atmospheres have been calculated with a variety of assumptions and in different contexts (e.g. Eddington 1916; Chandrasekhar 1935, 1960; King 1956; Matsui & Abe 1986; Weaver & Ramanathan 1995; Pujol & North 2003; Chevallier et al. 2007; Shaviv et al. 2011). However, the discovery of super-Earths, giant exoplanets, brown dwarfs, and low-mass stars close to a source of intense radiation has prompted the need for solutions that account for both an outside and an inside radiation field, and properly link low and high optical depths levels. Hubeny et al. (2003), Rutily et al. (2008), Hansen (2008), Guillot (2010), Robinson & Catling (2012), and Heng et al. (2012) provide these solutions in the framework of a semi-grey model, with one opacity for the incoming irradiation (generally mostly at visible wavelengths), and one opacity for the thermal radiation field (generally mostly at infrared wavelengths). These approximations have been used in hydrodynamical models of planetary atmospheres (e.g. Heng et al. 2011; Rauscher & Menou 2013), planetary evolution models (e.g. Miller-Ricci & Fortney 2010; Guillot & Havel 2011; Budaj et al. 2012), planet synthesis models (Mordasini et al. 2012a,b), retrieval methods (Line et al. 2012), and a variety of other applications.

Fig. 1 Optical depth vs. atmospheric temperature in units of the effective temperature. A numerical solution obtained from Fortney et al. (2008) (thick black line) is compared to the semi-grey analytical solutions of Guillot (2010) for values of the greenhouse factor ${\mathit{\gamma }}_{\mathrm{v}}^{-1}$ ranging from 0.01 to 100 (black to red lines). Low values of γv are redder. We used ${\mathit{\mu }}_{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$, Tirr = 1250K corresponding to the dayside average profile of a planet at 0.05  au from a sun-like star, and accounting for the albedo obtained in the numerical model. The internal temperature is Tint = 125  K and gravity 25ms-2. The effective temperature of the studied slice of atmosphere1 is defined by ${\mathit{T}}_{\mathrm{eff}\mathit{,}\hspace{0.17em}{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{+}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$. For the numerical solution, the relation between pressure and optical depth was calculated using Rosseland mean opacities; TiO and VO opacities were not included.

As shown in Fig. 1 for an atmosphere irradiated from above with a flux $\mathit{\sigma }{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$ and heated from below with a flux $\mathit{\sigma }{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}$, while semi-grey models provide solutions that are well-behaved when compared to full numerical solutions at optical depths larger than about unity, the temperatures at low-optical depths appear to be systematically hotter than in the numerical solutions. Most importantly, this occurs regardless of the choice of the two parameters of the problem, i.e. the thermal (infrared) opacity κ_th and the ratio of the visible to infrared opacity γ_v ≡ κ_v/κ_th. For hot Jupiters, as in the example in Fig. 1, the real temperature profiles at low optical depths can be several hundreds of Kelvins cooler than predicted by the semi-grey solutions.

The levels probed both by transit spectroscopy and by the observations of secondary eclipses of exoplanets often correspond to low-optical depth levels (e.g. Burrows et al. 2007; Fortney et al. 2008; Showman et al. 2009), i.e. where semi-grey models seem to systematically overestimate the temperatures. Furthermore, because the problem persists regardless of the main parameters, this implies that the functional form of the semi-grey solutions is probably not appropriate for inversion models. Non-grey effects are known to facilitate the cooling of the upper atmosphere (see Pierrehumbert 2010, for a qualitative explanation) so they must be included. This is the purpose of the present paper.

We first describe previous analytical methods used to solve the radiative transfer problem analytically. In Sect. 3, we then derive an analytical non-grey line model and compare it to previous models in Sect. 4. In Sect. 5 we study the role of non-grey effects in shaping the atmospheric thermal structure. Eventually we apply our model to the structure of irradiated giant planets in Sect. 6. We note here that while we focus the discussion on exoplanets, we believe that this model is applicable to a much wider variety of problems, as long as an atmosphere is irradiated both from above and below. Our method can also be used to solve the radiative transfer equations in other geometries, such as the thermal structure of protoplanetary disk. We provide our conclusions in Sect. 7.

2. Assumptions and previous analytical models

2.1. Setting

2.1.1. The equation of radiative transfer

Following Guillot (2010), we will consider the problem of a plane-parallel atmosphere in local thermodynamic equilibrium that receives from above a collimated flux $\mathit{\sigma }{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$ at an angle θ_∗ = cos^-1(μ_∗) from the vertical, and from below an isotropic flux $\mathit{\sigma }{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}$. The total energy budget of the modelled atmosphere is then set by ${\mathit{T}}_{\mathrm{eff}\mathit{,}\hspace{0.17em}{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{+}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$, which defines the effective temperature in this paper¹.The irradiation and intrinsic fluxes are generally characterized by very different wavelengths. Although this is not required in the solution that we propose, it is convenient to think of them as being emitted preferentially in the visible and in the infrared, respectively. Scattering processes can influence both the thermal and visible radiation. As shown by Heng et al. (2012), the solution including symmetrical scattering for the incoming radiation with the Eddington approximation is equivalent to the one obtained without scattering when the irradiation flux is reduced by a factor of (1 − A) and the visible opacity is reduced by a factor 1/$\sqrt{\mathit{\xi }}$, where A is the Bond albedo of the planet and ξ is the ratio of the absorption to the extinction opacity (see also Meador & Weaver 1980, for a review of the different two stream methods including scattering). In this paper, the irradiation flux and the visible opacity are treated as parameters, thus we implicitly take into account symmetrical scattering of the incoming radiation (e.g. Rayleigh scattering). Scattering of the thermal radiation, however, is neglected.

In order to solve the radiative transfer problem for a plane parallel atmosphere in local thermodynamic equilibrium, one has to solve the following equation for all frequencies ν and all directions μ (Chandrasekhar 1960), $\mathrm{-}\mathit{\mu }\frac{\mathrm{d}{\mathit{I}}_{\mathit{\mu \nu }}}{\mathrm{d}\mathit{m}}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{I}}_{\mathit{\nu \mu }}\mathrm{-}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{,}$(1)where I_μν is the specific intensity at the wavelength ν propagating with an angle θ = cos^-1(μ) with the vertical, κ_ν is the opacity at a given wavelength, B_ν is the Planck function, and dm = ρdz is the mass increment along the path of the radiation. As usual, T, ρ, and z are the atmospheric temperature, density and height, respectively. The main difficulty in solving Eq. (1)lies in its triple dependence on μ, ν, and T and its additional dependence on m. An analytical solution requires simplifications of the opacities and of the dependence of the radiation intensity on angle.

2.1.2. Opacities and optical depth

The need for simplification implies that mean opacities must be used. The most common one is the Rosseland mean, defined as $\frac{\mathrm{1}}{{\mathit{\kappa }}_{\mathrm{R}}}\mathrm{\equiv }{\left({\int }_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\nu }}}{\mathit{\partial T}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\right)}^{-1}{\int }_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{{\mathit{\kappa }}_{\mathit{\nu }}}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\nu }}}{\mathit{\partial T}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathit{.}$(2)When at all wavelengths the mean free path of photons is small compared to the scale height of the atmosphere, the radiative gradient obeys its well-defined diffusion limit and (unless convection sets in) the temperature gradient become that obtained from a grey atmosphere in which the opacity is set to the Rosseland mean (Mihalas & Mihalas 1984, p. 350). We hence define the optical depth τ on the basis of the Rosseland mean opacity, such that, along the vertical direction $\mathrm{d}\mathit{\tau }\mathrm{\equiv }{\mathit{\kappa }}_{\mathrm{R}}\mathrm{d}\mathit{m}\mathit{.}$(3)Assuming hydrostatic equilibrium, the relation between pressure and optical depth can be found by integrating Eq. (3): $\mathit{\tau }\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{=}{\int }_{\mathrm{0}}^{\mathit{P}}\frac{{\mathit{\kappa }}_{\mathrm{R}}\mathrm{\left(}{\mathit{P}}^{\mathrm{\prime }}\mathit{,T}\mathrm{\right(}{\mathit{P}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}}{\mathit{g}}\mathrm{d}{\mathit{P}}^{\mathrm{\prime }}\mathit{.}$(4)The optical depth thus becomes the natural variable to account for the dependence on depth in the radiative transfer problem. For any strictly positive Rosseland mean opacities, Eq. (4)is a bijection relating pressure and optical depth. Thus, solution of the radiative transfer equations in terms of optical depth can be converted to a solution in term of pressure for any functional form of the Rosseland mean opacities.

The second mean opacity that is traditionally used for radiative transfer is the so-called Planck mean: ${\mathit{\kappa }}_{\mathrm{P}}\mathrm{\equiv }{\mathrm{(}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{B}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}{\mathit{\nu }}^{\mathrm{)}}}^{-1}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathit{.}$(5)We use the ratio of the Planck and Rosseland mean opacities to quantify the non-greyness of the atmosphere: ${\mathit{\gamma }}_{\mathrm{P}}\mathrm{\equiv }\frac{{\mathit{\kappa }}_{\mathrm{P}}}{{\mathit{\kappa }}_{\mathrm{R}}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(6)While the value of the Rosseland mean opacity is dominated by the lowest values of the opacity function κ_ν, the Planck mean opacity is dominated by its highest values. Thus, it can be shown that γ_P = 1 for a grey atmosphere and γ_P > 1 for a non-grey atmosphere (King 1956).

In irradiated atmospheres, a collimated flux coming from the star is absorbed at different atmospheric levels. We name κ_v the opacity relevant to the absorption of the stellar flux. As will be shown in Sect. 3.8, the absorption of the visible flux appears linearly in the radiative transfer equations. Thus, a solution can be found using multiple visible opacity bands κ_v1, κ_v2, etc.

We further define the ratio of the visible opacity to the mean (Rosseland) thermal opacity: ${\mathit{\gamma }}_{\mathrm{v}}\mathrm{\equiv }{\mathit{\kappa }}_{\mathrm{v}}\mathit{/}{\mathit{\kappa }}_{\mathrm{R}}\mathit{.}$(7)In order to solve the radiative transfer problem analytically, we suppose that γ_v is constant with optical depth. Once γ_v is chosen, we can solve the equations for the visible radiation independently from the final thermal structure of the atmosphere. Of course, purely grey models are such that γ_v = 1.

2.1.3. The picket-fence model

Fig. 2 Simplified thermal opacities for the picket-fence model. β = δν/Δν is the equivalent bandwidth (see text).

It is important to note at this point that two sets of opacities with different wavelength dependences may have the same Rosseland and Planck means. We must constrain the problem further, and to this intent, we now consider the simplest possible line model, known as the picket-fence model (Mihalas 1978), where the thermal opacities can take two different values κ₁ and κ₂ (see Fig. 2) such that ${\mathit{\kappa }}_{\mathit{\nu }}\mathrm{=}\{\begin{array}{c}\\ {\mathit{\kappa }}_{\mathrm{1}}& \mathrm{for}\mathit{\nu }\mathrm{\in }\mathrm{[}\mathit{n}\mathrm{\Delta }\mathit{\nu ,n}\mathrm{\Delta }\mathit{\nu }\mathrm{+}\mathit{\delta \nu }\mathrm{]}\\ {\mathit{\kappa }}_{\mathrm{2}}& \mathrm{for}\mathit{\nu }\mathrm{\in }\mathrm{[}\mathit{n}\mathrm{\Delta }\mathit{\nu }\mathrm{+}\mathit{\delta \nu ,}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{\left)}\mathrm{\Delta }\mathit{\nu }\mathrm{\right]}\end{array}\mathit{n}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,N}\mathrm{\right]}\mathit{.}$(8)We define an equivalent bandwidth by: $\mathit{\beta }\mathrm{=}\frac{\mathrm{1}}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{B}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }}\sum _{\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}{\int }_{\mathit{n}\mathrm{\Delta }\mathit{\nu }}^{\mathit{n}\mathrm{\Delta }\mathit{\nu }\mathrm{+}\mathit{\delta \nu }}{\mathit{B}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathit{.}$(9)The characteristic width of the Planck function can be defined as $\mathrm{\Delta }{\mathit{\nu }}_{\mathit{P}}\mathrm{\equiv }{{}^{\mathrm{(}}\frac{\mathrm{1}}{\mathit{B}}{\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\nu }}}{\mathit{\partial \nu }}}^{\mathrm{)}}}^{-1}$. When choosing Δν ≪ Δν_P, the Planck function can be considered constant over Δν and we get β = δν/Δν. The Planck and Rosseland mean opacities then become (see Eqs. (2)and (5)) $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{R}}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\kappa }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\kappa }}_{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}{\mathit{\kappa }}_{\mathrm{1}}}\\ {\mathit{\kappa }}_{\mathrm{P}}& \mathrm{=}& \mathit{\beta }{\mathit{\kappa }}_{\mathrm{1}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}{\mathit{\kappa }}_{\mathrm{2}}\mathit{.}\end{array}$We also define the following ratios: $\begin{array}{ccc}& & {\mathit{\gamma }}_{\mathrm{1}}\mathrm{\equiv }{\mathit{\kappa }}_{\mathrm{1}}\mathit{/}{\mathit{\kappa }}_{\mathrm{R}}\\ & & {\mathit{\gamma }}_{\mathrm{2}}\mathrm{\equiv }{\mathit{\kappa }}_{\mathrm{2}}\mathit{/}{\mathit{\kappa }}_{\mathrm{R}}\\ & & \mathit{R}\mathrm{\equiv }{\mathit{\kappa }}_{\mathrm{1}}\mathit{/}{\mathit{\kappa }}_{\mathrm{2}}\mathrm{=}{\mathit{\gamma }}_{\mathrm{1}}\mathit{/}{\mathit{\gamma }}_{\mathrm{2}}\mathit{.}\end{array}$Following Chandrasekhar (1935), we also define a limit optical depth ${\mathit{\tau }}_{\lim }\mathrm{\equiv }\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}}\sqrt{\frac{{\mathit{\gamma }}_{\mathrm{P}}}{\mathrm{3}}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(15)The role of τ_lim in shaping the final temperature profile is discussed in Sect. 5 and its variations with R and β is pictured in Fig. 8.

2.2. The method of discrete ordinates for the non-irradiated problem

2.2.1. The grey case

An approximate method to solve Eq. (1) including the angular dependency has been developed by Chandrasekhar (1960) in the case of a non-irradiated atmosphere (T_irr = 0). The idea is to replace the integrals over angle in Eq. (1) by a Gaussian sum over μ. It can then be solved to an arbitrary precision by increasing the number of terms in the sum. The boundary condition at the top of the atmosphere is simply given as I_{μ    <    0}(0) = 0. The expansion to the fourth term yields the temperature profile $\mathit{T}\mathrm{(}\mathit{\tau }{\mathrm{)}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\mathrm{(}\mathit{\tau }\mathrm{+}\mathit{Q}\mathrm{+}{\mathit{L}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathrm{1}}\mathit{\tau }}\mathrm{+}{\mathit{L}}_{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathrm{2}}\mathit{\tau }}\mathrm{+}{\mathit{L}}_{\mathrm{3}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathrm{3}}\mathit{\tau }}\mathrm{)}\mathit{,}$(16)with Q = 0.706920, L₁ = −0.083921, L₂ = −0.036187, L₃ = −0.009461, k₁ = 1.103188, k₂ = 1.591778, and k₃ = 4.45808 (Chandrasekhar 1960, Table VIII)². One of the important results from this formalism is that the skin temperature of the planet (defined as the temperature at zero optical depth) is independent of the order of expansion and therefore corresponds to the exact value ${\mathit{T}}_{\mathrm{skin}}^{\mathrm{4}}\mathrm{=}\sqrt{\mathrm{3}}\frac{{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(17)This expression is exact only in the limit of a grey, non-irradiated atmosphere.

2.2.2. The non-grey case

Chandrasekhar (1960) also developed a perturbation method in order to include non-grey thermal opacities. This method was improved by Krook (1963). However, these perturbation methods either work for small departures from the grey opacities or involve a fastidious iterative procedure (e.g. Unno & Yamashita 1960; Avrett & Krook 1963) and are no longer fully analytical. However, considering that the variations in the opacities are small compared to the variations of the Planck function, analytical solutions can be found for an arbitrarily large departure from the grey opacities. Noting the similar role of μ and κ_ν in Eq. (1), King (1956), following Münch (1946), used the method of discrete ordinates in order to turn the integrals over frequency into Gaussian sums. For the picket-fence model defined in Sect. 2.1.3 and the second approximation for the angular dependency, King’s method leads to the following temperature profile ${\mathit{T}}^{\mathrm{4}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\left[\frac{\mathrm{1}}{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\mathrm{+}\mathit{\tau }\mathrm{+}\frac{\mathrm{(}\sqrt{{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{)}\mathrm{(}\sqrt{{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}}{{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{-}\mathrm{1}\mathrm{)}\right]\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(18)As in the grey case, the method of discrete ordinates leads to an exact relation for the skin temperature, whatever the dependency of κ_ν on frequency (but no dependence on pressure or temperature): ${\mathit{T}}_{\mathrm{skin}}^{\mathrm{4}}\mathrm{=}\sqrt{\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{P}}}}\frac{{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(19)In the grey limit, γ_P = 1 and we recover Eq. (17). Otherwise, γ_P > 1, implying that for a non-irradiated atmosphere, non-grey effects always tend to lower the atmospheric skin temperature.

2.3. Moment equation method

2.3.1. Equations for the momentum of the radiation intensity

A simpler way to solve the radiative transfer equation has been carried out by Eddington (1916). The idea is to solve the equation using the different momentum of the intensity defined as $\mathrm{\left(}{\mathit{J}}_{\mathit{\nu }}\mathit{,}{\mathit{H}}_{\mathit{\nu }}\mathit{,}{\mathit{K}}_{\mathit{\nu }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\int }_{-1}^{\mathrm{1}}{\mathit{I}}_{\mathit{\mu \nu }}\mathrm{\left(}\mathrm{1}\mathit{,\mu ,}{\mathit{\mu }}^{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathit{.}$(20)Then, integrating over μ Eq. (1)and μ times Eq. (1)one gets the momentum equations $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{H}}_{\mathit{\nu }}}{\mathrm{d}{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{=}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{.}\\ & & \frac{\mathrm{d}{\mathit{K}}_{\mathit{\nu }}}{\mathrm{d}{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{=}{\mathit{H}}_{\mathit{\nu }}\mathit{.}\end{array}$Assuming the atmosphere to be in radiative equilibrium, we can write ${\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{(}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{=}\mathrm{0.}$(23)For a grey atmosphere (κ_ν = κ_R  ∀ν), Eqs. (21)−(23) can be integrated over frequency, leading to an equation for J, H, K, and B, the frequency-integrated versions of J_ν, H_ν, K_ν, and B_ν. The radiative equilibrium equation becomes $\mathit{J}\mathrm{=}\mathit{B}\mathit{.}$(24)The frequency integrated versions of Eqs. (21)−(23) are a set of three equations with four unknowns. The system is not closed because by integrating Eq. (1) over all angles we have lost the information on the angular dependency of the irradiation. A closure relationship that contains this angular dependency is therefore needed. A common closure relationship, known as the Eddington approximation is $\mathit{J}\mathrm{=}\mathrm{3}\mathit{K}\mathit{.}$(25)This relationship is exact in two very different cases: when the radiation field is isotropic (I_μ independent of μ), and in the two-stream approximation (${\mathit{I}}_{\mathit{\mu }\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathrm{0}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{+}}$ and ${\mathit{I}}_{\mathit{\mu }\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{0}}\mathrm{=}{\mathit{I}}_{\mathrm{0}}^{\mathrm{-}}$). Although this seems to be a very restrictive approximation, it is relevant for the deep layers of the atmosphere because of the quasi-isotropy of the radiation field there. It is also good for the top of the atmosphere, where the flux comes mainly from the τ ≈ 1 layer. Indeed, the exact solution gives a ratio J/K that differs by no more than 20% from the 1/3 ratio over the whole atmosphere and leads to a temperature profile that is correct to 4% in the grey case (see the plain blue line in Fig. 4).

2.3.2. Top boundary condition

Although in the method of discrete ordinates the boundary condition at the top of the atmosphere is intuitive, in the momentum equations method it is less obvious and different choices have been made by different authors. Usually, the expression for J is known and some integration constants need to be found. Two equations are needed, one for J(0) and one for H(0). Four possibilities are widely used in the literature from which one has to choose two:

• 1.

  The radiative equilibrium equation that relates the emergent flux at the top of the atmosphere to the internal flux from the planet and the incident flux from the star.

• 2.

  An ad-hoc relation between H(0) and J(0) at τ = 0: H(0) = f_HJ(0), where f_H is often called the second Eddington coefficient.

• 3.

  A calculation of H(0) from the second moment equation (Eq. (22)) and the Eddington approximation.

• 4.

  A calculation of H(0) from the integration of the source function through the entire atmosphere, known as the Milne equation (Mihalas & Mihalas 1984, p. 347): $\mathit{H}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathit{B}{}^{\mathrm{\right(}}\mathit{\tau }^{\mathrm{\left)}}{\mathit{E}}_{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\tau }$ .

For grey and semi-grey models, the first condition is natural and so it was used by Hansen (2008) and Guillot (2010). For the other part of the top boundary condition, Guillot (2010) chose to use the second and Hansen (2008) the fourth condition (see Appendix A of Guillot 2010 for a comparison of the two expressions).

In the case of a non-grey model, the first condition cannot be implemented (at least directly) because it is a constraint on the total thermal flux, but it provides no information on how the thermal flux is split between the opacity bands that are considered. Chandrasekhar (1935) therefore uses conditions 2 and 3 in each of the opacity bands for his non-grey, non-irradiated model. He also notes that using condition 4 instead of condition 3 should yield better results, but it leads to more complex expressions. In this work, because an accurate treatment of the flux is needed for the non-grey irradiated model, we will use conditions 2 and 4 in each of the opacity bands. All these models are discussed in the next sections and summarized in Table 1.

2.3.3. Non-irradiated, grey case

In this section we consider the case T_irr = 0. Under the grey approximation, using the conditions 1 and 4, the temperature profile is given by (Mihalas & Mihalas 1984, p. 357) ${\mathit{T}}^{\mathrm{4}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\left(\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\mathit{\tau }\right)\mathit{,}$(26)which leads to the same solution as assuming condition 2 with $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$. The skin temperature is then ${\mathit{T}}_{\mathrm{skin}}^{\mathrm{4}}\mathrm{=}\mathrm{2}\frac{{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}$(27)which differs from the exact solution (Eq. (17)) by a factor of $\sqrt{\mathrm{3}}$/2. Assuming $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}$/$\sqrt{\mathrm{3}}$ is thus tempting, as it leads to the correct skin temperature. Unfortunately, it leads to a temperature profile which is less accurate around τ ≈ 1.

2.3.4. Non-irradiated non-grey picket-fence model

Chandrasekhar (1935) provides solutions to the moment equations for the picket-fence model presented in Sect. 2.1.3. He assumes that the relation $\mathit{H}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\stackrel{\mathrm{h}}{.}}\mathit{J}\mathrm{\left(}\mathrm{0}\mathrm{\right)}$ with $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$, valid in the grey case under the Eddington approximation (see Eq. (27)), holds for the two thermal channels separately. Using this condition together with condition 3, he obtains a temperature profile $\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\left[\mathit{\tau }\mathrm{+}\frac{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\sqrt{\frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\right]\\ & & \mathrm{+}\frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\left(\frac{{\mathit{\gamma }}_{\mathrm{P}}\mathrm{-}\mathrm{1}}{\sqrt{{\mathit{\gamma }}_{\mathrm{P}}}}\right)\frac{\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}\mathrm{+}\sqrt{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{\tau }}_{\lim }}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{)}\end{array}$(28)and an equation for the skin temperature ${\mathit{T}}_{\mathrm{skin}}^{\mathrm{4}}\mathrm{=}\mathrm{2}\left(\frac{\mathrm{2}\mathrm{+}\sqrt{\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{P}}}}}{\mathrm{2}\mathrm{+}\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\right)\frac{{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\mathit{,}$(29)As expected, this equation reduces to Eq. (27) in the limit γ_P = 1. For high values of γ_P, this relation differs by a factor of 4/3 from the exact one derived with the method of discrete ordinates [19]. As happens for the grey case, using $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}$/$\sqrt{\mathrm{3}}$ would lead to the exact solution for the skin temperature, but at the expense of the accuracy of the profile at deeper levels. Again, we note that, in the non-irradiated case, the temperature at the top of the atmosphere is determined by a single parameter, γ_P, representing the non-greyness of the atmosphere. A comparison of the different expressions for the skin temperature is provided in Sect. 4.

2.3.5. Irradiated semi-grey model

In the case of irradiated atmospheres, the presence of an incoming collimated flux at the top of the atmosphere breaks the angular symmetry of the equations. The radiative transfer problem thus can no longer be solved analytically (at least not in a simple way) through the discrete ordinates technique. The momentum method is thus required.

To solve the problem, the radiation field is split into two parts: the incoming, collimated radiation field on one hand, the thermal radiation field on the other. The radiative equilibrium equation (Eq. (23)) links the two streams, as can be seen in Sect. 3.1 (see also Hansen 2008; Guillot 2010; Robinson & Catling 2012). As mentioned in Sect. 2.1.2, when the incident radiation is at a much shorter wavelength than the thermal emission of the atmosphere, the two streams correspond to different characteristic wavelengths and may often be labelled as visible and infrared. This is not a requirement, however: the solutions apply if the radiation field corresponds to other wavelengths or if they overlap.

As discussed previously (Sect. 2.3.2), the boundary condition at the top of the model can be chosen in several ways. When using condition 2, Guillot (2010) lets the value of $\stackrel{\mathrm{h}}{.}$ be either 1/2 or $\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$. Those values are based on the non-irradiated case: $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$ is the value that arises from the calculation of the angle dependence between H(0) and J(0) in the isotropic case, but $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ provides a skin temperature that agrees with the exact value. The two solutions differ by ≈3% at most (see Fig. 5 hereafter), and choosing one over another is not crucial. In any case, for an easier comparison, we provide here the solution of Guillot (2010) for $\stackrel{\mathrm{h}}{.}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}$$\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}& \mathrm{=}& \frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\left[\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\mathit{\tau }\right]\\ & & \mathrm{+}\frac{\mathrm{3}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}}{\mathrm{4}}{\mathit{\mu }}_{\mathrm{\ast }}\left[\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\frac{{\mathit{\mu }}_{\mathrm{\ast }}}{{\mathit{\gamma }}_{\mathrm{v}}}\mathrm{+}\left(\frac{{\mathit{\gamma }}_{\mathrm{v}}}{\mathrm{3}{\mathit{\mu }}_{\mathrm{\ast }}}\mathrm{-}\frac{{\mathit{\mu }}_{\mathrm{\ast }}}{{\mathit{\gamma }}_{\mathrm{v}}}\right){\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}\mathit{\tau }\mathit{/}{\mathit{\mu }}_{\mathrm{\ast }}}\right]\mathit{,}\end{array}$(30)where μ_∗ is the cosine of the angle of the incident radiation. The skin temperature is ${\mathit{T}}_{\mathrm{skin}}^{\mathrm{4}}\mathrm{=}\mathrm{2}\frac{{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{+}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}}{\mathrm{4}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{v}}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathit{.}$(31)For γ_v → 0, the incident radiation is absorbed in the deep layers of the atmosphere and the skin temperature converges to the skin temperature of a grey model with an effective temperature ${\mathit{T}}_{\mathrm{eff}\mathit{,}\hspace{0.17em}{\mathit{\mu }}_{\mathrm{\ast }}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{+}{\mathit{\mu }}_{\mathrm{\ast }}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$. The semi-grey model depends only on the parameter γ_v.

As discussed in the introduction, the semi-grey model predicts minimum temperatures that are generally higher than the numerical solutions for irradiated exoplanets, independent of the choice of γ_v (see Fig. 1). In fact, similarly to the skin temperature, the minimum temperature of a semi-grey atmosphere, shown in Fig. 3, depends only on the values of T_{eff,   μ∗} and γ_v. It is the lowest and equal to T_{eff,   μ∗}/2^1/4 both in the γ_v → 0 and γ_v → ∞ limits. This lower bound for the semi-grey temperature profile is hotter than the one obtained by numerical calculations taking into account the full set of opacities. The discrepancy is much larger than the variations resulting from the approximation of the momentum method. Clearly, non-grey effects must be invoked to explain the low temperatures obtained by numerical models at low optical depths.

Fig. 3 Minimum temperature of the semi-grey model in terms of the effective temperature as a function of γv/μ∗.

3. An analytical irradiated non-grey picket-fence model

3.1. Equations

We now derive the equations for an irradiated atmosphere in local thermodynamic equilibrium with infrared line opacities as described in Sect. 2.1.2. Thus, our model contains three different opacities: κ₁ and κ₂ for the thermal radiation and κ_v for the incoming radiation of the star. As explained before, the difference between the thermal and the visible channel depends on the angular dependency of the radiation and not on the frequency. Although the method of discrete ordinates is shown to lead to more exact results, it is difficult to adapt to the irradiated case. Therefore, following Chandrasekhar (1935) and Guillot (2010), we solve the radiative transfer equations using the momentum equations. Integrating Eqs. (21) and (22) over each thermal band we obtain

$\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{H}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}}{\mathrm{d}\mathit{\tau }}\mathrm{=}{\mathit{\gamma }}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}{\mathit{J}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{(}\mathit{\beta ,}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\mathit{B}\mathrm{\left(}\mathit{T}\mathrm{\right)}\\ & & \frac{\mathrm{d}{\mathit{K}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}}{\mathrm{d}\mathit{\tau }}\mathrm{=}{\mathit{\gamma }}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}{\mathit{H}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathit{,}\end{array}$where the subscript indicates the integrated quantities over the given thermal band. Thus, for a quantity X_ν we have:

$\begin{array}{ccc}& & {\mathit{X}}_{\mathrm{1}}\mathrm{=}\sum _{\mathit{n}}{\int }_{\mathit{n}\mathrm{\Delta }\mathit{\nu }}^{\mathit{n}\mathrm{\Delta }\mathit{\nu }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\delta \nu }}{\mathit{X}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\\ & & {\mathit{X}}_{\mathrm{2}}\mathrm{=}\sum _{\mathit{n}}{\int }_{\mathit{n}\mathrm{\Delta }\mathit{\nu }\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{\delta }}_{\mathit{\nu }}}^{\mathrm{(}\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathrm{\Delta }\mathit{\nu }}{\mathit{X}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathit{.}\end{array}$The Planck function is considered constant over each bin of frequency Δν and therefore B₁ = βB and B₂ = (1 − β)B.

We now assume that the Eddington approximation is valid in the two bands separately: J_(1,2) = 3K_(1,2). Equations (32)and (33)can be combined into $\begin{array}{ccc}& & \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{\mathrm{1}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{=}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{(}{\mathit{J}}_{\mathrm{1}}\mathrm{-}\mathit{\beta B}\mathrm{)}\\ & & \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{\mathrm{2}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{=}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{(}{\mathit{J}}_{\mathrm{2}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{\left)}\mathit{B}\mathrm{\right)}\end{array}$and the radiative equilibrium equation becomes ${\mathit{\gamma }}_{\mathrm{1}}{\mathit{J}}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}{\mathit{J}}_{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{v}}{\mathit{J}}_{\mathrm{v}}\mathrm{=}{\mathit{\gamma }}_{\mathrm{P}}\mathit{B,}$(38)where the quantities with subscript v are the momentum of the incident stellar radiation. Assuming that the incoming stellar radiation arrives as a collimated flux and hit the top of the atmosphere with an angle θ_∗ = cos^-1μ_∗ they can are given by $\mathrm{\left(}{\mathit{J}}_{\mathrm{v}}\mathit{,}{\mathit{H}}_{\mathrm{v}}\mathit{,}{\mathit{K}}_{\mathrm{v}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{1}\mathit{,}{\mathit{\mu }}_{\mathrm{\ast }}\mathit{,}{\mathit{\mu }}_{\mathrm{\ast }}^{\mathrm{2}}\mathrm{\right)}\mathit{I}\ast$(39)with ${\mathit{I}}_{\mathrm{\ast }}\mathrm{=}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{I}}_{\mathrm{\ast }\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }$ the total incident intensity.

The absorption of the stellar irradiation can be treated separately from the thermal radiation and J_v is given by Eq. (13) of Guillot (2010), ${\mathit{J}}_{\mathrm{v}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{{\mathit{\mu }}_{\mathrm{\ast }}}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{,}$(40)where we have simplified the notation by introducing the parameter ${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{\equiv }{\mathit{\gamma }}_{\mathrm{v}}\mathit{/}{\mathit{\mu }}_{\mathrm{\ast }}$.

Equations (36) to (38) are a set of three coupled equations with three unknowns J₁, J₂, and B. In order to decouple these equations we define two new variables: $\{\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{J}}_{\mathit{\gamma }}\mathrm{\equiv }\frac{{\mathit{J}}_{\mathrm{1}}}{{\mathit{\gamma }}_{\mathrm{1}}}\mathrm{+}\frac{{\mathit{J}}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{2}}}\\ \begin{array}{c}\end{array}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{\equiv }\frac{{\mathit{J}}_{\mathrm{1}}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{3}}}\mathrm{+}\frac{{\mathit{J}}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}\end{array}$(41)Conversely, we can come back to the original variables: $\{\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{J}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{3}}\mathrm{(}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{-}{\mathit{J}}_{\mathit{\gamma }}\mathrm{)}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}\\ \begin{array}{c}\end{array}{\mathit{J}}_{\mathrm{2}}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{3}}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{-}{\mathit{J}}_{\mathit{\gamma }}\mathrm{)}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}\end{array}$(42)Using the combination of Eqs. $\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}}}$(36) + $\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{2}}}$(37) and (38)we get $\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{\mathit{\gamma }}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\mathrm{3}{\mathit{\gamma }}_{\mathrm{v}}{\mathit{J}}_{\mathrm{v}}\mathit{.}$(43)The combination of Eqs. $\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{3}}}$(36)+$\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{3}}}$(37) yields $\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{=}\mathrm{3}\mathrm{(}{\mathit{J}}_{\mathit{\gamma }}\mathrm{-}\mathit{B}\mathrm{)}\mathit{.}$(44)Noting that ${\mathit{\gamma }}_{\mathrm{1}}{\mathit{J}}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}{\mathit{J}}_{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{)}{\mathit{J}}_{\mathit{\gamma }}\mathrm{-}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}$, Eq. (38) becomes $\mathit{B}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{J}}_{\mathit{\gamma }}\mathrm{-}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{v}}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{J}}_{\mathrm{v}}\mathit{.}$(45)Equations (43)−(45) are now a set of two uncoupled differential equations and a linear equation.

3.2. Boundary conditions

To solve the differential equations we need to specify the boundary conditions. When τ →  + ∞ we want to fulfill the diffusion approximation J_ν = B_ν (Mihalas & Mihalas 1984, p. 350). In our case this translates to J₁ = βB and J₂ = (1 − β)B. Furthermore, at these levels, the gradient of B should also obey the diffusion approximation (Mihalas & Mihalas 1984) $\frac{\mathrm{d}\mathit{B}}{\mathrm{d}\mathit{\tau }}\underset{\mathit{\tau }\mathrm{\to }\mathrm{+}\mathrm{\infty }}{\mathrm{~}}\mathrm{3}{\mathit{H}}_{\mathrm{\infty }}\mathit{,}$(46)where $\mathrm{4}\mathit{\pi }{\mathit{H}}_{\mathrm{\infty }}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}$ is the thermal flux coming from the interior of the planet. Using the system of Eqs. (41) and noting that $\frac{\mathit{\beta }}{{\mathit{\gamma }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}\mathit{\beta }}{{\mathit{\gamma }}_{\mathrm{2}}}\mathrm{=}\mathrm{1}$, we can derive a condition on J_γ and J_γ³: $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{J}}_{\mathit{\gamma }}}{\mathrm{d}\mathit{\tau }}\underset{\mathit{\tau }\mathrm{\to }\mathrm{+}\mathrm{\infty }}{\mathrm{~}}\mathrm{3}{\mathit{H}}_{\mathrm{\infty }}\\ & & \frac{\mathrm{d}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}}{\mathrm{d}\mathit{\tau }}\underset{\mathit{\tau }\mathrm{\to }\mathrm{+}\mathrm{\infty }}{\mathrm{~}}\left(\frac{\mathit{\beta }}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{3}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}\mathit{\beta }}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{3}}}\right)\mathrm{3}{\mathit{H}}_{\mathrm{\infty }}\hspace{0.17em}\mathit{.}\end{array}$For τ → 0 we specify the geometry of the intensity by setting: ${\mathit{J}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathit{H}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(49)Furthermore, we calculate the flux at the top of the atmosphere in each band using Eq. (79.21) from Mihalas & Mihalas (1984). From the assumption of local thermodynamic equilibrium, the source function in the two bands is S₁(τ₁) = βB(τ/γ₁) and S₂(τ₂) = (1 − β)B(τ/γ₂). The upper boundary condition on the flux of the two bands thus becomes ${\mathit{H}}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{\left(}\mathit{\beta ,}\mathrm{\right(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{\left)}\mathrm{\right)}\mathit{B}\left(\frac{\mathit{\tau }}{{\mathit{\gamma }}_{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}}\right){\mathit{E}}_{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\tau }\hspace{0.17em}\mathit{.}$(50)

3.3. Solution

The solution of a second-order differential equation with constant coefficient is the sum of the solutions of the homogeneous equation and a particular solution of the complete equation. Thus, solutions of Eq. (43) must be of the form: ${\mathit{J}}_{\mathit{\gamma }}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{1}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}\mathit{\tau }\mathrm{+}\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{.}$(51)Applying the boundary condition Eq. (47), we get C₂ = 3H. For τ = 0 we obtain ${\mathit{J}}_{\mathit{\gamma }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{1}}\mathrm{+}\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(52)Using Eq. (45)to eliminate B and replacing J_γ by its solution, Eq. (44) becomes: $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{-}\mathrm{3}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{=}& & \mathrm{3}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}\right)\mathrm{(}{\mathit{C}}_{\mathrm{1}}\mathrm{+}\mathrm{3}\mathit{H\tau }\mathrm{)}\\ & & \mathrm{+}\mathrm{3}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}\right)\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\\ & & \mathrm{+}\mathrm{3}\frac{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{.}\end{array}$(53)Again, solutions of this differential equation must be the sum of the solutions of the homogeneous equation and one solution of the complete equation. The homogeneous solution must have the form ${\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}\mathrm{H}}\mathrm{=}{\mathit{C}}_{\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{+}{\mathit{C}}_{\mathrm{4}}{\mathrm{e}}^{\mathrm{+}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathit{,}$(54)where C₃ and C₄ are constants of integration to be determined using the boundary conditions. We look for a particular solution formed by the superposition of an exponential and an affine function. The affine function must then be a solution of Eq. (53) with H_v(0) = 0${\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}\mathrm{P}\mathrm{1}}\mathrm{=}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{P}}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}\right)\mathrm{(}{\mathit{C}}_{\mathrm{1}}\mathrm{+}\mathrm{3}\mathit{H\tau }\mathrm{)}\mathit{,}$(55)and the exponential function must be solution of Eq. (53), keeping only the exponential part on the right-hand side ${\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}\mathrm{P}\mathrm{2}}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{P}}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}{\mathit{\tau }}_{\lim }{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\left[\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}\right)\frac{\mathrm{3}}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{{\mathit{\gamma }}_{\mathrm{P}}}\right]{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{.}$(56)Applying the boundary condition defined by Eq. (48) to the full solution J_γ³ = J_γ³P1 + J_γ³P2 + J_γ³H, we find C₄ = 0. The full solution of Eq. (44)is therefore given by ${\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}\mathrm{P}\mathrm{1}}\mathrm{+}{\mathit{J}}_{{\mathit{\gamma }}^{\mathrm{3}}\mathrm{P}\mathrm{2}}\mathrm{+}{\mathit{C}}_{\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathit{.}$(57)We can get an expression for the source function by replacing J_γ and J_γ³ in the radiative equilibrium equation (Eq. (45)): $\begin{array}{ccc}\mathit{B}& \mathrm{=}& {\mathit{C}}_{\mathrm{1}}\mathrm{+}\mathrm{3}\mathit{H\tau }\mathrm{-}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{C}}_{\mathrm{3}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\\ & & \mathrm{+}\frac{\mathrm{[}\mathrm{3}\mathrm{-}\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{/}{\mathit{\gamma }}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}\mathrm{]}\mathrm{[}\mathrm{3}\mathrm{-}\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{/}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{]}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{)}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{.}\end{array}$(58)To get the complete solution of the problem, we need to determine the two remaining integration constants C₁ and C₃ using the boundary condition (49). For that we need to calculate J₁(0), J₂(0), H₁(0) and H₂(0). The first two quantities can be evaluated by inserting the values of J_γ(0) and J_γ³(0) from Eqs. (51) and (57) into system (42): $\begin{array}{ccc}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}\left({\mathit{\gamma }}_{\mathrm{2}}\mathrm{-}\mathrm{1}\right)}{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{C}}_{\mathrm{1}}\\ & & \mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{(}\mathrm{3}\left({\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\right)\left({\mathit{\gamma }}_{\mathrm{2}}\mathrm{-}\mathrm{1}\right)\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{{\mathit{\tau }}_{\lim }^{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)} \\ & & \mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{3}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{\mathit{C}}_{\mathrm{3}}\\ {\mathit{J}}_{\mathrm{2}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\gamma }}_{\mathrm{2}}\left({\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}\mathrm{1}\right)}{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{C}}_{\mathrm{1}}\\ & & \mathrm{-}\frac{{\mathit{\gamma }}_{\mathrm{2}}\mathrm{(}\mathrm{3}\left({\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\right)\left({\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}\mathrm{1}\right)\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{{\mathit{\tau }}_{\lim }^{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)} \\ & & \mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{3}}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}{\mathit{C}}_{\mathrm{3}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}\end{array}$Noting that

$\begin{array}{ccc}& & {\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{E}}_{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{d}\mathit{\tau }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\\ & & {\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\tau }{\mathit{E}}_{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{d}\mathit{\tau }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\\ & & {\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}\mathit{\alpha \tau }}{\mathit{E}}_{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\tau }\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathit{\alpha }}\mathrm{-}\frac{\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{)}}{{\mathit{\alpha }}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}\end{array}$we can evaluate H₁(0) and H₂(0) by inserting Eqs. (58) into (50). Then, Eq. (49) is a linear system of two equations with two unknowns. After some calculations we get the expressions for C₁ and C₃$\begin{array}{ccc}{\mathit{C}}_{\mathrm{1}}& \mathrm{=}& \mathrm{(}{\mathit{a}}_{\mathrm{0}}\mathrm{+}{\mathit{a}}_{\mathrm{1}}{\mathit{b}}_{\mathrm{0}}\mathrm{)}\mathit{H}\mathrm{+}\mathrm{\left(}{\mathit{a}}_{\mathrm{1}}{\mathit{b}}_{\mathrm{0}}{\mathit{b}}_{\mathrm{1}}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{3}}\mathrm{)}\mathrm{+}{\mathit{a}}_{\mathrm{2}}\mathrm{+}{\mathit{a}}_{\mathrm{3}}\mathrm{)}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)} \\ {\mathit{C}}_{\mathrm{3}}& \mathrm{=}& {\mathit{b}}_{\mathrm{0}}\mathit{H}\mathrm{+}{\mathit{b}}_{\mathrm{0}}{\mathit{b}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{3}}\mathrm{)}{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{,}\end{array}$where we have $\begin{array}{ccc}{\mathit{a}}_{\mathrm{0}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{2}}}\\ {\mathit{a}}_{\mathrm{1}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{3}{\mathit{\tau }}_{\lim }^{\mathrm{2}}}\left[\frac{{\mathit{\gamma }}_{\mathrm{p}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{p}}}\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{-}\mathrm{2}}{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}}\mathrm{+}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}{\mathit{\tau }}_{\lim }\mathrm{-}\mathrm{(}{\mathit{A}}_{\mathrm{t}\mathit{,}\mathrm{1}}\mathrm{+}{\mathit{A}}_{\mathrm{t}\mathit{,}\mathrm{2}}\mathrm{)}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\right] \\ {\mathit{a}}_{\mathrm{2}}& \mathrm{=}& \frac{{\mathit{\tau }}_{\lim }^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}}\\ \mathrm{\times }& & \frac{\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}}^{\mathrm{)}}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}}^{\mathrm{)}}\left({\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\right)\mathrm{-}\mathrm{3}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{(}\mathrm{6}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{)}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}} \\ {\mathit{a}}_{\mathrm{3}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{A}}_{\mathrm{v}\mathit{,}\mathrm{2}}\mathrm{+}{\mathit{A}}_{\mathrm{v}\mathit{,}\mathrm{1}}\mathrm{)}}{{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{)}}\\ {\mathit{b}}_{\mathrm{0}}& \mathrm{=}& {\left(\frac{{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}\frac{{\mathit{A}}_{\mathrm{t}\mathit{,}\mathrm{1}}\mathrm{-}{\mathit{A}}_{\mathrm{t}\mathit{,}\mathrm{2}}}{\mathrm{3}}\mathrm{-}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}}\mathrm{-}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}}\right)}^{-1} \\ {\mathit{b}}_{\mathrm{1}}& \mathrm{=}& \frac{{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}{\mathit{\tau }}_{\lim }^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{(}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}{\mathit{\tau }}_{\lim }^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\\ {\mathit{b}}_{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{3}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{3}}}{\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}\mathrm{(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{)}}\\ {\mathit{b}}_{\mathrm{3}}& \mathrm{=}& \frac{{\mathit{A}}_{\mathrm{v}\mathit{,}\mathrm{2}}\mathrm{-}{\mathit{A}}_{\mathrm{v}\mathit{,}\mathrm{1}}}{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}\mathrm{)}}\mathit{,}\end{array}$ where we defined $\begin{array}{ccc}{\mathit{A}}_{\mathrm{t}\mathit{,i}}& \mathrm{=}& {\mathit{\gamma }}_{\mathit{i}}^{\mathrm{2}}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\tau }}_{\lim }{\mathit{\gamma }}_{\mathit{i}}}\right)\\ {\mathit{A}}_{\mathrm{v}\mathit{,i}}& \mathrm{=}& {\mathit{\gamma }}_{\mathit{i}}^{\mathrm{2}}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}}{{\mathit{\gamma }}_{\mathit{i}}}\right)\hspace{0.17em}\mathrm{·}\hspace{0.17em}\end{array}$

3.4. Atmospheric temperature profile

Using the relations B = σT⁴/π, $\mathit{H}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{/}\mathrm{4}\mathit{\pi }$, and ${\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{\mu }}_{\mathrm{\ast }}\mathit{\sigma }{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathit{/}\mathrm{4}\mathit{\pi }$ and Eq. (58), we can derive the equation for the temperature at any optical depth $\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}& \mathrm{=}& \frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\mathrm{(}\mathit{\tau }\mathrm{+}\mathit{A}\mathrm{+}\mathit{B}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{)}\\ & & \end{array}$(76)with $\begin{array}{ccc}& & \mathit{A}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}{\mathit{a}}_{\mathrm{0}}\mathrm{+}{\mathit{a}}_{\mathrm{1}}{\mathit{b}}_{\mathrm{0}}\mathrm{)}\\ & & \mathit{B}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{p}}}{\mathit{b}}_{\mathrm{0}}\\ & & \mathit{C}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\left[{\mathit{b}}_{\mathrm{0}}{\mathit{b}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{3}}\mathrm{)}{\mathit{a}}_{\mathrm{1}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}\mathrm{+}{\mathit{a}}_{\mathrm{3}}\right]\\ & & \mathit{D}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{\mathrm{(}{\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{p}}}{\mathit{b}}_{\mathrm{0}}{\mathit{b}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{3}}\mathrm{)}\\ & & \end{array}$

3.5. Grey limit

In the grey limit, γ_P → 1 (as γ₁ and γ ₂ ) and we obtain $\begin{array}{ccc}& & \mathit{A}\mathrm{\to }\mathrm{2}\mathit{/}\mathrm{3}\\ & & \mathit{B}\mathrm{\to }\mathrm{0}\\ & & \mathit{C}\mathrm{\to }\mathrm{2}\mathit{/}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{/}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{/}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{+}\mathrm{2}\log \mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{)}\mathrm{(}\mathrm{1}\mathit{/}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }\mathrm{3}}\mathrm{-}\mathrm{1}\mathit{/}{\mathrm{\left(}\mathrm{3}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{\right)}}^{\mathrm{)}} \\ & & \mathit{D}\mathrm{\to }\mathrm{0}\\ & & \mathit{E}\mathrm{\to }{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{/}\mathrm{3}\mathrm{-}\mathrm{1}\mathit{/}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathit{.}\end{array}$ If we also assume that ${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{\to }\mathrm{0}$ we obtain $\mathit{C}\mathrm{\to }\mathrm{2}\mathit{/}\mathrm{3}\mathrm{+}\mathrm{1}\mathit{/}{\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}$ and E → −1/${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}$ and the solution converges towards that of Guillot (2010) (see Eq. (30)). For other values of ${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}$ our model differs from the solutions of Guillot (2010, see also Hansen 2008) because of the different boundary conditions used in the two models (see Sect. 2.3.2). However, calculations show that the value of C obtained here differs from the same coefficient extracted from Eq. (30) by at most 12% and that the two solutions also converge for ${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}\mathrm{\to }\mathrm{\infty }$. As seen in Fig. 5, in the semi-grey limit, and when calculating the full temperature profile, our model differs by at most 2% from the Guillot (2010) model. The difference between the various solutions must be attributed to the Eddington approximation.

3.6. Using the model

The temperature vs. optical depth profile for our irradiated picket-fence model is given by Eq. (76). The profile has been derived using the Rosseland optical depth as vertical coordinate. It is therefore valid for any functional form of the Rosseland opacities. Equation (4)allows us to switch from τ to P as the vertical coordinate. Although, for convenience, this expression contains four different variables, γ₁,γ₂,γ_P, and τ_lim, it must be kept in mind that, besides the Rosseland mean opacity, there are only two independent variables in the problem. The variables β and R ≡ γ₁/γ₂ = κ₁/κ₂ are the ones to consider to control the shape of the thermal opacities. The variables γ_P and τ_lim are the ones to consider to control the profile itself. The variable γ_P is directly related to the skin temperature of the planet (see Sect. 4.3) whereas τ_lim is the optical depth at which the irradiated picket-fence model differs from the semi-grey model. The steps to use our model are as follow:

• 1)

  choose the pair of variables suitable for the problem: (R, β) or (γ_P, τ_lim), for example;

• 2)

  using Eqs. (87) to (95), calculate the values of γ_P, γ₁, γ₂, and τ_lim;

• 3)

  using Eqs (77) to (81) and Eqs. (66) to (75), calculate the coefficients A, B, C, D, and E.;

• 4)

  using Eq. (76), calculate the temperature/optical depth profile;

• 5)

  using Eq. (4), calculate the pressure/optical depth relationship and therefore the pressure/temperature profile.

For Rosseland opacities depending on the temperature, step 5) can be iterated until convergence. Given the apparent complexity of the solution, we provide a ready-to-use code³ in different languages that gives the temperature/optical depth profile (steps 1 to 4) or the temperature/pressure profile given a Rosseland mean opacity.

The relationship between the different variables are listed below: $\begin{array}{ccc}& & {\mathit{\gamma }}_{\mathrm{P}}\mathrm{=}\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}\mathrm{+}\frac{\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}}{\mathit{R}}\mathrm{-}\frac{{\left(\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}\right)}^{\mathrm{2}}}{\mathit{R}}\\ & & {\mathit{\tau }}_{\lim }\mathrm{=}\frac{\sqrt{\mathit{R}}\sqrt{\mathit{\beta }{\mathrm{(}\mathit{R}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}{\mathrm{(}\mathit{R}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}\mathit{R}}}{\sqrt{\mathrm{3}}{\left(\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}\right)}^{\mathrm{2}}}\\ & & \mathit{R}\mathrm{=}\frac{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\tau }}_{\lim }\mathrm{+}\sqrt{\mathrm{\Delta }}}{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\tau }}_{\lim }\mathrm{-}\sqrt{\mathrm{\Delta }}}\\ & & \mathit{\beta }\mathrm{=}\frac{\sqrt{\mathrm{\Delta }}\mathrm{-}\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\tau }}_{\lim }}{\mathrm{2}\sqrt{\mathrm{\Delta }}}\\ & & \mathrm{\Delta }\mathrm{=}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}\mathrm{+}\mathrm{3}\sqrt{{\mathit{\gamma }}_{\mathrm{P}}}{\mathit{\tau }}_{\lim }\left(\mathrm{2}\sqrt{\mathrm{3}}{\mathit{\gamma }}_{\mathrm{P}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\tau }}_{\lim }\mathrm{-}\mathrm{4}\sqrt{\mathrm{3}}\right)\\ & & {\mathit{\gamma }}_{\mathrm{1}}\mathrm{=}\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}\\ & & {\mathit{\gamma }}_{\mathrm{2}}\mathrm{=}\frac{\mathit{\beta }\mathrm{+}\mathit{R}\mathrm{-}\mathit{\beta R}}{\mathit{R}}\\ & & {\mathit{\gamma }}_{\mathrm{1}}\mathrm{=}\frac{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\tau }}_{\lim }\mathrm{+}\sqrt{\mathrm{\Delta }}}{\mathrm{6}{\mathit{\tau }}_{\lim }}\\ & & {\mathit{\gamma }}_{\mathrm{2}}\mathrm{=}\frac{\sqrt{\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}}\mathrm{+}\mathrm{3}{\mathit{\gamma }}_{\mathrm{P}}{\mathit{\tau }}_{\lim }\mathrm{-}\sqrt{\mathrm{\Delta }}}{\mathrm{6}{\mathit{\tau }}_{\lim }}\\ & & \mathit{R}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{P}}\mathrm{-}\mathrm{1}}{\mathrm{2}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}}\mathrm{+}\sqrt{{\left(\frac{{\mathit{\gamma }}_{\mathrm{P}}\mathrm{-}\mathrm{1}}{\mathrm{2}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}}\right)}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{P}}\mathrm{-}\mathrm{1}}{\mathrm{2}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}\end{array}$

3.7. About averaging

Equation (76) can thus be considered to depend on κ_R, γ_P ≡ κ_P/κ_R, and β. While κ_R can be considered a function of pressure and temperature (e.g. extracted from a known Rosseland opacity table) when deriving the atmospheric temperature profile, it is important to realize that the analytical solution remains valid only if γ_P and β are held constant. This analytical solution therefore cannot accommodate consistent Rosseland and Planck opacities as a function of depth. A solution consisting of atmospheric slices with different values of γ_P have been derived by Chandrasekhar (1935) for the non-irradiated case, but it becomes too complex to be handled easily.

Furthermore, the solution is provided only for one fixed direction of the incoming irradiation. When considering the case of a non-resolved planet around a star, any information acquired on its atmosphere will have been averaged over at least a fraction of its surface. Solving this problem for the particular case of Eq. (76) goes beyond the scope of the present work, but it can be approximated relatively well on the basis of the study by Guillot (2010). This work shows that given an irradiation flux at the substellar point $\mathit{\sigma }{\mathit{T}}_{\mathrm{sub}}^{\mathrm{4}}\mathrm{\equiv }\mathit{\sigma }\mathrm{(}{\mathit{R}}_{\mathrm{\ast }}\mathit{/}\mathit{D}{\mathrm{)}}^{\mathrm{2}}{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}$, where T_∗ is the star’s effective temperature, R_∗ its radius, and D the star-planet distance, the average temperature profile of the planet will be very close to that obtained from the one-dimensional solution with an average angle ${\mathit{\mu }}^{\mathrm{\ast }}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ and an average irradiation temperature T_irr = (1 − A)^1/4f^1/4T_sub, where A is the (assumed) Bond albedo of the atmosphere and f is a correction factor, equal to 1/4 when averaging on the entire surface of the planet and equal to 1/2 when averaging on the dayside only. This corresponds to the so-called isotropic approximation. In the semi-grey case, it is found to be within 2% of the actual average for a typical hot-Jupiter (see Fig. 2 of Guillot 2010).

For the interpretation of spectroscopic and photometric data of secondary eclipses, the dayside average is often used (f = 1/2). For the calculation of evolution models, the global average is the correct physical quantity to be used when the composition and opacity variations in latitude and longitude are not precisely known (see Guillot 2010). In that case, f = 1/4 which is equivalent to setting the irradiation temperature equal to the usual equilibrium temperature defined as T_eq ≡ T_∗(R_∗/2D)^1/2 (Saumon et al. 1996). Obviously detailed interpretations must use an approach mixing three-dimensional dynamical and radiative transfer models (see Guillot 2010; Heng et al. 2012; Showman et al. 2009).

3.8. Adding several bands in the visible

Although, for the simplicity of the derivation, our model used only one spectral band in the visible channel, it can be easily extended to n visible bands. The most important point is that our equations, and in particular Eq. (43), are linear in the visible. Thus, the equations can be solved for any linear combination of visible bands. In that case the first momentum of the visible intensity (see Eq. (40)) writes ${\mathit{J}}_{\mathrm{v}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{H}}_{\mathrm{v}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{{\mathit{\mu }}_{\mathrm{\ast }}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{\beta }}_{\mathrm{v}\mathit{i}}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}\mathit{i}}^{\mathrm{\ast }}\mathit{\tau }}\mathit{,}$(97)where β_vi is the relative spectral extent of the ith band and γ_vi = κ_vi/κ_R with κ_vi the opacity in the ith visible band. Equation (76)then becomes $\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}& \mathrm{=}& \frac{\mathrm{3}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}}{\mathrm{4}}\mathrm{(}\mathit{\tau }\mathrm{+}\mathit{A}\mathrm{+}\mathit{B}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{)}\\ & & \mathrm{+}\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{n}}\frac{\mathrm{3}{\mathit{\beta }}_{\mathrm{v}\mathit{i}}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}}{\mathrm{4}}{\mathit{\mu }}_{\mathrm{\ast }}\mathrm{(}{\mathit{C}}_{\mathit{i}}\mathrm{+}{\mathit{D}}_{\mathit{i}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\tau }}_{\lim }}\mathrm{+}{\mathit{E}}_{\mathit{i}}{\mathrm{e}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{v}\mathit{i}}^{\mathrm{\ast }}\mathit{\tau }}\mathrm{)}\end{array}$(98)where C_i, D_i, and E_i are the coefficients C, D, and E given by Eqs. (79)to (81)where ${\mathit{\gamma }}_{\mathrm{v}}^{\mathrm{\ast }}$ have been replaced by ${\mathit{\gamma }}_{\mathrm{v}\mathit{i}}^{\mathrm{\ast }}$.