# Anomaly Detection in Gravitational Waves data using Convolutional
AutoEncoders

Filip Morawski<EMAIL_ADDRESS>Nicolaus Copernicus Astronomical Center,
Polish Academy of Sciences, Bartycka 18, 00-716, Warsaw, Poland  Michał Bejger
Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka
18, 00-716, Warsaw, Poland  Elena Cuoco  European Gravitational Observatory
(EGO), I-56021 Cascina, Pisa, Italy  Scuola Normale Superiore, Piazza dei
Cavalieri 7, I-56126 Pisa, Italy  INFN, Sezione di Pisa, Largo Bruno
Pontecorvo, 3, I-56127 Pisa, Italy  Luigia Petre  Department of Computer
Science, Faculty of Science and Engineering  Åbo Akademi University,
Tuomiokirkontori 3, 20500 Turku, Finland

###### Abstract

As of this moment, fifty gravitational waves (GW) detections have been
announced, thanks to the observational efforts of the LIGO-Virgo
Collaboration, working with the Advanced LIGO and the Advanced Virgo
interferometers. The detection of signals is complicated by the noise-
dominated nature of the data. Conventional approaches in GW detection
procedures require either precise knowledge of the GW waveform in the context
of matched filtering searches or coincident analysis of data from multiple
detectors. Furthermore, the analysis is prone to contamination by instrumental
or environmental artifacts called glitches which either mimic astrophysical
signals or reduce the overall quality of data.

In this paper, we propose an alternative generic method of studying GW data
based on detecting anomalies. The anomalies we study are transient signals,
different from the slow non-stationary noise of the detector. Presented in the
manuscript anomalies are mostly based on the GW emitted by the mergers of
binary black hole systems. However, the presented study of anomalies is not
limited only to GW alone, but also includes glitches occurring in the real
LIGO/Virgo dataset available at the Gravitational Waves Open Science Center.
To search for anomalies we employ deep learning algorithms, namely
convolutional autoencoders, which are trained on both simulated and real
detector data. We demonstrate the capabilities of our deep learning
implementation in the reconstruction of injected signals. We study the
influence of the GW strength, defined in terms of matched filter Signal-to-
Noise Ratio, on the detection of anomalies. Moreover, we present the
application of our method for the localization in time of anomalies in the
studied time-series data. We validate the results of anomaly searches on real
data containing confirmed gravitational wave detections; we thus prove the
generalization capabilities of our method, towards detecting GWs unknown to
our deep learning models during training.

##  1 Introduction

The first gravitational wave (GW) detection on September 14, 2015 [1]
inaugurated a new era in astrophysics. The joint observational effort of the
LIGO and Virgo Collaborations, working with the Advanced LIGO [2] and the
Advanced Virgo [3] interferometers in a global network has provided, until the
suspension of the LIGO-Virgo observing run O3 in Spring 2020, fifty GW
candidate signal detections [4, 5, 6, 7]. Such an impressive number of
statistically-significant signal candidates allows for the verification of
many theoretical models, describing various sources of GW radiation, like
binary systems of black holes (BH) and neutron stars (NS), as well as the very
nature of gravity [8]. As a result of continuous work to improve their
sensitivity, the Advanced LIGO and the Advanced Virgo interferometers will
soon probe a much larger volume of space and expand the capability of
discovering new GW sources. According to theoretical predictions, tens of BH
mergers and a few NS mergers will be soon routinely registered every year [9].
Such a large number of events will deliver extraordinary information about the
nature of those objects, phenomena, as well as the space-time itself.

What can be observed through the GW window strongly depends on the fidelity of
the data analysis methods. The GW data is always noise dominated:
astrophysical signals are buried deep into the detectors’ noise, caused by
various sources. Some of them are associated with the environment [8], e.g.,
the seismic and environmental activity; others originate from the detector
itself [9], such as the thermal fluctuations of the mirrors and the laser beam
photon shot noise [10] \- time-dependent fluctuations in the laser interacting
with the mirrors of the interferometer. In order to reveal the hidden GWs, the
traditional approach is based on matched filtering algorithms [11]. The
examples of existing pipelines utilising a matched filtering approach are
PyCBC [12] and GstLAL [13]. Assuming the model gravitational waveform is
known, which is not always the case for all astrophysical sources, the
algorithm scans the data to match an optimal template using a template bank.
Matched filtering is computationally expensive and in case of large template
banks requires a lot of computational resources. As an optimal method of
filtering in stationary Gaussian noise, it is prone to contamination by non-
stationary instrumental artifacts, the so-called glitches, which either mimic
GW signals or reduce the quality of the collected data. Existing alternatives
to matched filtering methods are based on unmodeled searches of GW burst
signals. Such GWs are of short duration with an unknown or partially modelled
waveform morphology related to the complicated or unknown astrophysics. These
signals are searched by measuring an excessive power in the time-frequency
domain that occurs coherently between multiple detectors. An example of GW
burst pipeline is coherent Wave Burst (cWB) [14, 15]. The sensitivity of this
method is affected by short duration glitches that may occur coincident in
time between multiple detectors. Versatile and rapid pipelines are needed to
deal with non-stationary noise (and instrumental glitches) and to perform
preliminary analysis of large amounts of data from multiple detectors.

Deep learning (DL) [16] fits in that role perfectly. DL has commenced a new
era of machine learning (ML), a field of computer science based on specially-
designed algorithms that can generalize (‘learn’) from examples in order to
solve problems and make predictions, without the need of being explicitly
programmed [17]. DL algorithms, based on the concept of neural networks -
models of neuron connectivity in the brain - are able to analyse different
representations of data with varied dimensionality, like images (spectrograms,
Q-transforms, Wavelet-transform) or time series. Moreover, they can quickly
process large amounts of data - a requirement for the real-time (low latency)
analysis at the GW interferometers.

The purpose of this work is to propose a generic, model-independent data
analysis method that searches for anomalies in signals recorded by the GW
detectors, based on performant DL models. The innovation of the proposed
method is that our DL model ‘learns’ the features of the noise and detects if
anomalies occur in the detected signal. We define an anomaly as an
extraordinary, sparse transient data feature, outstanding with respect to the
‘normal’ background noise of the detector. The anomaly signal could therefore
either be represented as a GW or an instrumental glitch. In particular, the
GWs studied in the presented work are based on the signals emitted by the
binary BH systems (BBH).

The term ‘anomaly’ is a well defined concept in statistics and data analysis
(hence also in machine learning), sometimes also referred to as ‘outlier’. The
DL methods we apply are fitting perfectly in their detection and analysis role
here, being especially suited to detect, compress, and reconstruct non-
linearities in the input signal data. The analysis assessing the capabilities
of the method is performed with the BBH signals, injected (added) to both
simulated and real detector data. Once trained and tested on injections, the
DL models are validated on the confirmed GW detections - real astrophysical
signals registered by the detectors. Our results provide a proof-of-concept
for the advantages of using DL methods in the GW data analysis, namely the
processing speed and ability to capture complicated non-linear relationships
in the data. These features enable our method to be potentially used as an
Event-Trigger-Generator (ETG). In this context, the advantage of our method is
the computational speed - ML algorithms once trained are extremely efficient
in the processing of data, in particular when used on computers equipped with
a Graphics Processing Units (GPU). By searching for anomalies in the
gravitational waves’ data via ML algorithms, our method could support the
currently existing ETG such as Omicron [18] and Q-transform based Omega [19].

In the GW astronomy, while DL is being actively researched, it is still quite
a novel method. Therefore this research fits very well in the early adoption
scheme of the modern state-of-the-art development in the field, a knowledge of
which will become indispensable in the nearest future. In the following we
mention a few interesting test cases. George et al. [20] developed the Deep
Filtering method for signal processing, based on a system of two deep
convolutional neural networks, designed to detect and estimate parameters of
compact binary coalescence signal in highly noisy time-series data streams.
The same authors have been involved in a group working on denoising
gravitational waves with autoencoders [21]. Dreissigacker et al. [22] have
been using DL as a search method for the continuous GW emitted by spinning
neutron stars. Similar work has been conducted by Morawski et al. [23] on the
application of convolutional neural networks for classifying continuous GW
signal candidates. Furthermore, DL has been used by Beheshtipour and Papa [24]
for clustering continuous gravitational wave candidates. DL has also been
successfully used in the classification of glitches by Razzano and Cuoco [25].
Finally DL has been used in searches for GW emitted by core-collapse supernova
explosions by Iess et al. [26].

Within the GW astrophysics, application of DL in anomaly detection is a fairly
new concept (see however a recent application in [27]). In different fields of
astrophysics the ML anomaly detection has already been successfully
implemented [28, 29] (see [30] for a recent general overview of ML in
astronomy). Outside astronomy, anomaly detection has also been proposed in
search for signatures of new physics in the Large Hadron Collider data [31],
where it was demonstrated to discover a specific class of highly-energetic
particle jets, without the prior knowledge of their specific features.

The outline of this work is as follows. In Section 2 we briefly discuss the DL
architecture applied in our work and in Section 3 we describe our data
generation procedures. In Section 4 we explain the results of the anomaly
detection studies performed on both simulated and real data. Finally, in
Section 5 we further discuss our results and draw some conclusions.

##  2 Deep Learning algorithms

In this paper we employ a combination of two deep learning methods (CNNs and
AEs) for distinguishing GW from noise signals. We first briefly describe these
learning methods and then we explain how we applied them to our problem.

###  2.1 CNN and AE

A _convolutional neural network_ [16] (CNN) is a deep, feed-forward artificial
neural network (processes the information one-way, from the input to the
output), designed for processing structured arrays of data, e.g., for
classifying images. The core feature of a CNN is the convolution operator,
that differentiates them from regular (linear) neural networks employing
simple matrix multiplication. The convolution operation envisions the input
structured array (say, of dimension $m\times n$) as a sequence of overlapping
elements of dimension, say, $p\times p$, often with $p<min(m,n)$. The
convolution operator inner-multiplies each such element with a so-called
kernel (or filter), of the same dimension ($p\times p$), that ‘slides’ over
each element of the original array. Each $p\times p$ element of the input
array is thus replaced by the result of the convolution operation (a number),
and we thus have dimensionality reduction. This somewhat simplified
convolution concept is illustrated in Figure 1.

Figure 1: The 2-dimensional 4$\times$5 array on the left is seen as 6
overlapping 3$\times$3 arrays; each of these smaller arrays is inner-
multiplied with the 3$\times$3 kernel, resulting in the 2$\times$3 array on
the right.

Multiple layers are typically sequentially used, each with varying numbers and
types of kernels. The kernels are chosen so that the network learns specific
features in one convolution (e.g., edges, corners, etc). Convolution layers
are alternated with other specialised layers, all further reducing the
dimensionality, until a final activation function maps the last layer to a
vector whose elements correspond to the desired options (classes) for
classifying the input structured array. Each class in this vector is usually a
probability, so that the class with the highest value is indicated as the
recognised class.

The overlapping of the input array elements together with the sliding kernel,
as well as the sequencing of applying convolution step-wise simulates the
structure and operation of the visual human cortex that processes incoming
images in a series of layers, by identifying progressively more complex
features. CNNs have been shown to work extremely well at picking up patterns
in the input structured arrays, for instance in images, and thus have rapidly
become the state-of-the-art in image classification and computer vision.

In our approach we use all the advantages brought by CNNs, but we further
embed them into an _autoencoder_ (AE) architecture [32]. An AE [16] is a
special type of deep artificial neural network that step-wise encodes and
compresses the input and then it (re)constructs an output based only on the
most compressed encoding, called hidden layer, latent representation, or
bottleneck. The main AE hypothesis is that some structure exists in the input
data, for instance some form of correlation between the input features; such
structure can be and is learnt by an AE and consequently leveraged when
forcing the input through the latent layer. (If the input features are
independent of each other, then the compression and subsequent reconstruction
are very difficult.)

An ideal AE is sensitive enough to the inputs to accurately build a
reconstruction and insensitive enough to the same inputs so that the model
does not memorize/overfit the training data. Differently formulated, this
forces the model to maintain in the latent representation only the variations
in data needed for reconstruction, without holding on to redundancies within
the input. There are different types of AE, but in our work we use the
undercomplete AE, where the dimension of the latent representation is strictly
smaller than the input dimension. In this way we take care of avoiding
overfitting, since the model will not be able to copy the input to the output.

Figure 2: An AE with one hidden layer.

The simple AE in Figure 2 (with only one hidden layer) works by first encoding
the input $y\in\mathbb{R}^{d}$ to the element $h\in\mathbb{R}^{p}$, where
$p<d$. Assuming we have the encoder activation function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$, then $h=f(Wy+b)$, where $f$ can
be sigmoid, ReLU, or something else, $W$ is a weight matrix and $b$ a bias
factor (initialised randomly and updated iteratively during training). The
decoder activation function $g:\mathbb{R}^{p}\rightarrow\mathbb{R}^{d}$ will
map $h$ to $\tilde{y}=g(\tilde{W}h+\tilde{b})$, where the activation function
$g$, the weight matrix $\tilde{W}$ and the bias factor $\tilde{b}$ may be
unrelated to $f$, $W$, and $b$. If only ReLU is used for activation and we
have only one hidden layer, then we have a linear AE; if we have more hidden
layers, or non-linear activation(s), then the AE becomes non-linear, which is
better at detecting abstract features. Thus, in general, the encoder and
decoder are proper neural networks in themselves, not simply activation
functions.

The training of the AE works by attempting to minimize the reconstruction
loss, most often using the mean squared error (MSE) formula [16]

${\cal
L}(y,\tilde{y})=\left\|{y-\tilde{y}}\right\|^{2}=\left\|{y-g(\tilde{W}f(Wy+b)+\tilde{b})}\right\|^{2}.$
(1)

The training of the network works by updating the parameters $W$, $b$,
$\tilde{W}$, $\tilde{b}$ until ${\cal L}(y,\tilde{y})$ is sufficiently small
and further training does not decrease it anymore - in that case the network
has converged. There are numerous algorithms for updating the parameters, but
here we use the ADAM algorithm (adaptive moment estimation) [33], as it adapts
the learning rate during training and has been empirically shown superior to
other methods for large datasets, large number of parameters, as well as non-
stationary input. The learning rate (training ‘step’ of updating) together
with the batch size (number of training examples that use the same learning
rate), the update method (also called as optimizer), the number and size of
the layers and the loss function are the hyperparameters of the AE.

Since the AE is reducing dimensionality for encoding, it has often been used
for only that - for instance, for feature learning. However, when compared to
other dimensionality reduction techniques, such as Principal Component
Analysis (PCA), AE appears as a powerful generalization, since it is able to
learn non-linear relationships in input data. While PCA attempts to discover a
lower dimensional hyperplane describing the original data, AE is capable of
learning non-linear manifolds of the least possible size, as illustrated in
Figure 3. Essentially, the AE learns a vector field for mapping input data
towards lower dimensional manifolds, that describe the high density region
where input data concentrates. If the manifold accurately describes input
data, then the AE has effectively learnt the input data.

Figure 3: Example of the manifold concept. Non linear vs linear dimension
reduction.

###  2.2 Applying the CNN AE

Although CNNs were designed for the analysis of 2D data (i.e., images) [16],
here we apply them for a simpler, 1D implementation, as our analysed signals
are time series. The CNN AE is designed to learn two things. In case of input
data instances containing only detector background noise (no anomalies), the
AE is trained to reconstruct the noise as closely as possible. However, in
case of instances containing an anomalous signal (GW or glitch in case of real
data), the AE is trained to disregard the anomaly and reconstruct the data as
if the signal was not present. By comparing the input to the output
reconstructed by the AE, the anomaly present in the instance of data time
series is recovered and further studied. Using the AE network is thus
instrumental, since it will reconstruct only the signals it has been trained
for, in our case the detector noise, and disregard anything else in its
reconstruction, in our case the GW or the glitch, as noise.

Figure 4: Diagram shows the networks’ layer structure and architecture. The
size below convolutional layers correspond to the Number of filters times
Kernels size of every layer. The second dimension of all layers is equal to
unity since the initial input data were 1D time-series vectors.

Depending on the presence of anomalies in the input data, the loss value
${\cal L}(y,\tilde{y})$ is expected to vary. The ${\cal L}(y,\tilde{y})$
computed for an ‘anomalous’ input reaches higher values than in case of
‘anomaly-free’ input, since the difference between $y$ (noise or noise with
anomaly) and $\tilde{y}$ (only noise) is larger. The difference between these
signals is proportional to the amplitude of the anomaly. As a result, the AE
trained on data containing stronger anomalies is expected to converge during
training towards higher ${\cal L}(y,\tilde{y})$.

The final architecture111We implement the algorithm in Python [34] using the
Keras/TensorFlow library [35, 36] with the GPU support. The development was
performed on the NVidia Quadro P6000 (sponsorship via the NVidia GPU seeding
grant). The production runs were deployed on the Prometheus cluster (Academic
Computer Centre CYFRONET AGH) equipped with Tesla K40 GPU nodes, running CUDA
10.0 [37] and the cuDNN 7.3.0 [38]. used in the Section 4 was chosen based on
empirical tests on the data. We tested architectures ranging from 1 to 8
hidden layers. The final layout of the architecture is presented in Figure 4.
The chosen architecture allowing to reach the minimum value for ${\cal
L}(y,\tilde{y})$ for the fixed training data set was the network containing 3
hidden convolutional layers: an encoding layer, a decoding layer, and a latent
representation layer in between them, with 256, 128, and 256 neurons,
respectively. The kernel size was fixed for all layers to 3$\times$1\. All but
the final output layer use the ReLU as the activation function, whereas the
final layer reconstructing the initial signal uses a sigmoid activation
function. The other hyperparameters used for training were the ADAM optimizer
[33] with learning rate of 0.0005 and batch size of 32.

In the following section, we detail our datasets used for training, validation
and testing.

##  3 Training data sets and data flow

We prepared two kinds of training datasets: a simplified one by means of
simulated detector strain time series (based on the colored normal
distribution of noise), denoted DataSet 1 (DS1), and a realistic one based on
the real LIGO-Virgo O2 observing run [39], publicly available at the
Gravitational Waves Open Science Center (GWOSC) [40] and denoted DataSet 2
(DS2). In both cases we use the same general data flow presented in Figure 5.

Figure 5: Data flow of the project. The first step consists of data
generation. GW waveforms were injected into the detector strain, either
simulated or real (e.g. public data from the GWOSC platform). Since the raw
time-series signal varied with frequency, the whitening procedure was applied
to simplify data for further training. The training of the DL algorithms aimed
to recover as many injected anomalies as possible, while limiting the false
positive rate of the pipeline (noise samples incorrectly classified as
anomaly).

The whitening mentioned in the workflow diagram removes the contribution of
the stationary detector noise and re-weights the sensitivity at different
frequencies [41]. As a result, the amplitude spectral density of the data
becomes uniform and GW signals buried in the data are easier to search for and
compare with each other. The whitening filter was re-computed separately for
DS1 and DS2 as well as for every interferometer to take into account the
differences in the sensitivity. The whitening procedure used for the analysis
of both studied datasets was conducted in the frequency domain following pyCBC
Python library modules [42].

To simulate an astrophysical GW signal emitted by a BBH we used the
IMRPhenomv4 waveform model [43], which includes the binary inspiral, merger of
the components and the final BH ringdown. The component BH masses $m_{1}$ and
$m_{2}$ of the waveform model were chosen to be compatible with the first
detected GW150914 [1]. We selected $m_{1}$, $m_{2}$ based on the Initial Mass
Function (IMF) in ranges associated with uncertainties of mass estimation for
GW150914: $m_{1}:32.5-40.3$ $M_{\odot}$, $m_{2}:26.2-33.6$ $M_{\odot}$. For
the index of the IMF power law, we chose the value of $\alpha=-2.35$ [44]. The
luminosity distances were chosen uniformly from the range between $200$ to
$800$ $Mpc$, to cover a realistic range of the matched filter Signal-to-Noise
Ratio (SNR) - from 4 to 40 varying for different interferometers as shown in
the bottom right plot in Figure 6. The position in the sky was chosen to be
optimal for every detector in a given moment of time. The examples of a few
simulated GW signals are presented on the top right plot in Figure 6.

The DS1 dataset was created for the two assumed sensitivity curves of the GW
detectors. Each curve described the level of the detector sensitivity with
respect to the frequency in such a way that the generated strain was mimicking
the realistic time series output. In our analysis we used the designed
sensitivity for the advanced Virgo (aVirgo) from O3 run (version without
squeezing) [45, 46] and the advanced LIGO (aLIGO) [47] interferometers.
‘Designed’ means that the interferometers were expected to reach this level of
sensitivity after all the planned upgrades. Band-pass filtering was then
applied to the generated noise to remove high frequency (above 1 kHz) and low
frequency (below 30 Hz, corresponding to the seismic noise) components from
the data, as current interferometers are not sensitive enough to detect GWs
outside that frequency range. The data was then resampled from 4096 Hz to 1024
Hz. An example of the output time series from the DS1 is shown on the bottom
left plot in Figure 6. Prepared in advance GW signals were injected into the
generated strain and subjected to the procedure of whitening.

Figure 6: Top left: Designed sensitivity of aVirgo and aLIGO interferometers.
The detectors were expected to reach this level of sensitivity over broad
range of frequencies after all planned upgrades. Top right: Examples of
generated BBH GW waveforms injected into the strain data as anomalies. The GW
were generated for the following parameters: blue signal - $m1=29$ $M\odot$,
$m2=24$ $M\odot$, $distance=440$ $Mpc$; orange signal - $m1=33$ $M\odot$,
$m2=27$ $M\odot$, $distance=380$ $Mpc$; green signal - $m1=28$ $M\odot$,
$m2=23$ $M\odot$, $distance=600$ $Mpc$. Bottom left: Examples of the simulated
data using the above sensitivity curves. Bottom right: Distributions of
matched filter SNR of simulated GW injected into the real data for detectors:
LIGO Hanford (blue), LIGO Livingston (orange) and Virgo (green) as well
simulated data for: aLIGO (violet) and aVirgo (pink).

The DS2 (realistic) dataset was created based on the publicly available LIGO-
Virgo O2 observing run, using the data stored at the GWOSC platform [40] for
three interferometers: LIGO Hanford, LIGO Livingston and Virgo described in
the text with the respective abbreviations $H1$, $L1$ and $V1$. For each
detector, we chose six hours of data to train the DL models. For $L1$, we took
segments between 1187270656 and 1187295232 (in GPS time units); for $H1$,
between 1174958080 and 1174982656; whereas for $V1$, between 1187672064 and
1187696640. We injected into the strain the same GW signals as for the
simulated data. As a result, we obtained (for the same set of injections)
three different distributions of matched filter SNR, since every detector had
a different sensitivity. The distributions are illustrated in the bottom right
plot in Figure 6. For the comparison, we included additional SNR distributions
of the simulated datasets. The obtained real data strain with injected
anomalies was then subjected to the procedure of whitening and resampled from
4096 Hz to 1024 Hz.

In total, we generated five datasets (two simulated for aVirgo and aLIGO, and
three based on real LIGO Livingston, LIGO Hanford and Virgo O2 data) which
were further split into one second segments and divided for the training,
validation and testing datasets (65%, 10% and 25% respectively). An additional
test set, containing confirmed GW detections, has been created by using one
hour of data around GPS time for each confirmed GW, whitened and resampled as
described above. As confirmed GWs, we chose three BBH detections with the
highest network SNR from O2 run: GW150914, GW170608 and GW170814 [48].

##  4 Results

The results presented below are split into subsections. The first presents the
results of the anomaly searches on the simulated dataset containing injected
GW. The second subsection covers searches of anomalies in the real data of
Virgo and LIGO interferometers. The last subsection presents capabilities of
anomaly searches on the real data containing confirmed detections of GW.

###  4.1 Anomaly searches on simulated data

The CNN-AE described in Section 2 was first trained on the whitened, simulated
data containing GWs. The convergence of the model was achieved after 100
epochs, during which the MSE loss function reached a value of around $6\cdot
10^{-5}$ for the aVirgo data and $10^{-4}$ for the aLIGO data. We extended the
training for another 100 epochs to investigate the onset of overfitting.
However, the overfitting didn’t appear and MSE fluctuated around the values
mentioned above. The learning history of the AE, trained on both simulated
datasets, with results for both the training and the validation sets is shown
on the left plot in Figure 7. For aVirgo, the same set of gravitational
waveforms covered a SNR range of smaller values than for aLIGO, as a result of
the worse detector sensitivity (see the top left plot in Figure 6 for
comparison). This, in turn, resulted in the convergence towards a lower MSE
for the aVirgo dataset, since the difference between the ‘anomalous’ input and
‘anomalous-free’ reconstruction were smaller than in the case of aLIGO data
(see Section 2.2 for more details).

Figure 7: The learning history of the AE trained on the whitened, simulated
aVirgo and aLIGO datasets (left plot) as well as on the whitened, real data
for $L1$, $H1$ and $V1$ datasets (right plot). The convergence was achieved
after 100 epochs. We trained for 200 epochs, to investigate the onset of
overfitting: it did not appear.

Furthermore, as described in Section 2.2, the correctly trained CNN-AE
reconstructed the pure detectors’ noise, regardless of the anomalies presence
in the input data. Then, by subtracting the input from the reconstructed
output we aimed to recover the underlying signal. We computed the mentioned
differences between the input initial data and the output AE reconstruction.
Examples of the results are presented in Figure 8, where blue, dashed lines in
the background correspond to the expected results, whereas red lines
correspond to values obtained with the AE. In the majority of cases, we were
able to correctly reconstruct the GW waveforms for aLIGO data. The anomaly was
significantly different from the surrounding noise, with its recovered part
mostly associated with the merger part of the gravitational waveforms,
together with a small contribution from the inspiral part. In contrast, for
the aVirgo dataset, the reconstruction was significantly worse, often
dominated by the surrounding noise. In rare cases, as presented on the right
plot in Figure 8, the merger part was reconstructed. In the Appendix A we
present the summary of the match between the injected and reconstructed
waveforms using $<x_{1}|x_{2}>$ metric performed in the time domain.

Figure 8: The examples of signal reconstruction using the CNN-AE on aLIGO
(left plot) and aVirgo data (right plot). The plots (red curves) were
generated after subtracting the input signal from the AE predictions. For
comparison, we added the expected signal using blue, dashed curve (difference
between the input signal and the ground-truth output signal).

Since the AE was able to correctly reconstruct the majority of anomalies
buried in the data (at least for the aLIGO dataset), we aimed then to define a
metric allowing automatic anomaly detection. We chose the MSE as such a
metric, computed between the input data and the AE output as detailed in
Section 2.2. Figure 9 show histograms of MSE for the two different signal
types present in the studied data: noise and injected GW. As expected, values
of MSE for the noise were much smaller than for the latter case and close to
zero. There was a range of MSE in which histograms of both types overlapped
(denoted with the burgundy colour in Figure 9). Nevertheless, the majority of
noise with injected GW instances in the aLIGO dataset and almost half of these
instances in the aVirgo dataset had values larger than the noise. Moreover, we
added the detection threshold (defined in the following paragraph) to the
histograms to stress out how many of initially injected GW were correctly
detected as anomalies (hatched area in Figure 9).

We defined the threshold for the anomaly detection using the relation between
the false positive rate (FPR) and MSE. Then, by fixing FPR at a particular
value, we set the detection threshold ($DT$) on the corresponding MSE. In the
presented analysis, we fixed FPR at $5\%$, resulting in the following
thresholds: $DT_{simV}=1.6\cdot 10^{-5}$ for aVirgo and $DT_{simL}=3.1\cdot
10^{-5}$ for aLIGO. The results of the anomaly searches at FPR=5% are shown in
Table 1 in the form of confusion matrix. Additionally, the comparison in the
anomaly detection efficiency for both interferometers is shown on the left
plot in the Figure 10 in the form of Receiver Operating Characteristic (ROC)
curves. Over all ranges of FPR, the aLIGO detector achieved a significantly
higher detection efficiency, or True Positive Rate (TPR).

Figure 9: Distribution of MSE between the AE predictions and the input strain
for two types of studied signals: noise (blue histogram) and injected GW (red
histogram). In the broad range of MSE for aLIGO dataset the injected GW had
significantly higher MSE than noise, allowing a definite distinction between
these signal types. Additional black vertical line representing anomaly $DT$
was added to emphasize amount of detected anomalies (hatched areas).

Presented in Tab. 1 values quantitatively represent the results from both
panels in Fig. 9. Anomaly row relates to the hatched area from those panels.
In the case of aLIGO dataset 96 % of detected anomalies are correctly related
to the injected GW with minor contamination of noise samples. In case of
aVirgo 59 % of samples characterized by low SNR (around 10 and less) samples
did not exceed the $DT$ and contributed to the non-anomalous group. Details
presenting the relation between the SNR of the injected GW and MSE of the
reconstructed waveform can be found in the Appendix B.

| aLIGO | aVirgo
---|---|---
Injected GW | Noise | Injected GW | Noise
Anomaly | 96 % | 5 % | 41 % | 5 %
Non-anomaly | 4 % | 95 % | 59 % | 95 %
Table 1: Results of anomaly detection of CNN-AE at FPR=5% for aLIGO and aVirgo
dataset in the form of confusion matrix. Columns relate to the ground-truth
values whereas rows to the predictions. For aLIGO dataset significant majority
of detected anomalies corresponded to the data samples with injected GW.
However in case of aVirgo more than a half of data samples with injected GW
did not exceed the $DT_{simV}$ as a result of low SNR (see Fig. 6 for
comparison between aLIGO and aVirgo GW SNR distributions).

Figure 10: Receiver operating characteristic curves for simulated data (left
plot) and real data (right plot). Black vertical line corresponds to the
$FPR=5\%$ chosen as the criterion for the anomaly detection threshold.

###  4.2 Anomaly searches on real data

Later on, the AEs were trained on whitened, real data from the O2
observational run, collected for three existing interferometers: $V1$, $L1$
and $H1$ with injected BBH gravitational waveforms. The right plot in Figure 7
presents the learning history of the AE trained on the dataset for every
detector. As in the case of the simulated data, AE reached convergence after
around 100 epochs. We again prolonged the training to investigate the onset of
overfitting, which didn’t appear. Since the difference between the ‘anomalous’
input and the ‘anomalous-free’ reconstruction for $V1$ was the smallest among
the considered datasets (as a result of the SNR smallest range), the AE during
training converged towards the lowest MSE. Adequately, a smaller range of SNRs
for the $H1$ dataset with respect to the $L1$ dataset resulted in the
difference between the respective values of MSE (for $L1$ being higher, and
for $H1$ being lower).

The ability of the AE to reconstruct the detectors’ noise was investigated as
previously in the case of simulated data. The differences between the input
strain and AE reconstructions were compared with the expected values. Examples
of this comparison are presented in Figure 11. The summary of the match
between the injected and the reconstructed waveforms is presented in Appendix
A. The AE trained on $L1$ data achieved the best results, since the GW
injected into the strain were extracted to the greatest extent in the merger
part and partially in the inspiral part. On the other hand, the AE trained on
the other datasets reconstructed mainly the merger part. Overall, the AE
seemed to fail to reconstruct lower amplitudes and frequencies of the GW
signal (the inspiral and ringdown part).

Figure 11: Three examples of GW reconstruction using AE. Presented plots were
generated after subtracting from the AE predictions, input signal (red curve).
For the comparison we added the expected signal using blue, dashed curve
(difference between the input signal and the ground-truth output signal). From
left to right: results for $L1$, $H1$ and $V1$ detectors.

To compute the threshold for the anomaly detection, we generated histograms of
the MSE for every detectors’ dataset and compared its values with the FPR. The
results are shown in Figure 12. The anomalies covered a greater range of MSE
values than the noise, with an overlapping range varying for different
datasets (coloured in burgundy in Figure 12). In the case of $L1$ data, this
region was the smallest whereas for $V1$, it was the largest. As in the case
of simulated data, we defined a detection threshold for anomalies assuming
FPR=5%, which resulted in the following thresholds: $DT_{L1}=1.3\cdot 10^{-5}$
for $L1$, $DT_{H1}=2.2\cdot 10^{-5}$ for $H1$ and $DT_{V1}=4.3\cdot 10^{-6}$
for $V1$. The results of the anomaly searches on real datasets are present in
Table 2 in the form of confusion matrix. Anomaly row relates to the hatched
area on panels in Figure 12. In the case of $L1$ and $H1$ datasets, around
half of detected anomalies are correctly related to the injected GW. Samples
that did not exceed corresponding $DT$ had low SNR. That was also the case of
the $V1$ dataset where only 27 % of samples exceeded $DT_{V1}$. Details
presenting the relation between the SNR of the injected GW and MSE of the
reconstructed waveform for the real datasets can be found in the Appendix B.

Figure 12: Distribution of MSE between the AE predictions and the input strain for two types of studied signals: injected GW (red histograms) and noise (blue histograms). Plots correspond to the $L1$ (left plot), $H1$ (middle plot) and $V1$ (right plot) datasets. Additional black vertical line representing anomaly $DT$ was added to emphasize amount of detected anomalies (hatched areas).  | $L1$ | $H1$ | $V1$
---|---|---|---
Inj. GW | Noise | Inj. GW | Noise | Inj. GW | Noise
Anomaly | 52 % | 5 % | 50 % | 5 % | 27 % | 5 %
Non-anomaly | 48 % | 95 % | 50 % | 95 % | 73 % | 95 %
Table 2: Results of anomaly detection of CNN-AE at FPR=5% for $L1$, $H1$ and
$V1$ datasets in the form of confusion matrix. Columns relate to the ground-
truth values whereas rows to the predictions. For all datasets significant
majority of detected anomalies correctly corresponded to the data instances
with injected GW. However, more than a third part of non-anomalous class
(samples that did not exceed DT for a given detector) related to the low SNR
injected GW (see Fig. 6 for comparison between $L1$, $H1$ and $V1$ GW SNR
distributions).

Figure 13: Anomaly localisation based on the difference in peak positions
between reconstructed and injected signals. Left plot: The results for the
simulated datasets. Right plot: The results for the real datasets.

Additionally, we performed tests on the localisation in time of anomalies
detected by our method. We compared the known times of the GW injection into
the data with the time of the reconstructed signal. Specifically, we
subtracted the times corresponding to the maximum amplitude peaks of both
signals, and plotted histograms of the resulting differences for all the
anomalies that exceeded computed previously DT. For comparison, we performed
this procedure also for the simulated datasets. The results are shown in
Figure 13. In all the studied cases, around $95\%$ of detected anomalies were
localised within $0.05$ s intervals around the injection times. We conclude
that this feature may be useful not only for detection and reconstruction, but
also potentially for other applications, such as the localisation of signals
in the data from many detectors, and subsequent sky localisation of the
sources.

###  4.3 Anomaly searches on confirmed GW detections

Using a selection of real GW detections provided by the LIGO-Virgo
collaboration on the GWOSC platform [40], we tested the AE on three relatively
strong signals from the GWTC-1 O1-O2 catalog [39]: GW150914 [1], GW170608 [49]
and GW170814 [50]. The reported network SNR (square root of a sum of squares
of SNRs from individual detectors $\rho_{i}$)
$\rho_{net}=\sqrt{\sum_{i}\rho^{2}_{i}}$ is ${\simeq}24$ for GW150914,
${\simeq}15$ for GW170608 and ${\simeq}16$ for GW170814 [40]. Assuming two
equally sensitive detectors, each of them would measure SNR of ${\simeq}17$
for GW150914, ${\simeq}10$ for GW170608 and and ${\simeq}11$ for GW170814;
note that the GW170814 was a three-detector event, with the single-detector
SNRs in $H1$, $L1$ and $V1$ equal to 7.3, 13.7, and 4.4, respectively [50]. In
reality, due to differences in sensitivity the detectors registered the
signals with different SNRs, which were nevertheless near the single-detector
SNR detection threshold, established by the FPR=5% condition.

Figure 14: Test of the CNN-AE on the dataset containing confirmed detections
using $L1$ and $H1$ datasets for GW150914 (four upper plots), GW170608 (four
middle plots) and GW170814 (four bottom plots). Left plots: spectrograms
presenting the relation between frequency and time for the segment of real
data containing GW. Right plots: the difference between the CNN-AE predictions
and the input data.

After whitening, test data were fed into the AE. Then the reconstructed values
were subtracted from the input data. The results of subtraction for GW150914
are shown on the two upper, right plots in Figure 14 for the both LIGO
detectors. The presented signals resulted in the $MSE_{L1}=3.2\cdot 10^{-4}$
and $MSE_{H1}=1.0\cdot 10^{-3}$. For both detectors, MSE had significantly
higher values than the corresponding detection thresholds of FPR=5%.

In case of GW170608, AE detected the event; however, the reconstructed signal
was limited to the merger part for both the LIGO detectors as shown on the two
middle, right plots in Figure 14. A potential explanation for this weaker
reconstruction of this particular GW is related to the different mass ranges
of GW170608 and the GWs used for the training of the AE. GW170608 had BBH
component masses corresponding to around $m_{1}=11.0$ and $m_{2}=7.6$
$M_{\odot}$ which were substantially smaller than the masses used in the data
generation (see Section 3 for more details). We note that the H1 detector was
nominally outside of observing mode at the time of this event. This data
release includes H1 data around the time of this event, by using a modified
segment list, as was done for the published analysis [49, 48, 40].
Nevertheless, the AE detected GW170608, proving its generalisation
capabilities towards recognising gravitational waveforms it was not trained
for. This is a sure advantage of our proposed method, when compared to the
matched filtering method. Furthermore, the MSE values for both detectors had
slightly higher values above the detection thresholds at FPR=5%:
$MSE_{L1}=5.3\cdot 10^{-5}$ and $MSE_{H1}=4.1\cdot 10^{-5}$.

The last studied test case, GW170814, was detected in both LIGO detectors with
a substantial part of the waveform being recovered, as shown on the two
bottom, right plots in Figure 14. For the $V1$ data, AE failed to detect the
event evidently due to low reported SNR of ${\simeq}4$ [50]. The MSE for the
$L1$ and $H1$ datasets were above the data threshold, at $5\%$ and equal to
$MSE_{L1}=2.2\cdot 10^{-4}$, $MSE_{H1}=2.2\cdot 10^{-4}$, whereas for $V1$ the
MSE was below the threshold: $MSE_{V1}=1.8\cdot 10^{-6}$.

These results confirm that our proposed method of using CNN-AE for anomaly
searches is able to detect real GWs, even though the deep learning model was
trained on relatively uncomplicated datasets, based on the information related
to a particular GW waveform models, varying with respect to the limited range
of masses and distances.

##  5 Summary

In this paper, we proved that AEs are a potentially powerful method for
anomaly searches in the GW data. A relatively simple AE, consisting of only
three hidden layers, was capable of detecting anomalies, defined in terms of
transient BBH GWs (as well glitches in real data), and trained on either
simulated or real data. Moreover, our proposed method was able to detect as
anomalies all three confirmed GWs used as a test case and even partially
reconstruct the waveforms of the underlying signals.

In the proposed method we introduced a metric allowing for the automatic
detection of anomalies. The chosen metric was MSE. Using this metric, we
defined the threshold for the anomaly detection by associating MSE with FPR.
Presented in the manuscript results referred to FPR=5%. At such a threshold we
were able to detect in case of LIGO data, almost all of the injected GW into
the simulated dataset, as well as around 50% for real detector data. In the
case of Virgo, for simulated data around half of injected signals were
detected and quarter for real data. The reason for worse results on real data
taken from O1-O2 observational run, was substantial difference in the
sensitivity with respect to the simulated counterpart. The sensitivity of real
detectors were overall worse than the designed one. However, the real
sensitivity significantly improved for O3 observational run after detector
upgrades. Once the data from O3 becomes public, the anomaly searches of
proposed CNN-AE method are expected to improve.

Our method also proved to be useful in the precise time localisation of
anomalies. Anomalies exceeding the detection threshold for all the studied
datasets (real and simulated) turned out to be localised with the accuracy
below 0.05 second. Localisation of the peak amplitude may be particularly
useful for the localisation of GW sources on the sky in the multi-detector
analysis. However, such an application of our method requires a separate
study.

Finally, the successful detection of GW170608 proved the generalisation
capabilities of our AE, towards the detection of GWs with parameters that are
different than those of the GWs used for the AE training and in the data
nominally outside of the observing mode.

Among projects we are considering are applications of recurrent neural
networks instead of CNN in the anomaly searches as an alternative suited for
the time-series data as well as anomaly searches of different GW types such as
core collapse supernova signals.

The work was partially supported by the Polish NCN grants no.
2016/22/E/ST9/00037, 2017/26/M/ST9/00978 and 2020/37/N/ST9/02151, as well as
the European Cooperation in Science and Technology COST action G2Net
(CA17137). The Quadro P6000 used in this research was donated by the NVIDIA
Corporation. The production computations were supported in part by PL-Grid
Infrastructure and by the MNiSW grant for expansion of the strategic IT
infrastructure for science. This research has made use of data, software
and/or web tools obtained from the Gravitational Wave Open Science Center
(https://www.gw-openscience.org/ ), a service of LIGO Laboratory, the LIGO
Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and
Advanced LIGO are funded by the United States National Science Foundation
(NSF) as well as the Science and Technology Facilities Council (STFC) of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * [1] Abbott B P, Abbott R, Abbott T D, Abernathy M R et al. 2016 Phys. Rev. Lett. 116 061102 (Preprint 1602.03837)
  * [2] Aasi J, Abbott B P, Abbott R, Abbott T et al. 2015  32 074001 (Preprint 1411.4547)
  * [3] Acernese F, Agathos M, Agatsuma K, Aisa D et al. 2015  32 024001 (Preprint 1408.3978)
  * [4] O3 summary , see more: https://www.ligo.caltech.edu/WA/news/ligo20200326
  * [5] Virgo status , see more: https://www.virgo-gw.eu/status.html
  * [6] Gravitational-Wave Candidate Event Database, see more: https://gracedb.ligo.org/superevents/public/O3/
  * [7] Abbott R, Abbott T D, Abraham S and et al F A 2020 Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run (Preprint 2010.14527)
  * [8] Abbott B P, Abbott R, Abbott T D, Abraham S, Acernese F, Ackley K, Adams C, Adhikari R X, Adya V B, Affeldt C and et al 2019 The Astr. Jour. 882 L24 ISSN 2041-8213 URL http://dx.doi.org/10.3847/2041-8213/ab3800
  * [9] Acernese F, Agathos M, Agatsuma K, Aisa D et al. 2016 Phys. Rev. D 93 122003 (Preprint 1602.03839)
  * [10] Kumar A, Saminadayar L, Glattli D C, Jin Y and Etienne B 1996 Phys. Rev. Lett. 76(15) 2778–2781 URL https://link.aps.org/doi/10.1103/PhysRevLett.76.2778
  * [11] Owen B and Sathyaprakash B 1999 Phys. Rev. D 60(2) 022002 URL https://link.aps.org/doi/10.1103/PhysRevD.60.022002
  * [12] Usman S A, Nitz A H, Harry I W, Biwer C M, Brown D A, Cabero M, Capano C D, Dal Canton T, Dent T, Fairhurst S, Kehl M S, Keppel D, Krishnan B, Lenon A, Lundgren A, Nielsen A B, Pekowsky L P, Pfeiffer H P, Saulson P R, West M and Willis J L 2016 Classical and Quantum Gravity 33 215004 (Preprint 1508.02357)
  * [13] Sachdev S et al. 2019 (Preprint 1901.08580)
  * [14] Klimenko S et al. 2016 Phys. Rev. D 93 042004 (Preprint 1511.05999)
  * [15] Klimenko S et al. 2008 Classical and Quantum Gravity 25 114029 URL http://stacks.iop.org/0264-9381/25/i=11/a=114029
  * [16] Goodfellow I, Bengio Y and Courville A 2016 Deep Learning (The MIT Press) ISBN 0262035618, 9780262035613
  * [17] Samuel A L 1959 IBM Journal of Research and Development 3 210–229
  * [18] Robinet F, Arnaud N, Leroy N, Lundgren A, Macleod D and McIver J 2020 SoftwareX 12 100620 ISSN 2352-7110 URL https://www.sciencedirect.com/science/article/pii/S2352711020303332
  * [19] Chatterji S, Blackburn L, Martin G and Katsavounidis E 2004 Classical and Quantum Gravity 21 S1809–S1818 URL https://doi.org/10.1088/0264-9381/21/20/024
  * [20] George D and Huerta E 2018 Phys. Lett. B 778 64 – 70 ISSN 0370-2693 URL http://www.sciencedirect.com/science/article/pii/S0370269317310390
  * [21] Shen H, George D, Huerta E A and Zhao Z 2019 Denoising gravitational waves with enhanced deep recurrent denoising auto-encoders ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) pp 3237–3241
  * [22] Dreissigacker C, Sharma R, Messenger C, Zhao R and Prix R 2019 Phys. Rev. D 100(4) 044009 URL https://link.aps.org/doi/10.1103/PhysRevD.100.044009
  * [23] Morawski F, Bejger M and Ciecieląg P 2020 Machine Learning: Science and Technology URL http://iopscience.iop.org/10.1088/2632-2153/ab86c7
  * [24] Beheshtipour B and Papa M A 2020 Phys. Rev. D 101(6) 064009 URL https://link.aps.org/doi/10.1103/PhysRevD.101.064009
  * [25] Razzano M and Cuoco E 2018 Clas. Quant. Grav. 35 095016 URL https://doi.org/10.1088%2F1361-6382%2Faab793
  * [26] Iess A, Cuoco E, Morawski F and Powell J 2020 Machine Learning: Science and Technology 1 025014 URL https://doi.org/10.1088%2F2632-2153%2Fab7d31
  * [27] Corizzo R, Ceci M, Zdravevski E and Japkowicz N 2020 Expert Systems with Applications 151 113378 ISSN 0957-4174 URL http://www.sciencedirect.com/science/article/pii/S0957417420301986
  * [28] Giles D and Walkowicz L 2019 MNRAS 484 834–849 (Preprint 1812.07156)
  * [29] D’Addona M, Riccio G, Cavuoti S, Tortora C and Brescia M 2020 Anomaly detection in astrophysics: a comparison between unsupervised deep and machine learning on kids data (Preprint 2006.08235)
  * [30] Baron D 2019 arXiv e-prints arXiv:1904.07248 (Preprint 1904.07248)
  * [31] Farina M, Nakai Y and Shih D 2020 Phys. Rev. D 101(7) 075021 URL https://link.aps.org/doi/10.1103/PhysRevD.101.075021
  * [32] Baldi P 2011 Autoencoders, unsupervised learning and deep architectures Proceedings of the 2011 International Conference on Unsupervised and Transfer Learning Workshop - Volume 27 UTLW’11 (JMLR.org) p 37–50
  * [33] Kingma D P and Ba J 2014 arXiv e-prints arXiv:1412.6980 (Preprint 1412.6980)
  * [34] Van Rossum G and Drake F L 2009 Python 3 Reference Manual (Scotts Valley, CA: CreateSpace) ISBN 1441412697
  * [35] Chollet F et al. 2015 Keras https://keras.io
  * [36] Abadi M, Agarwal A, Barham P, Brevdo E et al. 2015 TensorFlow: Large-scale machine learning on heterogeneous systems software available from tensorflow.org URL https://www.tensorflow.org
  * [37] Nickolls J, Buck I, Garland M and Skadron K 2008 Queue 6 40–53 ISSN 1542-7730 URL http://doi.acm.org/10.1145/1365490.1365500
  * [38] Chetlur S, Woolley C, Vandermersch P, Cohen J, Tran J, Catanzaro B and Shelhamer E 2014 ArXiv e-prints (Preprint 1410.0759)
  * [39] Abbott B P, Abbott R, Abbott T D, Abraham S, Acernese F, Ackley K, Adams C et al. 2019 Physical Review X 9 031040 (Preprint 1811.12907)
  * [40] Collaboration T L S, the Virgo Collaboration, Abbott R, Abbott T D et al. 2019 (Preprint 1912.11716)
  * [41] Cuoco E, Calamai G, Fabbroni L, Losurdo G, Mazzoni M, Stanga R and Vetrano F 2001 Classical and Quantum Gravity 18 1727–1751 URL https://doi.org/10.1088%2F0264-9381%2F18%2F9%2F309
  * [42] Nitz A, Harry I, Brown D, Biwer C M, Willis J, Canton T D, Capano C, Pekowsky L, Dent T, Williamson A R, Davies G S, De S, Cabero M, Machenschalk B, Kumar P, Reyes S, Macleod D, dfinstad, Pannarale F, Massinger T, Tápai M, Singer L, Kumar S, Khan S, Fairhurst S, Nielsen A, SSingh087, shasvath, Gadre B U V and Dorrington I 2020 gwastro/pycbc: Pycbc release v1.16.10 URL https://doi.org/10.5281/zenodo.4063644
  * [43] Hannam M, Schmidt P, Bohé A, Haegel L, Husa S, Ohme F, Pratten G and Pürrer M 2014 Phys. Rev. Lett. 113(15) 151101 URL https://link.aps.org/doi/10.1103/PhysRevLett.113.151101
  * [44] Salpeter E E 1955 ApJ 121 161
  * [45] The Virgo Collaboration 2009 _Advanced Virgo Baseline Design, note VIR-027A-09, May 16 2009_ https://tds.virgo-gw.eu/?call_file=VIR-0027A-09.pdf
  * [46] Virgo Interferometer Monitoring webpage https://vim-online.virgo-gw.eu
  * [47] Updated Advanced LIGO sensitivity design curve https://dcc.ligo.org/LIGO-T1800044/public
  * [48] Abbott B P, Abbott R, Abbott T D and Abraham S e a (LIGO Scientific Collaboration and Virgo Collaboration) 2019 Phys. Rev. X 9(3) 031040 URL https://link.aps.org/doi/10.1103/PhysRevX.9.031040
  * [49] Abbott B P, Abbott R, Abbott T D, Acernese F, Ackley K, Adams C, Adams T, Addesso P, Adhikari R X, Adya V B, Affeldt C, Afrough M, Agarwal B, Agathos M, Agatsuma K, Aggarwal N, Aguiar O D, Aiello L, Ain A, Ajith P, Allen B, Allen G, Allocca A et al. 2017 ApJ 851 L35 (Preprint 1711.05578)
  * [50] Abbott B P, Abbott R, Abbott T D, Acernese F, Ackley K, Adams C, Adams T, Addesso P, Adhikari R X, Adya V B, Affeldt C, Afrough M, Agarwal B, Agathos M, Agatsuma K, Aggarwal N, Aguiar O D, Aiello L, Ain A, Ajith P, Allen B, Allen G, Allocca A, Altin P A, Amato A, Ananyeva A, Anderson S B, Anderson W G, Angelova S V, Antier S, Appert S, Arai K, Araya M C, Areeda J S, Arnaud N et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2017 Phys. Rev. Lett. 119(14) 141101 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.141101

##  Appendix A Match/overlap between the injected and reconstructed waveforms

###  A.1 Simulated dataset

To measure the match between the injected and the reconstructed waveforms we
used the normalized scalar product $<x_{1}|x_{2}>$ in the time domain
resulting in values in a range $(0,1)$. Zero related to no match, whereas one
to full match between the waveforms. Presented results in Fig. 15 corresponds
to the whole studied dataset for aLIGO and aVirgo detectors (left panel) as
well as samples exceeding the anomaly detection thresholds (right panel).
Applying respective $DT$ allowed to substantially reduce number of samples
with no match between the waveforms as a result of low SNR of injected GW.
Samples exceeding $DT$ were reconstructed to a greater extent which resulted
in $<x_{1}|x_{2}>$ closer to one.

Figure 15: Distribution of a match between the injected ($x_{1}$) and the
reconstructed ($x_{2}$) waveforms as a function of normalized scalar product
$<x_{1}|x_{2}>$. Scalar product equals to one relates to the full match
between waveforms whereas value of zero - no match. Left plot: results for the
whole studied dataset of particular detector; right plot: results of a match
after applying the anomaly detection threshold at FPR=5%.

###  A.2 Real dataset

The same metric as in case of simulated dataset was used to study the match
between the injected and the reconstructed waveforms for the real datasets.
Presented in Fig. 16 results show the similarities in the match between
consecutive datasets ($L1$, $H1$ and $V1$) as well as the effect of applying
the anomaly detection threshold. Samples with $<x_{1}|x_{2}>$ close to zero
had low SNR. As a result they were poorly reconstructed which in turn
translated into low value of MSE. Samples exceeding respective $DT$ were
reconstructed to a greater extent as in the case of simulated data. However,
overall match was worse - the mean values of $<x_{1}|x_{2}>$ for real datasets
were around 0.6, whereas for simulated data 0.8.

Figure 16: Distribution of a match between the injected ($x_{1}$) and the
reconstructed ($x_{2}$) waveforms as a function of normalized scalar product
$<x_{1}|x_{2}>$. Scalar product equals to one relates to the full match
between waveforms whereas value of 0 - no match. Left plot: results for the
whole studied dataset of particular detector; right plot: results of a match
after applying the anomaly detection threshold at FPR=5%.

##  Appendix B Signal-to-Noise Ratio vs Mean Squared Error

###  B.1 Simulated dataset

For the aLIGO dataset above SNR=20, the MSE-SNR relation was almost linear,
with a small spread of individual data instances along MSE. Whereas with the
decline of SNR, the spread of MSE significantly increased, characterized by
the non-linearity in the MSE-SNR relation. Anomalies around the same SNR for
values below $10$ varied up to an order of magnitude in MSE. Manual inspection
of data samples containing anomalies of low SNR provided an explanation of
this behaviour. In the analysed samples, only the merger part of the
gravitational waveform was recovered. For lower SNRs (below $10$) the recovery
was partial and dependent on the variability of the noise. If the amplitude of
the noise in a given data segment was small comparing to the injected GW
signal, the resulting MSE was higher than in the case of noise samples with
larger amplitudes. Overall, for the aLIGO dataset, the susceptibility of AE to
the local variability of the noise was inversely proportional to the SNR of
the injected anomaly.

In case of the aVirgo dataset, the mentioned susceptibility was more
significant. Matched filter SNRs for all injected GWs covered a smaller range
of values than for aLIGO (compare SNR ranges on the bottom right plot in
Figure 6 for aLIGO and aVirgo datasets). In the majority of studied cases, the
recovery of the anomaly was partial and limited to the merger part.

Figure 17: Relation between MSE and matched filter SNR for aVirgo and aLIGO
datasets colored with a distance to the GW source. Black horizontal line
corresponds to the detection threshold of anomalies at FPR=5%.

###  B.2 Real dataset

The relation between SNR and MSE presented similar features as for the
simulated data discussed above. The variability of MSE for samples of similar
SNR was largest among the weakest anomalies. The increase of the matched
filter SNR led to the decrease of the spread in MSE, and thus their relation
became more linear. This was the case for the AEs trained on the $L1$ and $H1$
datasets. In contrast, the AE trained on the $V1$ dataset was more susceptible
to the local variability of the noise. Since the injected anomalies had a
small SNR for $V1$ (majority of injected GW had SNR below $10$), their
amplitude was significantly lower than the detectors noise. As a result, the
anomaly detection depended on the variability of the noise, which had a random
character.

Figure 18: Relation between the MSE and the matched filter SNR for $L1$ (left
plot), $H1$ (middle plot) and $V1$ (right plot) datasets colored with a
distance to the GW source. Black horizontal line corresponds to the detection
threshold of anomalies at FPR=5%.


# Quasiparticle band structures of bulk and few-layer PdSe2 from first-
principles GW calculations

Han-gyu Kim  Hyoung Joon Choi<EMAIL_ADDRESS>Department of Physics,
Yonsei University, Seoul 03722, Korea

###### Abstract

We performed first-principles density functional theory (DFT) and GW
calculations to investigate electronic structures of bulk and few-layer PdSe2.
We obtained the quasiparticle band structure of bulk PdSe2, and the obtained
energy gap agrees excellently with the reported experimental value. For
monolayer and bilayer PdSe2, we obtained quasiparticle band structures with
respect to the vacuum level. We analyzed DFT and GW band structures in detail,
finding $k$-space positions of valence band maxima and conduction band minima,
effective masses, the quasiparticle density of states, work functions,
ionization potentials, electron affinities, and $k$-space shapes of electron
and hole pockets. These results provide a foundation for development of basic
studies and device applications.

##  I INTRODUCTION

Two-dimensional (2D) materials have attracted much attention due to their
versatile electronic and magnetic properties and fabrication of
heterojunctions by mechanical exfoliation and stacking. While graphene is
metallic without any band gap, transition-metal dichalcogenides (TMDs) are
semiconducting and applied to various device applications due to their finite
band gaps [1, 2, 3, 4, 5, 6, 7]. Recently, PdSe2, a layered TMD material with
stacks of a puckered pentagonal lattice, was adopted for device fabrication
and showed electron-dominant transport with high mobility, a large on-off
ratio of current [8], and air stability [9]. PdSe2 is also interesting for the
possibility of Stoner-type ferromagnetism in hole-doped 2D layers [10].

Bulk PdSe2 is a semiconductor with a band gap of 0.4 eV from the resistivity
measurement [11]. A previous density functional theory (DFT) calculation using
the local density approximation, however, showed that bulk PdSe2 is a
semimetal with no band gap [12], and a DFT calculation using the modified
Becke-Johnson potential showed that bulk PdSe2 is a semiconductor with a very
small band gap of 0.03 eV [13]. For accurate study of the band gap, the GW
method [14, 15, 16, 17, 18] is required. The GW method directly computes
quasiparticle (QP) energies, and it is very successful in describing the band
gap accurately [19]. For bulk PdSe2, so far, there is no report on band
structures using the GW method.

For monolayer and bilayer PdSe2, previous DFT calculations [9, 21, 20] showed
that their band gaps are about 1.3 and 0.8 eV, respectively. GW calculations
were also preformed previously, reporting that monolayer and bilayer PdSe2
have QP band gaps of 2.55 and 1.89 eV, respectively [21]. Meanwhile, band gaps
of monolayer and bilayer PdSe2 estimated from optical absorption experiments
are close to DFT results [9], so the difference between experimental band gaps
and GW band gaps is expected to be related to excitonic binding energies,
substrate effects, and structural defects such as Se vacancies [21, 22]. For
device applications, the accurate study of band-edge energies with respect to
the vacuum level is required, which is related to work functions, ionization
potentials, and electron affinities. These properties are not reported yet,
and can be obtained accurately by the GW method with careful convergence in 2D
geometries.

In the present work, we study electronic structures of bulk and few-layer
PdSe2 using DFT and GW calculations. We obtain quasiparticle band structures,
band gaps, the projected density of states, work functions, ionization
potentials, and electron affinities for the monolayer, bilayer, and bulk. We
also obtain $k$-space positions of valence band maxima (VBMs) and conduction
band minima (CBMs), effective masses at VBMs and CBMs, and $k$-space shapes of
electron and hole pockets.

This paper is organized as follows. In Sec. II, we describe our calculation
method including atomic structures of bulk and few-layer PdSe2 and parameters
for DFT and GW calculations. In Sec. III, we present and discuss electronic
structures of monolayer, bilayer, and bulk PdSe2 obtained from DFT and GW
calculations. Finally, we summarize our work in Sec. IV.

##  II METHODOLOGY

###  II.1 Atomic structures

We use the experimental atomic structure of bulk PdSe2 shown in Fig. 1.
Experimental lattice constants are $a$ = 5.7457, $b$ = 5.8678, and $c$ =
7.6946 Å [23], and in our work the lattice vectors $a,b,$ and $c$ are along
the $x,y,$ and $z$ axes, respectively. We construct atomic structures of
monolayer and bilayer PdSe2 by cutting the bulk atomic structure and including
a large enough vacuum region of about 20 Å. To simulate the vacuum level for
bulk PdSe2, we construct 12-layer PdSe2.

The monolayer PdSe2 is a triple layer consisting of three atomic layers of Se,
Pd, and Se, where each Pd atom is neighbored by four Se atoms and each Se atom
is neighbored by two Pd atoms and one Se atom, as shown in Figs. 1(b) and
1(c). The triple layer has a thickness of 1.45 Å and contains two Pd atoms and
four Se atoms in its unit cell, as shown in Fig. 1. The bulk geometry contains
two triple layers in the unit cell, as shown in Fig. 1(a), where triple layers
are separated by 2.40 Å. Atomic positions are from experimental values [23];
that is, one of the Pd atoms is at (0, 0, 0), and one of the Se atoms is at
(0.11125, 0.11799, 0.40573) in terms of fractional coordinates with respect to
lattice vectors, and other atomic positions are determined by symmetries.

The bulk geometry belongs to space group 61 ($Pbca$), which has inversion
symmetry and three glide-mirror planes consisting of $xy$, $yz$, and $zx$
mirror planes and translation by a half unit cell along the $x,y,$ and $z$
directions, respectively. The monolayer geometry belongs to space group 14
($P2{{}_{1}}/c$), which has inversion symmetry and one glide-mirror plane
consisting of the $yz$ mirror plane and translation by one half of the lattice
vector $b$. The bilayer geometry belongs to space group 7 ($Pc$), which has no
inversion symmetry but one glide-mirror plane consisting of the $yz$ mirror
plane and translation by one half of the lattice vector $b$.

Figure 1: Atomic structure of bulk PdSe2: (a) side view, (b) top view of the
first triple layer, and (c) top view of the second triple layer. Gray and
light green dots are Pd and Se atoms, respectively. In (a), distances between
layers are experimental ones. The lattice vectors $a,b,$ and $c$ are along the
$x,y,$ and $z$ axes, respectively. In (b) and (c), Se atoms above (below) the
Pd plane are represented with larger (smaller) light green dots.

###  II.2 DFT calculations

We perform DFT calculations [24, 25] using the QUANTUM ESPRESSO code [26] with
the Perdew-Burke-Ernzerhof-type generalized gradient approximation [27] to the
exchange-correlation energy. We use a kinetic-energy cutoff of 100 Ry for the
plane wave and norm-conserving pseudopotentials. The norm-conserving
pseudopotential for Pd is constructed by including semicore $4s$, $4p$, and
$4d$ electrons as valence electrons. For Se, we include $4s$, $4p$, and $4d$
electrons as valence electrons. For self-consistent calculations, we use an
$8{\times}8{\times}6$ $k$-point sampling in the three-dimensional (3D)
Brillouin zone (BZ) for bulk PdSe2 and an $8{\times}8{\times}1$ $k$-point
sampling in the 2D BZ for monolayer and bilayer PdSe2. The density of states
(DOS) is calculated with a finer $k$ grid of $20{\times}20{\times}15$ for bulk
and $40{\times}40{\times}1$ for 2D layers.

###  II.3 GW calculations

We performed GW calculations for QP band structures in monolayer, bilayer, and
bulk PdSe2 using the BERKELEYGW code [28, 18, 29]. We used the one-shot GW
method (G0W0) which constructs the dielectric matrix and the Green’s function
from DFT eigenvalues and wave functions and calculates the self-energy once.
The generalized plasmon-pole model is used to consider the frequency
dependence of the inverse dielectric function $\epsilon^{-1}(\omega)$.

In GW calculations, summation over unoccupied bands is needed to obtain the
self-energy. Generally, it is known that a large number of unoccupied bands
and a large energy cutoff for the dielectric matrix are required in GW
calculations involving semicore $d$ states, while simple $sp$ semiconductors
show fast convergence [30]. For calculations of PdSe2, Pd has fully occupied
$4d$ states, so its semicore $4s$ and $4p$ states should be considered as
valence electrons in pseudopotentials. We include 3900 unoccupied bands in our
GW calculations for bulk PdSe2, which is about 65 times the number of occupied
bands, and 1050 (2100) unoccupied bands for monolayer (bilayer) PdSe2, which
is about 35 times the number of occupied bands. For better convergence, we
used the static remainder method that considers the contribution of unincluded
unoccupied bands to the expectation value of the static Coulomb-hole operator
[31]. Kinetic-energy cutoffs of 40 and 20 Ry are used to calculate dielectric
matrices of the bulk and layer, respectively. A $4{\times}4{\times}3$
$q$-point sampling is used for the bulk, and for monolayer and bilayer PdSe2,
a $6{\times}6{\times}1$ uniform $q$ grid is used along with additional 10 $q$
points determined by the nonuniform neck subsampling method [32], which is
equivalent to a $1143{\times}1143{\times}1$ uniform $q$ grid. With these
calculation parameters, QP band gaps are converged within a few tens of meV.

In general, VBMs and CBMs may not be at a high-symmetry point in BZ. Thus, to
find the QP band gap, which can be indirect, we interpolate the QP band
structure using the dual-grid method [29]. In the dual-grid method, the
expectation value of an observable for a wave function in a fine $k$ grid is
calculated by summing its expectation values for wave functions in a coarse
$k$ grid weighted by coefficients expanding the fine-$k$-grid wave function
with coarse-$k$-grid wave functions. QP DOS is also calculated by
interpolating QP energies from a coarse $k$ grid to a fine $k$ grid. In our
present work, the coarse $k$ grid is $4{\times}4{\times}3$ for bulk and
$6{\times}6{\times}1$ for 2D layers, and the fine $k$ grid is
$20{\times}20{\times}15$ for bulk and $40{\times}40{\times}1$ for 2D layers

##  III RESULTS AND DISCUSSIONS

We performed DFT and GW calculations of monolayer, bilayer, and bulk PdSe2 to
find direct and indirect band gaps, effective masses, work functions,
ionization potentials, and electron affinities. We also obtained the projected
density of states to analyze contributions of Pd and Se orbitals to electronic
states near VBMs and CBMs. Finally, we obtained equienergy lines at 0.1 eV
below the VBM and 0.1 eV above the CBM to find the shape of hole and electron
pockets in the $k$-space.

###  III.1 DFT band structures

Figure 2:  Band structures of (a) monolayer, (b) bilayer, and (c) bulk PdSe2
from DFT calculations, and (d) 3D BZ of bulk PdSe2 and projected 2D BZ of
monolayer and bilayer PdSe2. The $M$ point is not a high-symmetry point, but
it is $(t,0.50,0.00)$, where $t$ is $0.46$, $0.40$, and $0.42$ for monolayer,
bilayer, and bulk PdSe2, respectively. In (a)-(c), green arrows connect the
VBM and CBM for direct and indirect band gaps.

Figure 2 shows electronic band structures for monolayer, bilayer, and bulk
PdSe2 obtained by DFT calculations. Some VBMs and CBMs are located away from
high-symmetry lines in the BZ. Band gaps in monolayer, bilayer, and bulk PdSe2
are all indirect, unlike other TMDs in which the band gap is direct or
indirect, depending on the number of layers [1, 4].

Monolayer PdSe2 has an indirect band gap of 1.33 eV, with the VBM at $(\pm
0.34,0,0)$ and the CBM at $(\pm 0.44,\pm 0.48,0)$ [Fig. 2(a)]. Here, $k$-point
coordinates are with respect to reciprocal lattice vectors. The work function
is 4.71 eV, which we estimate as the energy between the vacuum level and the
center of the indirect band gap because the chemical potential of an undoped
semiconductor at low temperature is equal to the average of CBM and VBM
energies [5]. The ionization potential is 5.37 eV, and the electron affinity
is 4.04 eV. Near the VBM, the effective mass is $0.8m_{e}$ along the $x$
direction and $1.3m_{e}$ along the $y$ direction. Here, $m_{e}$ is the bare
mass of an electron in vacuum. Near the CBM, the effective mass is $4.6m_{e}$
along the $x$ direction and $0.5m_{e}$ along the $y$ direction. The direct
band gap is 1.50 eV, located at $(\pm 0.36,0,0)$ [Fig. 2(a)]. Our DFT indirect
band gap of monolayer PdSe2 agrees well with results of previous DFT
calculations [20, 13].

Bilayer PdSe2 has the VBM at $(\pm 0.22,0,0)$ and the CBM at $(\pm 0.36,\pm
0.43,0)$, exhibiting an indirect band gap of 0.79 eV [Fig. 2(b)]. Near the
VBM, the valence band shows almost flat band dispersion along the $x$
direction. The VBM at $(\pm 0.22,0,0)$ is 5 $\mu$eV higher than the energy of
the highest valence band at the $\Gamma$ point. In the line from $(\pm
0.22,0,0)$ to the $\Gamma$ point, the highest valence band has a local
minimum, and it is about 4 meV lower than the VBM. Near the VBM, the effective
mass is $2.2m_{e}$ along the $x$ direction and $1.4m_{e}$ along the $y$
direction. Near the CBM, the effective mass is $0.6m_{e}$ along the $x$
direction and $0.3m_{e}$ along the $y$ direction. The work function is 4.65
eV, the ionization potential is 5.06 eV, and the electron affinity is 4.27 eV.
The direct band gap is 1.05 eV at $(\pm 0.28,0,0)$ [Fig. 2(b)].

Bulk PdSe2 is semimetallic with a negative indirect band gap of $-$0.10 eV
[Fig. 2(c)]; that is, the VBM is higher than the CBM by 0.10 eV. The VBM is at
the $\Gamma$ point, and the CBM is at $(\pm 0.37,\pm 0.44,0)$. The effective
mass near the VBM is $4.1m_{e}$ along the $x$ direction and $2.0m_{e}$ along
the $y$ direction. The effective mass near the CBM is $0.3m_{e}$ along the $x$
direction and $0.3m_{e}$ along the $y$ direction. The direct band gap is 0.56
eV at $(\pm 0.25,0,0)$ [Fig. 2(c)]. We obtained band energies of bulk PdSe2
with respect to the vacuum level by extracting the vacuum level from DFT
calculation of 12-layer PdSe2, which is thick enough to have a band gap very
close to that of bulk. This method was used in the study of bulk TMDs [5].
From band energies with respect to the vacuum level, the work function of bulk
PdSe2 is estimated to be 4.55 eV.

From the resistivity measurement, bulk PdSe2 is a semiconductor with a band
gap of 0.4 eV [11], while our DFT calculation yields a negative indirect band
gap of $-0.1$ eV. This underestimation of the band gap in PdSe2 in DFT
calculations may originate from overestimation of the interlayer interaction.
Previous theoretical calculations reported that, without van der Waals forces,
the relaxed out-of-plane lattice constant $c$ is 10% larger than the
experimental value [9], and with this lattice constant, the DFT calculation
gives a band gap of 0.4 eV [9], which is close to the experiment although the
atomic structure is not. A similar tendency has also been reported for black
phosphorus [33].

###  III.2 QP band structures from GW calculations

We obtained QP band structures of PdSe2 from GW calculations. Since the CBM is
not in a high-symmetry line in the BZ, it is necessary to calculate band-edge
energies through the dual-grid interpolation in order to determine the
indirect band gap. We found the VBM and CBM of QP band structures by
interpolating the QP correction near locations of conduction and valence band
edges in DFT band structures. Then we determined the chemical potential using
interpolated QP band structures. In valence and conduction bands of all the
considered cases, obtained GW corrections are fairly smooth with energy, and
QP bandwidths are slightly larger than DFT bandwidths. Thus, effective masses
in QP band structures are slightly smaller than those in DFT band structures,
as shown in Table 1. GW corrections are almost constant shifts near conduction
and valence band edges, and band gaps are corrected mainly by GW corrections
for valence bands.

Table 1:  Indirect band gaps $E^{(i)}_{g}$, direct band gaps $E^{(d)}_{g}$, effective masses, ionization potentials (IPs), electron affinities (EAs), and work functions $W$ of monolayer (ML), bilayer (BL), and bulk PdSe2 from DFT and GW calculations. Effective masses $m^{*}_{e,x}$, $m^{*}_{e,y}$, $m^{*}_{h,x}$, and $m^{*}_{h,y}$ are for electrons and holes along the $x$ and $y$ axes, respectively. Band gaps, IP, EA, and $W$ are in eV. Effective masses are in units of the bare electron mass in vacuum.  |  | $E^{(i)}_{g}$ | $E^{(d)}_{g}$ | $m^{*}_{e,x}$ | $m^{*}_{e,y}$ | $m^{*}_{h,x}$ | $m^{*}_{h,y}$ | IP | EA | $W$
---|---|---|---|---|---|---|---|---|---|---
ML | DFT | 1.33 | 1.50 | 4.6 | 0.5 | 0.8 | 1.3 | 5.37 | 4.04 | 4.71
GW | 2.37 | 2.54 | 4.3 | 0.5 | 0.8 | 1.3 | 6.44 | 4.06 | 5.26
BL | DFT | 0.79 | 1.05 | 0.6 | 0.3 | 2.2 | 1.4 | 5.06 | 4.27 | 4.65
GW | 1.55 | 1.81 | 0.6 | 0.3 | 2.1 | 1.4 | 6.04 | 4.48 | 5.26
bulk | DFT |  $-$0.10 | 0.56 | 0.3 | 0.3 | 4.1 | 2.0 |  |  | 4.55
GW | 0.45 | 1.10 | 0.3 | 0.3 | 3.9 | 1.9 | 5.41 | 4.96 | 5.18

Figure 3(a) shows the QP band structure of monolayer PdSe2 obtained from GW
calculations. Monolayer PdSe2 has an indirect QP band gap of 2.37 eV, with the
VBM at $(\pm 0.34,0,0)$ and the CBM at $(\pm 0.39,\pm 0.47,0)$. GW correction
is much larger in monolayer PdSe2 than in the bulk. For monolayer PdSe2, a
large band-gap correction leads to a large work-function change. While the
work function from DFT calculation is 4.71 eV for the monolayer, GW
calculations increase this value to 5.26 eV for the monolayer. The ionization
potential is 6.44 eV, and the electron affinity is 4.06 eV. The GW correction
for band-edge energies is prominent in valence bands. Near the VBM, the
effective mass is $0.8m_{e}$ along the $x$ direction and $1.3m_{e}$ along the
$y$ direction. Near the CBM, the effective mass is $4.3m_{e}$ along the $x$
direction and $0.5m_{e}$ along the $y$ direction. The effective masses are
slightly reduced by GW corrections when compared with DFT results.

Figure 3: QP band structures of (a) monolayer, (b) bilayer, and (c) bulk PdSe2
from GW calculations and (d) 3D BZ of bulk PdSe2 and projected 2D BZ of
monolayer and bilayer PdSe2. The $M$ point is not a high-symmetry point, but
it is $(t,0.50,0.00)$, where $t$ is $0.42$, $0.40$, and $0.40$ for monolayer,
bilayer, and bulk PdSe2, respectively. In (a)-(c), green arrows connect the
VBM and CBM for direct and indirect band gaps.

Figure 3(b) shows the QP band structure of bilayer PdSe2 obtained from GW
calculations. Bilayer PdSe2 has an indirect QP band gap of 1.55 eV, with the
VBM at $(\pm 0.24,0,0)$ and the CBM at $(\pm 0.34,\pm 0.43,0)$. GW correction
is much larger in bilayer PdSe2 than in the bulk. For bilayer PdSe2, a large
band-gap correction leads to a large work-function change. While the work
function from the DFT calculation is 4.65 eV for the bilayer, the GW
calculation increases this value to 5.26 eV for the bilayer, which is the same
value as that for the monolayer. The ionization potential is 6.04 eV, and the
electron affinity is 4.48 eV. The GW correction for band-edge energies is
prominent in the valence bands as in the monolayer. Near the VBM, the
effective mass is $2.1m_{e}$ along the $x$ direction and $1.4m_{e}$ along the
$y$ direction. Near the CBM, the effective mass is $0.6m_{e}$ along the $x$
direction and $0.3m_{e}$ along the $y$ direction. The effective masses are
slightly reduced by GW corrections when compared with DFT results as in the
monolayer.

Figure 3(c) shows the QP band structure of bulk PdSe2 obtained from GW
calculations. Bulk PdSe2 has an indirect QP band gap of 0.45 eV, with the VBM
at $\Gamma$ and the CBM at $(\pm 0.35,\pm 0.44,0)$. The calculated QP band gap
of 0.45 eV is in excellent agreement with the experimental band gap of 0.4 eV
from the resistivity measurement [11]. The GW correction in bulk PdSe2 is
about $-0.4$ and $-0.9$ eV for conduction and valence bands, respectively,
increasing the band gap by about 0.5 eV and increasing the work function. The
GW correction in bulk PdSe2 is much smaller than in the monolayer and bilayer,
but the correction is large enough to change the semimetallic DFT band
structure to the semiconducting one. The work function is 5.18 eV, the
ionization potential is 5.41 eV, and the electron affinity is 4.96 eV. Near
the VBM, the effective mass is $3.9m_{e}$ along the $x$ direction and
$1.9m_{e}$ along the $y$ direction. Near the CBM, the effective mass is
$0.3m_{e}$ along the $x$ direction and $0.3m_{e}$ along the $y$ direction. The
effective masses are slightly reduced by GW corrections when compared with DFT
results as in the monolayer and bilayer.

Figure 4: QP DOS of (a) monolayer, (b) bilayer, and (c) bulk PdSe2 from GW
calculations. QP energies are interpolated to a dense $k$ grid to calculate
DOS. States near the CBM and VBM are mainly from Pd $4d$ and Se $4p$ orbitals.

###  III.3 Density of states

Figure 4 shows the QP DOS for monolayer, bilayer, and bulk PdSe2 obtained from
GW calculations. We obtained QP band energies for a coarse $k$-point grid by
GW calculations and then interpolated them to a fine $k$-point grid to
calculate the QP DOS. In Figs. 4(a)-4(c), the VBM is set to zero in energy.
Projected QP DOSs show that states near the VBM and CBM are mainly from Pd
$4d$ and Se $4p$ orbitals, and contributions of Se $4p$ orbitals to states
near the VBM and CBM are greater than those of Pd $4d$ orbitals. This implies
that each Se or Pd atom contributes similarly to the QP DOS near the VBM and
CBM, for Se atoms are twice as numerous as Pd atoms in the material. Thus,
hybridization of Pd $4d$ orbitals and Se $4p$ orbitals is important for
describing the band gap in PdSe2.

In Fig. 4(a), the QP DOS of monolayer PdSe2 shows a very sharp peak right
below the VBM originating from almost flat bands near the VBM shown in Fig.
3(a). In Fig. 4(b), the QP DOS of bilayer PdSe2 shows a peak near the VBM
which is not as strong as that in the monolayer but closer to the VBM in
energy. In contrast, the QP DOS of bulk PdSe2, shown in Fig. 4(c), does not
have a peak near the VBM. Thus, with hole doping, ferromagnetic instability
may be less effective in the bulk than in the monolayer and bilayer.

Figure 5: Equienergy lines at 0.1 eV below the VBM in (a) monolayer, (c)
bilayer, and (e) and (g) bulk, and at 0.1 eV above the CBM in (b) monolayer,
(d) bilayer, and (f) and (h) bulk obtained with GW calculations. In (a)-(f),
equienergy lines are plotted in the $k_{x}k_{y}$ plane. In (g), equienergy
lines are plotted in the $k_{y}=0$ plane. In (h), equienergy lines are plotted
in the $k_{y}=0.44$ plane. In each plot, red dots represent locations of the
VBM or CBM.

###  III.4 Equienergy lines in the k space

Figure 5 shows equienergy lines in the BZ obtained from GW calculations of
monolayer, bilayer, and bulk PdSe2 at energies which are 0.1 eV below the VBM
and 0.1 eV above the CBM. These lines correspond to hole or electron pockets
in the BZ if hole or electron doping is done.

Monolayer PdSe2 has the VBM at $(\pm 0.34,0,0)$ and the CBM at $(\pm 0.39,\pm
0.47,0)$, as marked with red dots in Figs. 5(a) and 5(b), respectively. In
monolayer PdSe2, equienergy lines 0.1 eV below the VBM form four close loops,
two of which (marked with $\alpha$) are around the VBM and larger than the
other two (marked with $\beta$) [Fig. 5(a)]. Meanwhile, equienergy lines 0.1
eV above the CBM form loops near four corners of the BZ [Fig. 5(b)], which are
two intersecting slanted ellipses (marked with $\gamma_{1}$ and $\gamma_{2}$)
when all BZ corners are put together.

Bilayer PdSe2 has the VBM at $(\pm 0.24,0,0)$ and the CBM at $(\pm 0.34,\pm
0.43,0)$, as marked with red dots in Figs. 5(c) and 5(d), respectively. In
bilayer PdSe2, the equienergy line 0.1 eV below the VBM is a large loop
(marked with $\alpha$) surrounding two VBM points [Fig. 5(c)], indicating that
its valence band is quite flat near the VBM. Meanwhile, equienergy lines 0.1
eV above the CBM are loops near four corners of the BZ, which form four closed
loops (marked with $\gamma_{1}$, $\gamma_{2}$, $\gamma_{3}$, and $\gamma_{4}$)
when horizontal BZ boundaries are put together [Fig. 5(d)]. Compared with the
monolayer, VBM points are shifted a little bit toward the $\Gamma$ point, and
CBM points are shifted a little bit away from BZ corners.

Bulk PdSe2 has the VBM at $\Gamma$ and the CBM at $(\pm 0.35,\pm 0.44,0)$, as
marked with red dots in Figs. 5(e) and 5(f), respectively. In bulk PdSe2, the
equienergy line 0.1 eV below the VBM is a closed loop (marked with $\alpha$)
surrounding the $\Gamma$ point in the $k_{x}k_{y}$ plane [Fig. 5(e)],
indicating that the in-plane effective mass is large. Meanwhile, equienergy
lines 0.1 eV above the CBM are loops near four corners of the BZ, which form
four closed loops (marked with $\gamma_{1}$, $\gamma_{2}$, $\gamma_{3}$, and
$\gamma_{4}$) when horizontal BZ boundaries are put together [Fig. 5(f)],
similar to the case of the bilayer.

Since the BZ of the bulk is three-dimensional, we also plot the equienergy
line 0.1 eV below the VBM in the $k_{x}k_{z}$ plane, that is, the $k_{y}=0$
plane [Fig. 5(g)], and equienergy lines 0.1 eV above the CBM in the
$k_{y}=0.44$ plane [Fig. 5(h)]. In the $k_{y}=0$ plane, the equienergy line
0.1 eV below the VBM is a very narrow ellipselike loop corresponding to
$\alpha$ [Fig. 5(g)], indicating that the effective mass is much smaller along
the $z$ direction than the $x$ direction. In the $k_{y}=0.44$ plane which
passes through two CBM points, equienergy lines 0.1 eV above the CBM are two
ellipselike loops corresponding to $\gamma_{1}$ and $\gamma_{3}$ [Fig. 5(h)],
indicating that the effective mass is larger along the $z$ direction than the
$x$ direction.

##  IV SUMMARY

We studied the electronic structures of monolayer, bilayer, and bulk PdSe2
using first-principles DFT and GW calculations. First, we calculated the
electronic structure of bulk PdSe2 using DFT, obtaining a semimetallic
electronic structure due to underestimation of the band gap. To calculate the
band gap of PdSe2 correctly, we performed GW calculations. We obtained the QP
band structure of bulk PdSe2 with an indirect band gap of 0.45 eV, which is
consistent with the reported experimental band gap of 0.4 eV. Then we
calculated the electronic structure of 2D PdSe2 with respect to the vacuum
level using DFT and the GW method. Our DFT calculations produced band gaps of
1.33 and 0.79 eV for monolayer and bilayer PdSe2, respectively. From GW
calculations, QP band gaps of monolayer and bilayer PdSe2 are 2.37 and 1.55
eV, respectively. We also obtained the DOS, effective masses, work functions,
ionization potentials, electron affinities, and $k$-space shapes of electron
and hole pockets for monolayer, bilayer, and bulk PdSe2 from DFT and GW
calculations. These results provide basic information for the development of
material properties and device applications.

ACKNOWLEDGMENTS

This work is supported by NRF of Korea (Grant No. 2020R1A2C3013673).
Computational resources were provided by the KISTI Supercomputing Center
(Project No. KSC-2020-CRE-0335).

APPENDIX: DOS FROM DFT CALCULATIONS

Figure 6: DOS of (a) monolayer, (b) bilayer, and (c) bulk PdSe2 from DFT
calculations. States near the CBM and VBM are mainly from Pd $4d$ and Se $4p$
orbitals.

Figure 6 shows the electronic DOS of monolayer, bilayer, and bulk PdSe2
obtained from DFT calculations. In Figs. 6(a) and 6(b), VBM is set to zero in
energy. In Fig. 6(c), the chemical potential is set to zero. Overall shapes of
the electronic DOS are not much different from those of the QP DOS from GW
calculations, except for sizes of band gaps. Compared with the QP DOS from GW
calculations, the DOSs from DFT calculations show smaller band gaps in the
monolayer and bilayer [Figs. 6(a) and 6(b)] and a negative indirect band gap
in the bulk [Fig. 6(c)]. Because of the negative indirect band gap in DFT
calculations, DOS in Fig. 6(c) does not show any band gap explicitly, although
bulk PdSe2 has a finite energy difference between valence and conduction bands
throughout the BZ, as shown in Fig. 2(c). Meanwhile, similar to the QP DOS
from GW calculations, the DOSs from DFT calculations show sharp peaks near the
VBM in the monolayer and bilayer, as shown in Figs. 6(a) and 6(b), but no peak
near the VBM in the bulk [Fig. 6(c)].

## References

  * [1] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically Thin MoS2: A New Direct-Gap Semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * [2] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Emerging photoluminescence in monolayer MoS2, Nano Lett. 10, 1271 (2010).
  * [3] L. Britnell, R. M. Ribeiro, A. Eckmann, R. Jalil, B. D. Belle, A. Mishchenko, Y.-J. Kim, R. V. Gorbachev, T. Georgiou, S. V. Morozov, A. N. Grigorenko, A. K. Geim, C. Casiraghi, A. H. Castro Neto, and K. S. Novoselov, Strong light-matter interactions in heterostructures of atomically thin films, Science 340, 1311 (2013).
  * [4] Y. Zhang, T.-R. Chang, B. Zhou, Y.-T. Cui, H. Yan, Z. Liu, F. Schmitt, J. Lee, R. Moore, Y. Chen, H. Lin, H.-T. Jeng, S.-K. M., Z. Hussain, A. Bansil, and Z.-X. Shen, Direct observation of the transition from indirect to direct bandgap in atomically-thin epitaxial MoSe2. Nat. Nanotechnol. 9, 111 (2014).
  * [5] H.-g. Kim and H. J. Choi, Thickness dependence of work function, ionization energy, and electron affinity of Mo and W dichalcogenides from DFT and GW calculations, Phys. Rev. B 103, 085404 (2021).
  * [6] S. Oh, J. Y. Lim, S. Im, and H. J. Choi, Stability, efficiency, and mechanism of $n$-type doping by hydrogen adatoms in two-dimensional transition metal dichalcogenides, Phys. Rev. B 100, 085416 (2019).
  * [7] J.-S. Park and H. J. Choi, Band-gap opening in graphene: A reverse-engineering approach, Phys. Rev. B 92, 045402 (2015).
  * [8] W. L. Chow, P. Yu, F. Liu, J. Hong, X. Wang, Q. Zeng, C.-H. Hsu. C. Zhu, J. Zhou, X. Wang, J. Xia, J. Yan, Y. Chen , D. Wu, T. Yu, Z. Shen, H. Lin, C. Jin, B. K. Tay, and Z. Liu, High Mobility 2D Palladium Diselenide Field effect Transistors with Tunable Ambipolar Characteristics, Adv. Mater. 29, 1602969 (2017).
  * [9] A. D. Oyedele, S. Yang, L. Liang, A. A. Puretzky, K. Wang, J. Zhang, P. Yu, P. R. Pudasaini, A. W. Ghosh, Z. Liu, C. M. Rouleau, B. G. Sumpter, M. F. Chisholm, W. Zhou, P. D. Rack, D. B. Geohegan, and K. Xiao, PdSe2: Pentagonal Two-Dimensional Layers with High Air Stability for Electronics, J. Am. Chem. Soc. 139, 14090 (2017).
  * [10] S.-H. Zhang, and B.-G. Liu, Hole-doping-induced half-metallic ferromagnetism in a highly-air-stable PdSe2 monolayer under uniaxial stress, J. Mater. Chem. C. 6, 6792 (2018).
  * [11] F. Hulliger, Electrical properties of some nickel-group chalcogenides J. Phys. Chem. Solids 26, 639 (1965).
  * [12] A. Hamidani, B. Bennecer, and K. Zanat, Structural and electronic properties of the pseudo-binary compounds PdX2 (X = P, S and Se), J. Phys. Chem. Solids 71, 42 (2010).
  * [13] J. Sun, H. Shi, T. Siegrist, and D. J. Singh, Electronic, transport, and optical properties of bulk and mono-layer PdSe2, Appl. Phys. Lett. 107, 153902 (2015).
  * [14] L. Hedin, New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem, Phys. Rev. 139, A796 (1965).
  * [15] G. Strinati, H. J. Mattausch, and W. Hanke, Dynamical Correlation Effects on the Quasiparticle Bloch States of a Covalent Crystal, Phys. Rev. Lett. 45, 290 (1980).
  * [16] G. Strinati, H. J. Mattausch, and W. Hanke, Dynamical aspects of correlation corrections in a covalent crystal, Phys. Rev. B 25, 2867 (1982).
  * [17] M. S. Hybertsen and S. G. Louie, First-Principles Theory of Quasiparticles: Calculation of Band Gaps in Semiconductors and Insulators, Phys. Rev. Lett. 55, 1418 (1985).
  * [18] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * [19] M. van Schilfgaarde, T. Kotani, and S. Faleev, Quasiparticle Self-Consistent GW Theory, Phys. Rev. Lett. 96, 226402 (2006).
  * [20] S. Lebègue, T. Björkmank, M. Klintenberg, R. M. Nieminen, and O. Eriksson, Two-Dimensional Materials from Data Filtering and Ab initio Calculations, Phys. Rev. X 3, 031002 (2013).
  * [21] A. V. Kuklin and H. Ågren, Quasiparticle electronic structure and optical spectra of single-layer and bilayer PdSe2: Proximity and defect-induced band gap renormalization, Phys. Rev. B 99, 245114 (2019).
  * [22] M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, F. Wang, S. G. Louie, and M. F. Crommie, Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor, Nat. Mater. 13, 1091 (2014).
  * [23] C. Soulard, X. Rocquefelte, P.-E. Petit, M. Evain, S. Jobic, J.-P. Itié, P. Munsch, H.-J. Koo, and M.-H. Whangbo, Experimental and Theoretical Investigation on the Relative Stability of the PdS2\- and Pyrite-Type Structures of PdSe2, Inorg. Chem. 43, 1943 (2004).
  * [24] P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136, B864 (1964).
  * [25] W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140, A1133 (1965).
  * [26] P. Giannozzi et al., QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * [27] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * [28] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures, Comput. Phys. Commun. 183, 1269 (2012).
  * [29] M. Rohlfing and S. G. Louie, Electron-hole excitations and optical spectra from first principles, Phys. Rev. B 62, 4927 (2000).
  * [30] B.-C. Shih, Y. Xue, and P. Zhang, Quasiparticle Band Gap of ZnO: High Accuracy from the Conventional $G^{0}W^{0}$ Approach, Phys. Rev. Lett. 105, 146401 (2010).
  * [31] J. Deslippe, G. Samsonidze, M. Jain, M. L. Cohen, and S. G. Louie, Coulomb-hole summations and energies for GW calculations with limited number of empty orbitals: A modified static remainder approach, Phys. Rev. B 87, 165124 (2013).
  * [32] F. H. da Jornada, D. Y. Qiu, and S. G. Louie, Nonuniform sampling schemes of the Brillouin zone for many-electron perturbation-theory calculations in reduced dimensionality, Phys. Rev. B 95, 035109 (2017).
  * [33] S. S. Baik, K. S. Kim, Y. Yi, and H. J. Choi, Emergence of Two-Dimensional Massless Dirac Fermions, Chiral Pseudospins, and Berry’s Phase in Potassium Doped Few-Layer Black Phosphorus, Nano Lett. 15, 7788 (2015).