Uppsala University

(1)

Uppsala University

Master’s Thesis

Supervisor: Prof. Dr. Nikolai Piskunov, Dr. Andreas Korn

Planning Observations of Terrestrial Exoplanets

Around M Type Stars with CRIRES+

Author: Jonas Zbinden

(2)

Svensk Sammanfattning

(3)

Abstract

(4)

1 Introduction

Today, the atmospheric composition of the planets in our Solar System has been roughly determined (de Pater and Lissauer,2010). Besides, discoveries of exoplanets have revealed that the architecture of

our Solar System is not typical at all (Wright et al., 2011). Therefore, we may expect very different

planetary compositions, and atmospheric properties than can be found in our Solar System. On the other side, investigating their compositions will bring us closer to answer one of the most interesting questions about our own existence. Are we alone in the universe? However, we shall see, before we are able to answer this question we still have a long way to go. The next step will be to characterize

planetary atmospheres and learn about the habitability of these remote and different worlds. To

shed light onto the question of atmospheric composition of exoplanets, the CRyogenic high-resolution

InfraRed Echelle Spectrograph (CRIRES)(Kaeufl et al., 2004) was developed. An upgraded version

of CRIRES, now called CRIRES+, is mounted at Paranal observatory on the very large telescope, (VLT), unit telescope 3 (UT3). To plan observations with CRIRES+ we developed a software tool called (CRIRES-planning-tool). Although it was developed to plan transit observations, it may easily be extended to other kinds of observations with CRIRES+. With the upgrade CRIRES+ has now enhanced observation efficiency, increased wavelength coverage, an upgrade of the focal plane detector array to increase the number of accessible diffraction orders, new wavelength calibration methods to

reach wavelength precision of 5 m/s and a new spectropolarimetric unit (Follert et al.,2014). Previous

scientific highlights with CRIRES were measuring the length of an exoplanet day for the first time (Snellen et al.,2014), creating the first weather map for the nearest brown dwarf to Earth (Crossfield et al., 2014) or finding the first superstorm on an exoplanet, HD209458 b and measuring its mass (Snellen et al.,2010a) to name but a few.

To plan transit observations of exoplanets we used the python libraries astropy, (Price-Whelan et al.,

2018) and astroplan, (Morris et al.,2018), which provide useful methods to assess times and

coordi-nates for observability at observatories all over the world. To calculate the signal-to-noise ratio (SNR) for an observation we use the Exposure Time Calculator (ETC), developed by the European Southern

Observatory (ESO) in Garching, Munich, Germany together with Uppsala University (UU)1.

CRIRES-planning-tool is written in python and selects candidates from the NASA Exoplanet Archive (Akeson et al.,2013) to assess which candidates can be observed when, what signal-to-noise ratio can be reached for CRIRES+, and ranks the planets after criteria deduced from these properties. The output of the software tool is used for decision making, which targets should be observed and when, and provides an important basis to plan for best efficiency in time and data quality. This manuscript provides an overview of the structure of the software tool and describes and discusses some of the used methods. The software comes with an extensive documentation. The README file provides information about the different functionalities, accessibility, structure, installation, and further development possibilities

of the CRIRES-planning-tool. The tool is available on GitHub2.

The thesis contains two major parts. The first part (chapter 2) gives an overview of the present state

of exoplanetary astrophysical techniques, both in observation and theory. In section 2.1we introduce

the reader to the technique of transit photometry, highlighting its strengths and limitations, and derive

the basic equations for transmission spectroscopy. In section2.2we present two important methods for

transmission spectra retrieval, and two techniques to compare retrieved spectra with synthesized spectra 1_{https://etctestpub.eso.org/observing/etc/crires}

2

(6)

from planetary models. In section 2.3 we provide an extensive overview about the different physics treated by a planetary model, such as radiative transfer, thermodynamics, convection, chemistry, clouds, and hydrodynamics and present examples to the different kinds of planetary models developed to study exoplanetary atmospheres. At the end of this part we will give a short summary and conclusion of the reviewed topics.

In the second part (chapter4), we shall provide the reader with a brief introduction to the CRIRES-planning-tool.

In sections4.1, and4.2, we explain the functionalities and used methods to select candidates for

obser-vations and assess their observability from the Paranal observatory, in section 4.3we give an overview

about the ETC and the computation of signal-to-noise ratio, and in section 4.4explain each parameter

used for the ETC. In section4.5we compare our tool with two other transit observation planning tools,

highlighting differences between these tools and verifying our own findings. Finally in 4.6 we present

preliminary results of possible transit observations over the course of one year. A discussion of the

CRIRES-planning-tool can be found in chapter5, and a short conclusion of the findings is given in6.

2 Theory and literature review

One idea to determine the composition of an exoplanet’s atmosphere is to retrieve the transmission spectrum of the star light passing through the atmosphere and investigate the spectrum for known molecular transitions. Since planetary atmospheres consist of molecules, and in case of hot Jupiters, atoms, one aims at observing only molecular bands. Subtracting the stellar spectrum, and if the observatory is based on Earth subtracting the telluric spectrum, the exoplanet’s transmission spectra should reveal its atmospheric composition and allow the deduction of the prevailing conditions at the surface of the planet. However, one is faced with several challenges. We will present these challenges in the following chapter one by one, and solutions or methods to circumnavigate these issues and retrieve as much information as possible about the composition and structure of an exoplanet’s atmosphere. We will discuss these solutions and point out where more research or improvement of the mentioned techniques is necessary.

2.1 Exoplanet transit

Transits of exoplanets in front of different stars have been observed thousands of times by now and are increasing by the month. The technique is based on the photometric observation of a star and investigation of its lightcurve for a periodically occurring dip in the received brightness or flux. The number of exoplanets that can be observed this way is limited by the inclination of their orbit with

respect to the observer. Assuming a planet of mass Mp and radius Rp transiting in front of its host star

with M_? and R_? for stellar mass and radius the orbit can be projected onto the plane of the observer.

Any point of the orbit of the exoplanet has a projected distance rsky to the barycenter of the orbit.

The distance r of a planet from the barycenter is given by

r = a(1 − e

2₎

1 + e cos f, (1)

and can be projected onto the plane of the observer, such that the projected distance r_sky takes the

form

rsky= r

q

(7)

a is the semi-major axis of the planet’s orbit, e the eccentricity, f the true anomaly, i the orbital inclination with respect to the projected plane of the observer, and ω is the argument of periapse. The term under the square root represents the projection of the orbit onto the observer’s plane. The

projected distance expressed in Cartesian coordinates takes the form r_sky =√X2_{+ Y}2_{, where the X}

axis is pointing in the direction of the descending node of the planet’s orbit. Following the conventions

of Winn(2010), X can be inferred directly from orbital mechanics and the projection term in rsky:

X = −r · cos (ω + f ) (3)

Applying X = 0 we get the condition for transit mid-time f_tra = +π₂ − ω. Similarly the condition for

occultation of the planet by the star is given by focc= −π₂− ω. From a photometric measurement of a

transit and occultation one can constrain the parameter space for ω and the eccentricity e. If possible, the transit observations can be combined with the radial velocity method to measure e directly and in this way also infer ω. With the exception of a and the inclination i, the parameter space of an orbit is

covered. Winn(2010) also deduce the probability of an observable transiting planet from the penumbra

and antumbra cones of a transiting planet:

ptra= R?± Rp a 1 + e sin ω 1 − e2 . (4)

The penumbra is represented here by the + sign and includes grazing orbits (not full planetary disk

occulting the star) and − represents fully enclosed transits. The probability p_tra to observe a transit is

therefore determined by the orbital parameters a, e, and ω. Using population studies for e (Xie et al.,

2016) and integrating out ω, we can estimate how many exoplanets we should be able to observe. For

our purpose we do not want to have an estimate of how many exoplanets might be observable in the

future (Brakensiek and Ragozzine, 2016), however, we would like to calculate the expectancy to find

another Earth-twin around a Sun-like star, i. e. a G-dwarf. Plugging in the same values as for Earth

in (4), we can compute the expected number of systems we must observe to find an Earth-sized planet

at 1 AU around a G-dwarf. With the fraction η of systems with an Earth like planet we must observe

N > 200 · η−1 systems to find a transiting Earth-twin. Using estimations for the relative abundance

of G-dwarfs (Mamajek,2016), the distribution of semi-major axis a around G-dwarfs (Han et al.,2014

(accessed July 20, 2020) (see Figure1), and the distribution of eccentricities e (Xie et al.,2016) of all exoplanets, we estimate η:

η = fa× fG× fe ≈ 0.00054, (5)

where f_a= ₁₄₉₇27 stands for the fraction of exoplanets found at 1 AU restricted to G-dwarf host stars,

fG ≈ 0.06 the fraction of G-dwarf stars in our galactic neighbourhood and fe ≈ 0.5 the fraction of

planets with eccentricity e ≈ 0. Plugging in the numbers from above we get N ≈ 3700000. This means

on average out of 3700000 stars, 1 star is a G-dwarf with a transiting planet the size of Earth at 1 AU

and a low eccentricity e below 0.1. Filtering the exoplanets archive (Han et al.,2014 (accessed July 20,

2020) for transiting planets around G-dwarfs with 0.8 < a < 1.2 and e < 0.1, we have not found any

(8)

0.01 0.1 1 10 200 150 100 50 0

Semi-Major Axis [Astronomical Units (AU)]

Distr

ib

ution

exoplanets.org | 7/20/2020

Figure 1: Distribution of the semi-major axis of exoplanets around G-dwarfs,

pro-duced with (Han et al.,2014 (accessed July 20, 2020).

Measuring the lightcurve and observing the transit of an exoplanet we can infer two other important parameters of an eclipsing system, namely the period and the transit duration. Both are important to

predict transit observability. We denote the period with P and the transit duration with Ttot. This

should not be confused with the full transit duration T_{f ull}, which stands for the time the complete

planetary disk is occulting the star. The time difference from T_tot to the beginning of T_{f ull} and vice

versa is called the ingress and egress time during a transit (see Figure 2). Using these properties we

can predict future transits of an exoplanet in question:

tc(n) = tc(0) + P n, (6)

where t_c is the time of conjunction or for our purposes the transit mid-time. The uncertainty in t_c(n)

is mostly sensitive to the uncertainty in P. Using Gaussian error propagation we get

∆tc(n) =

p

∆tc(0)2+ n2∆P2 ≈ n∆P. (7)

Other errors that may influence the predictions and planning of transit observations are influences from other bodies in the observed system, which are hard to constrain. However, provide an excellent way of discovering the presence of other unobserved objects. If we also would like to observe occultations we can get another estimate for the eccentricity e and the observer’s celestial longitude ω from the time

difference between transit and occultation ∆T_c

(9)

Additionally to observe exoplanets from an observatory on Earth, we must include the barycentric correction, which accounts for the light travel time from the observed system to the barycenter of the Solar System compared to the travel time to the observer. In the following expression U T C stands for universal time coordinated and if not mentioned otherwise, we always refer to UTC. BJD stands for

the barycentric julian date and JD for the julian date. ~r is the position of the observer with respect to

the barycenter of the Solar System, ˆn the unit vector from the observer to the observed system and c

is the speed of light. Then the barycentric correction for objects outside of our Solar System can be

written as (Eastman et al.,2010)

BJ DU T C ' J DU T C +

~r · ˆn

c . (9)

The correction term is interchangeable and we can compute J D_{U T C}, to compute our predicted transit

times. This method was used to convert the transit observation times t_cfrom the barycentric frame to

the observers frame, the Paranal observatory.

If one measures the radial velocity of the host star and thereby knows the eccentricity e, the velocity

semi-amplitude K?and the period P , one can estimate the mass of the planet Mp(Murray and Correia,

2010;Winn,2010): Mp (Mp+ M?)2/3 = K? √ 1 − e2 sin i P 2πG 1/3 . (10)

G is Newton’s gravitational constant. Using the transit method together with the radial velocity method the sin i degeneracy is approximately broken since in case of a transit sin i ≈ 1. However, normally

we are faced with the situation that Mp M? and we can only measure Mp/M?2/3. To measure the

mass of the planet the mass of the star must be determined by other means, such as spectral type,

luminosity, and spectroscopically determined properties, such as effective temperature Tef f, surface

gravity g?, and metallicity [Fe/H]. The surface gravity gp at the surface of the planet can be derived by

using the expression for the radial velocity semi amplitude of the star K_? (Murray and Correia,2010)

K? = Mp Mp+ M? 2πa/P sin i √ 1 − e2 ≈ Mp M? 2πa/P sin i √ 1 − e2 (11)

and Kepler’s third law. We obtain the surface gravity in terms of Rp/a, P, K?, and i:

gp= G Mp R2 p = 2π P √ 1 − e2_K ? (Rp/a)2sin i . (12)

Measuring the photometric lightcurve during a transit we can obtain the depth of the lightcurve δtra,

which provides us with further information about the observed planet star system, such as the ratio

between the radius of the planet and the star k = R_p/R?:

δtra ≈ k2 1 −Ip(ttra) I? . (13)

In case of negligible nightside emission Ip = 0 compared to the disk averaged stellar intensity I?, δtra

can be set equal to k2, and the stellar radius R_? can be computed from stellar models or approximately

by R_?∝ M_?−1/3. Thereby one can also solve the degeneracy between R_p and a, since by using Kepler’s

third law again we can compute a using M?. Notice that without making assumptions about the stellar

properties of the host star, the presented treatment to obtain Rp is always degenerate with either R?,

(10)

Winn(2010) present another way to approximate the stellar radius and the orbital inclination from the shape of the transit lightcurve. However, in a real observation of a lightcurve we are faced with limb darkening, which means the blocked fraction of the stellar flux decreases at the edges of the stellar disk, resulting in an increase of the received flux. To obtain information about the shape of the lightcurve such that one may determine transit, ingress and egress time, one needs to correct the lightcurves for limb darkening.

Knutson et al. (2007) propose a way to fit the lightcurves in different bandpasses simoultanously for the shape of a limb darkening corrected lightcurve, using the non-linear limb darkening law presented by Mandel and Agol(2002):

ˆ I(r) = 1 − 4 X n=1 cn(1 − µn/2) (14) ˆ

I(r) is the specific intensity at position r on the stellar disk, where ˆI(0) = 1 and 0 ≤ r ≤ 1, and 0

stands for the center of the stellar disk, while 1 stands for the outer edge of the stellar disc. r is related

to µ through µ = cos θ =√1 − r2_{, where θ represents the angle between the line of sight and and the}

line of emergence of the flux, standing perpendicular on the stellar surface. c_n are the four coefficients

to fit the limb darkening from stellar models or from photometric transit lightcurves. A more frequently

encountered limb darkening law is the quadratic limb darkening with c1 = c3 = 0, c2 = γ1+ 2γ2 and

c4 = −γ2, and the new fitting coefficients γ1, 2. The choice of using the quadratic over the full treatment

is due to higher efficiency and less complexity in fitting lightcurves. Mandel and Agol (2002) present

the full treatment of possible transit scenarios and the following lightcurve geometry with application

of limb darkening. Knutson et al.(2007) apply the non-linear limb darkening for small planets k ≤ 0.1

which returns a set of equations to fit the limb darkening coefficients simultaneously with the decoupled parameters k and i. The calculations are lengthy and complicated and we forgo presenting them here

and only show the expressions derived byWinn (2010). Another technique developed to fit the specific

intensity dependent on limb darkening was proposed by Aronson and Piskunov (2018) and basically

follows a similar idea as described by Mandel and Agol(2002).

The impact parameter b is depicted in Figure 2, and is fundamentally related to the geometry of the

orbit of the transiting planet:

btra= a cos i R? 1 − e2 1 + e sin ω , (15)

where the + in the denominator must be replaced with − for occultations b_occ. If we can retrieve

the limb darkening corrected shape of the lightcurve, we can decouple R?/a from b using the time of

ingress and egress (Ttot− Tf ull) and computing the scaled stellar radius R?/a for non-grazing orbits

Rp R? a and small eccentricities in terms of Ttot and Tf ull:

b2 = (1 −

√

δtra)2− (Tf ull/Ttot)2(1 +

√ δtra)2 1 − (Tf ull/Ttot)2 (16) R? a = π 2δ1/4_tra q T_tot2 − T2 f ull P 1 + e sin ω √ 1 − e2 (17)

These expressions can be derived by computing Ttotand Tf ull, also depicted in Figure2, in terms of the

shape of the lightcurve and the orbit and inverting for b2and R?/a. Note that we have now independent

(11)

Figure 2: Illustration of the different quantities introduced to describe the

lightcurve of a transiting planet in front of its host star. T represents T_{f ull}, T_tot

would cover the times from t_I until t_IV, t_I−II is the ingress time, and t_III−IV the

egress time. b is called the impact parameter of the transit and is closely related

to the orbital parameters of the transiting planet (Winn,2010)

.

some input parameters for the star, either from fitted observational data, with stellar models, or from other theoretical assumptions for the star.

Using photometric lightcurves together with radial velocity observations enables us to constrain the

parameter space of the planet’s orbit, make good estimates for the radius of the planet R_p and the

surface gravity gp and can be used as inputs for planetary climate and atmospheric models. However,

the observed planetary radius Rpfrom a lightcurve is not well defined since cloudy or hazy atmospheres

can alter the fraction of blocked stellar flux to total flux and defining the surface of a gas planet might be

ambiguous. Nevertheless, the different values for Rpwe may get by observing the lightcurve at different

wavelengths is at the base of the idea of transmission spectroscopy of exoplanetary atmospheres, which is further elaborated in the next section.

2.1.1 Photometric lightcurve and transmission spectroscopy

Using spectroscopy to probe transmission spectra of a planetary atmosphere needs some assumptions

first. We define the radius of the planet Rp as the opaque part of the planet in all wavelengths, deduced

from a photometric lightcurve as derived above. This means that the planetary radius Rp may be set

(12)

The photometric lightcurve, observed at frequency ν is defined by the fraction of the received flux to the stellar flux:

f (t)ν ≡ Fν(t) Fν?(t) = 1 + Ip(t) I? −      k2_α tra(t), transits 0, outside eclipse k2_α occ(t)Iνp_I_ν?(t), occultations (18)

where we have used the definition for the monochromatic flux (Pradhan and Nahar,2011)

Fν = Z ∂Σ Iνcos θdΩ = 2π Z 1 −1 Iνµdµ (19)

with Iν = Ip(t)+I?the received intensity in [ergs/(cm2 s str Hz)]. Here, the polar angle θ stands for the

angle under which the light was emitted from the surface dA on the star. The monochromatic flux F_ν

is given in [ergs/(cm2 s Hz)] or [W/(m2s Hz)]. αtra and αoccdescribe the geometry of the area blocking

the star light during a transit or the area adding to emitted light from the planet during occultations

and out of transit. Mandel and Agol (2002) give the full analytic expressions for all possible cases

that might be encountered doing transit observations. The different cases are derived from geometrical considerations.

Measuring the radius of the planet Rp at a certain wavelength might return a different radius than

obtained from photometry. We rewrite the new radius as Rν = Rp+ δRν depending on the wavelength,

respectively the frequency ν. Then the part adding to the new occulting planetary disk can be described

by an annulus 2πRpδRν+ πδR2ν with thickness δRν. The thickness of the annulus δRν can be related

to the pressure scale height of the atmosphere observed at wavelength ν via the optical depth τν at ν.

The optical depth is the integrated opacity along a path C weighted with the density ρ:

τν =

Z

C

κν(s)ρ(s)ds. (20)

The opacity can be related to the atomic absorption coefficient αν: κν = ανn, where n is the relative

number density of a certain species. The expression for the absorption coefficient α is often seen to describe a single absorption line, whereas the opacity κ is used in general for radiative transfer with its continuous contributions from absorption and emission as well. Depending on the geometry and

the physics we are investigating, κν or αν can be dependent on the integration path C. We will now

show how the thickness of the annulus δR_ν can be related to the absorption at ν. For more details, see

section 2.3.3and (Gray,2005).

Assume we have an atmopshere consisting of a single species. For a grazing light ray traveling through

the annulus of the semi-transparent atmosphere at an impact parameter xν, we can compute the

tangential optical depth τ_ν according toWebb and Wormleaton(2001) andBrown(2001). The geometry

is depcited in Figure 3. In this derivation the geometry does not take into account any effects of

refraction.

The tangential optical depth in this geometry is given by

(13)

Figure 3: Geometry for a grazing light ray traveling through the semi-transparent atmosphere of a planet at impact parameter x, maximum atmospheric height H

and planetary radius R_p (Webb and Wormleaton,2001)

.

with the molecular mass µ of the species in question. In case of hydrostatic equilibrium we can use a barometric law for the number density of the species with height r:

n(r) = n0e−r/h, h =

kBT

µg , (22)

where the pressure scale height h is dependent on the surface gravity g, under the assumption of

H Rp, such that g does not change much with altitude, the temperature T of the atmosphere, and

kB is the Boltzmann factor. This may be generalized to non-isothermal atmospheres with different

temperature dependencies in different layers and non-homogeneous processes such as strong irradiation effects on one side of the planet in case of tidally locked planets. Assuming the absorption coefficient to be relatively constant along a tangential line of sight, the integral above can be evaluated using the

modified Bessel function K1((x + Rp)/h) (Abramowitz and Stegun,1965):

τν(x) = 2n0κν(x + Rp)eRp/hK1

(x + Rp)

h

(23) In the absence of absorption lines we can estimate the overall pressure scale height h to one order of

magnitude accuracy (Benneke and Seager, 2012; Miller-Ricci et al., 2008) by measuring the impact

parameter x at two different wavelengths under the assumption τ (x1,2) = τ (δR1,2) = 1. Equating both

expressions and solving for h, we get

h = x1− x2 lnα1 α2 q x1 x2 , (24)

where κ1 and κ2 are the opacity or scattering cross section at the two different wavelengths λ1 and λ2

(14)

To investigate spectral features, we need to relate the measured impact parameter x to the decrease

of the stellar flux by the semi-transparent atmosphere: δF?,ν. The specific photometric light curve fν

(18) in terms of the blocked stellar flux δF?,ν can be expressed as:

fν = 1 −

δF?,ν

F?,ν

, (25)

where the second term is given by (Brown,2001)

δF?,ν F?,ν = 1 πR2 ? Z ∞ 0 2π(Rp+ x0)[1 − exp (−τν(x0))]dx0, (26)

and one integrates over all impact parameters. τ is the integral along every line of sight at each impact parameter. The factor 2π arises from integrating around the whole annulus visible by the observer, neglecting effects of rotation of the planet, the star or wind in the atmosphere of the planet.

Following the treatment of Brown(2001) and the method to measure the transmission spectrum of an

exoplanetary atmosphere by Aronson and Waldén (2015), we may equate the expression for δF_ν with

ia· P , where the latter describes the spectrum of the semi-transparent atmosphere expressed in terms

of the measured impact parameter xν:

ip&a(·2πRpxν+ x2ν+ R2p) = Z ∞ 0 2π(Rp+ x0)Iν(µ(φ))[1 − exp (−τν(x0))]dx0 (27) ip&a(φ) = Z Ap&a(φ) Iν(µ(φ))dµ(φ) (28)

where ip&a denotes the specific intensity from the star covered by the planet and its atmosphere. For

I(µ(φ)) one can use a quadratic Limb darkening law or spectral models. The dependency of µ on the

orbital parameter φ can be looked up in Figure 4. It describes the relation between the position of

the planet on its orbit around the star to the observer’s view of the position of the planet. What we have not treated in our derivation of δF/F is the influence of the refractive index and the night side

emission of the planet. Brown (2001) point out that the path of a light ray through the atmosphere

is not perfectly straight. However to compute the complete obscurred area by the atmosphere of the planet they explain that the minimum height of a tangential light ray not altered by refraction is

relevant. They do not discuss night side emission, andAronson and Waldén(2015) do not discuss the

influence of refraction or night side emission and neglect it. Using the previous arguments we may now present their novel retrieval method.

2.2 Atmospheric retrieval methods for transmission spectroscopy

2.2.1 A novel planetary model-independent transmission spectrum recovery method

Aronson and Waldén (2015) describe a method to recover transmission spectra of exoplanetary atmo-spheres in the near-infrared from Earth based observatories with high-resolution spectrometry. The method should be applied on data from CRIRES+ in the future. The main idea is that even in presence of strong telluric absorption CRIRES+ will be able to measure in between the telluric lines and observe features shifting with time due to a time dependent Doppler shift, while the planet moves across the surface of the star. The full signal received by a high-resolution spectrograph can be expressed as received flux per exposure normalized to the exposure time:

(15)

with F (ν, t) the stellar flux and T (ν, t) the telluric transmission and f (v_p, vs) a Doppler correction

factor for the movement of the planet with respect to the star, the rotation of the star, the movement

of the Earth, and the rotation of the Earth. vs and vp are the velocities of the star and the planet

relative to Earth depending on their present position during the transit. They are dependent also on the orbital parameters φ at t and the shape of the orbit (a, ω, e, i). The geometry is depicted in

Figure 4. Additional Doppler shifts can be introduced by planetary winds across the terminator (the

only part of the atmosphere observable with transmission spectroscopy) and the rotation of the planet, however, a reasonable assumption is that the planets are tidally locked and so the influence of the planets rotation is small. On the other hand, horizontal winds arise from the strong stellar irradiation

onto a tidally locked planet and dissipate some of the energy to the night side. Snellen et al.(2010b)

used CRIRES, the predecessor of CRIRES+, to measure the high altitude winds of the gas giant HD 209458 b with transmission spectroscopy and a cross correlation technique that we will present in the

next section 2.2.2. They showed that these winds can reach speeds of several 100 ms−1, close to the

sound speed of the transported gases and can therefore alter the spectrum through line broadening and overall wavelength shifts significantly. To obtain the complete received signal one needs to convolve the spectrum with the instrument profile Γ. It can be approximated by a Gaussian with full-width at half maximum equal to the spectral resolution of the spectrograph and a noise function for each pixel. The

noise function follows a Poisson distribution with a width of λ = _S/N1 , where S/N is the signal-to-noise

ratio. The signal is then finally multiplied with a continuum normalization η. Applying these factors we get Sn,syn(ν, t) =h ˜Sn(ν, t) ⊗ Γ · (1 − N (ν, φ)) i · η. (30) 𝜙 𝜇! 𝜇" 𝐴#&%(𝜙)

Figure 4: In this figure we depict the geometry of an orbit in relation to the transit observation. The dependency of µ on φ is dependent on the shape of the orbit given by the orbital param-eters (a, e, i, ω) and the orientation of the orbit observed from Earth.

The continuum normalization is necessary due to the limited wavelength range of high-precision spectrographs and generally to correct for any in-strumental changes in the signal over time. Con-ventionally the continuum is normalized via fit-ting of the recovered spectrum or removing all broadband signal contributions and computing the relative line strengths to the total stellar

flux (de Kok et al.,2013). Aronson and Waldén

(2015) use sensitivity curves, derived from

mea-suring the Solar spectrum before and after the science exposure to compare with a theoretical Solar spectrum and normalize the measured spec-tra accordingly. Additionally they present a tech-nique to correct the measured spectra with a brightness correction parameter ν(t), to control the continuum contribution to the retrieved spec-trum. Computing this parameter gives not much insight into the physics at play and we forego its presentation. To compute the impact

parame-ter xν from Sn, which gives rise to the recovered

(16)

Aronson and Waldén(2015) formulate the recovery of the spectrum as a minimization problem, similar to a least square minimization between synthetic and measured signal, respectively spectra:

Φν =

X

n

(ν(t) · Sn(ν, t) − Sn,syn(ν, t))2 = min, (31)

where Sn,syn(ν, t) denotes the synthetic spectra according to (30). Notice, in Sn the parameter xν is

free, while all the other contributions must be constrained prior to retrieving x_ν. The minimization

problem is then rewritten in the form

Ω = Φ + Λ · R = min (32) R =X ν dP dν !2 (33) P = δRν · f (vs, vp) πR2 ? (34) where P denotes the part of the absorption by the atmosphere of the planet normalized to the projected stellar surface. The parameter Λ is used to balance the relative importance of the regularization term R, which balances the smoothness of the solution with the precision of the fit. By imposing the condition of all derivatives vanishing at the minimum and linearizing Ω in R, we get a system of equations for

xν that is only dependent on the synthetic input spectra F (ν, t, vs), Iν(µ(φ), vs) and T (ν, t) of the star

and the Earth. The flux spectrum for F (ν, t, v_s) is generated using the stellar atmosphere model code

MARCS (Gustafsson, B. et al., 2008) with the stellar parameters effective temperature Tef f, surface

gravity log g? and metallicity [F e/H] and is corrected with out of transit observations. The telluric

absorption spectrum T (ν, t) is generated from the synthetic atmosphere spectrum algorithm to model

transmission spectra and radiance, LinePak (Gordley et al., 1994). However, Aronson and Waldén

(2015) recommend for the future to use tools for synthesizing high-precision telluric transmission spectra

fitted to observation. I_ν(µ(φ), vs) can be generated the same way as F or by using limb darkening

equations, as mentioned earlier. However, Aronson and Waldén (2015) report higher accuracy using

stellar models, althoughAronson and Piskunov(2018) andCzesla et al.(2015) argue that using stellar

models to compute the limb darkened specific intensity fails at reproducing center to limb variations even

for the Sun. Aronson and Waldén(2015) propose to find the optimal value for Λ by comparing simulated

exoplanet spectra to the synthetic one, generating a simulated spectra with the known parameters of the exoplanet as the one observed. The sky emission received by the telescope can either be included

through models or can be corrected for by nodding (Chromey, 2010). To test their method, Aronson

and Waldén(2015) generated synthetic planetary atmospheres with LinePak and used the transmission

spectra to use their recovery method on. In this way, they could directly compare the recovered

spectrum with the synthetic spectrum. They tested three types of exoplanetary atmospheres: super-Earth with super-Earth-like atmosphere, super-super-Earth with Venus-like atmosphere, both orbiting M-dwarfs and a hot Jupiter orbiting a Sun-like star. The transmission spectrum of the Venus-like atmosphere could be recovered within one transit of 75 min, the Earth-like atmosphere could be recovered by combining 4 transits of each 90 min and the hot Jupiter case could be recovered in a single transit. However, recent analysis including instrumental stability and instrumental drift from tests of CRIRES+ revealed that for the Earth like case about 10-15 transits, for the Venus-like case about 5-10 transits, and for the hot Jupiter case about 5 transits will be necessary, (private communication N. Piskunov). To test their

method under different circumstances Aronson and Waldén (2015) introduced systematic errors, such

(17)

inserting a T_{ef f} 5% too high. They also comment on stellar activity through star spots, the influence of clouds, Rayleigh scattering, or no atmosphere. Constant over- or underestimation of the telluric transmittance leads to complete loss of the transmission spectrum. Starspots and stellar flux/intensity variability increases the uncertainty in the received spectrum. Young active M-dwarfs with significant number of starspots are not well understood and may make spectrum recovery challenging, if not impossible. Tests with Solar like stars and including spectra of starspots showed the possibility of complete recovery of the transmission spectrum. However, they point out that star spots or stellar activity altering the flux of the observed star significantly could not be properly simulated by stellar atmospheric models and might lead to loss of the transmission spectral features. No atmosphere, a very thin atmosphere, or a Rayleigh scattering dominated atmosphere could not be distinguished from each other. Measuring the Rayleigh scattering slope using ground based spectroscopy turns out to be difficult to achieve, due to Rayleigh scattering in Earth’s atmosphere in the visual band. Also errors in planetary size can lead to over or underestimation of the atmospheric height, however, this can be corrected for by setting the continuum to zero, using the technique to compute ν(t) mentioned before.

Finally they present the application on data from Snellen et al.(2010b) and showed the possibility to

retrieve the same CO lines asSnellen et al.(2010b) have. To reach the specified recovery performance, a

certain signal-to-noise ratio S/N is required. The requirements can be found in (Aronson and Waldén,

2015) and are applied in our planning tool to determine if the signal from a transit is strong enough

for atmospheric spectra retrieval. The minimum S/N requirement can be different for other methods,

such as the cross correlation method, described in section2.2.2.

What is important to notice with this method, is that it recovers the transmission spectrum prior to assuming any knowlegde about the atmospheric composition or structure of the exoplanet, however, at the cost of more observation time to obtain higher S/N and higher sensitivity to data quality, resolution and wavelength coverage than other methods.

2.2.2 Cross correlation methods to retrieve the transmission spectrum of an exoplanetary

atmosphere

The cross correlation between two continuous functions f and g is defined as: (f ? g)(τ ) ≡

Z ∞

−∞

f (t)g(t + τ )dt, (35)

where f (t) is the complex conjugate of f (t). In our case we can denote the signal fν(t) as the received

signal Sn, adjusted for all the instrument and continuum contributions and gν(t) as the synthetic signal

Sn,syn, generated from an atmospheric modeling code for transmission spectroscopy, such as LinePak

(Gordley et al., 1994). What the cross correlation then measures is the similarity between S_n and

Sn,syn, as a function of velocity shift from the transit midpoint: ∆ν_ν ∼ ∆v_v . This technique has been

widely used for previous studies in transmission spectroscopy (Brogi et al., 2013, 2012;de Kok et al.,

2013; Snellen et al., 2010b). Following the procedure used by de Kok et al. (2013), first the telluric and stellar contributions need to be removed. To align the spectra on a common wavelength grid, they corrected the spectra with a second order polynomial describing a wavelength shift for each pixel, obtained by the best cross correlated synthetic telluric and stellar spectra with the measured spectra.

Next, they propose to use single value decomposition (Kalman,1996). The single value decomposition

of any matrix A containing all the spectra of a single detector is given by:

(18)

The matrix U contains all the right singular vectors, which are the eigenvectors of A, V are the left singular vectors, which represent the eigenvectors scaled by the singular values of A. The matrix W is then a diagonal matrix containing the singular values of A. The main idea in applying singular value decomposition normally is that A can be approximated by the highest singular values of W. However, here the procedure is turned around and one is only interested in the lower singular values, while the higher singular values, containing the major contribution to the measured spectra, are set to zero. The major contribution comes from the telluric absorption lines and to some degree from the host star. The key to retrieve the planetary transmission spectra relies again on the fact that the shift of the absorption lines from the planetary atmosphere due to the change in radial velocity is greater than the shift in the telluric and stellar lines. Thus, the time varying information of A can be isolated to compare with model spectra by using cross correlation to obtain the best fit from a model grid. The

model spectra used by de Kok et al. (2013) were generated with the parametric forward model from

(Madhusudhan and Seager,2009) (see chapter2.3). Madhusudhan and Seager(2009) report constraints on the shape of the temperature-pressure profile (P-T) of HD 189733 b and the relative abundances of the major chemical constituents as well. However, they had not been able to successfully constrain the

presence or abundance of CO.de Kok et al.(2013) report the presence of CO in the atmosphere of HD

189733 b at a 5σ significance. Still, due to the degeneracy of pressure, temperature and abundance, the abundance of CO could not be constrained from the line cores. Although constraints on P-T

profiles exist from (Madhusudhan and Seager,2009) and others,de Kok et al.(2013) report that there

is still an ongoing discussion about the exact P-T profiles. Previous low-resolution observations of HD 189733 b have not provided strong enough constraints on P-T profiles at low pressure, which is where the observed spectral line cores are most sensitive. Nevertheless, they mention that with a wider

wavelength coverage, observing other species such as H₂O and CH₄, and by measuring the H₂ induced

Rayleigh scattering, it would be possible to infer more exact P-T profiles and determine the absolute

abundances of gas species. Benneke and Seager(2012) for instance present a retrieval method with their

parametric forward model to distinguish hydrogen H₂ dominated atmospheres from water vapour H₂O

dominated atmospheres and to constrain absolute abundance of the gas species assumed to dominate

exoplanetary atmospheres. Brogi et al. (2012) use a similar method as de Kok et al. (2013) to prove

the presence of CO in the atmosphere of τ Boötis b at a 6σ confidence level, and Brogi et al. (2013)

find absorbtion by CO and water in the atmosphere of 51 Pegasi b at a 5.9σ confidence level. All the above mentioned observations applying the cross correlation method were conducted at the VLT with

CRIRES (Kaeufl et al.,2004).

2.2.3 Optimal estimation and bayesian inference

Traditionally for atmospheric retrieval a spectrum of parametric forward models, such as (Madhusudhan

and Seager, 2009) would be generated with up to ∼ 107 models, (∼ 107 points in a ten-dimensional parameter space) and compared via least square fit or cross correlation, to find the model closest to the

recovered spectrum. Madhusudhan and Seager(2009)’s model includes 6 parameters for the

pressure-temperature profile P-T and 4 parameters for the chemical composition. Another model from Guillot

(2010) used by Benneke and Seager (2012) uses one parameter less for the P-T profile. The models

are both described in detail in chapter 2.3. On the other side, self-consistent planetary atmospheric

(19)

by parameters and does not necessarily require complete physical consistency. On the other hand self-consistent models give raise to their P-T profile naturally from the modelled physical processes. A standard method to compare recovered spectra to self-consistent models is the method of optimal Estimation, which is described next.

Optimal estimation (OE) is suitable to apply in high-resolution spectroscopy with a wide wavelength coverage and high signal-to-noise ratio. High-resolution spectroscopy allows to identify single features in the transmission spectrum and compare to detailed self-consistent model atmospheres. Optimal estima-tion optimizes the likelihood funcestima-tion with least square minimizaestima-tion using the Levenberg-Marquardt

algorithm. Shulyak, D. et al.(2019) developed their code for emission spectroscopy. Nevertheless, the

method holds also for transmission spectroscopy and they present the OE method to be used with data

from CRIRES+ and the τ -REx software package (Waldmann et al., 2015). (As it happens, τ -REx

is also equipped with Bayesian inference methods in form of Markov-Chain Monte Carlo MCMC and nested sampling), The OE parameter fitting is used to obtain the global minimum of the cost function

φ containing the χ2 sum and the deviation from the initial guess state vector x_i, which contains the

parameters of the model fitted to observation (Rodgers,2004). The deviation of a measured signal y

to the synthetic signal yi = F (xi) at iteration i generated from the parameters in xi, and the deviation

of x_i from an initial guess x_a are assumed to be Gaussian and follow a normal distribution with the

co-variance matrices Sy and Sa:

P (y|x) = 1 (2π)n/2 _{mod S}1/2 y exp −1 2(y − yi) T_S−1 y (y − yi) (37) P (x) = 1 (2π)n/2 _{mod S}1/2 a exp −1 2(xi− xa) T_S−1 a (xi− xa) . (38)

The OE method is then based on Bayes’ theorem and the above assumptions. Bayes’ theorem is also the key to Bayesian inference and therefore we quickly recall it here:

P (y|x) = P (x|y)P (y)

P (x) . (39)

We can rewrite this expression in the form

− 2 ln P (y|x) = −2 ln P (x|y) − 2 ln P (y) − (−2 ln P (x)) (40)

and using P(y|x), P(x) we can define the cost function φ

φ = (y − yi)TSy−1(y − yi) + (xi− xa)TSa−1(xi− xa). (41)

P (y) is fixed by the measured signal y and can therefore simply be treated as a normalization constant.

The co-variance matrix Sa is generally given by the expectation value E of the difference of xi from

the initial guess x_a:

S_amn = E((xm_i − x_am)(xn_i − xn_a)). (42)

For the application in planetary atmospheric forward models Shulyak, D. et al. (2019) assume the

off-diagonal parameters to be correlated by the correlation length lcorr which depicts the number of

layers between layer m and n and p_m, pn represent the pressure in those layers

S_amn=pSmm

a Sann exp (

−| ln (pm/pn)|

lcorr

(20)

Sy contains the measurement errors and the estimated errors for the parameters in x. The final errors

after convergence can be calculated from the requirement of lim

xi→ˆx

∇xiφ = 0, ˆx is the vector containing

the parameters at the minimum of the cost function: lim

xi→ˆx

−K_iTS_y−1(y − yi) + Sa−1(xi− xa) = 0. (44)

Ki = ∇xiF (xi) is the Jacobian of the forward model at iteration step xi and arises from acting with ∇

on (y − yi). Expressing yi as

yi ≈ F (xa) + ∇xiF (xi)(xi− xa), (45)

where xais assumed to be close to the minimum, (44) then becomes

lim xi→ˆx −K_iTS_y−1y + K_iTS_y−1F (xa) + KiTS −1 y Kixi− KiTS −1 y Kixa+ Sa−1xi− Sa−1xa= 0. (46)

Collecting all the terms quadratic in xi in (41), all other terms remain constant varying xi and we can

compare (46) with the alternative form of lim

xi→ˆx ∇xiφ = 0, ˆx: lim xi→ˆx ∇xiφ = lim xi→ˆx ∇xi(xi− ˆx) T_S_ˆ−1_(x i− ˆx) = lim xi→ˆx ( ˆS−1xi− ˆS−1x) = 0,ˆ (47)

and thereby identify the co-variance matrix of the final result ˆx at the global minimum of φ as ˆS,

containing the co-variance of the model parameters, the measurement errors and the pre-assumed errors of the forward model:

ˆ

S = (S_a−1+ K_iTS−1_y Ki)−1. (48)

The construction of the iteration step from xi to xi+1 can be deduced from setting (44)=(47) and

solving for ˆx, and is given by:

xi+1= xi+(1 + γi)Sa−1+ KiTSy−1Ki

−1

KT

i Sy−1(y − yi) − Sa−1(xi− xa) . (49)

γ is a fine-tuning parameter controlling the balance between the measurement and the initial guess.

Notice that in case of x_i+1= xi = ˆx and γ = 0 the first term in brackets becomes ˆS and the second term

becomes zero due to the minimization condition of the cost function φ. Although it might seam like a well-motivated assumption that the posterior probability distribution P (y|x) are following a normal

distribution, Benneke and Seager (2012) criticise this assumption. Using a parametric forward model

from (Guillot,2010), they present the method of Bayesian inference, described next.

The method of Bayesian inference itself is also based on Bayes’theorem (39) as is the optimal Estimation

method and combined with a suitable parameter iteration method investigates the full parameter distribution around the minimum of the likelihood function

L = L₀exp (−χ2_n/2), (50) χ2_n= N X j=1 (yj− F (xnj))2 σ2 j , (51)

where σj denotes the standard deviation of the j-th measured datapoint yj. Looking at the different

terms in (39), we have P (y|x) = L the likelihood function, the probability of the measured data y under

(21)

xn at iteration step n under the data set y. P (x) = Π is the prior parameter probability distribution and P (y) = Z is called the evidence, independent of n. Z is the likelihood of the data marginalized over the parameter space and is normally treated as a normalization factor for the employed model. It

can be used to compare different models with each other (Madhusudhan,2018). The evidence Z can

be rewritten as the integral over the likelihood L and the prior parameter probability distribtution Π: P (y) =

Z

P (y|xn)P (xn)dNxn=

Z

LnΠndNxn= Z. (52)

This integral may be challenging to compute depending on the number of parameters N and the parameter range. The posterior probability distribution can be calculated nevertheless by combining the posterior probability distribution of the iteration step n + 1 and the step n, such that

P (xn+1|y)

P (xn_|y) =

P (y|xn+1)P (xn+1)

P (y|xn_{)P (x}n₎ , (53)

and we can eliminate P (y). The next step n + 1 is randomly chosen using for instance a Gaussian distribution around the current step n with a "jump-length" σ. The new step is accepted in case of the probability p is greater than a random number m between 0 and 1:

p = minnP (x

n+1_|y)

P (xn_|y) , 1

o

> m, m ∈ Uniform[0, 1]. (54)

This is called the Markov Chain Monte Carlo (MCMC) algorithm, under the Metropolis-Hasting

algo-rithm to generate the iteration steps xn. The method depends highly on the "jumping-length" σ and

must be chosen such that the MCMC actually probes the parameter space towards higher values of the Likelihood function with increasing n. On the other hand σ needs to be chosen big enough such that

a reasonable convergence speed can be reached (∼ 105 evaluations according to Benneke and Seager

(2012)). To omit sampling of only one local minimum while there exist several, Benneke and Seager

(2012) use the method of parallel tempering (Gregory, 2005). By replacing P (y|xn) with a tempered

distribution profile P (y|xn)β and 0 < β < 1. Thereby, one can flatten the likelihood function and

probe a much wider parameter space. Since the tempered profile leads to a higher acceptance of steps

xn+1, the parameter space gets explored on a wider scale. β must be chosen as reasonalbe ladder βk

increasing towards β = 1, which represents the original likelihood function. Modifying (54) to

p = minnP (x

n+1_{|y, β}k_{)P (x}n_{|y, β}k+1₎

P (xn_{|y, β}k_{)P (x}n+1_{|y, β}k+1₎, 1

o

> m, m ∈ Uniform[0, 1]. (55)

enables the cross influence of iteration step acceptance and the β-ladder, also called the "temperatures", to run in parallel. The behaviour of this modified MCMC algorithm depends strongly on the choice of

the βk’s and must be chosen by probing. Additionally one can run several simulations with different

initial parameter distributions x0. A proof why the Metropolis-Hastings works can be found in (Gregory,

2005).

2.3 Planetary atmospheric models

(22)

radiative transfer, cloud physics, photochemistry, non-equilibrium chemistry, atmospheric escape, and many other processes, and calculating transmission spectra thereof. Following from that, the stellar input as energy input into the atmosphere as well as the light source to create transmission spectra would have to be as accurately and precisely modelled as the atmosphere of the transiting exoplanet itself. The reality shows that such models, at least for now, are not efficient to answer present scientific questions. Additionally, the instruments and telescopes available or under development could not satisfy the same need of accuracy in data retrieval to actually justify the effort in producing exclusively such models. Obviously, for this reason a great variety of ideas how to model exoplanetary atmospheres and how to retrieve transmission spectra have emerged and are specialized on certain science cases. One needs to be very careful in choosing which model to use to interpret their data. Naturally, that has

lead to conflicts between retrieval methods and models, applied to the same data sets. See (de Kok

et al. (2013) Discussion) and (Sotzen et al.(2019) Conclusion, and literature therein).

Historically, the first planetary atmospheric models were adapted from atmospheric models of brown

dwarfs (Burrows et al., 1997) or irradiated stellar binary atmospheric models (Seager and Sasselov,

1998). Those were used to investigate close in extraSolar giant planets, highly irradiated by their host

stars. The atmospheric models were based on the assumptions of local thermodynamic equilibrium (LTE) and convection was treated using the method of the mixing length theory. The radiative transfer

equation is solved using Feautrier method or two-stream approximation (Gandhi and Madhusudhan,

2017;Heng et al.,2011a,2014;Mihalas,1978). Seager and Sasselov (1998) used an equation of state consisting of the elements up to two ionization stages H, He, C, N, O, Ne, Na, Mg, Al, Si, S, K, Ca,

Cr, Fe, Ni, and Ti; the ions H₂, H; and the molecules H₂, H₂O, OH, CH, CO, CN, C₂, N₂, O₂, NO,

NH, C₂H₂, HCN, C₂H, HS, SiH, C₃H, C₃, CS, SiC, SiC₂, NS, SiN, SiO, SO, S₂, SiS, and TiO.

Equation of states describe the thermodynamic balance of the different constituents in the atmosphere and can be encountered in two ways. The type of E. o. S. we might be more familiar with, is the

mechanical equation of state (Mihalas,1978):

P = P (ρ, T, chemical composition), (56)

and accounts for the relation between pressure, temperature and density, whereas density can be related to the sum of the densities of the different constituents. The other type important in astrophysical simulations, especially allowing hydrodynamical motion, is the caloric equation of state:

e = e(ρ, T, chemical composition), (57)

and accounts for the connection between the internal energy of the gas to the temperature and again the chemical constituents, for instance through the specific heat coefficients. In most cases for the equation of state we simply use the ideal gas law. However, this might not hold for instance in non-LTE cases or in case of adiabatic or dissipative heating or cooling and we need the caloric equation of state as well. Nevertheless, we will see that the ideal gas law is almost always valid for the internal processes in a planetary atmosphere.

Seager and Sasselov(1998) also included bound-free and free-free atomic transition opacities, Thomson

scattering, and Rayleigh scattering by H₂ and H, and straight means opacities of H₂O and TiO. The

(23)

locally. The incident flux is treated in the full manner of radiative transfer, computing the heat contribution from bound-bound, bound-free, and free-free transitions in different layers. Scattering does not add to the heat budget, however changes the radiation field towards a black body spectrum. The entropy conservation describes the redistribution of heat in the atmosphere. A parameter f describing a qualitative redistribution is defined such that f = 1 complete redistribution of heat and f = 2, night and day side are thermodynamically completely disconnected. Additionally, they incorporate the influence of dust as an extra parameter, without generally regarding processes of dust grain formation. The influence of scattering through dust is treated using Mie theory, to compute scattering opacities

averaged over all angles. However, Seager and Sasselov (1998) do not treat cloud opacities, cloud

formation, or winds in exoplanetary atmospheres. After Seager et al. (2000) found high dependency

of the pressure temperature profile (P-T) and the emergent flux from the cloud formation depth, and

the importance of the alkali metal Na I and K I resonance doublets and the He I 23S - 23P triplet,

Seager and Sasselov (2000) update the code from Seager and Sasselov (1998) with Gibbs free energy minimization to calculate solids and gases in chemical equilibrium and condensate opacities for three solid species. This gives rise to cloud formation with cloud top and cloud base of different condensation species.

Complementary, others proposed other ways of modelling exoplanetary atmospheres, and focus on

certain aspects differently. Brown (2001) for instance created a non-strict self-consistent model to

describe the terminator region of an exoplanetary atmosphere. Thereby, restricting himself in the number of chemical species or the calculation of He opacities. However, he includes photoionization of the alkali metals Na and K, the dynamics of emergent horizontal winds due to the one sided heating at different incident angles, and differential rotation of the planet. Another model referred to before

called LinePak, developed byGordley et al. (1994) treats radiative transfer under similar assumptions

as Seager and Sasselov (1998), however does not treat scattering, clouds, or dust. The atmospheric model can be generated in a customized way with different layers with different mass functions and

accounts for refraction effects cell by cell of light rays passing through. At the same timeFortney et al.

(2008,2010) and references therein present their way of adopting and testing models for brown dwarfs

as proxys for highly irradiated hot Jupiters. Gandhi and Madhusudhan (2017) review some of the

work done until 2017. A table extracted from their work on the new self-consistent planetary model

(24)

Figure 5: The table is copied fromGandhi and Madhusudhan(2017) and presents the previous self consistent models and their major differences considering the treatment of convection, clouds and scattering. Some of the models presented here we have also covered in our text.

In this chapter we will give an outline of the different types of atmospheric models and present the most important physical inputs and assumptions these models are based on.

2.3.1 Parametric forward models

Madhusudhan and Seager (2009) present a parametric forward model based on a simple atmospheric three-layer model with holistic assumptions for the general structure of the P-T profile, and chemical composition. The model is generated with 10 parameters in a computationally fast and simple manner,

and has for instance been applied byde Kok et al.(2013). These kind of models are called parametric

forward models, since they are generated based on parameters rather than physics. The method was designed for hydrogen-dominated hot Jupiter atmospheres. Each layer is described through a pressure-temperature relation approximated by

P = P0eα(T −T0)

β

(58)

with the free parameters P0, T0, α and β for each layer. Layer number three is assumed to be convective

and has a constant temperature T3. Eliminating two more parameters at the boundary of layer 2-3 and

1-2, 6 parameters remain to describe the P-T profile. The approach of a three-layer parametric P-T profile is motivated by studies of the atmosphere of Solar System bodies and self-consistent planetary

(25)

Figure 6: General temperature and pressure profile presented by Madhusudhan and Seager (2009) that can be adjusted by parameters for each layer. Layer 3 is isothermal. Notice that the profile allows temperature inversions to occur.

Al-though it is a completely parameter dependent model, Madhusudhan and Seager

(2009) enforce energy balance at the top of the atmosphere and discard any

pa-rameter combinations violating this requirement.

The radiative transfer through the atmosphere is solved using a line-by-line code under the assumption of LTE, no scattering, and the requirement of energy balance at the top of the atmosphere. The density is given by the ideal gas law from the P-T profile and the radial dependence is determined by hydrostatic equilibrium. Stellar heating is incorporated into the model through the P-T profile, which in turn is indirectly dependent on the incident, absorbed, and scattered radiation. An important point to mention is that the emergent transmission spectrum in the infrared is decoupled from the opacities in the visual, responsible for some of the heating. Species in planetary atmospheres can have strong opacity contributions in the different wavelength regimes. Therefore, a parametric P-T profile can account for the overall heating, treating the underlying physics as a black box. On top of that the partial independence between the P-T profile and atmospheric composition allows for probing P-T profiles of atmospheres with loosely known chemical composition and does not need to be described self-consistently with the constituents in chemical equilibrium. This is called the dual-band approximation (Guillot,2010;Heng et al.,2011a), which we will explain more later.

Molecules included in the model are H₂, H₂O, CO, CO₂, CH₄, and NH₃, and it also incorporate

H2 - H2 collision-induced absorption opacities. CO2 and NH3 are chosen to have a fiducial arbitrary

concentration of Solar abundance, since NH3has only minor absorption features in the bands considered

(26)

The concentrations and radiative transfer is calculated in 100 layers uniformely spaced in log P . The molecular abundances are controlled by the ratio between the concentration and a fiducial equilibrium

concentration, equal to Solar abundances, of a molecule x in layer i: fx = cix/cix,eq. Although the

con-centration of each molecule is different for each layer, dividing by the fiducial equilibrium concon-centration

they get a constant abundance parameter fx. The abundance parameter for the main constituents

H2O, CO, CH4, and CO2 are free parameters of the model. Thus, the parametric forward model of

Madhusudhan and Seager(2009) comes with a total of 10 parameters, 6 parameters for the P-T profile and 4 parameters for the abundance of the main constituents of the atmosphere.

Clouds are not treated in this version, however, were planned to be included in future versions. Since

Madhusudhan and Seager(2009) focused on the study of hot Jupiters and previous literature suggested them to be cloud free, they had no need to include treatment of clouds initially. This assumption was

shown to be invalid in more recent transmission observations in 2013, described inBarstow et al.(2014)

and literature therein.

The energy balance is controlled via the incoming wavelength-integrated flux from the star F? and the

emitted flux from the top of the atmosphere F_p. It is scaled via the Bond albedo A_B, describing the

ratio of the incident radiation flux reflected back into space, and a redistribution factor f_r accounting

for the redistribution of heat to the night side:

Fp= (1 − AB)(1 − fr)F?, (59)

and the actual energy balance is constrained by η = (1 − AB)(1 − fr) ≤ 1. To use the model, they

propose to calculate a goodness-of-fit dependent on the difference between the model calculated flux and the observed flux at each wavelength. The goodness-of-fit is then computed for each model on the

parameter grid of about 107 grid points.

Guillot(2010) followed a different approach compared toMadhusudhan and Seager(2009) and derived a T − τ model using the dual band approximation under the assumption of constant opacities in two wavelength bands. The radiation contributions are treated again as independent distinct bands of incoming shortwave or visual and outgoing longwave or thermal radiation. The radiative transfer

equation to solve is (Heng et al.,2011a)

µ∂Iν ∂m = κνIν ξ − κνBν− κν(1 − ξ)Jν ξ . (60)

The parameter ξ carries information about scattering, which is not treated byGuillot(2010).

ξ will be taken up again in section 2.3.2. Using the three moments of the specific intensity (Guillot,

2010) (Jν, Hν, Kν) ≡ 1 2 Z 1 −1 Iµν(1, µ, µ2)dµ (61)

the radiative transfer equation (60) can be recast in the form:

dHν

dm = κν(Jν − Bν) (62)

dKν

dm = κνHν, (63)

where m is the column mass and carries the information about the height profile z, through dm =

(27)

with the flux, 4π/cK_ν represents the radiation pressure, and B_ν is the Planck function or Black-body radiation. The upper boundary condition is given by

Hv(0) = −µ?Jv(0), (64)

for the irradiation from the star. µ_? = cos θ?, where θ? denotes the angle of the incoming radiation

relative to the direction perpendicular to the atmosphere, and Hν is evaluated at m = 0. Guillot(2010)

introduces the parameter γ ≡ κv

κth, which is the ratio between the mean opacity in the visual range and

the mean opacity of the infrared or thermal radiation. The mean opacities are defined as the average

taken over the mean intensity Jν in the respective band:

κv = Jv−1 Z visual κνJνdν (65) κth= J_th−1 Z thermal κνJνdν. (66)

However, we can approximate κ_v by using the Planck function B_ν for the mean intensity J_ν at the

temperature Tirr, which is the temperature of the radiation field from the star, received at the top of

the atmosphere. Then the upper boundary condition becomes

Hv(0) = −µ?Jv = −

µ?

4πσT

4

irr (67)

The irradiation temperature can be computed by relating the stellar flux with the received stellar flux

at the substellar point µ_? = 0. Using the Stefan-Boltzmann law

σT_irr4 = 4π Z ∞ 0 Hν(0)dν = L? 4πa2 = σT_?4R2_? a2 (68) ⇒ T_irr = T? R? a 1/2 , (69)

where a is the average distance between the star and the planet and σ the Stefan-Boltzmann constant.

Tirr can be related to the equilibrium temperature for complete redistribution of energy around the

planet: σT_eq4 = σT 4 irrπR2p 4πR2 p ⇒ T_eq= T√irr 2 = T? R? 2a 1/2 . (70)

More accurately Teq would be dependent on the Bond albedo AB, carrying the information about the

efficiency of the radiative redistribution of energy around the planet. κth is averaged the same way as

in the visual band, but using the equilibrium temperature of the planet T_eq for the Planck function

Bν. Thereby, the mean opacities can be calculated from first principles and do not need to be iterated

by computing Jν. Generally, the calculation of the equilibrium temperature is also dependent on the

reflected energy back into space and we can not assume the complete uniform redestribution throughout the entire atmosphere. The Bond albedo and its relation to scattering properties in the atmosphere is

treated later in section 2.3.2as well.

The argument for the approximation of the mean opacities from above is that the contribution in the visual is dominant at low τ 1 (optically thin), where the visual radiation field is dominated by

the incoming stellar radiation, controlled by Tirr, and absorption effects are sparse. In the region of

(28)

radiation, which can be approximated by B_ν at T_eq. This of course tailors the crude assumption of

Bν ∼ 0 for ν in the visual, at high τ 1. Vice-versa the longwave stellar contribution to the radiation

field must be negligible, which might not hold for instance for M-dwarfs. Similarly the moments of the specific intensity are integrated over the desired wavelength range in the visual v and thermal

th regime. In the visual regime, following Eddington (Eddington, 1916) µ? can be approximated by

µ2_? = Kv

Jv, which is µ

2

? = 13 in case of an isotropic incoming radiation field. Thus, the radiative transfer

equation in the visual becomes

d2_(H v, Jv) dm2 = κ2 v µ2 ? (Hv, Jv) (71)

and Bv, the thermal emission of the atmosphere in the visual at m can be neglected. A general solution

for the visual band is then

(J (m)v, Hv(m)) = (Jv(0), Hv(0)) exp (−κvm/µ?). (72)

The optical depth τ is defined as the thermal opacity κ_th multiplied with m

τ = κthm. (73)

Guillot(2010) uses the first and second Eddington coefficients

fKth = Kth Jth , fHth = Hth(0) Jth(0) , (74)

where Hth(0) means that Hth is evaluated at the outer boundary τ = 0 of the atmosphere. The set of

equations for the thermal regime is then given by

dHth

dm = κth(Jth− B) (75)

dKth

dm = κthHth (76)

κth(Jth− B) + κvJv = 0 (77)

with B =R_thermalBνdν. We can see through the radiative equilibrium condition (77) how the thermal

and the visual band connect. The parameters of the Guillot (2010) profile are Tint, Teq, γ ≡ _κκv

th, fKth,

and f_H_th. The equations (75), (76), and (77) can be combined and integrated in m and using the

appropriate boundary conditions, such as

Z ∞

0

Hνdν = H, (78)

Hth(0) = H + µ?Jv(0), (79)

and (64), to express the source function B in terms of H, Hv(0), fHth, fKth, γ, µ?, and τ :

B = H 1 fHth + τ fKth − H_v(0) 1 fHth + µ? γfKth + γ µ? − µ? γfKth exp (−γτ /µ?) . (80)

H can be identified with H = σT_int4 , with Tint the temperature associated with the deeper layers, and

internal heating and radiative equilibrium matches the sum of the flux in the visual and the thermal

Uppsala University