# Model-independent time-delay interferometry based on principal component
analysis

Quentin Baghi<EMAIL_ADDRESS>IRFU, CEA, Université Paris-Saclay,
F-91191 Gif-sur-Yvette, France  John Baker  Jacob Slutsky  James Ira Thorpe
Goddard Space Flight Center
Mail Code 663
8800 Greenbelt Rd, Greenbelt, Maryland 20771, USA

###### Abstract

With a laser interferometric gravitational-wave detector in separate free
flying spacecraft, the only way to achieve detection is to mitigate the
dominant noise arising from the frequency fluctuations of the lasers via
postprocessing. The noise can be effectively filtered out on the ground
through a specific technique called time-delay interferometry (TDI), which
relies on the measurements of time-delays between spacecraft and careful
modeling of how laser noise enters the interferometric data. Recently, this
technique has been recast into a matrix-based formalism by several authors,
offering a different perspective on TDI, particularly by relating it to
principal component analysis (PCA). In this work, we demonstrate that we can
cancel laser frequency noise by directly applying PCA to a set of shifted data
samples, without any prior knowledge of the relationship between single-link
measurements and noise, nor time-delays. We show that this fully data-driven
algorithm achieves a gravitational-wave sensitivity similar to classic TDI.

##  I Introduction

Gravitational waves can be detected through a constellation of satellites
forming a network of laser interferometers that probe phase differences in
beams exchanged between spacecraft. This concept is the core of future space-
based detectors like the Laser Interferometer Space Antenna (LISA) [1] and
TianQuin [2]. It was soon realized [3] that the noise caused by the stochastic
fluctuations of laser frequencies cannot be canceled at the photometric
detectors onboard the satellites. The reason lies in the difference between
the delays they experience when being transmitted along two interferometric
arms. Rather, interferometry must be performed once the telemetry is retrieved
on the ground. By taking advantage of the fact that several delayed versions
of the same noise component appear in different measurements, it is possible
to cancel laser noise by combining and interpolating the received data streams
at the right timestamps. We call this post-processing technique time-delay
interferometry (TDI) [4]. TDI has been undergoing a constant development to
account for increasing degrees of complexity, including optical bench motion
noise [5], varying armlengths [6, 7], noise correlations [8], on-board
filtering [9], clock jitters [10] and the preferred representative units of
the processed data (phase or fractional-frequency) [11]. Other efforts focused
on finding new noise-canceling combinations [12].

Romano and Woan [13] showed that we could think of TDI combinations as low-
variance principal components of the laser noise covariance matrix. This idea
relied on a matrix representation of the interferometric measurements and was
also explored in Ref. [14]. Computing TDI variables is equivalent to
projecting the data onto the null space of the covariance. Recently, this
concept regained interest in the community. Vallisneri et al. [15] presented a
similar algebraic computation (called TDI-$\infty$) where they perform the
singular value decomposition (SVD) of the design matrix $\bm{M}$, i.e., the
matrix translating individual laser noise sources $\bm{c}$ into measured
noises $\bm{y}$ so that $\bm{y}=\bm{M}\bm{c}$. They approximate the model
likelihood by restricting it to the dominating terms, using the singular
vectors that correspond to the null space of $\bm{M}$. Concomitantly, we
derived a frequency-domain version of the covariance eigenvector decomposition
[16] called PCI (for principal component interferometry), where we focused on
the orthogonalization of the data with respect to secondary (i.e., non-laser)
noises. We showed that this method generalizes the A, E, and T TDI variables
[8] and can bring an improvement in signal-to-noise ratio (SNR) when the
levels of the noises affecting the single-link measurements are unequal.
Although the implementation details are still under debate [17, 18, 19], the
matrix representation of TDI is a powerful tool for data analysis, as it
allows one to express the model likelihood directly in terms of single-link
measurements and readily account for data gaps.

TDI and the extensions cited above rely on accurate modeling of the dependence
of the measurements on the various laser noise sources. TDI stems from a set
of equations modeling every interferometer, from which one derives the
sequence of combination and delaying operations needed to precisely cancel
laser noise terms. Equivalently, TDI-$\infty$ assumes perfect knowledge of the
design matrix, while PCI assumes a model of the noise covariance. In addition,
these methods need precise estimates of the time-delays, which can be obtained
independently from the telemetry [20] or directly estimated from the data
using noise power minimization [21] or, as it was recently shown, by Bayesian
inference [22, 16]. However, LISA is a complex apparatus involving many
control laws, couplings, and frequency planning inducing correlations among
lasers. Inaccurate modeling can lead to imperfect cancellation or even no
laser noise suppression at all [9, 10]. In this work, we present a fully data-
driven approach to cancel laser noise in space interferometry, with no prior
assumption on the noise structure. The method (that we name aPCI for automatic
PCI) relies on the principal component analysis (PCA) of a matrix directly
constructed from the measurements.

PCA is usually used to reduce the dimensionality of a problem by decomposing
the data into an orthogonal basis whose vectors are called principal
components. Then, one keeps only the vectors that contain the most information
about the data by restricting the analysis to the components having the
largest variance. However, when the noise overwhelmingly dominates the signal
of interest, as this is the case for gravitational waves (GW) by about a
factor $10^{8}$, we can use PCA with the opposite objective: to retain the
components having the lowest variance.

In most problems, the decomposition is performed on a collection of $n$
different realizations of the same $p$-dimensional measurement, forming a
matrix of size $n\times p$. However, in a gravitational-wave observatory like
LISA, the data is instead made of different channels affected by different
noises and measuring different projections of the same signal. To build a data
matrix that we can process with PCA, we can draw inspiration from the
technique used to form the TDI combinations [23], which relies on
interpolators like Lagrange filters to produce delayed versions of the
measured data. The measurement estimate at a given point in time is
interpolated by linearly combining nearby data samples. Thus, we build a
$n\times p$ data matrix whose rows are versions of the same data shifted by an
integer number of samples. Then, we compute the PCA of this matrix to obtain
components organized in decreasing order of variance. The lowest variance
components correspond to combinations where laser noise is significantly
reduced. Projecting the data onto these components yields data streams that
are readily usable for gravitational-wave searches.

We first detail the principle of the aPCI method in Section II, and present
the theoretical framework. In Section III, we provide a means to evaluate the
method performance by deriving the sensitivity of aPCI combinations to
monochromatic GWs in the case of LISA. In Section IV, we use a numerical
simulation of LISA data to demonstrate laser noise suppression, and also
compare aPCI with classic TDI by computing their sensitivities averaged over
sky location and orientation. Finally, in Section V we discuss important
aspects of this approach and pave the way for future studies.

##  II Principal component analysis theoretical framework

###  II.1 Construction of the data matrix

Space-based gravitational-wave detector data include the outputs of several
interferometers, which can be expressed in relative frequency deviations. We
denote the corresponding measurements as $s_{ij}(t)$ to conform to the
official convention (see, for example, Ref. [11]), where $i$ is the index of
the satellite hosting the optical bench, and $j$ is the index of the sending
spacecraft. In practice, we measure a discrete version of $s_{ij}(t)$ with $N$
samples spaced every $\tau_{s}$ seconds. Therefore, we can form a vector
$\bm{s}_{ij}$ with $N$ entries given by $s_{ij}(t_{n})$, where the time stamps
are $t_{n}=n\tau_{s},\,n\in\left[0,\,N-1\right]$. Hence, we have
$\bm{s}_{ij}\equiv\left(s_{ij}(t_{0}),\,\ldots,\,s_{ij}(t_{N-1})\right)^{\dagger}$.
We will stick to the LISA example, where the detector is made of a
constellation of 3 satellites, each of them carrying 2 optical benches on
board. If we consider the science measurements only, then we have 6
measurements that we can stack into a $N\times 6$ matrix $\bm{Y}$ as

$\displaystyle\bm{Y}\equiv\left(\bm{s}_{12},\,\bm{s}_{23},\,\bm{s}_{31},\,\bm{s}_{13},\,\bm{s}_{21},\,\bm{s}_{32}\right).$
(1)

Classical TDI algorithms usually operate in two steps. First, they compute
delayed versions of discrete signals, interpolated at specific times depending
on the light travel time delays along the constellation arms. In a second
step, they combine these interpolated time series in such a way that laser
noise terms vanish.

We would like to form laser-noise canceling combinations of the data without
knowing the delays nor the way they enter the measurements. We know that the
first step of TDI is usually performed through fractional-delay filter
techniques which use linear combinations of samples at times before and after
the target date [23]. In this process, one must choose the number of samples
$n_{h}$ to use before and after the desired interpolation time. In principle,
the larger $n_{h}$, the better the interpolation accuracy, but the longer the
computation time. To provide a fully data-driven algorithm with the
ingredients necessary to automatically find the correlations among the data,
we can therefore start from shifted versions of the measurement vector
$\bm{s}$. Let us assume that we have measured data over a time span from $t=0$
to $t=T_{\mathrm{obs}}$, sampled every $\tau_{s}$. We can build a matrix
$\bm{X}$ containing $\bm{Y}$ itself, along with $n_{h}$ backward-shifted
versions of it and $n_{h}$ forward-shifted versions, resulting in a matrix of
size $N\times 6(2n_{h}+1)$ such that:

$\displaystyle\bm{X}\equiv\left(\bm{D}_{-n_{h}}\bm{Y},\,\bm{D}_{-n_{h}+1}\bm{Y},\,\ldots,\,\bm{D}_{n_{h}}\bm{Y}\right),$
(2)

where we defined the integer delay matrix $\bm{D}_{m}$ of size $N\times N$
which acts on any vector $\bm{s}$ as

${\bm{D}_{m}\bm{s}}_{n}=\begin{cases}\bm{s}_{n-m}&\text{if
}n\in\left[n_{\mathrm{min}},\,n_{\mathrm{max}}\right]\\\ 0&\text{ otherwise,
}\end{cases}$

with $n_{\mathrm{min}}\equiv\mathrm{max}(0,m)$ and
$n_{\mathrm{max}}\equiv\mathrm{min}(N-1,N-1+m)$. This way, the matrix $\bm{X}$
is formed of sub-vectors including 6 interferometric measurements, where a
sub-vector is shifted by one time sample with respect to the previous one. In
this form, the data is ready to be processed through PCA. In the following, we
will set $p\equiv 2n_{h}+1$ as the number of subchannels.

###  II.2 PCA decomposition of the data matrix

PCA aims at finding a transformation of the data matrix into a new set of
variables which are as independent as possible with respect of the information
they carry, measured in terms of variance. This is done by computing a set of
vectors called principal components (PCs) representing orthogonal directions
that are ordered from the one that best fits the data, to the one that is less
representative. In the problem we study, the data matrix variance
overwhelmingly represents laser noise. The first PCs will therefore better fit
the laser noise. On the contrary, the last PCs should be more independent of
laser noise, representing the remaining information present in the data.

We can compute the PCs through singular value decomposition (SVD):

$\bm{X}=\bm{U}\bm{S}\bm{V}^{\intercal},$ (3)

where $\bm{U}$ is a $N\times N$ unitary matrix of left-singular vectors,
$\bm{S}$ is a $N\times 6p$ rectangular diagonal matrix whose non-zero elements
are the singular values, $\bm{V}$ is a $6p\times 6p$ unitary matrix of right-
singular vectors, and $\intercal$ is the matrix transpose operator. The PCs
are given by the vectors forming the columns of $\bm{V}$.

Note that, in regular PCA, the columns of the matrix we decompose correspond
to different features (like different sensor outputs) and the rows correspond
to different realizations of the measurement. In the way we apply it to laser
noise processing, the matrix columns include not only different interferometer
outputs, but also different shifts, while the rows correspond to particular
times at which the data is measured. Therefore, to draw a comparison between
TDI and what we do here, the elements of particular PC vector could be
interpreted as the kernels of several fractional delay filters.

The principal component decomposition of $\bm{X}$ is therefore the matrix
$\bm{T}$ resulting in the transformation

$\bm{T}\equiv\bm{X}\bm{V}.$ (4)

The rows of $\bm{T}$ are usually called the scores, but we will refer to them
as the PCA (or aPCI) combinations in the following.

In this problem, we are interested in the $q$ lowest-variance components of
the decomposition, below some threshold where the laser noise is significantly
suppressed. Let us rewrite $\bm{V}$ in two parts: a part
$\bm{V}_{q}=\bm{V}\bm{\Pi}_{q}$, where
$\bm{\Pi}_{q}\equiv\left(\bm{0},\,\bm{I}_{q}\right)^{\intercal}$ is the
$6p\times q$ matrix that projects the last $q$ low-variance PCs, and the
complementary part $\bm{V}_{\setminus{q}}$ containing the first $2p-q$ PCs
related to large variance:

$\begin{pmatrix}\bm{T}_{\setminus{q}}\\\
\bm{T}_{q}\end{pmatrix}\equiv\bm{X}\begin{pmatrix}\bm{V}_{\setminus{q}}&\bm{V}_{q}\end{pmatrix},$
(5)

where we split $\bm{V}=\left(\bm{V}_{\setminus{q}},\,\bm{V}_{q}\right)$. To a
good approximation, the low-variance transformation $\bm{T}_{q}$ should be all
we need to search for the signals of interest, i.e. gravitational waves.

A key effect of the PCA is that the transformation from $\bm{X}$ to $\bm{T}$
diagonalizes the data matrix sample covariance
$\bm{\hat{\Sigma}}_{\bm{X}}=\frac{1}{N}\bm{X}^{\intercal}\bm{X}$ and
transforms it to the diagonal matrix
$\bm{\hat{\Sigma}}_{\bm{T}}=\frac{1}{N}\bm{T}^{\intercal}\bm{T}=\frac{1}{N}\bm{S}^{\intercal}\bm{S}$.
Here we are primarily interested in the _inverse_ covariance, which the low-
variance decomposition allows us to write as:

$\displaystyle\bm{\hat{\Sigma}}^{-1}_{\bm{T}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\bm{V}_{\setminus{q}}^{\intercal}\bm{\hat{\Sigma}}^{-1}_{\bm{X}}\bm{V}_{\setminus{q}}&\bm{V}_{\setminus{q}}^{\intercal}\bm{\hat{\Sigma}}^{-1}_{\bm{X}}\bm{V}_{q}\\\
\bm{V}_{q}^{\intercal}\bm{\hat{\Sigma}}^{-1}_{\bm{X}}\bm{V}_{\setminus{q}}&\bm{V}_{q}^{\intercal}\bm{\hat{\Sigma}}^{-1}_{\bm{X}}\bm{V}_{q}\end{pmatrix}$
(6) $\displaystyle\approx$ $\displaystyle\begin{pmatrix}\bm{0}&\bm{0}\\\
\bm{0}&\bm{V}_{q}^{\intercal}\bm{\hat{\Sigma}}^{-1}_{\bm{X}}\bm{V}_{q}\end{pmatrix}$
$\displaystyle=$
$\displaystyle\bm{\Pi}^{\intercal}_{q}\bm{\Pi}_{q}\bm{\hat{\Sigma}}^{-1}_{\bm{T}}\bm{\Pi}^{\intercal}_{q}\bm{\Pi}_{q}$

where in the last line we used the diagonality of $\bm{\hat{\Sigma}}_{\bm{T}}$
and neglected the first $6p-q$ inverse variances compared to the $q$ last
ones. Here our reverse-ordered preference for the nominally least-significant
principal components has the effect of selecting the most significant parts of
the inverse covariance matrix.

###  II.3 Relation of the PCA process with the model covariance

As is common in gravitational-wave data analysis, we assume our noise is
Gaussian and stationary (over the relatively short timescales relevant for
laser noise.) Stationarity allows us to assume that frequency bins are
approximately uncorrelated in the Fourier domain. For any vector $\bm{z}$, we
will label as $\bm{\tilde{z}}(f)$ its discrete Fourier transform (DFT) element
at frequency $f$. Similarly, we will denote by $\bm{\tilde{\Sigma}}_{z}(f)$
its Fourier-domain covariance, which fully characterizes the (zero-mean)
noise.

The construction of the data matrix in Eq. (2) can be rewritten in the
frequency domain as

$\bm{\tilde{X}}(f)=\bm{\tilde{y}}^{\intercal}(f)\begin{pmatrix}\bm{\tilde{D}}_{-n_{h}}(f),\,\ldots,\,\bm{\tilde{D}}_{+n_{h}}(f)\end{pmatrix},$
(7)

where $\bm{\tilde{y}}(f)$ is now a $6\times 1$ column vector, and
$\bm{\tilde{X}}(f)$ is a $1\times 6p$ row vector. Here we approximate the
$m^{\mathrm{th}}$ integer delay operator involved in Eq. (2) by a
multiplicative complex exponential phase in the Fourier domain, encoded by a
$6\times 6$ diagonal matrix $\bm{\tilde{D}}_{m}(f)$ with entries

$\tilde{D}_{m,pq}(f)=\delta_{pq}e^{-2i\pi fm\tau_{s}},$ (8)

where $i=\sqrt{-1}$ is the complex number and $\delta_{pq}$ is the Kronecker
symbol ($\delta_{pq}=1\text{ if }p=q$ and $\delta_{pq}=0$ otherwise).

Similarly, we can write a frequency-domain version of the projection in Eq.
(4) by defining $\bm{\tilde{t}}(f)$ as the $6p\times 1$ vector formed by
stacking the $6p$ PCA combinations at frequency $f$. Then, we have

$\displaystyle\bm{\tilde{t}}(f)$
$\displaystyle=\bm{V}^{\intercal}\bm{\tilde{X}}^{\intercal}(f)=\bm{W}(f)\bm{\tilde{y}}(f),$
(9)

where we set
$\bm{\tilde{W}}(f)\equiv\bm{V}^{\intercal}\big{(}\bm{\tilde{D}}_{-n_{h}}(f),\,\ldots,\,\bm{\tilde{D}}_{+n_{h}}(f)\big{)}^{\intercal}$
the $6p\times 6$ matrix transformation which encodes, in the Fourier domain,
the construction of the data matrix and its PCA transformation.

The transformation matrix $\bm{\tilde{W}}$ has useful properties:

$\displaystyle\bm{\tilde{W}}^{{\dagger}}\bm{\tilde{W}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\bm{\tilde{D}}_{-n_{h}}^{{\dagger}}&\ldots&\bm{\tilde{D}}_{+n_{h}}^{{\dagger}}\end{pmatrix}\bm{V}\bm{V}^{\intercal}\begin{pmatrix}\bm{\tilde{D}}_{-n_{h}}\\\
\vdots\\\ \bm{\tilde{D}}_{+n_{h}}\end{pmatrix}$ (10) $\displaystyle=$
$\displaystyle\sum_{m=-n_{h}}^{+n_{h}}\bm{\tilde{D}}_{m}^{\dagger}\bm{\tilde{D}}_{m}=p\bm{I}_{6},$

where we used the fact that $\bm{V}$ is unitary. Note that here we took the
Hermitian transpose ${\dagger}$, because $\bm{\tilde{W}}$ is a complex-valued
matrix.

We denote by
$\bm{\tilde{\Sigma}}_{\bm{x}}\equiv\operatorname{E}\left[\bm{\tilde{x}}{\bm{\tilde{x}}}^{\dagger}\right]$
the covariance of any zero-mean complex random vector $\bm{\tilde{x}}$. Then
from Eq. (9) we can write the $6p\times 6p$ covariance of $\bm{\tilde{t}}$ as

$\bm{\tilde{\Sigma}}_{\bm{t}}=\bm{\tilde{W}}\bm{\tilde{\Sigma}}_{\bm{y}}\bm{\tilde{W}}^{\dagger},$
(11)

where we dropped the dependence on $f$ to lighten the expressions.

We seek a similar expression relating the inverse covariances. Because the
transformation $\bm{\tilde{t}}=(1/p)\bm{\tilde{W}}^{\dagger}\bm{\tilde{y}}$
involves a projection from $6p$ dimensions down to $6$, we do not have a
useful inverse identity involving $\bm{\tilde{W}}\bm{\tilde{W}}^{\dagger}$.
Consequently it is not straightforward to find an expression for the inverse
covariance by solving $\bm{A}\bm{\tilde{\Sigma}}_{\bm{t}}=\bm{I}_{6p}$.
However, if we set
$\bm{A}=(1/p^{2})\bm{\tilde{W}}\bm{\Sigma}_{\bm{y}}^{-1}\bm{\tilde{W}}^{\dagger}$,
we can verify that
$\bm{\tilde{\Sigma}}_{\bm{t}}\bm{A}\bm{\tilde{\Sigma}}_{\bm{t}}=\bm{\tilde{\Sigma}}_{\bm{t}}$
and $\bm{A}\bm{\tilde{\Sigma}}_{\bm{t}}\bm{A}=\bm{A}$ and that both
$\bm{A}\bm{\tilde{\Sigma}}_{\bm{t}}$ and $\bm{\tilde{\Sigma}}_{\bm{t}}\bm{A}$
are Hermitian, thus establishing that A is the pseudoinverse of
$\bm{\tilde{\Sigma}}_{\bm{t}}$, which we simply denote as

$\bm{\Sigma}_{\bm{t}}^{+}\bm{=}(1/p^{2})\bm{\tilde{W}}\bm{\Sigma}_{\bm{y}}^{-1}\bm{\tilde{W}}^{\dagger}.$
(12)

Using the explicit properties of $\bm{\tilde{W}}$, we can then invert this
expression to find

$\bm{\Sigma}_{\bm{y}}^{-1}=\bm{\tilde{W}}^{\dagger}\bm{\Sigma}_{\bm{t}}^{+}\bm{\tilde{W}}.$
(13)

If the sample variance from the PCA data matrix provides a sufficient estimate
of the (generally distinct) data variance as regards the laser noise, then as
in Section II.2 we can estimate
$\bm{\tilde{\Sigma}}_{\bm{t}}^{+}\approx\bm{\Pi}^{\intercal}_{q}\bm{\Pi}_{q}\bm{\tilde{\Sigma}}_{\bm{t}}^{+}\bm{\Pi}^{\intercal}_{q}\bm{\Pi}_{q}.$
Thus

$\displaystyle\bm{\tilde{\Sigma}}_{\bm{y}}^{-1}$ $\displaystyle\approx$
$\displaystyle\bm{\tilde{W}}^{{\dagger}}_{q}\bm{\Pi}_{q}\bm{\tilde{\Sigma}}_{\bm{t}}^{+}\bm{\Pi}^{\intercal}_{q}\bm{\tilde{W}}_{q}.$
(14)

where we set $\bm{\tilde{W}}_{q}\equiv\bm{\Pi}_{q}\bm{\tilde{W}}$.

This approximation amounts to ignoring the negligible terms depending on laser
noise when computing the inverse covariance, through a suitable transformation
that allows us to restrict the analysis to the terms that are relevant for the
GW parameter inference.

###  II.4 Bayesian inference

Gravitational-wave signals are usually analyzed using Bayesian inference. We
consider the frequency-domain data $\bm{\tilde{y}}$ comprising six channels
evaluated at frequency bin $f$. We drop the explicit dependence on $f$ for
more clarity. For a model parameterized by $\bm{\theta}$, the posterior
distribution given $\bm{\tilde{y}}$ is

$\displaystyle
p\left(\bm{\theta}|\bm{\tilde{y}}\right)=\frac{p\left(\bm{\tilde{y}}|\bm{\theta}\right)p\left(\bm{\theta}\right)}{p\left(\bm{\tilde{y}}\right)},$
(15)

where $p\left(\bm{\tilde{y}}|\bm{\theta}\right)$ is the likelihood,
$p\left(\bm{\theta}\right)$ is the prior distribution of the parameters, and
the evidence $p\left(\bm{\tilde{y}}\right)$ provides normalization.

For zero-mean Gaussian noise, the likelihood has the form:

$\displaystyle
p\left(\bm{\tilde{y}}|\bm{\theta}\right)=\frac{\exp\left\\{-\frac{1}{2}\left(\bm{\tilde{y}}-\bm{\tilde{h}}\right)^{\dagger}\bm{\tilde{\Sigma}}_{\bm{y}}^{-1}\left(\bm{\tilde{y}}-\bm{\tilde{h}}\right)\right\\}}{\sqrt{(2\pi)^{6}\left|\bm{\tilde{\Sigma}}_{\bm{y}}\right|}},$
(16)

where $\dagger$ denotes the Hermitian conjugate, $\bm{\tilde{h}}$ is the
contribution to the data $\bm{\tilde{y}}$ expected from the gravitational wave
signal implied by parameters $\bm{\theta}$, and $\bm{\tilde{\Sigma}}_{\bm{y}}$
is the $6N\times 6N$ covariance. At this level the covariance includes all the
noise processes including dominant contributions from the laser noise.

Bayesian analysis is clearest, when the data represented in the analysis comes
directly from the basic measurements. In practice though some pre-treatment or
approximating assumptions can be important to make the analysis practical. For
LISA, the challenge is that the basic problem represented in Eq. (16) is
dominated by the uninteresting laser noise. TDI processing is one approach to
recast the problem without this dominant part.

Here we propose to apply our PCA treatment instead. To see how this works, we
first plug the approximation for the inverse covariance derived in Eq. (14)
into Eq. (16) and take the logarithm to obtain:

$\displaystyle\ln p\left(\bm{\tilde{y}}|\bm{\theta}\right)$
$\displaystyle\approx$
$\displaystyle-\frac{1}{2}\left(\bm{\tilde{y}}-\bm{\tilde{h}}\right)^{{\dagger}}\bm{\tilde{W}}^{\dagger}_{q}\bm{\Pi}_{q}\bm{\tilde{\Sigma}}_{\bm{t}}^{+}\bm{\Pi}^{{\dagger}}_{q}\bm{\tilde{W}}_{q}\left(\bm{\tilde{y}}-\bm{\tilde{h}}\right)$
(17)
$\displaystyle-\frac{1}{2}\ln{\left((2\pi)^{6}\big{|}\bm{\Pi}_{q}^{{\dagger}}\bm{\tilde{W}}\bm{\tilde{\Sigma}}_{\bm{y}}\bm{\tilde{W}}^{\dagger}\bm{\Pi}_{q}\big{|}\right)}.$

Then, we use the transformation in Eq. (9) to write the log-likelihood as a
function of $\bm{\tilde{t}}$:

$\displaystyle\ln p\left(\bm{\tilde{y}}|\bm{\theta}\right)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\bm{\tilde{t}}_{q}-\bm{\tilde{t}}_{hq}\right)^{{\dagger}}\bm{\tilde{\Sigma}}_{\bm{t}_{q}\bm{t}_{q}}^{-1}\left(\bm{\tilde{t}}_{q}-\bm{\tilde{t}}_{hq}\right)$
(18)
$\displaystyle-\frac{1}{2}\ln{\left((2\pi)^{6}\big{|}\bm{\tilde{\Sigma}}_{\bm{t}_{q}\bm{t}_{q}}\big{|}\right)},$

where we have also defined
$\bm{\tilde{\Sigma}}_{\bm{t}_{q}\bm{t}_{q}}^{-1}\equiv\bm{\Pi}_{p}\bm{\tilde{\Sigma}}_{\bm{t}}^{+}\bm{\Pi}_{p}$,
$\bm{\tilde{t}}_{q}\equiv\bm{\tilde{W}}_{q}\bm{\tilde{y}}$ and
$\bm{\tilde{t}}_{hq}\equiv\bm{\tilde{W}}_{q}\bm{h}$. As a result, decomposing
the data with PCA and keeping only the lowest variance components is
equivalent to restricting the likelihood to its most significant terms (i.e.,
the terms that are the most relevant for GW parameters inference). The data-
driven PCA process outlined here is similar to the TDI analysis, which has
already been shown to be equivalent to an approximation of the likelihood [13,
15, 16]. Depending on the choice of $q$, Eq. (18) can be considered as a
family of approximations for Eq. (16). We have argued that we expect that
small $q$ will provide a sufficient approximation, but the result become exact
in the limit that $q\rightarrow 6p$.

###  II.5 Orthogonalization of aPCI streams with respect to secondary noises

To further simplify the computation of the log-likelihood in Eq (17), it is
more convenient to manipulate streams which are orthogonal with respect to the
noise, including its non-laser part. Therefore, we need a process equivalent
to the generation of A, E, T combinations [8] in TDI. We already have
generalized this operation in a previous work [16], where we simply compute
the eigenvectors of the covariance matrix:

$\bm{\tilde{\Sigma}}_{\bm{t}}(f)=\bm{\Phi}(f)\bm{\Lambda}(f)\bm{\Phi}(f)^{\dagger},$
(19)

where $\bm{\Phi}(f)$ is the matrix of eigenvectors, and $\bm{\Lambda}(f)$ is
the diagonal matrix of eigenvalues. Note that we perform this decomposition
for each frequency. Ultimately, the orthogonal PCA combinations are given by

$\bm{\tilde{t}}_{\perp}(f)=\bm{\Phi}(f)^{\dagger}\bm{\tilde{t}}(f),$ (20)

and their PSDs are given by the diaogonal elements of $\bm{\Lambda}(f)$ . We
have now a process to compute the noise levels necessary to derive the
gravitational strain sensitivities.

##  III Evaluation of aPCI performance

The right quantity to assess the performance of the PCA process to analyse
gravitational-wave measurements is the sensitivity of the PC projections,
which we derive in this section.

###  III.1 Orbit assumptions

We assume that the LISA spacecraft follows Keplerian orbits approximated to
first order in eccentricity, with $e=4.8237\times 10^{-3}$. We assume that the
spacecraft form an equilateral triangle whose center follows a circular orbit
around the Sun, such that its position at any time $t$ is given by

$\displaystyle\bm{r}_{0}=\left(R\cos\phi_{T},\,R\sin\phi_{T},\,0\right)^{\intercal},$
(21)

where $R=1\,\mathrm{au}$ is the distance between the Sun and the center of the
constellation, and $\phi_{T}(t)$ is the flight path angle.

###  III.2 Derivation of gravitational-wave responses

Here we derive the aPCI variables response to gravitational-wave signals,
which is needed to compute their sensitivities. At the time of reception $t$,
the response of gravitational waves in a single link $ij$ is given in
fractional frequency by [24]:

$y^{\mathrm{GW}}_{ij}(t)=\frac{H_{ij}\left(t_{j}-\frac{\bm{k}\cdot\bm{r}_{j}(t_{j})}{c}\right)-H_{ij}\left(t-\frac{\bm{k}\cdot\bm{r}_{i}(t)}{c}\right)}{2\left(1-\bm{k}\cdot\bm{\hat{n}}_{ij}\right)},$
(22)

where we labeled the GW projection onto arm $ij$ as $H_{ij}$, the spacecraft
$i$ position vector as $\bm{r}_{i}$, the unit vector along arm $ij$ as
$\bm{\hat{n}}_{ij}=\left(\bm{r}_{i}-\bm{r}_{j}\right)/\lVert\bm{r}_{i}-\bm{r}_{j}\rVert$
and the time of emission as $t_{j}=t-L_{ij}/c$, which depends on $L_{ij}$, the
effective distance between the sending spacecraft $j$ and the receiving
spacecraft $i$. We compute the vectors $\bm{r}_{i}$ from the orbital motion
postulated in Section III.1. In this equation, we also defined $\bm{k}$ as the
direction of GW propagation, which can be expressed as a function of the
source’s ecliptic latitude $\beta$ and longitude $\lambda$:

$\displaystyle\bm{k}=-\left(\cos\lambda\cos\beta,\,\sin\lambda\cos\beta,\,\sin\beta\right)^{\intercal}.$
(23)

The GW projection can be written as a function of antenna pattern functions
$F_{+}$ and $F_{\times}$ such that

$H_{ij}(\tau)=h_{+}(\tau)F_{+}\left(\psi,\bm{k},\bm{\hat{n}}_{ij}\right)+h_{\times}(\tau)F_{\times}\left(\psi,\bm{k},\bm{\hat{n}}_{ij}\right),$
(24)

where $h_{\alpha}$ is the strain deformation in the source frame for the
polarization mode $\alpha$. The antenna pattern functions depend on the GW
polarization angle $\psi$, the wave propagation vector $\bm{k}$ and the link
vector $\bm{\hat{n}}$ as

$\displaystyle F_{+}\left(\psi,\bm{k},\bm{\hat{n}}\right)$
$\displaystyle=\cos(2\psi)\xi_{+}\left(\bm{k},\bm{\hat{n}}\right)-\sin(2\psi)\xi_{\times}\left(\bm{k},\bm{\hat{n}}\right);$
$\displaystyle F_{\times}\left(\psi,\bm{k},\bm{\hat{n}}\right)$
$\displaystyle=\sin(2\psi)\xi_{+}\left(\bm{k},\bm{\hat{n}}\right)+\cos(2\psi)\xi_{\times}\left(\bm{k},\bm{\hat{n}}\right),$
(25)

where $\xi_{+}$ and $\xi_{\times}$ are modulation functions defined as

$\displaystyle\xi_{+}\left(\bm{k},\bm{\hat{n}}\right)$ $\displaystyle=$
$\displaystyle\left(\bm{u}\cdot\bm{\hat{n}}\right)^{2}-\left(\bm{v}\cdot\bm{\hat{n}}\right)^{2};$
$\displaystyle\xi_{\times}\left(\bm{k},\bm{\hat{n}}\right)$ $\displaystyle=$
$\displaystyle
2\left(\bm{u}\cdot\bm{\hat{n}}\right)\left(\bm{v}\cdot\bm{\hat{n}}\right).$
(26)

The vectors $\bm{u}$ and $\bm{v}$ in the above equation form the source
localization basis $\left(\bm{u},\,\bm{v},\,\bm{k}\right)$, with

$\displaystyle\bm{u}$
$\displaystyle=\left(\sin\lambda,\,-\cos\beta,\,0\right)^{\intercal};$
$\displaystyle\bm{v}$
$\displaystyle=\left(-\sin\beta\cos\lambda,\,-\sin\beta\sin\lambda,\cos\beta\right)^{\intercal}.$
(27)

To write the response as a function of the source frequency, it will be useful
to express the GW strain in terms of its Fourier transform
$\tilde{h}_{\alpha}(f)$:

$h_{\alpha}(t)=\int_{-\infty}^{+\infty}\tilde{h}_{\alpha}(f^{\prime})e^{2i\pi
f^{\prime}t}df^{\prime}.$ (28)

We plug Eq. (28) into Eq. (22) to get

$\displaystyle y^{\mathrm{GW}}_{ij}(t)$
$\displaystyle=\sum_{\alpha=+,\times}\int_{-\infty}^{+\infty}\tilde{h}_{\alpha}(f^{\prime})e^{2i\pi
f^{\prime}t}df^{\prime}\frac{F_{\alpha}\left(\psi,\bm{k},\bm{\hat{n}}_{ij}\right)}{2\left(1-\bm{k}\cdot\bm{\hat{n}}_{ij}\right))}$
$\displaystyle\times\left[e^{-2i\pi
f^{\prime}\left(L_{ij}+\bm{k}\cdot\bm{r}_{j}(t_{j})\right)/c}-e^{-2i\pi
f^{\prime}\bm{k}\cdot\bm{r}_{i}(t)/c}\right].$ (29)

In the following, we assume that the gravitational wave is monochromatic, with
frequency $f$. Choosing a null phase, the Fourier transform of $h_{\alpha}(t)$
can then be written as
$\tilde{h}_{\alpha}(f^{\prime})=\frac{1}{2}A_{\alpha}\delta\left(f^{\prime}-f\right)$
for positive frequencies. Inserting this into Eq. (III.2) gives

$\displaystyle y^{\mathrm{GW}}_{ij}(t,f)$
$\displaystyle=\sum_{\alpha=+,\times}\frac{1}{2}A_{\alpha}e^{2i\pi
ft}\frac{F_{\alpha}\left(\psi,\bm{k},\bm{\hat{n}}_{ij}\right)}{2\left(1-\bm{k}\cdot\bm{\hat{n}}_{ij}\right))}$
$\displaystyle\times\left[e^{-2i\pi
f\left(L_{ij}+\bm{k}\cdot\bm{r}_{j}(t_{j})\right)/c}-e^{-2i\pi
f\bm{k}\cdot\bm{r}_{i}(t)/c}\right].$ (30)

The modes amplitudes $A_{\alpha}$ can be expressed as a function of the
inclination $\iota$ of the system with respect to the line of sight:

$\displaystyle A_{+}$ $\displaystyle=A\left(1+\cos^{2}\iota\right);$
$\displaystyle A_{\times}$ $\displaystyle=-2A\cos\iota,$ (31)

where $A$ is the amplitude of the source, that we consider constant in our
monochromatic assumption.

We can now derive the response of PCA projections to GWs using the time-domain
equivalent of Eqs. (9):

$\bm{t}^{\mathrm{GW}}(f)=\bm{V}^{\dagger}\begin{pmatrix}\bm{y}^{\mathrm{GW}}(t+n_{h}\tau_{s},f)\\\
\vdots\\\ \bm{y}^{\mathrm{GW}}(t-n_{h}\tau_{s},f)\end{pmatrix}^{\dagger},$
(32)

where $\bm{y}^{\mathrm{GW}}(t,f)$ is simply the vector of single-link GW
responses

$\displaystyle\bm{y}^{\mathrm{GW}}(t,f)\equiv\begin{pmatrix}y^{\mathrm{GW}}_{12}(t,f)\\\
y^{\mathrm{GW}}_{23}(t,f)\\\ y^{\mathrm{GW}}_{31}(t,f)\\\
y^{\mathrm{GW}}_{13}(t,f)\\\ y^{\mathrm{GW}}_{21}(t,f)\\\
y^{\mathrm{GW}}_{32}(t,f)\end{pmatrix},$ (33)

whose entries are given by Eq. (III.2).

###  III.3 Computation of sky-averaged sensitivity

In Sec. II.3 we derived the frequency-domain covariance of the aPCI
combinations as a function of the single-link data covariance
$\bm{\tilde{\Sigma}}_{\bm{y}}$, while in Sec. III.2 we derived their response
to gravitational waves. Following Babak et al. [25], we define the
instantaneous strain sensitivity of a given component $m$ as its noise PSD
divided by its response averaged over all sky locations and orientations:

$S_{h,m}(f)=A^{2}\frac{S_{\perp,m}(f)}{\langle\lvert
t^{\mathrm{GW}}_{\perp,m}(f)\rvert^{2}\rangle},$ (34)

where
$\bm{t}^{\mathrm{GW}}_{\perp}\equiv\bm{\Phi}(f)^{\dagger}\bm{t}^{\mathrm{GW}}$
is the GW response in the orthogonalized aPCI combinations,
$S_{\perp,m}(f)\equiv\Lambda_{mm}(f)$ is the noise PSD given by the diagonal
elements of matrix $\Lambda(f)$ in Eq. (19), and the brackets denote the
averaging over sky location angles, polarization and inclination:

$\displaystyle\langle\lvert
t^{\mathrm{GW}}_{\perp,m}\rvert^{2}\rangle=\frac{1}{4\pi}\int_{\bm{k}}\frac{1}{2\pi}\int_{0}^{2\pi}\frac{1}{2}\int_{-1}^{+1}\lvert
t^{\mathrm{GW}}_{\perp,m}\rvert^{2}d^{3}\bm{k}d\psi d\cos\iota.$ (35)

In this definition, the sensitivity corresponds to the level of strain noise
within an infinitesimal frequency band $df$ (as stated in, e.g., [26]). Note
that we include the factor $A^{2}$ in front of Eq. (34) to cancel the one
present in the denominator via Eqs. (III.2), making the sensitivity
independent of the amplitude, as it precisely determines the variance of $A$.
We can sum the sensitivity over a subset of components between
$m_{\mathrm{min}}$ and $m_{\mathrm{max}}$:

$S_{h}(f)=\left[\sum_{m=m_{\mathrm{min}}}^{m_{\mathrm{max}}}S^{-1}_{h,m}(f)\right]^{-1},$
(36)

where the terms in the sum are the inverse sensitivities of each component,
which can be viewed as their SNR per unit time and unit amplitude. In the
following, we compute this sensitivity on a given example and compare it with
its TDI analog.

To evaluate the sensitivity in Eq. (34), we need an analytical model for the
PSD of the single-link data, at least under some assumptions. This will then
allow us to derive an expression for the PSD of the PCA combinations
$S_{\perp,m}(f)$. We assume that LISA’s single-link measurements include two
main sources of noise: a large term due to laser frequency fluctuations and a
much smaller term due to other instrumental noises. We assume that the signals
measured by the two optical benches on-board spacecraft 1 are given by the
equations

$\displaystyle s_{12}(t)$ $\displaystyle=$ $\displaystyle
p_{21}(t-c^{-1}L_{12})-p_{12}(t)+n_{12}(t)$ $\displaystyle s_{13}(t)$
$\displaystyle=$ $\displaystyle p_{31}(t-c^{-1}L_{13})-p_{13}(t)+n_{13}(t),$
(37)

where $p_{ij}$ are laser noises, $n_{ij}$ are secondary noises, and $c$ is the
speed of light, labeling by the convention noted in Section II.1. We obtain
the other measurements on-board spacecraft 2 and 3 by permuting the indices
(modulo 3) in Eq. (III.3). To simplify the analysis, we assume that the two
lasers onboard each spacecraft $i$ are identical, so that $p_{ij}=p_{ik}$.
This amounts to assuming that the $\bm{s}_{ij}$ are equal to the vectors of
intermediate variables $\bm{\eta}_{ij}$ in the TDI pre-processing (see, e.g.,
Ref. [27]). In this model, the measured laser noises can be seen as the result
of the mixing (and delaying) of a vector of laser noise sources
$\bm{p}=\left(\bm{p}_{12}^{\intercal},\,\bm{p}_{23}^{\intercal},\,\bm{p}_{31}^{\intercal}\right)^{\intercal}$.
We consider the $6N$ column-vector version of the $N\times 6$ matrix $\bm{Y}$
in Eq. (1), which stacks all the measurements as

$\bm{y}\equiv\left(\bm{s}_{12}^{\intercal},\,\bm{s}_{23}^{\intercal},\,\bm{s}_{31}^{\intercal},\,\bm{s}_{13}^{\intercal},\,\bm{s}_{21}^{\intercal},\,\bm{s}_{32}^{\intercal}\right)^{\intercal}.$
(38)

This mixing is represented by a $6N\times 3N$ design matrix $\bm{M}$, such
that the model for $\bm{y}$ is

$\bm{y}=\bm{M}\bm{p}+\bm{n},$ (39)

where we set
$\bm{n}\equiv\left(\bm{n}_{12}^{\intercal},\,\ldots,\,\bm{n}_{32}^{\intercal}\right)^{\intercal}$,
using the same ordering as in Eq. (38). Using this assumption, we have

$\bm{M}=\begin{pmatrix}-\bm{1}&\bm{D}_{12}&\bm{0}\\\
\bm{0}&-\bm{1}&\bm{D}_{23}\\\ \bm{D}_{31}&\bm{0}&-\bm{1}\\\
-\bm{1}&\bm{0}&\bm{D}_{13}\\\ \bm{D}_{21}&-\bm{1}&\bm{0}\\\
\bm{0}&\bm{D}_{32}&-\bm{1}\end{pmatrix},$ (40)

where $\bm{D}_{lm}$ represents the fractional delay operator acting on any
discretized signal $x(t_{n})$ as $\bm{D}_{lm}\bm{x}|_{n}\approx
x\left(t-L_{lm}/c\right)$.

In the following, we assume that $\bm{p}_{ij}$ and $\bm{n}_{ij}$ are zero-mean
stationary Gaussian noises vectors of size $N$ described by PSDs $S_{p,ij}(f)$
and $S_{n,ij}(f)$, respectively. We further assume that these entries are
uncorrelated. Since the laser noise level is much larger than the other
noises, we have $S_{p,ij}(f)\gg S_{n,ij}\,\forall i,j$. For this reason, we
will refer to $\bm{n}$ as the vector of secondary noises.

###  III.4 Practical implementation

For a number of aPCI variables, we compute the orientation and sky-averaged
sensitivity following Eqs. (36). We perform the integration with respect to
the polarization angle analytically, which is equivalent to setting
$\psi_{\mathrm{eq}}=\pi/4$ (as it corresponds to equal contributions between
the $+$ and $\times$ modes). We also use an analytical result for the
integration over inclination, as it amounts to multiplying the squared
response by $4/5$ [28]. Finally, we average over sky locations numerically by
randomly drawing 3000 values of the source ecliptic latitude $\beta$ and
longitude $\lambda$ from a uniform distribution in the interval
$\left[0,\,\pi\right]\times\left[0,\,2\pi\right]$.

Suppose we have computed the SVD decomposition. Among the $6p$ PCs provided by
the algorithm, let us consider only $q$ of them. For each sky location
$(\lambda,\,\beta)$ and each GW frequency $f$, we compute the response by
following 3 steps:

  1. 1.

Compute the arm responses $y^{\mathrm{GW}}_{ij}$ using Eq. (III.2);

  2. 2.

Form the $q$ PCA projections $\bm{t}^{\mathrm{GW}}_{m}$ with Eq. (32);

  3. 3.

Orthogonalize the projections using Eq. (20);

To get a baseline for comparison, we also compute the sensitivity of first-
generation TDI, as defined in, e.g., Ref. [29]. The computation steps follow a
similar scheme as above. The first step is the same: we calculate the arm
response to harmonic GWs exactly as before. Step 2 is replaced by the
computation of the Michelson TDI combinations $X$, $Y$ and $Z$. This is done
by applying the TDI 1.5 transfer function to the single-link response vector
$\bm{t}^{\mathrm{GW}}_{\mathrm{TDI}}=\bm{V}_{\mathrm{TDI}}^{\dagger}\bm{y}^{\mathrm{GW}}$,
where $\bm{V}_{\mathrm{TDI}}$ can be written as

$\displaystyle\bm{V}_{\mathrm{TDI}}$
$\displaystyle\equiv\begin{pmatrix}D_{13}D_{31}-1&\left(1-D_{23}D_{32}\right)D_{21}&0&\\\
0&D_{21}D_{12}-1&\left(1-D_{31}D_{13}\right)D_{32}\\\
\left(1-D_{12}D_{21}\right)D_{13}&0&D_{32}D_{23}-1\\\
1-D_{12}D_{21}&0&\left(D_{32}D_{23}-1\right)D_{31}\\\
\left(D_{13}D_{31}-1\right)D_{12}&1-D_{23}D_{32}&0\\\
0&\left(D_{21}D_{12}-1\right)D_{23}&1-D_{31}D_{13}\\\ \end{pmatrix}.$ (41)

Then, Step 3 is replaced by forming the $A$, $E$ and $T$ combinations from
$X$, $Y$ and $Z$, as defined in Ref. [4]. Finally, we compute the TDI
sensitivity over all channels with a calculation similar to Eq. (36).

##  IV Numerical simulations and results

In this section, we verify that the aPCI algorithm yields a significant noise
subtraction. Then we check the PSD model derived in Section II.3 against
realizations of the noise. Finally, we compute averaged sensitivities to
gravitational waves, and compare them with Sagnac TDI observables A, E and T.

###  IV.1 Laser noise subtraction

We test the aPCI algorithm on simulated LISA data streams. We adopt the same
assumptions as in Ref. [16]. To generate the noises, we assume that the laser
sources of two optical benches belonging to the same spacecraft are the same,
and that the light travel time delays between spacecraft are constant and
equal to $L/c$, where $L=2.5\times 10^{9}$ m is the mean inter-spacecraft
distance.

We simulate the science measurement time series $s_{ij}$ from Eq. (III.3). We
generate noises from zero-mean Gaussian stationary processes in the frequency
domain described by PSDs $S_{p}(f)$ and $S_{n}(f)$ given by:

$\displaystyle S_{p}(f)$
$\displaystyle=\left(\frac{28.2\,\mathrm{Hz.Hz^{-1/2}}}{\nu_{0}}\right)^{2};$
$\displaystyle S_{n}(f)$ $\displaystyle=\frac{S_{a}(f)}{\left(2\pi
fc\right)^{2}}+\left(\frac{2\pi f}{c}\right)^{2}S_{\mathrm{oms}}(f),$ (42)

where $\nu_{0}=c\cdot(1064\times 10^{-9}\mathrm{m})^{-1}$ is the central
frequency of the lasers, $S_{a}(f)$ is the test-mass acceleration noise PSD
(in $\rm m^{2}s^{-4}Hz^{-1}$) and $S_{\mathrm{oms}}(f)$ is the PSD
representing all other pathlength noises in the measurement (in $\rm
m^{2}Hz^{-1}$). Their expressions are given in Appendix A, where we adopted
assumptions made for the LISANode simulator [30]. Likewise, we apply light
travel time delays in the frequency domain. To prevent any artificial
periodicity, we generate $2N$ samples corresponding to Fourier frequencies
$\pm kf_{s}/(2N)$ before transforming the generated noise to the time domain,
and discarding the second half of the data to get a time series of size $N$.

We choose the duration of the time series to be 12 hours, with a sampling rate
of $f_{s}=2$ Hz, which allows us to cover the frequency band between 0.1 mHz
and 1 Hz. We choose the half-stencil size $n_{h}$ so that the stencil spans
the duration of one complete turn around the constellation $n_{h}=\lfloor
3Lf_{s}/c\rfloor=25$. Note that we did not find significant changes in the
results presented in the following when choosing larger values for $N$ and
$n_{h}$. A full assessment of the impact of data and stencil sizes is left for
future work. We form the data matrix following Eq. (2) and compute its PCA as
in Eq. (4) with the algorithm provided by the Python package Scikit-learn
[31]. In this example, the decomposition yields 306 PCA components.

We can analyze the relative amount of noise present in each component $m$ by
calculating the variance quantity $S^{2}_{m,m}/(N-1)$, where $\bm{S}$ is the
singular value matrix in Eq. (3). We report these values for all $m$ on Fig.
(1). We see that the variance decreases with increasing values of $m$, with a
drop spanning 15 orders of magnitudes. This plot is a good hint about the
ability of the algorithm to separate large noise from low noise components.

Figure 1: Variance associated to PCA components in decreasing order. Instead
of selecting the largest variance components like in classic PCA, here we are
interested in the lower end of the plot, i.e. the less noisy components at the
right-hand side of the red dashed line.

Now let us order the aPCI variables from the lowest variance to the largest
variance, so that index $m=1$ corresponds to the rightmost end of Fig. 1, and
$m=6p$ corresponds to the leftmost end. We project the data vector $\bm{y}$
onto the last 3 singular vectors using Eq. (4), obtaining 3 data projections
whose periodograms are plotted in Fig. 2 in red for $m=1$, in orange for
$m=2$, in blue for $m=3$. We define the cross-periodogram of two time series
$x$ and $y$ of size $N$ as

$\displaystyle P_{xy}(f)=\tilde{x}(f)\tilde{y}^{\ast}(f),$ (43)

where $\tilde{x}$ refers to the windowed discrete Fourier transform of $x$
normalized as

$\displaystyle\tilde{x}(f)\equiv\sqrt{\frac{2\tau_{s}}{N}}\sum_{n=0}^{N-1}w_{n}x_{n}e^{-2i\pi
fn\tau_{s}}.$ (44)

As a comparison, we plot the periodogram of link 12 (gray curve). The
difference between the gray and colored curves shows that the noise is
significantly suppressed by several orders of magnitude. In addition, we
verify that the noise model derived in Eq. (11) accurately predicts the noise
levels in the periodograms. To do so, we include the added contributions of
laser noise and secondary noises to the aPCI residuals.

Figure 2: Periodograms (light colors) of the three least noisy PCA projections
as a function of frequency. We also plot the theoretical PSDs (dark colors)
obtained with Eq. (11), expressed in units of relative frequency deviation. To
compare with the input noise, we plot in gray the periodogram of the raw
measurement from laser link 12, which is overwhelmed by laser noise.

Besides, we analyze the correlations among the aPCI variables by comparing the
cross-periodograms $P_{{t_{m}}{t_{l}}}(f)$ with the cross spectral density
(CSD) model provided by the non-diagonal terms of the covariance in Eq. (11).
As an example, we plot in Fig. 3 the CSD of the pairs $(m,l)=(2,3)$ in blue,
$(3,1)$ in brown, $(2,1)$ in turquoise. They are all consistent with the
cross-periodogram (lighter colors) averaged over 25 realizations of the noise.
We see that the aPCI variables are strongly correlated, with CSD levels
comparable to the PSDs. There is no reason why this should not be the case, as
the PCA decomposition is built mainly from the dominant laser noise
statistics. In addition, the matrix $\bm{X}$ we decompose is formed from the
data itself, not from the actual covariance of the data $\bm{\Sigma}_{y}$. As
a result, the process does not give exactly orthogonal data streams by
construction. That is why an extra diagonalization step can be applied, as
described in Section II.5, to account for secondary noises.

Figure 3: Averaged cross-periodograms of the 3 aPCI variables of lowest
variance (thin light lines), along with the corresponding analytical CSD model
(darker solid lines). Correlations among components are comparable to their
individual PSDs (plotted in Fig. 2).

While the noise levels in the aPCI streams $\bm{t}_{m}$ give us an indication
about the amount of laser noise that is suppressed, there is little sense in
comparing them directly to the TDI PSDs. It is rather preferable to compute GW
sensitivities, which we do in the next section.

###  IV.2 Comparison of aPCI and TDI sensitivities

In this section, we arrange the aPCI variables in increasing order of
variance, i.e. we start from the lower end of the curve in Fig. 1. We define
the sensitivity of the $q$ variables of lowest variance as the sensitivity
obtained by accumulating the SNRs of their orthogonal transformations, as
outlined in the steps of Section III.4.

We plot the sensitivities for a number of components ranging from $q=1$ to
$q=6$ in Fig. 4. Each blue dashed line corresponds to a value of $q$, with
darker colors for larger $q$. To better distinguish the levels, we fill the
areas in between two levels with related shades of blue. For example, the
uppermost, lightest curve is the sensitivity obtained with a single aPCI
variable ($q=1$), i.e. the one with lowest variance. The curve that comes just
below ($q=2$) is the sensitivity obtained by using the two variables of lowest
variance, and so on. As expected, we can see that the larger the number of
components we include, the better the sensitivity.

Figure 4:  Comparison of aPCI (blue dashed lines) and TDI (thick orange line)
sensitivities summed over combinations. The lightest, uppermost blue dashed
line corresponds to the sensitivity obtained with a single aPCI component,
chosen to be the one with the smallest singular value. The darker curve just
below corresponds to the accumulated sensitivity obtained by including an
additional component (corresponding to the next smallest singular value).
Darker blue dashed curves beneath are constructed similarly, by adding more
components. The lowermost dashed blue curve is the sensitivity reached with 6
components, which almost perfectly fits the orange line representing the
accumulated AET TDI sensitivity.

We also plot in Fig. 4 the summed sensitivity of TDI channels A, E and T with
a thick, solid orange line. Remarkably, the figure shows that we reach the
same sensitivity as orthogonal TDI with 6 aPCI variables, as we can see from
the superimposition of the lowest dashed blue line and the orange solid line.
We have checked that adding more aPCI variables does not improve the
sensitivity any further, as they are too much dominated by laser noise.

###  IV.3 SNR distribution and statistical behavior

Fig. 4 suggests that we need 6 aPCI variables to reach the exact same
sensitivity as the 3 orthogonal TDI combinations. To analyse the extent to
which this is true, we can assess the fraction of sensitivity, or rather
equivalent squared SNR, contained in each component. We do so by defining the
squared SNR as simply the inverse sensitivity integrated over the frequency
band:

$\mathrm{\overline{SNR}}_{m}^{2}\equiv\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}S^{-1}_{h,m}(f)df,$
(45)

where $S_{h,m}(f)$ is given by Eq. (34), $f_{\mathrm{min}}=10^{-4}$ Hz and
$f_{\mathrm{max}}=1$ Hz. We are not using the term “SNR” in the usual way
here, as the inverse sensitivity is just a proxy for the SNR of a
monochromatic source. It can be thought of as a SNR per unit time and unit
strain. Then, the squared SNR fraction of component $m$ relative to TDI is
simply $\mathrm{\overline{SNR}}_{m}^{2}$ divided by the total SNR2 achieved by
variables A, E, T.

Besides, since the PCA is performed on noisy data, singular vectors are random
variables with a variability depending on the particular realization of the
noise present in $\bm{y}$. Therefore, to study how the sensitivity is
distributed between different components, we need to consider several noise
realizations and examine this distribution on average.

For this purpose, we generate 100 realizations of the 12-hour long vector
$\bm{y}$, and compute the aPCI variables for each of them, along with their
theoretical PSDs. We gather the results in Fig. 5, where we plot the mean PSDs
for the three variables with lowest singular values (the same as in Fig. 2),
along with their 68% confidence intervals. This plot reveals the significant
variability of the relative noise levels within each component.

Figure 5: Impact of noise realizations on the three least noisy aPCI variables
noise PSDs. The solid lines correspond to the mean PSDs accross 100 noise
realizations, while the colored areas around each line represent the
1-$\sigma$ region. The indexing $m$ is ordered along increasing variance, from
$m=1$ (lowest variance) to $m=3$ (third lowest variance).

For each noise realization, we compute the sensitivity of the 6 lowest
variance aPCI variables, and the SNR2 fraction (relative to TDI) carried by
each variable. In Fig. 6, we plot the distribution of cumulative SNR2 fraction
as a function of the number of components in the form of box plots (blue),
using the seaborn package [32]. We range the components from largest to
smallest SNR2. As a reference, we plot the TDI combinations’ cumulative SNR2
when adding A, E and T in this order (horizontal yellow dashed lines).

Figure 6:  Distributions of cumulative fraction of SNR2 (in percent) as a
function of the number of variables, ordered in decreasing SNR2 fractions. The
blue boxes correspond to aPCI variables and are made of 3 line showing the 3
quartiles, while the vertical lines on both sides extend to the extreme values
of the distribution. The yellow horizontal lines show the cumulative SNR2
fraction of TDI variables A, E, T. The distributions are obtained from 100
realizations of the data vector $\bm{y}$.

We see that in average, most of the SNR2 (about 94%) is concentrated in the 3
most sensitive aPCI variables (boxes on the left-hand side of the figure).
While the first component always gather about 50% of the SNR2, the heights of
the boxes indicate a significant spread of the SNR2 fraction around the median
when adding components 2 and 3. The 3 last components (4, 5, 6) have marginal
SNR concentrations, since their total fraction is about 6% on average. It
should be noted that the last box is almost flat, which means that the total
SNR (or sensitivity) hardly changes from one realization to another. It
fluctuates around 99.99751(7)% of the TDI total SNR. Therefore, the SNR
contained in the last 6 aPCI variables is virtually equal to the cumulative
SNR of the orthogonal TDI variables. The model independence results in
negligible SNR loss, however the gravitational wave signal is diluted across
more than the traditional 3 TDI variables.

###  IV.4 Impact of gravitational waves

The PCA process outlined in this work acts directly on the measurements, hence
the combinations that we derive depend on the content of the data. A natural
question that arises is the impact of the presence of gravitational-wave
signals onto the PC vectors, and thus on the overall sensitivity of the aPCI
variables.

The core property which distinguishes the laser noise as we isolate it for
cancellation is just that it is loud. If a sufficiently large gravitational-
wave signal is present, we should expect that our process will cancel that as
well, but how loud does the signal have to be for such losses to be
significant?

Here, we take a first step in answering that question by simulating a quasi-
monochromatic signal from a ultra-compact galactic binary, whose amplitude $A$
is allowed to vary. For the purpose of this assessment, we consider arbitrary
and unphysical amplitude values chosen to give the source particular SNRs,
ranging from 1 to $10^{9}$. We choose a fixed GW frequency equal to 1 mHz, as
well as a fixed sky location. For each value of the amplitude, we create a
data set by injecting the GW signal into the phasemeter measurements, in
addition to a noise realization which remains the same for all sets. The
duration of the time series is the same as in the previous sections, i.e. 12
hours.

For each measurement set containing the GW signal, we build the data matrix,
compute its PCA decomposition, and derive the corresponding averaged
sensitivity that we denote by $S^{\mathrm{GW}}_{h}$. We compare these
sensitivities to the signal-free sensitivity (derived with data containing
only noise) denoted by $S_{h}$ by computing their relative error maximized
over the frequency range:

$\gamma=\underset{f}{\mathrm{max}}\left\\{\frac{\left|S^{\mathrm{GW}}_{h}(f)-S_{h}(f)\right|}{S_{h}(f)}\right\\}.$
(46)

We report the results in Fig. 7 where we plot $\gamma$ as a function of the
source SNR. The relative error in sensitivity remains below 1 % for SNR values
below or equal to $10^{7}$. Beyond $\mathrm{SNR}=10^{8}$, the error becomes
significant, reflecting the fact that the PCA decomposition starts to be
affected by the presence of the injected signal. This demonstrates the
robustness of the method in the presence of deterministic signals, as the
chosen SNRs are extremely large relative to typical GW sources, especially
given the short observation time we consider. As a comparison, the total GW
SNR accumulated over the mission duration should be of order $\sim 10^{4}$,
which is way below the threshold where it starts to affect the PCA
decomposition.

Figure 7: Impact of a 1 mHz monochromatic GW source onto the aPCI sensitivity:
relative difference of the sensitivity obtained with the presence of the
source relative to the source-free sensitivity.

##  V Discussion

We outlined a method designed to process space-borne interferometer telemetry
and minimize laser frequency noise for data analysis purposes. This method,
aPCI, is completely blind: it does not rely on any model describing the
measurements. The only implicit assumption is that there exists some unknown
linear relationship over short times between the various measurements. We
begin by building a data matrix by stacking integer-shifted versions of all
single-link measurements; then we compute its PCA decomposition. Projecting
the data onto the components with lowest variances yields data streams that
are almost free from laser noise, and are directly usable for gravitational-
wave searches. We demonstrate its equivalence with Bayesian inference based on
TDI. This approach is a robust, complementary alternative to existing TDI
techniques which depend on a specific model of the data (the laser noise
equations) and on physical parameters (the time delays). This model-
independent technique could be an important tool to ensure the robustness of
the data analysis in missions such as LISA.

We demonstrated the effectiveness of the method to remove laser noise using
simulated LISA single-link measurements, assuming fixed armlengths. We showed
that the accumulated GW sensitivity of the 6 lowest-variance PCA projections
is always effectively equal to the sensitivity of first-generation Michelson
TDI combinations, and that the major part of the SNR concentrates in the 2
most sensitive channels.

This work lays the foundation for future investigations to build a fully data-
driven approach to space-based interferometer data-analysis. Many aspects
remain to be studied in more detail. First, the ability of the method to
process measurements obtained with time-varying armlengths and other realistic
details should be tested. Second, even though we showed that the method is
robust against loud monochromatic gravitational-wave signals, we should extend
this assessment to the presence of various sources across the full sensitive
bandwidth, as expected in LISA. Third, one should investigate how aPCI can be
further improved towards a better compression of the GW information, using
orthogonality with respect to the sensitivity. Fourth, similar to TDI, aPCI
would have to face data artifacts like instrumental transients and data gaps
expected in LISA-like detector measurements. We should study the impact of
these perturbations on aPCI performance and test its combination with
mitigation techniques such as imputation [33]. Finally, an important step to
take is to fully develop and demonstrate Bayesian inference analysis using
this approach.

##  Appendix A Pathlength noises

We provide here the analytical forumlas for the secondary noise PSDs used in
this study. The acceleration noise PSD is given by

$\displaystyle
S_{a}(f)=a_{\mathrm{TM}}^{2}\left[1+\left(\frac{f_{1}}{f}\right)^{2}\right]\left[1+\left(\frac{f}{f_{2}}\right)^{4}\right]$
(47)

where $a_{\mathrm{TM}}=3\times 10^{-15}\,\mathrm{ms^{-2}}$, $f_{1}=4\times
10^{-4}$ Hz and $f_{2}=8\times 10^{-3}$ Hz. The OMS noise PSD is

$\displaystyle
S_{\mathrm{oms}}(f)=a_{\mathrm{oms}}^{2}\left[1+\left(\frac{f_{3}}{f}\right)^{4}\right],$
(48)

where $a_{\mathrm{oms}}=15\times 10^{-12}\,\mathrm{mHz^{-1/2}}$ and
$f_{3}=2\times 10^{-3}$ mHz.

###### Acknowledgements.

We would like to thank Tyson Littenberg and Jessica Page for their interesting
feedback. We also thank Martin Staab for his thorough reading and remarks.
This work is partly supported by an appointment to the NASA Postdoctoral
Program at the Goddard Space Flight Center, administered by Universities Space
Research Association under contract with NASA.

## References

  * [1] Karsten Danzmann, Pau Amaro-Seoane, Heather Audley, Stanislav Babak, John Baker, et al.  LISA laser interferometer space antenna.  arXiv, pages 1–41, 2017.
  * [2] V. K. Milyukov.  TianQin Space-Based Gravitational Wave Detector: Key Technologies and Current State of Implementation.  Astronomy Reports, 64(12):1067–1077, 2020.  Publisher: Pleiades journals.
  * [3] Massimo Tinto and J. W. Armstrong.  Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation.  Physical Review D - Particles, Fields, Gravitation and Cosmology, 59(10):1–11, 1999.  Publisher: American Institute of Physics Inc.
  * [4] Massimo Tinto and Sanjeev V. Dhurandhar.  Time-Delay Interferometry.  Living Reviews in Relativity, 17(1):6, 2014.  Publisher: Springer International Publishing.
  * [5] F. B. Estabrook, Massimo Tinto, and J. W. Armstrong.  Time-delay analysis of LISA gravitational wave data: Elimination of spacecraft motion effects.  Physical Review D - Particles, Fields, Gravitation and Cosmology, 62(4):8, 2000.  Publisher: American Physical Society.
  * [6] Neil J. Cornish and Ronald W. Hellings.  The Effects of Orbital Motion on LISA Time Delay Interferometry.  Classical and Quantum Gravity, 20(22):4851–4860, 2003.
  * [7] Daniel A. Shaddock, Massimo Tinto, Frank B. Estabrook, and J. W. Armstrong.  Data combinations accounting for LISA spacecraft motion.
  * [8] Thomas A. Prince, J. W. Armstrong, Massimo Tinto, and Shane L. Larson.  LISA optimal sensitivity.  Physical Review D - Particles, Fields, Gravitation and Cosmology, 66(12), 2002.
  * [9] Jean Baptiste Bayle, Marc Lilley, Antoine Petiteau, and Hubert Halloin.  Effect of filters on the time-delay interferometry residual laser noise for LISA.  Physical Review D, 99(8), 2019.
  * [10] Olaf Hartwig and Jean Baptiste Bayle.  Clock-jitter reduction in LISA time-delay interferometry combinations.  arXiv, May 2020.
  * [11] Jean-Baptiste Bayle, Olaf Hartwig, and Martin Staab.  Adapting time-delay interferometry for LISA data in frequency.  2021\.
  * [12] Martina Muratore, Daniele Vetrugno, and Stefano Vitale.  Revisitation of time delay interferometry combinations that suppress laser noise in LISA.  Classical and Quantum Gravity, 37(18):185019, aug 2020.
  * [13] J. D. Romano and G. Woan.  Principal component analysis for LISA: The time delay interferometry connection.  Physical Review D - Particles, Fields, Gravitation and Cosmology, 73(10), 2006.
  * [14] Michele Dawn Leighton.  A principal component approach to space-based gravitational wave astronomy.  Technical report, 2016.
  * [15] Michele Vallisneri, Jean Baptiste Bayle, Stanislav Babak, and Antoine Petiteau.  TDI-infinity: time-delay interferometry without delays.  arXiv, 2020.
  * [16] Quentin Baghi, James Ira Thorpe, Jacob Slutsky, and John Baker.  Statistical inference approach to time-delay interferometry for gravitational-wave detection.  Physical Review D, 103(4), 2021.
  * [17] Massimo Tinto, Sanjeev Dhurandhar, and Prasanna Joshi.  Matrix Representation of Time-Delay Interferometry.  2021\.
  * [18] Jean-Baptiste Bayle, Michele Vallisneri, Stanislav Babak, and Antoine Petiteau.  On the matrix formulation of time-delay interferometry.  2021\.
  * [19] Massimo Tinto, Sanjeev Dhurandhar, and Prasanna Joshi.  Reply to the Bayle et al. gr-qc document dated June 7, 2021.  2021\.
  * [20] Yan Wang, Gerhard Heinzel, and Karsten Danzmann.  First stage of LISA data processing: Clock synchronization and arm-length determination via a hybrid-extended Kalman filter.  Physical Review D, 90(6):064016–064016, September 2014.  Publisher: American Physical Society.
  * [21] Massimo Tinto, Michele Vallisneri, and J. W. Armstrong.  Time-delay interferometric ranging for space-borne gravitational-wave detectors.  Physical Review D, 71(4):041101, February 2005.  Publisher: American Physical Society.
  * [22] Jessica Page and Tyson Littenberg.  Bayesian Time Delay Interferometry.  arXiv:2107.08174 [gr-qc], July 2021.  arXiv: 2107.08174.
  * [23] D A Shaddock, B Ware, R E Spero, and M Vallisneri.  Postprocessed time-delay interferometry for LISA.  Physical Review D - Particles, Fields, Gravitation and Cosmology, 70(8):081101–081101, 2004.
  * [24] Frank B. Estabrook and Hugo D. Wahlquist.  Response of Doppler spacecraft tracking to gravitational radiation.  General Relativity and Gravitation, 6(5):439–447, October 1975\.
  * [25] Stanislav Babak, Martin Hewitson, and Antoine Petiteau.  LISA Sensitivity and SNR Calculations.  arXiv:2108.01167 [astro-ph, physics:gr-qc], August 2021.  arXiv: 2108.01167.
  * [26] Michele Vallisneri and Chad R. Galley.  Non-sky-averaged sensitivity curves for space-based gravitational-wave observatories.  Classical and Quantum Gravity, 29(12):124015, June 2012.  Publisher: IOP Publishing.
  * [27] Markus Otto, Gerhard Heinzel, and Karsten Danzmann.  TDI and clock noise removal for the split interferometry configuration of LISA.  Classical and Quantum Gravity, 29(20):205003–205003, October 2012\.  Publisher: IOP Publishing.
  * [28] Travis Robson, Neil J Cornish, and Chang Liug.  The construction and use of LISA sensitivity curves.  Classical and Quantum Gravity, 36(10):105011–105011, May 2019.  Publisher: IOP Publishing.
  * [29] Massimo Tinto, Frank B. Estabrook, and adn J. W. Armstrong.  Time Delay Interferometry with Moving Spacecraft Arrays.  Physical Review D, 69(8):082001, April 2004.  arXiv: gr-qc/0310017.
  * [30] Jean-Baptiste Bayle.  Simulation and Data Analysis for LISA (Instrumental Modeling, Time-Delay Interferometry, Noise-Reduction Performance Study, and Discrimination of Transient Gravitational Signals).  Theses, Université de Paris ; Université Paris Diderot ; Laboratoire Astroparticules et Cosmologie, October 2019.
  * [31] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, et al.  Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research, 12(85):2825–2830, 2011.
  * [32] Michael L. Waskom.  seaborn: statistical data visualization.  Journal of Open Source Software, 6(60):3021, 2021.
  * [33] Quentin Baghi, James Ira Thorpe, Jacob Slutsky, John Baker, Tito Dal Canton, et al.  Gravitational-wave parameter estimation with gaps in LISA: A Bayesian data augmentation method.  Physical Review D, 100(2):022003–022003, July 2019.  Publisher: American Physical Society.


[d]Carlos Tamarit

# The limits of the strong $CP$ problem

Wen-Yuan Ai  Juan S. Cruz  Björn Garbrecht

###### Abstract

While $CP$ violation has never been observed in the strong interactions, the
QCD Lagrangian admits a $CP$-odd topological interaction proportional to the
so called $\theta$ angle, which weighs the contributions to the partition
function from different topological sectors. The observational bounds are
usually interpreted as demanding a severe tuning of $\theta$ against the
phases of the quark masses, which constitutes the strong $CP$ problem. Here we
report on recent challenges to this view based on a careful treatment of
boundary conditions in the path integral and of the limit of infinite
spacetime volume, which leads to $\theta$ dropping out of fermion correlation
functions and becoming unobservable, implying that $CP$ is preserved in QCD.

KCL-PH-TH/2022-31,TUM-HEP-1403-22

##  1 Introduction

The Lagrangian of quantum chromodynamics (QCD) with gauge fields $A_{\mu}$ in
the Lie Algebra of SU(3) and $N_{f}$ flavours of Dirac fermions $\psi_{i}$
admits $CP$-noninvariant terms. After a rotation to Euclidean spacetime, used
throughout the rest of the paper, these terms read

$\displaystyle{\cal
L}\supset\sum_{i=1}^{N_{f}}\bar{\psi}_{i}\left(\mathfrak{m}_{i}P_{\rm
R}+\mathfrak{m}_{i}^{*}P_{\rm L}\right)\psi_{i}-\frac{{\rm
i}}{16\pi^{2}}\,\theta\,{\rm tr}\,F\widetilde{F}.$ (1)

Above, $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}-{\rm
i}[A_{\mu},A_{\nu}]$ is the non-Abelian field strength, while
$\tilde{F}_{\mu\nu}=\frac{1}{2}\,{\epsilon_{\mu\nu\rho\sigma}}F_{\rho\sigma}$.
$P_{\rm R/L}$ are right/left handed projectors, and $\mathfrak{m}_{i}\equiv
m_{i}{\rm e}^{{\rm i}\alpha_{i}},m_{i},\alpha_{i}\in\mathbb{R}$ are complex
fermion masses. $CP$-violating physical observables may be sensitive to
combinations of the $CP$-odd parameters $\alpha_{n},\theta$ that remain
invariant under chiral field redefinitions. As follows from chiral anomalies,
the relevant combination is

$\displaystyle\bar{\theta}=\theta+\sum_{i}\alpha_{i}\equiv\theta+\bar{\alpha}.$
(2)

The $\theta$-dependent contribution in Eq. (1) can be seen to be a total
derivative. As such, it is a boundary term whose effects can only be captured
with nonperturbative computations. Moreover, whenever the fields go to pure
gauge configurations at the boundary $\partial\Omega$ of the four-dimensional
spacetime $\Omega$, the $\theta$-term is proportional to an integer $\Delta n$
known as “topological charge”,

$\displaystyle\left.A_{\mu}\right|_{\partial\Omega}={\rm
i}{}U(x)\partial_{\mu}U^{\dagger}(x),\,U(x)\in{\rm
SU}(3)\Rightarrow\int_{\Omega}d^{4}x\frac{1}{16\pi^{2}}\,\theta\,{\rm
tr}\,F\widetilde{F}=\theta\Delta n,\,\Delta n\in\mathbb{Z}.$ (3)

It is known that semiclassical expansions about saddle points of the Euclidean
QCD action known as instantons [5, 8] are sensitive to boundary effects, and
hence it is generally expected that nonzero values of the $\bar{\theta}$ in
Eq. (2) will induce nonvanishing predictions for $CP$-violating observables
such as the neutron electric dipole moment $d_{n}$. However, experiments
constrain the latter to [3]

$\displaystyle d_{n}<1.28\times 10^{-26}\,{e\cdot{\rm cm}}.$ (4)

With traditional theoretical calculations predicting
$d_{n}\propto\bar{\theta}$, the stringent bounds on $d_{n}$ are thought to
imply a severe tuning of $\bar{\theta}<10^{-10}$, which is known as the strong
$CP$ problem [1, 2].

Here we report on recent calculations that challenge the above conclusion [7],
and show a lack of $CP$ violation in the strong interactions with and without
relying on instantons methods. For brevity, in these proceedings we will not
cover the instanton-based calculations and will focus instead on arguments
using a minimal set of assumptions related to mathematical consistency in the
choice of boundary conditions, cluster decomposition [4], and the index
theorem [6].

##  2 The importance of boundary conditions in the path integral

The Euclidean QCD partition function is meant to correspond to a transition
amplitude from the vacuum onto itself after a time interval $T$,

$\displaystyle Z[T]=\langle 0|{\rm e}^{-HT}|0\rangle={\rm e}^{-E_{0}T}.$ (5)

On the other hand, path integrals correspond to transition amplitudes between
eigenstates of the field operators in the Heisenberg picture. Denoting the
different fields collectively as $\phi$, one has

$\displaystyle{\int_{\phi(-T/2)=\bar{\phi}}^{\phi(T/2)=\bar{\phi}^{\prime}}}{\cal
D}\phi{\rm e}^{-{S_{E}[\phi;\,T]}}={}_{T/2}\langle\bar{\phi}^{\prime}|{\rm
e}^{-HT}|\bar{\phi}\rangle_{-T/2}=\sum_{n}{}_{T/2}\langle\bar{\phi}^{\prime}|n\rangle\langle
n|\bar{\phi}\rangle_{-T/2}{\rm e}^{-E_{n}T}.$ (6)

Above, the path integration is such that the field trajectories are assumed to
tend to $\bar{\phi}^{\prime},\bar{\phi}$ for $t\rightarrow-T/2,T/2$,
respectively. $|\bar{\phi}\rangle_{t}$ are eigenstates of the Heisenberg-
picture field operators $\hat{\phi}(t)$, with eigenvalue $\bar{\phi}$. Before
the last equality in Eq. (6) we have inserted the spectral decomposition of
the identity into eigenstates of the Hamiltonian,
$H|n\rangle=E_{n}|n\rangle,\sum_{n}|n\rangle\langle n|=\mathbb{I}$. Clearly,
the path integration in Eq. (6) is a very different object than the vacuum
partition function in Eq. (5). To recover the latter, one can consider two
options, as detailed next.

First, one might take $T\rightarrow\infty$, such that for finite
$\langle\bar{\phi}^{\prime}|n\rangle,\langle n|\bar{\phi}\rangle$ the
contributions of the excited states in Eq. (6) become infinitely suppressed
with respect to the vacuum contribution. With $T\rightarrow\infty$, the
spacetime volume $VT$ goes to infinity, and the boundary conditions remain
arbitrary. It is generally expected from Picard-Lefschetz theory that the full
path integral with open boundary conditions can be expressed as a sum over
path integrations over complexified steepest descent field paths, each of them
passing through a saddle point of the Euclidean action [10]. Over each path,
$\exp(-S_{E})$ is maximal at the saddle, so that if the saddles have $S_{E}>0$
one only needs to consider those with finite action. For
$VT\rightarrow\infty$, finite action requires the fields to approach pure
gauge configurations asymptotically. In this case the condition in Eq. (3)
applies, and the saddle points fall into equivalence classes with integer
topological charge $\Delta n$. As the steepest descent field paths can be
described by means of a continuous deformation from the saddle points, we
expect each path to fall into a single equivalence class. Hence for
$VT\rightarrow\infty$ one can write

$\displaystyle Z=\sum_{\Delta n=-\infty}^{\infty}Z_{\Delta n}\text{ \,\,\, for
infinite spacetime }(T\rightarrow\infty).$ (7)

As the boundary conditions for infinite $T$ are arbitrary, there is nothing
wrong a priori with choosing fixed boundary conditions. In particular, one
could restrict the path integration to a single topological sector $\Delta m$
at infinite volume. As $\Delta m$ is gauge invariant, so is $Z_{\Delta m}$. In
this case, $\theta$ would always enter all correlators through a global phase
$\exp({\rm i}\Delta m\theta)$, which cancels out when normalizing the
correlators by the partition function. Hence there can’t be any $CP$-violation
in a fixed topological sector. As will be seen, this conclusion still holds
for $Z$ as in Eq. (7).

The alternative to an infinite spacetime volume is significantly more
challenging. In order to ensure that one projects into the vacuum state, one
should integrate over the boundary conditions $\bar{\phi},\bar{\phi}^{\prime}$
in Eq. (6) with weights given by the vacuum wave functionals $\langle
0|\bar{\phi}^{\prime}\rangle$, $\langle 0|\bar{\phi}\rangle$. Indeed, by
appropriately inserting spectral resolutions of the identity in terms of field
eigenstates, one has [11]

$\displaystyle\begin{aligned} \langle 0|{\rm e}^{-HT}|0\rangle=&\,\int{\cal
D}\bar{\phi}{\cal D}\bar{\phi}^{\prime}\langle
0|\bar{\phi}^{\prime}\rangle_{T/2}\,\,{}_{T/2}\langle\bar{\phi}^{\prime}|{\rm
e}^{-HT}|\bar{\phi}\rangle_{-T/2}\,\,{}_{-T/2}\langle\bar{\phi}|0\rangle\\\
=&\,\int{\cal D}\bar{\phi}{\cal D}\bar{\phi}^{\prime}\langle
0|\bar{\phi}^{\prime}\rangle_{T/2}\,\,{}_{-T/2}\langle\bar{\phi}|0\rangle{\int_{\phi(-T/2)=\bar{\phi}}^{\phi(T/2)=\bar{\phi}^{\prime}}}{\cal
D}\phi{\rm e}^{-{{S_{E}[\phi;\,T]}}}.\end{aligned}$ (8)

The computation of all path integrations in Eq. (8) is a daunting task, not
least because the vacuum wave functionals are unknown. Moreover, with the
boundary conditions $\bar{\phi},\bar{\phi}^{\prime}$ applied at a finite
surface, there is no reason that $\bar{\phi},\bar{\phi}^{\prime}$ correspond
to pure gauge configurations, and as such the classification into equivalence
classes with integer topological classes does not apply. Still, such an
assumption is usually made in the literature. In this spirit, starting without
a $\theta$ term in the action, the said term and the formal sum of Eq. (7)
could be recovered for finite $\Omega$ by assuming that the vacuum wave-
functionals are a linear combination of $\delta$ functionals with support on
static pure gauge configurations (i.e. classical vacua). The latter can be
classified with integer Chern-Simons numbers $k$ (so that in an appropriate
gauge $\Delta n=k(T/2)-k(-T/2)$) and when demanding gauge invariance the
ensuing wave-functional corresponds to the usual $\theta$ vacuum in the
literature, $|0\rangle=|\theta\rangle=\sum_{k}\exp({\rm i}k\theta)|k\rangle$.
Such a wave functional with no support outside classical vacua has no clear
justification, and certainly goes against the usual behaviour in quantum
mechanics.

In conclusion, the only practical way to ensure that one is computing the
vacuum partition function relies on taking the infinite spacetime volume
limit, more explicitly taking $T\rightarrow\infty$, regardless of whether the
spatial dimensions are infinite or compact without boundaries. Only after
taking this limit one can assume that the field configurations fall into
classes with integer $\Delta n$.

##  3 Correlation functions from the index theorem

The aim of this section is to show how the functional dependence of the
partition functions $Z_{\Delta n}$ on $\theta$, the spacetime volume $\Omega$
and the complex fermion masses $\mathfrak{m}_{i},\mathfrak{m}_{i}^{*}$ is
constrained. Knowing the dependence of $Z_{\Delta n}$ on
$\mathfrak{m}_{i},\mathfrak{m}_{i}^{*}$ allows one to obtain spacetime-
averaged fermion correlators—which should be sensitive to constant $CP$-odd
phases coming from $\theta,\alpha_{i}$—from the partition function with open
boundary conditions $Z$. This is because $\mathfrak{m}_{i}$ and
$\mathfrak{m}_{i}^{*}$ can be seen as sources for integrated two-point
functions. As the topological charge $\Delta n$ is only enforced to be an
integer for an infinite $\Omega$, we will take the limit
$\Omega\rightarrow\infty$ before summing over topological sectors. Within a
fixed topological sector $\Delta n$, the spacetime averages of the fermionic
correlators can be obtained as

$\displaystyle\begin{aligned} \frac{1}{\Omega}\int{\rm
d}^{4}x\,\langle\bar{\psi}_{i}P_{\rm R}\psi_{i}\rangle_{\Delta
n}=&\,-\frac{1}{\Omega}\frac{\partial}{\partial\mathfrak{m}_{i}}Z_{\Delta
n},\\\ \frac{1}{\Omega}\int{\rm d}^{4}x\,\langle\bar{\psi}_{i}P_{\rm
L}\psi_{i}\rangle_{\Delta
n}=&\,-\frac{1}{\Omega}\frac{\partial}{\partial\mathfrak{m}^{*}_{i}}Z_{\Delta
n},\end{aligned}$ (9)

as should be clear from Eq. (1) and the fact that the path integrals are
weighted by $\exp(-\int d^{4}x{\cal L})$. Let’s consider now a partition of
$\Omega$ into subvolumes $\Omega_{1}$ and $\Omega_{2}$. Following Ref. [4],
noting that $\Delta n$ is a surface integral allows to write the partition
function as

$\displaystyle Z_{\Delta n}(\Omega)=\sum_{\Delta
n_{1}=-\infty}^{\infty}Z_{{\Delta n}_{1}}(\Omega_{1})Z_{{\Delta n}-{\Delta
n}_{1}}(\Omega_{2}).$ (10)

Next, we aim to express $Z_{\Delta n}$ as a product of phases times a real
function. First, there is the phase $Z_{\Delta n}(\Omega)\propto{\rm e}^{{\rm
i}\Delta n\theta}$ that follows from $Z_{\Delta n}$ being a path integration
over $\exp(-{\cal L})$ and from Eqs. (1), (3). Other complex phases can only
come from the $\alpha_{i}$ that enter through the fermionic path integration.
At one-loop order about each saddle-point, the fermionic path integration
leads to a determinant of the massive Dirac operator
$\not{D}+\mathfrak{m}_{i}P_{R}+\mathfrak{m}_{i}^{*}P_{L}$ at the saddle for
each flavour. The phase of the determinant is fixed by parity and the Atiyah-
Singer index theorem [6]. Indeed, parity relates eigenvalues of eigenfunctions
not annihilated by $\not{D}$, which come in mutually conjugate pairs and hence
do not contribute to the phase [7]. The latter is then fully determined by the
zero modes of $\not{D}$. Given the structure of the operator
$\not{D}+\mathfrak{m}_{i}P_{R}+\mathfrak{m}_{i}^{*}P_{L}$, right/left-handed
zero modes have phases ${\rm e}^{\pm{\rm i}\alpha_{i}\theta}$, respectively,
and the total phase per fermion flavour $i$ is equal to $\alpha_{i}$ times the
difference in the number of right and left-handed zero modes. But according to
the index theorem, this quantity coincides with $\Delta n$. Accounting for all
flavours, one gets a phase ${\rm e}^{{\rm i}\Delta n\bar{\alpha}}$. Together
with the $\theta$-dependent phase and using Eq. (1), this gives

$\displaystyle Z_{\Delta n}(\Omega)={\rm e}^{{\rm i}\Delta
n\bar{\theta}}\tilde{g}_{\Delta n}(\Omega)\quad\text{ with
}\quad\tilde{g}_{\Delta n}(\Omega)\in\mathbb{R}.$ (11)

Equation (10) then yields the relations

$\displaystyle\tilde{g}_{{\Delta
n}}(\Omega_{1}+\Omega_{2})=\\!\\!\sum_{{\Delta
n}_{1}=-\infty}^{\infty}\tilde{g}_{{\Delta
n}_{1}}(\Omega_{1})\tilde{g}_{{\Delta n}-{\Delta n}_{1}}(\Omega_{2}).$ (12)

Setting $\Omega_{2}=0$ gives $\tilde{g}_{\Delta n}(0)=\delta_{\Delta n,0}$.
Under parity, $\Delta n$ changes sign, and since the real $\tilde{g}_{\Delta
n}(0)$ are insensitive to $CP$-odd phases, one must have $\tilde{g}_{-\Delta
n}(\Omega)=\tilde{g}_{\Delta n}(\Omega)$. This motivates the ansatz

$\displaystyle\tilde{g}_{\Delta n}(\Omega)=\Omega^{|{\Delta n}|}f_{|{\Delta
n}|}(\Omega^{2}),\quad f_{|{\Delta n}|}(0)\neq 0.$ (13)

Under the assumption of analyticity in $\Omega$, the solution turns out to be
unique, giving in the end [7]

$\displaystyle Z_{\Delta n}(\Omega)=I_{\Delta n}(2\beta\Omega)\,{\rm e}^{{\rm
i}\bar{\theta}\Delta n}\,,$ (14)

where $I_{\Delta n}$ are modified Bessel functions and $\beta$ is a real
parameter that can only depend on the moduli of the $\mathfrak{m}_{i}$:
$\beta=\beta(\mathfrak{m}_{j}\mathfrak{m}_{j}^{*})$. From Eqs. (9), (14),
normalizing by the full partition function and summing over $\Delta n$ after
taking the limit $\Omega\rightarrow\infty$, we obtain the following
correlators,

$\displaystyle\begin{aligned} \frac{1}{\Omega}\int
d^{4}x\,\langle\bar{\psi}_{i}P_{\rm
R}\psi_{i}\rangle\equiv\lim_{N\rightarrow\infty}\lim_{\Omega\rightarrow\infty}\frac{\sum_{|\Delta
n|<N}\int d^{4}x\langle\bar{\psi}_{i}P_{\rm R}\psi_{i}\rangle_{\Delta
n}}{\Omega\sum_{|\Delta m|<N}Z_{\Delta
m}}={-}2\mathfrak{m}^{*}_{i}\,\partial_{\mathfrak{m}_{i}\mathfrak{m}^{*}_{i}}\beta(\mathfrak{m}_{k}\mathfrak{m}^{*}_{k}),\\\
\frac{1}{\Omega}\int
d^{4}x\,\langle\bar{\psi}_{i}P_{L}\psi_{i}\rangle\equiv\lim_{N\rightarrow\infty}\lim_{\Omega\rightarrow\infty}\frac{\sum_{|\Delta
n|<N}\int d^{4}x\langle\bar{\psi}_{i}P_{\rm L}\psi_{i}\rangle_{\Delta
n}}{\Omega\sum_{|\Delta m|<N}Z_{\Delta
m}}={-}2\mathfrak{m}_{i}\,\partial_{\mathfrak{m}_{i}\mathfrak{m}^{*}_{i}}\beta(\mathfrak{m}_{k}\mathfrak{m}^{*}_{k}).\end{aligned}$
(15)

To arrive to the previous expressions, we used $\lim_{x\to\infty}I_{\Delta
n}(x)/I_{\Delta n^{\prime}}(x)=1$. The results remain $\theta$-independent and
the correlators have phases which, being aligned with those of the tree-level
masses, can be rotated away with suitable chiral field redefinitions. By
taking additional derivatives with respect to the masses
$\mathfrak{m}_{j}$,$\mathfrak{m}_{j}^{*}$, the results can be extended to
correlation functions involving more fermion fields. While the field
redefinitions that yield real correlators change $\theta$, this has no impact
on correlators/physical observables as they are $\theta$-independent. Hence,
there is no $CP$ violation. This conclusion depends crucially on the order of
the limits in Eqs. (15). Using the (mathematically incorrect) opposite
ordering, one obtains a dependence on both $\theta$ and the $\alpha_{i}$, so
that observables would be predicted to depend on $\bar{\theta}$, recovering
the traditional picture.

##  4 Discussion

Our results for spacetime-averaged fermion correlators in QCD imply that there
cannot be any $CP$-violation coming from nonperturbative corrections to
fermionic interactions. Both the form of the partition functions in Eq. (14)
and the alignment of correlators with the phases $\alpha_{i}$ can be recovered
as well with instanton calculus [7]. When matching our results for fermion
correlators with effective operators, in particular the ’t Hooft interactions
with $2N_{f}$ fermion fields that are expected from anomalies, one infers that
the phases of the ‘t Hooft vertices are fixed by the $\alpha_{i}$:

$\displaystyle{\cal L}_{\rm eff}\supset-\Gamma_{N_{f}}{\rm e}^{{\rm
i}\xi}\prod_{j=1}^{N_{f}}(\bar{\psi}_{j}P_{\rm L}\psi_{j})-\Gamma_{N_{f}}{\rm
e}^{-{\rm i}\xi}\prod_{j=1}^{N_{f}}(\bar{\psi}_{j}P_{\rm
R}\psi_{j})\,,\quad\xi=-\bar{\alpha}.$ (16)

The former result applies regardless of whether one sums over topological
sectors in the partition function or not. The only other QCD computation of
the phase in the ‘t Hooft vertices that we are aware of is that of ‘t Hooft
himself, who used instanton calculus in the dilute gas approximation, but
assumed the opposite order of limits than the one used in Eq. (15) and
concluded that the phase $\xi$ in Eq. (16) was $\xi=\theta$ [8, 9]. Again,
mathematical consistency demands the ordering of limits used in eq. (15),
leading to $\xi=-\bar{\alpha}$ (which complies with chiral symmetries) and no
$CP$ violation.

The effective interactions in Eq. (16) can be matched to the determinant terms
in the chiral Lagrangian, whose phases are then predicted to be fixed by
$\bar{\alpha}$. This leads to no $CP$-violation in the low-energy effective
theory, and in particular no dipole moment for the neutron, while still
leading to an enhancement of the $\eta^{\prime}$ mass [7].

## References

  * [1] C. G. Callan, Jr., R. F. Dashen and D. J. Gross, Phys. Lett. B 63 (1976), 334-340 doi:10.1016/0370-2693(76)90277-X
  * [2] R. Jackiw and C. Rebbi, Phys. Rev. Lett. 37 (1976), 172-175 doi:10.1103/PhysRevLett.37.172
  * [3] C. Abel et al. Phys. Rev. Lett. 124 (2020) no.8, 081803 doi:10.1103/PhysRevLett.124.081803 [arXiv:2001.11966 [hep-ex]].
  * [4] S. Weinberg, The quantum theory of fields. Vol. 2: Modern applications
  * [5] A. A. Belavin, A. M. Polyakov, A. S. Schwartz and Y. S. Tyupkin, Phys. Lett. B 59 (1975), 85-87 doi:10.1016/0370-2693(75)90163-X
  * [6] M. F. Atiyah and I. M. Singer, Bull. Am. Math. Soc. 69 (1969), 422-433 doi:10.1090/S0002-9904-1963-10957-X
  * [7] W. Y. Ai, J. S. Cruz, B. Garbrecht and C. Tamarit, Phys. Lett. B 822 (2021), 136616 doi:10.1016/j.physletb.2021.136616 [arXiv:2001.07152 [hep-th]].
  * [8] G. ’t Hooft, Phys. Rev. D 14 (1976), 3432-3450 [erratum: Phys. Rev. D 18 (1978), 2199] doi:10.1103/PhysRevD.14.3432
  * [9] G. ’t Hooft, Phys. Rept. 142 (1986), 357-387 doi:10.1016/0370-1573(86)90117-1
  * [10] E. Witten, AMS/IP Stud. Adv. Math. 50 (2011), 347-446 [arXiv:1001.2933 [hep-th]].
  * [11] A. D. Plascencia and C. Tamarit, JHEP 10 (2016), 099 doi:10.1007/JHEP10(2016)099 [arXiv:1510.07613 [hep-ph]].