# Assessment of a Physics-based Retrieval of Exoplanet Atmospheric
Temperatures from Infrared Emission Spectra

Franz Schreier,1 J. Lee Grenfell,2 Fabian Wunderlich,2 and Thomas Trautmann1
1DLR — Deutsches Zentrum für Luft- und Raumfahrt, Institut für Methodik der
Fernerkundung, 82234 Oberpfaffenhofen, Germany
2DLR — Deutsches Zentrum für Luft- und Raumfahrt, Institut für
Planetenforschung, Rutherfordstr. 2, 12489 Berlin, Germany E-mail:
<EMAIL_ADDRESS>(FS)

(Accepted 2023 May 16. Received 2023 May 15; in original form 2023 February
22)

###### Abstract

Atmospheric temperatures are to be estimated from thermal emission spectra of
Earth-like exoplanets orbiting M-stars as observed by current and future
planned missions. To this end, a line-by-line radiative transfer code is used
to generate synthetic thermal infrared (TIR) observations. The range of
“observed” intensities provides a rough hint of the atmospheric temperature
range without any a priori knowledge. The equivalent brightness temperature
(related to intensities by Planck’s function) at certain wavenumbers can be
used to estimate the atmospheric temperature at corresponding altitudes. To
exploit the full information provided by the measurement we generalize
Chahine’s original approach and infer atmospheric temperatures from all
spectral data using the wavenumber-to-altitude mapping defined by the
weighting functions. Chahine relaxation allows an iterative refinement of this
“first guess”. Analysis of the $4.3\rm\,\mu m$ and $15\rm\,\mu m$ carbon
dioxide TIR bands enables an estimate of atmospheric temperatures for rocky
exoplanets even for low signal to noise ratios of 10 and medium resolution.
Inference of Trappist-1e temperatures is, however, more challenging especially
for $\mathrm{CO_{2}}$ dominated atmospheres: the “standard” $4.3\rm\,\mu m$
and $15\rm\,\mu m$ regions are optically thick and an extension of the
spectral range towards atmospheric window regions is important. If atmospheric
composition (essentially $\mathrm{CO_{2}}$ concentration) is known
temperatures can be estimated remarkably well; quality measures such as the
residual norm provide hints on incorrect abundances. In conclusion,
temperature in the mid atmosphere of Earth-like planets orbiting cooler stars
can be quickly estimated from thermal IR emission spectra with moderate
resolution.

###### keywords:

Astrobiology – Radiative transfer – Techniques: spectroscopic – Planets and
satellites: atmospheres – Infrared: planetary systems; Methods: data analysis

††pubyear: 2023††pagerange: Assessment of a Physics-based Retrieval of
Exoplanet Atmospheric Temperatures from Infrared Emission Spectra–Assessment
of a Physics-based Retrieval of Exoplanet Atmospheric Temperatures from
Infrared Emission Spectra

##  1 Introduction

A quarter century after the detection of the first planet orbiting a main
sequence star (Mayor & Queloz, 1995) exoplanetary science has developed
astonishingly quickly and is likely to continue on this path (Mérand et al.,
2021). Some 5000 exoplanets are known today111http://exoplanet.eu/, several
dedicated space missions are already in orbit for detection (Kepler/K2,
Transiting Exoplanet Survey Satellite (TESS), …), and characterisation
(CHaracterising ExOPlanet Satellite (CHEOPS, Benz et al., 2018)), and some
others are in development, e.g. PLAnetary Transits and Oscillations of stars
(PLATO, Rauer et al., 2014), Atmospheric Remote-sensing Infrared Exoplanet
Large-survey (ARIEL, Tinetti et al., 2018) or the proposed Large
Interferometer For Exoplanets (LIFE, Defrère et al., 2018; Quanz et al.,
2022a, b; Konrad et al., 2022). The successfully launched James Webb Space
Telescope (JWST) with its several infrared (IR) instruments will considerably
expand our knowledge for a wide range of exoplanets (Greene et al., 2016).

Atmospheric characterisation by means of microwave, IR, or ultraviolet
spectroscopy is performed routinely for Earth as well as for the Solar System
planets and moons (Hanel et al., 2003) and its feasibility has also been
demonstrated for extrasolar planets, mostly hot Jupiters (e.g. Madhusudhan &
Seager, 2009). For atmospheric remote sensing “optimal estimation” (Rodgers,
1976, 2000) is by far the most common inversion technique and this Bayesian
method has also been used in an exoplanet context, e.g. Irwin et al. (2008);
Lee et al. (2012); Barstow et al. (2013, 2016) or more recently Shulyak et al.
(2019). Bayesian methods such as optimal estimation (OE) heavily rely on a
priori knowledge that is in general readily available for Earth’s atmosphere
as well as many Solar System bodies, but which is less well-known for
exoplanets. Grid based search methods have become a standard approach for
exoplanet characterisation; because of the need to perform millions to
billions of radiative transfer forward model runs, sophisticated methods based
on Monte Carlo Markow Chains and Nested Sampling are exploited to speed up the
search (for recent reviews see e.g. Madhusudhan, 2018; Barstow & Heng, 2020;
MacDonald & Batalha, 2023).

Line et al. (2013) have performed an intercomparison of three retrieval codes
utilising OE and two Monte Carlo methods and concluded that for good
measurements (high spectral resolution and little noise) the estimates agree
quite well. More recently, Barstow et al. (2020) compared three retrieval
codes CHIMERA (Line et al., 2012), NEMESIS (Irwin et al., 2008), and Tau-REx
(Waldmann et al., 2015) and reported mostly consistent results but emphasised
the important role of radiative transfer modeling since the same inverse
problem solver (nested sampling MultiNest, Feroz et al., 2009) is used by all
three codes; hence the testing of additional inverse solvers is desirable.
Likewise, the critical impact of the forward model was also emphasised by
Barstow et al. (2022) who presented an intercomparison of five codes (Min et
al., 2020, ARCiS,, NEMESIS, Cubillos & Blecic, 2021, Pyrat BAY,, Tau-REx, and
MacDonald & Madhusudhan, 2017, POSEIDON,).

Pressure, temperature, and molecular concentrations depend on space (and
time), and a discretisation is mandatory for numerical analysis. The
latitudinal and longitudinal dependencies are largely ignored (for a
discussion of 2D or 3D effects on emission and transmission see e.g. Feng et
al. 2016; Blecic et al. 2017; Taylor et al. 2020 and Caldas et al. 2019;
MacDonald et al. 2020, respectively; see also Pluriel 2023 for a recent
review). Whereas assuming molecular concentrations which are constant in
altitude is likely acceptable as a first step, the assumption of isothermal
temperature profiles (e.g. Barstow et al., 2022) is problematic. Layer-by-
layer representations (standard for Earth remote sensing) have been criticised
as troublesome due to the limited information content (Line et al., 2013;
Parmentier & Guillot, 2014), and parameterised representations of the vertical
dependence of temperature are quite common (e.g. Madhusudhan & Seager, 2009;
von Paris et al., 2013; Morley et al., 2017a). A novel function expansion
approach along with a standard nonlinear least squares solver has been studied
by Schreier et al. (2020). Barstow & Heng (2020) recommended to “conduct
retrievals … with a variety of temperature structures and … to investigate
alternative approaches.”

Here we examine the feasibility of temperature sounding by means of an
iterative relaxation developed in the late sixties by Chahine (1968, 1970,
1972) that is presented in several textbooks on atmospheric radiation (e.g.
Goody & Yung, 1989; Liou, 1980; Hanel et al., 2003; Zdunkowski et al., 2007)
but rarely used today. According to “Subsection 6.5.2 — A physical approach to
retrieval” in Goody & Yung (1989) it is a “simple idea easy to visualize and
to extend to new circumstances” and exploits the properties of the weighting
functions (the derivatives of the transmission w.r.t. altitude). Extensions
and/or refinements of this approach were presented by Smith (1970) and Twomey
et al. (1977). These relaxation methods have been used for analysis of IR and
microwave observations of Earth’s atmosphere, and for temperature sounding of
Venus (e.g. Taylor et al., 1980), Mars (Lellouch et al., 1991a; Lellouch et
al., 1991b; Haus & Titov, 2000), and Jupiter (Gautier et al., 1977, 1979).

The organisation of the paper is as follows: The following section describes
our methodology (radiative transfer, Chahine relaxation, …) along with the
code and data. We demonstrate the feasibility to estimate the temperature for
Earth-like exoplanets in Section 3: Computationally fast estimates exploiting
a “mapping” of the wavenumber space to the altitude space as well as iterative
refinements based on a comparison of observed and model spectra are presented.
We continue with a discussion in Section 4 and give our conclusions in Section
5.

Figure 1: Transmission and weighting functions in the longwave TIR
($\mathrm{CO_{2}}$ $\nu_{2}$ band at $15\rm\,\mu m$, nadir view, Earth’s
midlatitude summer (MLS) atmosphere (Anderson et al., 1986), Gaussian response
function of half width $\Gamma=0.25\rm\,cm^{-1}$ corresponding to a resolution
$R\approx 2800$.). The upper panels show the transmission as a function of
wavenumber for several atmospheric paths (right) and transmission vs. altitude
(related to path length $s=z_{\text{ToA}}-z$) for selected wavenumbers (left).
The lower panel shows individual weighting functions for selected wavenumbers
(left) and a contour plot. Numbers in the left legends indicate the wavenumber
[cm-1] and the corresponding peak altitude [km] (bottom-left).

##  2 Theory

###  2.1 Forward model — infrared radiative transfer

In a gaseous atmosphere with local thermodynamic equilibrium the upwelling
intensity (radiance) $I$ at wavenumber $\nu$ is described by the Schwarzschild
equation of radiative transfer Goody & Yung (1989); Hanel et al. (2003)

$\displaystyle I(\nu)~{}$
$\displaystyle=~{}{\mathcal{T}}(\nu,0)\;B(\nu,T_{\text{surf}})~{}+~{}\int_{0}^{\tau(\nu,z)}B(\nu,T(\tau^{\prime}))~{}\exp{\bigl{(}-\tau^{\prime}(\nu)\bigr{)}}\,\mathrm{d}\tau^{\prime}$
(1)
$\displaystyle=~{}{\mathcal{T}}(\nu,0)\;B(\nu,T_{\text{surf}})~{}-~{}\int_{0}^{\infty}B(\nu,T(z^{\prime}))~{}{\partial{\mathcal{T}}(\nu,z^{\prime})\over\partial
z^{\prime}}\,\mathrm{d}z^{\prime}$ (2)

where $B=2hc^{2}\nu^{3}\big{/}\left[\exp{(hc\nu/k_{\text{B}}T)}-1\right]$ is
Planck’s function at temperature $T$ ($h$, $c$, $k_{\text{B}}$ are the Planck
constant, speed of light, and Boltzmann constant, respectively). For the
surface contribution
$I_{\text{surf}}(\nu)~{}\equiv~{}{\mathcal{T}}(\nu,0)\;B(\nu,T_{\text{surf}})$
we assume a Planck black body emission attenuated by the intermediate
atmosphere and a temperature identical to the bottom of the atmosphere (BoA)
temperature, i.e. $T_{\text{surf}}=T_{\text{BoA}}=T{(z=0)}$.

The monochromatic transmission ${\mathcal{T}}$, closely related to the optical
depth $\tau$, between observer and altitude $z$ is given by Beer’s law

$\displaystyle{\mathcal{T}}(\nu,z)~{}$
$\displaystyle=~{}\exp\bigl{(}-\tau(\nu,z)\bigr{)}$ (3)
$\displaystyle=~{}\exp\left[-\int_{z}^{\infty}\sum_{m}k_{m}\bigl{(}\nu,p(z^{\prime}),T(z^{\prime})\bigr{)}\>n_{m}(z^{\prime})\;dz^{\prime}\right]~{},$

with $n_{m}$ the density of molecule $m$, and $k_{m}$ the pressure and
temperature dependent absorption cross section obtained by summing over the
contributions from many lines. For simplicity a vertical path is assumed; for
a slant path with angle $\theta$ in a plane-parallel atmosphere replace
$z^{\prime}\longrightarrow z^{\prime}/\cos(\theta)$. The finite spectral
resolution of the instrument is taken into account by convolution of the
monochromatic intensity (1) (or transmission (3)) with a spectral response
function (SRF, e.g. Gaussian).

The upper atmosphere has a low abundance of absorbers and is therefore almost
transparent (i.e., transmission close to one); with an increasingly longer
atmospheric path (decreasing $z$) attenuation becomes stronger especially at
wavenumbers with strong absorption in the band or line center (Fig. 1 top-
right). Here path length $s$ refers to the distance to the observer,
essentially at “infinity”, in practice at top-of-atmosphere (ToA), and is
linked to altitude via $s=z_{\text{ToA}}-z$, cf. Eq. (3). Viewed as a function
of altitude the transmission decays rapidly to zero for these wavenumbers
(Fig. 1 top-left), i.e. photons from the lower atmosphere cannot penetrate to
space and the ToA radiation arises mainly from the upper atmospheric layers.
The so-called weighting function,222A note on terminology: different terms are
used for the derivative $\partial{\mathcal{T}}/\partial z$: weighting
functions Liou (1980); Hanel et al. (2003); Zdunkowski et al. (2007), kernel
functions Goody & Yung (1989), or contribution functions. Moreover, “Jacobian”
and “weighting function” are often used interchangeably; however, for
temperature retrievals using nonlinear least squares the Jacobian is the
partial derivative of the radiance w.r.t. the state vector $\vec{x}$,
$\partial I/\partial x$, where $\vec{x}$ is a discrete representation of the
temperature profile. the partial derivative
$K(\nu,z)\equiv\partial{\mathcal{T}}(\nu,z)/\partial z$ in (2), quantifies the
dominant contribution (or weight) of an altitude layer to the outgoing
radiation (Fig. 1 bottom panels, see also Liou (1980, Fig. 7.6), Goody & Yung
(1989, Fig. 6.17), and Hanel et al. (2003, Fig. 8.2.1)). Note that the
weighting function is similar, but not identical to the temperature Jacobian
$\partial I(\nu)/\partial T(z)$ that measures the radiation’s sensitivity to
changes of temperature (see Schreier et al., 2020, Fig. 10). Both the
weighting functions and the Jacobian clearly demonstrate that the radiation
carries little information of the lowermost and upper atmospheric layers.

Figure 2: “Observed” equivalent brightness temperatures $T_{\text{B}}$ (5) of
hypothetical Earth-like planets orbiting the M-dwarfs indicated in the legend.
Earth’s MLS atmosphere is shown for comparison (see subsection 2.4). Gaussian
response function of HWHM $\Gamma=0.25\rm\,cm^{-1}$ (longwave, left) and
$\Gamma=1.0\rm\,cm^{-1}$ (shortwave, right). For clarity no noise has been
added. The TIR-LW and TIR-SW spectra comprise 1441 and 281 pixels,
respectively. The numbers in the legend list the minimum and maximum
atmospheric temperature [K]. Wavelengths are given at the top (numerically
$\lambda[\rm\,\mu m]=10^{4}/\nu[\rm\,cm^{-1}]$).

###  2.2 Inversion — Chahine relaxation

According to Rodgers (1976) “the intensity to be measured is … a weighted mean
of the Planck function profile with the weighting function”. The bell-shape of
the weighting function (Fig. 1 bottom-left) can be exploited for analysis of
TIR spectra. Assuming a delta-function-like weighting function the
Schwarzschild equation (1) reduces (Hanel et al., 2003) to

$I(\tilde{\nu})~{}\approx~{}I_{\text{surf}}(\tilde{\nu})~{}+~{}B\bigl{(}\tilde{\nu},T(\tilde{z}_{\nu})\bigr{)}$
(4)

with $\tilde{z}_{\nu}$ the altitude where $K(\nu,z)$ has a maximum for a given
$\nu$. A first approximation of the atmospheric temperature can thus be
inferred from the observed Equivalent Brightness Temperature (EBT), i.e.

$\displaystyle T(\tilde{z}_{\nu})~{}\approx~{}T_{\text{B}}(\tilde{\nu}){}$
$\displaystyle\equiv~{}B^{-1}\Bigl{(}I_{\text{obs}}(\tilde{\nu})\,-\,I_{\text{surf}}(\tilde{\nu})\Bigr{)}$
(5)
$\displaystyle=~{}{hc\tilde{\nu}/k_{\text{B}}\over\log\Bigl{(}{2hc^{2}\tilde{\nu}^{3}/\bigl{(}I_{\text{obs}}(\tilde{\nu})\,-\,I_{\text{surf}}(\tilde{\nu})\bigr{)}}\Bigr{)}}$

This estimate can be iteratively improved using a relaxation scheme originally
proposed by Chahine (1968, 1970)

$T_{i+1}(\tilde{z}_{\nu})~{}\approx~{}B^{-1}\left({I_{\text{obs}}(\tilde{\nu})\,-\,I_{\text{surf}}(\tilde{\nu})\over
I_{\text{mod}}(\tilde{\nu},T_{i})\,-\,I_{\text{surf}}(\tilde{\nu})}~{}B\bigl{(}\tilde{\nu},T_{i}(\tilde{z}_{\nu})\bigr{)}\right)$
(6)

where $T_{i}$ denotes the temperature for iteration $i$ and
$I_{\text{mod}}(\tilde{\nu},T_{i})$ the corresponding modelled radiance
according to (1).

Figure 3: ToA intensities for M-Earths. The gray-shaded area indicates the
Planck function $B(\nu,T)$ for the minumum and maximum atmospheric
temperatures $T_{\text{min}}\approx 171\rm\,K$ for GJ 832 and
$T_{\text{max}}\approx 286\rm\,K$ by “construction”. (Computed with GARLIC
using atmospheric data from 1D-TERRA (Wunderlich et al., 2020); Gaussian
response function with resolution $R=1000$.)

###  2.3 Implementation

For our forward model we use Python for Computational ATmospheric Spectroscopy
(Py4CAtS, Schreier et al., 2019, available at
https://atmos.eoc.dlr.de/tools/py4cats/), a Python re-implementation of the
Generic Atmospheric Radiation Line-by-line Infrared Code (Schreier et al.,
2014). GARLIC has been thoroughly verified by intercomparison with other codes
(e.g. Schreier et al., 2018b) and validated by comparison to effective height
spectra (Schreier et al., 2018a) generated from Earth observations of the
Atmospheric Chemistry Experiment — Fourier transform spectrometer (ACE-FTS,
Bernath, 2017).

Py4CAtS and GARLIC compute molecular absorption cross sections $k_{m}$
assuming a (default) Voigt line shape (Schreier, 2018), where the wavenumber
grid point spacing is adjusted automatically for each molecule, pressure and
temperature to a fraction of the typical line width. Next, cross sections
scaled by molecular number densities are summed up to absorption coefficients;
then standard quadrature schemes are used to compute optical depths and
radiances. Both observed and modeled spectra are convolved with a Gaussian
spectral response function of constant width (with a default sampling of 4
points per half width at half maximum (HWHM)). Py4CAtS makes heavy use of
Numpy (van der Walt et al., 2011; Harris et al., 2020) (and occasionally SciPy
(Virtanen et al., 2020) and MatPlotLib (Hunter, 2007)).

The synthetic measurement spectrum is generated by adding generic Gaussian
noise (generated by the Numeric Python numpy.random.randn function)
independent of wavenumber,

The equivalent brightness temperature spectrum is obtained by “inversion” (5)
of Planck’s function and depicted in Fig. 2 for noise-free simulated
observations (see next subsection).

Figure 4: Weighting function maxima in the longwave TIR ($\mathrm{CO_{2}}$
$\nu_{2}$ band at $15\rm\,\mu m$ with Gaussian response function of half width
$\Gamma=0.25\rm\,cm^{-1}$, left) and shortwave TIR ($\mathrm{CO_{2}}$
$\nu_{3}$ band at $4.3\rm\,\mu m$ with $\Gamma=1.0\rm\,cm^{-1}$, right). The
length of these spectra (i.e. number of data points) is identical to those of
Fig. 2. The black dots indicate manually selected $\nu\leftrightarrow z$ pairs
(as also shown in Fig. 1 lower left).

###  2.4 Data

Atmospheric data for Earth are taken from the “AFGL Atmospheric Constituent
Profiles” atmospheres (Anderson et al., 1986) providing pressure and
temperature vs. altitude along with concentration profiles for 28 gases
including water vapor, carbon dioxide (with 360 ppm volume mixing ratio
(VMR)), ozone, methane.

Atmospheric data for the assumed Earth-like planets around M-dwarfs
(henceforth called “M-Earths”) are taken from Wunderlich et al. (2019). These
temperatures and concentrations are inferred from a 1D photochemistry model
(Gebauer et al., 2018) coupled to a climate model and are defined on 64 levels
with a ToA pressure of about $0.08\rm\,mb$. The former work assumed
hypothetical planets with Earth’s properties orbiting different M-dwarf stars
placed at the location where the planet receives modern Earth’s instellation.
For all M-Earths the surface temperature is approximately 288 K by appropriate
selection of the orbital radius, and the $\mathrm{CO_{2}}$ VMR is about 355
ppm in the lower atmosphere.

Data for Trappist-1e from Wunderlich et al. (2020) are derived from the Berlin
1D steady-state, cloud-free, radiative-convective photochemical model
1D-TERRA. This dataset comprises dry&dead, wet&dead, and wet&live scenarios
for $\mathrm{CO_{2}}$ surface partial pressures of $10^{-3}$, $10^{-2}$,
$10^{-1}$, and $10^{0}\rm\,bar$ (corresponding to VMRs of approximately
$10^{-3}$, $10^{-2}$, $10^{-1}$, and $0.5\cdot 10^{0}$ (see Table 11 in
Wunderlich et al. (2020) and is indicated by the exponent in subsection 3.7).
For comparison, the Trappist-1 planet of the M-Earth dataset (Wunderlich et
al., 2019) has a VMR of 355 ppm in the lower atmosphere. Note that this planet
is a purely hypothetical “Earth” orbiting Trappist-1, whereas the Trappist-1e
scenarios are based on orbital and stellar data (see Section 3.2 and Table 9
in Wunderlich et al., 2020).

The Chahine approach delivers temperatures only at a small set of altitudes,
but data on a moderately dense grid from BoA to ToA are required for the
radiative transfer modeling and we will use function expansion for
inter/extrapolation (see subsection 3.5). For the generation of the singular
vectors to be used as base vectors for this expansion we use the set of 42
Earth atmospheric profiles collected by Garand et al. (2001) augmented with
the eleven M-Earth temperatures regridded to a uniform altitude grid with 2 km
steps (Schreier et al., 2020, subsection 3.2). Note that the first six Garand
atmospheres correspond to the AFGL data.

Molecular line parameters are taken from the Hitran database; instead of the
most recent data (Gordon et al., 2022) (clearly mandatory for analysis of real
observations) we use data from the initial 1986 release (Rothman et al., 1987)
to speed-up the computations. (See the further discussion in subsection 4.3.)
Only the main IR absorbers are considered, i.e. $\mathrm{CO_{2}}$ and the
interfering species $\mathrm{H_{2}O}$, $\mathrm{CH_{4}}$, and
$\mathrm{O_{3}}$.

See also Schreier et al. (2020) for more details on atmospheric and molecular
data and a discussion of some of our approximations and assumptions.

##  3 Results

###  3.1 First preliminary constraints

Inspection of the Schwarzschild equation (1) can be used for a first estimate
of the range of atmospheric temperatures. Replacing the height-dependent
temperature $T(z)$ in the Planck function by the minimum atmospheric value
$T_{\text{min}}$, the equation simplifies to
$I(\nu)=B_{BoA}\mathrm{e}^{-\tau}+\int
B\mathrm{e}^{-\tau^{\prime}}\mathrm{d}\tau^{\prime}\geq
B_{min}\left(\mathrm{e}^{-\tau}+(1-\mathrm{e}^{-\tau^{\prime}})\right)=B_{\text{min}}$,
hence the ToA radiance is greater (or equal) to the minimum Planck function.
The upper limit can be derived in a similar manner, hence $B_{\text{min}}\leq
I(\nu)\leq B_{\text{max}}$.

The intensities shown in Fig. 3 confirm these constraints, i.e. the lower and
upper “Planck envelope” can be used as a preliminary estimate of the
temperature range. The minimum and maximum atmospheric temperatures can then
be readily estimated from the corresponding EBT minima and maxima.

The EBT spectra in Fig. 2 are, except for resolution, essentially a zoom-in on
the intensity spectra of Fig. 3 transformed via Eq. (5). In addition to the
strong absorption at $15\rm\,\mu m$ due to $\mathrm{CO_{2}}$ the ozone
fundamental band at $1042\rm\,cm^{-1}$ ($9.6\rm\,\mu m$) is clearly visible in
Fig. 3.

###  3.2 Mapping Wavenumbers to Altitudes

Before estimating atmospheric temperatures from IR spectra according to the
recipe of Eq. (5) several issues have to be addressed. First a set of
appropriate wavenumber-altitude pairs has to be identified. Suitable
wavenumbers are conveniently searched for in absorption band(s) of a molecule
with well-known concentration and ideally little variability, e.g. the
shortwave or longwave TIR bands of $\mathrm{CO_{2}}$. The corresponding
altitudes are given by the location of the weighting function maxima, hence
depend on the properties of the transmission ${\mathcal{T}}$ and as a
consequence depend on atmospheric temperature, pressure, and composition.
Obviously these data are unknown for exoplanets, but fortunately weighting
function of Earth-like exoplanets (with $\mathrm{N_{2}}$-$\mathrm{O_{2}}$
dominated atmospheres) are rather similar and closely resemble weighting
functions of typical Earth climates, see Fig. 4. Accordingly we will use the
$\nu\leftrightarrow z$ mapping of Earth’s atmosphere, e.g. for midlatitude and
subarctic summer or winter (MLS, MLW, SAS, SAW) in the $\nu_{2}$ longwave (LW,
wavelength $\lambda\approx 15\rm\,\mu m$) and $\nu_{3}$ shortwave (SW,
$\lambda\approx 4.3\rm\,\mu m$) bands of $\mathrm{CO_{2}}$. Note that these
spectra also depend on resolution (compare Figures 1 and 2 of Schreier et al.
(2020)), hence a lower resolution leads to smoother spectra and reduces the
sensitivity to upper atmospheric layers (See the discussion in subsection
4.1).

Figure 5: Estimate of AD Leo atmospheric temperature using a
$\nu\leftrightarrow z$ mapping with 8 pairs for MLS (see Fig. 4). The main
plot shows the ideal noise-free “observed” EBT (blue), the noise contaminated
EBT spectrum ($S/N=10$, cyan), and selected EBT values taken from the noise-
free spectrum; the inset compares the estimated temperatures with the true
profile. (TIR-LW $660\text{--}750\rm\,cm^{-1}$; Gaussian with
$\Gamma=0.25\rm\,cm^{-1}$.)

###  3.3 First Guesses — Selected Data in TIR-LW

Having identified the translation from wavenumber to altitude space
(henceforth called “mapping”) it appears to be straightforward to infer the
temperatures from the observed spectrum using Eq. (5). However, both real
spectra as well as our theoretical spectra simulating planned instrumental
measurements (see Fig. 2) are contaminated by noise, and exploiting single
data pairs $(\tilde{\nu},T_{\text{B}})$ is likely to lead to a “noisy”
temperature profile. This is confirmed by Fig. 5 where the temperature
retrieval for a hypothetical Earth-like planet orbiting AD Leo (Wunderlich et
al., 2019) is illustrated: the inset shows the eight temperature values taken
from the observed noise-free observation (intensity spectrum $I(\nu)$
converted to EBT $T_{\text{B}}(\nu)$ according to (5)). However, estimating
these temperatures from the noisy EBT spectrum would clearly lead to a zigzag
temperature profile.

Figure 6: First estimate of M-Earths’ atmospheric temperatures using a
$\nu\leftrightarrow z$ mapping with 8 pairs for MLS (see Fig. 4). The red and
magenta curves labeled “avg=20” are identical except for the different noise
vector. (Spectra as in Fig. 5)

To compensate for the noise the average of some neighboring pixels from the
observed EBT spectrum can be used instead. Figure 6 depicts temperature
estimates for the eleven M-Earth atmospheres (taken from Wunderlich et al.
(2019) and already used in Schreier et al. (2020)) and clearly demonstrates
that averaging EBTs from larger windows (e.g. 10 or 20 pixels) leads to the
inferred temperature profile becoming smoother. Averaging 20 pixels leads to a
significantly reduced zigzag, however, profiles estimated from two different
measurements (model spectrum contaminated by two randomly generated noise
vectors) are still distinct. For the M-Earths these first guess temperatures
are encouraging, for some cases almost “perfect” (e.g. AD Leo, GJ 644, and GJ
832 even without averaging), but in a few other cases (e.g. GJ 551, GJ 876,
and Trappist-1) deviations to the true profile are clearly visible. Moreover,
temperatures of Trappist-1 in the lowermost altitudes are significantly
underestimated.

Figure 7: Estimate of AD Leo atmospheric temperatures using all
$\nu\leftrightarrow z$ pairs for MLS. The colors and symbols for the $S/N=10$
observation indicate the altitude range where radiation is mainly coming from
according to the weighting functions. (For example green and cyan triangles
for the lowest and highest altitude. Inset similar to Fig. 5, spectra as
above.)

The profiles of GJ 551, GJ 876, and Trappist-1 show oscillations of up to a
few Kelvin. The cold trap of Trappist-1 is not very strong, i.e. the minimum
temperature is only 20 K cooler than the maximum temperature (288 K at BoA).
The tropopause temperatures of GJ 551 and GJ 876 are slightly cooler (258 K,
and 251 K, respectively). Moreover, the local temperature maximum of
Trappist-1 in the mid atmosphere at about 38 km is only modest and not
reproduced by the EBT estimate. Overall however, the results are rather
encouraging.

###  3.4 First Guesses — Exploiting the entire TIR-LW

Despite the promising results there are some caveats. Exploiting just a few
data points implies that the majority of data remain unused. Moreover, the
hand-picked selection of $\nu\leftrightarrow z$ mapping pairs is somewhat
arbitrary, and the estimated atmospheric temperatures are therefore likely to
change with different mappings. Consider for example the mapping pairs
$(677\rm\,cm^{-1},~{}18\,km),~{}(698\rm\,cm^{-1},12\,km),~{}(711\rm\,cm^{-1},~{}7\,km)$
(cf. Fig. 4): One issue is whether the mapping wavenumber should be chosen to
lie in the valley or on the peak of the absorption line. Choosing a point
somewhere in-between is less sensitive to resolution and less likely sensitive
to noise. Furthermore, the number of pixels which are used for averaging to
account for the noise is limited, in order to avoid overlap of spectral
regions. This is especially the case near the band center.

Figure 4 indicates that (except for the bottom and top altitudes) several
wavenumbers are sensitive to a particular altitude. For example, the
atmospheric layer around 32 km influences the radiance at several wavenumbers
in the $660\,\text{--}\,670\rm\,cm^{-1}$ interval and near $720\rm\,cm^{-1}$
(for some planets only). In order to address this issue, we therefore use the
mean EBT of all pixels contributing to a particular altitude according to the
weighting function peak height spectrum (Fig. 4). Because these peak altitudes
rarely coincide exactly to a given grid point we accept all pixels within a
given _tolerance interval $\delta z$_, i.e. to estimate the atmospheric
temperature at an altitude $z$ we take the average of all EBTs at wavenumbers
with a weighting function peak height in the interval $[z-\delta z,z+\delta
z]$. This concept is illustrated in Fig. 7 for a 5 km tolerance: Altitudes
beyond 40 km are seen only in a narrow interval around $667\rm\,cm^{-1}$, and
the average of all EBT’s in this interval is interpreted as atmospheric
temperature in the $40\text{--}50\rm\,km$ altitude range. Likewise,
temperatures for altitudes below 10 km are estimated from EBT’s at wavenumbers
beyond $700\rm\,cm^{-1}$ (except for the peak at $720\rm\,cm^{-1}$).

Figure 8: Estimate of M-Earths’ atmospheric temperatures using all MLS
$\nu\leftrightarrow z$ pairs. For the 1.0 km and 5.0 km tolerance mappings 18
and 4 altitudes have been found with $5\leq z\leq 39\rm\,km$ and $5\leq z\leq
35\rm\,km$, respectively. (Spectra as above.)

Figure 8 compares results for all M-Earths and for various $\delta
z$-intervals. In accordance with Fig. 4 no mappings are found for the upper
atmosphere, and the highest altitude point estimated depends on the magnitude
of the tolerance $\delta z$. In the “middle atmosphere” the estimated
temperatures are roughly equivalent. In some cases (e.g. GJ 664 or GJ 832)
where temperatures were slightly underestimated using 8 $(\nu,z)$ pairs only,
the deviation is reduced or eliminated. For the generalised Chahine case the
retrieved temperatures are often smoother compared to Fig. 6. Moreover, the
extended retrievals are less sensitive to noise: temperatures estimated with
the 5 km tolerance from two observations with different noise vectors are
largely identical except for the highest values (at 35 km) of GJ 551 and
Trappist-1 (compare the two “avg=20” estimates in Fig. 6). The warm M-Earths
(GJ 551, GJ 876, and Trappist-1) remain problematic with clear oscillations
above about 30 km and considerably underestimate near BoA especially for the 1
and 2 km tolerances, but otherwise the retrieved and true temperature profiles
are in good agreement with deviations to $T_{\text{true}}$ less than ten
Kelvin.

The oscillations of the profiles estimated with the 1 or 2 km tolerances can
be interpreted as follows: A closer look to Fig. 1 shows that the bell-shaped
weighting functions have a finite width of several kilometers in the lower
atmosphere and almost 10 kilometers in mid to the upper atmosphere. Hence
retrieving temperatures with one kilometer resolution is questionable.

###  3.5 Iterative Refinements — TIR-LW

The temperature estimates presented in the previous subsection are not always
satisfactory, but they can be used as initial guesses for the Chahine
relaxation scheme (6) or other iterative solvers like nonlinear least squares
(e.g. Schreier et al., 2020). Clearly the main advantage is the computational
speed, i.e. the whole “retrieval” is simply the inversion (5) of the Planck
function and therefore can be performed within fractions of a second. However,
because the peak heights of the weighting functions do not cover the entire
altitude range, this initial guess temperature values cannot readily be used
as an input for radiative transfer modeling. Moreover, with the $4$ or
$5\rm\,km$ tolerances only few temperature values are estimated in the mid
atmosphere.

Extrapolation appears to be a tempting solution, but this is well known to be
problematic. In Schreier et al. (2020) (Fig. 5) we have shown that Earth-like
temperature profiles can be represented as a linear combination of some base
vectors resulting from a singular value decomposition (SVD) of a large matrix
comprising “representative” temperatures (comprising the 42 Garand atmospheres
and the 11 M-Earth atmospheres introduced in subsection 2.4). Hence we will
use here a linear least squares fit to determine the expansion coefficients
for the “Chahine initial guess” profile and then “extrapolate” this profile to
the entire altitude range. In addition this is used for interpolation to a
dense altitude grid in the lower and mid atmosphere appropriate for radiative
transfer modeling. (For brevity this will be called “extrapolation”
henceforth.)

Before starting the iterations according to (6) one more issue has to be
discussed: how to stop the process. For nonlinear least squares solvers such
as MINPACK (More, 1978) or NL2SOL (Dennis, Jr. et al., 1981a, b) two
convergence criteria are usually employed: firstly the change of the estimated
state vector (here temperature) is small and the change of the residual norm
(the norm of the model minus observed signal vector, here the radiance
spectrum) is small (where “small” should be related to $S/N$). Inspection of
the radiance residual norm is clearly a natural choice for least squares.
Exploiting the deviation of the fitted temperature to the true temperature or
the deviation of the model to the “true” spectrum is clearly impossible for
analysis of true observations and is hence not used as a convergence
criterions.

Figure 9: Chahine relaxation for AD Leo starting with SAW pressure,
temperature (green), and concentrations ($S/N=100$, TIR-LW). The cyan crosses
show the initial guess (5) (similar to Fig. 7) and the cyan dotted line its
extrapolation. Intermediate temperature profiles (6) are shown in yellow, the
numbers in the legend indicate the maximum temperature change
($T_{i+1}-T_{i}$) and maximum EBT deviation (observed - model). The red curve
shows the final temperature profile with the maximum EBT deviation in the
legend.

For the analysis of the TIR-LW spectra using iterative Chahine relaxation (6)
we start with, e.g., Earth’s SAW atmospheric data; the SAW temperature is
needed to compute the atmospheric transmission ${\mathcal{T}}$ and surface
emission $I_{\text{surf}}$ in Eq. (5). The following stopping criteria have
been used: a maximum change of the updated temperature less than 5 K or a
maximum change of the EBT of less than 5 K. For all M-Earths the relaxation
stops after two or three iterations (see Fig. 9 for an illustrative example).
The results (Fig. 10) are consistent with those of Fig. 8: for most planets
the temperature is retrieved quite well for altitudes below about $40\rm\,km$,
but GJ 551 and Trappist-1 (and to a lesser extent GJ 1214) appear to be
problematic. Nevertheless, for all planets the cold trap temperatures closely
resemble the true one, and even for Trappist-1 the absolute temperature
differences are small.

Figure 10: Comparison of atmospheric temperatures using Chahine relaxation for
all Earth data (specified in the legend). The red solid line indicates the
best fit according to the norm $\|\Delta I\|$ of the residual radiance
spectrum (in units $\rm erg/s/(cm^{2}\,sr\,cm^{-1})$, see legend). The green
temperature shows the fit using the correct atmospheric densities. The last
subplot shows the temperature differences (true - fit) for the best Earth
model. (All $\nu\leftrightarrow z$ pairs with 5 km tolerance, extrapolation
with 4 base vectors, TIR-LW spectra as above.)

Further runs with other Earth model data (MLS, MLW, SAS, and tropical)
essentially confirm these findings, see Fig. 10. For planets with very low
minimum temperatures (AD Leo, GJ 644, or GJ 832) all models lead to almost
identical temperatures with small differences only in the upper atmosphere.
Differences are clearly visible for GJ 1214, GJ 551, GJ 876, and Trappist-1,
i.e. planets with relatively high minimum temperatures.

Although the most favorable model is unknown for real observations we can
nevertheless use the minimum radiance residuum norm $\|\Delta
I\|=\|I_{\text{obs}}-I_{\text{mod}}\|$ as a hint for selection of the “best”
fit. In all but two cases this selected solution corresponds to the optimum
solution according to the minimum norm of the EBT difference spectrum $\Delta
T_{\text{B}}$: for GJ 551 and GJ 876 MLS give the smallest EBT difference,
whereas SAS and TRO yield the smallest radiance difference (however, both
norms are identical within four digits). Using the mean EBT difference as
criterion gives also different solutions for GJ 581.

Table 1: Comparison of the radiance residuum norms (in radiance units $\rm erg/s/(cm^{2}\,sr\,cm^{-1})$) for iterative Chahine relaxation with the correct atmospheric pressure and concentrations with fits using one of Earth’s model atmospheres. The “min” and “max” columns give the range of residuum norms for these fits. The second column lists the range of atmospheric temperatures, i.e. $\max(T)-\min(T)$ in Kelvin. |  | $S/N=10$ | $S/N=100$
---|---|---|---
| $\Delta T$ | corr. | min | max | corr. | min | max
AD Leo | 100.6 | 249.7 | 245.3 | 249.5 | 42.43 | 41.2 | 49.2
GJ1214 | 35.1 | 330.1 | 328.2 | 331.7 | 58.94 | 60.4 | 71.5
GJ176 | 62.0 | 283.3 | 282.9 | 285.5 | 40.73 | 43.8 | 51.6
GJ436 | 56.7 | 286.8 | 286.4 | 288.7 | 32.14 | 36.6 | 47.4
GJ551 | 27.3 | 348.4 | 350.1 | 351.3 | 71.32 | 72.8 | 83.1
GJ581 | 51.5 | 304.7 | 302.7 | 305.1 | 34.03 | 39.0 | 48.7
GJ644 | 106.9 | 249.7 | 245.9 | 250.1 | 47.77 | 45.0 | 54.7
GJ667c | 69.8 | 281.3 | 276.4 | 282.1 | 36.33 | 34.8 | 48.5
GJ832 | 114.2 | 255.8 | 248.4 | 258.7 | 60.31 | 47.8 | 67.6
GJ876 | 35.2 | 313.2 | 312.4 | 313.7 | 55.08 | 57.8 | 67.7
Trappist1 | 20.2 | 371.1 | 372.0 | 376.1 | 90.76 | 92.8 | 101.8

In addition to the retrievals using the five Earth models Fig. 10 also shows
the temperature fitted with the correct exoplanet pressure and composition
(clearly unknown in case of real observations); the initial guess temperature
required to compute the weighting functions and total atmospheric transmission
${\mathcal{T}}(\nu,0)$ is set to the mean EBT. Interestingly the norm of the
radiance difference with the correct exoplanet composition is smaller than the
norm with the best Earth composition only for some of the planets. However,
Table 1 indicates that the variability of these residual norms is always
small, the difference of the largest and smallest norm is less than a few
percent, and the residual of all fits is very large because of the significant
noise of the synthetic intensities (i.e. the residual is dominated by noise).

For synthetic measurements with less noise the residual norm of fits with
different atmospheric models shows larger variations up to 41% for $S/N=100$
(because of the smaller noise the stopping criterium has been tightened to 3 K
temperature change). Note that the variation of the norms is especially large
for planets with a large range of atmospheric temperatures. Again the correct
atmosphere does not always deliver the best fit, but the best Earth model and
the correct densities always give similar temperature profiles. Table 1 shows
that the correct densities do not yield a model spectrum closer to the
observed spectrum for the four M-Earths with the largest temperature gradients
$\Delta T$ (AD Leo, GJ 644, GJ 667c, and GJ 832).

At first glance this “failure” appears to be quite disturbing and we interpret
this as follows: The residual norm is clearly the essential number
characterising least squares, i.e. the quantity to be minimised. However, it
is not the decisive quantity for Chahine relaxation; obviously the ratio of
observed to model spectra in (6) should finally be close to one, but there is
no uniquely defined number for the progress of the relaxation and quality of
the solution.

These results may also be considered as a hint that Chahine relaxation relies
on a strong simplification, i.e. the intensity (1) at a particular wavenumber
can be approximated by the Planck emission at a particular altitude according
to (4). This assumption clearly ignores the shape of the weighting functions
(compare Fig. 1 bottom-left) and is apparently problematic when the
temperature differences are large. Furthermore, large temperature differences
make the extrapolation more difficult. Nevertheless, the last subplot of Fig.
10 demonstrates that in the mid atmosphere (about 15 to 35 km) the temperature
can be estimated within $\pm 5\rm\,K$.

Figure 11: Estimate of atmospheric temperatures using all $\nu\leftrightarrow
z$ pairs for MLS and TIR-SW (2350 – 2420$\rm\,cm^{-1}$, $S/N=10$, Gaussian
with $\Gamma=1.0\rm\,cm^{-1}$.)

###  3.6 TIR-SW

The shortwave TIR appears to be less favourable for temperature retrieval for
several reasons: The radiance values in the TIR-LW are higher compared to the
TIR-SW (Fig. 11 in Schreier et al. (2020)), the TIR-LW is more favourable
because of the higher star-planet contrast, and the TIR-LW is also less
affected by scattering. Moreover, the TIR-SW weighting function peak heights
do not cover altitudes above 30 km (Fig. 4), and for shorter wavelengths
reflection of thermal radiation at the surface is becoming increasingly
important. On the other hand, the TIR-SW weighting functions indicate some
more sensitivity to the lowest atmosphere, and in fact IR instruments of
meteorological satellites used for sounding of Earth’s temperature (Menzel et
al., 2018), e.g. AIRS (Chahine et al., 2006) and IASI (Hilton et al., 2012),
exploit both regions.

Figure 12: Comparison of atmospheric temperatures estimated from the TIR-SW,
the TIR-LW, and the concatenated spectrum. (All $\nu\leftrightarrow z$ pairs
with a 4 km tolerance for MLS, $S/N=10$)

Figure 11 depicts the temperature inferred from the equivalent brightness
temperatures (5) in the TIR-SW using all $\nu\leftrightarrow z$ pairs for MLS.
Considering altitudes within one or two kilometers around a weighting function
peak height gives zigzag temperature profiles as in the longwave analysis, cf.
Fig. 8. For the $5\rm\,km$ tolerance the profiles are smoother, however only
four temperature values for the upper troposphere and lower stratosphere are
estimated. Compared to the TIR-LW estimates, Fig. 8, the profiles appear to be
somewhat smoother and closer to the true temperature. Trappist-1, however, is
reasonable only in the $10\text{\,--\,}20\rm\,km$ range. For other Earth
mappings (MLW, SAS, SAW, tropical) temperatures are almost identical (not
shown).

In Fig. 12 results from temperature estimates using the SW and LW interval
individually are compared with those from the combined spectrum (with 1722
data points, cf. Fig. 2). The ”data fusion” product clearly benefits from the
sensitivity of the longwave spectrum beyond $30\rm\,km$; however, it also
inherits the underestimated temperature in the lowest atmospheric levels,
especially for GJ 1214, GJ 551, GJ 876, and Trappist-1. A more sophisticated
data fusion approach might possibly be able to avoid the shortcomings of the
TIR-SW in the upper atmosphere and TIR-LW in the lower atmosphere.

###  3.7 Trappist-1e

The Trappist-1 planetary system (Gillon et al., 2017) orbiting a nearby
M-dwarf has attracted considerable attention because several terrestrial-type
exoplanets lie in the circumstellar habitable zone (e.g. Barstow & Irwin,
2016; Grimm et al., 2018; Krissansen-Totton et al., 2018; Lustig-Yaeger et
al., 2019; Fauchez et al., 2020; Krissansen-Totton & Fortney, 2022). Planets e
and f have been the object of numerous studies (e.g. Barstow & Irwin, 2016;
Morley et al., 2017b; Mikal-Evans, 2021). Recently, Wunderlich et al. (2020)
have used the newly developed radiation–convection–photochemistry model
1D-TERRA to study the feasibility of atmospheric characterisation, in
particular the possibility of finding any evidence for an ocean or biosphere.

Here we use these dry&dead, wet&dead, and wet&live scenarios with varying
$\mathrm{CO_{2}}$ levels (see subsection 2.4) for further tests of the Chahine
methodology. We generate synthetic observations in the same way as for the
M-Earths (Fig. 2), i.e. monochromatic intensity spectra according to (1) are
convolved with a Gaussian and contaminated with noise.

Figure 13: Trappist-1e temperatures from TIR-LW (660 – 750$\rm\,cm^{-1}$ with
$\Gamma=0.25\rm\,cm^{-1}$, cyan diamonds) and TIR-SW (2350 –
2420$\rm\,cm^{-1}$ with $\Gamma=1.0\rm\,cm^{-1}$, green squares) with entire
$\nu\leftrightarrow z$ MLS map. Third estimate with an extended TIR-LW (LWX,
660 – 850$\rm\,cm^{-1}$, red circles) and the correct map (i.e. correct
$\mathrm{CO_{2}}$). (All spectra with $S/N=10$). The number $n$ immediately
following the planetary scenario description in each subplot title indicates
the $\mathrm{CO_{2}}$ surface partial pressure, i.e. $p_{\text{CO2}}\approx
10^{-n}\rm\,bar$. (Regarding notation see subsection 2.4.) Minimum and maximum
atmospheric temperatures are listed in the title, the range of the EBTs of the
noise-free LW, SW, and LWX spectra are given inside the plot. The gray shaded
area shows the EBT range for the SW and LW combination.

Figure 13 is similar to Fig. 12 and shows atmospheric temperatures retrieved
from the TIR-LW and TIR-SW equivalent brightness temperatures (5) using the
MLS weighting function peak altitudes (cf. Fig. 4) with a 4 km tolerance
(without iteration). The LW and SW estimates are quite similar and relatively
smooth. For the $\mathrm{CO_{2}}$ $\text{VMR}=10^{-3}$ atmospheres (bottom
row) the temperature near BoA is somewhat underestimated, but in the middle
atmosphere up to about 40 km the estimates are close to the truth. However,
for the atmospheres with larger $\mathrm{CO_{2}}$ concentrations the
temperatures are clearly too cool.

Due to the large mixing ratios the atmospheres are optically thick in the
spectral windows considered here, and radiation (photons) from the lower warm
altitudes cannot propagate upwards to ToA (and the observer). This
interpretation is confirmed by the range of equivalent brightness temperatures
given in the plot, i.e. the warm low atmosphere does not show up in the EBT
spectra. Analysis of effective height spectra (Fig. 14) that might be
available from primary transit spectroscopy also indicates that the TIR-LW
interval (660 – 750$\rm\,cm^{-1}$) is insensitive to the lower atmosphere.

Figure 14: Trappist-1e effective heights for a Gaussian response function with
$R=1000$. (Computed with GARLIC.)

Furthermore, the success of the temperature estimates for the
$\text{VMR}=10^{-3}$ cases suggests that the $\nu\leftrightarrow z$ mappings
of Earth’s atmospheres are not adequate for Trappist-1e atmospheres with more
carbon dioxide. Fig. 15 shows that with higher $\mathrm{CO_{2}}$
concentrations the location of the weighting function maxima move upwards by 5
or even $10\rm\,km$, whereas the “nature” of the planet (wet vs. dry, dead vs.
alive) does not have a big impact.

Figure 15: Trappist-1e altitudes of weighting function maxima in the extended
LW-TIR (Gauss $\Gamma=0.25\rm\,cm^{-1}$). The wet&dead weighting function
maxima are largely similar to the wet&live maxima.

In summary, Figures 14 and 15 suggest that atmospheres with abundant
$\mathrm{CO_{2}}$ are not well reproduced from TIR-LW and TIR-SW spectra
especially in the lower regions (Fig. 13). Therefore we also considered an
extended spectral range along with the correct weighting function peak heights
(i.e. weighting functions in the extended TIR-LW computed with the true
concentrations) which considerably improved the estimate (esp. in the mid
atmosphere, Fig. 13 red dots; however, the sensitivity to the lower atmosphere
is lost, i.e. no mappings with altitudes below 10 or $20\rm\,km$ for the 10%
and 50% $\mathrm{CO_{2}}$ planets).

Having demonstrated that atmospheric temperatures can be estimated from
equivalent brightness temperatures (5) using an extended TIR-LW interval if
the composition is known, we now examine iterative Chahine relaxation (6). The
assumption of known abundances appears reasonable since some objects will have
information from transmission spectra. In particular quantifying
$\mathrm{CO_{2}}$ might be possible due to e.g. its dominant absorption bands
(clearly visible around 660$\rm\,cm^{-1}$ and 2400$\rm\,cm^{-1}$ in Fig. 14);
concentration estimates for other gases however might be more difficult.

First we examine whether we can then distinguish the temperature profiles for
the three scenarios (wet vs. dry, dead vs. live), i.e. we compare retrievals
using different surface and lower atmosphere scenarios (for example, for the
wet&live planet “observation” (right column in Fig. 16) we also investigate
fits using boundary parameters of the dead scenarios with the correct
$\mathrm{CO_{2}}$).

Figure 16 (top) shows that the effect of changing the scenario is largest for
the high $\mathrm{CO_{2}}$ abundance: in the lower atmosphere temperature is
significantly underestimated for the two dead scenarios which is likely
related to the fact that the effective heights (Fig. 14) and weighting
functions (Fig. 15) do not show any sensitivity to the lowest atmosphere (esp.
for the wet cases). For the observation of the wet&live $\text{VMR}\\!=\\!0.5$
planet (top right) both wet models yield an only moderately underestimated
temperature. The difficulty of the high $\mathrm{CO_{2}}$ atmospheres is also
demonstrated by the increased number of iterations required for convergence.
Fits with large residuum norm show the largest deviations between true and
retrieved temperature in the mid atmosphere. For moderate and low
$\mathrm{CO_{2}}$ (other rows in Fig. 16) the impact of the scenarios is only
weak, and only one or two iterations are necessary. Results here are close to
the true profile (blue line) in the low and mid atmosphere, although
deviations are clearly evident in the upper atmosphere where the weighting
function exhibits little sensitivity.

Figure 16: Chahine iterative estimate of Trappist-1e temperatures with correct
$\mathrm{CO_{2}}$ abundance and varying scenarios: wet & live (red), wet &
dead (cyan long dashed), and dry & dead (green dashed). SVD “extrapolation”
with four base vectors. Extended TIR-LW spectra as above. The legend lists the
number of iterations and the radiance residual norm.

A second set of runs has been conducted assuming MLS molecular profiles scaled
to the correct Trappist-1e column density, see Fig. 17. Similar to Earth the
Trappist-1e atmospheres have $\mathrm{CO_{2}}$ VMR profiles almost constant in
altitude, and the $\mathrm{H_{2}O}$ VMRs are strongly decreasing for pressures
below $100\rm\,mb$ (see Fig. 7 in Wunderlich et al., 2020). (In contrast the
$\mathrm{CO_{2}}$ VMR of all M-Earths increases by 10 – 20% for $p<100\rm\,mb$
(see Fig. 4 in Wunderlich et al., 2019).) However, at least for Earth, the
total optical depth is essentially dominated by the $\mathrm{CO_{2}}$
contribution (see Fig. A3 in Schreier et al., 2020), which suggests that for
these Trappist-1e atmospheres with even higher $\mathrm{CO_{2}}$
concentrations this dominance will be even stronger. Hence, the assumption of
isoprofiles does not have a strong impact on the quality of the retrievals,
i.e. temperature in the mid atmosphere can be estimated with little deviation
from the truth.

Figure 17: Chahine iterative estimate of Trappist-1e temperatures starting
with MLS $\nu\leftrightarrow z$ mapping and molecular profiles with correct
column density. The retrieved temperatures correspond to SVD “extrapolation”
with two (blue dash-dotted), three (red dashed), four (cyan long dashed), and
five (green dotted) base vectors. (Numbers in legend list the radiance
residual norm, spectra as above.)

Figure 17 also shows the importance of the “extrapolation” scheme based on the
expansion using singular vectors. Although a temperature profile
representation using only two base vectors delivers temperatures close to the
true even for the lowest atmosphere (including the $\text{VMR}=0.5\cdot
10^{0}$ cases, top row), a profile expansion with three or four base vectors
appears to be more reasonable. Using five base vectors works well for
$\text{VMR}\leq 10^{-2}$, but clearly fails for very high $\mathrm{CO_{2}}$
concentrations. The radiance residual norms shown in the plot also indicate
that three or four base vectors can be used reliably except for the high
$\mathrm{CO_{2}}$ concentrations.

Further tests have also been conducted with reduced or increased
$\mathrm{CO_{2}}$ concentrations. For moderate changes with quarter, half,
double and quadruple isoprofiles the retrieved temperature profiles are still
close to the true temperature. (See the further discussion in subsection 4.7.)

##  4 Discussion

###  4.1 Spectral resolution

The retrievals reported so far have been conducted assuming moderate
resolution TIR-LW and TIR-SW spectra. The motivation for $R>2500$ has been
discussed in Schreier et al. (2020, Figure 1) where we showed that for
decreasing resolution the sensitivity to upper atmospheric layers decreases:
for a Gaussian response function with HWHM $\Gamma=0.25\rm\,cm^{-1}$
(corresponding to a resolution $R=2800$ at $700\rm\,cm^{-1}$ (wavelength
$14.3\rm\,\mu m$)) the maximum of the weighting function peak height spectrum
(cf. Fig. 4) reaches altitudes above $40\rm\,km$; a coarser resolution reduces
the peak height (for $\Gamma=1.0\rm\,cm^{-1}$ ($R=700$),
$\Gamma=2.0\rm\,cm^{-1}$ ($R=350$), and $\Gamma=7.0\rm\,cm^{-1}$ ($R=100$) the
maxima lie at 40 km, 24 km, and 18 km, respectively).

The TIR-LW and TIR-SW $\mathrm{CO_{2}}$ bands can be observed by the Medium
Resolution Spectrometer (MRS) of the JWST Mid Infrared Instrument (MIRI) with
a resolving power $R$ of about 2500 (Rieke et al., 2015) (the Low Resolution
Spectrometer (LRS) only sees the TIR-SW band with a resolution $R=100$).
However, Morley et al. (2017b) caution that “for temperate planets
spectroscopy with JWST-MIRI will likely be unrealistically expensive”. In the
JWST Guaranteed Time Observations (GTO) program several exoplanets have
already been observed with the Near Infrared Imager and Slitless Spectrograph
(NIRISS), NIR Spectrograph (NIRSPEC), and NIR Camera (NIRCam) instruments,
that can deliver spectra with wavelengths up to $5\rm\,\mu m$ at low and
medium resolution.

For an assessment of the impact of resolution on the retrieval, synthetic
observations have been generated with different resolutions for selected
exoplanets. Noise has been adjusted assuming a square root relationship
between resolution and $S/N$. The results depicted in Fig. 18 confirm the
expectations discussed above, i.e. with decreasing resolution information on
the upper atmosphere is diminishing. In particular, for $R=500$ Chahine
estimates according to (6) are only available for altitudes below $20\rm\,km$
(shown as cyan circles) and “extrapolation” is clearly problematic.

Figure 18: Impact of resolution on Chahine iterative temperatures estimates.
MLS initial guess atmosphere, SVD “extrapolation” with 4 base vectors. The red
diamonds and cyan circles show the final update according to the relaxation
equation (6) for $R=2800$ and $R=500$, respectively. TIR-LW (M-Earths) or TIR-
LWX (Trappist-1e). The solid, dashed etc. lines show the corresponding
extrapolated temperature profiles.

###  4.2 Surface emission

Apart from the type of inversion (least squares vs. Chahine relaxation) the
methodology used here is largely identical to that of our previous feasibility
study. However, there is one large difference worth discussing: in Schreier et
al. (2020) we have ignored surface emission. Obviously the first term in (1)
is mandatory for modeling IR spectra in atmospheric window regions where the
atmosphere is relatively transparent such as the
$800\text{\,--\,}1200\rm\,cm^{-1}$ interval for Earth (except for
$\mathrm{O_{3}}$ absorption around $9.6\rm\,\mu m$). For the terrestrial
exoplanets, transmission ${\mathcal{T}}\approx 0$ in the $\mathrm{CO_{2}}$
bands considered in Schreier et al. (2020) (LW:
$660\text{\,--\,}720\rm\,cm^{-1}$), hence the surface contribution has no
impact on mid atmospheric temperature estimates.

However, for wavenumbers beyond $720\rm\,cm^{-1}$ corrections cannot be
neglected, so an estimate of the surface temperature is required. If
wavenumbers with transmission close to one are observed, the corresponding
radiance can be used to estimate $T_{\text{surf}}$. Here we approximate
surface temperature by the largest equivalent brightness temperature observed,
assuming that the maximum atmospheric temperature corresponds to the BoA
temperature and is approximately equal to the surface temperature. (Although
warmer temperatures are possible at the stratopause, these do not contribute
strongly to the spectrum as indicated by the weighting functions (cf. Fig.
1)).

Due however to noise, the largest EBT may not always be the best estimate of
the surface temperature. Alternatively the maximum of a smoothed EBT spectrum
(e.g. by running averages) or the mean of the largest ten etc. values could be
used.

###  4.3 Molecular spectroscopy data

Similar to Schreier et al. (2020) we have used the very first edition of the
Hitran database (Rothman et al., 1987) rather than the latest version to speed
up the computations. This simplification appears to be justified for the
feasibility study presented here, but is clearly inadequate for analysis of
real observations. In our previous study we discussed the then latest 2016
version (Gordon et al., 2017), here we provide a brief update for the current
Hitran 2020 (Gordon et al., 2022).

For a realistic modeling of IR emission spectra in the
$660\text{\,--\,}750\rm\,cm^{-1}$ interval line data in an enlarged interval
have to be considered to properly account for line wing contributions and
convolution with the spectral response function. Hence lines in the range
$648.75\text{\,--\,}761.25\rm\,cm^{-1}$ are read: Hitran 2020 knows 449445
lines of 29 molecules, which reduces to 79394 lines of the five main IR
absorbers ($\mathrm{CO_{2}}$, $\mathrm{O_{3}}$, $\mathrm{N_{2}O}$,
$\mathrm{CH_{4}}$, $\mathrm{H_{2}O}$). In contrast, Hitran 86 returns 16003
lines of 4 molecules; furthermore removing weak lines results in 4255 lines
total including 2827 lines of carbon dioxide. Despite this drastic reduction
of active lines the final spectra do not change significantly: For the MLS
atmosphere the peak height of the weighting function is modified by a few
hundred meters only (except for $\Delta z\approx 0.8\rm\,km$ at the high
wavenumber end of the interval); Likewise, the difference of the equivalent
brightness temperature spectra is usually less than $1\rm\,K$ (with a maximum
$|\Delta T_{\text{B}}|\approx 1.8\rm\,K$ at $718\rm\,cm^{-1}$).

Air-broadened half-widths are listed in the Hitran database since its
beginnings half a century ago, whereas self-broadening was introduced later
with the 1986 release (Rothman et al., 1987). Collision (pressure) broadening
parameters for perturbers other than “air” (i.e. Earth’s $\mathrm{N_{2}}$,
$\mathrm{O_{2}}$) can be important for modeling spectra of other planets and
were included only recently in Hitran (Gordon et al., 2017). For the
terrestrial $\mathrm{N_{2}-O_{2}}$ dominated atmospheres considered here the
Hitran 86 broadening parameters are clearly appropriate (with the possible
exception of the 50% $\mathrm{CO_{2}}$ Trappist-1e scenarios). In this context
it might also be worth noting that the line strengths in Hitran are also tuned
to Earth, i.e. the strengths are scaled by the relative natural abundance of
the isotopes in Earth’s atmosphere.333https://hitran.org/docs/iso-meta/

In addition to line transitions continua can be also important, especially
water self and foreign continua (Shine et al., 2012) and collision induced
absorption (CIA, Richard et al., 2012; Karman et al., 2019). Comparison of EBT
spectra (MLS atmosphere, $R=100$) computed with and without $\mathrm{H_{2}O}$,
$\mathrm{CO_{2}}$, $\mathrm{O_{2}}$ and $\mathrm{N_{2}}$ corrections (Clough
et al., 1989) indicates negligible differences in the TIR-LW regions
considered here; however, in the TIR-SW differences up to almost 8 K show up
around $2400\rm\,cm^{-1}$ in the right wing (short wavelength) of the
$\mathrm{CO_{2}}$ band. (For the entire TIR
($500\,\text{--}\,2500\rm\,cm^{-1}$) maximum differences up to $\approx
10\rm\,K$ at $1600\rm\,cm^{-1}$ can be seen.)

###  4.4 ToA 60 vs 120 km

In Schreier et al. (2020) all runs have been performed with a ToA at 60 km,
mainly because the Garand et al. (2001) atmospheres are defined only for
pressures down to $0.08\rm\,mb$ (about 60 to 65 km), but also to speed up the
computations. Likewise, the M-Earth atmospheres of Wunderlich et al. (2019)
have ToA altitudes in the range 61 to 74 km. However, for IR radiative
transfer modeling altitudes up to about 100 km might be important (the AFGL
data (Anderson et al., 1986) are given up to 120 km). For an assessment of the
importance of the upper atmosphere we have compared radiances for 60 and 120
km ToA and the MLS atmosphere. In the TIR-LW region the EBT difference is
usually less than 1 K except for the strong radiance peak at 668$\rm\,cm^{-1}$
and in the right wings of the $\mathrm{CO_{2}}$ $\nu_{2}$ band.

###  4.5 Geometry

Observations of exoplanet thermal emission will deliver disk averaged spectra
only, whereas we have used a single line-of-sight assuming a strict nadir
view. A common approximation for disk averaged spectra is to model a slant
path with about $35^{\circ}$ from nadir. The equivalent brightness spectra for
the vertical and slant paths differ by less than 2 K in the center of the TIR-
LW band (sensitive to the mid atmosphere) with somewhat larger differences for
$\nu>700\rm\,cm^{-1}$. On the other hand, the peak altitudes of the weighting
functions are shifted downwards by some 10 km (see Fig. 1 lower left panel).
As a consequence, the magnitude of the retrieved temperature will not change
significantly, but the associated altitudes could change by several
kilometers. It is also important to note that the preliminary constraints
discussed in subsection 3.1 are independent of the viewing angle.

###  4.6 Auxiliary data, initial guess and a priori

For the solution of “real world” inverse problems a large variety of auxiliary
data are required, e.g. instrument parameters or observation geometry. In the
case of atmospheric IR spectroscopy these auxiliary data also comprise
molecular optical properties (e.g. line lists or k-distributions); for
temperature sounding these data also include pressure and molecular
concentrations. These auxiliary data are often denoted “a priori”.

However, for optimal estimation (Rodgers, 1976, 2000) this solely refers to
the knowledge of the unknown state vector (e.g. temperature) prior to any
measurement (all other auxiliary parameters are treated as “model
parameters”). OE provides an estimate close to the a priori state vector where
the relative weight of observation and a priori is determined by the
respective measurement and a priori covariance matrices. (Note that this
weighting is independent of the weighting function
$\partial{\mathcal{T}}/\partial z$ defined in subsection 2.1.) Obviously a
priori knowledge of exoplanetary atmospheric properties is scarce (see e.g.
Shulyak et al. 2019 for a thorough discussion). Barstow (2020) noted that “the
dependence of OE on an informative prior means that it is less appropriate for
exoplanets”, which motivated the upgrade of the NEMESIS code (Irwin et al.,
2008) with MultiNest (Feroz et al., 2009). Monte Carlo type retrievals are
less sensitive to a priori, but for robust retrievals thousands to millions of
time expensive forward model evaluations are required (Fortney et al., 2021).

Some prior idea of the atmospheric state is required for Chahine-type
inversions to compute the weighting functions and the total atmospheric
transmission attenuating the surface emission (cf. Eqs. (1) and (5), (6)).
However, our simulations have demonstrated that the inferred M-Earth
temperatures are largely independent of the assumed Earth atmospheric model
and/or initial guess temperature profile.

Parameterisations as proposed by Madhusudhan & Seager (2009); von Paris et al.
(2013); Morley et al. (2017a) or Fossati et al. (2020) are frequently used for
temperature retrievals (for a detailed discussion see Section 3 of Barstow &
Heng, 2020). Obviously, these parameterisations are motivated by physical
insight and can also be considered as a kind of a priori knowledge.

For the iterative solution of nonlinear inverse problems by least squares
(e.g. Schreier et al., 2020), OE (e.g. Irwin et al., 2008) (technically a
constrained least squares) or Chahine relaxation, an initial guess is
mandatory. For the Earth-like exoplanets considered here a climatological
Earth temperature profile or an isoprofile defined by the mean of the observed
EBT can be used to start the iteration. Alternatively, an atmospheric
characterisation resulting from coupled photochemistry-climate codes such as
1D-TERRA (Wunderlich et al., 2020) can be used.

Figure 19: Impact of $\mathrm{CO_{2}}$ concentration. MLS initial temperature,
SVD “extrapolation” with 4 base vectors, $S/N=20$, TIR-LW (M-Earths) or TIR-
LWX (Trappist-1e) with Gaussian $0.25\rm\,cm^{-1}$. The legend indicates the
factor used to scale the $\mathrm{CO_{2}}$ VMR isoprofile of the model
atmosphere and the radiance residual norm.

###  4.7 Atmospheric composition

Thermal emission IR spectra are sensitive to temperature but provide little
information on atmospheric (molecular) composition. However these data are
mandatory for radiative transfer modeling. In particular knowledge of the
$\mathrm{CO_{2}}$ concentration is important for temperature sounding
exploiting its strong TIR-LW and TIR-SW bands. Actually the concentrations of
all molecules absorbing in the spectral region to be analysed are required: in
the case of the two TIR bands, $\mathrm{H_{2}O}$, ozone ($\mathrm{O_{3}}$,
longwave only) and methane ($\mathrm{CH_{4}}$, shortwave) absorption are
relevant. Note that $\mathrm{H_{2}O}$ does not have a large impact on the TIR-
LW weighting functions.

Assuming $\mathrm{CO_{2}}$ concentrations which are too low or too high likely
has an impact on the quality of the retrievals. Further test runs have been
performed with a setup similar to Fig. 17, but with $\mathrm{CO_{2}}$
isoprofiles scaled by factors from 0.25 up to 10 in the model atmosphere.
Figure 19 demonstrates that in most cases the residual norm is larger for
incorrect $\mathrm{CO_{2}}$ mixing ratios; in some cases the iteration fails
(e.g. for the 0.1% $\mathrm{CO_{2}}$ Trappist-1e atmospheres, usually because
of negative temperatures), the number of iterations becomes larger, or the
temperature profile shows stronger oscillations (zigzag profiles are
unrealistic because the observed radiance is an integral (1) that is
insensitive to small-scale “perturbations” of temperature and/or molecular
densities).

An independent estimate of concentrations is therefore desirable. Carbon
dioxide has several strong bands throughout the IR: in addition to the TIR
bands there are further strong rotation-vibration bands in the shortwave and
near IR (SWIR, NIR) at $2.7\rm\,\mu m$, $2.0\rm\,\mu m$, and $1.6\rm\,\mu m$.
In fact the NIR bands as well as the TIR-SW band enabled the identification of
$\mathrm{CO_{2}}$ in the atmosphere of the Saturn-mass exoplanet WASP-39b with
various JWST instruments (Ahrer et al., 2023; Alderson et al., 2023;
Rustamkulov et al., 2023) and the latter two SWIR bands are used operationally
for monitoring of Earth’s $\mathrm{CO_{2}}$ budget by several satellite
missions such as OCO-2/3 or GOSAT (Crisp et al., 2004; Kuze et al., 2009).
Regarding exoplanets the analysis of effective height / transit depth spectra
provided by primary transit observations is therefore valuable. Moreover, the
ratio of observed signals in two appropriate filter bands is a sensitive
indicator of $\mathrm{CO_{2}}$ atmospheric concentration (Rieke et al., 2015).
For the joint analysis of primary and secondary transits see also Griffith
(2014).

##  5 Summary and Conclusions

Temperature profiles of Earth-like exoplanets orbiting M-dwarfs have been
retrieved from synthetic thermal IR emission spectra using Chahine-type
relaxation methods along with a line-by-line radiative transfer code. The
essential assumption is that for a particular wavenumber the outgoing
radiation arises from a corresponding, well-defined altitude: In the band
centre absorption is strong and only radiation from the upper atmosphere will
be seen remotely; in the band wings absorption is weak and even photons from
the lower atmosphere can traverse the entire atmosphere to the ToA and beyond.
The feasibility of this method has been demonstrated using synthetic noise-
contaminated observations of various Earth-like planets with
$\mathrm{N_{2}}$-$\mathrm{O_{2}}$ dominated atmospheres orbiting M-dwarfs and
for Trappist-1e planets of different surface conditions (wet/dry and
dead/alive) and different carbon dioxide concentrations up to about 50%. (Note
that the assumption of $\mathrm{N_{2}}$-$\mathrm{O_{2}}$ dominance is
questionable for the high $\mathrm{CO_{2}}$ Trappist-1e scenarios.)

The equivalent brightness temperature (EBT) spectrum corresponding to the
observed intensity can be used to deliver temperature estimates extremely
quickly: The minimum and maximum EBT values provide first constraints on the
range of atmospheric temperatures independent of any a priori knowledge. The
EBT in the carbon dioxide absorption bands can be readily “translated” (within
seconds) to mid atmospheric temperatures. Using a handful of manually selected
data points is problematic because of the noise, hence averaging of some
neighboring spectral pixels is required (cf. Fig. 5). However, exploiting the
entire intensity spectrum and the corresponding wavenumber-altitude mapping,
as defined by the weighting functions, is preferably and clearly advantageous
(Fig. 7). Furthermore, using both TIR-LW and TIR-SW data can be helpful to
overcome the limited altitude sensitivity range of one region alone (Fig. 12).
In any case, the quality of this guess is however related to an appropriate
knowledge of the $\mathrm{CO_{2}}$ concentrations, in particular in the case
of the $\mathrm{CO_{2}}$-rich Trappist-1e planets.

Iterative relaxation allows a refinement of the first guess, however, the
success of the improvement relies on the inter/extrapolation used to complete
the limited set of $T$ data. Note that an update of the temperature will
change cross sections, transmission, and weighting functions and hence also
the wavenumber-altitude mapping.

Compared to classical nonlinear least squares fitting a clear advantage of the
Chahine relaxation (or the related Smith and Twomey schemes) is the fact that
no Jacobian (derivatives of the intensity $I$ w.r.t. to the state vector
elements) are required. This is clearly beneficial when finite differences are
used to approximate the Jacobian as for (nonlinear) optimal estimation or in
the SVEEEETIES study (Schreier et al., 2020) using Py4CAtS (GARLIC exploits
algorithmic differentiation, where the overhead for temperature Jacobians is
only about a factor 2, see Schreier et al. (2015)). On the other hand, the
time required for nonlinear least squares fitting (proportional to the number
of iterations) can be reduced if a good initial guess is provided. Of course
the “Chahine first guess” can be used as initial guess for any other iterative
optimisation solver.

In conclusion, we have used an extension of the classical Chahine approach in
an exoplanet context for the first time to our knowledge. This approach can
deliver stable and reasonable temperature estimates for terrestrial-type
exoplanets quickly (first guess in seconds, iterative refinements in minutes),
even for challenging cases such as atmospheres with weak inversions and large
$\mathrm{CO_{2}}$ abundances. It is also attractive in view of the growing
awareness on the “carbon footprint of large scale computing” and green
computing (e.g. Jahnke et al., 2020). Hence it allows interesting new insight
and provides a valuable addition to existing methods.

## Acknowledgements

This foundations for this research were the DFG projects SCHR 1125/3-1,
RA-714/7-1, and RA 714/9-1. We acknowledge the support of the DFG priority
programme SPP 1992 “Exploring the Diversity of Extrasolar Planets (GO
2610/2-1)”. J.L.G. thanks ISSI Team 464 for useful discussion. We thank the
German Research Foundation (DFG) for financial support via the project The
Influence of Cosmic Rays on Exoplanetary Atmospheric Biosignatures (Project
number 282759267). Finally we would like to thank Joanna Barstow for the
constructive review.

## Data Availability

Atmospheric data have been taken from Wunderlich et al. (2019, 2020) and
molecular spectroscopy data have been obtained from Hitran
(https://hitran.org/); the data analysis software is an extension of the
Py4CATS package (Schreier et al., 2019) available at
https://atmos.eoc.dlr.de/tools/Py4CAtS/.

## References

  * Ahrer et al. (2023) Ahrer E.-M., et al., 2023, Nature, 614, 649
  * Alderson et al. (2023) Alderson L., et al., 2023, Nature, 614, 664
  * Anderson et al. (1986) Anderson G., Clough S., Kneizys F., Chetwynd J., Shettle E., 1986, Technical Report TR-86-0110, AFGL Atmospheric Constituent Profiles (0 – $\rm 120\,km$). AFGL, Hanscom AFB, MA, USA
  * Barstow (2020) Barstow J. K., 2020, MNRAS, 497, 4183
  * Barstow & Heng (2020) Barstow J., Heng K., 2020, Space Sci. Rev., 216, 82
  * Barstow & Irwin (2016) Barstow J., Irwin P., 2016, MNRAS Letters, 461, L92
  * Barstow et al. (2013) Barstow J., Aigrain S., Irwin P., Fletcher L., Lee J.-M., 2013, MNRAS, 434, 2616
  * Barstow et al. (2016) Barstow J. K., Aigrain S., Irwin P. G. J., Kendrew S., Fletcher L. N., 2016, MNRAS, 458, 2657
  * Barstow et al. (2020) Barstow J., Changeat Q., Garland R., Line M., Rocchetto M., Waldmann I., 2020, MNRAS, 493, 4884
  * Barstow et al. (2022) Barstow J., Changeat Q., Chubb K. L., Cubillos P. E., Edwards B., MacDonald R. J., Min M., Waldmann I., 2022, Exp. Astron., 53, 447
  * Benz et al. (2018) Benz W., Ehrenreich D., Isaak K., 2018, in Deeg H., Belmonte J., eds, , Handbook of Exoplanets. Springer, doi:10.1007/978-3-319-55333-7_84
  * Bernath (2017) Bernath P., 2017, JQSRT, 186, 3
  * Blecic et al. (2017) Blecic J., Dobbs-Dixon I., Greene T., 2017, APJ, 848, 127
  * Caldas et al. (2019) Caldas A., Leconte J., Selsis F., Waldmann I. P., Bordé P., Rocchetto M., Charnay B., 2019, A & A, 623, A161
  * Chahine (1968) Chahine M., 1968, J. Opt. Soc. Am., 58, 1634
  * Chahine (1970) Chahine M., 1970, J. Atmos. Sci., 27, 960
  * Chahine (1972) Chahine M., 1972, J. Atmos. Sci., 29, 741
  * Chahine et al. (2006) Chahine M., et al., 2006, Bull. Am. Met. Soc., 87, 911
  * Clough et al. (1989) Clough S., Kneizys F., Davies R., 1989, Atm. Res., 23, 229
  * Crisp et al. (2004) Crisp D., et al., 2004, Adv. Space Res., 34, 700
  * Cubillos & Blecic (2021) Cubillos P. E., Blecic J., 2021, MNRAS, 505, 2675
  * Defrère et al. (2018) Defrère D., et al., 2018, Exp. Astron., 46, 543
  * Dennis, Jr. et al. (1981a) Dennis, Jr. J., Gay D., Welsch R., 1981a, ACM Trans. Math. Soft., 7, 348
  * Dennis, Jr. et al. (1981b) Dennis, Jr. J., Gay D., Welsch R., 1981b, ACM Trans. Math. Soft., 7, 369
  * Fauchez et al. (2020) Fauchez T. J., et al., 2020, Geosci. Model Dev., 13, 707
  * Feng et al. (2016) Feng Y., Line M., Fortney J., Stevenson K., Bean J., Kreidberg L., Parmentier V., 2016, APJ, 829, 52
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Fortney et al. (2021) Fortney J. J., Barstow J. K., Madhusudhan N., 2021, Atmospheric Modelling and Retrieval (arXiv:2109.01976)
  * Fossati et al. (2020) Fossati L., et al., 2020, A & A, 643, A131
  * Garand et al. (2001) Garand L., et al., 2001, JGR, 106, 24017
  * Gautier et al. (1977) Gautier D., Lacombe A., Revah I., 1977, J. Atmos. Sci., 34, 1130
  * Gautier et al. (1979) Gautier D., Marten A., Baluteau J. P., Lacombe A., 1979, Icarus, 37, 214
  * Gebauer et al. (2018) Gebauer S., Grenfell J., Lehmann R., Rauer H., 2018, ACS Earth Space Chem., 2, 1112
  * Gillon et al. (2017) Gillon M., et al., 2017, Nature, 542, 456
  * Goody & Yung (1989) Goody R., Yung Y., 1989, Atmospheric Radiation — Theoretical Basis, second edn. Oxford University Press
  * Gordon et al. (2017) Gordon I., et al., 2017, JQSRT, 203, 3
  * Gordon et al. (2022) Gordon I., et al., 2022, JQSRT, 277, 107949
  * Greene et al. (2016) Greene T. P., Line M. R., Montero C., Fortney J. J., Lustig-Yaeger J., Luther K., 2016, APJ, 817, 17
  * Griffith (2014) Griffith C. A., 2014, Philos. Trans. Roy. Soc. London, Ser. A, 372, 20130086
  * Grimm et al. (2018) Grimm S. L., et al., 2018, A & A, 613, A68
  * Hanel et al. (2003) Hanel R., Conrath B., Jennigs D., Samuelson R., 2003, Exploration of the Solar System by Infrared Remote Sensing, 2 edn. Cambridge University Press
  * Harris et al. (2020) Harris C., et al., 2020, Nature, 585, 357
  * Haus & Titov (2000) Haus R., Titov D., 2000, PSS, 48, 473
  * Hilton et al. (2012) Hilton F., et al., 2012, Bull. Am. Met. Soc., 93, 347
  * Hunter (2007) Hunter J., 2007, Comp. Sci. & Eng., 9, 90
  * Irwin et al. (2008) Irwin P., et al., 2008, JQSRT, 109, 1136
  * Jahnke et al. (2020) Jahnke K., et al., 2020, Nat. Astron., 4, 812
  * Karman et al. (2019) Karman T., et al., 2019, Icarus, 328, 160
  * Konrad et al. (2022) Konrad B. S., et al., 2022, A & A, 664, A23
  * Krissansen-Totton & Fortney (2022) Krissansen-Totton J., Fortney J. J., 2022, APJ, 933, 115
  * Krissansen-Totton et al. (2018) Krissansen-Totton J., Garland R., Irwin P., Catling D., 2018, AJ, 156, 114
  * Kuze et al. (2009) Kuze A., Suto H., Nakajima M., Hamazaki T., 2009, Appl. Opt., 48, 6716
  * Lee et al. (2012) Lee J.-M., Fletcher L., Irwin P., 2012, MNRAS, 420, 170
  * Lellouch et al. (1991a) Lellouch E., Encrenaz T., Phillips T., Falgarone E., Billebaud F., 1991a, PSS, 39, 209
  * Lellouch et al. (1991b) Lellouch E., Paubert G., Encrenaz T., 1991b, PSS, 39, 219
  * Line et al. (2012) Line M., Zhang X., Vasisht G., Natraj V., Chen P., Yung Y., 2012, APJ, 749, 93
  * Line et al. (2013) Line M. R., et al., 2013, APJ, 775, 137
  * Liou (1980) Liou K.-N., 1980, An Introduction to Atmospheric Radiation. Academic Press, Orlando
  * Lustig-Yaeger et al. (2019) Lustig-Yaeger J., Meadows V., Lincowski A., 2019, AJ, 158, 27
  * MacDonald & Batalha (2023) MacDonald R. J., Batalha N. E., 2023, Res. Notes AAS, 7, 54
  * MacDonald & Madhusudhan (2017) MacDonald R., Madhusudhan N., 2017, MNRAS, 469, 1979
  * MacDonald et al. (2020) MacDonald R. J., Goyal J. M., Lewis N. K., 2020, APJ, 893, L43
  * Madhusudhan (2018) Madhusudhan N., 2018, in Deeg H., Belmonte J., eds, , Handbook of Exoplanets. Springer, doi:10.1007/978-3-319-30648-3_104-1
  * Madhusudhan & Seager (2009) Madhusudhan N., Seager S., 2009, APJ, 707, 24
  * Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
  * Menzel et al. (2018) Menzel W. P., Schmit T. J., Zhang P., Li J., 2018, Bull. Am. Met. Soc., 99, 583
  * Mérand et al. (2021) Mérand A., et al., 2021, The Messenger, 184, 8
  * Mikal-Evans (2021) Mikal-Evans T., 2021, MNRAS, 510, 980
  * Min et al. (2020) Min M., Ormel C. W., Chubb K., Helling C., Kawashima Y., 2020, A & A, 642, A28
  * More (1978) More J., 1978, in Watson G., ed., , Numerical Analysis. Springer, Lecture Notes in Mathematics 630
  * Morley et al. (2017a) Morley C., Knutson H., Line M., Fortney J., Thorngren D., Marley M. S., Teal D., Lupu R., 2017a, AJ, 153, 86
  * Morley et al. (2017b) Morley C., Kreidberg L., Rustamkulov Z., Robinson T., Fortney J., 2017b, APJ, 850, 121
  * Parmentier & Guillot (2014) Parmentier V., Guillot T., 2014, A & A, 562, A133
  * Pluriel (2023) Pluriel W., 2023, Remote Sens., 15, 635
  * Quanz et al. (2022a) Quanz S., et al., 2022a, Exp. Astron.
  * Quanz et al. (2022b) Quanz S. P., et al., 2022b, A & A, 664, A21
  * Rauer et al. (2014) Rauer H., Catala C., et al., 2014, Exp. Astron., 38, 249
  * Richard et al. (2012) Richard C., et al., 2012, JQSRT, 113, 1276
  * Rieke et al. (2015) Rieke G. H., et al., 2015, PASP, 127, 584
  * Rodgers (1976) Rodgers C., 1976, Rev. Geophys. Space Phys., 14, 609
  * Rodgers (2000) Rodgers C., 2000, Inverse Methods for Atmospheric Sounding: Theory and Practise. World Scientific, Singapore
  * Rothman et al. (1987) Rothman L., et al., 1987, Appl. Opt., 26, 4058
  * Rustamkulov et al. (2023) Rustamkulov Z., et al., 2023, Nature, 614, 659
  * Schreier (2018) Schreier F., 2018, MNRAS, 479, 3068
  * Schreier et al. (2014) Schreier F., Gimeno García S., Hedelt P., Hess M., Mendrok J., Vasquez M., Xu J., 2014, JQSRT, 137, 29
  * Schreier et al. (2015) Schreier F., Gimeno García S., Vasquez M., Xu J., 2015, JQSRT, 164, 147
  * Schreier et al. (2018a) Schreier F., Städt S., Hedelt P., Godolt M., 2018a, Molec. Astrophys., 11, 1
  * Schreier et al. (2018b) Schreier F., Milz M., Buehler S., von Clarmann T., 2018b, JQSRT, 211, 64
  * Schreier et al. (2019) Schreier F., Gimeno García S., Hochstaffl P., Städt S., 2019, Atmosphere, 10, 262
  * Schreier et al. (2020) Schreier F., Städt S., Wunderlich F., Godolt M., Grenfell J., 2020, A & A, 633, A156
  * Shine et al. (2012) Shine K., Ptashnik I., Rädel G., 2012, Surv. Geophys., 33, 535
  * Shulyak et al. (2019) Shulyak D., Rengel M., Reiners A., Seemann U., Yan F., 2019, A & A, 629, A109
  * Smith (1970) Smith W., 1970, Appl. Opt., 9, 1993
  * Taylor et al. (1980) Taylor F., et al., 1980, JGR, 85, 7963
  * Taylor et al. (2020) Taylor J., Parmentier V., Irwin P., Aigrain S., Lee G., Krissansen-Totton J., 2020, MNRAS, 493, 4342
  * Tinetti et al. (2018) Tinetti G., et al., 2018, Exp. Astron., 46, 135
  * Twomey et al. (1977) Twomey S., Herman B., Rabinoff R., 1977, J. Atmos. Sci., 34, 1085
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Waldmann et al. (2015) Waldmann I., Tinetti G., Rocchetto M., Barton E., Yurchenko S., Tennyson J., 2015, APJ, 802
  * Wunderlich et al. (2019) Wunderlich F., et al., 2019, A & A, 624, A49
  * Wunderlich et al. (2020) Wunderlich F., et al., 2020, APJ, 901, 126
  * Zdunkowski et al. (2007) Zdunkowski W., Trautmann T., Bott A., 2007, Radiation in the Atmosphere — A Course in Theoretical Meteorology. Cambridge University Press
  * van der Walt et al. (2011) van der Walt S., Colbert S., Varoquaux G., 2011, Comp. Sci. & Eng., 13, 22
  * von Paris et al. (2013) von Paris P., Hedelt P., Schreier F., Trautmann T., 2013, A & A, 551, A120