# Impact of the noise knowledge uncertainty for the science exploitation of
cosmological and astrophysical stochastic gravitational wave background with
LISA

Martina Muratore  contact<EMAIL_ADDRESS>Jonathan Gair  Lorenzo
Speri

###### Abstract

This paper investigates the impact of a lack of knowledge of the instrumental
noise on the characterisation of stochastic gravitational wave backgrounds
with the Laser Interferometer Space Antenna (LISA). We focus on constraints on
modelled backgrounds that represent the possible backgrounds from the mergers
of binary black holes of stellar origin, from primordial black hole
generation, from non-standard inflation, and from sound wave production during
cosmic fluid phase transitions. We use splines to model generic, slowly
varying, uncertainties in the auto and cross-spectral densities of the LISA
time delay interferometry channels. We find that allowing for noise knowledge
uncertainty in this way leads to one to two orders of magnitude degradation in
our ability to constrain stochastic backgrounds, and a corresponding increase
in the background energy density required for a confident detection. We also
find that to avoid this degradation, the LISA noise would have to be known at
the sub-percent level, which is unlikely to be achievable in practice.

##  I Introduction

The Laser Interferometer Space Antenna (LISA) is part of the European Space
Agency Cosmic Vision program and is due to be launched in the mid-2030s. LISA
will be the first observatory in space to study gravitational waves (GWs) at
mHz frequencies. It will consist of a constellation of three satellites
forming a quasi equilateral triangle and continuously exchanging laser beams
LISA:2017pwj . LISA is expected to observe a large variety of sources, such as
galactic binaries (GBs), massive black hole binaries (MBHBs) (Klein:2015hvg,
), stellar-origin black hole binaries (SOBHB) (PhysRevLett.116.231102, ;
Gerosa:2019dbe, ; Moore:2019pke, ), extreme-mass-ratio inspirals (EMRIs)
(Babak:2017tow, ) and possibly stochastic backgrounds arising from
astrophysical and cosmological processes Caprini_2019 .

When considering the science that can be done with LISA, it is typical to
assume a known model for the instrumental noise in the detector data channels.
However, these noise levels will not be known in practice. This is also true
for ground-based gravitational wave detectors, but in that context spectral
density estimation is easier because signals are rare and short-lived,
allowing the spectral density to be estimated from data in the vicinity of
observed events. LISA signals, by contrast, are typically long lived, which
means that noise and signal properties must be simultaneously estimated by
fitting a suitable model. While such methods and models are still under
development, it is expected that the characterisation of deterministic signals
will not be significantly affected by lack of instrumental noise knowledge
(see Appendix B.4). The case of stochastic GW backgrounds (SGWB) is different,
however, as these are intrinsically of the same character as the stochastic
instrumental noise. Searches for stochastic signals in ground-based
interferometers rely on the cross-correlation of data from independent
detectors 1999PhRvD..59j2001A . This would only be possible if there is
another space-based interferometer in operation concurrently, such as Taiji
2018arXiv180709495R , but this is not certain at the moment. Here we explore
the challenge of distinguishing between the stochastic instrumental noise and
a stochastic GW signal.

One approach is to use a model for the instrumental noise. It is possible to
derive analytical models that describe how different known noise sources
propagate into the LISA data stream. However, not all noise sources will be
known in advance, so we will not be able to strictly rely on the models as we
will not be able to perform full tests and directly measure the noise. In the
LISAPathfinder mission PhysRevLett.116.231101  it was seen that at low
frequency the analytical models couldn’t fully explain the measured noise.
Therefore, when we plan for LISA data analysis, we must be prepared for
uncertainty in the noise models.

The goal of this paper is to assess the impact of lacking a noise model for
LISA in parameter estimation of SGWBs. We consider four different models of
cosmological and astrophysical SGWBs: a power law to model signals from
stellar origin binary black hole inspirals, a Gaussian bump to model a
background from primordial black hole generation, a power law with running to
model a background from non-standard inflation and finally a first order phase
transition model, representing GW production from sound waves in the cosmic
fluid generated by colliding phase transition bubbles Caprini_2019 . For each
model, we will take a reference amplitude that corresponds to a relatively low
signal to noise ratio (SNR) that is close to the boundary for detection. These
are the backgrounds that will be most difficult to distinguish from
instrumental noise. We will also explore what happens as the background energy
density is varied in each model.

We represent our lack of knowledge of the LISA instrumental noise by
multiplying a set of reference auto- and cross-spectral densities with cubic
splines. For the reference spectral densities we use the noise model from
LISA:2017pwj , which includes only the so-called secondary noises
Muratore:2022nbh , the test mass (TM) acceleration and optical metrology noise
(OMS). This noise model assumes that the laser noise Armstrong_1999 , clock
noise Tinto:2018kij , tilt to length coupling Hartig_2022 ;
PhysRevD.106.042005  have been suppressed by the initial noise reduction
pipeline Hartwig:2021dlc ; Muratore:2021rwq . To represent the fact that we
will have some amount of information from noise modelling before launch, we
place a Gaussian prior on the weights of the cubic spline. By varying the
Gaussian variance we explore the effect of having more or less knowledge of
the noise.

Several previous studies have tackled the problem of detecting a SGWB with
LISA and distinguishing it from the noise, but these have used different
methods than the one we employ in this paper. In Flauger:2020qyi  it was shown
that SGWB reconstruction was possible for generic SGWB models, if the LISA
instrumental noise can be represented by just two parameters, representing the
level of TM and OMS noise, assumed equal for all arms of the interferometer.
The authors of Adams:2010vc ; Hartwig:2023pft  allowed the TM and OMS noises
to differ from arm to arm, but still assumed that these noises had a known
spectral shape as a function of frequency. In Baghi:2023qnq  an arbitrary
noise shape was allowed, described by a spline, but using a simplified noise
model for the single link. Finally, Muratore:2022nbh  derived an upper bound
on the detectable SGWB amplitude when being agnostic on both the signal and
noise shape and discussed limitations of the utility of the null channel for
distinguishing between instrumental noise and a stochastic GW background.

The paper is organised as follows: in Section II we introduce the general data
model that we use in the analysis and we describe the Fisher matrix formalism
that will be used for this work. In Section II.3 we describe the spline model
that we use to represent the uncertainties in the power spectral density (PSD)
and cross-spectral density (CSD) of the instrumental noise. In Section II.4 we
give the analytical noise model for a single LISA link that is used as the
reference model, and the corresponding PSDs and CSDs for the time delay
interferometry (TDI) channels, $A$, $E$ and $\zeta$. In Section II.5 we
describe how a stochastic signal appears in the three TDI channels and their
cross-correlations, while in Section II.6 we describe the models for the
cosmological and astrophysical SGWBs that we use in this paper. In Section
III.1 we show how well we can estimate the parameters of the different SGWB
models when we allow for uncertainty in our knowledge of the instrumental
noise. For each model, we compare the precision of parameter estimation to
that when noise knowledge is perfect and show how the parameter precisions
vary as a function of the background energy density, $\Omega$, evaluated at
1mHz. In Section III.2 we show how the results change as we vary our priors
uncertainty on the instrumental noise. We conclude our results in Section
III.3 by showing how well the signal, noise and galactic foreground can be
reconstructed for a power law SGWB background. Section IV summarises our
conclusions and future perspectives.

##  II Methods

###  II.1 Likelihood

We assume that the output of a gravitational wave detector, $s(t)$, is
expressed as a linear combination of a signal, $h(t|\vec{\mu})$, determined by
a finite set of (unknown) parameters, $\vec{\mu}$, and instrumental noise,
$n(t)$. If we ignore the presence of calibration errors Savalle_2022 , the
content of a single data stream, i.e., one output channel from one detector,
can be written in the frequency domain as:

$\tilde{s}(f)=\tilde{h}(f|\vec{\mu})+\tilde{n}(f),$ (1)

where the tilde indicates the Fourier transform. The likelihood for the
observed data can be written as
$p(\tilde{s}(f)|\vec{\mu})=p(\tilde{n}(f)=\tilde{s}(f)-\tilde{h}(f|\vec{\mu}))$.
In a gravitational wave context it is usual to further assume that the
instrumental noise follows a Gaussian distribution characterized by a one-
sided PSD, $S_{n}(f)$, defined such that

$\mathbb{E}[\tilde{n}^{*}(f)\tilde{n}(f^{\prime})]=\frac{1}{2}S_{n}(f)\delta(f-f^{\prime}),$
(2)

for $f$, $f^{\prime}>0$, where the expectation value $\mathbb{E}$ is taken
over the data generating process. The delta function in the previous equation
implies that different frequencies are not correlated.

In reality the noise model is not known perfectly and could vary from the
assumption above in a number of ways. For example, the PSD might have a
different shape from the reference one babak2021lisa , the probability
distribution of the noise might not be Gaussian, or the noise might not be
stationary, leading to correlations between frequencies.

In this work we will continue to assume that the noise is Gaussian and
stationary, but we will allow the power spectral density to vary using a
parametrized spectral density, $S_{n}(f)\rightarrow S_{n}(f|\vec{\lambda})$,
described by parameters $\vec{\lambda}$. Then the log-likelihood depends on
both sets of parameters, $\vec{\mu}$ and $\vec{\lambda}$, and can be written
as:

$\displaystyle l:=\ln
p(\tilde{s}|\vec{\mu},\vec{\lambda})=-\sum_{k=1}^{n}\ln\left[2\pi\frac{S_{n}(f_{k}|\vec{\lambda})}{4\Delta
f}\right]-$ (3)
$\displaystyle\frac{1}{2}\sum_{k=1}^{n}\frac{|\tilde{s}(f_{k})-\tilde{h}(f_{k}|\vec{\mu})|^{2}}{\frac{1}{4\Delta
f}S_{n}(f_{k}|\vec{\lambda})}$ (4)

where the sum is performed over $n$ frequencies and
$\tilde{n}(f_{k})=\tilde{s}(f_{k})-\tilde{h}(f_{k}|\vec{\mu})$ are the
discrete Fourier component at frequency $f_{k}=k\Delta f$, of the data minus
signal model. The frequency bin width, $\Delta f$, is related to the total
observation time as $T=1/\Delta f$. The first term does not include the $1/2$
factor because the real and imaginary parts of $\tilde{n}(f_{k})$ are
independent random variables. This follows from the fact that, for a real time
series, $\tilde{n}^{*}(f)=\tilde{n}(-f)$, which combined with Eq. (2) means
that $\langle\tilde{n}(f)\tilde{n}(f^{\prime})\rangle=0$ for $f$,
$f^{\prime}>0$. This allows equation 2 to be rewritten as

$\langle\Re[\tilde{n}(f_{k})]^{2}\rangle=\langle\Im[\tilde{n}(f_{k})]^{2}\rangle=\frac{\langle|\tilde{n}(f_{k})|^{2}\rangle}{2}=\frac{S_{n}(f_{k}|\vec{\lambda})}{4\Delta
f}$ (5)

for a discrete set of frequencies.

Stochastic gravitational wave backgrounds are not deterministic signals and
can be treated on the same footing as the instrumental noise by defining the
total variance at frequency $f_{k}$ as

$S_{t}(f_{k}|\vec{\theta},\vec{\lambda})=S_{\rm GW}(f_{k}|\vec{\theta})+S_{\rm
n}(f_{k}|\vec{\lambda})\,.$ (6)

If we assume that all the deterministic sources have been correctly subtracted
from the datastream $s$ the log-likelihood becomes:

$\displaystyle l(\vec{\theta},\vec{\lambda})$
$\displaystyle=-\sum_{k=1}^{n}\ln\left[T\pi\frac{S_{t}(f_{k}|\vec{\theta},\vec{\lambda})}{2}\right]$
$\displaystyle-\frac{1}{2}\sum_{k=1}^{n}\frac{|\tilde{s}(f_{k})|^{2}}{\frac{T}{4}S_{t}(f_{k}|\vec{\theta},\vec{\lambda})}$
(7)

The derivation of this likelihood can be found in Appendix A.

###  II.2 Fisher matrix

We are interested in understanding the impact of noise knowledge uncertainties
on the parameter measurement precision of stochastic gravitational wave
backgrounds. The Fisher information matrix provides a lower bound on the
covariance of an unbiased estimator of the model parameters and provides a
good approximation to the precision of parameter estimation in the high
signal-to-noise ratio limit. We will therefore use it to quantify our ability
to measure both the noise parameters, $\vec{\lambda}$, and the background
parameters, $\vec{\theta}$.

In a general context the Fisher matrix is defined by

$\Gamma_{ij}=\mathbb{E}\left[\frac{\partial
l}{\partial\upsilon^{i}}\frac{\partial
l}{\partial\upsilon^{j}}\right]=-\mathbb{E}\left[\frac{\partial^{2}l}{\partial\upsilon^{i}\partial\upsilon^{j}}\right]$
(8)

where the expectation value $\mathbb{E}$ is taken over the data generating
process, and the partial derivatives are taken with respect to the parameters,
$\vec{\upsilon}$, on which the likelihood depends. We want to compute the
Fisher matrix on the extended parameter space
$\vec{\upsilon}=\\{\vec{\theta},\vec{\lambda}\\}$.

It can be shown that the expectation value of the product between the
derivative of the log-likelihood with respect to deterministic and stochastic
parameters is zero. Therefore, at the level of the Fisher matrix approximation
it can be shown that the estimation of the noise and deterministic signal
parameters is independent (see Appendix B.4).

For SGWBs, we can compute the Fisher matrix in continuous domain as:

$\Gamma_{ij}=T\int_{0}^{\infty}(\Sigma^{-1})_{lr}\frac{\partial\Sigma^{rp}}{\partial\upsilon^{i}}(\Sigma^{-1})_{pm}\frac{\partial\Sigma^{ml}}{\partial\upsilon^{j}}\,\rm{d}f\,.$
(9)

with

$\displaystyle\Sigma(f|\vec{\upsilon}=\\{\vec{\theta},\vec{\lambda}\\})=\frac{1}{2}\begin{pmatrix}S^{AA}_{t}&S^{AE}_{t}&S^{A\zeta}_{t}\\\
S^{AE*}_{t}&S^{EE}_{t}&S^{E\zeta}_{t}\\\
S^{A\zeta*}_{t}&S^{E\zeta*}_{t}&S^{\zeta\zeta}_{t}\end{pmatrix}\,.$ (10)

where each element of the matrix can be written as a sum of an instrumental
noise component and a stochastic gravitational wave component as indicated in
Eq. 6. A complete derivation of this formula can be found in Appendix B.

Prior knowledge on the noise can be incorporated by imposing a prior on the
instrumental parameters, $\vec{\lambda}$. When doing numerical marginalisation
any prior can be imposed, but in the Fisher matrix formalism it is easiest to
work with a Gaussian prior Savalle_2022 . The posterior covariance is then
given by the inverse of the modified Fisher matrix:

$\Gamma=\begin{pmatrix}\Gamma^{\theta\theta}&\Gamma^{\theta\lambda}\\\
(\Gamma^{\theta\lambda})^{T}&\Gamma^{\lambda\lambda}+\Theta^{\lambda\lambda}\end{pmatrix}$
(11)

with normal prior on the instrumental noise parameters with zero mean and
covariance given by $(\Theta^{\lambda\lambda})^{-1}$. The diagonal elements of
the inverse of this matrix provide estimates for the precision with which the
corresponding parameters can be measured. The estimated precision of
measurement of the SGWB parameters accounting for noise model uncertainty is
thus given by the diagonal elements of the matrix:

$\sigma_{\theta}=\sqrt{{\rm
diag}[(\Gamma^{\theta\theta}-\Gamma^{\theta\lambda}(\Gamma^{\lambda\lambda}+\Theta^{\lambda\lambda})^{-1}(\Gamma^{\theta\lambda})^{T})^{-1}]}$
(12)

Note that in the limit in which the instrumental noise parameters are
perfectly known $\Theta\rightarrow\infty$ and the measurement precision of the
SGWB parameters is given by $\sigma_{\theta}=\sqrt{{\rm
diag}[(\Gamma^{\theta\theta})^{-1}]}$.

###  II.3 Modeling noise knowledge uncertainties

To model noise uncertainties, we allow the PSD and CSD of the instrumental
channels to deviate from the design specification. However, we assume that
such deviations vary smoothly over a relatively wide range of frequency and
model the noise uncertainties as fractional deviations from the design PSD/CSD
that are described by natural cubic splines. We write the PSD of the
instrumental noise in each channel as:

$S_{n}(f|\lambda)=S_{\rm des}(f)\,10^{C(f|\vec{\lambda})}\,,$ (13)

where $C(f|\vec{\lambda})$ is a natural cubic spline. The parameters
$\vec{\lambda}$ specify the values of the spline at the knots, labelled by
$i$. In this study we use knots evenly spaced in $\log_{10}(f)$ between
$\log_{10}(f)=-4$ and $\log_{10}(f)=0$ and we fix the number of knots to 13.
Noise curves corresponding to this model, with the weights at each knot drawn
randomly from a $\log_{10}(f)\sim U[-1,1]$ distribution, are shown in fig. 1.
We note that this choice of prior means we are allowing approximately one
order of magnitude variation in the PSD. When we evaluate the Fisher matrix we
will always do so at the reference point where the weights of the spline are
zero, i.e., where the PSD is equal to the reference value shown in Fig. 2.

We cannot follow the same procedure for specifying the CSD, because the
reference model is smaller by 1 or even 2-3 order of magnitude a low frequency
with respect to the PSD (see Fig. 3 and Fig. 2). It was shown in
Hartwig:2023pft  that when the LISA response is constructed allowing for
unequal noises in the different laser links, the CSD can be much larger and
become comparable to the PSD. Since our goal is to allow the splines to vary
in such a way to mimic un-expected and un-modelled noise components with
respect to the simplified scenario (three unequal but fixed-length arms) we
model the CSD as111 In principle our model does not force the matrix to be
positive definite. We are forcing the reference spectral density matrix to be
positive definite, but in principle we could have a factor of 10 variation in
the CSD while the PSD is unchanged. It doesn’t matter for the Fisher matrix
because this is a local approximation and we are evaluating it at a point
where the matrix is positive definite. The CSD at the central point is $0.1$
of its maximum value, so in an open set around that point it will be positive
definite and thus all derivatives are well defined. The conclusion is that the
model used here is fair for what we want to demonstrate but would not be a
suitable model to use when analysing the data. :

$\displaystyle S_{n}(f|\\{\lambda_{i}\\})=$ $\displaystyle\sqrt{S_{\rm
des,i}(f)S_{\rm
des,i}(f)}\,\sigma_{R}\,10^{C(f|\\{\log_{10}(f_{i})\\},\\{\lambda_{i}\\})}$
$\displaystyle+$ $\displaystyle\sqrt{S_{\rm des,i}(f)S_{\rm
des,i}(f)}\,i\sigma_{I}\,10^{C(f|\\{\log_{10}(f_{i})\\},\\{\lambda_{i}\\})}\,$
(14)

where we fix $\sigma_{R}=0.1,\sigma_{I}=0.8\sigma_{R}$. We do not expect that
varying the relationship between $\sigma_{I}$ and $\sigma_{R}$ to
significantly change the conclusions; although we fix $\sigma_{I}$ slightly
smaller than $\sigma_{R}$ accordingly to Fig. 3. There at lower frequencies
the imaginary components are about 1 order of magnitude smaller than the real
components. The indexes $i$ and $j$ run over the number of detectors or
channels with $i\neq j$. The additional factors $\sigma_{I}$ and $\sigma_{R}$
are used to limit the amplitude, and allow us to model the CSD as a sum of
splines times the geometric mean of the square root-PSDs. Using the
(scaled-)geometric mean of the PSDs as a reference for the CSD rather than the
CSD of the reference equal-noise configuration, allows for much larger CSD
variations. This is consistent with the results presented in Hartwig:2023pft .

It is important to state that our model is not completely general since we are
imposing a certain amount of smoothness in the PSD variation, and consequently
in the CSD, when we specify the number and spacing of the knots. Thus we are
not able to fit for all possible noise scenarios. In particular, this model
does not attempt to reproduce the zeros of the TDI transfer functions
faithfully. This will become important above $f\sim 0.05$Hz, but this should
not affect our results as the SGWBs we consider do not have much power at
those frequencies, as can be seen from Fig. 6. Other models could be
considered, for example by imposing the spline variations at the level of the
noise in individual laser links (generalising the approach taken in
Baghi:2023qnq ), before applying the TDI transfer function. This should be
explored in the future, but this would increase the number of parameters
further so we might expect there to be additional degeneracies, which would
lead to practical difficulties in fitting noise and signal simultaneously.
However, for the purpose of the current study, the model we use is adequate to
represent generic, slowly varying, fluctuations in the PSD and CSD.

Figure 1:  Deviations from the design power spectral density obtained using
the cubic spline model, with $\lambda_{i}\sim U[-1,1]$, and with knots equally
spaced between $\log_{10}(f)=-4$ and $\log_{10}(f)=0$. The plot shows the
ratio of the total PSD and the design one for different parameter realizations
$\lambda_{i}$.

###  II.4 Noise at the TDI input and outputs

Here, we present the instrumental noise model used to define the reference PSD
in this work. Among the different noise sources for LISA, the laser noise is
the main source of noise, which must be reduced by eight orders of magnitude
by applying a post-processing technique called time delay interferometry (TDI)
tintoD . TDI synthesises an equal arm-length interferometer by appropriately
delaying and combining the interferometric measurements in many different ways
to form TDI channels free from laser noise. The standard second generation TDI
channels (unequal and time varying arm-length) are the Michelson
interferometer channels, $X$, $Y$ and $Z$, from which we form the more GW
sensitive channels, $A$ and $E$ PhysRevD.66.122002 . Together with the GW
sensitive channels we consider a null channel, the $\zeta$ channel
Muratore:2021uqj  that is less sensitive to GWs and can in principle be used
as a noise monitor.

In the current work we will assume that laser noise has already been reduced
thus we can work directly with the first generation TDI PhysRevD.105.062006 ,
allowing us to consider three unequal but fixed arm-lengths. This choice
should not significantly affect the conclusions of the analysis but simplifies
the transfer function derivations.

We also assume that all known calibrated and measured instrumental noise
sources have been subtracted, such as the optical tilt to length cross-
coupling to spacecraft motion and clock noise (Hartig_2022 ,
PhysRevApplied.14.014030  and PhysRevD.103.123027 ). The remaining noises, for
which we have neither a measurement for coherent subtraction nor a high
precision a priori model Muratore:2022nbh , fall into two broad categories,
the acceleration noise of each individual test-mass (TM) and an overall
optical metrology system (OMS) noise term for each single link measurement
(see nam2022tdi  for the case of multiple OMS noise terms).

We represent the TM acceleration noise PSD of a single TM by $S_{g_{ij}}$. To
directly compare the OMS and TM contributions we can directly convert the
acceleration noise of a single TM to an equivalent displacement, whose PSD is
given by

$S_{g_{ij}}^{\rm disp}=S_{g_{ij}}/(2\pi f)^{4}\,$ (15)

where $f$ is the Fourier frequency. We denote the time series associated with
this displacement as $x^{g}_{ij}(t)$. We also define the PSD of the OMS noise
as $S_{{\rm oms}_{ij}}(f)$ and we denote the time series of the single OMS as
$x^{m}_{ij}(t)$. All TDI combinations can be constructed from a combination of
single link TM to TM measurements. Such measurements are represented by the
intermediary variables Hartwig:2021dlc :

$\tilde{\eta}^{N}_{ij}(\omega)=\tilde{x}(\omega)^{g}_{ji}e^{-i\omega
L_{ji}}+\tilde{x}(\omega)^{g}_{ij}+\tilde{x}(\omega)^{m}_{ij},$ (16)

where $\tilde{\eta}^{N}_{ij}(\omega)$ is the noise in a single link
measurement, the first index $i$ indicates the spacecraft where the
measurement is performed at time $t$, and the second index $j$ indicates the
distant spacecraft from which light was emitted at time $t-\tau$, and
$\omega=2\pi f$. Equation 16 implies that each single link measurement
contains TM noise terms from the distant and local spacecraft, such that the
TM noise appearing in the measurements on the two ends of the same arm is
correlated (between the two links):

$\langle\tilde{\eta}^{N}_{ij}(\omega)\tilde{\eta}^{N}_{ji}(\omega)\rangle\neq
0$ (17)

From these measurements it is possible to build any TDI channels
Armstrong_1999 ; muratore2020revisitation  and therefore the corresponding
first generation orthogonal channels $A_{1}$ and $E_{1}$ PhysRevD.105.062006
that will be used in this work:

${\rm A_{1}}=\frac{{\rm Z_{1}}-{\rm X_{1}}}{\sqrt{2}}\;,\qquad{\rm
E_{1}}=\frac{{\rm X_{1}}-2{\rm Y_{1}}+{\rm Z_{1}}}{\sqrt{6}}.$ (18)

The $X_{1}$ variable is defined as:

$\displaystyle{\rm X_{1}}=$
$\displaystyle(D_{13}D_{31}-1)(\eta_{12}+D_{12}\eta_{21})$
$\displaystyle+(1-D_{12}D_{21})(\eta_{13}+D_{13}\eta_{31})\;,$ (19)

where the delays $D_{ij}$ corresponds to a constant time shift and thus in
frequency to $\mathcal{F}\\{D_{ij}\\}=e^{-i\omega L_{ij}}$. The $Y_{1}$ and
$Z_{1}$ are given from $X_{1}$ by cyclic permutations of the three satellites.
The fully symmetric channel, $\zeta_{1}$, is defined by:

$\zeta_{1}=D_{12}(\eta_{31}-\eta_{32})+D_{23}(\eta_{12}-\eta_{13})+D_{31}(\eta_{23}-\eta_{21})\;.$
(20)

The assumed model for the TM acceleration noise is

$\displaystyle\operatorname{E}\left\langle\tilde{x}^{g}_{ij}(f)\tilde{x}^{g*}_{lm}(f^{\prime})\right\rangle=\frac{1}{2}\delta_{il}\delta_{jm}\delta(f-f^{\prime})S_{g_{ij}}(f)$
$\displaystyle S_{g_{ij}}(f)=\left(3\times
10^{-15}~{}\frac{\textrm{m}}{\textrm{s}^{2}~{}\sqrt{\textrm{Hz}}}\right)^{2}$
(21)
$\displaystyle\times\left(1+\left(\frac{0.4~{}\textrm{mHz}}{f}\right)^{2}\right)\left(1+\left(\frac{f}{8~{}\textrm{mHz}}\right)^{4}\right),$

and for the OMS noise

$\displaystyle\operatorname{E}\left\langle\tilde{x}^{m}_{ij}(f)\tilde{x}^{m*}_{lm}(f^{\prime})\right\rangle=\frac{1}{2}\delta_{il}\delta_{jm}\delta(f-f^{\prime})S_{{\rm
oms}_{ij}}(f)$ $\displaystyle S_{{\rm
oms}_{ij}}(f)=\left(15~{}\textrm{pm}/\sqrt{\textrm{Hz}}\right)^{2}\times\left(1+\left(\frac{2~{}\textrm{mHz}}{f}\right)^{4}\right),$
(22)

This model assumes that individual noise components are uncorrelated. In
reality the test masses in the same satellite will share environmental noise,
such as temperature fluctuations, so this assumption might not hold. However,
this model serves as a reference one, and any variation is captured by the
flexible spline model previously presented.

To derive the PSDs and CSDs of the TDI channels $A$, $E$ and $\zeta$, as
illustrated in Muratore:2021uqj , one can express the arm-lengths $L_{ij}$ in
terms of the breathing modes of the LISA triangle, $\delta_{a}$ and
$\delta_{b}$ as:

$\displaystyle L_{12}(t)$
$\displaystyle=L\left[1+\frac{1}{2}\left(\sqrt{3}\,\delta_{a}-\delta_{b}\right)\right]\;,$
(23a) $\displaystyle L_{23}(t)$ $\displaystyle=L\left(1+\delta_{b}\right)\;,$
(23b) $\displaystyle L_{31}(t)$
$\displaystyle=L\left[1-\frac{1}{2}\left(\sqrt{3}\,\delta_{a}+\delta_{b}\right)\right]\;.$
(23c)

The full expressions are rather long thus we give them in a separate
$Mathematica$ notebook file (noise-analytical-model) and we plot them in Fig.
2. The Amplitude spectral density (ASD) are computed in terms of $\delta_{a}$
and $\delta_{b}$. Indeed, while
$L=(L_{12}+L_{23}+L_{31})/3\approx$8.3\text{\,}\mathrm{s}$$ is the average
arm-length, the small parameters $\delta_{a}$ and $\delta_{b}$ are typically
$\sim 0.005-0.009$ for realistic ESA orbits lisa_orbit . The case
$\delta_{a}=\delta_{b}=0$ corresponds to the equal-arm LISA scenario.

Figure 2: Reference Amplitude spectral density for the Time delay
interferometry channels A, E and $\zeta$ considering only test mass
acceleration and optical metrology noise and assuming a constellation of three
fixed unequal arm-lengths.  Figure 3:  Real and imaginary part of the
reference square root of the cross spectral density for the time delay
interferometry channels AE, E$\zeta$ and A$\zeta$ considering only TM
acceleration and OMS noise and assuming a constellation of three fixed unequal
arm-lengths.

We consider that the six TMs have the same PSD as well as the six OMS noise
terms, but this suppresses the contribution in the CSD. It was shown in
Hartwig:2023pft  that if the levels of the noises differ by $20\%$ then the
CSD can be $10\%$ of the PSD at low frequencies and several tens of percent at
high frequency. This motivates the particular choice of flexible CSD model
that we introduced in Eq. (14) and is illustrated in Fig. 4.

Figure 4:  Real and Imaginary part of splines-square root cross spectral
density for the time delay interferometry channels AE, E$\zeta$ and A$\zeta$
considering test mass acceleration and optical metrology noise assuming a
constellation of three fix unequal arm-lengths

###  II.5 Signal transfer function

The detector response to a stochastic background can be computed by expressing
a GW signal as a superposition of plane waves, and by assuming that the LISA
constellation has static arm lengths and is in a flat background spacetime.
Following Hartwig:2023pft , it is possible to show that the component of the
single link measurement $\eta_{ij}(t)$ due to a GW is given by:

$\displaystyle\eta^{GW}_{ij}(t)$
$\displaystyle=i\int_{-\infty}^{\infty}\left\\{\frac{f}{f_{ij}}{\rm e}^{2\pi
if(t-L_{ij})}\right.$ $\displaystyle\left.\int\left[{\rm e}^{-2\pi
if\hat{k}\cdot\vec{x}_{i}}\sum_{\mathcal{A}}\xi_{ij}^{\mathcal{A}}(f,\hat{k})\tilde{h}_{\mathcal{A}}(f,\hat{k})\right]{\rm
d}\Omega_{\hat{k}}\right\\}{\rm d}f,$ (24)

where $i$ stands for imaginary component, $f_{ij}=(2\pi L_{ij})^{-1}$,
$\vec{x}_{i}$ denotes the position of satellite $i$, $\mathcal{A}=+,\,\times$
denotes the GW polarization, $\tilde{h}_{\mathcal{A}}(f,\hat{k})$ is the
Fourier transform of the GW signal, $f$ is the GW frequency, $\hat{k}$ the
outward vector in the direction of the incoming GW and ${\rm
d}\Omega_{\hat{k}}$ is the infinitesimal solid angle.
The above expression quantifies the fractional frequency shift due to a
superposition of plane waves coming from different directions $\hat{k}$.
The term $\xi^{\mathcal{A}}_{ij}$ projects the incoming wave with polarization
$\mathcal{A}$ onto the detector, and its functional dependence is given by:

$\xi^{\mathcal{A}}_{ij}\left(f,\hat{k}\right)=e^{-2\pi
if\hat{k}\cdot\vec{L}_{ij}}\mathcal{M}_{ij}(f,\hat{k})\;\mathcal{G}^{\mathcal{A}}(\hat{k},\hat{l}_{ij})\,,$
(25)

where

$\mathcal{M}_{ij}(f,\hat{k})\equiv\mathrm{e}^{\pi
ifL_{ij}(1+\hat{k}\cdot\hat{l}_{ij})}\;\mathrm{sinc}\left(\pi
fL_{ij}(1+\hat{k}\cdot\hat{l}_{ij})\right)\\\ $ (26)

and

$\mathcal{G}^{\mathcal{A}}(\hat{k},\hat{l}_{ij})\equiv\frac{\hat{l}^{a}_{ij}\hat{l}^{b}_{ij}}{2}e^{\mathcal{A}}_{ab}(\hat{k})\,,$
(27)

where $\hat{l}_{ij}=(\vec{x}_{j}-\vec{x}_{i})/|\vec{x}_{j}-\vec{x}_{i}|$ is a
unit vector pointing from spacecraft $i$ to $j$ and
$e^{\mathcal{A}}_{ab}(\hat{k})$ denotes the GW polarization tensors.

For an a homogeneous, isotropic and non-chiral, stochastic background, the GW
signal is only specified statistically

$\langle\tilde{h}_{\mathcal{A}}(f,\hat{k})\,\tilde{h}_{B}^{*}(f^{\prime},\hat{k}^{\prime})\rangle=\delta(f-f^{\prime})\delta(\hat{k}-\hat{k}^{\prime})\delta_{\mathcal{A}B}\frac{P_{h}^{\mathcal{A}B}(f)}{16\pi}$

and

$\langle\tilde{h}_{\mathcal{A}}(f,\hat{k})\,\tilde{h}_{B}(f^{\prime},\hat{k}^{\prime})\rangle=0\,.$

Homogeneity and isotropy implies that $P_{h}^{\mathcal{A}B}(f)$ is diagonal,
whereas the non-chirality implies $P^{\times\times}_{h}=P^{++}_{h}$, so that
we can define $P_{h}:=\sum_{\mathcal{A}}P^{\mathcal{A}\mathcal{A}}_{h}$.

We characterise the response of the individual links to a stochastic
background statistically

$\langle\tilde{\eta}^{GW}_{ij}\,\tilde{\eta}^{GW}_{mn}\rangle=\frac{1}{2}S^{\eta,\mathrm{GW}}_{ij,mn}(f)\delta(f-f^{\prime})\,,$
(28)

where spectral densities for the link measurements are given by

$S^{\eta,\mathrm{GW}}_{ij,mn}(f)=\frac{f^{2}}{f_{ij}f_{mn}}e^{-2\pi
if(L_{ij}-L_{mn})}\sum_{\mathcal{A}}P_{h}^{\mathcal{A}\mathcal{A}}(f)\;\Upsilon_{ij,mn}^{\mathcal{A}}(f)\;,$
(29)

with:

$\Upsilon_{ij,mn}^{\mathcal{A}}(f)=\int\frac{\rm{d}\Omega_{\hat{k}}}{4\pi}\;\textrm{e}^{-2\pi
if\hat{k}\cdot(\vec{x}_{i}-\vec{x}_{m})}\;\xi^{\mathcal{A}}_{ij}(f,\hat{k})\,\xi^{\mathcal{A}}_{mn}(f,\hat{k})^{*}\,.$
(30)

The power spectral densities of the signal in the TDI variables described in
sec. II.4 can then be computed from

$\displaystyle\langle\tilde{U}(f)\tilde{V}^{*}(f^{\prime})\rangle$
$\displaystyle=\frac{1}{2}S^{\rm GW}_{UV}(f)\delta(f-f^{\prime})$
$\displaystyle S^{\rm GW}_{UV}(f)$
$\displaystyle=\sum_{ij,mn\in\mathcal{I}}c^{U}_{ij}(f)c^{V*}_{mn}(f)S^{\eta,\mathrm{GW}}_{ij,mn}(f)\;,$
(31)

where $\tilde{U}$ and $\tilde{V}$ denote any two TDI variables, which in our
case are TDI A, E and $\zeta$, and $\mathcal{I}=\\{12,13,23,21,31,32\\}$
denotes the set of pairs of indices that define the six inter-satellite links.
The coefficients $c^{U}_{ij/mn}$ map the single-link measurements onto the TDI
variable $U$. Refer to the Mathematica code for the computation of such
coefficients (noise-analytical-model).

Note that considering each polarization of the SGWB contributes equally to the
background, i.e. $P^{\times\times}_{h}=P^{++}_{h}$, we can rewrite Eq. (31) as
a product of the SGWB spectral density $P_{h}(f)$ and a transfer function
$\mathcal{T}^{GW}(f)$ which takes into account the LISA detector response,
i.e. $S^{\rm GW}_{UV}(f)=\mathcal{T}^{GW}(f)P_{h}(f)$.
$S^{\rm GW}_{UV}(f)$ would correspond to the first term on the right hand side
of Eq. (6). The transfer functions for the three TDI channels and their cross
correlation are shown in Fig.5.

Figure 5: Upper panel: Gravitational wave transfer functions
$\mathcal{T}^{GW}(f)$ of the three time delay interferometry channels A,E and
$\zeta$ assuming a constellation of three fixed unequal arm-lengths; lower
panel: real and imaginary components of the gravitational wave transfer
functions $\mathcal{T}^{GW}(f)$ of the time delay interferometry channels AE,
E$\zeta$,A$\zeta$.

###  II.6 SGWB Signal models

There are a large variety of models for stochastic gravitational wave
backgrounds that might manifest in the LISA band red_b . In this work we focus
on four models, which can be described by their energy density,
$h^{2}\Omega_{GW}$ Caprini:2018mtu , which is a function of some parameters,
$\theta$:

  * •

Power law:

$h^{2}\Omega_{GW}(f)\approx A\left(\frac{f}{f_{p}}\right)^{n},$ (32)

where $f_{p}$ is the pivot frequency, defined as the geometrical mean of the
LISA frequency interval ($10^{-4}$ Hz, $0.1$ Hz), $f_{p}=3$mHz. The model
parameters are the log-amplitude, $A$, and slope, $n$ babak2023stochastic ;
lehoucq2023astrophysical . We use reference values of $n=2/3$ and
$A=7.87\times 10^{-13}$, representing a SGWB from stellar origin black hole
binaries that has energy density at 1mHz of $h^{2}\Omega_{GW}(1mHz)=3.78\times
10^{-13}$. This value was chosen to be compatible with LIGO-VIRGO-KAGRA
constraints babak2023stochastic

  * •

Gaussian bump:

$h^{2}\Omega_{GW}=Ae^{-\frac{1}{2\sigma^{2}}\ln(\frac{f}{f_{p}})^{2}},$ (33)

where $f_{p}$ is the pivot frequency as before. The model parameters are the
log-amplitude, $A$, and width, $\sigma$. We use reference values of
$A=10^{-12.48}$ and $\sigma=0.3$ whose energy density at 1mHz is
$h^{2}\Omega_{GW}(1mHz)=4.05\times 10^{-16}$. This signal is chosen as a
simple way to mimic the one which might arise from particle production taking
place for a limited number of e-folds during inflation (as, for instance,
required by some models of primordial BH generation). (See e.g. Bartolo_2019 ;
Bartolo:2018evs ; franciolini2021primordial  )

  * •

Power law with running:

$h^{2}\Omega_{GW}=A\left(\frac{f}{f_{p}}\right)^{n+\alpha\ln(\frac{f}{f_{p}})},$
(34)

where $f_{p}$ is the pivot frequency as before. The model parameters are the
log-amplitude, $A$, the slope, $n$, and the running index, $\alpha$. We use
reference values of $A=10^{-12.65}$, $n=1$ and $\alpha=-0.1$. This signal is
motivated by non-standard inflationary models. For example, gravitational wave
generation can be enhanced by sustained particle production during inflation,
leading to a power law stochastic GW background, which would deviate from a
simple power law at higher frequency when back-reaction kicks in (see e.g.
Bartolo:2016ami ). The energy density at 1mHz is
$h^{2}\Omega_{GW}(1mHz)=6.61\times 10^{-14}$

  * •

First Order Phase Transition:

$h^{2}\Omega_{GW}(f)=h^{2}\Omega_{p}\left(\frac{f}{f_{p}}\right)^{3}\left(\frac{7}{4+3\big{(}\frac{f}{f_{p}}\big{)}^{2}}\right)^{n},$
(35)

where $f_{p}=2\cdot 10^{-4}$ Hz (note this is different to the reference
frequency in the previous models). The model parameters are the energy
density, $h^{2}\Omega_{p}$, and spectral index, $n$. We use reference values
of $A\equiv h^{2}\Omega_{p}=10^{-10}$ and $n=7/2$ whose energy density at 1mHz
is $h^{2}\Omega_{GW}(1mHz)=2.59\times 10^{-12}$. This signal is motivated by
the production of sound waves in the cosmic fluid from colliding phase
transition bubbles Caprini_2016 ; PhysRevD.96.103520 .

In our analysis, we also include the contribution to the spectral density from
the foreground of galactic binaries (GB). We use the following model for the
foreground red_b :

  * •

Foreground of Galactic Binaries:

$\displaystyle S_{GB}(f)=$ $\displaystyle
A_{GB}\left(\frac{f}{Hz}\right)^{-\frac{7}{3}}e^{-(f/f_{1})^{\alpha}}$
$\displaystyle\times\frac{1}{2}\left[1+\tanh\left(\frac{f_{knee}-f}{f_{2}}\right)\right]$
(36)

with:

$f_{1}=10^{a_{1}\log_{10}(T)+b_{1}},f_{knee}=10^{a_{k}\log_{10}(T)+b_{k}}$

.

setting $A=1.15\cdot
10^{-44};~{}\alpha=1.56;~{}a_{1}=-0.15;~{}b_{1}=-2.72;~{}a_{k}=-0.37;~{}b_{k}=-2.49;~{}f_{2}=6.7\times
10^{-4}$Hz. When considering the background in conjunction with other SGWBs we
allow the amplitude to vary, but keep the other parameters fixed.

The relation between the energy density $\Omega_{GW}$ and the stochastic GW
background power spectral density $P_{h}(f)$ is given by Caprini:2018mtu :

$\Omega_{GW}(f)=\frac{4\pi^{2}}{3H_{0}^{2}}f^{3}P_{h}(f),$ (37)

where $H_{0}$ is the Hubble constant fixed to be 67.8 km/s/Mpc, as a
consequence $h=0.678$. The conversion between the energy density
$\Omega_{GW}(f)$ and gravitational power spectral density $P_{h}(f)$ used to
compute Eq.(31) is then Caprini:2018mtu :

$P_{h}(f)=7.98\times
10^{-37}\left(\frac{Hz}{f}\right)^{3}h^{2}\Omega_{GW}(f)\frac{1}{Hz},$ (38)

.

We report in Fig. 6 the ASD of the four SGWB models together with the ASD of
the reference instrumental noise in TDI channel A. babak2021lisa

Figure 6:  Amplitude spectral density of the stochastic GW background models
and the amplitude spectral density of the reference test mass and optical
metrology noise in the time delay interferometry channel A.

We also provide the computation of the SNRs of these different backgrounds in
the TDI channel A using the following formula Thrane:2013oya :

$SNR_{A}=\sqrt{T}\left[\int_{0}^{\infty}\frac{S^{GW}_{AA}(f)^{2}}{S^{A}_{n}(f)^{2}}{\rm
d}f\right]^{1/2}$ (39)

with an observation time span of $T=4$ years. Here $S^{GW}_{AA}(f)$ is the
spectral density in channel $A$ that can be computed from Eq. (31) and
$S^{A}_{n}$ is the PSD of the A channel. The results222We note that this
formula is derived assuming that we have access to two independent channels
that have uncorrelated noise and perfectly correlated signals. This is not a
good approximation to LISA, so the SNR is not directly interpretable. However,
it still provides an indication of the relative detectability of different
backgrounds.  are shown in Table 1.

SGWB Model | SNR w/o GB | SNR w/ GB
---|---|---
Power law with running | 14.54 | 13.35
Power law | 48.70 | 42.89
Gaussian bump | 13.51 | 11.65
First order phase transition | 118.68 | 64.18
Table 1: Signal to noise ratio in TDI channel A for the four SGWB models, with
and without the presence of the galactic foreground as an additional noise
component. The galactic foreground here considered has an SNR of 1627.39

It is possible to notice that including the foreground as part of the noise
($S^{A}_{n}(f):=S^{A}_{GB}(f)+S^{A}_{n}(f)$) leads to a substantial decrease
of the SNR for the FOPT background but the SNR does not change very much for
the other models.

We plot in fig. 7 the value of the SNR in channel $A$ vs. the energy density
at 1mHz for the different models. The dotted lines assume no presence of the
foreground, whereas the continuous lines include the presence of the
foreground. As expected there is a direct correlation between increasing the
energy density and an increase in the SNR. Moreover, the presence of the
foreground mostly affects the SNR of the FOPT. In fact, in the presence of the
foreground the energy density must be two times larger to have the same SNR as
it would in the absence of a foreground.

Figure 7: Signal to noise ratio of four different SGWB signals in the TDI
channel A, Power law, Power law-with-running, Gaussian bump and First order
phase transition, versus the energy density at 1mHz, both considering
(continuous lines) or not (dashed lines) the presence foreground

##  III Results

###  III.1 Impact of instrumental noise knowledge uncertainty on SGWB recovery

Here, we explore how the measurement precision of the SGWB parameters changes
in the presence of instrumental noise knowledge uncertainty, for each of the
SGWB models described in Section II.6. We use the Fisher matrix formalism
described in Section II.2, which assumes that the noise is uncorrelated at
different frequencies. We assume we use three TDI channels in our analysis, A,
E and $\zeta$, as described in Section II.4. We model uncertainties in the PSD
and CSD at each frequency following the model described in Section II.3. To
build the Fisher matrix we need the following elements:

  1. 1.

The derivatives of the PSD and CSD at each frequency with respect to the
parameters of the SGWB model.

  2. 2.

The derivatives of the PSD and CSD at each frequency with respect to the
parameter (amplitude) of the Galactic Binaries. The addition of this parameter
extends the dimension of the Fisher matrix by one.

  3. 3.

The derivatives of the PSD and CSD at each frequency with respect to the
parameters of the instrumental noise model. The instrumental noise model is
based on 9 different splines: 3 splines to model the PSD of A, E and $\zeta$,
and 3 splines each for real and imaginary parts of the CSDs for AE, A$\zeta$
and E$\zeta$. Each spline has a number of parameters equal to the umber of
knots, which we take to be $13$. The total number of noise parameters is
therefore 9 x 13 = 117.

  4. 4.

The evaluation of the Fisher matrix from these elements using Eq. (9), which
is summed over frequency.

  5. 5.

The choice of a prior on the instrumental noise parameters. We use a Gaussian
prior, which is implemented in the Fisher matrix formalism by adding the prior
matrix to the Fisher matrix before computing its inverse (see Eq. (11)). For
this first study we take the priors on each noise parameter to be independent,
with zero mean and equal variance, $\sigma_{\rm inst}$. In this section we fix
$\sigma_{\rm inst}=1$, which means we are allowing for up to an order of
magnitude uncertainty in the instrumental noise at each frequency.

  6. 6.

We compute the inverse of the Fisher matrix after adding the prior to obtain
an estimate of the measurement uncertainty, from the square root of the
diagonal elements of the inverse as explained in Sec. II.2. We also compute
the inverse of the SGWB-parameter only sub-matrix of the Fisher matrix, which
represents the expected uncertainty in the absence of instrumental noise
uncertainties.

For each SGWB model we will present the results in two different ways.
Firstly, we will show the ratio of the uncertainties in the SGWB parameters in
the presence of instrumental noise uncertainties to those uncertainties when
perfect knowledge of the instrumental noise is assumed. These results
illustrate the impact of lack of noise knowledge on SGWB characterisation.
Secondly, we will show the actual uncertainties in the SGWB parameters, as
computed from the Fisher matrix. Of particular interest is the uncertainty in
the log-energy density of the background. As a rule of thumb, a background
will be detectable if the uncertainty $\Delta\ln(A)<1$ ( $\Delta\ln(A)=(\Delta
A)/A$). In both cases, we will plot results as a function of the background
amplitude/the background energy density at a reference frequency of 1mHz (the
logarithm of these quantities are linearly related, so they can be easily
represented using bottom/top axes in a single figure). For the second type of
plot, solid lines show results in the presence of noise knowledge uncertainty,
and dashed lines give results assuming perfect noise knowledge. In both the
analysis we consider the foreground amplitude to vary and we consider it as
additional source of noise together with the instrumental noise.

####  III.1.1 Power law

A power law SGWB is described by two parameters: the slope and the amplitude.
The full Fisher matrix, including instrumental noise and foreground
parameters, is $120\times 120$.

Figure 8:  Results for the power law SGWB model considering the foreground as
additional source of noise. The upper panel shows the ratio of the
uncertainties of the SGWB parameters (amplitude and slope) when including
instrumental noise uncertainties or assuming perfect noise knowledge. This
ratio is plotted versus the amplitude (bottom axis) and SGWB energy density at
1mHz (top axis). The lower panel shows the estimated parameter uncertainties
for the two cases. Once again this is as a function of amplitude/energy
density, but for a restricted range. The horizontal red dashed line
corresponds to an uncertainty of one, which is our threshold on the
uncertainty in log-energy density for deciding that a background is
detectable. The vertical red dashed line indicates the reference SGWB
amplitude given in Section II.6.

Figure 8 shows the results computed for this model. We see that in the
presence of instrumental noise uncertainties, the uncertainty in the SGWB
parameters increases by a factor of $\sim 55$–$60$, with the uncertainty in
the slope being slightly more affected than that of the amplitude. The
increase is lower for high background amplitudes, as expected, but only when
the background is one to two orders of magnitude brighter than the reference
value. Considering the raw uncertainties, we see that the uncertainty in the
log-energy density is typically a factor of $\sim 50$ larger and the
background energy density would have to to be a factor of $\sim 50$ times
higher to be characterised with the same measurement precision when there is
instrumental noise uncertainty as it could be without those uncertainties.
However, a background with amplitude equal to the reference value should
(just) be detectable even allowing for confusion with instrumental noise
mismodelling.

####  III.1.2 Power law with running

For the power law with running SGWB, the fisher matrix is $121\times 121$, as
the SGWB model depends on 3 parameters: slope, amplitude and running index,
$\alpha$.

Figure 9:  As Figure 8, but now for the power-law with running SGWB model

The results for this model are shown in fig. 9. In this case we see that the
uncertainties in the SGWB parameters increase by a factor of $\sim 30$–$75$,
with the uncertainty on the amplitude being most affected in this case. Once
again, the relative increase in the uncertainty is somewhat lower at higher
background amplitudes. The lower panel of fig. 9 shows that the background is
not detectable at the reference amplitude. An energy density $\sim 20$ times
higher would be required for a detection. In general, the background has to
have an energy density $\sim 60$ times higher to be characterised with the
same measurement precision when there is instrumental noise uncertainty as it
could be without those uncertainties and there is a similar increase in the
parameter measurement uncertainty at fixed background energy density.

####  III.1.3 Gaussian bump

As for the power law, the Fisher matrix is a $120\times 120$ as we have 2
signal parameters: the Gaussian width and the amplitude.

Figure 10:  As Figure 8, but now for the Gaussian bump SGWB model

The results for this model are shown in Figure 10. In this case, the
degradation in the precision of parameter measurement is a factor of $\sim
2$–$8$ when allowing for lack of knowledge of the instrumental noise. This
difference in behaviour is related to the different shapes of the SGWBs being
considered. A Gaussian is more distinct from the spline model being used to
represent the instrumental noise uncertainties than a power law, and hence the
degree of confusion between the two models is less in this case. From the
lower panel of Figure 10: we see that the energy density in a Gaussian bump
SGWB has to be a factor of $\sim 10$ times higher for it to be characterised
with the same measurement precision when there is instrumental noise
uncertainty as it could be in the absence of those uncertainties. A Gaussian
bump background at the reference amplitude would be detectable, and the width
of such a Gaussian could be measured to a few tens of percent precision. This
measurement precision improves approximately linearly with the background
energy density.

####  III.1.4 First order phase transition

The FOPT model is again characterised by two parameters, an amplitude and a
spectral index, and has a $120\times 120$ Fisher matrix. The results for this
model are shown in Figure 11. When allowing for instrumental noise knowledge
uncertainties, the precision with which the SGWB log-energy density can be
characterised degrades by a factor of $\sim 20$. The degradation in the
determination of the spectral index is even larger, $\sim 35$. Once again, to
achieve the same measurement precision, the background energy density would
have to be $\sim 20$ times larger than it would need to be in the absence of
noise knowledge uncertainties. Nonetheless, a FOPT background at the reference
amplitude would still be detectable and provide a measurement of the spectral
index at the level of $\sim\pm 0.8$.

Figure 11:  As Figure 8, but now for the first order phase transition SGWB
model

The previous results were computed considering the presence of the galactic
foreground. In the appendix C we report similar results, computed without
taking into consideration the Galactic foreground. Redoing these analyses
ignoring the foreground we do not see big differences in the uncertainty
ratio, nor in the absolute uncertainties, when these are compared at fixed
SNR, i.e., when the signal-to-noise ratio is recomputed without the galactic
binaries included in the spectral density. To illustrate this, we show in
Fig.12 the precision of the measurement of the log-energy density of the
background, as a function of the SNR in TDI channel A, for all SGWB models and
both including and not including the galactic binary foreground. We see that
the uncertainty is typically larger when the foreground is present, but this
is typically less than a factor of a few. The Gaussian bump and power law
background are most affected, with the uncertainty at fixed SNR and the SNR
required for detection both decreasing by a factor of a few when the galactic
binary background is removed from the spectral density. For the Gaussian bump,
the uncertainty decreases by a factor of a little more than two when the
Galactic background is excluded, and the SNR needed to reach the
$\Delta\ln(A)<1$ threshold for detection decreases by a similar factor. For
the power law, the uncertainty decreases by about a factor of $4$, and the
$\Delta\ln(A)<1$ threshold required for detection is reached at an SNR that is
a factor of $\sim 4$ smaller. For the power law with running and the FOPT
backgrounds, the uncertainty at fixed SNR is almost unchanged, and the
threshold SNR for detection is within a factor of $1.5$ and $2$, respectively.

This behaviour can be understood by looking at the shapes of the various SGWBs
in Figure 6. Figure 7 demonstrates that the removal of the foreground does not
affect the SNR very much. The only SGWB that shows a significant change is the
FOPT, for which most of the power is at frequencies where the foreground is
significant. However, in the region around $300\mu$Hz, where the majority of
the SNR is generated, the shape of the FOPT is very different to the
foreground. This is also true for the power-law-with running model around
$5$mHz, where the majority of its SNR is generated. The power-law model, on
the other hand, is quite parallel to the foreground at low frequency, and the
Gaussian bump is quite parallel to the foreground at a few mHz. This most
likely explains why the latter two backgrounds are more difficult to
distinguish from a galactic foreground, and therefore more affected by its
inclusion.

Figure 12: Signal to noise ratio of four different SGWB signals in the TDI
channel A, Power law, Power law-with-running, Gaussian bump and First order
phase transition, versus the error in log amplitude both considering
(continuous lines) or not (dashed lines) the presence foreground

###  III.2 Setting a noise knowledge requirement

In this section we will now explore how the amount of uncertainty in the
instrumental noise impacts the results. In practice we will not be completely
ignorant of the instrumental noise. Measurements on-board the satellites will
provide an indication of the size of certain noise components. In principle,
it might therefore be possible to place a requirement on how well the
instrumental noise must be known in order to not degrade the science output of
the mission. To assess this, we will recompute the results while changing the
variance of the Gaussian prior on the instrumental noise spline parameters. We
will vary the prior on the spline weights from very small values
($\log_{10}(\sigma_{inst})=-10$), representing near-perfect knowledge of the
noise, to very high values $\log_{10}(\sigma_{inst})=6$, representing no
knowledge of the noise.

We fix the amplitude of the background for each SGWB model, so that it
corresponds to an SNR of $\sim 120-140$ in each case 333Specifically the SNR
of the power law with running is 135, the SNR of the power law is 138, the SNR
of the FOPT is 142 and the SNR of the Gaussian bump is 120.. This choice of
SNR was motivated by Figure 12, which shows that an SNR greater than 100 is
required in order to ensure all types of SGWB are detectable. For the case of
the power law we also show results with the amplitude set to the reference
energy density, which shows that the exact choice of background amplitude does
not make a significant difference to the qualitative behaviour, only to the
absolute value of the uncertainty.

For all SGWB models we again present the results in two different ways — as a
ratio of the SGWB parameter measurement uncertainties when instrumental noise
uncertainties are considered to those assuming perfect noise knowledge, and as
the absolute measurement uncertainty. Results for the power law model are
shown in Figure 13, for the Power-law-with-running model in Figure 14, for the
Gaussian bump model in Figure 15 and for the FOPT in Figure 16.

Figure 13:  As Figure 8, but now for fixed background amplitude and varying
the variance of the Gaussian prior on the instrumental noise spline model.
This plot is for a power law background, and the amplitude has been fixed such
that the SNR in TDI channel A is 138 (continous line) and 43.7 (dashed lines).

The results for all four SGWB models are qualitatively similar. For very low
prior uncertainties the ratio of the uncertainties tends to unity. This is
expected as this limit corresponds to the limit in which the instrumental
noise is perfectly known. As the prior uncertainty is increased beyond $\sim
10^{-6}$ the measurement precision in the presence of noise knowledge
uncertainties starts to increase. When the noise knowledge uncertainty reaches
$\sim 10^{-2}$ for power law and Gaussian bump, $\sim 10^{-1}$ for power law-
with-running and 10 for FOPT, the measurement precision ratio saturates. This
final value reflects the expected uncertainty in the absence of any noise
knowledge. The results given in Section III.1 were all computed in this
regime.

The main conclusion from these results is that if we wanted to ensure that
there was no degradation in LISA science due to lack of noise knowledge, the
necessary requirement on the noise knowledge would be $<<10\%$. In the LISA
Pathfinder mission, which was designed to accurately characterize the free-
fall performance of test masses in a space-based environment, the observed
noise could only be explained within some margin: the physical origin of the
measured sub-mHz acceleration is only partially understood, as more than 50%
of its PSD is still unmodeled LorePhD ; lpf_noiseperf3 .

It is therefore unrealistic to expect that a noise requirement at the $\sim
1$–$10\%$ level could be met. At noise uncertainties above this threshold,
there is little difference between some and no noise-knowledge, at least
within the model for instrumental noise variations considered here. We
conclude that no useful and achievable noise knowledge requirement could be
implemented in practice.

While we will not be able to achieve the precision that would be possible
under ideal circumstances, it is important to emphasise that this does not
mean we will not be able to detect and characterise modelled SGWBs. In all
cases, at SNR of $\sim>100$, the amplitude can be constrained to a few tens of
percent, even without any knowledge of the instrumental noise.

Figure 14:  As Figure fig. 13 but now for the power-law-with-running model.
The background amplitude has been fixed to give an overall SNR of $135$.

Figure 15:  As Figure fig. 13 but now for the Gaussian bump model. The
background amplitude has been fixed to give an overall SNR of $122$.

Figure 16:  As Figure fig. 13 but now for the FOPT model. The background
amplitude has been fixed to give an overall SNR of $142$.

###  III.3 Signal reconstruction

To finish this section we will use our Fisher matrix results to illustrate how
well we can reconstruct the Power law, the foreground, and the instrumental
noise. To do this, we will approximate the posterior distribution on the model
parameters using a multi-variate Gaussian with covariance matrix equal to the
inverse of the Fisher matrix. We can then take random draws from this
fake444Even if the measurements $x$ are Gaussian the posterior is not because
is a function of $\theta$. A Gaussian with unknown variance $\sigma$ is not a
Gaussian on that variance,
$p(x|\theta)\propto\frac{1}{\sigma}\exp(-\frac{x^{2}}{2\sigma^{2}})$, is a
Gaussian in $x$, but not in $\sigma$ posterior distribution and plot the PSD
of the SGWB, the foreground and instrumental noise corresponding to the drawn
parameters. In fig. 17 we follow this procedure for a power law signal with an
SNR of 872. The choice of the SNR is driven to avoid the breakdown of the
Fisher matrix approximation for the SGWB parameters, because the SGWB
parameter uncertainties are large and no longer in the linear signal regime as
explained and, shown in details, in Appendix C.1.

The three panels show the reconstructed ASDs for the Power law, for the
foreground and for the instrumental noise and the total, which is the sum of
the three.

What we would expect is that our ability to measure the total spectral density
is roughly independent of the relative amplitudes of the two components, since
this is what we actually see and measure in the data. Our model attempts to
split that measurement into constituent components. If one of those components
is much weaker than the other we would not expect to recover it as well as
when the components are making comparable contributions to the data. Top panel
of Fig. 17 (and the Figure 26 in the appendix C.1) is consistent with this
expectation. We see indeed in Fig. 17 that the galactic binaries are well
recovered. Moreover, the only noise reconstruction for the TDI A and E suffers
from the presence of the GW background and foreground in the regime $0.4mHz$ –
$4mHz$ where these two signals have the majority of power. As a final point,
it is clear that the SGWB is being best constrained around a frequency of 4mHz
where the power of the GB is less and the uncertainty in the instrumental
noise is also largest (although still small) at this point. This can be
understood from Figure 6, which shows that the power law background is closest
to the instrumental noise ASD at that frequency and the Galactic binaries pick
at 1mHz, and so this frequency range dominates the SNR in the signal. We
expect to be able to measure the background best in the frequency range where
it is most dominant relative to the instrumental noise and distinguishable
from the galactic binaries.

Figure 17:  We show the Power law signal, Galactic binaries and noise ASDs
corresponding to random draws from the posterior, approximated using the
Fisher matrix as described in the text. In each panel the curves correspond to
the three TDI channels: A (blue), E (red) and $\zeta$ (green). Upper first
panel: reconstructed SGWB; upper second panel: reconstructed TM and OMS
instrumental noise; middle panel: reconstructed Galactic binaries; lower
panel: total reconstructed ASD (signal + noise + GB).

##  IV Discussion and conclusion

We have explored the impact of noise knowledge uncertainty on measuring the
parameters of various modelled stochastic gravitational wave backgrounds. This
was done by modelling instrumental noise uncertainties using cubic splines to
represent deviations away from the design PSDs and CSDs for the three TDI
channels $A$, $E$ and $\zeta$. We then used a Fisher matrix analysis to
evaluate the expected uncertainties in the measurements of the model
parameters when fitting a model including the instrumental noise uncertainties
and compared it to fitting a model without those uncertainties. The degree of
uncertainty was characterized by including a Gaussian prior on the
instrumental noise parameters, allowing us to quantify the impact of imposing
a requirement on our noise knowledge.

This analysis showed that, for all SGWB models, allowing for instrumental
noise uncertainties leads to a significant increase in the uncertainty in our
measurements of the background parameters. The increase in uncertainty was a
factor of $2-8$ for the Gaussian bump model, which reduces to $2-4$ when not
including GB as foreground, $55-60$ for the power law ($15-30$ without GB as
foreground), $20-35$ for the First order phase transition ($20-50$ without GB
as foreground) and $30-75$ for the power law with running ($20-75$ without GB
as foreground). These increased uncertainties correspond to the threshold
background energy density required for detection increasing by a factor of
$10$ ($5$ without GB) for the Gaussian bump model, a factor of $60$ with and
without GB for the Power-with-running and a factor of $20$ with and without GB
for all other models ($50$ for the power law when including GBs). The
threshold energy density at 1mHz at which the backgrounds start to be
detectable are $4\times 10^{-13}$, $4\times 10^{-13}$, $2\times 10^{-16}$ and
$10^{-12}$ ($10^{-13}$, $2.5\times 10^{-13}$, $8\times 10^{-17}$ and $5\times
10^{-13}$ if we do not include the GB foreground) for the power law, power law
with running, Gaussian bump and FOPT models respectively. Comparing these to
the reference background amplitudes introduced in Section II.6 we see that the
power law, the Gaussian bump and FOPT backgrounds are detectable at the
reference amplitudes, while the power law with running is not as the threshold
energy density is $1.5$ times higher than the reference. However, for this
latter background the amplitudes were specified based on the SNR and not on a
physical model. The reference amplitudes for the power law and FOPT
backgrounds are based on physical model predictions, so it is more important
that these backgrounds are still detectable. We note that this result does
depend on the particular choice of model we used for representing the
instrumental noise uncertainty. If this model is made even more flexible, for
example by increasing the number of spline knots used, the threshold would
increase further and potentially also make the reference backgrounds
undetectable.

When we vary our assumed level of knowledge of the instrumental noise, we find
that the uncertainties on the SGWB parameters show a similar trend for all
models, starting to degrade at relative small uncertainties, increasing and
then saturating after a certain point. The point at which the sensitivity
starts to degrade is when the uncertainty in the log-spectral density of the
noise reaches $\sim 10^{-6}$–$10^{-5}$, depending on the SGWB model. The
uncertainty saturates at log-spectral density uncertainties of $\sim
10^{-2}/10^{-1}/10^{-2}/10^{1}$ ($10^{-2}/10^{-1}/10^{-3}/10^{-1}$ in the
absence of the GB foreground) for the power law/power-law with
running/Gaussian bump/FOPT backgrounds respectively. This means that if we
wanted to limit the degradation in the science that arises from lack of noise
knowledge we would have to impose a very stringent requirement on our
knowledge of the noise. This is likely to be impossible to implement in
practice, so we will have to accept that our ability to resolve SGWBs will not
be as good as calculations that assume perfect noise knowledge predict.

It is important to note that these results are based on some assumptions which
might not hold in practice. In particular, we have considered only modelled
SGWBs and we have assumed a particular form for variations in the PSD that
forces variations to be smoothly varying as a function of frequency. If the
number of knots is increased to obtain a more flexible instrumental noise
model, with potentially faster variations of the PSD as a function of
frequency, we already see a degradation of 2 or 3 orders of magnitude in the
estimation of the log-energy density.
It is the distinguishability of the models that allows us to measure the
parameters of the SGWBs. In the extreme picture where we do not want to make
any assumption at all about the form of the instrumental or SGWB spectral
densities, then spectral separation will not be possible. We will be able to
report measured power spectral density in all channels, and cross-spectral
densities between them, and translate these into upper limits on the SGWB
amplitude; but any interpretation of this as an actual detection will require
independent confirmation from another detector Muratore:2022nbh .

All previous studies of the separation of instrumental noise and stochastic
backgrounds have required assumptions: in Baghi:2023qnq  it was assumed that
the instrumental uncertainty is a spline and that the SGWB has a power law
spectrum; in Hartwig:2023pft  it was assumed that the instrumental noise is
determined by 12 individual noise levels; and in this paper we are assuming
something similar to Baghi:2023qnq , although with a bit more flexibility, a
wider variety of SGWB models and a different noise model for single satellite
links. The SGWBinner Caprini_2019  is agnostic on the spectrum of the
background but it can only work because it assumes a specific model for the
instrumental noise. That is not going to be possible in practice. SGWBinner
could be adapted to use a more flexible noise model, similar to the model used
here, but the precision on the background recovery will be degraded. If we
have a completely general instrumental noise model and a completely general
SGWB model then we won’t be able to separate them. In that case, the only hope
would be that the SGWB is above the design sensitivity and we trust that the
instrumental noise meets the mission requirements, in which case the best
interpretation of such an observation would be a SGWB. However, even then an
assumption would be made that the mission had met the design sensitivity
requirements. An exploration of how our ability to separate instrumental noise
as the spectral models of the SGWB and the instrumental noise are made more
complicated should be the focus of future work.

##  V Acknowledgement

M.M gratefully acknowledges the support of the German Space Agency, DLR. The
work is supported by the Federal Ministry for Economic Affairs and Climate
Action based on a resolution of the German Bundestag (Project Ref. No. FKZ 50
OQ 2301). The authors thank Chiara Caprini for her review and suggestions
regarding the use of SGWB templates in the text. Additionally, we acknowledge
Mauro Pieroni for providing feedback on the presented results and the
implementation of spline models. Furthermore, we thank Olaf Hartwig for
discussions on the accurate computation of the GW transfer function and on the
noise models and Marc Lilley for comparisons on the results obtained with the
Fisher matrix formalism. We thank Lorenzo Sala for feedbacks on the recent
analysis on the LISA Pathfinder performance results. The scientific
discussions with Quentin Baghi, Jean Baptiste Bailey, Germano Nardini,
Nikolaos Karnesis, Jesus Torrado, Nam Dam Quam, Henry Hinshauspe, Antoine
Petiteau, and Riccardo Buscicchio have also been highly appreciated. Their
feedback on the results and methodology presented in the paper has been
valuable for the final outcome

##  Appendix A Likelihood derivation

We derive the likelihood starting from the noise properties and explain why it
takes the form shown in Sec. II. If we assume that the real time series $n(t)$
is a stationary, zero-mean, Gaussian and ergodic random process, then the
Fourier transform of the noise $\tilde{n}_{k}=\tilde{n}(f_{k})$ at each
frequency $f_{k}$ is normally distributed with zero-mean and variance
$\sigma_{k}^{2}$. Thus, the natural log-likelihood at each frequency takes the
form of a two-dimensional normal distribution

$\ln
p(\tilde{n}_{k})=-\frac{1}{2}\Big{(}\frac{\Re[\tilde{n}_{k}]^{2}}{\sigma^{2}_{k}}+\frac{\Im[\tilde{n}_{k}]^{2}}{\sigma^{2}_{k}}\Big{)}-\frac{1}{2}\ln[(2\pi\sigma_{k}^{2})^{2}]$
(40)

where we assumed that the real $\Re$ and imaginary $\Im$ part of the noise are
not correlated and have the same variance. If we further assume that the
variance of the noise at different frequencies follow a one-sided power
spectral density $S_{n}(f)$ then :

$\displaystyle<\tilde{n}^{*}(f^{\prime})\tilde{n}(f)>$
$\displaystyle=\frac{1}{2}S_{n}(f)\delta(f-f^{\prime})$ (41)
$\displaystyle<\tilde{n}(f^{\prime})\tilde{n}(f)>$ $\displaystyle=0$ (42)

where we used the expectation value of $<>$ over the data generating process.
For a set of frequencies the first relation can be written as

$\displaystyle<\tilde{n}^{*}_{k}\tilde{n}_{j}>$
$\displaystyle=\frac{T}{2}S_{n}(f_{k})\delta_{jk}$ (43)

Therefore the variance of the real and imaginary part of the noise is given by

$<\Re[\tilde{n}_{k}]^{2}>=<\Im[\tilde{n}_{k}]^{2}>=\frac{T}{4}S_{n}(f_{k})=\sigma_{k}^{2}\,.$
(44)

We can write the natural log-likelihood for all the measured frequencies as:

$\sum_{k}\ln
p(\tilde{n}_{k})=-\sum_{k=1}^{n}\ln[2\pi\frac{T}{4}S_{n}(f_{k})]-\frac{1}{2}\sum_{k=1}^{n}\frac{|\tilde{n}(f)|^{2}}{\frac{T}{4}S_{n}(f_{k})}$
(45)

which becomes in the continuum limit:

$\ln
p(\tilde{n})=-\ln\bigg{[}2\pi\frac{1}{4}\det[S_{n}(f)\delta(f-f^{\prime})]\bigg{]}-\frac{1}{2}4\int_{0}^{\infty}\frac{|\tilde{n}(f)|^{2}}{S_{n}(f)}\rm{d}f$
(46)

Note that in the continuum limit the variance of the noise can be thought as
an operator. In fact, one can define the inner product

$(a(t)|b(t))=4\Re\int_{0}^{\infty}\int_{0}^{\infty}\tilde{a}^{*}(f)\,\Sigma^{-1}(f,f^{\prime})\tilde{b}(f^{\prime})\,\rm{d}f\,\rm{d}f^{\prime}$
(47)

with $\Sigma^{-1}$ defined through the relation:

$\int_{0}^{\infty}\Sigma^{-1}(f,f^{\prime})\Sigma(f^{\prime},f^{\prime\prime})\rm{d}f^{\prime}=\delta(f-f^{\prime\prime})\,,$
(48)

where if we set in Eq.48 that
$\Sigma(f^{\prime},f^{\prime\prime})=\delta(f^{\prime}-f^{\prime\prime})S_{n}(f^{\prime})$
we obtain

$\Sigma^{-1}(f,f^{\prime\prime})S_{n}(f^{\prime\prime})=\delta(f-f^{\prime\prime})$
(49)

and the inner product becomes

$\displaystyle(a(t)|b(t))$
$\displaystyle=4\Re\int_{0}^{\infty}\int_{0}^{\infty}\frac{\tilde{a}^{*}(f)\,\delta(f-f^{\prime})\tilde{b}(f^{\prime})}{S_{n}(f^{\prime})}\,\rm{d}f\,\rm{d}f^{\prime}$
$\displaystyle=4\Re\int_{0}^{\infty}\frac{\tilde{a}^{*}(f)\tilde{b}(f)}{S_{n}(f)}\,\rm{d}f$
(50)

##  Appendix B Fisher matrix derivation

###  B.1 Single detector

To compute the Fisher matrix of Eq. 8 we need the first derivative of the log-
likelihood $l$ with respect to the parameters of the power spectral density.

Here we present the derivation of the Fisher matrix for the noise parameters
$\vec{\lambda}$ affecting the one-sided spectral density
$S_{n}(f|\vec{\lambda})$, but this can be easily extended also including the
gravitational wave background parameters $S_{n}(f|\vec{\lambda})\rightarrow
S_{n}(f|\vec{\lambda})+S_{\rm GW}(f|\vec{\theta})$. We differentiate the log-
likelihood of equation 45 with respect the parameters $\vec{\lambda}$:

$\displaystyle\frac{\partial l}{\partial\lambda^{i}}=$
$\displaystyle\sum_{k=1}^{n}[-\frac{1}{S_{n}(f_{k})}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{i}}+\frac{1}{2}\frac{|\tilde{n}(f_{k})|^{2}}{\frac{T}{4}S_{n}(f_{k})^{2}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{i}}],$ (51)

where we have omitted the dependency from $\vec{\lambda}$ to have a lighter
notation. The second derivative of the likelihood is then

$\displaystyle\frac{\partial^{2}l}{\partial\lambda^{i}\partial\lambda^{j}}=$
$\displaystyle\sum_{k=1}^{n}[\frac{1}{S^{2}_{n}(f_{k})}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{i}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{j}}$
$\displaystyle-\frac{1}{S_{n}(f_{k})}\frac{\partial^{2}S_{n}(f_{k})}{\partial\lambda^{i}\partial\lambda^{j}}$
$\displaystyle-\frac{1}{2}\frac{2|\tilde{n}(f_{k})|^{2}}{\frac{T}{4}S_{n}(f_{k})^{3}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{j}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{i}}$
$\displaystyle+\frac{1}{2}\frac{|\tilde{n}(f_{k})|^{2}}{\frac{T}{4}S_{n}(f_{k})^{2}}\frac{\partial^{2}S_{n}(f_{k})}{\partial\lambda^{j}\partial\lambda^{i}}]$
(52)

using the definition of eq. 43 the second and last term cancels and we get:

$\displaystyle\Gamma_{ij}$
$\displaystyle=\sum_{k=1}^{n}\frac{1}{S_{n}(f_{k})^{2}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{i}}\frac{\partial
S_{n}(f_{k})}{\partial\lambda^{j}}\,.$ (53)

If we want to get the continuum limit we need to recast a factor of
$T\rm{d}f$:

$\Gamma_{ij}=T\int_{0}^{\infty}\frac{1}{S_{n}(f)^{2}}\frac{\partial
S_{n}(f)}{\partial\lambda^{i}}\frac{\partial
S_{n}(f)}{\partial\lambda^{j}}\rm{d}f$ (54)

###  B.2 Multiple detectors: real and imaginary part as separate random
variables

If we want to generalize our derivation to multiple detectors or channels we
need to define the noise properties of each channel. For simplicity let us
consider two channels $A$ and $E$ with 4 independent random variables
$X(f_{k})=X_{k}=\\{\Re[\tilde{X}_{k}^{A}],\Im[\tilde{X}_{k}^{A}],\Re[\tilde{X}_{k}^{E}],\Im[\tilde{X}_{k}^{E}]\\}$
at each frequency $f_{k}$. Since the final likelihood will be given by the
product over all the frequencies, we consider only one frequency and we drop
the subscript ”k”. We can specify the spectral densities of each channel and
the cross-spectral densities with:

$\displaystyle<\tilde{X}^{c*}(f^{\prime})\tilde{X}^{c}(f)>$
$\displaystyle=\frac{1}{2}S_{c}(f)\delta(f-f^{\prime})$ (55a)
$\displaystyle<\tilde{X}^{c}(f^{\prime})\tilde{X}^{c}(f)>$ $\displaystyle=0$
(55b) $\displaystyle<\tilde{X}^{E*}(f^{\prime})\tilde{X}^{A}(f)>$
$\displaystyle=\frac{1}{2}S^{*}_{AE}(f)\delta(f-f^{\prime})$ (55c)
$\displaystyle<\tilde{X}^{A*}(f^{\prime})\tilde{X}^{E}(f)>$
$\displaystyle=\frac{1}{2}S_{AE}(f)\delta(f-f^{\prime})$ (55d)
$\displaystyle<\tilde{X}^{A}(f^{\prime})\tilde{X}^{E}(f)>$ $\displaystyle=0$
(55e)

where the first two rows are valid for both channels $c=A,E$, and $S_{c}$ is
real and $S_{AE}$ is complex. From the above expression we can deduce:

$\displaystyle<\Re[\tilde{X}_{c}]\Re[\tilde{X}_{c}]>+<\Im[\tilde{X}_{c}]\Im[\tilde{X}_{c}]>$
$\displaystyle=\frac{S_{c}}{2}$ (56a)
$\displaystyle<\Re[\tilde{X}_{c}]\Re[\tilde{X}_{c}]-\Im[\tilde{X}_{c}]\Im[\tilde{X}_{c}]>$
$\displaystyle=0$ (56b) $\displaystyle<\Re[\tilde{X}_{c}]\Im[\tilde{X}_{c}]>$
$\displaystyle=0$ (56c)
$\displaystyle<\Re[\tilde{X}_{A}]\Re[\tilde{X}_{E}]>+<\Im[\tilde{X}_{A}]\Im[\tilde{X}_{E}]>$
$\displaystyle=\frac{\Re[S_{AE}]}{2}$ (56d)
$\displaystyle<\Re[\tilde{X}_{A}]\Re[\tilde{X}_{E}]>-<\Im[\tilde{X}_{A}]\Im[\tilde{X}_{E}]>$
$\displaystyle=0$ (56e)
$\displaystyle<\Re[\tilde{X}_{A}]\Im[\tilde{X}_{E}]>-<\Im[\tilde{X}_{A}]\Re[\tilde{X}_{E}]>$
$\displaystyle=\frac{\Im[S_{AE}]}{2}$ (56f)
$\displaystyle<\Re[\tilde{X}_{A}]\Im[\tilde{X}_{E}]>+<\Im[\tilde{X}_{A}]\Re[\tilde{X}_{E}]>$
$\displaystyle=0\,,$ (56g)

where in the last four rows $c=A,E$. Note that these are in total 10
independent conditions (3 equations for A, 3 equations for E and 4 equations
for AE) that specify uniquely the 10 independent elements of a symmetric
covariance matrix.
For a single frequency we can generalize the likelihood to two channels as:

$p(X)=\frac{1}{\sqrt{(2\pi)^{2\times
N_{c}}\det(\Sigma)}}e^{-\frac{1}{2}X^{T}\Sigma^{-1}X}$ (57)

where $N_{c}=2$ are the number of channels, X is a quadri-dimensional vector
define above, $\Sigma$ is the multiple-channels covariance matrix:

$\Sigma=\left(\begin{array}[]{cccc}\frac{{S_{A}}}{4}&0&\frac{\Re({S_{AE}})}{4}&\frac{\Im({S_{AE}})}{4}\\\
0&\frac{{S_{A}}}{4}&\frac{\Im({S_{AE}})}{4}&\frac{\Re({S_{AE}})}{4}\\\
\frac{\Re({S_{AE}})}{4}&\frac{\Im({S_{AE}})}{4}&\frac{{S_{E}}}{4}&0\\\
\frac{\Im({S_{AE}})}{4}&\frac{\Im({S_{AE}})}{4}&0&\frac{{S_{E}}}{4}\\\
\end{array}\right)$ (58)

where here each element is evaluated at the fixed frequency. It can be shown
that the expectation value of $X^{T}\Sigma^{-1}X$ equals the degrees of
freedom, in this case 4. We have two channels, where each one has two degrees
of freedom associated with the real and imaginary part of $\tilde{X}$.

We can then derive the Fisher matrix for the multiple channel case. Taking the
first derivative of the log-likelihood

$\displaystyle\frac{\partial\ln p(X)}{\partial\lambda^{i}}$
$\displaystyle=-\frac{1}{2}\frac{1}{\det(\Sigma)}\frac{\partial\det(\Sigma)}{\partial\lambda^{i}}-\frac{1}{2}X^{T}\frac{\partial\Sigma^{-1}}{\partial\lambda^{i}}X$
(59)

where we can use the following property of the determinant:

$\displaystyle\frac{\partial\det(\Sigma)}{\partial\lambda^{i}}$
$\displaystyle=\det(\Sigma)\;\rm{Tr}[\Sigma^{-1}\frac{\partial\Sigma}{\partial\lambda^{i}}]$
$\displaystyle=\det(\Sigma)\;\Big{[}\Sigma^{-1}\Big{]}_{lm}\Big{[}\frac{\partial\Sigma}{\partial\lambda^{i}}\Big{]}^{ml}$
(60)

to obtain

$\displaystyle\frac{\partial\ln p(X)}{\partial\lambda^{i}}$
$\displaystyle=-\frac{1}{2}\Big{[}\Sigma^{-1}\Big{]}_{lm}\Big{[}\frac{\partial\Sigma}{\partial\lambda^{i}}\Big{]}^{ml}-\frac{1}{2}X^{T}\frac{\partial\Sigma^{-1}}{\partial\lambda^{i}}X\,.$
(61)

Then, the second derivative of the log-likelihood takes the form:

$\displaystyle\frac{\partial^{2}\ln
p(X)}{\partial\lambda^{i}\partial\lambda^{j}}$
$\displaystyle=-\frac{1}{2}\frac{\partial(\Sigma^{-1})^{lm}}{\partial\lambda^{i}}\frac{\partial\Sigma_{ml}}{\partial\lambda^{j}}$
$\displaystyle-\frac{1}{2}\Sigma_{lm}^{-1}\frac{\partial^{2}\Sigma^{ml}}{\partial\lambda^{i}\partial\lambda^{j}}-\frac{1}{2}X^{T}\frac{\partial^{2}\Sigma^{-1}}{\partial\lambda^{i}\lambda^{j}}X$
(62)

.

We can finally compute the Fisher matrix for a single frequency with:

$\displaystyle\Gamma_{ij}$
$\displaystyle=\frac{1}{2}\Big{[}\frac{\partial\Sigma_{lm}^{-1}}{\partial\lambda^{i}}\frac{\partial\Sigma^{ml}}{\partial\lambda^{j}}+\Sigma_{lm}^{-1}\frac{\partial^{2}\Sigma^{ml}}{\partial\lambda^{i}\partial\lambda^{j}}+\Sigma_{ml}\frac{\partial^{2}(\Sigma^{-1})^{lm}}{\partial\lambda^{i}\lambda^{j}}\Big{]},$
(63)

where we have considered
$<X_{l}^{T}\frac{\partial^{2}\Sigma_{lm}^{-1}}{\partial\lambda^{i}\lambda^{j}}X_{m}>=\frac{\partial^{2}\Sigma_{lm}^{-1}}{\partial\lambda^{i}\lambda^{j}}\Sigma_{ml}$.
If we use the property

$\frac{\partial(\Sigma)^{-1}_{lm}}{\partial\lambda}=-(\Sigma)^{-1}_{ln}\frac{\partial(\Sigma)_{nq}}{\partial\lambda}(\Sigma)^{-1}_{qm}$
(64)

we obtain the following expression:

$\displaystyle\Gamma_{ij}=$
$\displaystyle\frac{1}{2}\rm{Tr}\Big{[}-\Sigma^{-1}\frac{\partial\Sigma}{\partial\lambda^{i}}\Sigma^{-1}\frac{\partial\Sigma}{\partial\lambda^{j}}+\Sigma^{-1}\frac{\partial^{2}\Sigma}{\partial\lambda^{i}\partial\lambda^{j}}+\Sigma\frac{\partial^{2}\Sigma^{-1}}{\partial\lambda^{i}\partial\lambda^{j}}\Big{]}$
(65)

which can be further simplified if we use the following properties:

$\displaystyle\partial(\Sigma\Sigma^{-1})=0$ (66a)
$\displaystyle\partial\Sigma\Sigma^{-1}+\Sigma\partial\Sigma^{-1}=0$ (66b)
$\displaystyle\partial^{2}\Sigma\Sigma^{-1}+2\partial\Sigma\partial\Sigma^{-1}+\Sigma\partial^{2}\Sigma^{-1}=0$
(66c)
$\displaystyle\partial^{2}\Sigma\Sigma^{-1}-2\partial\Sigma\,\Sigma^{-1}\,\partial\Sigma\,\Sigma^{-1}+\Sigma\partial^{2}\Sigma^{-1}=0$
(66d)

Note that in the above expression the first and last term of Eq.66d correspond
to the last two terms in the Fisher matrix expression (65). The final
expression for all frequencies can be easily obtained by taking the sum over
all the frequencies:

$\Gamma_{ij}=\frac{1}{2}\sum_{k=1}^{n}[(\Sigma_{k}^{-1})_{lr}\frac{\partial\Sigma_{k}^{rp}}{\partial\lambda^{i}}(\Sigma_{k}^{-1})_{pm}\frac{\partial\Sigma_{k}^{ml}}{\partial\lambda^{j}}]\,.$
(67)

Note that there is an additional factor of $1/2$ with respect to Eq. 53. If we
insert only the first 2 columns and rows of $\Sigma$ we obtain the previous
equation for the single channel as expected.

###  B.3 Multiple detectors: complex random variables

Equivalently, the likelihood can be written in terms of complex variables
$\tilde{X}_{A}$ and $\tilde{X}_{E}$ (complex_rnd_distr, )

$p(\tilde{X}_{A},\tilde{X}_{E})=\frac{e^{-[\tilde{X}_{A},\tilde{X}_{E}]^{\rm
H}\Sigma^{-1}[\tilde{X}_{A},\tilde{X}_{E}]}}{(2\pi)^{N_{c}}\det(\Sigma_{k})}$
(68)

where ”H” indicates the Hermitian conjugate, the factor of $1/2$ disappeared
because it is a complex distribution and must match Eqs (57)–(58), and the new
complex covariance matrix is defined as

$\Sigma=\frac{1}{2}\begin{pmatrix}S_{A}&S_{AE}\\\
S^{*}_{AE}&S_{E}\end{pmatrix}\,,$ (69)

where $\Sigma$ is now an Hermitian matrix and can be obtained from the
conditions imposed in Eqs. 55. The expectation value of
$[\tilde{X}_{A},\tilde{X}_{E}]^{\rm
H}\Sigma^{-1}[\tilde{X}_{A},\tilde{X}_{E}]$ over complex variable realizations
$[\tilde{X}_{A},\tilde{X}_{E}]$ is now 2. However, since the exponential does
not have any factor of $1/2$ for a complex distribution, we recover the same
number of degrees of freedom in the argument of the exponent as in the
previous derivation, i.e. we got $\exp{[\frac{1}{2}\,4]}$ for the case of
multiple detectors with real and imaginary part as separate random variables’
and $\exp{[2]}$ for the case considered here.

The derivation of the Fisher matrix differs from the previous one (Eq. 67)
only by the factor $1/2$:

$\Gamma_{ij}=\sum_{k=1}^{n}\left[(\Sigma_{k}^{-1})_{lr}\frac{\partial\Sigma_{k}^{rp}}{\partial\lambda^{i}}(\Sigma_{k}^{-1})_{pm}\frac{\partial\Sigma_{k}^{ml}}{\partial\lambda^{j}}\right],$
(70)

where the matrix $\Sigma_{k}$ is given by $\Sigma$ with spectral densities
evaluated at given frequency $f_{k}$.
Note that we can recover the single channel realization by using the first
element of $\Sigma$.
The continuum limit of the Fisher matrix in this formulation is given by

$\Gamma_{ij}=T\int_{0}^{\infty}(\Sigma_{k}^{-1})_{lr}\frac{\partial\Sigma_{k}^{rp}}{\partial\lambda^{i}}(\Sigma_{k}^{-1})_{pm}\frac{\partial\Sigma_{k}^{ml}}{\partial\lambda^{j}}\,\rm{d}f\,.$
(71)

###  B.4 Deterministic sources and noise cross-correlation

In the presence of a deterministic source, $\tilde{h}(f_{k}|\vec{\mu})$, the
derivative of the log-likelihood in Eq. (7) with respect to the source
parameters, $\vec{\mu}$, is:

$\frac{\partial
l}{\partial\mu^{i}}=-\sum_{k=1}^{n}\frac{|\tilde{s}(f_{k})-\tilde{h}(f_{k}|\vec{\mu})|}{\frac{T}{4}S_{n}(f_{k}|\vec{\lambda})}\frac{\partial\tilde{h}(f_{k}|\vec{\mu})}{\partial\mu^{i}}\,,$
(72)

and the derivative with respect to the parameters characterising the spectral
density, $\vec{\lambda}$, is:

$\displaystyle\frac{\partial l}{\partial\lambda^{i}}=$
$\displaystyle\sum_{k=1}^{n}[-\frac{1}{2}\frac{T}{\det(S_{n}(f_{k}|\vec{\lambda}))}+$
$\displaystyle+\frac{|\tilde{s}(f_{k})-\tilde{h}(f_{k}|\vec{\mu})|^{2}}{\frac{T}{2}S_{n}(f_{k}|\vec{\lambda})^{2}}]\frac{\partial
S_{n}(f_{k}|\vec{\lambda})}{\partial\lambda^{i}}\,.$ (73)

The first of these expressions is odd in the noise component,
$\tilde{n}(f_{k})=\tilde{s}(f_{k})-\tilde{h}(f_{k}|\vec{\mu})$, while the
second term is even. Since $\mathbb{E}[\tilde{n}(f_{k})]=0$, from this we
deduce that

$\mathbb{E}_{\cal L}\left[\frac{\partial l}{\partial\mu^{i}}\frac{\partial
l}{\partial\lambda^{j}}\right]=0,$ (74)

i.e., at this level of approximation the terms in the Fisher matrix that mix
signal and noise parameters vanish. We conclude that the estimation of the
noise parameters and of the signal parameters is, at leading order,
independent. Lack of knowledge of the noise should therefore not significantly
affect measurements of the parameters of deterministic signals, except
indirectly through the change in the spectral density that enters the
likelihood for the deterministic sources.

##  Appendix C Impact of instrumental noise knowledge uncertainty on SGWB
recovery in absence of galactic foreground

Below we show the same computations we did in Sec. III but in case we do not
consider the presence of the foreground.

####  C.0.1 Power law

Figure 18 shows the results computed for the power law model. We see that in
the presence of instrumental noise uncertainties, the uncertainty in the SGWB
parameters increases by a factor of $\sim 19$–$36$, with the uncertainty in
the slope being more affected than that of the amplitude. Considering the raw
uncertainties, to achieve the same measurement precision, the background
energy density would have to be $\sim 33$ times larger than it would need to
be in the absence of noise knowledge uncertainties. However, a background with
amplitude equal to the reference value should be detectable even allowing for
confusion with instrumental noise mis-modelling.

Figure 18:  Results for the power law SGWB model without foreground

####  C.0.2 Power law with running

Figure 19:  Results for the power-law with running SGWB model without
foreground.

The results for this model are shown in Figure 19. In this case we see that
the uncertainties in the SGWB parameters increase by a factor of $\sim
21$–$72$, with the uncertainty on the log-amplitude being most affected in
this case. Once again, the relative increase in the uncertainty is somewhat
lower at higher background amplitudes. The lower panel of fig. 19 shows that
the background is not detectable at the reference amplitude. An energy density
$\sim 5$ times higher would be required for a detection. In general, the
background again has to have an energy density $\sim 100$ times higher to be
characterised with the same measurement precision when there is instrumental
noise uncertainty as it could be without those uncertainties.

####  C.0.3 Gaussian bump

Figure 20:  Results for Gaussian bump SGWB model without foreground

The results for this model are shown in Figure 20. In this case, the
degradation in the precision of parameter measurement is a factor of $\sim
2.3$–$4$ when allowing for lack of knowledge of the instrumental noise. From
the lower panel of Figure 20 we see that the energy density in a Gaussian bump
SGWB has to be just a small factor of $\sim 2.5$ times bigger to achieve the
same measurement precision when the instrumental noise is not known perfectly.
Moreover a Gaussian bump background at the reference amplitude can be measured
to percent precision at the reference amplitude. The width of the Gaussian can
be measured to a few tens of percent precision at the reference amplitude,
improving approximately linearly with the background energy density.

####  C.0.4 First order phase transition

Figure 21:  Results for first order phase transition SGWB model without
foreground.

The results for First order phase transition are shown in Figure 21. The
results for this SGWB model are quite similar to those for the power law
background. When allowing for instrumental noise knowledge uncertainties, the
precision with which the SGWB log-amplitude can be characterised degrades by a
factor of $\sim 18$. The degradation in the determination of the spectral
index is even larger, $\sim 50$. Once again, to achieve the same measurement
precision, the background energy density would have to be $\sim 50$ times
larger than it would need to be in the absence of noise knowledge
uncertainties. Nonetheless, a FOPT background at the reference amplitude would
still be detectable and provide a measurement of the spectral index at the
level of about 0.5 percent.

We can do the same analysis we did in section III.2 but without including the
foreground. The results for all four SGWB models are qualitatively similar
among themselves and also to the previous case with foreground in Sec. III.2.
Figure 22 consider the power law case, Fig. 23 the case of power law with
running, Fig. 24 the case of gaussian model and Fig. 25 the case of FOPT
model.

Figure 22:  As Figure 18, but now for fixed background amplitude and varying
the variance of the Gaussian prior on the instrumental noise spline model.
This plot is for a power law background, and the amplitude has been fixed such
that the SNR in TDI channel A is 136 (continous lines) and 43 (dashed lines)

The main conclusion from these results is that again if we wanted to ensure
that there was no degradation in LISA science due to lack of noise knowledge,
the necessary requirement on the noise knowledge would be $\ll 10\%$.

Figure 23:  As Figure 22 but now for the power-law-with-running model. The
background amplitude has been fixed to give an overall SNR of $145$.

Figure 24:  As Figure 22 but now for the Gaussian bump model. The background
amplitude has been fixed to give an overall SNR of $135$.

Figure 25:  As Figure fig. 22 but now for the FOPT model. The background
amplitude has been fixed to give an overall SNR of $149$.

###  C.1 Signal reconstruction without foreground of Galactic binaries

We consider a power law signal with an SNR of 48.70 in Fig. 26. The three
panels show the reconstructed ASDs for the SGWB and for the instrumental noise
and the total, which is the sum of the three. No foreground has been
considered in this case. In Fig 27 we show corresponding results for a power
law with a higher SNR, of 862.

Figure 26:  Signal and noise reconstruction for a power law SGWB with SNR
48.70. We show the signal and noise PSDs corresponding to random draws from
the posterior, approximated using the Fisher matrix as described in the text.
In each panel the curves correspond to the three TDI channels: A (blue), E
(red) and $\zeta$ (green). Upper panel: reconstructed SGWB; middle panel:
reconstructed instrumental noise; lower panel: total reconstructed ASD (signal
+ noise).

Figure 27:  As Figure 26, but for a power law background with higher SNR, of
862.

We see that our ability to reconstruct the signal component of the data stream
is poor when the SNR is low. However, we are able to obtain good measurements
of the instrumental noise and the total spectral density. We note that the
total ASD reconstruction in Figure 26 is somewhat poorer than the noise-only
component, which does not fit with the expectation that we are actually
measuring the total. This happens due to the breakdown in the Fisher matrix
approximation for the SGWB parameters in this case, because the SGWB parameter
uncertainties are large and no longer in the linear signal regime. At higher
SNR, we start to be able to reconstruct the SGWB more precisely, shown by a
reduction in the scatter in Figure 27. As the SNR is increased we would expect
the scatter to reduce further. The reconstruction of the noise spectral
density is comparable to what is seen in the lower SNR case, but we would
eventually expect it to degrade as the SGWB becomes more dominant in the data.
The reconstruction of the total spectral density is similar to the low SNR
case, as expected. However, this higher SNR case does not show the noise at
low frequency that arises from the breakdown of the Fisher matrix
approximation, presumably because the measurement uncertainties are within the
linear regime in this case. although, what is interesting to notice in the
total reconstruction is that the channels A and E are affected by the SGWB
where between $1mHz$ and $4mHz$ the noise level deviated from the expected one
shown in Fig 2.

## References

  * [1] Lisa red-book, in prep.
  * [2] Matthew R. Adams and Neil J. Cornish.  Discriminating between a Stochastic Gravitational Wave Background and Instrument Noise.  Phys. Rev. D, 82:022002, 2010.
  * [3] Bruce Allen and Joseph D. Romano.  Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities.  Phys. Rev. D, 59(10):102001, May 1999.
  * [4] Pau Amaro-Seoane et al.  Laser Interferometer Space Antenna.  2 2017.
  * [5] M. Armano, H. Audley, G. Auger, J. T. Baird, M. Bassan, P. Binetruy, M. Born, D. Bortoluzzi, N. Brandt, M. Caleno, L. Carbone, A. Cavalleri, A. Cesarini, G. Ciani, G. Congedo, A. M. Cruise, K. Danzmann, M. de Deus Silva, R. De Rosa, M. Diaz-Aguiló, L. Di Fiore, I. Diepholz, G. Dixon, R. Dolesi, N. Dunbar, L. Ferraioli, V. Ferroni, W. Fichter, E. D. Fitzsimons, R. Flatscher, M. Freschi, A. F. García Marín, C. García Marirrodriga, R. Gerndt, L. Gesa, F. Gibert, D. Giardini, R. Giusteri, F. Guzmán, A. Grado, C. Grimani, A. Grynagier, J. Grzymisch, I. Harrison, G. Heinzel, M. Hewitson, D. Hollington, D. Hoyland, M. Hueller, H. Inchauspé, O. Jennrich, P. Jetzer, U. Johann, B. Johlander, N. Karnesis, B. Kaune, N. Korsakova, C. J. Killow, J. A. Lobo, I. Lloro, L. Liu, J. P. López-Zaragoza, R. Maarschalkerweerd, D. Mance, V. Martín, L. Martin-Polo, J. Martino, F. Martin-Porqueras, S. Madden, I. Mateos, P. W. McNamara, J. Mendes, L. Mendes, A. Monsky, D. Nicolodi, M. Nofrarias, S. Paczkowski, M. Perreur-Lloyd, A. Petiteau, P. Pivato, E. Plagnol, P. Prat, U. Ragnit, B. Raïs, J. Ramos-Castro, J. Reiche, D. I. Robertson, H. Rozemeijer, F. Rivas, G. Russano, J. Sanjuán, P. Sarra, A. Schleicher, D. Shaul, J. Slutsky, C. F. Sopuerta, R. Stanga, F. Steier, T. Sumner, D. Texier, J. I. Thorpe, C. Trenkel, M. Tröbs, H. B. Tu, D. Vetrugno, S. Vitale, V. Wand, G. Wanner, H. Ward, C. Warren, P. J. Wass, D. Wealthy, W. J. Weber, L. Wissel, A. Wittchen, A. Zambotti, C. Zanoni, T. Ziegler, and P. Zweifel.  Sub-femto-$g$ free fall for space-based gravitational wave observatories: Lisa pathfinder results.  Phys. Rev. Lett., 116:231101, Jun 2016.
  * [6] J. W. Armstrong, F. B. Estabrook, and Massimo Tinto.  Time-delay interferometry for space-based gravitational wave searches.  The Astrophysical Journal, 527(2):814, dec 1999.
  * [7] Stanislav Babak, Chiara Caprini, Daniel G. Figueroa, Nikolaos Karnesis, Paolo Marcoccia, Germano Nardini, Mauro Pieroni, Angelo Ricciardone, Alberto Sesana, and Jesús Torrado.  Stochastic gravitational wave background from stellar origin binary black holes in lisa, 2023.
  * [8] Stanislav Babak, Jonathan Gair, Alberto Sesana, Enrico Barausse, Carlos F. Sopuerta, Christopher P. L. Berry, Emanuele Berti, Pau Amaro-Seoane, Antoine Petiteau, and Antoine Klein.  Science with the space-based interferometer LISA. V: Extreme mass-ratio inspirals.  Phys. Rev. D, 95(10):103012, 2017.
  * [9] Stanislav Babak, Antoine Petiteau, and Martin Hewitson.  LISA Sensitivity and SNR Calculations.  8 2021.
  * [10] Quentin Baghi, Nikolaos Karnesis, Jean-Baptiste Bayle, Marc Besançon, and Henri Inchauspé.  Uncovering gravitational-wave backgrounds from noises of unknown shape with LISA.  2 2023.
  * [11] N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso, and A. Riotto.  Primordial Black Hole Dark Matter: LISA Serendipity.  Phys. Rev. Lett., 122(21):211301, 2019.
  * [12] N. Bartolo, V. De Luca, G. Franciolini, M. Peloso, D. Racco, and A. Riotto.  Testing primordial black holes as dark matter with LISA.  Physical Review D, 99(10), may 2019.
  * [13] Nicola Bartolo et al.  Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves.  JCAP, 12:026, 2016.
  * [14] Jean-Baptiste Bayle, Aurélien Hees, Marc Lilley, and Christophe Le Poncin-Lafitte.  Lisa orbits, April 2022.
  * [15] Chiara Caprini and Daniel G. Figueroa.  Cosmological Backgrounds of Gravitational Waves.  Class. Quant. Grav., 35(16):163001, 2018.
  * [16] Chiara Caprini, Daniel G. Figueroa, Raphael Flauger, Germano Nardini, Marco Peloso, Mauro Pieroni, Angelo Ricciardone, and Gianmassimo Tasinato.  Reconstructing the spectral shape of a stochastic gravitational wave background with LISA.  Journal of Cosmology and Astroparticle Physics, 2019(11):017–017, nov 2019.
  * [17] Chiara Caprini, Mark Hindmarsh, Stephan Huber, Thomas Konstandin, Jonathan Kozaczuk, Germano Nardini, Jose Miguel No, Antoine Petiteau, Pedro Schwaller, Gé raldine Servant, and David J. Weir.  Science with the space-based interferometer eLISA. II: gravitational waves from cosmological phase transitions.  Journal of Cosmology and Astroparticle Physics, 2016(04):001–001, apr 2016.
  * [18] M. Chwalla, K. Danzmann, M. Dovale Álvarez, J.J. Esteban Delgado, G. Fernández Barranco, E. Fitzsimons, O. Gerberding, G. Heinzel, C.J. Killow, M. Lieser, M. Perreur-Lloyd, D.I. Robertson, J.M. Rohr, S. Schuster, T.S. Schwarze, M. Tröbs, G. Wanner, and H. Ward.  Optical suppression of tilt-to-length coupling in the lisa long-arm interferometer.  Phys. Rev. Appl., 14:014030, Jul 2020.
  * [19] LISA Pathfinder collaboration.  In-depth analysis of LISA Pathfinder performance results: time evolution, noise projection, physical models and implications for LISA.
  * [20] Raphael Flauger, Nikolaos Karnesis, Germano Nardini, Mauro Pieroni, Angelo Ricciardone, and Jesús Torrado.  Improved reconstruction of a stochastic gravitational wave background with LISA.  JCAP, 01:059, 2021.
  * [21] G. Franciolini.  Primordial black holes: from theory to gravitational wave observations, 2021.
  * [22] Davide Gerosa, Sizheng Ma, Kaze W. K. Wong, Emanuele Berti, Richard O’Shaughnessy, Yanbei Chen, and Krzysztof Belczynski.  Multiband gravitational-wave event rates and stellar physics.  Phys. Rev. D, 99(10):103004, 2019.
  * [23] Marie-Sophie Hartig, Sönke Schuster, and Gudrun Wanner.  Geometric tilt-to-length coupling in precision interferometry: mechanisms and analytical descriptions.  Journal of Optics, 24(6):065601, may 2022.
  * [24] Olaf Hartwig.  Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna.  PhD thesis, Leibniz U., Hannover, 2021.
  * [25] Olaf Hartwig and Jean-Baptiste Bayle.  Clock-jitter reduction in lisa time-delay interferometry combinations.  Phys. Rev. D, 103:123027, Jun 2021.
  * [26] Olaf Hartwig, Marc Lilley, Martina Muratore, and Mauro Pieroni.  Stochastic gravitational wave background reconstruction for a non-equilateral and unequal-noise LISA constellation.  3 2023.
  * [27] Olaf Hartwig and Martina Muratore.  Characterization of time delay interferometry combinations for the lisa instrument noise.  Phys. Rev. D, 105:062006, Mar 2022.
  * [28] Mark Hindmarsh, Stephan J. Huber, Kari Rummukainen, and David J. Weir.  Shape of the acoustic gravitational wave power spectrum from a first order phase transition.  Phys. Rev. D, 96:103520, Nov 2017.
  * [29] Antoine Klein et al.  Science with the space-based interferometer eLISA: Supermassive black hole binaries.  Phys. Rev. D, 93(2):024003, 2016.
  * [30] Leonard Lehoucq, Irina Dvorkin, Rahul Srinivasan, Clement Pellouin, and Astrid Lamberts.  Astrophysical uncertainties in the gravitational-wave background from stellar-mass compact binary mergers, 2023.
  * [31] Christopher J. Moore, Davide Gerosa, and Antoine Klein.  Are stellar-mass black-hole binaries too quiet for LISA?  Mon. Not. Roy. Astron. Soc., 488(1):L94–L98, 2019.
  * [32] Martina Muratore.  Time delay interferometry for LISA science and instrument characterization.  PhD thesis, Trento U., 2021.
  * [33] Martina Muratore, Olaf Hartwig, Daniele Vetrugno, Stefano Vitale, and William Joseph Weber.  On the effectiveness of null TDI channels as instrument noise monitors in LISA.  7 2022.
  * [34] Martina Muratore, Daniele Vetrugno, and Stefano Vitale.  Revisitation of time delay interferometry combinations that suppress laser noise in lisa, 2020.
  * [35] Martina Muratore, Daniele Vetrugno, Stefano Vitale, and Olaf Hartwig.  Time delay interferometry combinations as instrument noise monitors for LISA.  Phys. Rev. D, 105(2):023009, 2022.
  * [36] Dam Quang Nam, Yves Lemiere, Antoine Petiteau, Jean-Baptiste Bayle, Olaf Hartwig, Joseph Martino, and Martin Staab.  Tdi noises transfer functions for lisa, 2022.
  * [37] S. Paczkowski, R. Giusteri, M. Hewitson, N. Karnesis, E. D. Fitzsimons, G. Wanner, and G. Heinzel.  Postprocessing subtraction of tilt-to-length noise in lisa.  Phys. Rev. D, 106:042005, Aug 2022.
  * [38] B. Picinbono.  Second-order complex random vectors and normal distributions.  IEEE Transactions on Signal Processing, 44(10):2637–2640, 1996\.
  * [39] Thomas A. Prince, Massimo Tinto, Shane L. Larson, and J. W. Armstrong.  Lisa optimal sensitivity.  Phys. Rev. D, 66:122002, Dec 2002.
  * [40] Wen-Hong Ruan, Zong-Kuan Guo, Rong-Gen Cai, and Yuan-Zhong Zhang.  Taiji Program: Gravitational-Wave Sources.  arXiv e-prints, page arXiv:1807.09495, July 2018.
  * [41] Lorenzo Sala.  Residual test mass acceleration in LISA Pathfinder: in-depth statistical analysis and physical sources.  PhD thesis, University of Trento, 2023.
  * [42] Etienne Savalle, Jonathan Gair, Lorenzo Speri, and Stanislav Babak.  Assessing the impact of instrumental calibration uncertainty on LISA science.  Physical Review D, 106(2), jul 2022.
  * [43] Alberto Sesana.  Prospects for multiband gravitational-wave astronomy after gw150914.  Phys. Rev. Lett., 116:231102, Jun 2016.
  * [44] Eric Thrane and Joseph D. Romano.  Sensitivity curves for searches for gravitational-wave backgrounds.  Phys. Rev. D, 88(12):124032, 2013.
  * [45] M. Tinto and S.V Dhurandhar.  Time-delay interferometry.  Living Rev Relativ 24, 2021.
  * [46] Massimo Tinto and Olaf Hartwig.  Time-Delay Interferometry and Clock-Noise Calibration.  Phys. Rev. D, 98(4):042003, 2018.