Self-luminous and Irradiated Exoplanetary Atmospheres Explored with HELIOS

2.1.1. Definitions

In HELIOS's original design, the stellar shortwave flux and planetary thermal flux were treated equally using the two-stream approximation (Malik et al. 2017). With only the two-stream fluxes, it is not possible to set a stellar irradiation angle. Furthermore, the temperatures on the planetary dayside can only be adjusted globally by the f heat redistribution parameter, which sets the efficiency of heat exchange between the day- and nightside of the planet (Spiegel & Burrows 2010). However, it is not possible to model the local radiative balance for a given latitude and longitude. We account for these two insufficiencies by extending the two-stream equations to incorporate a direct stellar beam in addition to the two-stream flux expressions. This separation of stellar flux and planetary emission has been well known in the planetary science literature (Toon et al. 1989); however, in contrast to earlier studies, we do not differentiate between shortwave and longwave radiation. Our flux expressions hold irrespective of the explored wavelength.

In the following, we define the relevant physical quantities, which deviate slightly from the definitions given in Heng et al. (2018). Nonetheless, the subsequent derivations and final flux expressions are identical.

Let us start with the plane-parallel radiative transfer equation,

$$\begin{eqnarray}&&\mu \displaystyle \frac{{{dI}}^{\mathrm{tot}}}{d\tau }={I}^{\mathrm{tot}}-S,\end{eqnarray} \tag{ 1 }$$

which provides the change of the total intensity I^tot with optical depth τ. The intensity generally depends on the zenith and azimuth angles θ and ϕ, while τ simply measures the vertical change in optical depth, increasing in the downward direction. By convention, θ is measured from the upward-pointing plane normal vector, which makes $\mu =\cos \theta$ negative for downward-pointed and positive for upward-pointed radiation rays. The function S incorporates all radiation sources that add to the intensity I^tot. It is written as (Chandrasekhar 1960; Mihalas 1970, 1978)

$$\begin{eqnarray}&&S=(1-{\omega }_{0})B+\displaystyle \frac{{\omega }_{0}}{4\pi }{\int }_{0}^{4\pi }{ \mathcal P }(\theta ,\phi ;{\theta }^{{\prime} },{\phi }^{{\prime} }){I}^{\mathrm{tot}}({\theta }^{{\prime} },{\phi }^{{\prime} })\,d{{\rm{\Omega }}}^{{\prime} },\end{eqnarray} \tag{ 2 }$$

where the first term is the thermal blackbody emission with the Planck function B and the single-scattering albedo ω₀. The latter describes the relative strength of atmospheric extinction due to scattering only to the total extinction, which is scattering and absorption. The second term considers the incoming rays from the direction (θ', ϕ') and calculates, by multiplying them with the scattering phase function ${ \mathcal P }$, whether they are scattered into the line of sight (θ, ϕ). The integral over all possible incoming angles provides the total addition of scattered intensity to the source function S.

In addition to our previous theoretical excursion in Malik et al. (2017), here we divide the total intensity

$$\begin{eqnarray}&&{I}^{\mathrm{tot}}={I}^{\mathrm{dir}}+{I}^{\mathrm{diff}}\end{eqnarray} \tag{ 3 }$$

into the direct stellar beam I^dir and the diffuse intensity I^diff.⁷ The direct beam changes with the optical depth as

$$\begin{eqnarray}&&{I}^{\mathrm{dir}}={F}_{* ,\mathrm{TOA}}{e}^{\tau /\mu }\delta (\mu -{\mu }_{* })\delta (\phi -{\phi }_{* }),\end{eqnarray} \tag{ 4 }$$

propagating on a straight line through the atmosphere downward along the stellar angles ${\theta }_{* }={\cos }^{-1}{\mu }_{* }$ and ϕ_*.

The incoming stellar flux at the top of the atmosphere (TOA) is given by

$$\begin{eqnarray}&&{F}_{* ,\mathrm{TOA}}={\left(\displaystyle \frac{{R}_{* }}{a}\right)}^{2}{F}_{* },\end{eqnarray} \tag{ 5 }$$

where R_* is the stellar radius, a is the orbital distance, and F_* is the stellar surface flux, which can be approximated by either the blackbody Planck function B for the stellar temperature T_*, i.e., F_* = πB(T_*), or a synthetic model spectrum, e.g., PHOENIX  (Husser et al. 2013) or Kurucz/ATLAS (Kurucz 1979; Murphy & Meiksin 2004; Munari et al. 2005). The flux associated with the direct beam is

$$\begin{eqnarray}\begin{array}{rcl}{F}^{\mathrm{dir}} & \equiv & -{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-1}^{1}\mu {I}^{\mathrm{dir}}\,d\mu \,d\phi ,\\ & = & \ -{\mu }_{* }{F}_{* ,\mathrm{TOA}}{e}^{\tau /{\mu }_{* }},\end{array}\end{eqnarray} \tag{ 6 }$$

also simultaneously leading to the upper boundary condition of ${F}_{\mathrm{TOA},\downarrow }$ = −μ_*F_*,TOA, where τ = 0. The leading minus sign in Equation (6) is needed because we define the radiative flux to always be a positive quantity.

The diffuse up- and downward fluxes are defined as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\uparrow }^{\mathrm{diff}} & \equiv & {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{1}\mu {I}^{\mathrm{diff}}\,d\mu \,d\phi ,\\ {F}_{\downarrow }^{\mathrm{diff}} & \equiv & -{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-1}^{0}\mu {I}^{\mathrm{diff}}\,d\mu \,d\phi ,\end{array}\end{eqnarray} \tag{ 7 }$$

from which we immediately obtain the total up- and downward fluxes, ${F}_{\uparrow }={F}_{\uparrow }^{\mathrm{diff}}$ and ${F}_{\downarrow }={F}_{\downarrow }^{\mathrm{diff}}+{F}^{\mathrm{dir}}$. Note that the planetary thermal emission is included in the diffuse component of the flux.

2.1.2. Interface Flux Expressions

For the in-depth derivation of the diffuse fluxes, we refer the reader to Heng et al. (2018). Here we specifically present and write the diffuse flux expressions as included in HELIOS. From here on, we omit the superscript "diff" and subscript i for certain quantities, e.g., ω₀, for better readability. The fluxes at the ith interface in the staggered numerical grid (see Malik et al. 2017, Figure 2) read

$$\begin{eqnarray}\begin{array}{rcl}{F}_{i,\uparrow } & = & \displaystyle \frac{1}{\chi }\left(\psi {F}_{i-1,\uparrow }-\xi {F}_{i,\downarrow }+2\pi \epsilon {{ \mathcal B }}_{\uparrow }+\displaystyle \frac{1}{{\mu }_{* }}{{ \mathcal I }}_{\uparrow }\right),\\ {F}_{i,\downarrow } & = & \displaystyle \frac{1}{\chi }\left(\psi {F}_{i+1,\downarrow }-\xi {F}_{i,\uparrow }+2\pi \epsilon {{ \mathcal B }}_{\downarrow }+\displaystyle \frac{1}{{\mu }_{* }}{{ \mathcal I }}_{\downarrow })\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal B }}_{\uparrow } & \equiv & (\chi +\xi ){B}_{i}-\xi {B}_{i-1}+\displaystyle \frac{\epsilon }{1-{\omega }_{0}{g}_{0}}(\chi -\xi -\psi ){B}^{{\prime} },\\ {{ \mathcal I }}_{\uparrow } & \equiv & \psi {{ \mathcal G }}_{+}{F}_{i-1}^{\mathrm{dir}}-(\xi {{ \mathcal G }}_{-}+\chi {{ \mathcal G }}_{+}){F}_{i}^{\mathrm{dir}},\\ {{ \mathcal B }}_{\downarrow } & \equiv & (\chi +\xi ){B}_{i}-\psi {B}_{i+1}+\displaystyle \frac{\epsilon }{1-{\omega }_{0}{g}_{0}}(\xi -\chi +\psi ){B}^{{\prime} },\\ {{ \mathcal I }}_{\downarrow } & \equiv & \psi {{ \mathcal G }}_{-}{F}_{i+1}^{\mathrm{dir}}-(\chi {{ \mathcal G }}_{-}+\xi {{ \mathcal G }}_{+}){F}_{i}^{\mathrm{dir}},\end{array}\end{eqnarray} \tag{ 9 }$$

and B_i = B(T_i) is the Planck function evaluated for the temperature at interface i. Further, we use above

$$\begin{eqnarray}\begin{array}{rcl}\chi & \equiv & {\zeta }_{-}^{2}{{ \mathcal T }}^{2}-{\zeta }_{+}^{2},\\ \xi & \equiv & {\zeta }_{+}{\zeta }_{-}(1-{{ \mathcal T }}^{2}),\\ \psi & \equiv & ({\zeta }_{-}^{2}-{\zeta }_{+}^{2}){ \mathcal T },\\ {\zeta }_{\pm } & \equiv & \displaystyle \frac{1}{2}\left[1\pm {\left(\displaystyle \frac{1-{\omega }_{0}}{1-{\omega }_{0}{g}_{0}}\right)}^{1/2}\right],\\ {{ \mathcal G }}_{\pm } & \equiv & \displaystyle \frac{1}{2}\left[{ \mathcal L }\left(\displaystyle \frac{1}{\epsilon }\pm \displaystyle \frac{1}{{\mu }_{* }(1-{\omega }_{0}{g}_{0})}\right)\pm \displaystyle \frac{{\omega }_{0}{g}_{0}{\mu }_{* }}{1-{\omega }_{0}{g}_{0}}\right],\\ { \mathcal L } & \equiv & \displaystyle \frac{(1-{\omega }_{0})(1-{\omega }_{0}{g}_{0})-1}{1/{\mu }_{* }^{2}-1/{\epsilon }^{2}(1-{\omega }_{0})(1-{\omega }_{0}{g}_{0})},\\ {B}^{{\prime} } & \equiv & \displaystyle \frac{{\rm{\Delta }}B}{{\rm{\Delta }}\tau },\end{array}\end{eqnarray} \tag{ 10 }$$

where the last equation in Equation (10) shows our choice of expanding the Planck function linearly in optical depth, with ΔX representing the difference of any quantity X across a layer. Setting B' = 0 and ${B}_{i\pm 1}\to {B}_{i}$ in Equation (9) translates into a model with isothermal layers. Equation (8) is valid for nonisotropic coherent scattering with the direction of scattering given by the asymmetry parameter g₀ ∈ [−1, 1], as defined in, e.g., Goody & Yung (1989) or Pierrehumbert (2010). The transmission function is given as

$$\begin{eqnarray}&&{ \mathcal T }\equiv {e}^{-1/\epsilon \sqrt{(1-{\omega }_{0}{g}_{0})(1-{\omega }_{0})}{\rm{\Delta }}\tau }.\end{eqnarray} \tag{ 11 }$$

We make use of the first Eddington coefficient , which is defined as the assumed (constant) ratio between the first and second moments of the intensity and is the inverse of the commonly known diffusivity parameter ${ \mathcal D }$ (Armstrong 1968; Heng et al. 2014). The latter effectively determines how diffuse the radiation behaves (${ \mathcal D }=2$: fully isotropic radiation; ${ \mathcal D }=1$: fully vertically oriented radiation).

In the case of ω₀ = 1 (pure scattering), Equation (8) needs to be replaced by

$$\begin{eqnarray}\begin{array}{rcl}{F}_{i,\uparrow } & = & {F}_{i-1,\uparrow }+\displaystyle \frac{({F}_{i,\downarrow }-{F}_{i-1,\uparrow })(1-{g}_{0}){\rm{\Delta }}\tau }{(1-{g}_{0}){\rm{\Delta }}\tau +2\epsilon }+{{ \mathcal J }}_{\uparrow },\\ {F}_{i,\downarrow } & = & {F}_{i+1,\downarrow }+\displaystyle \frac{({F}_{i,\uparrow }-{F}_{i+1,\downarrow })(1-{g}_{0}){\rm{\Delta }}\tau }{(1-{g}_{0}){\rm{\Delta }}\tau +2\epsilon }+{{ \mathcal J }}_{\downarrow },\end{array}\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal J }}_{\uparrow } & \equiv & \displaystyle \frac{1}{(1-{g}_{0}){\rm{\Delta }}\tau +2\epsilon }{{ \mathcal K }}_{\uparrow },\\ {{ \mathcal K }}_{\uparrow } & \equiv & \left[{\mu }_{* }+(1-{g}_{0}){\rm{\Delta }}\tau +\epsilon \right]{F}_{i}^{\mathrm{dir}}-({\mu }_{* }+\epsilon ){F}_{i-1}^{\mathrm{dir}},\\ {{ \mathcal J }}_{\downarrow } & \equiv & \displaystyle \frac{1}{(1-{g}_{0}){\rm{\Delta }}\tau +2\epsilon }{{ \mathcal K }}_{\downarrow },\\ {{ \mathcal K }}_{\downarrow } & \equiv & \left[{\mu }_{* }-(1-{g}_{0}){\rm{\Delta }}\tau -\epsilon \right]{F}_{i}^{\mathrm{dir}}-({\mu }_{* }-\epsilon ){F}_{i+1}^{\mathrm{dir}}.\end{array}\end{eqnarray} \tag{ 13 }$$

Although in real atmospheres, the limit ω₀ = 1 is strictly academic, we include those expressions in the code for numerical stability reasons when ${\omega }_{0}\gt 1-| \epsilon |$, with $| \epsilon |$ ≪ 1.

We account for atmospheric convection by adjusting the superadiabatic lapse rates back to the dry adiabatic lapse rate (Manabe et al. 1965; Manabe & Wetherald 1967). In the model grid, we check for each pair of adjacent layers i and i + 1 (the latter being above the former) whether

$$\begin{eqnarray}&&{T}_{i+1}\geqslant {T}_{i}{\left(\displaystyle \frac{{P}_{i+1}}{{P}_{i}}\right)}^{\kappa }\end{eqnarray} \tag{ 14 }$$

is satisfied, where T is the temperature, P is the pressure, and κ is the adiabatic coefficient, for which we use the definition from the planetary climate literature

$$\begin{eqnarray}&&\kappa \equiv {\left(\displaystyle \frac{d\mathrm{ln}T}{d\mathrm{ln}P}\right)}_{S},\end{eqnarray} \tag{ 15 }$$

with S being the entropy. In stellar astrophysics, κ corresponds to ∇_ad ≡ (Γ₂ − 1)/Γ₂ with

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{2}\equiv {\left[1-\displaystyle \frac{P}{{c}_{P}\rho T}\displaystyle \frac{{\chi }_{T}}{{\chi }_{\rho }}\right]}^{-1},\end{eqnarray} \tag{ 16 }$$

where c_P is the specific heat capacity, ρ is the density, ${\chi }_{T}\equiv {(\partial \mathrm{ln}P/\partial \mathrm{ln}T)}_{\rho }$, and ${\chi }_{P}\equiv {(\partial \mathrm{ln}P/\partial \mathrm{ln}\rho )}_{T}$ (Hansen et al. 2004). Invoking Equation (14) is equivalent to testing for Schwarzschild's criterion. If the condition is not satisfied, we correct the temperatures back to the adiabatic lapse rate.

To find the correct adiabat in the convective region, we check for one additional criterion. As convection does not create or eliminate energy, we demand that the enthalpy, which is the total internal energy plus the work done to reach equilibrium with the surrounding medium, must be conserved within a convective region. In our one-dimensional case, the quantity of interest is the product of the specific enthalpy, i.e., enthalpy per unit mass, and m_col is the column mass,

$$\begin{eqnarray}&&{\int }_{\mathrm{conv}.\mathrm{region}}\widetilde{H}\,{{dm}}_{\mathrm{col}}={\int }_{\mathrm{conv}.\mathrm{region}}{\widetilde{H}}^{{\prime} }\,{{dm}}_{\mathrm{col}},\end{eqnarray} \tag{ 17 }$$

where we evaluate the integral over the whole convective region. The prime superscript denotes a quantity after the adjustment. From thermodynamics and the ideal gas law,

$$\begin{eqnarray}&&{dH}={dU}+d({PV})={C}_{V}{dT}+{nRdT},\end{eqnarray} \tag{ 18 }$$

with the volume V; the heat capacities at constant volume and pressure, C_V and C_P, respectively; and the number of gas particles N. With C_V + nR = C_P, we obtain dH = C_PdT. Integrating, dividing by mass, and inserting into Equation (17) leads to

$$\begin{eqnarray}&&{\int }_{\mathrm{conv}.\mathrm{region}}{c}_{P}T\,{dP}={\int }_{\mathrm{conv}.\mathrm{region}}{c}_{P}^{{\prime} }{T}^{{\prime} }\,{dP},\end{eqnarray} \tag{ 19 }$$

with c_P being the specific heat capacity per unit mass. We have further substituted m_col = P/g and used the common assumption of treating g as constant throughout the atmosphere. As the solution T' is following an adiabat, it is convenient to express it via the potential temperature Θ, which is constant along the adiabat, so that

$$\begin{eqnarray}&&{T}^{{\prime} }(P)={\rm{\Theta }}{\left(\displaystyle \frac{{P}_{0}}{P}\right)}^{-\kappa },\end{eqnarray} \tag{ 20 }$$

where we denote P₀ as a reference pressure (e.g., the bottom of the convective region). The potential temperature satisfying Equation (17) is

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{\int }_{\mathrm{bot}}^{\mathrm{top}}{c}_{P}T\,{dP}}{{\int }_{\mathrm{bot}}^{\mathrm{top}}{c}_{P}{\left({P}_{0}/P\right)}^{-\kappa }\,{dP}},\end{eqnarray} \tag{ 21 }$$

where "top" and "bot" mark the top and bottom boundaries of the convective region (Mendonça et al. 2018b).

Numerically, it is necessary to iterate between radiative convergence and convective adjustment. In practice, the following procedure is executed.

• 1.  
  HELIOS iterates until the atmosphere is in radiative equilibrium. See Malik et al. (2017) for more details.
• 2.  
  All layers are checked for unstable lapse rates with Equation (14).
• 3.  
  If unstable layers are found, they are corrected as described by Equations (20) and (21).
• 4.  
  Pairs of layers whose in-between lapse rate has been corrected are now called convective layers.
• 5.  
  One forward step with the radiative iteration is executed. Since this again changes the lapse rates, steps 2–4 need to be repeated.
• 6.  
  In general, steps 2–5 are repeated until a converged solution is found. Such a solution needs to satisfy (i) convective stability, (ii) local radiative equilibrium in all radiative layers, and (iii) global energy equilibrium. Whereas criteria (i) and (ii) are inherent to the convective adjustment mechanism, the satisfaction of the global energy criterion is somewhat more tricky. We discuss our method of how to solve it in Appendix D.

The net convective flux, which is nonzero in the convective layers, is a by-product of the convective adjustment and reads ${{ \mathcal F }}_{\mathrm{conv},-}={{ \mathcal F }}_{-}-{{ \mathcal F }}_{\mathrm{rad},-}$, where ${{ \mathcal F }}_{-}$ is the net flux and ${{ \mathcal F }}_{\mathrm{rad},-}$ is the net radiative flux. We denote a wavelength-independent flux by ${ \mathcal F }$.

In this work, we employ the equation of state of Saumon et al. (1995) to obtain the adiabatic coefficient. A simple ideal gas extension toward lower pressures (P < 10⁻² bar) and temperatures (T ≲ 100 K) was included by Alibert et al. (2005) and Mordasini et al. (2012). Since the species other than hydrogen and helium (the metals) represent a small contribution by mass, we include them only approximately in the calculation of the adiabatic gradient. As in Baraffe et al. (2008), we compute the effective helium fraction Y_eff = Y + Z, with X + Y + Z = X + Y_eff ≡ 1, where X, Y, and Z are the mass fractions of elemental hydrogen, helium, and the metal elements. We then use Y_eff to interpolate within the Saumon et al. (1995) tables.

Figure 1 shows κ in the temperature–pressure (T–P) regime of interest. The low values of κ ∼ 0.04 around 3000 K come from the dissociation of molecular hydrogen, and those near 10⁴ K come from the ionization of atomic hydrogen.

Zoom In Zoom Out Reset image size

To obtain the heat capacity, we use Equation (2.11) in Pierrehumbert (2010), c_P = k_B/(μm_uκ), with the Boltzmann constant k_B, the mean molecular weight μ, and the atomic mass unit m_u. In the last stages of this work, we realized that this expression is valid only in the limit of a perfect gas (i.e., with a constant number of degrees of freedom) or when the number of degrees of freedom changes independently of pressure. Since this is not given at dissociation or ionization of hydrogen, for instance, the heat capacity turns out to be underestimated by a factor of ∼3 in those narrow temperature regions. Explicitly, one can show easily (see Appendix E) that the correct equation is

$$\begin{eqnarray}&&{c}_{{\rm{P}}}=\displaystyle \frac{{k}_{{\rm{B}}}}{{m}_{{\rm{H}}}\kappa }{\left(-\displaystyle \frac{\partial \tilde{S}}{\partial \mathrm{ln}P}\right)}_{T},\end{eqnarray} \tag{ 22 }$$

where $\tilde{S}\equiv S/({k}_{{\rm{B}}}/{m}_{{\rm{u}}})$ is dimensionless. As the Sackur–Tetrode equation reveals and we verified numerically, the point is that the factor ${\left(-\partial \tilde{S}/\partial \mathrm{ln}P\right)}_{T}$ is equal to 1/μ everywhere in the perfect gas portion of the T–P plane, recovering the limit of Pierrehumbert (2010), except where, for instance, the dissociation and ionization of hydrogen and, to a lesser extent, the ionization of helium or metals occur. For our models, we verified that changing c_P does not lead to any noticeable differences in the obtained T–P profile (not shown).

Another addition to HELIOS is the output of the contribution function and the transmission weighting function. They both serve the purpose of finding the photosphere with and without consideration of the atmospheric emission, respectively. Considering a simple two-stream model with isothermal layers and pure absorption, the upward flux at the ith interface reads

$$\begin{eqnarray}&&{F}_{i,\uparrow }={{ \mathcal T }}_{i-1}{F}_{i-1,\uparrow }+2\pi \epsilon {B}_{i-1}(1-{{ \mathcal T }}_{i-1}).\end{eqnarray} \tag{ 23 }$$

Being interested in the TOA flux, we set the index i in the previous expression to the TOA and take all layers below into account. This leads to

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{TOA},\uparrow } & = & {F}_{\mathrm{BOA},\uparrow }\displaystyle \prod _{i=0}^{n-1}{{ \mathcal T }}_{i}\\ & & +\ \displaystyle \sum _{i=0}^{n-1}\left[2\pi \epsilon {B}_{i}(1-{{ \mathcal T }}_{i})\displaystyle \prod _{\mathop{j=i+1}\limits_{i\ne n-1}}^{n-1}{{ \mathcal T }}_{j}\right],\end{array}\end{eqnarray} \tag{ 24 }$$

with n being the number of layers. Using the altitude z as the vertical grid coordinate, Equation (24) is equivalent to

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{TOA}} & = & {F}_{\mathrm{BOA}}{ \mathcal T }({z}_{\mathrm{BOA}},{z}_{\mathrm{TOA}})\\ & & +\ {\displaystyle \int }_{{z}_{\mathrm{BOA}}}^{{z}_{\mathrm{TOA}}}B(z)\displaystyle \frac{d{ \mathcal T }(z,{z}_{\mathrm{TOA}})}{{dz}}\,{dz},\end{array}\end{eqnarray} \tag{ 25 }$$

where ${ \mathcal T }(z,{z}_{\mathrm{TOA}})$ is the total transmission between the atmosphere at altitude z and TOA. In a discrete layer model, the integral term directly translates to the bracket term in Equation (24). This term represents the contribution of a specific layer i to the TOA spectral emission, hence the name contribution function  (Irwin 2009). Physically, it corresponds to the atmospheric location where the transmission function ${ \mathcal T }$ exhibits the steepest gradient, weighted by the layer's emission. Without the consideration of the emission, one obtains the transmission weighting function Ψ, which simply reads

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{i}=(1-{{ \mathcal T }}_{i})\displaystyle \prod _{j=i+1}^{n-2}{{ \mathcal T }}_{j}.\end{eqnarray} \tag{ 26 }$$

Usually, the contribution function is the quantity analyzed to determine the origin of the planetary emission for a given wavelength. The according atmospheric layer is given by the location, where the contribution function peaks.

In the plane-parallel assumption, the stellar photons are assumed to travel along a straight path, which manifests itself in the use of a constant stellar zenith angle θ_* throughout the model grid. However, for large zenith angles,⁸ the stellar beam path length, proportional to $1/\cos ({\theta }_{* })$, is substantially overestimated; i.e., the path length exceeds the real atmospheric extent. If taking the spherical geometry into account, the stellar path length needs to be corrected downward. This is done by adjusting the zenith angle θ_* and, consequently, ${\mu }_{* }=\cos ({\theta }_{* })$, depending on the location in the atmosphere. In Figure 2, the situation is drawn for a vertical column model with three layers, located along the dotted line. First, the zenith angle should depend on the altitude, which means that α₀ < α₁ < α₂ < α₃ in Figure 2 for α₀ = 90° − θ_*. Second, each layer of the model is reached by stellar photons that traveled along a different path through the planetary hemisphere, i.e., ${\alpha }_{1}\ne {\alpha }_{1}^{{\prime} }$, ${\alpha }_{2}\ne {\alpha }_{2}^{{\prime\prime} }$, ${\alpha }_{3}\ne {\alpha }_{3}^{\prime\prime\prime }$, etc. We rewrite and implement Li & Shibata's (2006) Equation (2) as

$$\begin{eqnarray}&&{\mu }_{* ,{ij}}=-\sqrt{1-{\left(\displaystyle \frac{{R}_{\mathrm{pl}}+{z}_{i}}{{R}_{\mathrm{pl}}+{z}_{j}}\right)}^{2}(1-{\mu }_{* })}\forall j\geqslant i,\end{eqnarray} \tag{ 27 }$$

where z is the altitude and R_pl is the planetary radius. For each layer i, one must consider the stellar pathway through all overlying layers j separately with μ_*,ij (see the Appendix in Mendonça et al. 2018a). We set R_pl as the planet's measured white-light radius and associate it with a pressure of 10 bars (as done, e.g., in Kreidberg et al. 2015).⁹ The height difference from layer j to R_pl is given by z_j (and, analogously, for layer i). The quantity μ_* corresponds to the stellar zenith angle at the layer of interest i. As expected, μ_*,ij reduces to μ_*, if i = j.

Zoom In Zoom Out Reset image size

Since the stellar path length correction was envisioned and applied for Earth's atmosphere models, we explore in Section 4.4 the effect of such a correction to radiative transfer calculations for a case of a hot Jupiter and a super-Earth.

Two-stream radiative transfer models possess an accuracy deficit compared to multistream radiative transfer methods. In a meticulous study, Kitzmann et al. (2013) showed that the hemispheric two-stream radiative transfer method (e.g., Heng et al. 2014), as employed in HELIOS, overestimates the scattering and transmission of stellar radiation through cloud decks by up to 20%. More specifically, the two-stream method overestimates the net greenhouse effect of CO₂ clouds, which are strongly scattering in the thermal wavelengths. Thus, using a too-simplistic radiative transfer method may lead to wrong temperature estimates and corresponding atmospheric conditions, e.g., as happened for studies of the early Mars climate (Kitzmann 2016). Yet two-stream models are still prevalent, as more sophisticated radiative transfer methods often prove to be too computationally challenging to be efficiently used in large-scale climate simulations. Hence, there exist a number of tweaks to enhance the scattering behavior of two-stream models, like the δ-Eddington (Joseph et al. 1976) or the two-stream source function methods (Toon et al. 1989). Recently, Heng & Kitzmann (2017) introduced a correction factor E into the two-stream coupling coefficients ζ_± (see Section 2.1.2), which they calibrated to reproduce the results of a 32-stream DISORT code (Hamre et al. 2013). Later, Heng et al. (2018) developed a new and improved formalism consistently embedding the correction factor E into the two-stream equations. This correction factor depends only on the single-scattering albedo ω₀ and the scattering asymmetry parameter g₀. Conveniently, all scattering problems can be reduced to those two quantities, which means that a model equipped with E can treat any atmospheric condition imaginable. The latest version of HELIOS includes Heng et al.'s (2018) version of the correction factor E through the fitting function given by their Equation (31).

We test the two-stream method, including the scattering correction versus the standard method, in Section 4.5.

The atmospheric mixing ratios are obtained with the FastChem code, which reliably calculates the thermochemical equilibrium abundances of around 550 gas phase species for temperatures between 100 and 6000 K (Stock et al. 2018). We use the solar elemental abundances as stated in Table 1 of Asplund et al. (2009) throughout this study, unless otherwise stated. For nonsolar metallicities, we adjust all elements heavier than He, and for nonsolar C/O ratios, we keep O at solar value and adjust C accordingly.

We post-process the equilibrium abundances with the effects of condensation; i.e., we remove species in the parameter space where their gas abundance is expected to decrease drastically due their own condensation or participation in condensate dust species. See Appendix B for more details.

We have massively extended our list of opacity sources since Malik et al. (2017); see Table 1. The new version of our in-house opacity calculator HELIOS-K  (Grimm & Heng 2015) is able to automatically download and calculate the spectral line lists from the ExoMol, HITRAN/HITEMP, and Kurucz online databases, which facilitates the inclusion of new opacities significantly. In this study, we include the main spectroscopically active and atmospheric abundant molecules, as expected in atmospheres with temperatures from a few hundred to several thousand K. We divide our opacities into the main infrared absorbers, metal or other hydrides, metal oxides, Na and K,¹⁰ and H⁻. We use this division in Section 4.2 to investigate the spectroscopic impact of these absorbing species. For the atmospheric model grid, we use all of the opacity sources listed in Table 1. Figure 3 shows the included opacities for one temperature and pressure weighted with their respective equilibrium mixing ratio. The molecular opacities are computed at a resolution of 10⁻² cm⁻¹, the Na and K opacities at 10⁻¹ cm⁻¹, and H⁻ at 10 cm⁻¹. The calculation of the alkali metals is described in detail in Appendix A. For all other molecules, we use a Voigt profile, and, due to the absence of a first-principles theory for pressure broadening, we adopt the community procedure of an ad hoc truncation of the line wings at 100 cm⁻¹ from the line center. For the opacities with data from the ExoMol database, we include the default broadening as provided in their online library. For the HITRAN/HITEMP opacities, we use their self-broadening parameters.

Zoom In Zoom Out Reset image size

Lastly, we are aware that we are missing the opacities of collision-induced absorption (CIA) H2–H2 blueward of 1 μm  (Borysow et al. 2001; Borysow 2002). They play an insignificant role for hot, irradiated atmospheres because other opacity sources present, such as the alkali metals or the metal hydrides and oxides, assume the role of the dominant shortwave absorbers. Furthermore, CIA absorption generally only becomes important in the deep atmosphere, P ≳ 10 bars, as it scales with the pressure squared. However, at these depths, the convective region usually begins. That means that the atmosphere there is optically thick and thus not visible in the planetary emission, and the according temperatures are given by the adiabatic lapse rate and not by radiative equilibrium.

For the iterative runs, i.e., to find the temperature profile in radiative–convective equilibrium, we employ the κ-distribution method with a correlated-κ approximation as described in Malik et al. (2017). We use 300 wavelength bins over a range of 0.33–10⁵ μm. Each bin contains 20 Gaussian points where independent flux calculations are performed. With the addition of strong shortwave absorbers, occasional discontinuities in the converged temperature profile may occur. This is countered by a temperature-smoothing algorithm (see Appendix C). Furthermore, there is now the option to bypass the κ-distribution method and sample the opacities at the discrete flux wavelength values. This approach, commonly called opacity sampling, is used in this study to generate high-resolution spectra (8263 wavelength bins) by post-processing a given temperature profile. To elucidate on our terminology, we regard the difference between "opacity sampling" and "line-by-line" as a matter of sampling resolution. When the number of sampling points greatly exceeds the number of spectral lines considered, the computation qualifies as being line-by-line. When the number of sampling points is comparable to or smaller than the number of lines, as it is in our case, then one is performing opacity sampling.

4.1.1. Exo-REM, petitCODE, and ATMO

In a series of benchmark tests, Baudino et al. (2017) compared the three radiative transfer models Exo-REM  (Baudino et al. 2015), petitCODE  (Mollière et al. 2015), and ATMO  (Amundsen et al. 2014) to identify troublesome model parameters and how these may affect simulated atmospheres. For this purpose, they aligned the three codes in terms of chemistry, opacity line lists, and the treatment of spectral lines and investigated the effect of each of those aspects. As such an intercode comparison goes beyond the scope of this study, we limit our comparison to the self-consistent calculation of the four therein explored planets: the self-luminous GJ 504b and VHS 1256-1257b, the irradiated super-Earth GJ 436b, and the hot Jupiter WASP-12b.¹¹ In order to mimic their setup, we limit ourselves to the following opacities: H₂O, CO₂, CO, CH₄, NH₃, PH₃, the alkali metals Na and K, and CIA from H₂–H₂ and H₂–He. Furthermore, we include Rayleigh scattering for H₂ and H.

Using their planetary parameters and numerical setup, we self-consistently calculate the temperature profiles and emission spectra of the aforementioned planets; see Figure 4, top left and right panels, respectively. The temperature profiles of the self-luminous planets agree very well. Also, for GJ 436b, the temperatures are consistent among all models, although HELIOS leads to slightly warmer deep atmospheric layers (pressure <1 bar). In general, the temperature differences are minuscule, and thus it is not surprising that the corresponding emission spectra for these three planets also match very well. In fact, we are surprised about such a high level of agreement, as there are still differences in the treatment of opacities between HELIOS and the other models. Not only do we use different line lists but, particularly, the calculation of the far wings deviates. They employ a pure Voigt profile for the Na and K resonance lines, whereas our line wings are influenced by Burrows & Volobuyev's (2003) formalism. This discrepancy is particularly visible in the submicron wavelength range, where the far wings of the alkali resonance lines dominate the absorption. Note that the shown spectra have been downgraded to the same resolution to facilitate a comparison.

Zoom In Zoom Out Reset image size

A special case is WASP-12b, because they additionally include TiO and VO as absorbers for this planet. Although we also extend our opacity list with these two molecules, we obtain a substantially different temperature profile and emission spectrum. We cannot reproduce their large temperature inversion and are consequently missing their strong emission features in the optical. We speculate that the difference must be caused by differences in the TiO and VO line lists. Incomplete line lists or erroneous line strengths may both lead to stark differences in the shortwave heating of the upper atmosphere.

In Figure 4 (bottom panels), we show the corresponding vertical abundance profiles for the four planets used in the models. Even though the models are based on the same elemental abundances of Asplund et al. (2009), the exact molecular abundances depend on the choice of solving method, the utilized thermochemical data, and the size of the chemical network considered. Baudino et al. (2017) gave an overview of the chemistry in petitCODE, ATMO, and Exo-REM in their Section 2.2. We hypothesize three main reasons for the discrepancies in the chemical abundances between the models. First, some differences in the abundances stem directly from differences in the temperature profiles. This is best visible for WASP-12b, for which the found temperature profiles differ the most. Second, the removal of gas species due to condensation is modeled differently by each of the models. Critical to this issue is at which temperature condensation of the gas species occurs and whether the species in question might be removed indirectly by the formation of a dust compound. Unlike the other models, Exo-REM includes a "cold trap" by removing species located above the one where its condensation curve is crossed. We describe the method of how we apply condensation in Appendix B. Third, we include ions in our thermochemical network, which means that at higher temperatures, neutral species are depleted due to ionic transitions. Without this depletion, the number of neutral species may be overestimated. This explains, e.g., the discrepant abundances of K for WASP-12b between the models.

4.1.2. COOLTLUSTY and GENESIS

Next, we join in Gandhi & Madhusudhan's (2017) model comparison of their new radiative model GENESIS with the results of Burrows et al. (2008) for HD 189733b. The latter used the COOLTLUSTY model, which has a long-standing history of exoplanet applications. The most detailed description of the vast array of opacity sources included in COOLTLUSTY is given in Sharp & Burrows (2007). Many fewer opacities are included in GENESIS, namely, H₂O, CO₂, CO, CH₄, NH₃, HCN, C₂H₂, Na, K, and CIA H₂–H₂ and H₂–He. We use the latter batch of GENESIS's opacities for this particular comparison. Furthermore, we use their planetary parameters to obtain the same effective dayside temperature and include a stellar blackbody to avoid potential discrepancies due to a mismatched stellar model.

Figure 5 (left panel) shows the synthetic secondary eclipse spectra together with the corresponding temperature profiles in radiative–convective equilibrium according to each model. The shown spectra are downsampled to a common resolution to facilitate a qualitative comparison. In general, HELIOS's spectrum matches the other two well, with the largest difference being somewhat smaller peaks and troughs in the features between 2 and 10 μm. This is a consequence of the smaller temperature difference between the upper and lower atmosphere, which directly translates to the magnitude of spectral features. Also, in our case, the photosphere is located slightly deeper as compared to the other models. Yet the global difference between HELIOS and COOLTLUSTY appears only slightly larger than the difference between GENESIS and COOLTLUSTY. As in the last subsection, we argue that the most likely causes for the discrepancies between the models stem from the opacities. An in-depth comparison of the radiative transfer method, which exceeds the scope of this study, would need to eliminate the following problematic factors that are present here: different sets of employed opacities, different line lists, different spectral line-wing treatments (see Baudino et al. 2017), and differences in the stellar spectrum (which directly affect secondary eclipse spectra). Until those factors are brought to a common denominator, an absolute convergence of the models remains elusive or coincidental.

Zoom In Zoom Out Reset image size

In Figure 5 (right panel), the vertical abundance profiles as obtained in the HELIOS (solid lines) and GENESIS (dashed lines) models are shown. The main near-infrared absorbers are H₂O and CO, and CH₄ in the bottom atmosphere. The relatively abundant alkali metals, Na and K, provide the shortwave absorption and thus play an important role for the extinction of starlight. We attribute the main discrepancies in the abundance profiles between HELIOS and GENESIS to the size of the employed chemical network and differences in the atmospheric temperatures. The latter cause is particularly noticeable in the case of CH₄, as its abundance substantially decreases with temperature T if T ≳ 850 K. The discrepancy in the H₂O abundance may stem from using a slightly different value for the solar elemental abundances of carbon and oxygen, C/O_HELIOS = 0.55 and C/O_GENESIS = 0.5. Also, in the larger chemical network used in HELIOS, the oxygen atoms are divided among more oxygen-bearing species, which leads to a decrease in the H₂O abundance. Lastly, since the equilibrium chemistry formulae in GENESIS do not include the alkali metals, they appear to use the solar elemental abundances for Na and K. This approach returns somewhat smaller values than what we have.

4.1.3. BT-Cond

There are two reference models that stand out when it comes to atmospheres of self-luminous planets or brown dwarfs: the COOLTLUSTY models and the PHOENIX models. As we compared to COOLTLUSTY in last subsection, we turn here toward PHOENIX. With this code, numerous atmospheric grids have been published with different flavors of ingredients over the past 20 yr (e.g., Allard et al. 2001; Baraffe et al. 2003; Allard et al. 2011; Lothringer et al. 2018). The newest and most up-to-date members of this series are the BT-Cond (only gas opacities) and BT-Settl (includes a cloud model) atmospheric grids (Allard et al. 2012a). We choose BT-Cond as a comparison partner, since we do not model clouds. As BT-Cond includes a vast array of opacity sources, we too employ the full array of opacities at our disposal, as listed in Table 1. We take as test objects a colder and a hotter nonirradiated planet with effective temperatures of 800 and 2400 K, respectively. We set the surface gravity to 10⁴ cm s⁻² and compare the self-consistently calculated radiative–convective temperature profiles and the corresponding emission spectra; see Figure 6 (left panel). We further analyze the appearance of convective zones (shown as broader lines). We find that the temperatures in the hotter test case agree rather well with the only deviation in the upper atmosphere, where HELIOS predicts colder temperatures. However, being optically thin, this region plays only a minor role for the spectroscopic appearance of the planet and manifests itself as slightly weaker emission at wavelengths larger than 1.3 μm. The deep convective zone is almost identical for both models, which is the relevant aspect for the coupling to interior models. The colder planet exhibits a somewhat warmer photosphere and colder deeper layers (pressure <1 bar) in the HELIOS version. Interestingly, the BT-Cond model exhibits a significantly larger temperature gradient in the deep layers, which would be unstable and corrected with our employed equation of state. We attribute this discrepancy to the calculation of atmospheric entropy. On the other hand, the upper atmospheric temperatures are consistent, and both models even agree on the detached convective zone, although it is larger in the HELIOS model. Important for the coupling with interior models is the fact that the deep convective zones are somewhat shifted. According to the lower deep temperature in the HELIOS model, the spectral emission is somewhat weaker below 1.3 μm than in the BT-Cond model. Interestingly, we see strong CrH absorption features at 0.9 and 1 μm that are missing in the BT-Cond model. Indeed, observing the vertical abundance profiles of the two models (Figure 6, top right panel), the abundance of CrH is strongly muted for BT-Cond compared to HELIOS. We attribute this difference to a condensation effect implemented in BT-Cond but missing in HELIOS. Another interesting case is SiO. Namely, the temperatures of the upper atmosphere in the two models are located below (HELIOS) and above (BT-Cond) the condensation temperature of MgSiO₃, responsible for depleting SiO. Also, other species, like NH₃ and H₂S, show considerable differences in their abundance between the two models, which could be caused by differences in either the temperatures and/or thermochemical data and the chemical network utilized.

Zoom In Zoom Out Reset image size

We now present a few numerical tests to shed light on the impact of some of the new features included in HELIOS. In this section, we investigate the effects of different opacity sources on converged atmospheric temperature profiles after running the iteration for radiative–convective equilibrium.

For these tests, we use different opacity tables for our calculations. The tables consist of varying sets of opacities that are premixed using the same chemical network in order to single out the effect of missing opacities while keeping the molecular abundances the same.

• 1.  
  Sample A includes all opacities at our disposal, as listed in Table 1.
• 2.  
  Sample B is equal to sample A minus the H⁻ continuum opacity.
• 3.  
  Sample C is equal to sample B minus the metal and other hydrides and metal oxides, as shown in the bottom left and right panels of Figure 3. Effectively, that means we are removing SiH, CaH, MgH, NaH, AlH, CrH, TiO, VO, AlO, SiO, and CaO.
• 4.  
  Sample D is equal to sample C minus the Na and K opacities. This removes the last strong shortwave absorbers from the opacity pool.

The first setup consists of two self-luminous planets with effective temperatures of 1000 and 2500 K and a surface gravity of 10⁴ cm s⁻²; see Figure 7. For the colder planet, samples A and B lead to the identical temperature profile. Sample C leads to somewhat cooler temperatures, and sample D is markedly different. In the latter case, there are no shortwave absorbers left, so the hot bottom atmosphere cools too efficiently to be convectively unstable. Conversely, including all available opacities leads to the hottest temperatures, since having more absorbing species increases the optical depth and moves the photosphere upward. Subsequently, the lower atmosphere emits radiation less efficiently, which drives the temperatures up. The hotter planet shows similar trends, even more pronounced. Here the addition of H⁻ has a noticeable effect on the temperature. The hotter planet exhibits convective zones extending to lower pressures compared to the colder planet, which is likely due to more spectral lines being active, increasing the overall optical depth. As discussed in Section 1, the entropy found in the deep convective zones is of special interest for the coupling with interior models. Hence, we list the obtained values for the entropies S in the conducted runs. For the colder planet, the entropy in the model with sample A is S_A = S_B = 8.95 and S_C = 8.74. There is no deep convective zone when using sample D. The hotter planet models lead to S_A = 11.23, S_B = 10.99, S_C = 10.53, and S_D = 9.46, with a modeling accuracy of ΔS ≈ 0.02. The values for S are given in units of the Boltzmann constant.

Zoom In Zoom Out Reset image size

The second setup consists of the following irradiated planets covering a broad range of effective dayside temperatures: GJ 1214b (T_eff = 721 K), HD 189733b (T_eff = 1465 K), WASP-12b (T_eff = 3013 K), and KELT-9b (T_eff = 4528 K). We also include an internal temperature of 85 K. For each of these planets, we explore the diversity of dayside temperatures obtained when using the various samples of opacity sources. Our findings are shown in Figure 8. Two global statements can be made immediately. First, the temperature profiles vary substantially depending on the employed opacities. Second, this variance increases with stellar irradiation and effective planetary temperature. For GJ 1214b and HD 189733b, samples A and B lead to the same temperatures. Samples C and D lead to cooler upper and warmer lower atmospheres, resulting in a larger atmospheric greenhouse effect. Sample D additionally moves the stellar photosphere downward due to the lack of shortwave absorbers. With WASP-12b and KELT-9b, we move to effective temperatures beyond 3000 K, where H⁻ opacity starts to have a strong effect. Because H⁻ introduces a very strong continuum opacity throughout all wavelengths (see Figure 3), it makes the temperature profile more isothermal. In addition, the metal hydrides and oxides play a strong role in this regime, leading to minor (WASP-12b) or major (KELT-9b) temperature inversions. In the case of KELT-9b, because the stellar irradiation is so strong, even Na and K suffice to lead to two separate temperature inversions. In fact, we expect nonmonotonic temperature profiles with multiple inversions to be naturally occurring in very hot atmospheres (see Section 4.3).

Zoom In Zoom Out Reset image size

In conclusion, we find that missing opacities in models may lead to stark temperature discrepancies. We thus recommend employing a list of opacity sources as complete as possible to avoid erroneous temperature predictions, an error that directly propagates into the strength of the spectroscopic emission features. More specifically, we find the inclusion of metal hydrides, oxides, and H⁻ imperative to obtaining accurate temperatures. We remark that in very hot atmospheres, atomic or ion species like Fe, Mg, Al or Fe⁺, Mg⁺, Ca⁺ are expected to play a significant role as well, due to their high abundance and vast number of spectral lines in the optical and near-infrared (Kitzmann et al. 2018). We further remark that Fe, Fe⁺, and Ti⁺ were recently observationally confirmed in the atmosphere of the ultrahot Jupiter KELT-9b (Hoeijmakers et al. 2018), reinforcing the notion of the prevalent existence of metals in very hot atmospheres.

We will investigate the impact of more atoms and ions on self-consistent radiative transfer calculations in a future study.

As seen in the previous section, we find that for hotter planets with effective temperatures over 2000 K, nonmonotonic temperature profiles are not the exception but the norm. In this regime, sufficiently high temperatures activate a multitude of rovibronic transitions that generate billions of spectral lines belonging to many different species that cause complex extinction patterns for both stellar and planetary radiation. Even more, with planets orbiting cool M dwarfs or ultrahot Jupiters exhibiting temperatures like K stars, the wavelengths of incoming irradiation and local thermal emission may heavily overlap, eroding the classical shortwave and longwave paradigm. In this sense, we feel that the question of whether a planetary atmosphere exhibits a temperature inversion or not (e.g., Fortney et al. 2008) is ill-posed, as there may be small or large and one or several inversions appearing in higher or lower altitudes.

We conduct a numerical test to explore the origins of nonmonotonic temperature profiles. As a test setup, we choose to model the dayside of KELT-9b and employ our full opacity list without H⁻. As seen in Figure 8, this leads to a temperature profile with two large temperature inversions. In order to better explore the wavelength-dependent opacity effect, we use, for this test only, opacities calculated at a fixed pressure and temperature of P = 1 bar and T = 4500 K. With this simplification, the photospheric pressure at which the optical depth τ equals unity is given directly by

$$\begin{eqnarray}&&P=g/\kappa ,\end{eqnarray} \tag{ 28 }$$

where g is the surface gravity and κ is the wavelength-dependent opacity function. For the following demonstration, it is useful to know that the stellar irradiation of KELT-9 and the planetary thermal emission peak at 0.3 and 0.7 μm, respectively.

The converged temperature profile is shown in the left panel and the contribution function in the right panel of Figure 9. The latter is overlaid by the location of the photosphere, given by Equation (28). Due to the use of the κ-distribution method, each wavelength bin contains a range of unordered opacities that returns a range of photospheric pressures. Now, we ask the question of which opacity bands and spectral lines are responsible for which features of the temperature profile. To this end, we manually switch off or strongly decrease opacities in certain wavelengths and study the effect on the temperature profiles. We elucidate our findings by focusing on individual temperature regions, which are connected by lines A–E to the pressure layers and opacity bands they are most affected by.

• 1.  
  A. The temperature inversion above¹² 10⁻⁶ bars is caused by the Na resonance doublet. Stellar radiation is absorbed by the line center, resulting in an increase of temperature until the thermal radiation emits sufficiently at 0.59 μm to counter the stellar energy deposition.
• 2.  
  B. The dip in temperatures around 10⁻⁶ bars is mainly caused by efficient cooling through the SiO band at 9 μm, with the band at 4.5 μm playing a minor part. The layers at this height cool down until the cooling is countered by stellar heating through the other K and Na line peaks at ≲1 μm.
• 3.  
  C. The inversion around and below 1 mbar is due to the absorption of stellar radiation by strong shortwave bands of SiH and AlH between 0.4 and 0.5 μm. This region represents the stellar photon deposition depth, as the bulk of the stellar flux is absorbed here. The temperatures here increase until the atmosphere can efficiently cool through the SiH and AlH bands as well.
• 4.  
  D. At 1–10 bars lies the majority of the planetary photosphere. The temperatures are relatively low, as the gas is efficient at cooling, and only a small fraction of the stellar radiation is left at these depths. The maxima of the contribution function show that the bulk of the thermal emission originates from these pressures at wavelengths between 0.5 and 2 μm.
• 5.  
  E. Below 100 bars, the atmosphere is optically thick to both the thermal and the stellar flux. Radiative equilibrium under these isotropic conditions manifests through an isothermal temperature profile.

Zoom In Zoom Out Reset image size

In conclusion, we find that multiple temperature inversions are caused by a combination of the strong spectral bands and lines of metal oxides, metal hydrides, and the alkali metals. If these species are present in sufficient atmospheric abundances, as is expected for ultrahot planets such as KELT-9b, we expect complex temperature profiles with multiple inversions to be the common state.

However, we note once more that in this test, we used opacities calculated at one fixed value for the pressure and temperature, whereas in reality, the opacities would change throughout the atmosphere. Finally, we have neglected the H⁻ continuum opacity, which makes the temperatures substantially more isothermal than presented here (see Figure 8).

As stated in Section 2.4, the usual practice of one-dimensional models is to use a plane-parallel grid instead of the accurate spherical geometry of planets. In doing so, the path length of photons in the atmosphere and, consequently, the optical depth are overestimated. This error can be avoided with a geometrical correction that shortens the actual stellar path length by adjusting the zenith angle depending on the location in the atmosphere. In the following, we test our implementation of the geometric correction to the stellar path for its influence on modeling exoplanetary atmospheres. We employ our full opacity list and generate temperatures in radiative–convective equilibrium and the corresponding mock spectra for the hot Jupiter HD 189733b and the super-Earth GJ 1214b. We pick each of those two as poster children for the two prevalent classes of currently characterized exoplanets. We treat the atmosphere of GJ 1214b as hydrogen-dominated, so the effects presented here should be understood as an upper limit on the significance of the path length correction. Heavier atmospheres possess a smaller scale height, diminishing the effect of a path length correction.

As a first test, we use a precalculated atmospheric temperature profile and let the direct stellar beam propagate downward, once with and once without path length correction. We show in Figure 10 the ratio of the corrected to the uncorrected stellar fluxes with increasing pressure (or atmospheric depth). High up in the atmosphere, the ratio is very close to unity but increases as the stellar photons propagate deeper. As the effect of the geometric correction is relative to the path length, the discrepancy increases with pressure. Also, the error in the uncorrected path length increases with the size of the zenith angle. The region of interest is given by the shortwave photosphere (optical depth τ ∼ 1), where the stellar photons are deposited and contribute to the atmospheric heating. The layers below the photosphere are left with a vanishingly small stellar flux, making an error irrelevant. In Figure 10, the photospheres of HD 189733b and GJ 1214b are shown at the wavelength of the respective star's peak emission. At these depths, the ratio of fluxes for HD 18973b and GJ 1214b is ∼2 and ∼3 for zenith angles of 80° and ∼1.1 and ∼1.3 for zenith angles of 70°, respectively. We find the effect of the path length correction to be larger for GJ 1214b than for HD 198773b. This is because the path length error, when ignoring the correct geometry, scales with the relative width of the atmosphere, which is given by the ratio of the atmospheric scale height over the planetary radius H/R_pl. This ratio is larger for GJ 1214b, assuming a hydrogen-dominated atmosphere, than it is for HD 189733b. In fact, we estimate that (H/R_pl)_GJ1214b ∼ 0.02, whereas (H/R_pl)_HD189733b ∼ 0.003. For comparison, the Earth exhibits (H/R_pl)_Earth ∼ 0.001, meaning that the Earth's atmosphere is geometrically thinner compared to higher-mass planets.

Zoom In Zoom Out Reset image size

As a second test we run the fully self-consistent simulation with and without path length correction. We find that the impact of the correction on the temperature profile is negligible for zenith angles up to 70°; see Figure 11 (left panels). For zenith angles θ > 80°, the effect becomes noticeable, with the uncorrected models overestimating the upper atmospheric heating (pressure ≲0.1 bar). Consequently, due to the larger attenuation of the stellar flux above, the lower atmospheric layers (pressure ≳0.1 bar) are cooler in the uncorrected models. As in the previous test, the effect is significantly larger for GJ 1214b than for HD 189733b. In Figure 11 (right panels), we show the effect of the path length correction on the resulting emission spectra. Depicted is the wavelength-dependent percentage error in the outgoing TOA flux of the uncorrected versus the corrected model for different zenith angles, color-coded analogously to the left panels. For HD 189733b, the errors in the spectrum are exceeding the percent level for zenith angles higher than 80°. For GJ 1214b, even θ = 70° leads to errors of several percent for some wavelengths. For both planets, very close at the planetary limb (θ = 89°), the error in the flux is close to or exceeds 100%. It may be surprising that even though the temperatures in the different models are very similar, the error in the flux is significant for high zenith angles. Since the emitted flux depends on the temperature as T⁴, a small relative error in T translates to roughly four times the error in the flux.

Zoom In Zoom Out Reset image size

Based on these tests, we conclude that in most cases, the path length correction should play a secondary role for zenith angles up to around 70°. However, if investigating planetary limb regions with zenith angles around 80° or higher, we recommend using the path length correction for the accuracy of both the modeled temperatures and the spectra. If the atmosphere is expected to be very extended, it may be reasonable to include the correction for zenith angles from 70° upward.

We explore the effects of the scattering correction, which we have implemented into the hemispheric two-stream formalism of HELIOS (see Section 2.5). We self-consistently calculate the dayside temperatures and the corresponding spectral emission of HD 198733b, once with and once without the scattering correction. As we do not include clouds in this work, the only present scattering originates from hydrogen Rayleigh scattering. In this case, we would expect isotropic scattering, where g₀ = 0. However, to increase the tested parameter range, we artificially vary g₀ between 0 and 1 to account for different kinds of scattering particles. The value of ω₀ is given by

$$\begin{eqnarray}&&{\omega }_{0}=\displaystyle \frac{{\sigma }_{{\rm{H}},{{\rm{H}}}_{2}}}{{\sigma }_{{\rm{H}},{{\rm{H}}}_{2}}+{\sigma }_{\mathrm{mol}}},\end{eqnarray} \tag{ 29 }$$

where ${\sigma }_{{\rm{H}},{{\rm{H}}}_{2}}$ is the Rayleigh scattering cross section due to H and H₂, and σ_mol is the cross section due to molecular absorption.

In the first test, we use our full opacity list and find that the atmospheric temperatures do not change noticeably. Including the scattering correction, the deeper layers (pressure <1 bar) merely warm about 1–2 K compared to the nominal case, regardless of the used g₀ (not shown). However, turning our attention to the TOA emitted flux, we find that the correction decreases the backscattered light in the optical wavelengths by about 10%–30%; see Figure 12. For longer wavelengths, the discrepancy drops below 1%, as here the outgoing flux is dominated by thermal emission. In general, we also observe that the discrepancy, i.e., the error in the uncorrected model, becomes smaller as g₀ approaches unity. This may appear counterintuitive, as one could think that two-stream methods fare more poorly for larger aerosols, which exhibit forward scattering (g₀ ≈ 1). However, our results are consistent with those of Heng et al. (2018), well visualized in their Figure 1, who showed that the scattering correction vanishes, i.e., the factor $E\to 1$, as g₀ goes to unity. Apart from this, we find our reduction of scattered flux consistent with Kitzmann et al.'s (2013) result that two-stream methods overpredict scattering by up to 20%.

Zoom In Zoom Out Reset image size

In a second test, we include only the main infrared absorbers, shown in Figure 3 (top right panel). In this way, it is the Rayleigh scattering that becomes responsible for the bulk of the stellar flux extinction. This significantly increases the value of w₀ in the optical wavelengths (see Equation (29)), making the impact of the scattering correction larger. In this scenario, the bottom atmospheric temperatures diverge between 50 and 200 K depending on the g₀ parameter, and the corrected model is again warmer in the bottom (not shown). Hence, tuning the strength of the scattering directly affects the extent of the greenhouse warming. The thermal emission is subsequently more discrepant, caused by the larger temperature differences between the models. As can be seen in Figure 12, the error exceeds several percent in the near-infrared. The scattered flux discrepancy in the optical is with values around 10% on the same order of magnitude as in the previous test.

We conclude that the scattering correction provides a significant improvement for the accuracy of the stellar flux scattering. If scattering dominates over the absorption at shorter wavelengths, any potential errors on the stellar flux propagate to the temperatures and thermal emission.

In addition to the model comparisons and numerical tests presented, we also include self-consistently calculated cloud-free atmospheric models using our improved radiative transfer code HELIOS. We include the full opacity list, as given in Table 1, and the scattering correction, as described in Section 2.5. Each model output provides, among others, the atmospheric temperature profile in radiative–convective equilibrium and a post-processed high-resolution emission spectrum, as described in Section 3.2. Additionally, it contains information about any other calculated quantity, like the entropy, optical depth, contribution function, etc.

Below, we first describe the parameters of the grid for the self-luminous planets, and then we elaborate on the calculated sample of irradiated planets. Finally, we provide information on how to download the models.

4.6.1. Grid of Self-luminous Planets

The calculated atmosphere grids for self-luminous planets span the following dimensions. The effective temperature range is T_eff = [200, 3000 K] with a step size ΔT_eff = 100 K. The surface gravity, ${\mathrm{log}}_{10}g$ = [1.4, 5.0], with ${\rm{\Delta }}{\mathrm{log}}_{10}g$ = 0.2, [g] = cm s⁻². The metallicity range is [M/H] = [−1, 1] with Δ[M/H] = 0.5, where we write

$$\begin{eqnarray}&&\left[\displaystyle \frac{M}{{\rm{H}}}\right]={\mathrm{log}}_{10}\left(\displaystyle \frac{{n}_{M}}{{n}_{{\rm{H}}}}\right)-{\mathrm{log}}_{10}\left(\displaystyle \frac{{n}_{M,\odot }}{{n}_{{\rm{H}},\odot }}\right),\end{eqnarray} \tag{ 30 }$$

with n_X being the number density of element X and n_X,⊙ being the solar photospheric value (Asplund et al. 2009). We write M for every element heavier than helium. For each metallicity, we further calculate models with three values for the C/O ratio: 0.1, solar (0.55), and 1. The temperatures are calculated until a pressure of 1000 bars, so that every model atmosphere couples to the deep convective zone. The limits of the grid are chosen to accommodate low-mass young planets and brown dwarfs and will be extended in future to include older, high-mass brown dwarfs with ${\mathrm{log}}_{10}g\gt 5$.

Figure 13 shows an excerpt of the grid for self-luminous planets for various effective temperatures using ${\mathrm{log}}_{10}g$ = 4.0. Displayed are the temperatures in radiative–convective equilibrium and the corresponding TOA emission spectra for some of the models. Note that the spectra shown are downsampled from the native resolution for clarity.

Zoom In Zoom Out Reset image size

Figure 14 (left panel) shows the temperature at a pressure of 1000 bars, T₁₀₀₀, for the model grid calculated. As T₁₀₀₀ lies within the deep convective zone for all calculated models, it can be used directly to extrapolate to the interior adiabat (see Burrows et al. 1997). Otherwise, we also provide the entropy for direct coupling; see Figure 14 (right panel).

Zoom In Zoom Out Reset image size

In addition to the spectral fluxes, we provide a table with the expected magnitudes at the planetary surface for the NACO J, H, Ks, and L' filter bandpasses; see Figure 15 for an excerpt. They are calculated by convolving the filter profiles with the spectral flux. With a Vega spectrum,¹³ the zero-points were obtained. For more details, see Linder et al. (2019).

Zoom In Zoom Out Reset image size

4.6.2. Sample of Irradiated Planets

In addition to the self-luminous planets, we provide atmospheric models for a sample of irradiated planets of interest. Table 2 lists the planets currently included in our model library, together with the system parameters used. For each planet, we calculate 3 × 3 × 2 = 18 models, where we vary the metallicity, [M/H] = −1, 0, 1, C/O = 0.1, solar (0.55), 1, and the heat redistribution factor f = 1/2, 2/3. The latter parameter accounts for the heat redistribution efficiency from the dayside to the nightside (Burrows et al. 2008; Spiegel & Burrows 2010). Choosing f = 1/2 is computationally equivalent to setting the zenith angle to 60°. For these models, we conservatively use a blackbody for the stellar emission and an internal flux consistent with 75 K thermal emission.

Table 2.  System Parameters Used in This Study

Name	Rpl	${\mathrm{log}}_{10}$ g	a	R*	T*
GJ 1132ba	0.1058	3.068	0.0153	0.207	3270
GJ 436bb	0.38	3.106	0.04085	0.464	3684
HD 189733bc	1.216	3.29	0.03142	0.805	5050
Kepler-7bd	1.614	2.62	0.06246	1.966	5933
KELT-9be	1.891	3.30	0.03462	2.362	9560
WASP-8bf	1.038	3.74	0.0801	0.945	5600
WASP-12bg	1.776	3.066	0.02293	1.595	6300
WASP-14bh	1.281	4.01	0.036	1.306	6475
WASP-18bi	1.165	4.281	0.0247	1.23	6400
WASP-19bj	1.386	3.143	0.01655	0.99	5500
WASP-33bk	1.679	3.46	0.0259	1.509	7430
WASP-43bl	1.036	3.672	0.0152	0.667	4520

Notes. R_pl (R_Jup): planetary radius; g (cm s⁻²): surface gravity; a (au): orbital distance; R_* (R_⊙): stellar radius; T_* (K): stellar temperature.

^aBerta-Thompson et al. (2015); Southworth et al. (2017); Bonfils et al. (2018). ^bButler et al. (2004). ^cSouthworth (2010); de Kok et al. (2013); Boyajian et al. (2015). ^dLatham et al. (2010); Demory et al. (2011). ^eGaudi et al. (2017). ^fQueloz et al. (2010). ^gHebb et al. (2009); Chan et al. (2011). ^hJoshi et al. (2009). ⁱSouthworth et al. (2009); Hellier et al. (2009). ^jHebb et al. (2010); Hellier et al. (2011). ^kCollier Cameron et al. (2010); Kovács et al. (2013); Lehmann et al. (2015). ^lGillon et al. (2012).

Download table as:  ASCIITypeset image

Figure 16 displays the dayside temperatures and secondary eclipse spectra of the calculated planets for f = 1/2, solar metallicity, and three different C/O ratios. The planets are listed with increasing effective temperature.

Zoom In Zoom Out Reset image size

4.6.3. Download

The four resonance lines alone, i.e., the D doublet for Na at 0.589 μm and the doublet for K i at 0.77 μm, account for very strong absorption in the optical and are among the most detectable spectral lines in exoplanets and brown dwarfs (Seager & Sasselov 2000; Burrows et al. 2002a; Charbonneau et al. 2002). For modeling purposes, they pose a substantial challenge, as their line wings deviate significantly from a Lorentzian profile. Attempts to characterize the resonance line wings are motivated by theoretical quantum chemistry and enhanced with laboratory measurements (Burrows et al. 2000; Burrows & Volobuyev 2003; Allard et al. 2012b, 2016). We model the far wings as a composite of the Voigt profile in the line core, stitched together with a damping prescription for the far wings as described in Burrows et al. (2000).

While we are aware of an ongoing debate between Burrows et al. and Allard et al. regarding the correct shape of the sodium and potassium profiles in the line wings, we have chosen to adopt the Burrows et al. (2000) and Burrows & Volobuyev (2003) formulation, as it is available in analytical form and straightforwardly implementable. Since the current study is focused on the methodology of HELIOS, rather than the interpretation of measured spectra, we deem it reasonable to defer this issue to a future study.

The line profiles in their core follow a Voigt profile with the collisional broadening half-widths at half-maximum (HWHMs) Γ calculated from impact theory. The used values are Γ_Na = 0.071(T/2000 K)^−0.7 cm⁻¹ atm⁻¹ for sodium and Γ_K = 0.14(T/2000 K)^−0.7 cm⁻¹ atm⁻¹ for potassium. Using statistical theory, the outer line wings are modeled by a power law truncated with an exponential Boltzmann factor. The transition between those two regimes, called the location of detuning, happens at the frequency shift δν from the line center. For sodium, we use δν_Na = 30(T/500 K)^0.6 cm⁻¹, and for potassium, we use δν_K = 20(T/500 K)^0.6 cm⁻¹, where T is the temperature (Unsold 1955; Nefedov et al. 1999; Burrows et al. 2000; Iro et al. 2005). The total line profile at frequency ν is given by

$$\begin{eqnarray}\begin{array}{l}{\rm{\Phi }}(\nu )\left\{\begin{array}{cc}{{\rm{\Phi }}}_{{\rm{V}}}(\nu -{\nu }_{0}) & \mathrm{if}\ | \nu -{\nu }_{0}| \leqslant \delta \nu ,\\ {{\rm{\Phi }}}_{{\rm{V}}}(\delta \nu ){\left(\displaystyle \frac{\delta \nu }{\nu -{\nu }_{0}}\right)}^{3/2}\exp \left[-\displaystyle \frac{h(\nu -{\nu }_{0})}{{k}_{{\rm{B}}}T}\right] & \mathrm{if}\ | \nu -{\nu }_{0}| \gt \delta \nu ,\end{array}\right.\end{array}\end{eqnarray} \tag{ 31 }$$

where Φ_V is the Voigt profile, ν₀ is the line center frequency, δν is the distance from the detuning points to the line center, and h and k_B are the Planck and Boltzmann constants, respectively.

The individual impact of the weaker lines is smaller, and their far wings play a subdominant role. Hence, they are sufficiently modeled with a classical Voigt profile.

In the following, we briefly outline the expressions and parameters used in our nominal method to calculate atomic lines. The frequency ν dependent cross section is given by σ_ν = S · Φ_ν, where the line strength S and the profile Φ_ν are calculated independently. The line strength can be expressed as

$$\begin{eqnarray}&&S=\displaystyle \frac{\pi {e}^{2}{g}_{i}{f}_{{ij}}}{{m}_{e}c}\displaystyle \frac{{e}^{-{E}_{i}/{k}_{{\rm{B}}}T}}{Q(T)}\left(1-{e}^{-{\rm{\Delta }}E/{k}_{{\rm{B}}}T}\right),\end{eqnarray} \tag{ 32 }$$

where i and j denote the lower and upper quantum states, e is the electron charge, g is the statistical weight, f is the oscillator strength, m_e is the electron mass, c is the speed of light, E is the energy, Q is the partition function, and the energy difference ΔE = E_j − E_i. Also, ΔE = hν, with ν being the frequency of the absorbed or emitted photon.

The partition function is defined as

$$\begin{eqnarray}&&Q(T)=\displaystyle \sum _{k}{q}_{k}{e}^{-{E}_{k}/{k}_{{\rm{B}}}T},\end{eqnarray} \tag{ 33 }$$

where the sum leads over all allowed quantum states k. It can be either calculated directly from the energy levels, approximated by fitting functions, or read from precalculated tables.

As the ExoMol or HITRAN databases do not feature atomic species, we have taken the necessary spectral data from the Kurucz line list.¹⁴ With g_i, f_ij, E_i, and E_j provided in the database, the line strength is fully determined. Note that it is common to tabulate ${\mathrm{log}}_{10}({g}_{i}{f}_{{ij}})$ instead of the individual g_i and f_ij quantities.

The Voigt profile Φ_V at frequency ν is defined as the convolution between the Doppler Φ_D and Lorentz profiles Φ_L. Through

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{V}}}(\nu ;{{\rm{\Gamma }}}_{{\rm{D}}},{{\rm{\Gamma }}}_{{\rm{L}}})={\int }_{-\infty }^{\infty }{{\rm{\Phi }}}_{{\rm{L}}}(\nu -\nu ^{\prime} ;{{\rm{\Gamma }}}_{{\rm{L}}}){{\rm{\Phi }}}_{{\rm{D}}}(\nu ^{\prime} ;{{\rm{\Gamma }}}_{{\rm{D}}})\,d\nu ^{\prime} ,\end{eqnarray} \tag{ 34 }$$

the Voigt profile depends on the Doppler and Lorentz HWHMs Γ_D and Γ_L. For thermal velocities, the Doppler HWHM is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{D}}}=\displaystyle \frac{\nu }{c}\sqrt{\displaystyle \frac{2\mathrm{ln}2{k}_{{\rm{B}}}T}{m}},\end{eqnarray} \tag{ 35 }$$

where m is the particle mass. The Lorentz HWHM is a combination of natural, collisional, Stark, and Van der Waals broadening, i.e., Γ_L = Γ_nat + Γ_coll + Γ_Stark + Γ_VdW. Natural broadening arises due to the limited lifetime of an excited quantum state and the subsequent uncertainty on the energy. Using the Kurucz database, we write

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{nat}}=\displaystyle \frac{{10}^{{{\rm{\Gamma }}}_{{\rm{R}}}}}{2\pi },\end{eqnarray} \tag{ 36 }$$

where Γ_R is called the radiative damping coefficient and is a tabulated quantity. Collisional broadening is harder to quantify, as it depends on collisions with the surrounding gas and hence is a function of pressure and temperature. At the order-of-magnitude level, we estimate the HWHM due to collisions as (Heng et al. 2015)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{coll}}=P{\sigma }_{{{\rm{H}}}_{2}}\sqrt{\displaystyle \frac{2}{{{mk}}_{{\rm{B}}}T}},\end{eqnarray} \tag{ 37 }$$

where P is the pressure and ${\sigma }_{{{\rm{H}}}_{2}}\sim {10}^{-15}$ cm² is the cross section of molecular hydrogen, which we assume as the dominating collision partner. In this study, we neglect the remaining two broadening mechanisms due to Stark line splitting and Van der Waals perturbations. We refer to Sharp & Burrows (2007) for an excellent discussion of line profiles and the calculation of opacities.

On a final note, there is no closed analytical solution to calculate the Voigt integral (Equation (34)). However, the problem can be recast in the form of the real part of the complex Faddeeva function with sufficient accuracy (Zaghloul & Ali 2012). The Faddeeva function is included, e.g., in the scipy package for Python.