11institutetext:  ETH Zurich, Institute for Particle Physics & Astrophysics,
Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland 22institutetext: National
Center of Competence in Research PlanetS (www.nccr-planets.ch)
33institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D,
3001, Leuven, Belgium 44institutetext: STAR Institute, University of Liège,
19C allée du Six Août, 4000 Liège, Belgium 55institutetext: Centre Spatial de
Liège, Université de Liège, Avenue Pré-Aily, 4031 Angleur, Belgium
66institutetext: Institut de Ciències de l’Espai (ICE, CSIC), Campus UAB,
C/Can Magrans s/n, 08193 Bellaterra, Spain 77institutetext: Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
88institutetext: University of Oxford, Atmospheric, Oceanic and Planetary
Physics, Department of Physics, Sherrington Road, Oxford OX1 3PU, United
Kingdom 99institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17,
69117 Heidelberg, Germany 1010institutetext: Faculty of Aerospace Engineering,
Delft University of Technology, Kluyverweg 1, 2629 Delft, Netherlands
1111institutetext: www.life-space-mission.com

# Large Interferometer For Exoplanets (_LIFE_):

II. Signal simulation, signal extraction and fundamental exoplanet parameters
from single epoch observations

Felix Dannert Equal contribution, correspondence to:
<EMAIL_ADDRESS>Interferometer For Exoplanets (_LIFE_):Large
Interferometer For Exoplanets (_LIFE_):Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Maurice Ottiger Equal
contributionLarge Interferometer For Exoplanets (_LIFE_):Large Interferometer
For Exoplanets (_LIFE_): Sascha P. Quanz Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_):Large Interferometer For
Exoplanets (_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Romain
Laugier Large Interferometer For Exoplanets (_LIFE_):Large Interferometer For
Exoplanets (_LIFE_): Emile Fontanet Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Adrian Gheorghe Large
Interferometer For Exoplanets (_LIFE_):Large Interferometer For Exoplanets
(_LIFE_): Olivier Absil F.R.S.-FNRS Research AssociateLarge Interferometer For
Exoplanets (_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Colin
Dandumont Large Interferometer For Exoplanets (_LIFE_):Large Interferometer
For Exoplanets (_LIFE_): Denis Defrère Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Carlos Gascón Large
Interferometer For Exoplanets (_LIFE_):Large Interferometer For Exoplanets
(_LIFE_): Adrian M. Glauser Large Interferometer For Exoplanets (_LIFE_):Large
Interferometer For Exoplanets (_LIFE_): Jens Kammerer Large Interferometer For
Exoplanets (_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Tim
Lichtenberg Large Interferometer For Exoplanets (_LIFE_):Large Interferometer
For Exoplanets (_LIFE_): Hendrik Linz Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_): Jerôme Loicq Large
Interferometer For Exoplanets (_LIFE_):Large Interferometer For Exoplanets
(_LIFE_):Large Interferometer For Exoplanets (_LIFE_):Large Interferometer For
Exoplanets (_LIFE_): the _LIFE_ collaboration Large Interferometer For
Exoplanets (_LIFE_):Large Interferometer For Exoplanets (_LIFE_):

(Received: ¡date¿ / Accepted: ¡date¿)

###### Abstract

Context. The _Large Interferometer For Exoplanets (LIFE)_ initiative is
developing the science and a technology roadmap for an ambitious space mission
featuring a space-based mid-infrared (MIR) nulling interferometer in order to
detect the thermal emission of hundreds of exoplanets and characterize their
atmospheres.

Aims. In order to quantify the science potential of such a mission, in
particular in the context of technical trade-offs, an instrument simulator is
required. In addition, signal extraction algorithms are needed to verify that
exoplanet properties (e.g., angular separation, spectral flux) contained in
simulated exoplanet datasets can be accurately retrieved.

Methods. We present LIFEsim, a software tool developed for simulating
observations of exoplanetary systems with an MIR space-based nulling
interferometer. It includes astrophysical noise sources (i.e., stellar leakage
and thermal emission from local zodiacal and exo-zodiacal dust) and offers the
flexibility to include instrumental noise terms in the future. Here, we
provide some first quantitative limits on instrumental effects that would
allow the measurements to remain in the fundamental noise limited regime. We
demonstrate updated signal extraction approaches to validate signal-to-noise
ratio (SNR) estimates from the simulator. Monte-Carlo simulations are used to
generate a mock survey of nearby terrestrial exoplanets and determine to which
accuracy fundamental planet properties can be retrieved.

Results. LIFEsim provides an accessible way for predicting the expected SNR of
future observations as a function of various key instrument and target
parameters. The SNRs of the extracted spectra are photon-noise dominated, as
expected from our current simulations. Signals from multi-planet systems can
be reliably extracted. From single epoch observations in our mock survey of
small ($R<1.5R_{\mathrm{Earth}}$) planets orbiting within the habitable zones
of their stars, we find that typical uncertainties in the estimated effective
temperature of the exoplanets are $\lesssim$10%, for the exoplanet radius
$\lesssim$20%, and for the separation from the host star $\lesssim$2%. SNR
values obtained in the signal extraction process deviate less than 10% from
purely photon-counting statistics based SNRs.

Conclusions. LIFEsim has been sufficiently well validated so that it can be
shared with a broader community interested in quantifying various exoplanet
science cases that a future space-based MIR nulling interferometer could
address. Reliable signal extraction algorithms exist and our results underline
the power of the MIR wavelength range for deriving fundamental exoplanet
properties from single-epoch observations.

###### Key Words.:

Methods: data analysis – Techniques: interferometric – Techniques: high
angular resolution – Planets and satellites: detection – Planets and
satellites: terrestrial planets – Planets and satellites: fundamental
parameters

##  1 Introduction

Ever since the first detection of an exoplanet orbiting a Solar-like star
(Mayor & Queloz, 1995), thousands of exoplanets have been detected from the
ground and with dedicated space missions. Ongoing and future space missions,
such as the _CHaracterizing ExOPlanet Satellite (CHEOPS)_ (Benz et al., 2021),
the _James Webb Space Telescope (JWST)_ , and the _Atmospheric Remote-sensing
Infrared Exoplanet Large-survey (ARIEL)_ (Tinetti et al., 2018), will focus on
characterizing a subset of the known exoplanets and their atmospheres in
greater detail. All of these missions rely on investigating exoplanets that
transit in front of their host stars, but the vast majority of exoplanets do
not transit. Hence, the characterization space is limited and biased towards
close-in and large planets. In order to investigate the atmospheric properties
of hundreds of terrestrial exoplanets, including dozens that are potentially
habitable, optimized large-scale space missions that rely on a direct
detection technique will be required. This could be either a single aperture
UV/optical/NIR telescope to characterize exoplanets in reflected light (e.g.,
the _HabEx_ and _LUVOIR_ concepts; Gaudi et al., 2020; The LUVOIR Team, 2019),
or a mid-infrared (MIR) nulling interferometer to probe the thermal emission
of exoplanets. In Quanz et al. (2021) we have introduced the _Large
Interferometer For Exoplanets (LIFE)_ initiative that aims at developing the
science objectives, requirements and a concept (incl. a technology development
roadmap) for such a space-based MIR nulling interferometer mission. Of
significant importance in the current phase of the _LIFE_ initiative is the
availability of a simulator environment that, under well-defined assumptions,
can assess the performance of a certain interferometer architecture and
instrument concept in terms of scientific output. Similarly, a solid
understanding of how well a measured signal can be extracted from the
simulated data is crucial. The longer-term objective is to develop and
maintain a flexible framework that incorporates the fundamental properties of
various interferometer architectures and includes all relevant noise sources
(astrophysical and instrumental). Such a framework can then be used for even
more sophisticated scientific simulations and trade-offs (e.g., prioritization
of targets), but also to understand and quantify the impact of technical
trade-offs.

The structure of the paper is as follows: in Sect. 2, we give a brief
introduction to the basic principles of nulling interferometry needed to
describe the measuring process of _LIFE_. We further introduce the most
dominant astrophysical noise terms, explain how the signal-to-noise ratio
(SNR) of a measurement can be quantified and provide simulations of an Earth-
twin exoplanet seen at a distance of 10 pc as a specific example. In Sect. 3,
we revisit the assumption of the photon-based SNR calculations and present a
signal extraction method that enables locating exoplanets in the vicinity of
their host stars and estimating their fluxes. In Sect. 4 we investigate how
accurately physical properties (such as position and flux and, from this,
effective temperature and radius) of simulated exoplanets can be extracted
based on single-epoch observations. We do this by re-analyzing parts of the
Monte Carlo simulations presented in Quanz et al. (2021). We discuss the
results and conclude in Sect. 5. Throughout this work, we assume that the
measurement is primarily disturbed by the astrophysical noise sources and not
by noise emerging from instrumental effect. In Appendix A, we examine this
assumption by presenting estimates for the required instrument stability using
an explicit treatment of the instrumental noise.

##  2 Signal simulation

###  2.1 Nulling interferometry

Bracewell (1978) is the first to propose nulling interferometry as a method to
detect and characterize exoplanets. The general goal of a nulling
interferometer is to enable measurements of a faint source - i.e., of an
exoplanet - despite the presence of a much brighter source close by. The
simplest realization, a single Bracewell interferometer, consists of two
collector apertures which are separated by a baseline $b$. By combining the
collected light coherently, the response of the instrument is a sinusoidal
fringe pattern with spacing $\lambda/b$ projected onto the plane of the sky,
$\lambda$ being the observing wavelength (e.g., Bracewell, 1978; Lay, 2004).
By adding a $\pi$-phase shift to the beam of one of the collectors, the
central minimum of the fringe pattern coincides with the position of the star
so that its signal is effectively suppressed, i.e., “nulled”. Nulling the on-
axis stellar light is hence a way to optically separate the star light from
other nearby sources, such as exoplanets.

An interferometer has the advantage that it provides higher “spatial
resolution” than a single-aperture telescope with aperture size $D$: the first
positive transmission peak in the fringe pattern is located at $\lambda/2b$,
while the spatial resolution of a telescope is $\propto\lambda/D$.
Furthermore, in case of a free-flying space interferometer as foreseen for
_LIFE_ , the baselines are re-configurable and generally much larger (up to
hundreds of meters) than the apertures of a single-dish telescope (i.e.,
$b>>D$).

In a static, monochromatic configuration, the single Bracewell interferometer
provides only limited coverage of the uv-plane in the form of a single
baseline. Increasing this coverage will allow for the retrieval of more
complex signals from the measurements. One possible method for enhancing
coverage without changing the configuration of the array itself can be
achieved by rotation of the array (cf. rotation synthesis, Paresce & Crane
(1997)). By continuously rotating the interferometer around its line of sight,
the projected fringe pattern will also rotate, and while the on-axis star
remains nulled, any planet in its vicinity moves through the pattern of
varying transmission and its signal is modulated. The uv-plane coverage can be
additionally increased by adding collector apertures to the array and hence
raising the number of baselines, or by performing multi-wavelength
observations (see Fig. 6).

Nonetheless, the single Bracewell interferometer is insufficient for the
detection of Earth-like exoplanets (e.g., Angel & Woolf, 1997; Defrère et al.,
2010)111We refer the reader to Dandumont et al. (2020) of recent detection
yield estimates for Bracewell interferometers with various aperture sizes..
Furthermore, because the transmission pattern is symmetric, an ambiguity of
$180^{\circ}$ exists for the position angle of any detected exoplanet.

Driven by the ultimate goal to detect Earth-like exoplanets, extensive studies
of interferometer configurations consisting of more than two apertures led to
two major space mission concepts in the early 2000s: the _Darwin_ concept led
by the European Space Agency (ESA)(Cockell et al., 2009) and NASA’s
_Terrestrial Planet Finder - Interferometer (TPF-I)_ (Lawson et al., 2007).
_LIFE_ leverages the heritage of these concepts by rooting its design in
analyses performed in the context of the _Darwin_ and _TPF-I_ missions.

####  2.1.1 Beam combination in an X-array configuration

For _LIFE_ we assume an X-array architecture. It consists of four free-flying
collector spacecraft, arranged in a rectangular configuration (see Fig. 1) in
a plane perpendicular to the line of sight. Because the spacecraft are free-
flying, the baselines can be scaled depending on the desired spatial extent of
the fringe pattern. At the moment, the ratio between the length of the long
so-called “imaging baseline” and that of the short so-called “nulling
baseline” is chosen to be 6:1 (cf. Defrère et al., 2010). Compared to smaller
ratios, the 6:1 configuration is better suited for advanced post-processing
techniques aimed at removing a part of the instrumental effects on the
measurement called instability noise (Lay, 2006). Further investigations and
trade-offs are required to validate and confirm this choice for the _LIFE_
mission. The fifth spacecraft, in which the beams are combined coherently and
the signal is detected, is located in the geometric center of the array.
However, this beam combiner spacecraft could either be located in the same
plane as the collector spacecraft or, alternatively, it could fly several
hundreds of meters above this plane. This trade-off mostly affects stray light
and viewing zone considerations and while it is of no concern for the
following analyses, it does significantly affect the field of regard of the
mission (cf. Lay et al., 2007).

The beam combination for such complex array configurations need to be captured
by a suitable formalism. The following presents such a formalism as described
by Guyon et al. (2013).

Figure 1: Assumed beam combination scheme of the LIFE array. Four collector
spacecraft are located on a 6:1 rectangle with the combiner spacecraft located
in the geometric center. The shorter distance between the collectors $b$ is
called the nulling baseline, the longer distance is the imaging baseline $qb$.
The dual chopped Bracewell beam combination scheme, consisting of an
interconnection of two single Bracewell combiners, is shown. We note that the
collector spacecraft and the interferometric outputs are numbered for the
subsequent calculation and that the spacecraft diameters are not to scale.

To calculate the coherent beam combination, a nulling interferometer can be
fully described by $N$ apertures, each defined by its respective position and
radius $(x_{k},y_{k},r_{k})$ in a two-dimensional plane, where
$k\in[1,\dots,N]$. The complex field amplitude generated in the $k$-th
aperture by a point source with projected angular offset ($\alpha,\beta)$ from
the central line-of-sight is given by

$\displaystyle
V_{k}=r_{k}\,e^{\mathrm{i}2\pi(x_{k}\alpha+y_{k}\beta)/\lambda}\quad,$ (1)

where $\lambda$ is the wavelength. The interferometric combinations can then
be represented by a matrix $\@vec{\rm U}$, describing the linear mapping from
input amplitude vector $\vec{V}$ to the output amplitude vector

$\displaystyle\vec{W}=\@vec{\rm U}\vec{V}.$ (2)

Figure 2: Left and middle: Transmission maps $T_{\textit{m}}$ of output 3 and
4 calculated for a wavelength $\lambda=$10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$, a nulling baseline $b=$12\text{\,}\mathrm{m}$$ and a 6:1 baseline ratio.
These two outputs have destructive interference at their vertical center and
are used to create the differential map. The axes show the distance from the
center in milliarcseconds. The white star indicates the position of the star
in the center of the map. Right: Differential map $T_{\mathrm{dif}}$, which is
the difference between $T_{\mathrm{3}}$ and $T_{\mathrm{4}}$. The red dashed
line shows the path of a hypothetical planet at
$100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ through the differential map as
the interferometer array describes a full rotation.

Based on the heritage of the _Darwin_ and _TPF-I_ concepts, for this work we
assume that light picked up by the collector spacecraft is merged using a dual
chopped Bracewell combiner (Lay, 2004). The matrix

$\displaystyle\@vec{\rm U}=\frac{1}{\sqrt{4}}\pmatrix{0}$ 0 (3) (4) (5) (6)

describes a realizable (matrix is unitary, cf. Guyon et al. (2013)) and
lossless (matrix is orthonormal, cf. Laugier et al. (2020)) implementation of
this combiner.

The rows of this matrix reflect the two-stage approach of the beam combination
(Fig. 1). In a first stage, the collector apertures are grouped and combined
into two single Bracewell combiners. The respective constructive outputs,
denoted as Outputs 1 and 2, are captured by the first two rows of the matrix
$\@vec{\rm U}$, which do not induce any phase delays. The destructive outputs
of the first stage are expressed by the $e^{i\pi}$ phase difference between
entries one/two or three/four in the last two rows of the combiner matrix. In
a second stage, these destructive outputs are combined again, with one of the
destructive outputs receiving an additional phase shift of $\frac{\pi}{2}$.
This produces two complementary nulled outputs that can be used to apply phase
chopping to the final output, which reduces the susceptibility of the
measurement to instrumental noise effects (Lay, 2004). The outputs of this
second combination stage are denoted as the destructive Outputs 3 and 4 and
are captured by the rows 3 and 4 in the matrix $\@vec{\rm U}$. We note that
the analyses presented in this paper are solely based on these destructive
outputs.

####  2.1.2 Transmission and differential maps

The two-dimensional projected intensity transmission $T_{\mathrm{m}}$, which
describes the amount of signal originating from a point in the sky at position
$(\alpha,\beta)$ that transmits to the output m, can be described by a
combination of the input beams $V_{k}$ as

$\displaystyle
T_{\mathrm{m}}=|W_{\mathrm{m}}|^{2}=\bigl{\lvert}\sum_{k}U_{\mathrm{m},k}V_{k}\bigr{\rvert}^{2}.$
(7)

In the following, the transmission map of one of the interferometric outputs
is calculated analytically to illustrate the basic principle. In a coordinate
system that has its origin in the center of the interferometer array, the
aperture positions $x_{k}$ and $y_{k}$ are given as $\pm L$ and $\pm qL$
respectively, with $L$ representing half the nulling baseline $b$, and $q$
being the ratio between the imaging baseline and the nulling baseline. The
transmission of output 3 as a function of angular separation $(\alpha,\beta)$
is given by

$\displaystyle T_{3}=$
$\displaystyle\,|W_{3}|^{2}=\bigl{\lvert}\sum_{k}U_{3,k}V_{k}\bigr{\rvert}^{2}$
(8) $\displaystyle=$
$\displaystyle\,\dfrac{1}{4}\,\bigl{\lvert}\,e^{\mathrm{i}2\pi(L\alpha+qL\beta)/\lambda}-e^{\mathrm{i}2\pi(-L\alpha+qL\beta)/\lambda}$
(9) $\displaystyle-e^{\mathrm{i}2\pi(-L\alpha-
qL\beta)/\lambda+\mathrm{i}\pi/2}+e^{\mathrm{i}2\pi(L\alpha-
qL\beta)/\lambda+\mathrm{i}\pi/2}\,\bigr{\rvert}^{2}$ (10) $\displaystyle=$
$\displaystyle\,4\sin^{2}\left(\dfrac{2\pi
L\alpha}{\lambda}\right)\cos^{2}\left(\dfrac{2\pi
qL\beta}{\lambda}-\dfrac{\pi}{4}\right).$ (11)

To apply phase chopping, outputs 3 and 4 are subtracted from each other. This
yields the differential map

$\displaystyle T_{\mathrm{dif}}$ $\displaystyle=T_{3}-T_{4}$ (12)
$\displaystyle=4\sin^{2}\left(\dfrac{2\pi
L\alpha}{\lambda}\right)\sin\left(\dfrac{4\pi qL\beta}{\lambda}\right).$ (13)

This map is anti-symmetric with respect to its central point. Thus, any point
symmetric emission source (e.g., local zodiacal light, homogeneous exozodi
disks) will not transmit through the differential map, but will only
contribute to the statistical shot noise (e.g., Defrère et al., 2010). An
example of the two destructive output transmission maps $T_{3}$ and $T_{4}$ as
well as the differential map is shown in Fig. 2.

Figure 3: Top: Differential map in polar coordinates, normalized from $-1$
(blue) to 1 (yellow). The red dashed line at
$100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ angular separation corresponds to
the path of the planet shown in the right panel of Fig. 2 for a full rotation
of the array. Bottom: The normalized signal transmission intensity of the
planet as a function of rotation angle.

####  2.1.3 Signal modulation

To study the signal transmission of a potential planet, the differential map
is examined in polar coordinates as shown in the top panel of Fig. 3. Any
(static) source at some angular separation $\theta$ from the central star will
move through this map on a horizontal line as the telescope array is rotated
by an angle $\phi$. Thus, the normalized signal generated by such a source
over one array rotation corresponds to the value of the differential map on a
horizontal line through the map at the corresponding radial distance. The
normalized signal transmission as a function of rotation angle, i.e., the
modulated signal, of that point source is shown in the bottom panel of Fig. 3.

The modulation efficiency $\xi$ indicates the part of the incoming signal that
can be used for signal extraction and corresponds to the root-mean-square
(rms) of the differential map upon a complete array rotation (Lay, 2004). The
modulation efficiency over a full array rotation as a function of $\theta$ is
given by

$\displaystyle\xi(\theta)=\sqrt{\langle
T_{\mathrm{dif}}(\theta,\phi)^{2}\rangle_{\phi}}\quad.$ (14)

Because the spatial extent of the fringe pattern depends on the baseline
configuration of the interferometer and the considered wavelength, this is
also true for the modulation efficiency. Figure 4 shows $\xi(\theta)$ for
three different MIR wavelengths for the baseline configuration used in Fig. 2.
The curve for $\lambda=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, which is
based on the differential map shown in Figs. 2 and 3, has its first peak and
global maximum at $100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$. For shorter
wavelengths the pattern shrinks and the peak transmission is located closer to
the star; for longer wavelengths the pattern expands and the peak transmission
moves away from the star. The position of the maximum of $\xi(\theta)$ can be
calculated numerically as

$\displaystyle\theta_{\mathrm{\xi_{max}}}$
$\displaystyle=0.59\,\dfrac{\lambda}{2L}=0.59\,\dfrac{\lambda}{b}\quad,$ (15)

which is slightly more than the separation corresponding to the first maximum
in the transmission maps at $\lambda/2b$.

Figure 4: Modulation efficiency $\xi(\theta)$ as a function of angular
separation for the three wavelengths $\lambda=$
$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The vertical red line indicates an
angular separation of $100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$
corresponding to the separation of the planet shown in the right panel of Fig.
2.

###  2.2 Astrophysical sources

In the following we describe how various relevant astrophysical sources and
their photon flux are implemented in LIFEsim.

####  2.2.1 Host star and exoplanets

To first order, thermal radiation emitted by an exoplanet and its host star
can be approximated as black body radiation. Because we will base our SNR
calculation on the principles of photon statistics, we are interested in the
received spectral photon flux in units of
[$\mathrm{p}\mathrm{h}\,\mathrm{s}^{-1}\mathrm{m}^{-2}\mathrm{m}^{-1}$]

$\displaystyle
F_{\lambda}(\lambda,T)=\dfrac{2c}{\lambda^{4}}\dfrac{1}{e^{hc/\lambda
k_{\mathrm{B}}T}-1}\dfrac{\pi R^{2}}{d^{2}}\quad,$ (16)

where $h$ is Planck’s constant, $c$ is the speed of light, $k_{\mathrm{B}}$ is
the Boltzmann constant, $T$ and $R$ are the effective temperature and radius
of either the planet or the star, and $d$ is the distance from the observer.

Throughout this paper the star is modeled as a black body with effective
temperature $T_{\mathrm{s}}$ and radius $R_{\mathrm{s}}$. Limb darkening is
not taken into account, and thus the stellar brightness is distributed
uniformly across the projected circular disk. Because of the principle of
geometric stellar leakage, as described further below in Sect. 2.2.2, this is
thought to be a conservative approximation.

As one of the objectives of _LIFE_ is the atmospheric characterization of
(terrestrial) exoplanets by measuring wavelength-dependent absorption features
in their emission spectra, modeling their emission with black body radiation
is a simplification. Therefore, LIFEsim has the option to select either a
simple black body model for the planet or to upload a simulated / observed
emission spectrum. In the following subsections we make use of both options.
When using a black body, we denote the radius and the effective temperature of
the exoplanets with $R_{\mathrm{p}}$ and $T_{\mathrm{p}}$. We use a simulated
Earth’s spectrum for illustrative purposes when calculating the expected SNR
of an Earth-twin in Sect. 2.4.

####  2.2.2 Stellar leakage

Even with a perfect nulling interferometer, photons from the observed star
“leak” through the instrument. If the central null fringe is positioned on the
center of the star, the stellar flux is not completely suppressed because of
the star’s finite size. This phenomenon is described as geometric stellar
leakage (e.g., Lay, 2004). The photon noise contribution from geometric
stellar leakage is determined by applying the transmission maps to the surface
area of the star. As such, it depends on the surface brightness and angular
size of the star as well as on the baseline configuration.

####  2.2.3 Local zodiacal dust

Our Solar System contains diffuse dust that is mostly located near the
ecliptic plane. The dust adds a radiation background to any astronomical
observation either coming from sunlight that is scattered off the dust
particles or as thermal radiation emitted by the dust itself. As _LIFE_ works
in the MIR, we concentrate solely on the thermal emission. The local zodiacal
surface brightness is described by a sky position dependent model developed
for the DarwinSim science simulator (Den Hartog & Karlsson, 2005). The
parametrization, which is based on data from COBE (Kelsall et al., 1998),
gives the spectral surface brightness

$\displaystyle\begin{split}I_{\lambda}(\lambda,\lambda_{\mathrm{rel}},\beta)=&\,\tau\,\left[B_{\lambda}(\lambda,T_{\mathrm{eff}})+A\cdot
B_{\lambda}(\lambda,T_{\odot})\left(\dfrac{R_{\odot}}{$1.5\text{\,}\mathrm{A}\mathrm{U}$}\right)^{2}\right]\\\
&\cdot\left[\dfrac{\pi/\arccos(\cos(\lambda_{\mathrm{rel}})\cos(\beta))}{\sin^{2}(\beta)+0.36\cdot\left(\dfrac{\lambda}{$11\text{\,}\mathrm{\SIUnitSymbolMicro
m}$}\right)^{-0.8}\cos^{2}(\beta)}\right]^{\frac{1}{2}}\end{split}$ (17)

where $\lambda$ is the wavelength,
$\lambda_{\mathrm{rel}}=\lambda_{\mathrm{ecl}}-\lambda_{\mathrm{ecl,\sun}}$ is
the ecliptic longitude relative to the Sun, $\beta$ is the ecliptic latitude,
$\tau=$4\text{\times}{10}^{-8}\text{\,}$$ is the optical depth towards the
ecliptic poles, $B_{\lambda}(T)$ is the Plank function,
$T_{\mathrm{eff}}=$265\text{\,}\mathrm{K}$$ is the effective temperature of
the local zodiacal dust cloud at $1\text{\,}\mathrm{A}\mathrm{U}$,
$T_{\odot}=$5778\text{\,}\mathrm{K}$$ is the effective temperature of the Sun,
$A=0.22$ is the near-IR dust albedo, and $R_{\odot}$ is the radius of the Sun.
Since the local-zodiacal light received by the interferometer is diffuse, it
cannot be brought to interfere destructively. Therefore, it is not possible to
remove the local-zodi emission by virtue of the rotating nulling
interferometer. Subsequently, the most effective way of minimizing the noise
at long wavelengths is to observe in directions where the local-zodi emission
is weak (see Appendix C).

Compared to the original COBE data, the DarwinSim model overestimates the flux
from the zodiacal dust in the
$6\text{\,}\mathrm{-}$$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ range and for
a line of sight with a relative latitude of more than $90^{\circ}$ from the
Sun by $10-20\,\%$. For wavelengths below
$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ the difference increases to a factor
of 2-3 at $3\text{\,}\mathrm{\SIUnitSymbolMicro m}$. However, at these
wavelengths the photon noise contribution by the zodiacal dust is several
orders of magnitude lower than the photon noise from the stellar leakage as we
will show below in Sect. 2.3.

####  2.2.4 Exozodiacal dust

Exozodiacal dust clouds (“exozodis” for short) are the equivalent to the local
zodiacal dust in other stellar systems. Even if the dust density is as low as
in the Solar system, integrated over the full field-of-view the MIR radiation
emitted by the dust is two to three orders of magnitude higher than that of an
Earth-like planet (e.g., Defrère et al., 2010) and can thus significantly
affect the integration time required to detect (terrestrial) exoplanets.

The exozodiacal dust distribution and its surface brightness are simulated
according to the model described by Kennedy et al. (2015). It assumes a dust
distribution in exoplanetary systems similar to the one in the Solar system,
except for a global scaling factor $z$, the number of zodis222A probability
distribution for $z$ was derived by the HOSTS survey carried out at the Large
Binocular Telescope (LBTI) (Ertel et al., 2020).. The model assumes that the
dust emission is optically thin and the dimensionless, face-on surface density
$\Sigma(r)$ is approximated by a power law distribution. The emission is
modeled as black body radiation with the temperature-radius profile of the
disk following the equilibrium temperature

$\displaystyle T(r)=$278.3\text{\,}\mathrm{K}$\cdot
L_{\mathrm{s}}^{\frac{1}{4}}\,r^{-\frac{1}{2}}$ (18)

where $r$ is given in AU and $L_{\mathrm{s}}$ in Solar luminosities. The
spectral surface brightness of the disk as seen face-on is then given by

$\displaystyle I_{\lambda}(\lambda,r)=\Sigma(r)\,B_{\lambda}(\lambda,T(r)).$
(19)

Due to the high temperature and surface density close to the star, much of the
emitted radiation originates from the central regions.

It is important to note that the given definition of the zodi level $z$ is
coupled to the surface density and not to the spectral flux or the total
luminosity of the disk. For constant $z$, the total emitted radiation at a
given wavelength, as well as the luminosity of the disk, scale with the
surface area of the disk and thus with the stellar luminosity
$L_{\mathrm{disk}}\propto r_{0}^{2}\propto L_{\mathrm{s}}$. For high-
luminosity stars, which are expected to have larger exozodi disks in this
model, the total emitted flux by the disk is thus also larger than for low-
luminosity stars.

At the moment, the exozodi dust disks are assumed to be homogeneous and
symmetric around the star and always viewed face-on. Therefore, they only
contribute to the shot noise (see, Sect. 2.3 below). The influence of
asymmetric or clumpy exozodiacal disks viewed at an inclination $i>0^{\circ}$
was investigated in Defrère et al. (2010) and Defrère et al. (2012) and we
discussed these points and their potential impact on exoplanet detection yield
in Quanz et al. (2021).

###  2.3 Signal-to-noise ratio calculations

####  2.3.1 Astrophysical noise sources

The astrophysical sources located within the field-of-view of the collector
apertures create a scene that is described by the surface brightness
$I(\vec{\theta},\lambda)$. It depends on the position $\vec{\theta}$ within
the field-of-view and on the wavelength, and is here assumed to be time-
independent.

When the photons emitted by the different sources reach the interferometer,
they are detected in one of the interferometric outputs. The collected signal
$S_{\mathrm{m}}$ by output m, described by transmission map $T_{\mathrm{m}}$,
after some integration time is given by the integral of the total surface
brightness $I(\vec{\theta},\lambda)$ superimposed with the transmission map
over the complete field-of-view

$\displaystyle S_{\mathrm{m}}(\lambda)=\int
T_{\mathrm{m}}(\vec{\theta},\lambda)\,I(\vec{\theta},\lambda)\,t\,A\,\eta\,\mathrm{d\Omega}$
(20)

with $t$ the integration time, $A$ the total collecting area, and $\eta$ a
detection efficiency factor combining multiple effects. The efficiency factor
$\eta$ and the size of the field-of-view are generally also wavelength-
dependent. In this work, $\eta$ is the product of the quantum efficiency of
the detector $\eta_{\mathrm{QE}}$ and the overall instrument throughput
$\eta_{\mathrm{t}}$. Until a more concrete instrument design is available, it
is assumed that both parameters are wavelength-independent. As motivated in
more detail in Quanz et al. (2021), we chose $\eta_{\mathrm{QE}}=0.7$ and
$\eta_{\mathrm{t}}=0.05$ as default values. The effective field-of-view is
assumed to be $\mathrm{FoV}=\lambda/D$ in diameter as we assume the light will
be coupled into single-mode fibers.

The spectral signal generated by a single planet, approximated as a point
source at location $\vec{\theta}_{\mathrm{p}}$, is thus given by

$\displaystyle
S_{\mathrm{p,m}}(\lambda)=T_{\mathrm{m}}(\vec{\theta}_{\mathrm{p}},\lambda)\,F_{\mathrm{p}}(\lambda)\,t\,A\,\eta$
(21)

where $F_{\mathrm{p}}(\lambda)$ is the planet flux. Over a full array
rotation, the quadratic mean of the modulated signal
$S_{\mathrm{p}}=S_{\mathrm{3,p}}-S_{\mathrm{4,p}}$ detected from the planet is
given by

$\displaystyle\sqrt{\langle
S_{\mathrm{p}}^{2}(\lambda)\rangle}=\xi(\lambda,\theta_{\mathrm{p}})\,F_{\mathrm{p}}(\lambda)\,t\,A\,\eta$
(22)

where $\xi(\lambda,\theta)$ is the wavelength dependent modulation efficiency
at angular separation $\theta$ from the star (see, Eq. (14)).

If only rotationally symmetric sources, such as the central star, the local
zodiacal background or a homogeneous, face-on exozodiacal dust cloud are
considered, the detected signal does not depend on the rotation angle of the
array. If the source is only point-symmetric, such as a smooth exozodiacal
disk with some inclination $i>0^{\circ}$, the signal per output varies with
rotation angle, but is always equal for the two destructive outputs 3 and 4.
Thus, point-symmetric sources do not contribute to the modulated signal.

However, as Mugnier et al. (2006) pointed out, the linear combination of
detected signals is an incoherent combination, since it is performed
numerically after detection (i.e., in post processing). Removing the
contribution of symmetrically distributed sources by incoherent combination is
therefore only effective for the signal part of the data, but the sources
still contribute to the statistical noise. Hence, the SNR per spectral
wavelength bin, denoted as SNRλ and defined as the ratio between detected
exoplanet photons and detected photons from the various noise terms, is
calculated as

$\displaystyle\mathrm{SNR}_{\lambda}=\dfrac{\int\sqrt{\langle
S_{\mathrm{p}}^{2}(\lambda)\rangle}\,\mathrm{d}\lambda}{\sqrt{2\int\,(S_{\mathrm{sym,3}}(\lambda)+\,\sqrt{\langle
S_{\mathrm{p,3}}^{2}(\lambda)\rangle})\,\mathrm{d}\lambda}}\quad,$ (23)

where $S_{\mathrm{sym,3}}$ is the contribution from symmetric background
sources to the signal in output 3 and $S_{\mathrm{p,3}}$ is the signal of the
planet only333For the computation it does not matter whether we use output 3
or output 4 in the denominator. The factor of 2 in front of the integral
ensures that the contributions from both outputs are considered.. The integral
runs over the wavelength range covered by the wavelength bin.

Following Lay (2004), for a measurement across several wavelength bins and
under the assumption that the noise is uncorrelated between the bins, the
integrated SNR over the full wavelength range, or summed over all wavelength
bins, is given by

$\displaystyle\mathrm{SNR_{tot}}$
$\displaystyle=\sqrt{\sum_{\lambda}\mathrm{SNR_{\lambda}^{2}}}\propto\sqrt{t\,A\,\eta}\quad,$
(24)

which scales, as expected, with the square root of integration time, area and
detection efficiency.

In LIFEsim each noise term from the different sources is calculated
individually to avoid numerical discretization errors, as the angular extents
of the sources have different scales. For the local and exozodiacal dust, the
transmission map as well as the surface brightness distribution are calculated
and integrated over a two-dimensional artificial image covering the full
field-of-view defined by the instrumental parameters. For the stellar leakage,
it is only integrated over the solid angle covered by the stellar disk, but
with much higher resolution. The planet signal transmission efficiency is
calculated analytically using Eq. (21). All terms are calculated for each pre-
defined wavelength bin, with varying angular extent of the field-of-view and
possibly varying bin width, and are then combined to an integrated signal-to-
noise ratio as given in Eq. (24).

The calculated SNR depends strongly on the baseline of the array. The mean
modulated planet signal (Eq. (22)) scales with the modulation efficiency,
which in turn depends on the baseline as demonstrated in Sect. 2.1.3. For the
contribution from background sources, the amount of stellar leakage
contribution depends on the broadness of the central null. This broadness of
the null in the transmission map is also governed by the baseline
configuration. Thus, maximizing the SNR by changing the length of the
baselines poses a trade-off between the amount of planetary signal received
over the amount of stellar leakage affecting the measurement. Ideally, to
reduce the integration time required for detection, the baselines of the
interferometer should thus be chosen such that the SNR across the full
wavelength range is high for projected separations that are of major interest.
For instance, maximizing the chances of detecting an Earth-like planet
orbiting within the habitable zone (HZ) around nearby stars can be achieved by
maximizing the SNR across the projected habitable zone. As discussed in Quanz
et al. (2021), for _LIFE_ it is currently assumed that the baselines can
indeed be reconfigured depending on the target star where the minimum and
maximum separation between two spacecraft is $10\text{\,}\mathrm{m}$ and
$600\text{\,}\mathrm{m}$ respectively.

####  2.3.2 Instrumental noise terms

Instrument perturbations such as intensity variations and optical path
difference (OPD) errors can degrade the null and lead to additional stellar
leakage as well as instability noise (Lay, 2004; Defrere, 2009). Similarly,
detector noise (such as, e.g., dark current) and thermal background noise from
the instrument can also deteriorate the measurement. In the current version of
LIFEsim, it is assumed that the instrument will be designed such that photon
shot noise will dominate over these additional instrumental noise sources.
Appendix A presents a preliminary analysis of the aforementioned instrumental
effects and constrains the maximum level of allowed perturbations to operate
in this fundamental noise limited regime.

###  2.4 A simple case study: Earth-twin at 10 pc

We apply the above described SNR calculation to the example of an Earth-like
planet located around a Sun-like star at $10\text{\,}\mathrm{p}\mathrm{c}$,
with the aim to illustrate the effects of the different noise sources on the
SNR and their wavelength dependence. Instrument parameters and properties of
the simulated sources are given in Table 1.

Figure 5 (top panel) shows the total fluxes received from the different
sources as solid lines. The “leakage” terms, i.e., the part of the flux that
is not nulled, are shown as dotted lines. These leakage terms contribute to
the shot noise. It can be seen that at short wavelengths the stellar leakage
dominates strongly, while for longer wavelengths the emission from the local
zodiacal dust contributes most to the shot noise.

In the bottom panel of Fig. 5 we show the detected signal from the planet
(averaged over each spectral bin) as well as the shot noise and the SNR for
the parameters listed in Table 1. Integrating over the full wavelength range
gives an SNR$\approx$9.7. The maximum SNR per bin of $\approx$3 is reached
around $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$. For wavelengths shorter
than $8\text{\,}\mathrm{\SIUnitSymbolMicro m}$ the SNR decreases rapidly, due
to a decrease in planet signal and the increasing shot noise due to stellar
leakage.

Table 1: Instrumental parameters and properties of the simulated sources for
the SNR calculation of an Earth-like planet at 10 pc.

Parameter | Value | Description
---|---|---
$D$ | $2\text{\,}\mathrm{m}$ | Aperture diameter
$t$ | $200\,000\text{\,}\mathrm{s}$ | Integration time
$\eta_{QE}$ | 0.7 | Quantum efficiency
$\eta_{t}$ | 0.05 | Instrument throughput
FoV | $\lambda/D$ | Field-of-view ($\sim 1^{\prime\prime}@\lambda=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$)
| $4-18.5\,$\mathrm{\SIUnitSymbolMicro m}$$ | Wavelength range
$R$ | 20 | Spectral resolution ($\lambda/\Delta\lambda$)aa$a$Spectral resolution is assumed to be constant across the full wavelength range, such that the bin width $\Delta\lambda$ increases for larger wavelengths. For the parameters listed here, this results in 31 spectral bins.
$b$ | $15\text{\,}\mathrm{m}$ | Nulling baselinebb$b$The baseline is set according to Eq. (15) evaluated at $15\text{\,}\mathrm{\SIUnitSymbolMicro m}$.
$r$ | $6:1$ | Array baseline ratio
$d$ | $10\text{\,}\mathrm{p}\mathrm{c}$ | Target distance
$\theta_{\mathrm{p}}$ | $100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ | Planet-star separation
$R_{\mathrm{p}}$ | $1\,R_{\oplus}$ | Planet radius
$T_{\mathrm{p,0}}$ | $285\text{\,}\mathrm{K}$ | Planet surface temperaturecc$c$In Fig. 5, instead of assuming blackbody emission for the planet, we used the radiative transfer atmospheric model code petitRADTRANS to compute an MIR emission spectrum corresponding to an average cloud-free modern Earth spectrum (Konrad et al., 2021)
$R_{\mathrm{s}}$ | $1\,R_{\sun}$ | Stellar radius
$T_{\mathrm{s}}$ | $5778\text{\,}\mathrm{K}$ | Stellar effective temperature
$\lambda_{\mathrm{rel}},\,\beta$ | $135^{\circ},\,45^{\circ}$ | Ecliptic coordinates
$z$ | $3$ | Level of exozodi emissiondd$d$This value corresponds to the median of the best-fit nominal model derived from the HOSTS survey (Ertel et al., 2020).
444The table corresponds to Figure 5.

Figure 5: Top: Comparison of the various input signals in our example of an
Earth-like planet at $10\text{\,}\mathrm{p}\mathrm{c}$. The solid lines
indicate the total flux received from the respective sources: planet (red),
star (black), exozodiacal dust (grey) and local zodiacal dust (blue). The
dotted lines with the same color code indicate the fraction of the flux that
contributes to the shot noise. Bottom: Detected planet signal (red line) and
total shot noise (blue dashed line) averaged over each wavelength bin and for
an integration time of $200\,000\text{\,}\mathrm{s}$ (values refer to the
y-axis on the left). The resulting spectral SNR is shown as the black line and
is calculated as described in Sect. 2.3 (values refer to the y-axis on the
right). Integrated over the full wavelength range one obtains an
SNR$\approx$9.7.

##  3 Signal extraction

The approach and the simulated spectrum presented in the previous section are
based on photon counting statistics. It was implicitly assumed that the signal
extraction can be performed with sufficient accuracy.

In practice, data processing and signal extraction algorithms will be of great
importance, in particular if a mission like _LIFE_ features an initial search
phase555The baseline mission scenario features a 2.5 year search phase for
detecting previously unknown exoplanets, followed by a 2.5-3.5 year
characterization phase (Quanz et al., 2021).: in addition to separating the
planetary signal from the various noise sources, it will be crucial to
accurately derive physical characteristics of the exoplanet (e.g., radius,
effective temperature and separation from the host star) from single epoch
data in order to identify and rank-order the most interesting objects for in-
depth follow-up observations during the characterization phase. For both the
Darwin and the TPF-I missions, multiple signal analysis algorithms were
proposed. Draper et al. (2006) presented an algorithm which was based on an
advanced correlation process. For Darwin, Mugnier et al. (2006) and Thiébaut &
Mugnier (2005) presented a signal extraction scheme based on the maximum-
likelihood method (MLM). While in the meantime significant progress has been
made in many areas of data post-processing, especially thanks to machine
learning based approaches, the before-mentioned methods still offer a
straight-forward framework to deal with the unique characteristics of nulling
interferometry with respect to the modulation of the planet signal.

In the following subsections the maximum-likelihood method will be described,
modified, and applied to simulated _LIFE_ data. The goal is to find a metric
to quantify the robustness of an exoplanet detection and compare the extracted
SNR to that from the photon statistics used in the previous section.
Additionally, the signal extraction for multi-planet systems will be
investigated. Furthermore, we will revisit the Monte Carlo simulations
presented in Quanz et al. (2021) and quantify how accurately we can expect
exoplanet properties to be determined from single epoch data during the search
phase.

###  3.1 Maximum-likelihood method

As described in Sect. 2.3, the signal measurement with a nulling
interferometer can be described as a noisy time series. Following Mugnier et
al. (2006) and Thiébaut & Mugnier (2005), it is modeled as

$\displaystyle A_{t,\lambda}$
$\displaystyle=F_{\mathrm{p},\lambda}\,T_{t,\lambda}(\vec{\theta}_{\mathrm{p}})+n_{t,\lambda}$
(43)

where $A_{t,\lambda}$ is the recorded amplitude at time $t$ and for effective
wavelength $\lambda$, $T_{t,\lambda}(\vec{\theta}_{\mathrm{p}})$ is the
response of the instrument (i.e., the value of the differential map
$T_{\mathrm{dif}}$) at the location of the planet
$\vec{{\theta}}_{\mathrm{p}}$ as a function of time and wavelength,
$F_{\mathrm{p},\lambda}$ is the discretized planet flux and $n_{t,\lambda}$
denotes the noise, which is assumed to be spectrally independent normal noise
and whose variance $\sigma^{2}(t,\lambda)$ can be estimated from the data. In
a first step and as a simplification, only the case with one (detectable)
planet per system is considered to derive the signal analysis; afterwards, the
approach will be extended to multiple planets.

Figure 6: Left: Cross-correlation of the modulated signal of a noiseless point
source located at [-100, 0] mas with the transmission map at a single
wavelength. The resulting image can be understood as the demodulated signal
function of such a point source at the considered wavelength. The image is
anti-symmetric with respect to the image center. Right: The combined cross-
correlation across the full wavelength range considered for the ideal point
source. Considering a broader wavelength range reduces the side lobes and
allows for a better location of the point source.

The maximum likelihood method enables estimating the planet position
$\vec{\theta}_{\mathrm{p}}$ and spectral flux $F_{\mathrm{p}}(\lambda)$ from
the modulated signal by searching
$(\hat{\vec{\theta}}_{\mathrm{p}},\hat{F}_{\mathrm{p}})$ that maximize the
likelihood $L(\vec{\theta}_{\mathrm{p}},F_{\mathrm{p}})$. Maximizing the
likelihood is equivalent to minimizing the negative log-likelihood or cost
function

$\displaystyle
J(\vec{\theta}_{\mathrm{p}},\vec{F}_{\mathrm{p}})=\sum_{t,\lambda}\dfrac{|A(t,\lambda)-F_{\mathrm{p}}(\lambda)\,T(t,\lambda,\vec{\theta}_{\mathrm{p}})|^{2}}{\sigma^{2}(t,\lambda)}\quad,$
(44)

which measures the discrepancy between the data $A$ and the model of the data
for the estimated set of parameters
$(\vec{\theta}_{\mathrm{p}},\vec{F}_{\mathrm{p}})$, based on the differential
map $T$. $\vec{F}_{\mathrm{p}}$ denotes the set of flux values at all measured
wavelengths. It has been shown in the previously mentioned papers that for a
given planet position $\vec{\theta}_{\mathrm{p}}$, the optimal spectral flux
is obtained by maximizing $J$ (as defined in Eq. (65)) with respect to all
$F(\lambda$):

$\displaystyle\dfrac{\partial
J(\vec{\theta}_{\mathrm{p}},\vec{F}_{\mathrm{p}})}{\partial
F(\lambda)}=0\quad\forall\lambda.$ (45)

To ensure the result is physical, a positivity constraint has to be introduced
on the estimated flux. The most likely planet position is then found by
globally maximizing the cost function over a grid by inserting the optimized
flux values back into the cost function. The estimation of the planet position
and wavelength dependent flux can be further improved if an agreement of the
estimated parameters with a priori information is enforced (Thiébaut &
Mugnier, 2005). This can be done by introducing a regularization term
$J_{\mathrm{prior}}$, which estimates the roughness of a sampled planet
spectral flux:

$\displaystyle
J_{\mathrm{prior}}(\vec{\theta}_{\mathrm{p}},\vec{F}_{\mathrm{p}})=\mu\sum_{\lambda}\left(\dfrac{\partial^{m}F_{\mathrm{p}}(\lambda)}{\partial\lambda^{m}}\right)^{2}\quad,$
(46)

where $\mu$ is a parameter allowing to tune the relative weight of
regularization. The $m$-th order derivative is computed by finite differences
(Thiébaut & Mugnier, 2005). To smooth the estimated flux, only second order
derivatives are taken into account ($m$=2). The calculation of the optimal
flux values and the cost function with and without regularization is further
outlined in Appendix A. In accordance with the previous works, we denote by
$J^{\prime}$ the cost function that is optimized with respect to the fluxes,
and by $J^{\prime\prime}$ the cost function that has an additional positivity
constraint on the flux. The cost functions can be computed on a grid of
possible initial planet positions spanning the full or only a part of the
field-of-view. The most likely planet position
$\vec{\hat{\theta}}_{\mathrm{p}}$ is then the location of the maximum of
$J^{\prime\prime}(\vec{\theta}_{\mathrm{p}})$ on the grid.

It is noted here that the cost function is closely related to the cross
correlation of an ideal modulated signal generated by a noiseless point source
with the template functions (the transmission maps). The cross correlation
across a two-dimensional grid gives the demodulated signal function of the
interferometer (Defrère et al., 2010). This is shown in Fig. 6 for a single
wavelength and for the combination of multiple wavelength bins by addition of
the maps.

###  3.2 Detection criterion

To create a link between the outcome of the maximum likelihood signal
extraction method and a quantitative detection criterion we further
investigate the noise propagation and the probability of false alarm in the
following. As introduced in Sect. 2.3, and following the definition of Mugnier
et al. (2009), the SNR of the estimated planet signal at a single wavelength
is given by

$\displaystyle\mathrm{SNR_{\lambda}}(\vec{\theta}_{\mathrm{p}})=\mathrm{SNR}[\hat{F}_{\mathrm{p}}(\vec{\theta},\lambda)]$
$\displaystyle\hat{=}\dfrac{\hat{F}_{\mathrm{p}}(\vec{\theta},\lambda)}{\sigma(\hat{F}_{\mathrm{p}}(\vec{\theta},\lambda))}.$
(47)

with $\sigma(\vec{\hat{F}}_{\mathrm{p}}(\vec{\theta}))$ the standard deviation
of the estimated flux. Given the high photon flux values of the astrophysical
sources introduced in Sect. 2.2 (cf. Fig. 5), the noise in a measurement can
be approximated by a normal distribution. To get an estimated value for the
integrated SNR, the individual wavelength bins could be combined naively by

$\displaystyle\mathrm{SNR_{tot}}=\sqrt{\sum\mathrm{SNR_{\lambda}}^{2}}=\sqrt{\sum\left(\hat{F}_{\mathrm{\lambda}}/\hat{\sigma}_{\lambda}\right)^{2}},$
(48)

as is done to calculate the integrated SNR from the predicted photon noise in
Sect. 2.3. However, this does not take into account the random nature of the
measurement, i.e., in reality, even repeating the exact same observations will
result in a slightly different estimated SNR.

Flasseur et al. (2020) presents a method that allows us to derive a detection
criterion combining measurements at multiple wavelengths with uncorrelated
noise on the basis of the cost function for direct imaging methods. With some
modifications, this approach can also be applied to nulling interferometry. We
note that the dominant and normally distributed background noise in the
modulated signal propagates as a linear combination to the SNRλ maps. This
means that in the absence of point sources the distribution of SNRλ values
over a large enough field-of-view also follows a normal distribution. Thus,
the random values of the combined cost function $J^{\prime}$ map are expected
to follow a $\chi^{2}$-distribution $p(J^{\prime})=\chi^{2}_{L}$ with $L$
degrees of freedom, where $L$ is the number of wavelength bins (Flasseur et
al., 2020). For the simulated measurements presented above, a spectral
resolution of $R=20$ and a wavelength range from
$4\text{\,}\mathrm{\SIUnitSymbolMicro
m}18.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is used (unless otherwise
noted), which results in $L=31$ wavelength bins (with varying bin widths).

Figure 7 (top-left) shows the cost function $J^{\prime}$ calculated for a
simulated measurement without any point sources. The non-uniformity of the
image is a result of the statistical fluctuations of the background noise. The
bottom panel in Fig. 7 shows the corresponding distribution of the
$J^{\prime}$ values in the grey histogram. The simulated data follow the
$\chi^{2}$-distribution with $L=31$ degrees of freedom (red line) very
closely. Following the derivation by Flasseur et al. (2020), the probability
of false alarm (PFA) in a data set containing no point sources but only
background noise is given by $\mathrm{PFA}=\mathrm{p}(\chi^{2}_{L}>\eta)$,
where $\eta$ is the detection threshold, which has to be determined. To obtain
a PFA corresponding to a 5-$\sigma$ confidence level, one needs to solve

$\displaystyle F_{L}(\eta)=\dfrac{\gamma(L/2,\eta/2)}{\Gamma(L/2)}=\Phi(5)$
(49)

for $\eta$. Here, $F_{L}$ is the cumulative distribution function of the
$\chi^{2}_{L}$ distribution, which can be expressed in terms of the lower
incomplete gamma function $\gamma$ and the gamma function $\Gamma$, and $\Phi$
is the cumulative distribution function of the standard normal distribution.
Solving Eq. (49) for $\eta$ gives a detection threshold for $J^{\prime}$ of
$\eta\approx 87$. If the $J^{\prime}$ value of the most likely planet position
is greater than $\eta$ it can be considered a detection.

Using $J^{\prime}$ as a detection criterion is, however, unsatisfactory as
explained above, and $J^{\prime\prime}$ \- with positivity constraint on the
flux - is better suited. At a single wavelength, and without any point sources
being present, $J^{\prime\prime}_{\lambda}$ follows a $\chi^{2}$-distribution
with one degree of freedom - the distribution of squares of normally
distributed values - with a threshold at 0 due to the positivity constraint on
the estimated flux values:

$\displaystyle
p(J_{\lambda}^{\prime\prime})=\dfrac{1}{2}\delta_{0}(J^{\prime\prime}_{\lambda})+\dfrac{1}{2}\chi^{2}_{1}(J^{\prime\prime}_{\lambda}).$
(50)

The probability distribution of the sum of $L$ independent values following
the above distribution is given by

$\displaystyle
p(J^{\prime\prime})=\dfrac{1}{2^{L}}\delta_{0}(J^{\prime\prime})+\sum_{\lambda=0}^{L-1}\dfrac{L!}{2^{L}\lambda!(L-\lambda)!}\chi^{2}_{L-\lambda}(J^{\prime\prime}).$
(51)

In Fig. 7 the bottom panel shows as a blue histogram the distribution of the
$J^{\prime\prime}$ values found in the map on the top-right panel. It follows
the expected distribution from Eq. (50).

Thus, similar to Eq. (49), a detection criterion $J^{\prime\prime}>\eta$ can
be defined by finding $\eta$ such that

$\displaystyle\dfrac{1}{2^{L}}+\sum_{\lambda=0}^{L-1}\dfrac{L!}{2^{L}\lambda!(L-\lambda)!}F_{L-\lambda}(\eta)=\Phi(5)$
(52)

with the left hand side being the cumulative distribution function of the
probability distribution defined in Eq. (51).

Solving Eq. (52) with $L=31$ for $\eta$ gives a detection threshold of
$\eta\approx 65$ for the cost function with positivity constraint. A value of
$\eta=65$ in the cost map $J^{\prime\prime}$ thus corresponds to a 5-$\sigma$
detection.

Figure 7: Top: Different versions of the cost function $J$ for a simulated
observation without any point sources. $L=31$ wavelength bins are used. Top
left: $J^{\prime}$ without positivity constraint on the flux. The image is
point-symmetric with respect to the image center. Top right:
$J^{\prime\prime}$ with a positivity constraint on the planet flux. Compared
to the left image, the strict point-symmetry is lifted and the mean value is
lower. Bottom: Probability distribution of the cost function maps shown above.
The grey histogram shows the $J^{\prime}$ values obtained from the panel on
the top-left. The solid red line shows the probability density function of a
$\chi^{2}$-distribution with $L=31$ degrees of freedom. The blue histogram
shows the $J^{\prime\prime}$ values obtained from the panel on the top-right.
The solid orange line shows the theoretically expected distribution of the
$J^{\prime\prime}$ values given in Eq. (51). The dashed vertical lines
indicate the detection criteria of $\eta=87$ and $\eta=65$ as derived by Eq.
(49) and Eq. (52) respectively.

The detection measure presented here does not include spectral regularization.
While a regularization term improves the contrast of the detection maps (see,
Fig. 18 in the Appendix), it adds an interdependence between the wavelength
bins. Thus, the estimated signal, as well as the noise terms in the individual
wavelength bins, are no longer mutually independent. Because the derivation of
the detection criterion was based on the assumption of uncorrelated noise, it
cannot be applied directly in the case of regularization. However, it was
already pointed out by Mugnier et al. (2006), and is also indicated by the
example presented in Fig. 18, that regularization is mostly needed for the
detection of planets with a low underlying $\mathrm{SNR}\lesssim 5$.

##  4 Signal analysis

###  4.1 Single planet: Earth-twin at 10 pc

We apply the signal extraction method to the case study example of an Earth-
twin planet orbiting a Sun-like star at 10 pc (see, Table 1). We compare the
SNR of the extracted signal to predictions based on photon statistics and
determine the position of the signal. We simulate the planet at
$100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ projected angular separation from
the star. We remind the reader that for an assumed integration time of
$55\text{\,}\mathrm{h}$ one obtains a predicted integrated
$\mathrm{SNR_{pred}}=9.7$. The top panel of Fig. 8 shows the noisy time series
resulting from the simulation. The maximum likelihood analysis is performed on
a grid out to $200\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ angular separation
from the central star, with a resolution of
$\sim$0.5\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$$. The middle panel of Fig. 8
shows the cost function map $J^{\prime\prime}$ without regularization
($\mu=0$) for a single simulated measurement. The relatively high SNR allows
for a clear detection of the planet in the upper right corner. The detection
criterion as derived in Sect. 3.2 is fulfilled and the planet location is
estimated correctly within the resolution provided by the grid. The estimated
extracted spectral flux of the located planet is shown in the bottom panel in
Fig. 8. The extracted spectrum agrees with the underlying true spectrum within
the expected uncertainties.

Estimates for the location and the spectrum have been extracted from the
simulated data, but we still need to quantify the SNR and compare it to the
expected value. Applying Eq. (48) results in $\mathrm{SNR_{est}}=10.45$,
compared to the predicted $\mathrm{SNR_{pred}}=9.7$. The discrepancy is due to
the statistical nature of the simulated data as mentioned previously. If the
simulation is repeated, the outcome is slightly different as would be the case
for a real measurement. To better estimate the performance of the signal
extraction process in general, and the SNR estimation in particular, the
simulated measurement is repeated $1000$ times with the same set-up. For each
run the planet position and spectral flux are estimated from the simulated
data.

Figure 8: Top: Noisy time series signal from simulated data for an Earth-like
planet around a Sun-like star at $10\text{\,}\mathrm{p}\mathrm{c}$ (in blue).
The underlying noiseless planet signal is indicated by the red curve. Middle:
Detection map for simulated data of the Earth-like planet. The photon-noise
based SNR is $9.7$ (cf. Sect. 2.4). The planet is correctly found in the upper
right of the map. Bottom: Spectral flux of the Earth-like planet as estimated
from the simulated data (red data points with $1\sigma$ error bars). The blue
line shows the true flux of the simulated planet and the blue shaded region
indicates $\pm 1\sigma$ based on the photon statistics.

The top panel in Fig. 9 shows the distribution of the estimated angular
separation of the planets. The average of the estimated angular separation is
$99.7\pm 1.5\,\mathrm{mas}$. The azimuthal position was extracted correctly
for all simulated planets within the resolution of the grid ($\pm 1^{\circ}$).
While the extracted spectra are all noisy and are qualitatively similar to the
example shown in the bottom panel of Fig. 8, the bottom panel of Fig. 9 shows
the mean and the standard deviation of the extracted flux values per
wavelength bin over the $N=1000$ simulations in comparison with the input
spectrum. Apart from the short wavelength range, where the low flux and high
noise levels result in large uncertainties, the extracted mean spectrum agrees
very well with the input data.

Figure 9: Top: Distribution of the estimated angular separation from the host
star of 1000 simulated Earth-twins. The underlying true value is
$\theta=100\,\mathrm{mas}$, the average estimated value is $99.7\pm
1.5\,\mathrm{mas}$. The spacing of the bars corresponds to the angular
resolution of the grid on which the analysis was performed. Bottom: Mean (red
dots) with standard deviation (red error bars) of the distribution of the
estimated flux values over 1000 simulations of the same artificial planet. The
blue line shows the true spectrum, with the blue shaded region enclosing $\pm
1\sigma$ based on the photon statistics. Figure 10: Cost function detection
map for the iterative detection of three simulated planets (indicated by the
white arrows). The simulation parameters are given in Sect. 4.2. The left
image shows the detection map calculated from the full simulated data.
Subtracting the estimated signal contribution of the source detected in the
first map and recalculating the detection map results in the second image, and
after repeating in the third one. Removing a detected source from the data
also removes the side lobes created by the source, thus allowing for a clear
iterative detection.

###  4.2 Multi-planet systems

As most exoplanetary systems are expected to contain multiple planets, it is
important that extraction algorithms can deal with multiple objects. The
modulation signal recorded with a nulling interferometer is a superposition of
the signals from different sources within the field-of-view. The most
straightforward approach to differentiate individual planets is to find their
positions and the estimated fluxes one after another. To do so, the most
likely planet position, in absolute terms, is accepted as a detection provided
it fulfills the detection criterion described above, and the spectral flux is
estimated for that planet. The estimated planet signal is subtracted from the
data and the process is repeated until no other significant signal remains in
the data.

We applied the signal extraction process outlined above to a simulated
3-planet system. All planets have a radius $R_{\mathrm{p}}=1\,R_{\oplus}$ and
are located around a Sun-like star at $10\text{\,}\mathrm{p}\mathrm{c}$. Their
semi-major axes $a_{\mathrm{p}}$ are set at 0.8, 1.0, and 1.2 AU respectively.
All planets are assumed to radiate as black bodies with an equilibrium
temperature of $276\text{\,}\mathrm{K}$ scaled by the inverse square-root of
the semi-major axis, $T_{\mathrm{p}}\propto 1/\sqrt{a_{\mathrm{p}}}$. The
system is seen face-on, such that all planets appear at their respective
maximum angular separation. The azimuthal position of the planets is chosen
randomly. The system is assumed to harbor an exozodi disk with $z=3$. For an
assumed integration time of $55\text{\,}\mathrm{h}$ the predicted SNRs of the
three planets are 15.8, 10.4, and 6.9, respectively.

Figure 10 shows the three detection maps for the iterative signal extraction
process applied to the simulated data. In the left panel, the highest values
for the cost function indicate the correct position of the planet closest to
the host star, which has the highest expected SNR. Besides the side lobes of
the first planet, the positions of the other planets, as well as their side
lobes, are also visible. However, because the cost values of these features
are similar, it is unclear which signal corresponds to a true source and which
not. The subtraction of the highest SNR signal removes also its side lobes and
the second planet can be identified by the strongest remaining peak. Removing
the second signal allows the third source to be detected. After the third
source, no further point source is found above the detection criterion.
Extracting the signal of the three planets leads to estimated SNR values of
16.0, 10.1, and 7.1, respectively. To test the robustness of the extraction
process with respect to the planet position we iterate the signal extraction
over distinct realisations of the noise. We find that the angular separation
and the position angle of the planet is retrieved accurately within the given
uncertainties. These uncertainties scale in accordance to the respective
photon noise SNRs, meaning that the innermost planet is extracted to the
highest precision.

Figure 11: Temperature and radius estimates for the three simulated planets
with $R_{\mathrm{p}}=1\,R_{\oplus}$ and varying temperatures. The true values
are indicated by crosses, the estimated values by dots with corresponding
$1\sigma$ confidence ellipses. The values are estimated by fitting a black
body curve to the estimated spectra of the point sources, whose locations are
inferred from the detection maps shown in Fig. 10.

We fitted the extracted spectra of the three planets with a black body
spectrum to estimate their radii and effective temperatures as this will also
be done as one of the first steps for any planet candidates detected with
_LIFE_. Figure 11 shows the estimated values along with $1\sigma$ confidence
contours derived from the uncertainties in the extracted flux values. The true
values are indicated by crosses. The uncertainty in the estimated temperatures
and radii increases with decreasing SNR, from about $\sigma_{T}\approx$ 20 K
and $\sigma_{R}\approx$ 0.15 R⊕ for the brightest planet to
$\sigma_{T}\approx$ 35 K and $\sigma_{R}\approx$ 0.4 R⊕ for the faintest. As
the spectral shape of an emitting black body depends more strongly on the
temperature than on its size, the relative uncertainty in the derived
temperature is much smaller than the relative uncertainty in the radius.

###  4.3 Rocky, habitable zone exoplanets from _LIFE_ search phase

In the previous sections we investigated how well the signal of individual
exoplanets can be extracted and their basic properties, such as effective
temperature, radius and position, can be estimated. We now turn to a larger
sample of objects covering a broader range of properties and different noise
characteristics (e.g., spectral type of host star, exozodi level). The
extraction algorithm is applied to a subset of the simulated exoplanet
population that was used for the detection yield estimation in Paper 1 (Quanz
et al., 2021). Specifically, we focus on the rocky, habitable zone
planets666Defined as planets with radii $R_{\textrm{p}}$ in the range
$0.5\;R_{\oplus}\leq R_{\textrm{p}}\leq 1.5\;R_{\oplus}$ orbiting within the
spectral-type dependent empirical habitable zone (eHZ) of their host star. For
a Sun-like star the eHZ corresponds to an insolation range of
$0.32\;S_{\oplus}\leq S_{\textrm{p}}\leq 1.76\;S_{\oplus}$, where $S_{\oplus}$
corresponds to the present-day insolation level of the Earth. that were
detectable in a 2.5-year search phase with four $D=2$ m apertures and
investigate in detail how well their properties can be determined from single
epoch observations, assuming the planets radiate as black bodies at their
respective equilibrium temperature. From the 500 Monte Carlo runs that we
carried out in Quanz et al. (2021) in order to estimate the expected detection
yield of the _LIFE_ mission, we randomly select 100 runs. These runs contain
approximately 4400 rocky, HZ planets with predicted $\mathrm{SNR}\geq 7$,
which was considered as the detection threshold in Paper 1. Figure 12 shows
the distribution of the predicted SNR of these planets (based on photon-noise
statistics), their temperature, radius, and their projected angular separation
from their host star. The median of the SNR distribution is $\approx$13.3,
thus nearly twice the detection threshold. While, on average, the planet radii
are $>$1.0$\;R_{\oplus}$, the mean temperature of the detected planets is only
$T_{\mathrm{p}}\,$\approx 235\text{\,}\mathrm{K}$$, i.e., less than Earth’s
effective temperature.

Our signal extraction method is applied to each planet individually and the
detection maps are calculated without regularization of the estimated spectra.
The detection criterion as derived in Sect. 3.2 is applied. We find that
$\approx$1.5% of the simulated planets do not fulfill the detection criterion.
For a few additional planets the position was extracted incorrectly. Overall,
about 98% of the planets were detected at a separation deviating $\leq$15%
from the true value and with an error in position angle $\leq$10∘.

For the detected planets, the temperature and radius are estimated by fitting
a black body spectrum to the extracted flux. Additionally, the SNR is
estimated from the extracted flux as described in Sect. 4.1. Because each of
the estimated parameters has a broad distribution among the simulated planets,
the relative deviation of the estimated values from the underlying true values
is analyzed. The results of the parameter extraction are presented in Fig. 13.
The corner plot shows again the distribution of the predicted SNR and the
ratios between the estimated and the true values of the planetary parameters
(i.e., temperature, radius, and separation) as well as the ratio between the
estimated and the purely photon statistic based SNRs. From the first column it
is evident that the estimation of all parameters improves with increasing SNR
of the planet, as the statistical distribution of the estimates becomes
narrower. Overall, the SNR is slightly underestimated, even for relatively
high SNR values $>25$. This can be explained by an overestimated noise
variance from the simulated data, as it is calculated as the variance of the
full data including possible planet signals.

The deviation of the estimates of the temperature and the radius from the true
values are approximately symmetric in terms of over- or underestimation.
Across the full sample, the temperature is estimated to
$T_{\mathrm{est}}/T_{\mathrm{true}}=1.01^{+0.08}_{-0.06}$ and the radius to
$R_{\mathrm{est}}/R_{\mathrm{true}}=0.97\pm 0.18$. As expected, the estimate
of the temperature and the radius are highly correlated as both parameters are
positively correlated with the total emitted radiation of the planet. The
separation is estimated to
$\theta_{\mathrm{est}}/\theta_{\mathrm{true}}=0.99\pm{0.01}$.

Figure 12: Distribution of properties of detectable rocky, habitable zone
planets from the detection yield estimates presented in Quanz et al. (2021)
normalized over 100 Monte Carlo runs. From top to bottom: predicted SNR (based
on photon-noise), effective temperature, radius, and projected angular
separation. Figure 13: Corner plot showing the predicted SNR based on photon
statistics (cf. top panel in Fig. 12), and the estimated values for SNR,
temperature, radius, and angular separation relative to the true values based
on the signal analysis of $\approx 4400$ rocky, habitable zone planets. We
note that the discretization of the values for the estimated angular
separation corresponds to the resolution of the spatial grid on which the
analysis was performed. At the top of each column the mean and the standard
deviation of the distributions are given.

##  5 Discussion and Conclusions

We investigate the signal extraction process for simulated interferometric
measurements to study the detectability of potential exoplanets and the
ability to characterize them. We implement a maximum-likelihood algorithm
previously proposed for the _Darwin_ mission and apply it to simulated data.
This yields the three important results that (a) this method can be
successfully applied to multi-planet systems, (b) the signal extraction method
is able to reproduce SNR estimates that are in agreement with expectations
from photon-noise statistics and (c) planetary properties inferred by the
signal extraction are of high accuracy. These results are especially
interesting in the context of Quanz et al. (2021), since this publication
predicts the exoplanet yield of the LIFE mission based on a synthetic planet
population without considering the impact of signal extraction. Our results
indicate that an inclusion of the signal extraction in these simulations would
not fundamentally alter the results.

Yet, under certain conditions the signal extraction is especially demanding.
In the low-SNR regime, we have qualitatively shown how spectral regularization
can help to disentangle the planet signal from the background noise. However,
as spectral regularization relies on a priori assumptions about the planet
spectrum, it is not free of bias and adds a degree of freedom to the analysis.
To apply regularization to a large set of simulated planets, a detection
criterion depending on the amount of enforced spectral regularization would
have to be defined. Additionally it has to be investigated how the
regularization affects the estimated planet spectrum quantitatively.

To further assess potential limitations when dealing with multi-planet
systems, studies with more diverse planet properties or very small projected
angular separations could be envisaged. Potential improvements when dealing
with several signals were already presented by Thiébaut & Mugnier (2005).
These authors discussed an extrapolation of the maximum-likelihood method to
estimate the planet positions and spectra for multiple planets at once, which
avoids error propagation from the estimation of the first detected planet to
the other ones. This generalized approach is computationally more expensive,
but a simplified option estimates the planet positions one after the other but
recomputes the optimal flux of all detected planets at each iteration.

Furthermore, the ability to correctly extract multi-planet systems is also
governed by the effective angular resolution of the array. Therefore, a study
of diverse multi-planet extractions would be able to provide minimum
requirements for the angular resolution. Since the angular resolution has been
shown to be a dominant factor in the design of nulling interferometers (Lay,
2005), such a study will become vital in the ongoing development of _LIFE_.

Going forward, one could easily imagine that new approaches involving modern
machine-learning based signal extraction algorithms should be tested and their
performance compared. While already being explored in other fields (e.g.,
Cuoco et al., 2020; Gebhard et al., 2019, 2020), nulling interferometry was
not yet a focus point of these efforts. This will be particularly important as
soon as systematic noise terms can be simulated in a realistic way.

An additional factor degrading the quality of the measurement can be the
occurrence of instrumental noise. Within the _LIFE_ project, the current
working assumption is that the instrumental noise can be neglected based on
its contribution being significantly smaller than the fundamental noise. We
show, that this assumption places strict requirements on the instrumental
stability and prompts continued technical developments, for example in the
fields of low dark-current detectors or highly precise OPD control techniques.

Similarly, signal extraction has been shown to play an integral role in the
prediction of instrumental noise (Lay, 2004). Certain steps in the signal
retrieval, for example the explicit fitting and removal of systematic noise,
can significantly relax technical requirements (Lay, 2006). As LIFEsim
progresses to include the simulation of instrumental noise, the removal of the
same should be considered with the methods presented here.

In future work, we aim to shift from the fundamental noise limited regime into
a regime in which some instrumental noise is retained in the signal. A
prerequisite for this is that the instrumental noise is well understood,
tracked and considered in the interpretation of the measurements. Defrère et
al. (2010) offers a list of techniques targeted at facilitating this treatment
of instrumental noise. The goal of this future work will be to achieve the
scientific output postulated by the _LIFE_ mission with technology closer to
contemporary and near-future instrumentation.

In summary, we present a publicly available777www.life-space-mission.com tool
for simulating the exoplanet search phase of the _LIFE_ mission called
LIFEsim888Because the current version of LIFEsim has already been applied to a
number of science case simulations, we would like to enable the community to
verify and reproduce the results, and also to investigate additional science
cases.. Using the presented planetary signal extraction methods, we show that
the synthetic exoplanet population used for the mission yield evaluation in
Quanz et al. (2021) can be retrieved. We demonstrate the precise extraction of
planetary radius and temperature, which is of special interest since they are
fundamentally important parameters for characterizing the (atmospheric)
properties of exoplanets and for categorizing and prioritizing objects for
potential in-depth follow-up investigations. We remind the reader that these
two parameters are much more difficult to derive from reflected light
measurements as there is a degeneracy between planet size and geometric albedo
and no immediate information about the effective temperature can be obtained.
These results illustrate the high quality of information one can expect from
single-epoch observations with _LIFE_ and how they will complement the results
from future reflected light missions.

###### Acknowledgements.

We thank the anonymous referee for a critical and constructive review of the
original manuscript which helped improve the quality of the paper
significantly. This work has been carried out within the framework of the
National Centre of Competence in Research PlanetS supported by the Swiss
National Science Foundation. SPQ acknowledges the financial support of the
SNSF. TL was supported by the Simons Foundation (SCOL award No. 611576). This
research has made use of the following Python packages: astropy (Astropy
Collaboration et al., 2013, 2018), matplotlib (Hunter, 2007), numpy (Van Der
Walt et al., 2011), scipy (Virtanen et al., 2020).

_Author contributions:_ MO and FD contributed equally to this paper. MO wrote
the original LIFEsim tool, carried out the main analyses and wrote part of the
manuscript. FD contributed to the main analyses and the LIFEsim tool,
programmed the publicly available github version, led the instrumental noise
analyses and wrote part of the manuscript. SPQ initiated and guided this
project and wrote part of the manuscript. RL, EF and AGh contributed to the
LIFEsim tool. All authors discussed the results and commented on the
manuscript.

## References

  * Angel & Woolf (1997) Angel, J. R. P. & Woolf, N. J. 1997, The Astrophysical Journal, 475, 373
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Benz et al. (2021) Benz, W., Broeg, C., Fortier, A., et al. 2021, Experimental Astronomy, 51, 109
  * Bracewell (1978) Bracewell, R. N. 1978, Nature, 274, 780
  * Cabrera et al. (2019) Cabrera, M. S., McMurtry, C. W., Forrest, W. J., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 1
  * Cockell et al. (2009) Cockell, C. S., Herbst, T., Léger, A., et al. 2009, Experimental Astronomy, 23, 435
  * Cuoco et al. (2020) Cuoco, E., Powell, J., Cavaglià, M., et al. 2020, Machine Learning: Science and Technology, 2, 011002
  * Dandumont et al. (2020) Dandumont, C., Defrère, D., Kammerer, J., et al. 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 035004
  * Defrere (2009) Defrere, D. 2009, PhD thesis, University of Liege
  * Defrère et al. (2010) Defrère, D., Absil, O., Den Hartog, R., Hanot, C., & Stark, C. 2010, Astronomy and Astrophysics, 509, 1
  * Defrère et al. (2012) Defrère, D., Stark, C., Cahoy, K., & Beerer, I. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84420M
  * Den Hartog & Karlsson (2005) Den Hartog, R. & Karlsson, A. 2005, The DARWINsim science simulator, Tech. rep., ESA
  * Draper et al. (2006) Draper, D. W., Elias II, N. M., Noecker, M. C., et al. 2006, The Astronomical Journal, 131, 1822
  * Ertel et al. (2020) Ertel, S., Defrère, D., Hinz, P., et al. 2020, The Astronomical Journal, 159, 177
  * Flasseur et al. (2020) Flasseur, O., Denis, L., Thiébaut, E., & Langlois, M. 2020, Astronomy & Astrophysics, 9
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, Astronomy & Astrophysics, 616, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, Astronomy & Astrophysics, 595, A1
  * Gáspár et al. (2020) Gáspár, A., Rieke, G. H., Guillard, P., et al. 2020, Publications of the Astronomical Society of the Pacific, 133, 014504
  * Gaudi et al. (2020) Gaudi, B. S., Seager, S., Mennesson, B., et al. 2020, arXiv e-prints, arXiv:2001.06683
  * Gebhard et al. (2020) Gebhard, T. D., Bonse, M. J., Quanz, S. P., & Schölkopf, B. 2020, arXiv e-prints, arXiv:2010.05591
  * Gebhard et al. (2019) Gebhard, T. D., Kilbertus, N., Harry, I., & Schölkopf, B. 2019, Physical Review D, 100
  * Gheorghe et al. (2020) Gheorghe, A., Glauser, A. M., Ergenzinger, K. J., et al. 2020, in Optical and Infrared Interferometry and Imaging VII, ed. A. Mérand, S. Sallum, & P. G. Tuthill (Online Only, United States: SPIE), 111
  * Glasse et al. (2015) Glasse, A., Rieke, G. H., Bauwens, E., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 686
  * Guyon et al. (2013) Guyon, O., Mennesson, B., Serabyn, E., & Martin, S. 2013, Publications of the Astronomical Society of the Pacific, 125, 951
  * Hansen et al. (2022) Hansen, J. T., Ireland, M. J., & the LIFE Collaboration. 2022, Large Interferometer For Exoplanets (LIFE): VI. Ideal kernel-nulling array architectures for a space-based mid-infrared nulling interferometer
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Kelsall et al. (1998) Kelsall, T., Weiland, J., Franz, B., et al. 1998, The Astrophysical Journal, 508, 44
  * Kennedy et al. (2015) Kennedy, G. M., Wyatt, M. C., Bailey, V., et al. 2015, Astrophysical Journal, Supplement Series, 216
  * Konrad et al. (2021) Konrad, B. S., Alei, E., Angerhausen, D., et al. 2021, Large Interferometer For Exoplanets (LIFE): III. Spectral resolution, wavelength range and sensitivity requirements based on atmospheric retrieval analyses of an exo-Earth
  * Laugier et al. (2020) Laugier, R., Cvetojevic, N., & Martinache, F. 2020, arXiv e-prints, arXiv:2008.07920
  * Lawson et al. (2007) Lawson, P. R., Lay, O. P., Martin, S. R., et al. 2007, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6693, Techniques and Instrumentation for Detection of Exoplanets III, ed. D. R. Coulter, 669308
  * Lay (2004) Lay, O. P. 2004, Applied Optics, 43, 6100
  * Lay (2005) Lay, O. P. 2005, Applied Optics, 44, 5859
  * Lay (2006) Lay, O. P. 2006, in SPIE Astronomical Telescopes + Instrumentation, ed. J. D. Monnier, M. Schöller, & W. C. Danchi (SPIE), 62681A–14
  * Lay et al. (2007) Lay, O. P., Martin, S. R., & Hunyadi, S. L. 2007, Techniques and Instrumentation for Detection of Exoplanets III, 6693, 66930A
  * Mayor & Queloz (1995) Mayor, M. & Queloz, D. 1995, Nature, 378, 355
  * Mugnier et al. (2006) Mugnier, L., Thiébaut, E., Belu, A., et al. 2006, EAS Publications Series, 22, 69
  * Mugnier et al. (2009) Mugnier, L. M., Cornia, A., Sauvage, J.-F., et al. 2009, Journal of the Optical Society of America A, 26, 1326
  * Nagler et al. (2021) Nagler, P. C., Sadleir, J. E., & Wollack, E. J. 2021, Journal of Astronomical Telescopes, Instruments, and Systems, 7, 1
  * Paresce & Crane (1997) Paresce, F. & Crane, P., eds. 1997, Science with the VLT Interferometer: Proceedings of the ESO Workshop Held at Garching, Germany, 18–21 June 1996, ESO Astrophysics Symposia (Berlin, Heidelberg: Springer Berlin Heidelberg)
  * Perido et al. (2020) Perido, J., Glenn, J., Day, P., et al. 2020, Journal of Low Temperature Physics, 199, 696
  * Quanz et al. (2021) Quanz, S. P., Ottiger, M., Fontanet, E., et al. 2021, arXiv e-prints, arXiv:2101.07500
  * Rieke et al. (2015) Rieke, G. H., Ressler, M. E., Morrison, J. E., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 665
  * Roellig et al. (2020) Roellig, T. L., McMurtry, C., Greene, T., et al. 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6
  * The LUVOIR Team (2019) The LUVOIR Team. 2019, arXiv e-prints, arXiv:1912.06219
  * Thiébaut & Mugnier (2005) Thiébaut, E. & Mugnier, L. 2005, Proceedings of the International Astronomical Union, 1, 547
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, 46, 135
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

##  Appendix A Fundamental noise limit

Section 2.3.2 states that for the analysis presented in this paper a
measurement dominated by photon noise originating from the astrophysical
sources is assumed. In the following, this state will be referred to as the
fundamental noise limited regime. Considering the technical complexity of
nulling interferometers, it is important to substantiate this assumption and
offer proof that it is well founded. The aim of this appendix is the
implementation of an instrumental noise model and the subsequent prediction of
the maximum levels of perturbations to the instrument that allow the system to
stay in the fundamental noise dominated regime.

###  A.1 Assumptions

The term fundamental noise limited is readily used in literature as well as in
the _LIFE_ paper series. Under the goal of producing quantitative results,
this must be translated into a more tangible constraint. The perfect
instrument will only be subject to the photon noise by the astrophysical noise
sources $\sigma_{\mathrm{fund}}$999Square-root of the number of photons
contributing to photon noise per unit time. Perturbing the instrument will
give rise to two additional types of noise terms. First and foremost, a
systematic perturbation (e.g., a variation in the amplitude response of one of
the interferometric arms) of the instrument will introduce a systematic noise
term $\sigma_{s,\mathrm{inst}}$ arising from an additional photon rate
introduced by the perturbations. These additional photons can resemble a
signal as would be produced by a target exoplanet. Second, additional photon
noise sources like a detector dark current or thermal emission within the
instrument can be considered. Moreover, the systematic perturbations of the
central destructive fringe will average to a mean reduction in the effective
null-depth. This will increase the overall amount of photons received from the
target star and subsequently also increase the photon noise. We collect all
additional photon noise sources into the instrumental photon noise term
$\sigma_{p,\mathrm{inst}}$ and define the instrumental noise term

$\displaystyle\sigma_{\mathrm{inst}}=\sqrt{\sigma_{s,\mathrm{inst}}^{2}+\sigma_{p,\mathrm{inst}}^{2}}.$

The fundamental noise limited case is then interpreted as the measurement
being dominated by the photon noise arising in the perfect instrument over the
instrumental noise. Therefore, we define the fundamental noise limited case as

$\displaystyle\sigma_{\mathrm{fund}}\geq\sigma_{\mathrm{inst}},$ (53)

where the fundamental noise contribution is larger that the instrumental noise
contribution. Semantically, it could also be argued that the fundamental noise
being a multiple of the instrumental noise could be required. Therefore, we
also present the case of a five-times larger fundamental noise in the
following.

For the array and instrument configuration, we mirror the setup described in
Section 2.1 and Table 1. This includes the simulation of phase-chopping, which
is assumed to be facilitated by subtracting two simultaneous realisations of
phase-inverted outputs. This simplifies the calculation of instrumental noise
terms as, contrary to Lay (2004), we do not have to include errors at the
chopping frequency. All astrophysical sources are simulated according to the
models described in Section 2.2. If not stated otherwise, this appendix
analyses an observation of an Earth-twin around a Sun-like star located at 10
pc distance with a $z=3$ exozodi level (cf. Table 1).

In the following simulations, we assume a non-perfect instrument affected by
perturbations. Systematic perturbations are assumed to affect the amplitude
response $A$, the phase response $\phi$ and the polarization rotation $\theta$
of each of the four interferometric arms as well as the position $x$ and $y$
of the collector spacecraft. These dynamic perturbations follow a perturbation
spectrum centered around values for the configuration of the ideal
instrument101010We note the lack of DC-components in the perturbation spectra.
In this work, we omit any constant static offsets from the ideal instrument
configuration (e.g. a constant amplitude response offset in one of the arms).
This omission is based on static offsets being more easily calibrated and
removed compared to dynamic perturbations. Due to a lack of reliable empirical
data, we stay agnostic to how these perturbations are generated and assume the
general noise profiles listed in Table 11. Apart from the aforementioned
general systematic perturbations, we examine two more sources of perturbation
in an effort to illustrate the allowed magnitude of the remaining perturbation
terms. We consider the detector dark current $I_{D}$ and the thermal
background seen by a detector in an environment of temperature $T$. We assume
a detector with physical pixel dimensions of the JWST/MIRI detectors
($25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ pixel pitch, see Rieke et al.
2015) and a spectral sampling of 2.2 pixels per wavelength channel.

Table 2: Parameters of the instrument perturbations at a reference wavelength
of $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

Perturbation | Shape | Cutoff | rms | Wavelength Dependence
---|---|---|---|---
Amplitude $\delta A$ | pink noise, no DC | 10 kHz | 0.1 % | $\lambda^{-1.5}$
Phase $\delta\phi$ | pink noise, no DC | 10 kHz | 0.001 rad | $\lambda^{-1}$
Polarization $\delta\theta$ | pink noise, no DC | 10 kHz | 0.001 radaa$a$Decreased by an order of magnitude due to a typing error in the original publication (Dèfrere 2021, private communication) | none
Collector Position $\delta x$, $\delta y$ | white noise | 0.64 mHz | 1 cm | none

111111Values taken from Lay (2004).

###  A.2 Methods

The simulations are largely based on the approach presented in Lay (2004).
This paper was published in the context of the TPF-I mission and offers a
comprehensive framework for estimating the instrumental noise for nulling
interferometers. The following paragraphs will sparsely trace the calculations
performed in Lay (2004) to highlight fundamental properties.

First, the overall interferometric response is decomposed into a sum of
responses of the individual baselines

$\displaystyle
N=\sum_{j}\sum_{k}A_{j}A_{k}\left[\cos(\phi_{j}-\phi_{k})\overline{B}_{\mathrm{sym},jk}-i\sin(\phi_{j}-\phi_{k})\overline{B}_{\mathrm{asym},jk}\right].$
(59)

Here, $N$ is the detected photon rate, $A_{n}$ and $\phi_{n}$ are the
amplitude and phase responses of the $n$ interferometric arms and
$\overline{B}_{\mathrm{(a)sym},nm}$ is the Fourier transform of the (anti-)
symmetric sky brightness distribution under the baseline from the
$n^{\mathrm{th}}$ to the $m^{\mathrm{th}}$ collector spacecraft. Equation (59)
reflects an intrinsic symmetry property of nulling interferometers: baselines
with $\pi$-multiple phase difference respond to symmetric (with respect to the
line of sight) source while baselines with $\pi/2$-multiple phase differences
respond to anti-symmetric sources. The sensitivity of the photon rate against
perturbations of the instrument (here in phase $\delta\phi_{n}$ and relative
amplitude $\delta a_{n}=\delta A_{n}/A_{n}$) is captured in a second order
Taylor expansion

$\displaystyle\delta N\approx$ $\displaystyle\sum_{j}\left[\frac{\partial
N}{\partial A_{j}}A_{j}\delta a_{j}+\frac{\partial
N}{\partial\phi_{j}}\delta\phi_{j}\right]$
$\displaystyle+\sum_{j}\sum_{k}\left[\frac{1}{2}\frac{\partial^{2}N}{\partial
A_{j}\partial A_{k}}A_{j}A_{k}\delta a_{j}\delta a_{k}\right.$
$\displaystyle\left.+\frac{\partial^{2}N}{\partial
A_{j}\partial\phi_{k}}A_{j}\delta
a_{j}\delta\phi_{k}+\frac{1}{2}\frac{\partial^{2}N}{\partial\phi_{j}\partial\phi_{k}}\delta\phi_{j}\delta\phi_{k}\right].$
(60)

The partial derivatives in this equation represent the sensitivity that the
photon rate exhibits against a perturbation of respective kind and order. We
note the appearance of second order terms as well as amplitude-phase cross
terms. At this point in the calculation, the symmetry of the sources is used
to simplify the expressions and arrive at the sensitivity to the
perturbations. This will be demonstrated using the sensitivity to stellar
leakage. In first approximation, the brightness distribution of the stellar
disk is centrally symmetric, leading to $\overline{B}_{\mathrm{*,asym},nm}=0$.
This reduces Equation (59) to

$\displaystyle
N_{*}=\sum_{j}\sum_{k}A_{j}A_{k}\cos(\phi_{j}-\phi_{k})\overline{B}_{\mathrm{*},jk}.$

Now the sensitivity of the photon rate towards, e.g., first order amplitude
perturbations is given by

$\displaystyle C_{A_{j}}^{*}=A_{j}\frac{\partial N_{*}}{\partial
A_{j}}=2A_{j}\sum_{k}A_{k}\cos(\phi_{j}-\phi_{k})\overline{B}_{\mathrm{*},jk}.$

$C_{A_{j}}^{*}$ is called the first order stellar leakage sensitivity
coefficient and all remaining sensitivity coefficients in Equation (60) can be
determined in the same fashion. In a final step, the photon rates $N$ and
$\delta N$ are cross-correlated with planetary template functions, which
represents the signal extraction process. The photon noise of the instrument
then corresponds to the square-root of the mean photon rate and is therefore
given by

$\displaystyle\sigma_{\mathrm{p}}=\sqrt{N+\langle\delta N\rangle},$

where $\langle\rangle$ denotes the ensemble average over all perturbed states
of the system. In this description, we identify the following correspondences
to the variables defined in the previous Section A.1:
$N=\sigma_{\mathrm{fund}}^{2}t_{\mathrm{int}}$ and $\langle\delta
N\rangle=\sigma_{\mathrm{p,inst}}^{2}t_{\mathrm{int}}$ where
$t_{\mathrm{int}}$ is the integration time. The systematic noise is accessed
via

$\displaystyle\sigma_{\mathrm{s,inst}}=\sqrt{\langle\delta N^{2}\rangle}.$

Up to second order, the method of phase-chopping is able to remove all
systematic noise contributions except for one first order phase
($\delta\phi_{n}$) and one amplitude-phase cross term ($\delta
a_{n}\delta\phi_{m}$). We remind the reader that the order of the term is not
connected to it’s impact but refers to the proportionality between the
perturbation and the photon rate in Equation (60). Identifying the
perturbation with the noise spectra given in Table 11 enables the calculation
of the required noise sources. We refer the reader to Lay (2004) for the full
details of the methods described in this section.

###  A.3 Results

Figure 14: Wavelength dependence of noise contributions for observation of an
Earth-twin located at 10 pc with _LIFE_ in the baseline scenario. Panel a):
Division of the fundamental noise into its individual sources. Panel b):
Instrumental noise sources split into systematic noise sources (green) and
additional photon noise sources arising from instrumental effects. Panel c):
Planet signal, fundamental noise and instrumental noise contributions. A SNR
using only fundamental noise (orange bars) is representing the fundamental
noise limited case used in the _LIFE_ paper series. The green bars show the
full SNR additionally considering the instrumental noise as well as the
fundamental noise.

The signal and noise contributions using the assumptions laid out above for
every wavelength channel defined in the _LIFE_ baseline scenario are shown in
Figure 14. For the comparison of the fundamental noise sources in Panel a), we
reconfirm the behaviour described in Section 2.4: The short-wavelength end of
the spectrum is dominated by the stellar leakage term while the long-
wavelength end of the spectrum is dominated by the local-zodiacal dust
emission, which can not be suppressed by nulling due to its incoherent nature.
For the first and second order systematic noise terms shown in Panel b) we
stress again that the use of phase chopping removes all systematic noise terms
except for the first order phase deviations and the second order amplitude-
phase deviation cross terms. This removal helps to relax the technical
requirements imposed by the systematic noise. We reproduce the finding of Lay
(2004) that the second order term is generally larger than the first order
term. Panel b) also shows that the instrumental noise is either white or
increases towards shorter wavelengths. This is driven by the perturbation
scaling with wavelength described in Table 11 and by the increased amount of
stellar light emitted in this regime. Panel c) of Figure 14 presents the
impact of the additional noise sources on the SNR. Here, we assume a detector
environment temperature of $T=$11\text{\,}\mathrm{K}$$ and a dark current of
$I_{\mathrm{D}}=10^{-4}\,\,\mathrm{e^{-}\,px^{-1}\,s^{-1}}$. As expected, the
largest relative SNR reduction is induced towards the shorter wavelength end
of the spectrum where the systematic noise is most dominant. In the
fundamental noise limit, the condition in Equation (53) must be valid for all
wavelength bins. The relative trends of instrumental and fundamental noise in
Figure 14 indicate that this condition should be tested towards the shorter
wavelength end of the spectrum. We therefore select three discrete wavelength
bins below $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for evaluation121212Bins
at $4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$6.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$10.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$, with widths
$\Delta\lambda_{b}=$0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\Delta\lambda_{b}=$0.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$\Delta\lambda_{b}=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. For
completeness, we list the numerical values for the relative signal and noise
contribution in this shortest wavelength bin in Table 13.

Table 3: Values for the signal and noise contributions for the _LIFE_ baseline
case with perturbations as defined in Table 11 for the shortest wavelength bin
of $\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.

Source | Photon Rate ($\mathrm{s^{-1}}$)
---|---
Planet Signal | $1.1\cdot 10^{-4}$
Fundamental Noise | $5.8\cdot 10^{-2}$
Stellar geometric leakage | $5.7\cdot 10^{-2}$
Local zodi leakage | $5.0\cdot 10^{-4}$
Exo zodi leakage | $9.1\cdot 10^{-3}$
Instrumental systematic noise | $3.0\cdot 10^{-1}$
First order phase | $6.0\cdot 10^{-2}$
Second order amplitude-phase | $2.9\cdot 10^{-1}$
Instrumental photon noise | $1.8\cdot 10^{-2}$
Stellar null-floor leakage | $1.0\cdot 10^{-2}$
Polarisation angle | $2.2\cdot 10^{-3}$
Detector thermal background | $1.1\cdot 10^{-3}$
Detector dark current | $1.5\cdot 10^{-2}$
SNR, one rotation (56 h)aa$a$Calculated by integrating over all wavelength bins for the specified integration time and can be compared to the value computed in Figure 5 (SNR$\approx$9.7), where the same observing example was used, but without including instrumental noise terms. | $4.5$
131313

To infer the acceptable systematic noise levels, we run the simulations over a
grid of pink noise amplitude and phase perturbations with varying root-mean-
squared (rms) values141414We remind the reader that the use of phase chopping
removes most of the coupling to other sources of perturbation.. The rms is
specified at a reference wavelength of $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ and scaled according to Table 11. Phase errors are converted to optical
path difference (OPD) errors to allow for a better comparison with literature.
The ratio of instrumental to fundamental noise on this grid is shown in Figure
15.

Figure 15: The ratio of fundamental noise to instrumental noise on a grid of
amplitude and phase perturbations for the shortest-wavelength bin at
$\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ with width
$\Delta\lambda_{b}=$0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The
perturbations are assumed to be pink noise with the specified reference rms at
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and the phase perturbation is
converted to an optical path difference (OPD). The contours indicate where the
fundamental noise dominates the systematic noise by a factor of one, five or
ten respectively for the $\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$ case (in red) and for comparison for the
$\lambda_{b}=$6.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (in orange) and
$\lambda_{b}=$10.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (in cyan) cases.
The crosses indicate the reference levels specified by Table 11.

Firstly, the figure confirms that certain regions allow for a trade-off
between the OPD error rms and the amplitude error rms, where increasing one
allows for a reduction of the other. Outside of these regions, there exist
hard upper limits for the maximally allowed rms error for both amplitude and
OPD. It will be vital for the design of the instrument to stay clear of these
limits, since they imply that, if reached in one perturbation rms, any
improvement to the other perturbation rms does not improve instrument
performance. Additionally, we can reconfirm that the requirements are much
less strict in the longer wavelength regime. This grid can serve as a rough
predictor for the maximally acceptable noise levels: acceptable and too high
levels of noise are divided by the contour tracing a ratio of one between the
fundamental and systematic noise. Within the presented framework, this
produces sets of allowable noise levels shown in Table 15. These values can
only be seen as rough predictors, since we are not presenting a complete
instrument model and additional noise sources not considered in the present
analysis could further increase the noise requirements.

Table 4: Approximate predictions for the maximum allowed perturbation rms to
stay in the fundamental noise dominated case.

$\sigma_{\mathrm{fund}}=x\sigma_{\mathrm{inst}}$ | $x=1$ | $x=5$
---|---|---
OPD error rms | Relative amplitude error rms
(nm) | (in %)
0.75 | 0 | …
0.50 | 0.02 | …
0.16 | 0.08 | 0
0.10 | 0.13 | 0.02
0.05 | 0.24 | 0.05
0.01 | 0.62 | 0.12
0 | 0.67 | 0.14
151515rms perturbations values are given relative to the reference wavelength
of $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Presented values are for an
observation of an Earth-twin at 10 pc with LIFE configured in the baseline
scenario in the $\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$
wavelength bin, i.e. the given values follow the contours in Figure 15. The
asymptotic upper limits are indicated by the bold text.

A similar analysis can be done for the instrumental photon noise sources. As
representative terms for thermal background emission and white shot noise term
we analyse the thermal emission as seen by the detector and a dark current
term for the detector.

Figure 16: The ratio of fundamental noise to instrumental noise on a grid of
the detector environment temperature and the detector dark current for the
shortest-wavelength bin at
$\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ with width
$\Delta\lambda_{b}=$0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The
temperature induces black body thermal emission which is propagated directly
into the detector. The dark current is assumed to be constant with
temperature. The instrument is configured according to the reference state in
Table 11. The contours indicate where the fundamental noise dominates the
systematic noise by a factor of one or five. The cross indicates the reference
levels of $T=$11\text{\,}\mathrm{K}$$ and
$I_{\mathrm{D}}=10^{-4}\,\,\mathrm{e^{-}\,px^{-1}\,s^{-1}}$ used in Figure 14.

Figure 16 shows an expected behaviour. Since both noise sources are
independent of each other, there exists an absolute upper limit for the
acceptable level of both terms (see Table 16).

We note that while meeting the requirement for the thermal background
temperature of the detector does not appear too challenging, the requirements
on the dark current seem to be at least one order of magnitude more stringent
than what has been achieved so far for both Si:As IBC detector arrays and mid-
infrared MCT devices (e.g., Rieke et al. 2015; Cabrera et al. 2019; Roellig et
al. 2020; Gáspár et al. 2020). The challenges related to MIR detector
technology (and availability) has already been discussed in LIFE paper I
(Quanz et al. 2021). An interesting alternative approach could be kinetic
inductance detectors (KIDs) or Transition Edge Sensor detectors (TES) that may
also be applicable for the LIFE wavelength range (e.g., Perido et al. 2020;
Nagler et al. 2021).

Table 5: Approximate predictions for the asymptotic upper limits in detector
thermal background and detector dark current perturbation to stay in the
fundamental noise dominated case.

$n_{\mathrm{p}}=xn_{\mathrm{s}}$ | $x=1$ | $x=5$
---|---|---
Thermal background temperature (K) | 15.1 | 13.4
Detector dark current ($\mathrm{e^{-}\,s^{-1}\,px^{-1}}$) | $4.6\cdot 10^{-3}$ | $7.5\cdot 10^{-5}$
161616Presented values are for an observation of an Earth-twin at 10 pc with
LIFE configured in the baseline scenario in the
$\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ wavelength bin, i.e.
the given values correspond to the limits in Figure 16.

The acceptable perturbation levels will depend on the configuration of the
observed system. First and foremost, the properties of the star and exo-
zodical disk will influence both, the amount of fundamental noise received and
how the amount of systematic noise scales with a given level of perturbation.
Apart from this there are additional, more subtle effects that change the
noise levels. Prominent examples are the scaling of the array baseline with
the host star luminosity (see Section 2.3.1), relative shifts in main emission
wavelength between the astrophysical objects and scaling of absolute flux with
the target distance. A full study of all the aforementioned effects is vital,
but out of the scope of this paper. Nevertheless, a simulation varying only
the stellar type and the distance is performed below to confirm low impact on
the acceptable noise levels and find an order-of-magnitude estimate for the
relative change. To simulate the observation around a different stellar type,
an Earth-twin is placed around HD 232979171717Effective temperature
$T_{\mathrm{eff}}=$4047\text{\,}\mathrm{K}$$, effective radius
$R_{\mathrm{eff}}=0.57\,\mathrm{R_{\odot}}$, luminosity
$L=0.079\,\mathrm{L_{\odot}}$ (Gaia Collaboration et al. 2016, 2018), planet
separation scaled to $a=0.31$ according to an planetary equilibrium
temperature of $T_{\mathrm{eq}}=$265\text{\,}\mathrm{K}$$, an M0.5V-type star
at 10 pc distance. For this target as well as for an Earth-twin system at 5 pc
distance, the ratio of fundamental to instrumental noise over amplitude and
OPD perturbations is shown in Figure 17.

Figure 17: Impact of the target system configuration on the ratio of
fundamental noise to instrumental noise on a grid of amplitude and phase
perturbations for the shortest-wavelength bin at
$\lambda_{b}=$4.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ with width
$\Delta\lambda_{b}=$0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The system
parameters assumed for the M-star target are given in the text. The contours
indicate where the fundamental noise dominates the systematic noise by a
factor of one, five or ten respectively. The cross indicates the
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ reference level specified by Table
11.

It is apparent that the change in system configuration undertaken here changes
the acceptable perturbation levels by less than one order of magnitude. Away
from the amplitude and phase limits (towards the center of Figure 17), the
acceptable levels show almost no change. This can be taken as validation for
the approach of this appendix to retrieve order-of-magnitude estimates using
only a single representative Earth-twin target. However it must also be a
reminder that this will propagate systematic biases into the sample of
detectable planets, making an analysis of the full sample a top-priority for
future work. An additional interesting property found here is that decreasing
the distance to the target increases the strictness of the OPD error
requirement while decreasing the strictness of the amplitude error
requirement. The same effect can be seen when moving to hotter host stars.

###  A.4 Discussion

In Table 15 and Table 16, we present an estimate for the maximum allowed noise
level to stay in the fundamental noise limited regime. We have cross-validated
the simulations producing these results with Lay (2004). While the maximum
allowed noise levels can confidently be used as an order-of-magnitude
approximation, we briefly discuss the most prominent factors that will
influence these values as work on LIFE continues to progress. First, in Figure
17 we demonstrate three different configurations of the target systems which
do not produce significantly different requirements towards instrument
perturbations. However, the _LIFE_ mission aims at observing a much more
diverse sample than what is represented in Figure 17. Therefore, in order to
gain an understanding about the requirements to remain in the fundamental
noise limited regime for every target, we will apply the presented simulation
to the full synthetic planetary sample described in Quanz et al. (2021) in
future work. Secondly, the presented analysis assumes a specific type of beam
combination (double Bracewell). The _LIFE_ initiative is still in the process
of evaluating other combiner techniques (Hansen et al. 2022). An adoption of a
significantly different technique would certainly also change the requirements
on maximum allowed perturbations and the types of perturbations that dominate
the noise. Thirdly, the assumptions about the perturbations themselves are not
directly based on real data or a detailed instrument model. On the one hand,
an approach more strongly based on experiments could inform a more detailed
model for the shape of the perturbation spectra. As is shown in Defrère et al.
(2010) the shape can significantly influence the requirements for maximal
perturbations. On the other hand, it would enable the connection of
requirements as presented in this report to a component level breakdown. This
would enable an evaluation of the feasibility of reaching the required levels
of perturbation.

Meeting these requirements for amplitude and especially OPD error described in
Table 15 will be a considerable challenge. For comparison, currently planned
lab experiments in the context of the _LIFE_ mission (Gheorghe et al. 2020)
plan to achieve similar amplitude mismatch. The planned OPD error however is
one order of magnitude larger compared to what is presented in this work.
Based on these experiments it will be difficult to predict whether the
requirements in Table 15 are achievable in the mid-term future. In turn, this
raises the question if it is strictly required to operate the instrument in
the fundamental noise limited regime over the whole wavelength range to
accomplish the science goals proposed for the _LIFE_ mission. Allowing for
more OPD and amplitude perturbation will increase the instrumental noise at
short wavelengths. However, these wavelength bands are routinely dominated by
stellar geometric leakage (fundamental noise) and offer very little signal at
least for temperate exoplanets (see Figure 5). Accordingly, for those targets,
wavelengths below $6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ do not contribute
to the bulk detection nor the atmospheric retrieval (cf. Konrad et al. 2021).
We therefore expect that allowing for appropriately larger perturbations will
not significantly alter the results presented in this paper and in Quanz et
al. (2021).

##  Appendix B Implementation of maximum likelihood method

Figure 18: Different versions of the cost function $J$, which is calculated on
a grid up to $100\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ angular separation
from the star. An artificial planet with a photon-based SNR$\sim 5$ is located
in the lower right at $67\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ from the
center. Left: $J^{\prime}$ without positivity constraint on the flux. The
image is symmetric and the contrast is low such that the planet position
cannot be inferred. Middle: $J^{\prime\prime}$ with a simple positivity
constraint on the planet flux. The spatial degeneracy is removed and the most
likely planet position in the lower right indicates the correct planet
position. Two side lobes with similar cost value to the true planet position
are visible. Right: Cost function $J_{\mathrm{reg}}$ with positivity
constraint on the flux and spectral regularization. The overall contrast is
enhanced and the side lobes are suppressed compared with the unregularized
cost function.

We describe the implementation of the maximum likelihood method for the
setting described in Sect. 2.3, with symmetrical noise sources that contribute
by shot noise only. The time-dependence of the transmission map originates
only from the rotation of the array, such that
$T(t,\lambda,\vec{\theta}_{\mathrm{p}})=T(\phi(t),\lambda,\vec{\theta}_{\mathrm{p}})$.
The ideal modulation profiles for all possible planet positions
$\vec{\theta}_{\mathrm{p}}$ are already contained in the transmission map
shown earlier in Fig. 3 as horizontal lines. Additionally, the noise variance
is time-independent: $\sigma^{2}(t,\lambda)=\sigma^{2}(\lambda)$ and can be
estimated from the recorded data $A(t,\lambda)$.

For simplicity, here $J$ labels the total cost function including the
regularization term. To enforce smoothness on the spectrum, the second order
derivative $m=2$ is used. The discretized cost function is given by

$\displaystyle\begin{split}J(\vec{\theta},\vec{F})=\sum_{\lambda}\Bigg{[}&\sum_{t}-2F_{\lambda}\sum_{t}\dfrac{A_{t,\lambda}\,T_{t,\lambda}}{\sigma^{2}_{\lambda}}+F_{\lambda}^{2}\sum_{t}\dfrac{T^{2}_{t,\lambda}}{\sigma^{2}_{\lambda}}\\\
&+\mu\,\left(\dfrac{F_{\lambda+1}-2F_{\lambda}+F_{\lambda-1}}{(\Delta\lambda)^{2}}\right)^{2}\Bigg{]}\end{split}$
(61)

where the dependence of $A$ and $T$ on the considered position
$\vec{\theta}_{\mathrm{p}}$ is omitted in the notation on the right hand side.
To switch to matrix notation, the following parameters are defined:

$\displaystyle\vec{c}$
$\displaystyle:=\\{c_{\lambda}\\}\,;\qquad\quad\mathrm{where}\;c_{\lambda}=\sum_{t}A_{t,\lambda}\,T_{t,\lambda}/\sigma^{2}_{\lambda}$
(62) $\displaystyle\mathbf{B}$
$\displaystyle:=\operatorname{diag}\\{B_{\lambda}\\}\,;\quad\mathrm{where}\;B_{\lambda}=\sum_{t}T^{2}_{t,\lambda}/\sigma^{2}_{\lambda}$
(63)

$\mathbf{D}=\pmatrix{-}$ (64) &1⋯⋯⋯⋯ 1-21⋯⋯⋯ ⋮1-21⋯⋯ ⋮⋮⋱⋱⋱⋮ ⋮⋮⋮⋱-21 ⋮⋮⋮⋮1-1

where $\vec{c}$ corresponds to the unregularized flux estimate vector, and
$\mathbf{D}$ is the second order derivation matrix, which is symmetric.
Equation (61) can then be rewritten in matrix notation as

$\displaystyle J(\vec{\theta},\vec{F})$
$\displaystyle=-2\vec{F}\cdot\vec{c}+\vec{F}^{\intercal}\mathbf{B}\vec{F}+\mu\,\lVert\mathbf{D}\vec{F}\rVert^{2},$
(65)

The condition for minimizing the cost function with respect to the flux, Eq.
(45), then directly returns the optimized flux vector as

$\displaystyle\vec{\hat{F}}=(\mathbf{B}+\mu\,\mathbf{D}^{2})^{-1}\vec{c}$ (66)

which now also includes the regularization enforced by $\mathbf{D}$. If the
optimized flux is inserted back into Eq. (65), we get the cost function
$J^{\prime}(\vec{\theta})=J(\vec{\theta},\vec{\hat{F}})$, which is optimized
with respect to the flux:

$\displaystyle J^{\prime}(\vec{\theta},\vec{\hat{F}})$
$\displaystyle=-\vec{c}^{\intercal}(\mathbf{B}+\mu\,\mathbf{D}^{2})^{-1}\vec{c}=-\vec{\hat{F}}\cdot\vec{c}$
(67)

where only →c depends on the recorded data. As with the unregularized cost
function in Sect. 3.1, a positivity constraint can be applied on the
regularized estimated flux:

$\displaystyle J^{\prime\prime}(\vec{\theta})$
$\displaystyle=J[\vec{\theta},\vec{\hat{F}}_{\mathrm{pos}}(\vec{\theta})]=-\vec{\hat{F}}_{\mathrm{pos}}\cdot\vec{c}$
(68)

The computation of the likelihood maps can be performed on any grid of
possible planet positions and does not necessarily have to cover the whole
field-of-view of the telescope.

In Fig. 18 we show a concrete example for how the different cost functions
behave.

##  Appendix C Local zodiacal emission

Figure 5 indicates that the local-zodiacal thermal emission is the dominant
astrophysical noise source in _LIFE’s_ long-wavelength regime. Figure 19
depicts a visual representation this emission according to Formula (17). It
shows that with respect to local-zodi leakage it is generally favourable to
observe targets at high ecliptic latitudes and longitudes close to the anti-
sun direction.

Figure 19: Surface brightness distribution of local zodiacal dust cloud for
wavelengths of $5\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (from top to bottom). The flux is
given in MJy sr-1 to allow for a comparison with the DarwinSim technical
report (Den Hartog & Karlsson 2005) and other models for the local zodiacal
dust, e.g., the one used for sensitivity predictions for MIRI instrument
onboard the _JWST_ (Glasse et al. 2015). We note the different colorscale for
each wavelength.