# Terminus: A Versatile Simulator for Space-based Telescopes

Billy Edwards Blue Skies Space Ltd., 69 Wilson Street, London, EC2A 2BB, UK
Department of Physics and Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK  Ian Stotesbury  Blue Skies Space Ltd., 69 Wilson Street,
London, EC2A 2BB, UK

###### Abstract

Space-based telescopes offer unparalleled opportunities for characterising
exoplanets, Solar System bodies and stellar objects. However, observatories in
low Earth orbits (e.g. Hubble, CHEOPS, Twinkle and an ever increasing number
of cubesats) cannot always be continuously pointed at a target due to Earth
obscuration. For exoplanet observations consisting of transit, or eclipse,
spectroscopy this causes gaps in the light curve, which reduces the
information content and can diminish the science return of the observation.
Terminus, a time-domain simulator, has been developed to model the occurrence
of these gaps to predict the potential impact on future observations. The
simulator is capable of radiometrically modelling exoplanet observations as
well as producing light curves and spectra. Here, Terminus is baselined on the
Twinkle mission but the model can be adapted for any space-based telescope and
is especially applicable to those in a low-Earth orbit. Terminus also has the
capability to model observations of other targets such as asteroids or brown
dwarfs.

##  1 Introduction

To date, several thousand extra-solar planets have been discovered. With many
of these now being detected around bright stars, and with many more to come
from missions such as the Transiting Exoplanet Survey Satellite (TESS, Ricker
et al. (2014); Barclay et al. (2018)), the characterisation of these worlds
has begun and will accelerate over the next decade. Ground-based instruments
have detected absorption and emission lines in exoplanet atmospheres via high
resolution spectra (e.g. Hoeijmakers et al., 2018; Ehrenreich et al., 2020))
while the Hubble and Spitzer space telescopes have used lower resolution
spectroscopy or photometry to probe the chemical abundances and thermal
properties of tens of planets (e.g. Sing et al., 2016; Iyer et al., 2016;
Tsiaras et al., 2018; Garhart et al., 2020).

In the coming years several missions, some of which are specifically designed
for exoplanet research, will be launched to provide further characterisation.
While the James Webb Space Telescope (JWST, Greene et al. (2016)) and Ariel
(Tinetti et al., 2018) will be located at L2, observatories such as the
CHaracterising ExOPlanets Satellite (CHEOPS, Benz et al. (2020)), which was
launched in December 2019, and Twinkle (Edwards et al., 2019d) will operate
from a low Earth orbit and as such will have to contend with Earth
obscuration.

The orbit will cause gaps in some of the observations obtained by these
missions which will impact their information content due to parts of the
transit light curve being missed, decreasing the precision of the recovered
transit parameters. Additionally, the thermal environment of a low Earth orbit
and the breaks in observing can lead to recurring systematic trends such as
ramps in the recorded flux due to thermal breathing of the telescope and
detector persistence. Such gaps and systematics are experienced in all
exoplanet observations with Hubble (e.g. Deming et al., 2013; Kreidberg et
al., 2014). It should be noted, however, that Hubble is situated in an
equatorial orbit which is significantly different to the sun-synchronous
orbits of CHEOPS and Twinkle. Sun-synchronous orbits allow for certain areas
of sky, specifically those in the anti-sun direction, to be observed for
longer periods without interruption. Additionally, the thermal environment is
more stable due to the smaller variations in the spacecraft-Earth-Sun
geometry. Previous missions to have operated in sun-synchronous orbits include
the Convection, Rotation and planetary Transits (CoRoT, Bordé et al. (2003)),
Akari (Murakami et al., 2007) and WISE/NEOWISE (Wright et al., 2010; Mainzer
et al., 2014). Due to it’s Earth-trailing orbit, Spitzer (Werner et al., 2004)
did not experience gaps in its observations.

When designing future instrumentation, understanding the expected performance
for the envisioned science cases is paramount. Static models, often referred
to as radiometric or sensitivity models, are suitable for studying the
instrument performance over a wide parameter space (i.e. for many different
targets) as they are generally quick to run and require relatively minimal
information about the instrumentation. Radiometric models are a useful way to
understand the capabilities of upcoming exoplanet observatories and have been
widely used. The ESA Radiometric Model (ERM, Puig et al. (2015)) was used to
simulate the performance of the ESA M3 candidate EChO (Exoplanet
Characterisation Observatory, Tinetti et al. (2012)) and was subsequently used
for Ariel (Puig et al., 2018). A newer, python-based version, ArielRad, was
recently developed (Mugnai et al., 2020) while PandExo has been created for
simulating exoplanet observations with Hubble and JWST (Batalha et al., 2017)
and the NIRSpec Exoplanet Exposure Time Calculator (NEETC) was built
specifically for modelling transit and eclipse spectroscopy with JWST’s
NIRSpec instrument (Nielsen et al., 2016). These usually account for
efficiency of the optics and simple noise contributions such as photon, dark
current, readout and instrument/telescope emission.

More complex effects, such as jitter, stellar variability and spots and
correlated noise sources require models which have a time-domain aspect. These
tools usually also produce simulated detector images which can act as
realistic data products for the mission, accounting for detector effects such
as correlated noise between pixels or inter- and intra-pixel variations. For
example, ExoSim is a numerical end-to-end simulator of transit spectroscopy
which is currently being utilised for the Ariel mission (Pascale et al., 2015;
Sarkar et al., 2016, 2017). The tool has been created to explore a variety of
signal and noise issues that occur in, and may bias, transit spectroscopy
observations, including instrument systematics and the other effects
previously mentioned. By producing realistic raw data products, the outputs
can also be fed into data reduction pipelines to explore, and remove,
potential biases within them as well as develop new reduction and data
correction methods. End-to-end simulators such as ExoSim are therefore
powerful tools for understanding the capabilities of an instrument design.
Additional time-domain simulators of note include ExoNoodle (Martin-Lagarde et
al., 2019), which utilises MIRISim (Geers et al., 2019) to model time-series
with the JWST MIRI instrument, Wayne which models Hubble spatial scans of
exoplanets (Varley et al., 2017) and the simulators developed for the CHEOPS
and Colorado Ultraviolet Transit Experiment (CUTE) missions (Futyan et al.,
2020; Sreejith et al., 2019). While the complexity of these types of tools can
be hugely advantageous in understanding intricate effects it can also be their
biggest weakness; such sophisticated models require a great deal of time to
develop and run as well as an excellent understanding of all parts of the
instrument design. They can therefore only be applied to highly refined
designs and run for a small number of cases. The solution to the issue of
complexity versus efficiency is to use both types of models. For Ariel, ExoSim
is used to validate the outcomes of ArielRad for selected, representative
targets. ArielRad is then used as the workhorse for modelling the capability
of thousands of targets due to its superior speed (Edwards et al., 2019b;
Mugnai et al., 2020).

Here, we describe the Terminus tool which has been developed to model transit
(and eclipse) observations with Twinkle, to explore the impact of Earth
obscuration and allow for efficient scheduling methods to be developed to
minimise this impact. The simulator, however, is not mission specific and
could be adapted for other observatories, with a particular applicability for
satellites in low Earth orbit.

The Twinkle Space Mission111http://www.twinkle-spacemission.co.uk is a new,
fast-track satellite designed to begin science operations in 2024. It has been
conceived for providing faster access to spectroscopic data from exoplanet
atmospheres and Solar System bodies. Twinkle is equipped with a visible and
infrared spectrometer which simultaneously covers 0.5-4.5 $\mu$m with a
resolving power of R$\sim$20-70 across this range. Twinkle has been designed
with a telescope aperture of 0.45 m. Twinkle’s field of regard is a cone with
an opening angle of 40∘, centred on the anti-sun vector (Savini et al., 2016).

Previously the ESA Radiometric Model (ERM, Puig et al. (2015, 2018)), which
assumes full light curves are observed, has been used to model the
capabilities of Twinkle (see Edwards et al. (2019d)). Terminus includes a
radiometric model, built upon the concepts of the ERM, but it has been
upgraded to also have the capacity to simulate light curves. The code also
contains the ability to model the orbit of a spacecraft, thus allowing for the
availability of targets to be understood given solar, lunar and Earth
exclusion angles. The capability to model these gaps is not available in other
tools such as ArielRad or ExoSim and is one of the driving factors behind the
creation of Terminus. Additionally, the Twinkle mission will not be limited to
exoplanet characterisation and will also observe solar system bodies, brown
dwarfs and other astrophysical objects. As such, Terminus builds upon the work
of Edwards et al. (2019a, c) and can be used to calculate the predicted data
quality and observational periods for these objects, another feature which is
not present in other similar codes.

In this work we first describe the portion of the simulator which calculates
the target signal and noise contributions before comparing the outputs of
simulated light curve fitting to radiometric estimates. Next the orbital
module is detailed and validated against outputs from an orbital dynamics
software. Using this we explore the effect of gaps for observations of HD
209458 b and WASP-127 b with Twinkle. Finally, we discuss Twinkle’s ability to
observe asteroids by focusing on potential observations of the Near-Earth
Object (NEO) 99942 Apophis (2004 MN4).

##  2 Simulator Structure

Terminus has been constructed in Python and has several different stages. It
can be operated as a simple radiometric model, used to calculate expected
signal-to-noise ratio (SNR) on a given number of atmospheric scale heights, or
be utilised to create simulated light curves. A instrument file is loaded
(which includes parameters such as telescope aperture, quantum efficiency
etc.) and the star flux on the detector calculated. PSFs can be imported from
external sources. In Sections 2.1 to 2.4 we discuss the structure of the
simulator and an overview is given in Figure 1.

###  2.1 Target Parameters

A catalogue of planets has been created following the methodology of (Edwards
et al., 2019b) and data is taken from the NASA Exoplanet Archive (Akeson et
al., 2013).

Figure 1: Overview of the simulator structure. Generic parts are represented
by blue shapes while red indicates functions which are exoplanet specific.
Dotted lines indicate portions which are not compulsory.

###  2.2 Radiometric Model

Figure 2: Example detector images generated by Terminus for Twinkle Ch0 (top)
and Ch1 (bottom). These are used purely for the calculation of the saturation
time for each target.

The stellar flux at Earth is calculated using spectral energy distributions
(SEDs) from the PHOENIX BT-Settl models by Allard et al. (2012); Husser et al.
(2013); Baraffe et al. (2015). The spectral irradiance from a host star at the
aperture of the telescope is given by:

$E_{S}(\lambda)=S_{S}(\lambda)(\frac{R_{*}}{d})^{2}$ (1)

where S${}_{\rm S}(\lambda)$ is the star spectral irradiance from the Phoenix
catalogue (Wm${}^{-2}\mu$m-1) and d is the distance to the star. The effective
collecting area of the telescope is then accounted for before the flux is
integrated into the spectral bins of the instrumentation to give a photon flux
per bin. The signal is then propagated through the instrument to the detector
focal planes, taking into account the transmission of each optical component
and diffracting element as well as the quantum efficiency of the detector. The
final signal, in electrons per second, from the star in each spectral bin is
determined as a 1D flux rate before being convolved with 2D point spread
functions (PSFs) and the instrument dispersion to create a detector image. The
detector image, like the one shown in Figure 2, is utilised to calculate the
saturation time for the target while the 1D flux rate is used for all other
calculations. A variety of sources of noise are accounted for in each of the
models. In addition to photon noise, the simulator calculates the
contributions from dark current, instrument and telescope emission, zodiacal
background emission, and readout noise. Additionally, photometric
uncertainties due to spacecraft jitter can be imported and interpolated from
time-domain simulators such as ExoSim (see Section 2.2.3). Some of these noise
sources are wavelength dependent (e.g. zodiacal background) while others are
not (e.g. read noise).

####  2.2.1 Calculating Noise Per Exposure

In describing the acquisition of data we use the nomenclature of Rauscher et
al. (2007) in which a frame refers to the clocking and digitisation of pixels
within a specified area of the detector known as a sub-array. The size of sub-
array dictates the time required for it to read out. Here, given the footprint
of Twinkle’s spectrometer on the detector, we assume a fastest frame time of
0.5 seconds which is similar to that for the 1024A sub-array on JWST NIRSpec
(0.45 seconds, Pontoppidan et al. (2016)). A collection of frames then forms a
group although here, as with JWST time series, the number of frames is set to
one (i.e. tg = tf). A collection of non-destructively read groups, along with
a detector reset, forms an integration. Here, the detector reset time after a
destructive read is also assumed to be equivalent to the frame time. As the
duration of a transit/eclipse is generally orders of magnitude longer than the
saturation time of the detector, many integrations will be taken during an
observation. The total noise variance per integration, $\sigma^{2}_{\rm exp}$,
is given by:

$\sigma^{2}_{exp}=\frac{12(n_{g}-1)}{n_{g}(n_{g}+1)}n_{pix}\sigma^{2}_{read}+\frac{6(n_{g}^{2}+1)}{5n_{g}(n_{g}+1)}(n_{g}-1)t_{g}i_{total}$
(2)

from Rauscher et al. (2007) where ng is the number of groups (non-destructive
reads) per exposure, $\sigma_{\rm read}$ is the read noise in e-/pix rms, npix
is the number of pixels in the spectral bin, tg is the time for a single non-
destructive group read, and itotal is the total flux in e-/s. For JWST
observations, the standard practise for exoplanet observations is to maximise
the number of groups (Batalha et al., 2017). Meanwhile, Ariel will use a
variety of readout modes, depending upon the brightness of the target, with
correlated double sampling (CDS, ng =2) for brighter sources targets and
multiple up-the-ramp reads for fainter targets (Focardi et al., 2018).
Collecting several up-the-ramp reads can be useful in correcting for cosmic
ray impacts while also reducing the read noise. Additional reads, however,
increase the photon noise contribution and thus Terminus varies the number of
up-the-ramp reads according to the brightness of the target to attempt to
optimise noise. In each case, the maximum number of up-the-ramp reads is
calculated and Equation 2 used to selected the number of reads which yields
the lowest noise per transit observation (using Equations 2-7). npix can be
selected by specifying a required encircled energy but when importing jitter
simulations from ExoSim, npix is set to the values used in these simulations
as outlined in Section 2.2.3. In Equation 2, itotal is defined as:

$i_{total}=i_{sig}+n_{pix}(i_{dark}+i_{bdg})$ (3)

where isig is the total signal from the star in the spectral bin (e-/s) while
idark and ibdg are the dark current and background signals respectively (in
e-/s/pix). Currently, $i_{\rm bdg}$ is assumed to be from the emission of
optical elements and Zodiacal emission, as detailed in Section 2.2.2, but
future updates will include contributions from nearby stars. For exoplanet
spectroscopy, the total observational time is generally quantised in terms of
the duration of a transit/eclipse event, T. The model assumes the time spent
during ingress (T12) and egress (T34) is negligible to the primary transit
time (T23) and thus T = T23 = T14. The transit time can be calculated from:

$T_{14}=\sqrt{1-b^{2}}\frac{R_{*}P}{\pi a}$ (4)

for a given system where P is the orbital period. The fractional noise on the
star signal over one transit duration is then given by:

$\sigma_{Star}=\frac{1}{\sqrt{n_{int}}}\frac{\sigma_{exp}}{i_{sig}}$ (5)

where nint is the number of integrations over one transit duration which is
calculated from:

$n_{int}=\frac{T_{14}}{t_{r}+t_{g}n_{g}}$ (6)

where tr is the time taken to reset the detector. As a baseline we take tr to
be equivalent to the frame time, tf (0.5 seconds). As noted by Batalha et al.
(2017), if tg = tr = tf then the duty cycle (i.e. the efficiency) is given by
(ng $-$ 1)/(ng $+$ 1).

The measurement of the transit depth is differential and thus the error (i.e.
the 1$\sigma$ uncertainty) on the transit depth is given by:

$\sigma_{TD}=\sigma_{Star}\sqrt{1+\frac{1}{n_{{}_{T_{14}}}}}$ (7)

where n${}_{\rm T_{\rm 14}}$ is the number of transit durations observed out
of transit (i.e. the baseline). For all simulations presented here,
$n_{T_{14}}$ is set to 2 (i.e. 1 x T14 is spent both before/after the main
observation). The error is calculated in this way for every spectral bin.

####  2.2.2 Zodiacal Emission

We calculate the contribution of zodiacal emission using the prescription from
Pascale et al. (2015) and Sarkar et al. (2020a). The signal is composed of two
black bodies, with associated coefficients, to model the reflected and emitted
components. The spectral brightness is given by:

$\displaystyle Zodi(\lambda)=\beta(3.5\times 10^{-14}B_{\lambda}(5500K)$ (8)
$\displaystyle+3.58\times 10^{-8}B_{\lambda}(270K))$

where the coefficient $\beta$ modifies the intensity of the zodiacal light
based upon the declination of the target. At the ecliptic poles, $\beta$ = 1
provides a good fit to the intensity shown in Leinert et al. (1998). Sarkar et
al. (2020a) fitted a polynomial to data from this study, along with zodiacal
intensities from James et al. (1997); Tsumura et al. (2010), to provide a
measure of the increase in intensity at different latitudes. If d is the
ecliptic latitude, then the coefficient is given by:

$\displaystyle\beta=-0.22968868\zeta^{7}+1.12162927\zeta^{6}-1.72338015\zeta^{5}$
(9) $\displaystyle+1.13119022\zeta^{4}-0.95684987\zeta^{3}+0.2199208\zeta^{2}$
$\displaystyle-0.05989941\zeta+2.57035947$

where $\zeta$ = $log_{10}(d+1)$. This relation falls below 1 at d = 57.355 ∘
and so $\beta$ is fixed to 1 for latitudes greater than this (Sarkar et al.,
2020a).

####  2.2.3 Pointing Jitter

Directly modelling uncertainties due to spacecraft jitter is beyond the
capabilities of Terminus. Hence, ExoSim has been adapted to the Twinkle design
to study the effects of pointing jitter on science performance. ExoSim, first
conceived for EChO Pascale et al. (2015) and now used for the Ariel mission,
has previously been adapted for simulating observations with JWST (Sarkar et
al., 2020a) and the EXoplanet Climate Infrared TElescope (EXCITE, Tucker et
al. (2018); Nagler et al. (2019)). The modified version, christened
TwinkleSim, was run for a number of stellar types (TS = 3000, 5000, 6100 K)
and magnitudes (KS = 6, 9, 12) and the uncertainty due to jitter determined in
each case. Twinkle’s baseline pointing solution is based upon a high
performance gyroscope and a Power Spectral Density (PSD) was supplied by the
engineering team at the satellite manufacturer, Airbus. For each simulation, a
variety of different extraction apertures were trialled with larger apertures
reducing the jitter by ensuring clipping did not occur but increasing the
noise from other sources due to sampling more pixels (e.g. dark current).
After trialling a number of solutions, the aperture was set to be rectangular
with a width of 2.44 times the width of the Airy disk at longest wavelength of
each channel. In terms of pixels, this is equivalent to 12 and 22 in the
spatial direction, for Ch0 and Ch1 respectively, while the spectral pixels per
bin are set to 6 and 7.

When combining observations time-correlated noise may integrate down more
slowly than uncorrelated noise, which is assumed to decrease with the square
root of the number of observations, and thus can contribute more heavily to
the final noise budget. To account for this Allan deviations plots were
produced using TwinkleSim. A power law trend can be fitted to this and used to
derive a wavelength-dependent fractional noise term that jitter induces on the
photon noise. For more details on this process, we refer the reader to Sarkar
et al. (2020a).

####  2.2.4 Transit Signal

During transit, the critical signal is the fraction of stellar light that
passes through the atmosphere of the exoplanet. This signal is determined by
the ratio of the projected area of the atmosphere to that of the stellar disk
and thus is given by:

$\frac{2R_{p}\Delta z(\lambda)}{R_{*}^{2}}$ (10)

where $\Delta z$ is the height of the atmosphere. The size of the atmosphere
is taken to be equivalent to the height above the the 10 bar radius, at which
point the atmosphere is assumed to be opaque. The pressure of an atmosphere at
a height, z, is given by:

$p(z)=p_{0}e^{\frac{-z}{H}}$ (11)

where H is the scale height, the distance over which the pressure falls by
1/e. In the literature, 5 scale heights are often assumed for $\Delta z$ for a
clear atmosphere (at which point one is above 99.5$\%$ of the atmosphere)
while 3 would be more reasonable in the moderately cloudy case (Puig et al.,
2015; Tinetti et al., 2018; Edwards et al., 2019b). The scale height of the
atmosphere is calculated from:

$H=\frac{kT_{p}N_{A}}{\mu g}$ (12)

where k is the Boltzmann constant, $N_{A}$ is Avogadro’s number, $\mu$ is the
mean molecular weight of the atmosphere and g is the surface gravity
determined from:

$g=\frac{GM_{p}}{R_{p}^{2}}$ (13)

where $M_{p}$ and $R_{p}$ are the mass and radius of the planet and G is the
gravitational constant.

####  2.2.5 Eclipse Signal

During eclipse, the signal is calculated form two sources; reflected and
emitted light from the planet. Emission from the exoplanet day-side is
modelled as a black body and the wavelength-dependent surface flux density is
given by:

$S_{p}(\lambda,T_{p})=\pi\frac{2hc^{2}}{\lambda^{5}}\frac{1}{e^{\frac{hc}{\lambda
kT_{p}}}-1}$ (14)

where $T_{p}$ is the dayside temperature of the planet. The product of the
black body emission and the solid angle subtended by the exoplanet at the
telescope gives the spectral radiance at the aperture:

$E_{p}^{Emission}(\lambda,T_{p})=S_{p}(\lambda,T_{p})\left(\frac{R_{p}}{d}\right)^{2}$
(15)

in $Wm^{-2}\mu m^{-1}$. Additionally, a portion of the stellar light incident
on the exoplanet is reflected. The strength of this reflected signal is
strongly dependant on wavelength and can be significant at visible
wavelengths. The flux of reflected light at the telescope aperture is
calculated from:

$E_{p}^{Reflection}(\lambda)=\alpha_{geom}S_{s}(\lambda)\left(\frac{R_{*}}{d}\right)^{2}\left(\frac{R_{p}}{a}\right)^{2}$
(16)

where $S_{S}(\lambda)$ is the star spectral irradiance, $a$ is the star-planet
distance (i.e. the planet’s semi-major axis) and $\alpha_{geom}$ is the
geometric albedo, which is assumed to be that of a Lambertian sphere
($\frac{2}{3}\alpha_{bond}$), wavelength-independent and at a phase of $\phi$
= 1 (i.e. full disk illumination).

####  2.2.6 Signal-to-Noise Ratio

From these equations, and the error on the transit/eclipse depth, the signal-
to-noise (SNR) on the atmospheric signal can be obtained for a single
observation. Assuming the SNR increase with the square root of the number of
observations, the SNR after multiple transits/eclipses is given by:

$SNR_{N}=\sqrt{N}SNR_{1}$ (17)

where SNR1 is the SNR of a single observation and N is the total number of
observations. By setting a requirement on the SNR (SNRR), the number of
observations needed for a given planet can be ascertained from:

$N=\left(\frac{SNR_{R}}{SNR_{1}}\right)^{2}$ (18)

The current requirements are set to a median SNR $>$ 7 across 1.0-4.5 $\mu$m
for transit observations and 1.5-4.5 $\mu$m for eclipse measurements. In the
former of these the shorter wavelengths are excluded to avoid biasing against
planets around cooler stars while the latter is chosen as planetary emission,
even for relatively hot planets ($\sim$1500 K), is low at wavelengths shorter
than 1.5 $\mu$m. Using Equation 18, one can then determine the type(s) of
observation the planet is suited to.

###  2.3 Atmospheric Modelling

To simulate transmission (and emission) forward models, the open-source
exoplanet atmospheric retrieval framework TauREx 3 (Al-Refaie et al., 2019;
Waldmann et al., 2015a, b) is used. Within TauREx 3, cross-section opacities
are calculated from the ExoMol database (Yurchenko & Tennyson, 2012) where
available and from HITEMP (Rothman & Gordon, 2014) and HITRAN (Gordon et al.,
2016) otherwise. The H- ion is included using the procedure outlined in John
(1988); Edwards et al. (2020). For atmospheric chemistry, two options are
available within the Terminus infrastructure: chemical equilibrium, which is
achieved using the ACE code (Venot et al., 2012; Agúndez et al., 2012) and
takes the C/O ratio and metallicity as input, or free-chemistry which allows
the user to choose molecules and their abundances. Alternatively, a high-
resolution spectrum produced by another radiative transfer code can be read in
or, if a retrieval on actual data has been performed, the atmosphere can be
extrapolated from a TauREx 3 hdf5 file. Once the forward model is created at
high resolution, it is then binned to the instrument resolution using TauREx
3’s integrated binning function.

###  2.4 Light Curve Modelling and Fitting

For each spectral bin, PyLightCurve222https://github.com/ucl-
exoplanets/pylightcurve (Tsiaras et al., 2016a) is used to model a noise-free
transit/eclipse of the planet. The transits were all modelled with quadratic
limb darkening coefficients from Claret et al. (2013), calculated using
ExoTETHyS (Morello et al., 2020). The Twinkle spectrometer features a split at
2.43 $\mu$m, creating two channels. For each of these a white light curve is
also generated. The spectral light curves are created at the native resolution
of the instrument (R$\sim$20-70). A time-series is created with a cadence
equal to the time between destructive reads and the light curve integrated
over each of these exposures. The noise per integration, as calculated in
Section 2.2, is then used to create noisy light curves by adding Gaussian
scatter. Further updates will include the ability to add ramps due to detector
persistence as well as other time-varying systematics.

For the fitting of the light curves a Markov chain Monte Carlo (MCMC) is run
using emcee (Foreman-Mackey et al., 2013) via the PyLightCurve package, here
with 150,000 iterations, a burn-in of 100,000, and 100 walkers. For the
simulations shown here, both white light curves are individually fitted with
the inclination (i), reduced semi-major axis (a/R∗), transit mid-time (T0) and
planet-to-star radius ratio (Rp/Rs) as free parameters. A weighted average of
the recovered values for each of these parameters, except the planet-to-star
radius ratio, is then fixed for the fitting of the spectral light curves where
only the planet-to-star radius ratio is fitted. This provides a retrieved
transit/eclipse depth for each light curve, along with the error associated
with this parameter. If further complexity, such as ramps, is added to the
light curve, future iterations of the code will allow for multiple light curve
fits. In this case the uncertainties in the individual data points are
increased such that their median matches the standard deviation of the
residuals, a common technique when analysing Hubble observations of exoplanets
(e.g. Kreidberg et al., 2014; Tsiaras et al., 2016b).

For fainter targets, a spectrum with a reduced resolution can be requested and
Terminus will combine the light curves and provide a spectrum with a
resolution as close to the desired as possible. While the default cadence is
set by the saturation time of the detector it can lowered or exposures can be
combined. Additionally, multiple transits (or eclipses) can be individually
modelled, fitted and then combined. These functionalities are all controlled
by the input configuration file. Once a spectrum has been generated, an
automated interface with TauREx 3 can then be used to fit the data and
retrieve the atmospheric parameters.

To compare the errors predicted by the radiometric model to those from fitted
light curves, we model a single observation of HD 209458 b (Charbonneau et
al., 2000; Henry et al., 2000). For the atmosphere we model a composition
based loosely on that retrieved from the HST data of this planet (Tsiaras et
al., 2016b; MacDonald & Madhusudhan, 2017). We assume a plane parallel
atmosphere with 100 layers and include the contributions of collision-induced
absorption (CIA) of H2-H2 (Abel et al., 2011; Fletcher et al., 2018) and H2-He
(Abel et al., 2012), Rayleigh scattering and grey-clouds. In terms of
molecular line lists, we import the following: H2O (Polyansky et al., 2018),
NH3 (Yurchenko et al., 2011), CH4 (Yurchenko et al., 2017) and HCN (Barber et
al., 2014).

Figure 3 displays the errors on the transit depth predicted by the radiometric
portion of Terminus as well as the uncertainties recovered from the light
curve fits. While the agreement is generally good, within 10%, there appears
to be a wavelength-dependent effect on the accuracy of the radiometric tool.
The trend seen could be due to the limb darkening coefficients, which change
with wavelength and alter the shape of the light curve.

Figure 3: Comparison of error bars obtained from the radiometric model (black)
and light curve fitting for HD 209458 b. The wavelength dependent difference
between the models could be due to limb darkening coefficients.

##  3 Orbit Modelling

Observatories in low Earth orbits can experience interruptions in target
visibility due to Earth occultations. Additionally, instruments and spacecraft
usually have specific target-Sun, target-Moon or target-Earth limb
restrictions. To account for these, Terminus is capable of modelling the orbit
of a spacecraft and calculating angles between the target and the Earth limb,
the Sun or other celestial body, in a similar way to tools used for other
missions (e.g. for CHEOPS: Kuntzer et al. (2014)).

The tool operates within an Earth-centred frame and the positions of celestial
objects (the Sun, Moon etc.) are loaded from the JPL Horizons
service333https://ssd.jpl.nasa.gov/horizons.cgi. The spacecraft’s orbit is
defined by an ellipse which is subsequently inclined with respect to the X
plane. The right ascension of the ascending node (RAAN) is then used to rotate
this about the Z axis.

Twinkle will operate in a Sun-synchronous orbit and here we modelled the
following orbital parameters: altitude = 700 km, inclination = 90.4∘,
eccentricity = 0, RAAN = 190.4∘ (i.e. 6am). These are subject to change based
upon launch availability but provide an approximate description of the
expected operational state. The orbit of Twinkle during May 2024 is depicted
in Figure 4.

Figure 4: Modelled orbit of Twinkle (red) during May 2024. The yellow vector
indicates the direction of the Sun while the black represents the anti-sun
vector (i.e. the centre of Twinkle’s field of regard). The Earth is
represented by the sphere with the terminator between day and night roughly
shown.

Figure 5: Sky coverage of Twinkle given the specific exclusion angles. The
effects of individual constrains are shown for the Sun, Earth and Moon
alongside the combination of them all. Stars indicate known transiting
exoplanet hosts with HD 209458 and WASP-127 highlighted by light blue and
green stars respectively. We note that the colour bar axes differ between each
plot.

Figure 6: Sky coverage of JWST (left) and Ariel (right) which will have
continuous viewing zones at the ecliptic poles. These missions are unaffected
by Earth obscuration due to their L2 orbit.

As mentioned, the code can impose a number of exclusion angles to explore
their effects on target availability. Here we modelled Sun, Earth and Moon
exclusion angles of 140 ∘, 20 ∘ and 5 ∘ respectively. The first of these is
largely due to thermal constraints while the latter two are to reduce stray
light. The Earth and Moon exclusion angles for Twinkle are still under study
but the values chosen here are similar to those of other observatories
operating in sun-synchronous orbits or those proposed to do so (Kuntzer et
al., 2014; Deroo et al., 2012).

The effects of each exclusion angle on the sky coverage is shown in Figure 5
along with the effect of combining them all. In each case, the metric shown is
the total time the area of sky can be observed over the course of a year. The
plots highlight Twinkle’s excellent coverage of the ecliptic plane although
it, like CHEOPS, lacks the ability to study planets close to the ecliptic
poles. However, the JWST and Ariel missions will prefer the polar regions, as
shown in Figure 6, and thus both Twinkle and CHEOPS provide complimentary
coverage.

##  4 Partial Light Curves

Figure 7: Comparison of the predicted gap sizes for HD 209458 b (RA = 330, Dec
= 18) from Terminus and Freeflyer. The transit light curves are offset for
clarity and the gap sizes are seen to be highly similar. We note that these
gaps are due solely to physical obscuration by the Earth and no exclusion
angle is included. Figure 8: Effect of different Earth exclusion angles on the
percentage of time on target (black) and size of the gaps (red) for a transit
observation of HD 209458 b. Figure 9: The 17 transits of HD 209458 b that are
observable with Twinkle over the course of a single year. The gaps are due to
Earth obscuration plus an exclusion angle of 20∘. All light curves have gaps
of roughly 45 minutes which are comparable to those in the Hubble data of the
same planet and have been offset for clarity.

Figure 10: Recovered spectrum and error bars from different light curve fits
for HD 209458 b. In each case, red represents the fitting of a full light
curve (same as Figure 3), blue the fitting of the partial light curve (LC1
from Figure 9) and black represents the predicted error from the radiometric
model. The partial light curve results in far larger uncertainties due to the
reduction in the number of data points.

From an exoplanet modelling perspective, it has thus far it has been assumed
that a full light curve is observed. However, in reality, for space-telescopes
in a low-Earth orbit, sometimes only partial light curves will be obtained due
to Earth obscuration as discussed in Section 3. These gaps cannot be
completely accounted for in radiometric models and thus a time-domain code,
such as Terminus, is required.

To verify the orbital code created, and to explore the effect of partial light
curves, we check our results against those of Edwards (2019). In Edwards
(2019), the mission design, analysis and operation software
Freeflyer444https://ai-solutions.com/freeflyer/ was used to model the
obscurations of HD 209458 b by the Earth throughout a year. FreeFlyer has
previously been used to support planning for several missions including NASA’s
Solar Dynamics Observatory (SDO). We note that Freeflyer only models the
physical obscuration of the target star by Earth and thus for this comparison
we set the Earth exclusion angle to zero.

As mentioned, Twinkle’s field of regard means targets are not constantly
observable and in a year 17 transits of HD 209458 b would be observable by
Twinkle. Given the sky location of HD 209458, Right Ascension (RA): 330.79∘;
Declination (Dec): 18.88∘, the target will always be periodically obscured by
the Earth. In Figure 7, we show a comparison between the predicted gaps for
the first of these transits which are shown to be in excellent agreement.
Meanwhile, Figure 8 displays the increase in gap size that would be incurred
by various Earth exclusion angles. Going from an angle of 0 to 20 degrees
increases the gaps size from 20 minutes to 44 minutes. The latter case would
mean Twinkle could be on-target for over half an orbit (54 minutes). In
comparison, past Hubble observations featured gaps of 47 minutes, with 48
minutes on target per orbit (Deming et al., 2013; Tsiaras et al., 2016b).
Hence, Twinkle’s observing efficiency for HD 209458 b will probably be similar
to that of Hubble. All potential transit observations of HD 209458 b have gaps
or a similar size (see Figure 9).

Here we fit the first available light curve and the recovered spectrum, and
associated errors, is shown in Figure 10. As expected, the gaps increase the
uncertainties on the recovered transit depth. Using Equations 5 and 8, one
would expect the error to increase by 35%
($\sigma_{p}=\sigma_{f}\times\frac{1}{\sqrt{0.55}}=1.347\sigma_{f}$). We see
an increase of 20-40% and thus the radiometric model may also provide
reasonable errors for partial light curves.

Figure 11: The 18 transits of WASP-127 b that are observable with Twinkle in
2024 which have been offset for clarity. The gaps are due to Earth obscuration
plus an exclusion angle of 20∘. LC 9 has no gaps, highlighting the importance
of observational planning with Twinkle, or other LEO satellites, and the
benefit of a sun-synchronous orbit over the equatorial orbit of Hubble.

Figure 12: Recovered spectrum and error bars from different light curve fits
for WASP-127 b. In each case, red represents the fitting of a full light curve
(e.g. LC9 in Figure 11), blue the fitting of the partial light curve (LC1 from
Figure 11) and black represents the predicted error from the radiometric
model. The errors from the full light curve are found to agree with the
radiometric prediction, again with the exception of a slight, wavelength
dependent, variation. The partial light curve results in far larger
uncertainties.

However, some planets may have more variable gaps, due to their location in
the sky and a changing spacecraft-Earth-target geometry, and thus may be
affected more significantly. For these planets, the scheduling of observations
is likely to be highly important. Terminus is able to provide input into
studies exploring the effects of partial light curves.

As an initial step to understand the variability of Earth obscuration, we now
model observations of WASP-127 b (Lam et al., 2017). WASP-127 is located such
that Twinkle will potentially have a continuous, unobstructed view of the
target during a transit (RA: 160.56∘, Dec: -3.84∘). However, some potential
observations will incur Earth obscuration and the amount of time lost will be
dependent upon the Earth exclusion angle required. In the case of the 20∘
exclusion angle modelled here, Twinkle would have access to one complete
transit (i.e. no gaps due to Earth obscuration) in 2024 as shown in Figure 11.
The other available observation periods would incur interruption up to a
maximum of 45 minutes over a 98 minute orbit. In the case of the Hubble
observations of WASP-127 b (Skaf et al., 2020; Spake et al., 2020), the
spacecraft could only be pointed at the target for 40 minutes per orbit (55
minute gaps). Hence, through careful selection of observing windows, the
efficiency of Twinkle’s observations of WASP-127 b could be far greater than
that of Hubble’s for this target.

To understand the impact of these gaps, we simulate a set of light curves for
a single observation of WASP-127 b and compare the errors on the transit
depths when gaps are induced. Again we base the atmosphere off of current
observations which suggest a large water abundance and potentially the
presence of FeH (Skaf et al., 2020), which we model using the line lists from
Dulick et al. (2003); Wende et al. (2010).

The results of these fittings are shown in Figure 12. The full light curve
again has a wavelength dependent variation from the predicted radiometric
errors but this is again relatively small. As expected, the fitting of the
partial light curve results in larger uncertainties on the transit depth. In
the case modelled, LC1 from Figure 11, Twinkle only observes the target for
46% of the transit. Using Equations 5 and 8, one would expect the error to
increase by 48%
($\sigma_{p}=\sigma_{f}\times\frac{1}{\sqrt{0.46}}=1.476\sigma_{f}$). We see
the increase is wavelength dependent and generally between 20-40%, less than
predicted. Therefore the radiometric model may not always be capable of
providing accurate error estimations.

The recovered precision on different parameters is likely to be dependent upon
the location of the gaps in the light curve. In this case the central portion
of the transit is well sampled allowing for a precise recovery of the transit
depth. However, ingress/egress are less well sampled and thus orbital
parameters such as the inclination (i) and reduced semi-major axis (a/R∗) may
be less well determined.

Furthermore, the standard methodology of analysing transiting exoplanet data
is to fit to the light curves for planet-to-star radius ratio (Rp/Rs) to
achieve a spectrum with error bars before performing atmospheric retrievals on
said spectrum. This approach, which has essentially been followed here,
distils time-domain observations down to a single point and thus much
information about the orbital parameters of the system are lost. Fitting of
full light curves (no gaps) usually retrieves the orbital parameters
accurately but, as discussed, gaps can lead to less certainty. This potential
degeneracy is lost in the standard method and so, to bring the data analysis
one step closer to the raw data, retrievals with Terminus generated data could
be conducted using the light curves themselves and the methodology described
in Yip et al. (2019). The so called “L-retrieval” allows for the orbital
parameters (e.g. inclination, semi-major axis) to be free parameters in the
retrieval to ensure that orbital degeneracies are accounted for. Such a
methodology would be useful in the exploration of the effects of Earth
obscuration, particularly as these orbital elements have been shown to be
important in recovering the correct optical slope (Alexoudi et al., 2018). A
thorough analysis is needed to explore this fully and Terminus can feed vital
information into such an effort.

##  5 Availability of Solar System Bodies

Twinkle will also conduct spectroscopy of objects within our Solar System with
perhaps the most promising use of the mission in this regard being the
characterisation of small bodies. In particular, a diverse array of shapes for
the 3 $\mu$m hydration feature, which generally cannot be observed from the
ground, have been seen and used to classify asteroids (e.g. Mastrapa et al.,
2009; Campins et al., 2010; Rivkin & Emery, 2010; Takir & Emery, 2012; Takir
et al., 2013). Twinkle’s broad wavelength coverage will allow for studies of
this spectral feature, and many others, as outlined in Edwards et al. (2019a).

The times at which major and minor Solar System bodies are within Twinkle’s
field of regard has previously been studied in Edwards et al. (2019a, c).
These studies showed that the outer planets, and main belt asteroids, will
have long, regular periods within Twinkle’s field of regard. However, the
observation periods of Near-Earth Objects (NEOs) and Near-Earth Asteroids
(NEAs) are far more sporadic. Hence, we revisit this analysis with the
addition of considering Earth obscuration. For our example target, we choose
99942 Apophis (2004 MN4), a potentially hazardous asteroid (PHA). Apophis has
a diameter of around 400 m (Licandro et al., 2015; Müller et al., 2014) and
will have a close fly-by in 2029 (Figure 13). While it had been thought there
was potentially a high probability of impact during this fly-by, or one in
2036, this has now been significantly downgraded (Krolikowska & Sitarski,
2010; Chesley et al., 2010; Thuillot et al., 2015). Nevertheless, passing
around 31,000 km from the Earth’s surface, Apophis will come within the orbits
of geosynchronous satellites (see Figure 13).

By comparing the data to likely meteorite analogues, current spectral analyses
of Apophis have concluded it is an Sq-class asteroid that closely resembles
the LL ordinary chondrite meteorites in terms of composition (Binzel et al.,
2009; Reddy et al., 2018). This data was measured over 0.55-2.45 $\mu m$ and
similarities have been noted to that of the asteroid Itokawa which was visited
and studied by the Hayabusa mission (Abe et al., 2006).

Figure 13: Top: orbit of Earth and Apophis from June 2028 to June 2029. In the
period, Apophis crosses the orbit of Earth twice with the second of these
crosses occurring during April 2029. Bottom: the distance between Earth and
Apophis during the April 2029, highlighting that the minimum separation from
the Earth surface is closer than geosynchronous satellites. Data for these
plots was acquired via the NASA JPL Horizons service.

Figure 14: Visible magnitude (top) and rate of apparent motion (bottom) for
Apophis during it’s close fly-by in 2029. The availability of Apophis was
checked at a cadence of 1 minute with dark blue indicating it is unobstructed,
light blue showing times at which the Earth is occulting the target and black
representing times when it has left the field of regard (i.e. exclusion due to
Sun-target angle). The left-hand plots show these values for the week before
the closest approach while the right-hand plots display the Earth obscuration
more readily as Apophis approaches a rate of 30 mas/s. Figure 15: Average sky
coverage during the two weeks before the closest approach of Apophis and the
sky location of Apophis over that same period (white). It should be noted
that, for the plotted Apophis trajectory, the time spent outside the FOR is
only a few hours whereas the time spent within it equates to several days.
Figure 16: Simulated spectra for Apophis. The error bars are for a single
exposure with a 300 s integration time on an object at a visible magnitude of
12. The spectrum is of an LL6 ordinary chondrite meteorite, taken from the
RELAB database (bkr1dp015). We note that the reflectance shown here at shorter
wavelengths ($<0.8\mu m$) is slightly larger than found in actual studies of
Apophis (Binzel et al., 2009; Reddy et al., 2018).

Here, we analyse the availability of Apophis over the week before, and day
after, its closest approach to Earth. Terminus obtains asteroid ephemerides
using the astropy API to the JPL Horizons database (Astropy Collaboration et
al., 2018). In Figure 14 we show the visible magnitude and apparent rate of
motion of Apophis during this period. The interlaced dark and light blue
segments show the availability of the asteroid before it leaves the field of
regard soon after the closest point of its fly-by. The trajectory across the
sky of Apophis is depicted in Figure 15 along with the sky coverage of Twinkle
over this period.

The ability of spacecraft to accurately track non-sidereal objects is key for
their observation. The Spitzer Space Telescope was used extensively for
characterising small bodies (e.g. Trilling et al., 2007; Barucci et al., 2008)
and tracked objects moving at rates of 543 mas/s (Trilling et al., 2010).
Spitzer was oriented using a inertial reference unit comprising of several
high performance star trackers and observed asteroids using linear track
segments. These were commanded as a vector rate in J2000 coordinates, passing
through a specified RA and Dec at a specified time. The coordinates of the
target can be obtained from services such as Jet Propulsion Laboratory’s
Horizons System. JWST is expected be able to track objects moving at up to 30
mas/s (Thomas et al., 2016).

The maximum rate at which Twinkle can track non-sidereal objects is still
under definition but will be >30 mas/s which we take here as a conservative
maximum value. When this threshold is crossed, Apophis will have a visible
magnitude of approximately 11.8. During the day or so before this rate limit
is crossed, Apophis would be available for periods of 55 minutes, with 40
minute interruptions, again assuming a 20∘ Earth exclusion angle.

As demonstrated in Figure 16, such observation windows provide plenty of time
to achieve high quality spectra. Here we simulated spectra for Apophis at a
visible magnitude of 12 and an integration time of 5 minutes. We note that the
thermal emission from the asteroid has been subtracted, which was modelled as
a blackbody with a temperature of 300 K, to give the relative reflectance of
the asteroid. The input spectrum was taken from the RELAB
database555http://www.planetary.brown.edu/relab/ and is of an LL6 ordinary
chondrite meteorite.

Simulations have suggested the 2029 close encounter could cause landslides on
Apophis, if the structure of some parts of the structure are significantly
weak (Yu et al., 2014). The potential for resurfacing NEOs during terrestrial
encounters in discussed in e.g. Binzel et al. (2010); Nesvorný et al. (2010)
and spectral measurements can inform us on the freshness of the asteroid’s
surface, providing evidence for such mechanisms. Additionally, while an impact
in 2029 has been ruled out, the potential for a future collision cannot be
disregarded and further study of the object is needed to refine this. In
particularly, the Yarkovsky effect has been shown to significantly alter
predictions beyond 2029 and is sensitive to the physical parameters of
Apophis, such as its albedo, diameter and density (Farnocchia et al., 2013; Yu
et al., 2017).

By observing Apophis simultaneously from 0.5-4.5 $\mu m$, Twinkle could
significantly inform the debate surrounding the nature of Apophis and it’s
potential threat level to Earth. Therefore, Twinkle could have a role to play
in characterising known NEOs and NEAs, along with those predicted to be
discovered by Near-Earth Object Surveillance Mission (NEOSM, Mainzer et al.
(2019)) and Vera C. Rubin Observatory, previously known as the Large Synoptic
Survey Telescope (LSST, Jones et al. (2018)). The ability of Twinkle to
contribute to the study of NEOs and NEAs, and other specific asteroid
populations, will be thoroughly detailed in further work.

##  6 Conclusions and Future Work

Terminus, a simulator with some time-domain capabilities has been developed to
model observations with space-based telescopes. This model is especially
applicable to exoplanets and can incorporate gaps in the light curve, caused
by Earth obscuration, and be used to predict the potential impact on the
accuracy of the retrieved atmospheric composition. Here, Terminus is baselined
on the Twinkle Space Telescope but the model can be adapted for any space-
based telescope and is especially applicable to those in a low-Earth orbit.

The impact of gaps in exoplanet observations has not been fully explored and
further work is needed. Obtaining a full transit, or eclipse, light curve is
obviously the ideal case but when it is not possible, such as for HD 209458 b,
an optimisation of the location, and length, of the gaps is required. By being
able model when these gaps occur, it should be possible to begin to explore
this by running multiple fittings and comparing the retrieved transit depth
and atmospheric parameters.

The Earth exclusion angle considered here is identical for the lit and unlit
portions of the Earth. However, each will contribute different amounts of
stray light and thus likely have separate exclusion angles. Future work will
incorporate this capability, along with the capacity to quantitatively model
the expected stray light from the Earth and Moon to firmly establish the
exclusion angles required. The effect of different orbital parameters (e.g.
altitude, 6am vs 6pm RAAN) can also be explored. Terminus will be updated to
include the South Atlantic Anomaly (SAA) to model the impact in the event that
the spacecraft must limit scientific operations during its ingress into this
region. Other additional development aspects include satellite ground stations
and calculating potential accesses to these facilities. Such capabilities will
allow for the tool to serve wider concept of operations (CONOPS) concerns and,
in the event that spacecraft design for any reason limits operations during
downlink, this can then be accounted for in the scheduling. Additionally,
Terminus could also be used to model other effects such as stellar variability
or detector ramps such as those seen on Hubble and Spitzer.

Finally, Terminus will be incorporated into a web interface to provide the
community with simulations of Twinkle’s capabilities. Doing so will allow the
tool to be more widely used and facilitate in-depth studies of Twinkle’s
capabilities. These could include modelling various atmospheric scenarios for
each planet to judge its suitability for characterisation (e.g. Fortenbach &
Dressing, 2020), performing retrievals on populations of exoplanets (e.g.
Changeat et al., 2020), classifying groups of planets via colour-magnitude
diagrams (e.g. Dransfield & Triaud, 2020), testing machine-learning techniques
for atmospheric retrieval (e.g. Márquez-Neila et al., 2018; Zingales &
Waldmann, 2018; Hayes et al., 2020; Yip et al., 2020) or the exploration of
potential biases in current data analysis techniques (e.g. Feng et al., 2016;
Rocchetto et al., 2016; Changeat et al., 2019; Caldas et al., 2019; Powell et
al., 2019; MacDonald et al., 2020; Taylor et al., 2020). Additionally,
thorough analyses of Twinkle’s capabilities for specific scientific
endeavours, such as confirming/refuting the presence of thermal inversions and
identifying optical absorbers in ultra-hot Jupiters (e.g. Fortney et al.,
2008; Spiegel et al., 2009; Haynes et al., 2015; Evans et al., 2018;
Parmentier et al., 2018; von Essen et al., 2020; Edwards et al., 2020; Pluriel
et al., 2020; Changeat & Edwards, 2021), searching for an exoplanet mass-
metallicity trend (e.g. Wakeford et al., 2017; Welbanks et al., 2019), probing
the atmospheres of planets in/close to the radius valley to discern their true
nature (e.g. Owen & Wu, 2017; Fulton & Petigura, 2018; Zeng et al., 2019),
refining basic planetary and orbital characteristics (e.g. Berardo et al.,
2019; Dalba & Tamburo, 2019; Livingston et al., 2019), measuring planet masses
through accurate transit timings (e.g. Hadden & Lithwick, 2017; Grimm et al.,
2018; Petigura et al., 2018), verifying additional planets within systems
(e.g. Gillon et al., 2017; Bonfanti et al., 2021), studying non-transiting
planets by measuring the planetary infrared excess (Stevenson, 2020), or even
contributing to the search for exomoon candidates (e.g. Simon et al., 2015;
Heller et al., 2016; Teachey & Kipping, 2018), can also be undertaken.

##  7 Acknowledgements

This work has utilised data from FreeFlyer, a mission design, analysis and
operation software created by a.i. solutions. We thank Giovanna Tinetti,
Marcell Tessenyi, Giorgio Savini, Subhajit Sarkar, Enzo Pascale, Angelos
Tsiaras, Philip Windred, Andy Rivkin, Lorenzo Mugnai, Kai Hou Yip, Ahmed Al-
Refaie, Quentin Changeat and Lara Ainsman for their guidance, comments and
useful discussions. This work has been partially funded by the STFC grant
ST/T001836/1.

Software: TauREx3 (Al-Refaie et al., 2019), pylightcurve (Tsiaras et al.,
2016a), ExoTETHyS (Morello et al., 2020), ExoSim (Sarkar et al., 2020b),
Astropy (Astropy Collaboration et al., 2018), h5py (Collette, 2013), emcee
(Foreman-Mackey et al., 2013), Matplotlib (Hunter, 2007), Multinest (Feroz et
al., 2009; Buchner et al., 2014), Pandas (McKinney, 2011), Numpy (Oliphant,
2006), SciPy (Virtanen et al., 2020), corner (Foreman-Mackey, 2016).

## References

  * Abe et al. (2006) Abe, M., Takagi, Y., Kitazato, K., et al. 2006, Science, 312, 1334.  https://science.sciencemag.org/content/312/5778/1334
  * Abel et al. (2011) Abel, M., Frommhold, L., Li, X., & Hunt, K. L. 2011, The Journal of Physical Chemistry A, 115, 6805
  * Abel et al. (2012) —. 2012, The Journal of chemical physics, 136, 044319
  * Agúndez et al. (2012) Agúndez, M., Venot, O., Iro, N., et al. 2012, Astronomy & Astrophysics, 548, A73.  http://dx.doi.org/10.1051/0004-6361/201220365
  * Akeson et al. (2013) Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989
  * Al-Refaie et al. (2019) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2019, arXiv e-prints, arXiv:1912.07759
  * Alexoudi et al. (2018) Alexoudi, X., Mallonn, M., von Essen, C., et al. 2018, A&A, 620, A142
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, Astronomy & Astrophysics, 577, A42
  * Barber et al. (2014) Barber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828
  * Barclay et al. (2018) Barclay, T., Pepper, J., & Quintana, E. V. 2018, ArXiv e-prints 1804.05050, arXiv:1804.05050
  * Barucci et al. (2008) Barucci, M., Fornasier, S., Dotto, E., et al. 2008, A&A, 477, 665
  * Batalha et al. (2017) Batalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017, PASP, 129, 064501
  * Benz et al. (2020) Benz, W., Broeg, C., Fortier, A., et al. 2020, arXiv e-prints, arXiv:2009.11633
  * Berardo et al. (2019) Berardo, D., Crossfield, I. J. M., Werner, M., et al. 2019, AJ, 157, 185
  * Binzel et al. (2009) Binzel, R. P., Rivkin, A. S., Thomas, C. A., et al. 2009, Icarus, 200, 480 .  http://www.sciencedirect.com/science/article/pii/S001910350800434X
  * Binzel et al. (2010) Binzel, R. P., Morbidelli, A., Merouane, S., et al. 2010, Nature, 463, 331
  * Bonfanti et al. (2021) Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, arXiv e-prints, arXiv:2101.00663
  * Bordé et al. (2003) Bordé, P., Rouan, D., & Léger, A. 2003, A&A, 405, 1137
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Caldas et al. (2019) Caldas, A., Leconte, J., Selsis, F., et al. 2019, A&A, 623, A161
  * Campins et al. (2010) Campins, H., Hargrove, K., Pinilla-Alonso, N., et al. 2010, Nature, 464, 1320
  * Changeat et al. (2020) Changeat, Q., Al-Refaie, A., Mugnai, L. V., et al. 2020, AJ, 160, 80
  * Changeat & Edwards (2021) Changeat, Q., & Edwards, B. 2021, arXiv e-prints, arXiv:2101.00469
  * Changeat et al. (2019) Changeat, Q., Edwards, B., Waldmann, I. P., & Tinetti, G. 2019, ApJ, 886, 39
  * Charbonneau et al. (2000) Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45
  * Chesley et al. (2010) Chesley, S. R., Baer, J., & Monet, D. G. 2010, Icarus, 210, 158
  * Claret et al. (2013) Claret, A., Hauschildt, P. H., & Witte, S. 2013, A&A, 552, A16
  * Collette (2013) Collette, A. 2013, Python and HDF5 (O’Reilly)
  * Dalba & Tamburo (2019) Dalba, P. A., & Tamburo, P. 2019, ApJ, 873, L17
  * Deming et al. (2013) Deming, D., Wilkins, A., McCullough, P., et al. 2013, ApJ, 774, 95
  * Deroo et al. (2012) Deroo, P., Swain, M. R., & Green, R. O. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Proc. SPIE, 844241
  * Dransfield & Triaud (2020) Dransfield, G., & Triaud, A. H. M. J. 2020, MNRAS, 499, 505
  * Dulick et al. (2003) Dulick, M., Bauschlicher, C. W., J., Burrows, A., et al. 2003, ApJ, 594, 651
  * Edwards (2019) Edwards, B. 2019, PhD thesis, University College London
  * Edwards et al. (2019a) Edwards, B., Lindsay, S., Savini, G., et al. 2019a, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 034004
  * Edwards et al. (2019b) Edwards, B., Mugnai, L., Tinetti, G., Pascale, E., & Sarkar, S. 2019b, Astronomical Journal, 157, 242
  * Edwards et al. (2019c) Edwards, B., Savini, G., Tinetti, G., et al. 2019c, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 014006
  * Edwards et al. (2019d) Edwards, B., Rice, M., Zingales, T., et al. 2019d, Experimental Astronomy, 47, 29
  * Edwards et al. (2020) Edwards, B., Changeat, Q., Baeyens, R., et al. 2020, AJ, 160, 8
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, arXiv e-prints, arXiv:2003.05528
  * Evans et al. (2018) Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ, 156, 283
  * Farnocchia et al. (2013) Farnocchia, D., Chesley, S. R., Chodas, P. W., et al. 2013, Icarus, 224, 192
  * Feng et al. (2016) Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, ApJ, 829, 52
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Fletcher et al. (2018) Fletcher, L. N., Gustafsson, M., & Orton, G. S. 2018, The Astrophysical Journal Supplement Series, 235, 24
  * Focardi et al. (2018) Focardi, M., Pace, E., Farina, M., et al. 2018, Experimental Astronomy, 46, 1
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24.  https://doi.org/10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fortenbach & Dressing (2020) Fortenbach, C. D., & Dressing, C. D. 2020, PASP, 132, 054501
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fulton & Petigura (2018) Fulton, B. J., & Petigura, E. A. 2018, ArXiv e-prints, arXiv:1805.01453
  * Futyan et al. (2020) Futyan, D., Fortier, A., Beck, M., et al. 2020, A&A, 635, A23
  * Garhart et al. (2020) Garhart, E., Deming, D., Mandell, A., et al. 2020, AJ, 159, 137
  * Geers et al. (2019) Geers, V. C., Klaassen, P. D., Beard, S., & European Consortium, M. 2019, in Astronomical Society of the Pacific Conference Series, Vol. 523, Astronomical Data Analysis Software and Systems XXVII, ed. P. J. Teuben, M. W. Pound, B. A. Thomas, & E. M. Warner, 641
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456
  * Gordon et al. (2016) Gordon, I., Rothman, L. S., Wilzewski, J. S., et al. 2016, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 48, AAS/Division for Planetary Sciences Meeting Abstracts #48, 421.13
  * Greene et al. (2016) Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17
  * Grimm et al. (2018) Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68
  * Hadden & Lithwick (2017) Hadden, S., & Lithwick, Y. 2017, AJ, 154, 5
  * Hayes et al. (2020) Hayes, J. J. C., Kerins, E., Awiphan, S., et al. 2020, MNRAS, 494, 4492
  * Haynes et al. (2015) Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, The Astrophysical Journal, 806, 146.  http://stacks.iop.org/0004-637X/806/i=2/a=146
  * Heller et al. (2016) Heller, R., Hippke, M., Placek, B., Angerhausen, D., & Agol, E. 2016, A&A, 591, A67
  * Henry et al. (2000) Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000, ApJ, 529, L41
  * Hoeijmakers et al. (2018) Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018, Nature, 560, 453
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, Astronomy & Astrophysics, 553, A6.  http://dx.doi.org/10.1051/0004-6361/201219058
  * Iyer et al. (2016) Iyer, A. R., Swain, M. R., Zellem, R. T., et al. 2016, ApJ, 823, 109
  * James et al. (1997) James, J. F., Mukai, T., Watanabe, T., Ishiguro, M., & Nakamura, R. 1997, MNRAS, 288, 1022
  * John (1988) John, T. L. 1988, A&A, 193, 189
  * Jones et al. (2018) Jones, R. L., Slater, C. T., Moeyens, J., et al. 2018, Icarus, 303, 181
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69
  * Krolikowska & Sitarski (2010) Krolikowska, M., & Sitarski, G. 2010, arXiv e-prints, arXiv:1009.2639
  * Kuntzer et al. (2014) Kuntzer, T., Fortier, A., & Benz, W. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9149, Proc. SPIE, 91490W
  * Lam et al. (2017) Lam, K. W. F., Faedi, F., Brown, D. J. A., et al. 2017, A&A, 599, A3
  * Leinert et al. (1998) Leinert, C., Bowyer, S., Haikala, L. K., et al. 1998, A&AS, 127, 1
  * Licandro et al. (2015) Licandro, J., Müller, T., Alvarez, C., Alí-Lagoa, V., & Delbò, M. 2015, arXiv e-prints, arXiv:1510.06248
  * Livingston et al. (2019) Livingston, J. H., Crossfield, I. J. M., Werner, M. W., et al. 2019, AJ, 157, 102
  * MacDonald et al. (2020) MacDonald, R. J., Goyal, J. M., & Lewis, N. K. 2020, ApJ, 893, L43
  * MacDonald & Madhusudhan (2017) MacDonald, R. J., & Madhusudhan, N. 2017, MNRAS, 469, 1979
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30
  * Mainzer et al. (2019) Mainzer, A., Bauer, J., Cutri, R., et al. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–1049
  * Márquez-Neila et al. (2018) Márquez-Neila, P., Fisher, C., Sznitman, R., & Heng, K. 2018, Nature Astronomy, 2, 719
  * Martin-Lagarde et al. (2019) Martin-Lagarde, M., Lagage, P.-O., Gastaud, R., et al. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–1879
  * Mastrapa et al. (2009) Mastrapa, R., Sandford, S., Roush, T., Cruikshank, D., & Dalle Ore, C. 2009, ApJ, 701, 1347
  * McKinney (2011) McKinney, W. 2011, Python for High Performance and Scientific Computing, 14
  * Morello et al. (2020) Morello, G., Claret, A., Martin-Lagarde, M., et al. 2020, AJ, 159, 75
  * Mugnai et al. (2020) Mugnai, L. V., Pascale, E., Edwards, B., Papageorgiou, A., & Sarkar, S. 2020, Experimental Astronomy, 50, 303
  * Müller et al. (2014) Müller, T. G., Kiss, C., Scheirich, P., et al. 2014, A&A, 566, A22
  * Murakami et al. (2007) Murakami, H., Baba, H., Barthel, P., et al. 2007, PASJ, 59, S369
  * Nagler et al. (2019) Nagler, P. C., Edwards, B., Kilpatrick, B., et al. 2019, Journal of Astronomical Instrumentation, 8, 1950011
  * Nesvorný et al. (2010) Nesvorný, D., Bottke, W. F., Vokrouhlický, D., Chapman, C. R., & Rafkin, S. 2010, Icarus, 209, 510
  * Nielsen et al. (2016) Nielsen, L. D., Ferruit, P., Giardino, G., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99043O
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Owen & Wu (2017) Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110.  https://doi.org/10.1051/0004-6361/201833059
  * Pascale et al. (2015) Pascale, E., Waldmann, I. P., MacTavish, C. J., et al. 2015, Experimental Astronomy, 40, 601
  * Petigura et al. (2018) Petigura, E. A., Benneke, B., Batygin, K., et al. 2018, AJ, 156, 89
  * Pluriel et al. (2020) Pluriel, W., Whiteford, N., Edwards, B., et al. 2020, arXiv e-prints, arXiv:2006.14199
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2597
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, 991016
  * Powell et al. (2019) Powell, D., Louden, T., Kreidberg, L., et al. 2019, ApJ, 887, 170
  * Puig et al. (2015) Puig, L., Isaak, K., Linder, M., et al. 2015, Experimental Astronomy, 40, 393
  * Puig et al. (2018) Puig, L., Pilbratt, G., Heske, A., et al. 2018, Experimental Astronomy, doi:10.1007/s10686-018-9604-3
  * Rauscher et al. (2007) Rauscher, B. J., Fox, O., Ferruit, P., et al. 2007, PASP, 119, 7786
  * Reddy et al. (2018) Reddy, V., Sanchez, J. A., Furfaro, R., et al. 2018, AJ, 155, 140
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Proceedings of the SPIE, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, 914320
  * Rivkin & Emery (2010) Rivkin, A., & Emery, J. 2010, Nature, 464, 1322
  * Rocchetto et al. (2016) Rocchetto, M., Waldmann, I. P., Venot, O., Lagage, P. O., & Tinetti, G. 2016, ApJ, 833, 120
  * Rothman & Gordon (2014) Rothman, L. S., & Gordon, I. E. 2014, in 13th International HITRAN Conference, June 2014, Cambridge, Massachusetts, USA
  * Sarkar et al. (2020a) Sarkar, S., Madhusudhan, N., & Papageorgiou, A. 2020a, MNRAS, 491, 378
  * Sarkar et al. (2016) Sarkar, S., Papageorgiou, A., & Pascale, E. 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, Vol. 9904, International Society for Optics and Photonics (SPIE), 1243 – 1251.  https://doi.org/10.1117/12.2234216
  * Sarkar et al. (2017) Sarkar, S., Papageorgiou, A., Tsikrikonis, I. A., Vandenbussche, B., & Pascale, E. 2017, ARIEL Performance Analysis Report: ARIEL-CRDF-PL-AN-001
  * Sarkar et al. (2020b) Sarkar, S., Pascale, E., Papageorgiou, A., Johnson, L. J., & Waldmann, I. 2020b, arXiv e-prints, arXiv:2002.03739
  * Savini et al. (2016) Savini, G., Tessenyi, M., Tinetti, G., et al. 2016, in Preprint (SPIE astronomical telescopes + instrumentation 2016 paper 9904—175)
  * Simon et al. (2015) Simon, A. E., Szabó, G. M., Kiss, L. L., Fortier, A., & Benz, W. 2015, PASP, 127, 1084
  * Sing et al. (2016) Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59
  * Skaf et al. (2020) Skaf, N., Fabienne Bieger, M., Edwards, B., et al. 2020, arXiv e-prints, arXiv:2005.09615
  * Spake et al. (2020) Spake, J. J., Sing, D. K., Wakeford, H. R., et al. 2020, arXiv e-prints, arXiv:1911.08859
  * Spiegel et al. (2009) Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487
  * Sreejith et al. (2019) Sreejith, A. G., Fossati, L., Fleming, B. T., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 018004
  * Stevenson (2020) Stevenson, K. B. 2020, ApJ, 898, L35
  * Takir & Emery (2012) Takir, D., & Emery, J. 2012, Icarus, 219, 641
  * Takir et al. (2013) Takir, D., Emery, J., McSween, H., et al. 2013, Meteoritics and Planetary Science, 48, 1618
  * Taylor et al. (2020) Taylor, J., Parmentier, V., Irwin, P. G. J., et al. 2020, MNRAS, 493, 4342
  * Teachey & Kipping (2018) Teachey, A., & Kipping, D. M. 2018, Science Advances, 4, eaav1784
  * Thomas et al. (2016) Thomas, C., Abell, P., Castillo-Rogez, J., et al. 2016, PASP, 128, 018002
  * Thuillot et al. (2015) Thuillot, W., Bancelin, D., Ivantsov, A., et al. 2015, A&A, 583, A59
  * Tinetti et al. (2012) Tinetti, G., Beaulieu, J. P., Henning, T., et al. 2012, Experimental Astronomy, 34, 311
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, in press.  https://doi.org/10.1007/s10686-018-9598-x
  * Trilling et al. (2007) Trilling, D., Bhattacharya, B., Blaylock, M., et al. 2007, in Bulletin of the American Astronomical Society, Vol. 39, AAS/Division for Planetary Sciences Meeting Abstracts #39, 484
  * Trilling et al. (2010) Trilling, D., Mueller, M., Hora, J., et al. 2010, ArXiv e-prints, arXiv:1007.1009
  * Tsiaras et al. (2016a) Tsiaras, A., Waldmann, I., Rocchetto, M., et al. 2016a, ascl:1612.018
  * Tsiaras et al. (2016b) Tsiaras, A., Waldmann, I. P., Rocchetto, M., et al. 2016b, ApJ, 832, 202
  * Tsiaras et al. (2018) Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156
  * Tsumura et al. (2010) Tsumura, K., Battle, J., Bock, J., et al. 2010, ApJ, 719, 394
  * Tucker et al. (2018) Tucker, G. S., Nagler, P., Butler, N., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107025G
  * Varley et al. (2017) Varley, R., Tsiaras, A., & Karpouzas, K. 2017, ApJS, 231, 13
  * Venot et al. (2012) Venot, O., Hébrard, E., Agúndez, M., et al. 2012, Astronomy & Astrophysics, 546, A43.  http://dx.doi.org/10.1051/0004-6361/201219310
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * von Essen et al. (2020) von Essen, C., Mallonn, M., Hermansen, S., et al. 2020, arXiv e-prints, arXiv:2003.06424
  * Wakeford et al. (2017) Wakeford, H. R., Sing, D. K., Kataria, T., et al. 2017, Science, 356, 628
  * Waldmann et al. (2015a) Waldmann, I. P., Rocchetto, M., Tinetti, G., et al. 2015a, ApJ, 813, 13
  * Waldmann et al. (2015b) Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015b, ApJ, 802, 107
  * Welbanks et al. (2019) Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJ, 887, L20
  * Wende et al. (2010) Wende, S., Reiners, A., Seifahrt, A., & Bernath, P. F. 2010, A&A, 523, A58
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Yip et al. (2020) Yip, K. H., Changeat, Q., Nikolaou, N., et al. 2020, arXiv e-prints, arXiv:2011.11284
  * Yip et al. (2019) Yip, K. H., Waldmann, I. P., Tsiaras, A., & Tinetti, G. 2019, arXiv e-prints, arXiv:1811.04686
  * Yu et al. (2017) Yu, L.-L., Ji, J., & Ip, W.-H. 2017, Research in Astronomy and Astrophysics, 17, 070
  * Yu et al. (2014) Yu, Y., Richardson, D. C., Michel, P., Schwartz, S. R., & Ballouz, R.-L. 2014, Icarus, 242, 82
  * Yurchenko et al. (2017) Yurchenko, S. N., Amundsen, D. S., Tennyson, J., & Waldmann, I. P. 2017, A&A, 605, A95
  * Yurchenko et al. (2011) Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828
  * Yurchenko & Tennyson (2012) Yurchenko, S. N., & Tennyson, J. 2012, in EAS Publications Series, Vol. 58, EAS Publications Series, ed. C. Stehlé, C. Joblin, & L. d’Hendecourt, 243–248
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723
  * Zingales & Waldmann (2018) Zingales, T., & Waldmann, I. P. 2018, AJ, 156, 268


.tifpng.pngconvert #1 .tif

# Towards Explainable LiDAR Point Cloud Semantic Segmentation via Gradient
Based Target Localization

Abhishek Kuriyal1 and Vaibhav Kumar1 1Abhishek Kuriyal and Vaibhav Kumar are
with the GeoAI4Cities Lab, Department of Data Science and Engineering, IISER
Bhopal.
{abhishek19<EMAIL_ADDRESS>
###### Abstract

Semantic Segmentation (SS) of LiDAR point clouds is essential for many
applications, such as urban planning and autonomous driving. While much
progress has been made in interpreting SS predictions for images, interpreting
point cloud SS predictions remains a challenge. This paper introduces pGS-CAM,
a novel gradient-based method for generating saliency maps in neural network
activation layers. Inspired by Grad-CAM, which uses gradients to highlight
local importance, pGS-CAM is robust and effective on a variety of datasets
(SemanticKITTI, Paris-Lille3D, DALES) and 3D deep learning architectures
(KPConv, RandLANet). Our experiments show that pGS-CAM effectively accentuates
the feature learning in intermediate activations of SS architectures by
highlighting the contribution of each point. This allows us to better
understand how SS models make their predictions and identify potential areas
for improvement. Relevant codes are available at
https://github.com/geoai4cities/pGS-CAM.

##  I INTRODUCTION

Deep Learning (DL) techniques, such as Convolutional Neural Networks (CNNs)
employed for feature extraction, have gained substantial traction in various
applications, including instance classification, semantic segmentation, object
detection, and visual tracking [1, 2]. Nonetheless, a significant challenge
lies in the interpretability of their outcomes. This challenge can benefit
from the effective use of various Explainable AI (XAI) methods, with extensive
resources available in recent reviews by Van et al. [3] and Linardatos et al.
[4]. Over the years, Grad-CAM [5] has been one of the most widely used
techniques that uses the gradients of the output with respect to the final
convolutional layer to produce localization maps (saliency maps) highlighting
important regions. These important regions are most influential in the model’s
decision-making process. However, the implementation of Grad-CAM is
predominantly observed in studies focusing on tasks related to images.

Figure 1: pGS-CAM for single point (white dot) and class Road, points
corresponding to subset of class Person, points corresponding to all points of
class Car. All heatmaps are obtained with respect to the final activation
layers.

In the past few years, there has been a notable rise in interest towards
Semantic Segmentation (SS) applied to Light Detection and Ranging (LiDAR)
point clouds. This trend is attributed to LiDAR’s capacity for precise 3D
environmental perception. The reader can refer to the following review papers
[2, 6, 7] highlighting deep learning-based point cloud semantic segmentation
research in recent years.

Recently, efforts are made to explain feature extraction of point cloud DL
architectures for instance classification and object detection. Zhang et al.
[8] explains PointNet by visualizing what each point function has learnt
followed by generating class-attentive response maps. However, their approach
involved modifying the network’s architecture to utilize class-attentive
features. Zheng et al. [9] proposed a differentiable point-shifting process
that computes contribution scores to input points based on the loss value, but
it sacrifices global interpretability for highly accurate individual point
contributions. The resulting heatmaps resemble fine details similar to LRP
[10] methods used in images. Other methods such as [11, 12] utilize Grad-CAM
for instance classification and object detection of point clouds,
respectively. Needless to say, there is a lack of dedicated research on
interpretable SS of point clouds. Additionally, no existing study has proposed
and implemented a method utilizing gradients in architectures for
interpretation of large-scale outdoor point segmentation. Our study aims to
fill these gaps by making the following major contributions:

  * •

We propose pGS-CAM (point Grad-Seg Class Activation Mapping), a novel method
inspired by Grad-CAM. It creates saliency maps for each activation layer in
LiDAR Semantic Segmentation (LSS) architectures, revealing influential points.
Our approach marks the first to offer visual explanations for LSS process in
3D DL architectures.

  * •

Heatmaps generated via pGS-CAM can provide valuable purpose in comprehending
architecture’s working. One such utility is to verify that the network does
not fixate on peculiar intricacies unique to the training scans, which may not
extend to unseen scans.

  * •

We provide comprehensive qualitative analysis on two well-known point-based
architectures, KPConv [13] and RandLA-Net [14]. We utilize widely recognized
MLS datasets, including SemanticKITTI [15], Paris-Lille3D [13], and the DALES
ALS dataset [16]. We emphasize distinctions in the inner workings of
architectures, present counterfactual explanations, provide quantitative
analysis of pGS-CAM explanations etc., to prove method’s robustness.

##  II Methodology

Let ${\\{A^{k}\\}}^{K}_{k=1}$ denote the chosen feature maps of interest,
representing the K kernels from the final convolutional layer of the
classification network, and $l^{c}$ represents the logit for a selected class
$c$. Grad-CAM calculates the importance weight $\alpha^{k}_{c}$ for each
feature map $A^{k}$ by averaging the gradients of $l^{c}$ with respect to all
$N$ pixels (indexed as $u,v$) (see Equation 2). The heatmap

$L^{c}=ReLU\left(\sum_{k}\alpha^{k}_{c}A^{k}\right)$ (1)

with

$\alpha^{k}_{c}=\frac{1}{N}\sum_{u,v}\frac{\partial l^{c}}{\partial
A^{k}_{uv}}$ (2)

is formed by aggregating the feature maps using weights $\alpha^{k}_{c}$. This
summation is followed by the application of the ReLU function at the pixel
level to zero out negative values, thus emphasizing areas that positively
contribute to the classification decision for class $c$.

A fundamental challenge arises from the distinction that, while a image
classification network generates a single logit $l^{c}$ for each input image,
a LSS network provides logits $l^{c}_{i}$ for each of the $N$ points within
the point cloud $P$. For any intermediate activation layer $A$ in $R^{(M\times
k)}$ (consisting of M downsampled points with k feature dimensions), the
gradient of the logit $l_{i}$ for the target class $c$ with respect to $A_{k}$
is expressed by Equation 3.

$\frac{\partial l^{c}_{i}}{\partial A^{k}}$ (3)

To capture the collective gradient impact across all logits $l^{c}_{i}$, where
$i\in N$, we employ an aggregation operation. Our observations indicate that
summation is particularly effective for this purpose. Furthermore, there may
be cases where we require explanations for a limited subset of points, rather
than the entire point cloud. To address this need, we introduce pGS-CAM, which
focuses on a subset of points denoted as $P^{\prime}$, consisting of
$N^{\prime}$ points. We aim to provide individual explanations for each point
within this subset. $P^{\prime}$ may encompass a single point, points
associated with a specific class instance, or all $N$ points within the point
cloud (refer to Fig. 1 for a visual demonstration). As a result, the overall
pointwise gradient influence for this particular subset is defined by Equation
4.

$\frac{\partial l^{c}_{1}}{\partial A^{k}}+\frac{\partial l^{c}_{2}}{\partial
A^{k}}+...+\frac{\partial l^{c}_{N^{\prime}}}{\partial
A^{k}}=\frac{\partial\sum_{i=1}^{N^{\prime}}l_{i}^{c}}{\partial
A^{k}}=\frac{\partial L_{P^{\prime}}}{\partial A^{k}}$ (4)

Thereafter, the gradient influence $G^{k}$ similar to $\alpha^{k}_{c}$ in
Grad-CAM is approximated. $G^{k}$ (see Equation 5) is calculated by
aggregating over M points in activation layer $A_{k}$.

$G^{k}=\sum_{j=1}^{M}\left[\frac{\partial L_{P^{\prime}}}{\partial
A^{k}_{j}}\right]$ (5)

Heatmap at activation layer $A^{k}$ is generated via matrix multiplication of
$G^{k}$ with $A^{k}$ followed by ReLU operation to highlight positive
contributions only and MinMax normalization to align floating point values
between 0 and 1 (see Equation 6). The normalization is done to maintain
heatmap consistency over multiple point cloud instances aiding qualitative and
quantitative comparisons.

$H=MinMax(ReLU(\sum_{k}G^{k}A^{k}))$ (6)

The standard pGS-CAM produces saliency maps for each neural network activation
layer. However, in many semantic segmentation architectures with Encoder-
Decoder structures, the resolution of these maps decreases as we move through
the layers. To overcome this, we propose KDTree upsampling to maintain the
original point cloud resolution (see Algorithm 1). Note that unlike image
classification networks, where intermediate activation maps can be used as
coarse heatmaps, this isn’t the case for LSS networks. Localization is
apparent only after applying gradient influence $G_{k}$, highlighting a
significant difference between these architectural approaches.

1 Function  _Main(_Point cloud $P$, Segmentation Network $f$_)_:

Params: Point cloud $P$, Segmentation Network $f$

Data: Dataset $D$ with $N$ points

// $P^{M}$: Point resolution with resolution $M$

// $P^{N}$: Point resolution with resolution $N$

2 $l^{c}\leftarrow f(P)$

3 $G^{k}\leftarrow\sum_{j=1}^{M}\left[\frac{\partial L_{P^{\prime}}}{\partial
A^{k}_{j}}\right]$

4 $H_{A}\leftarrow\text{MinMax}(\text{ReLU}(\sum_{k}G^{k}A^{k}))$

5 tree $\leftarrow$ KDTree($P^{M}$, leaf_size=40)

6 upscale_idx $\leftarrow$ tree.query($P^{N}$)

7 $H_{A}\leftarrow H_{A}$[upscale_idx]

8 return _$H_{A}$_

9

Algorithm 1 pGS-CAM algorithm

##  III Experimentation and Discussion

The subsequent sections will illustrate the robustness of pGS-CAM. In Section
A, we offer a qualitative analysis of saliency maps and compare how these
architectures function at their core. In Section B, we connect pGS-CAM
heatmaps with dimensionality reduction techniques like PCA and t-SNE to
validate the heatmap results and offer diverse perspectives for heatmap
analysis. Section C presents intriguing counterfactual explanations,
highlighting points that significantly influence the network’s performance.
Lastly, in Section D, we provide a quantitative analysis of pGS-CAM
explanations through high-drop and low-drop experiments.

###  III-A Qualitative Analysis of Saliency Maps

Fig. 2 demonstrates pGS-CAM explanations for the car class in the
SemanticKITTI dataset using the RandLA-Net architecture. In the absence of
KDTree upsampling, the heatmaps exhibit varying resolutions due to differences
in activation layer resolutions. This is particularly challenging in
bottleneck layers with reduced point density. To address this issue, Fig. 2
presents full-resolution heatmaps using Algorithm 1 and nearest neighbor
upsampling, offering a close visual representation of decision-making at each
layer. A consistent pattern emerges, with edge-like structures in the initial
activation layer A1 and localization starting from A5 onwards. A2-A4 heatmaps
highlight road and terrain points, which influence car class predictions, and
A8 closely resembles per-point predictions for the car class. Localization
generally begins around A5 in RandLA-Net.

Fig. 3 displays pGS-CAM explanations for the KPConv architecture for the
Paris-Lille3D dataset. Similar to Fig. 2, initial activation layers A1-A2 show
edge-like structures. However, localization in KPConv starts from the A3
layer, indicating its superior representation capabilities. Fig. 4 presents
pGS-CAM heatmaps for the KPConv architecture for the DALES aerial dataset. The
pattern is similar to the Paris-Lille3D dataset, with early localization in
the encoder layers and potential redundancy in later layers (A6, A7, and A8).
Poor gradient highlights at encoder 3 suggest a high number of negative inputs
affecting information representation in that layer.

It is important to highlight that pGS-CAM heatmaps consistently exhibit
reproducibility across all datasets and architectures, making them valuable
for benchmarking purposes.

Figure 2: pGS-CAM heatmaps for SemanticKITTI car class and RandLA-Net
architecture. Figure 3: pGS-CAM heatmaps for Paris-Lille3D car class and
KPConv architecture. Figure 4: pGS-CAM heatmaps for DALES building class and
KPConv architecture.

###  III-B Interpreting pGS-CAM via Dimensionality Reduction

We conducted a comparative analysis between pGS-CAM heatmaps and 2D
dimensionality reduction techniques applied to corresponding activation
layers. Specifically, we employed Principal Component Analysis (PCA) and
t-Distributed Stochastic Neighbor Embedding (t-SNE). In the initial activation
layer, pGS-CAM revealed the presence of clustered points and overlapping
features among car, vegetation, and building points across all three methods.
However, in terms of interpretability and result differentiation, t-SNE
exhibited more intuitive and distinct patterns when compared to PCA. This
superiority is attributed to t-SNE’s advanced dimensionality reduction
capabilities (see Fig. 5). Transitioning to the final activation layer, we
observed a localized pGS-CAM heatmap, consistent with the PCA and t-SNE plots
(see Figure 6). In a specific case, we observed a slight highlighting of
vegetation points alongside car points in the heatmap, aligning with the t-SNE
plots, where red points (representing cars) were surrounded by green points
(representing vegetation). These findings highlight the coherence between pGS-
CAM and PCA, t-SNE observations and suggest that pGS-CAM offers more intuitive
per-point explanations compared to both techniques.

Figure 5: pGS-CAM heatmaps comparison with PCA and TSNE plots for initial
activation layer and target class $car$. Notice the highlights of vegetation
points in heatmap which is reflected by the cluttering of car and vegetation
in PCA and TSNE plots. Figure 6: pGS-CAM heatmaps comparison with PCA and TSNE
plots for final activation layer. Notice the presence of highlighted points
corresponding to the ”vegetation” class in close proximity to points belonging
to the ”car” class within the heatmap. This finding is further supported by
the t-SNE plot, which demonstrates a slight overlap between these two classes.

###  III-C Counterfactual pGS-CAM Explanations

An intriguing analysis, applicable to any explanation method, involves
identifying points that have a negative impact on the performance of a
specific class. This valuable insight can be obtained through a simple
adjustment in the gradient influence $G_{k}$. By negating the gradient
influence (as described in Equation 7), we can spotlight regions that
adversely affect the decision of the target class.

$G^{k}=\sum_{j=1}^{M}\left[-\frac{\partial L_{P^{\prime}}}{\partial
A^{k}_{j}}\right]$ (7)

Fig. 7 provides illustrative examples of these scenarios. The first example
presents a counterfactual explanation for the car class, suggesting that the
removal of fence points might enhance the performance of this class.
Similarly, the second example yields a similar result for the road class,
highlighting the influence of vegetation points. However, the most intriguing
case is the third example, which suggests that the performance of the person
class could be improved by eliminating all ground points. This observation
aligns with intuition, as LSS architectures often face challenges when dealing
with substantial differences in object sizes. Approaches such as KPConv employ
adaptive kernel convolutions to address these peculiarities.

Figure 7: Using negative gradient influence to highlight counterfactual pGS-
CAM explanations. These findings offer insightful inferences into the
contributing factors behind the adverse performance of the target class.

###  III-D Quantitative Point-Drop Experiments

Point Drop serves as a robustness assessment method for explanation
techniques. It involves removing points with high saliency scores (high drop),
causing a significant performance drop, or dropping points with low saliency
scores (low drop), resulting in a minor performance decrease. In our
experiments, we computed class-specific Intersection over Union (IoU)
stability and mean IoU stability after each point drop attack. Interestingly,
in certain scenarios, the model correctly identified objects as ”car” even
after the removal of 5000 points. Fig. 8 illustrates the IoU difference after
each drop attack, averaged across point clouds from the SemanticKITTI
validation set. High drop attacks led to a substantial class-specific IoU
decrease, averaging at 1.36, while low drop attacks maintained relatively
stable class-specific IoU, with an average decrease of 0.15. This underscores
the effectiveness of pGS-CAM. When calculating mean IoU for all classes, not
limited to the target class, the model’s performance remained reasonably
stable. The average decrease in mean IoU was 0.43 for high drop attacks and
0.34 for low drop attacks. Consequently, the overall performance decline
across all classes was less severe compared to focusing solely on the target
class.

Figure 8: Quantitative Point-Drop removal for car class and affect on mean IoU
performance for all class. Notice a significant decline in the performance of
the car class following a high-drop scenario.

##  IV Conclusion

We introduced pGS-CAM, a gradient based approach for explaining point cloud
semantic segmentation models. This method reveals discernible patterns across
network layers and consistently performs on popular datasets. When comparing
architectures, we found differences in localization initiation, with KPConv
displaying early localization compared to RandLA-Net. pGS-CAM effectively
aligns with 2D dimensionality reduction techniques, highlighting class
overlaps. We also proposed a gradient influence modification for
counterfactual explanations that improved class-specific performance. Our
Point-Drop experiment established pGS-CAM’s robustness, making it a valuable
benchmark. In summary, pGS-CAM provides fresh insights into the decision
making of point cloud semantic segmentation models.

## References

  * [1] J. Chai, H. Zeng, A. Li, and E. W. Ngai, “Deep learning in computer vision: A critical review of emerging techniques and application scenarios,” _Machine Learning with Applications_ , vol. 6, p. 100134, 2021.
  * [2] S. Dong, P. Wang, and K. Abbas, “A survey on deep learning and its applications,” _Computer Science Review_ , vol. 40, p. 100379, 2021.
  * [3] B. H. Van der Velden, H. J. Kuijf, K. G. Gilhuijs, and M. A. Viergever, “Explainable artificial intelligence (xai) in deep learning-based medical image analysis,” _Medical Image Analysis_ , vol. 79, p. 102470, 2022.
  * [4] P. Linardatos, V. Papastefanopoulos, and S. Kotsiantis, “Explainable ai: A review of machine learning interpretability methods,” _Entropy_ , vol. 23, no. 1, p. 18, 2020.
  * [5] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-cam: Visual explanations from deep networks via gradient-based localization,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 618–626.
  * [6] Y. Guo, H. Wang, Q. Hu, H. Liu, L. Liu, and M. Bennamoun, “Deep learning for 3d point clouds: A survey,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 43, no. 12, pp. 4338–4364, 2020.
  * [7] P. K. Vinodkumar, D. Karabulut, E. Avots, C. Ozcinar, and G. Anbarjafari, “A survey on deep learning based segmentation, detection and classification for 3d point clouds,” _Entropy_ , vol. 25, no. 4, p. 635, 2023.
  * [8] B. Zhang, S. Huang, W. Shen, and Z. Wei, “Explaining the pointnet: What has been learned inside the pointnet?” in _CVPR Workshops_ , 2019, pp. 71–74.
  * [9] T. Zheng, C. Chen, J. Yuan, B. Li, and K. Ren, “Pointcloud saliency maps,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 1598–1606.
  * [10] A. Binder, S. Bach, G. Montavon, K.-R. Müller, and W. Samek, “Layer-wise relevance propagation for deep neural network architectures,” in _Information science and applications (ICISA) 2016_. Springer, 2016, pp. 913–922.
  * [11] Y. Souai, G. Rouhafzay, and A.-M. Cretu, “A deep-learning-based approach for saliency determination on point clouds,” _Engineering Proceedings_ , vol. 27, no. 1, p. 17, 2022.
  * [12] D. Dworak and J. Baranowski, “Adaptation of grad-cam method to neural network architecture for lidar pointcloud object detection,” _Energies_ , vol. 15, no. 13, p. 4681, 2022.
  * [13] H. Thomas, C. R. Qi, J.-E. Deschaud, B. Marcotegui, F. Goulette, and L. J. Guibas, “Kpconv: Flexible and deformable convolution for point clouds,” in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 6411–6420.
  * [14] Q. Hu, B. Yang, L. Xie, S. Rosa, Y. Guo, Z. Wang, N. Trigoni, and A. Markham, “Randla-net: Efficient semantic segmentation of large-scale point clouds,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 11 108–11 117.
  * [15] J. Behley, M. Garbade, A. Milioto, J. Quenzel, S. Behnke, C. Stachniss, and J. Gall, “Semantickitti: A dataset for semantic scene understanding of lidar sequences,” in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 9297–9307.
  * [16] N. Varney, V. K. Asari, and Q. Graehling, “Dales: A large-scale aerial lidar data set for semantic segmentation,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops_ , 2020, pp. 186–187.