# Prototypes as Explanation for Time Series Anomaly Detection

Bin Li<EMAIL_ADDRESS>0000-0002-9707-4596 TU Dortmund
UniversityDortmundGermany ,  Carsten Jentsch<EMAIL_ADDRESS>TU Dortmund UniversityDortmundGermany and  Emmanuel Müller
<EMAIL_ADDRESS>0000-0002-5409-6875 TU Dortmund
UniversityDortmundGermany

(2022)

###### Abstract.

Detecting abnormal patterns that deviate from a certain regular repeating
pattern in time series is essential in many big data applications. However,
the lack of labels, the dynamic nature of time series data, and unforeseeable
abnormal behaviors make the detection process challenging. Despite the success
of recent deep anomaly detection approaches, the mystical mechanisms in such
black-box models have become a new challenge in safety-critical applications.
The lack of model transparency and prediction reliability hinders further
breakthroughs in such domains. This paper proposes ProtoAD, using prototypes
as the example-based explanation for the state of regular patterns during
anomaly detection. Without significant impact on the detection performance,
prototypes shed light on the deep black-box models and provide intuitive
understanding for domain experts and stakeholders. We extend the widely used
prototype learning in classification problems into anomaly detection. By
visualizing both the latent space and input space prototypes, we intuitively
demonstrate how regular data are modeled and why specific patterns are
considered abnormal.

Anomaly detection, Anomaly explanation, Prototypes, Time series

††copyright: acmcopyright††journalyear: 2022††conference: The 2nd KDD Workshop
on Anomaly and Novelty Detection, Explanation and Accommodation; August 14,
2022; Washington, DC††ccs: Computing methodologies Anomaly detection

##  1\. Introduction

Anomaly detection in time series data is gaining traction in several current
big data research areas. The vast development of data processing and analysis
facilitates real-time data monitoring applications in different branches,
e.g., health care, manufacturing, astronomy (Hundman et al., 2018; Ren et al.,
2019; Ming et al., 2019). Anomaly detection in time series data is an
important application area in those scenarios. Unlike traditional outlier
detection in the database system, detecting single abnormal points or even
whole abnormal time periods in time series data is much more challenging due
to its serial dependent and dynamic nature. On the one hand, time series data
is usually collected from sensor networks in real-time, leaving less time for
domain experts to generate labels for the model training. On the other hand,
anomalies often only appear on specific dimensions (subspace anomalies) and
specific temporal contexts (contextual anomalies), making detecting and
interpreting such anomalies even harder.

In recent years, deep reconstruction-based approaches have been developed for
time series anomaly detection tasks, for example, the autoencoders (Malhotra
et al., 2016; Zong et al., 2018). An autoencoder is a symmetric neural network
trained to reconstruct the regular111To prevent ambiguity, normal data in the
anomaly detection context is referred to as regular data in this paper. data
and is supposed to produce significantly larger reconstruction errors by
unseen anomalies. Specifically for time series data, recurrent neural networks
(RNNs) are used to construct the autoencoders to capture the temporal
information in the input sequences. Thanks to their model capability, deep
autoencoders have shown convincing performance in detecting high dimensional
and temporal anomalies from time series (Malhotra et al., 2016).

Despite the performance advantage, autoencoders also suffer the criticisms of
other deep neural networks - the lack of transparency. Though the detected
anomalies fit well with standard evaluation criteria, their reliability is not
well-proven. The lack of human interpretable information makes it hard to tell
why an anomaly is abnormal, especially in high-dimensional time series data or
long input sequences. In safety-crucial applications, for instance, a machine
learning model assists in detecting abnormal arrhythmia in patient ECG data.
Without a verifiable interpretation, it might lead to catastrophic
consequences by automatic treatment. Therefore, interpretability is an
essential demand of such anomaly detection applications.

To address the black-box issue, various interpretable machine learning
approaches recently gained attention in the fields that require high
transparency and human-understandability (Schlegel et al., 2020). We propose
to use example-based prototypes as an intuitive and explainable solution for
interpreting anomalies in time series data. Prototypes are widely used for
case-based reasoning in computer vision (Li et al., 2018; Chen et al., 2018;
Hase et al., 2019), graph learning (Zhang et al., 2021) and sequential data
learning (Ni et al., 2021; Ming et al., 2019; Gee et al., 2019). Embedding a
prototype layer is one of the widely used approaches to learn prototypes from
neural networks in an end-to-end fashion. With the build-in prototype layer in
the neural networks, prototype-based models are commonly efficient to train
without requiring extra investigation in the interpretability functionalities.
Moreover, the prototypes are usually self-contained and straightforward to
understand, e.g., representative animal faces, sentences, and sensor data
patterns. However, prototypes are still understudied in the anomaly detection
field.

In this paper, we propose to use prototypes to interpret the regular data
during anomaly detection using autoencoders. In this context, data showing a
certain repeating regular pattern is generated by several latent
distributions, while the anomaly is any data point or period that deviates
from from this regular pattern. We model the regular patterns of the time
series data with prototypes in the latent space of the autoencoder and learn
multiple prototypes to discover the latent components of the regular data
distribution. Anomaly patterns that lie distantly from the data generated
during the regular pattern state can be explained by comparing them with the
constructed prototypes.

Moreover, we further explain the anomaly patterns that lie distantly from
regular data in the latent space by comparing them with the prototypes. To our
knowledge, this is the first prototype-based explanation application in the
time series anomaly detection domain.

Our contribution to the paper can be summarized as follows:

  1. (1)

we propose ProtoAD, an end-to-end LSTM-Autoencoder for anomaly detection with
prototype learning

  2. (2)

we develop latent space prototype-based explanations for the understanding of
the regular state of the studied data

  3. (3)

we evaluate our method with synthetic and real-world time series data.
Moreover, we visually demonstrate examples of prototypes to show the benefit
of our model qualitatively

##  2\. Related works

This section briefly reviews the existing work in autoencoder-based anomaly
detection models and prototype-based explanations methods.

###  2.1. Reconstruction-based anomaly detection

Autoencoders have been used as an unsupervised anomaly detection approach for
years. Feed-forward autoencoders (Zong et al., 2018) and Variational
Autoencoders (Xu et al., 2018) are used for time-independent data. In
contrast, RNN-based autoencoders (Malhotra et al., 2016) show their strength
in detecting contextual anomalies in time series data. Based on the
reconstruction error, a standard approach for estimating anomaly likelihood is
to assume the reconstruction error following a normal distribution and measure
the Mahalanobis distance between the reconstruction error of unknown data and
the estimated distribution (Malhotra et al., 2016). In addition to
reconstruction error, the hidden representation in the latent space can also
be used for likelihood estimation (Zong et al., 2018). Gaussian Mixture Model
(GMM) (Zong et al., 2018) and energy-based model (Zhai et al., 2016) are also
used for the likelihood estimation. Common thresholding techniques over the
anomaly likelihood are based on maximizing the performance on a validation
set, which requires labels in advance (Malhotra et al., 2016).

Other approaches, including the hierarchical temporal memory (HTM) (Ahmad et
al., 2017) and temporal convolutional network (TCN) (He and Zhao, 2019) are
also adopted in time series anomaly detection concerning different use cases
and data properties. However, they are not directly relevant to the
reconstruction-based models.

###  2.2. Explanation with prototypes

Due to the complex properties of both feature and time dimensions of time
series data, prototypes are considered an intuitive explanation. Common
prototype learning approaches for neural networks follow a three-step
paradigm. 1) Representation learning, 2) prototype learning in the latent
space, and 3) class prediction. The objective commonly includes 1) minimizing
classification error, 2) minimizing the distances between each hidden
representation and one of the prototypes, and 3) maximizing the distances
between prototypes.

In the existing prototype learning literature, (Li et al., 2018) employs a
multi-layer convolutional neural network to construct the autoencoder, which
learns hidden representations for image data. They rely on the decoder to
project the learned prototypes in the human-understandable space, sometimes
producing unrealistic reconstructions. Using a single encoder to replace the
autoencoder is considered as a reduction of training effort in (Chen et al.,
2018; Ming et al., 2019), and they use the nearest neighborhood of each
prototype in the latent space as the corresponding realistic patterns in the
input space. Under different problem settings, (Chen et al., 2018) and (Hase
et al., 2019) build up the encoder with convolutional neural networks to
encode image data, (Ming et al., 2019) uses RNNs for sequential data, (Ni et
al., 2021) use a convolutional layer to learn time series representations,
(Zhang et al., 2021) employs graph neural networks for the encoder. In our
work, we use the single LSTM-Autoencoder for both reconstruction-based time
series anomaly detection and hidden space representation learning.

The standard objective functions of existing prototype learning approaches
consist of multiple regularisation terms that are trained jointly. To ensure
the representation ability of the prototypes, (Li et al., 2018; Ming et al.,
2019; Gee et al., 2019; Zhang et al., 2021; Chen et al., 2018) all minimize
the distance between each prototype and every nearby hidden representation as
well as every hidden representation to each prototype. Furthermore, the
learned prototypes are supposed to be diverse from each other (Ming et al.,
2019; Zhang et al., 2021; Gee et al., 2019). In the objective, we will follow
the standard design of the regularization terms above. However, different from
most existing works, which use cross-entropy for their classification task to
minimize the classification error (Li et al., 2018; Ming et al., 2019; Hase et
al., 2019), in our unsupervised setting, we use the reconstruction-error to
regularize the reconstruction process of regular data.

Besides prototypes, other techniques are also used for explaining time series
data. The representative subsequences Shapelets (Kidger et al., 2020; Li et
al., 2020) can be similarly used for explanation. Instead of finding the
representative pattern as prototypes, counterfactuals (Ates et al., 2021)
explain the instance towards the opposite class. Recently, the attention
mechanism is also used for explaining time series data (Li et al., 2020;
Vinayavekhin et al., 2018).

##  3\. Preliminaries

###  3.1. Terminology

Let $X=\\{X_{t}\\}_{t\in\mathbb{Z}}$ be a $d$-dimensional time series process
that shows a regularly repeating pattern over time periods of some length $L$.
These repeating patterns are contained in sliding windows
$W_{t}=\\{X_{t+1},\ldots,X_{t+L}\\}$ of $L$ consecutive elements of the time
series. Often, the window size $L$ can be selected based on prior knowledge on
the dataset, which is known to show seasonality over e.g., one day or one
week.

Anomalies in time series data are commonly divided into three categories (Choi
et al., 2021), point anomaly, contextual anomaly, and collective anomaly. In
this work, we consider point anomalies (e.g., abrupt peaks in the data window)
and contextual anomalies (e.g., the appearance of one or more points makes a
temporal context unusual). Formally, we assume that the data points are
generated by a periodically stationary time series process $X$ with
periodicity $L$ (Wang et al., 2015; Ursu and Duchesne, 2009). That is, the
time series consists of regularly repeating patterns of length $L$ which
evolve over time without distributional changes, i.e., we do not consider
concept drifts (Gama et al., 2014).

Let ${(W_{t},y_{t})}_{t\in\mathbb{Z}}$ be the dataset after applying the
sliding window and $y_{i}\in\\{0,1\\}$ is the label of the window $W_{t}$ ($0$
for regular data and $1$ for anomaly). The anomaly detection is conducted on
the window level. A window is considered abnormal if at least one point or
sub-window with multiple points that shows a significantly different behavior
to the window during the regular pattern state. The significance is determined
by comparing the window anomaly score predicted by the model and a user-
defined threshold.

###  3.2. Problem definition

Given the multi-dimensional time series data with applied sliding window, the
target is to train an autoencoder-based end-to-end anomaly detector that

  1. (1)

detect anomaly windows in an unsupervised manner

  2. (2)

deliver representative prototypes of regular data in the latent space

  3. (3)

leverage interpretation of anomalies based on the prototypes of regular data

##  4\. Methodology

In this section, we propose ProtoAD, an LSTM-Autoencoder with an additional
prototype layer, which can be trained end-to-end in an unsupervised manner.

###  4.1. ProtoAD architecture

The architecture of ProtoAD is in line with the existing prototype neural
networks (Ming et al., 2019; Li et al., 2018). We use an LSTM-Autoencoder to
learn time series hidden representations in the latent space and feed the
representations to the prototype layer for similarity-based prototype
comparison. Specifically, we designed the architecture and training procedure
for unsupervised anomaly detection, while only data consisting of regularly
repeating patterns is used for the training and prototype learning. An
overview of the ProtoAD architecture is shown in Figure 1.

EncoderDecoder$W_{t}$$W_{t}^{{}^{\prime}}$$\begin{array}[]{l}Anomaly\\\ \ \ \
score\end{array}$$\begin{array}[]{l}Reconstruction\\\ \ \ \ \ \ \ \ \
error\end{array}$ Figure 1. ProtoAD overview: ProtoAD consists of a encoder
and a decoder, which reconstruct the input data window $W_{t}$ to
$W^{{}^{\prime}}_{t}$. The anomaly score is calculated based on the
reconstruction error. In the latent space, prototypes (stars) are learned as
representative points of the encoded normal windows (white dots). In the test
phase, latent representations of the abnormal data windows (red dots) can also
be explained by comparing with the learned prototypes.

We construct the LSTM-Autoencoder in the fashion of (Malhotra et al., 2016).
More specifically, the $d$ dimensional input window
$W_{t}=\\{X_{t+1}^{t+L}\\}$ is feed into the encoder

$\textbf{f}\colon\mathbb{R}^{L\times d}\to\mathbb{R}^{m}$

The last hidden state of the encoder LSTM unit $h_{i}=\textbf{f}(W_{t})$
($h_{i}\in\mathbb{R}^{m}$) is used as the hidden representation of the input
window in latent space. A same-structured decoder

$\textbf{g}\colon\mathbb{R}^{m}\to\mathbb{R}^{L\times d}$

target at reconstructing the window from the hidden representation
$W_{t}^{{}^{\prime}}=\\{X_{t+1}^{{}^{\prime}t+L}\\}$. The decoder LSTM unit
takes $h_{i}$ as the initial hidden state while takes the real data from
previous timestamp as input. We train the autoencoder to minimize the
reconstruction error of regular windows, i.e., no anomaly data will be used
during training.

The reconstruction error at timestamp $t$ is defined as

$e_{t}=|X_{t}-X_{t}^{{}^{\prime}}|$

The training set is used to estimate a normal distribution
$\mathcal{N}(\mu,\Sigma)$ ($\mathcal{N}(\mu,\sigma)$ for univariate data) of
the reconstruction error for multivariate input data. And the likelihood of a
data point being abnormal is defined by the anomaly score

$a_{t}=\begin{cases}\frac{1}{{\sigma\sqrt{2\pi}}}e^{{{-\left({e_{t}-\mu}\right)^{2}}\mathord{\left/{\vphantom{{-\left({x-\mu}\right)^{2}}{2\sigma^{2}}}}\right.\kern-1.2pt}{2\sigma^{2}}}}&d=1\\\
(e_{t}-\mu)^{T}{\Sigma}^{-1}(e_{t}-\mu)&d>1\end{cases}$

The largest anomaly score is picked up to represented the window anomaly score

$a_{t+1}^{t+L}=\underset{i=1,...,L}{max}(a_{t+i})$

In our work, we do not specify a threshold over the window anomaly scores to
get a binary prediction. Instead, we directly evaluate the AUC score based on
the real-valued anomaly scores. Different existing thresholding techniques can
be applied to get a binary prediction in such a situation (Malhotra et al.,
2016; Su et al., 2019; Hundman et al., 2018).

Based on the anomaly detection model above, we introduce a prototype layer
between the encoder and decoder, which leverages interpretable prototypes of
the regular data during the end-to-end training process. The prototype layer
does not influence the information flow from the encoder to the decoder, i.e.,
the only information the decoder gets from the encoder is the last encoder
hidden state. The prototypes are learned by regularizing the objective
function.

The prototype layer contains $k$ prototypes $p_{i}\in\mathbb{R}^{m}(i=1...k)$
to be learned, where $k$ is a user-defined parameter and $k$ vectors are
randomly initialized within the range $[-1,1]$. As introduced in subsection
4.2, several regularization terms are employed in the objective function to
get the expected prototypes.

In most existing prototype-based models for classification task (Ming et al.,
2019; Kim et al., 2016; Gee et al., 2019), the prototype layer is followed by
some linear layers and a Softmax layer for the production of prediction, which
increases the complexity of the network and requires additional regularization
to enforce interpretability. As an anomaly detection model, our outputs are
derived directly from the autoencoder reconstruction errors. Therefore, we
omit the other common output layers after the prototype layer to simplify the
network structure.

###  4.2. Objective function

The objective in the training phase is to train the autoencoder with regular
windows such that the reconstruction error is minimized and to learn a batch
of prototypes from the regular data. The reconstruction error loss of the
autoencoder is given by

$\mathcal{L}_{e}=\frac{1}{n}\sum_{i=1}^{n}\sum_{l=1}^{L}e_{t+l}$

where $n$ is the number of sliding windows. To ensure the learned prototypes
are informative and diverse enough to each other, we use the diversity loss
designed by (Ming et al., 2019),

$\mathcal{L}_{d}=\sum_{i=1}^{k}\sum_{j=i+1}^{k}max(0,d_{min}-||p_{i}-p_{j}||_{2}^{2})^{2}$

which defines the threshold $d_{min}$ that only apply this penalize to nearby
prototype pairs. Finally, to ensure the prototypes are representative to the
local hidden representations, we define the following representation
regularization term

$\mathcal{L}_{r}=\frac{1}{k}\sum_{j=1}^{k}\underset{i\in[1,n]}{min}||p_{j}-h_{i}||^{2}+\frac{1}{n}\sum_{i=1}^{n}\underset{j\in[1,k]}{min}||h_{i}-p_{j}||^{2}$

The first term ensure that each prototype is close to at least one hidden
representation, which the second term ensures each hidden representation has
one prototype to be represented.

The overall objective function is

$\mathcal{L}=\lambda_{e}\mathcal{L}_{e}+\lambda_{d}\mathcal{L}_{d}+\lambda_{r}\mathcal{L}_{r}$

where $\lambda_{e}$, $\lambda_{d}$ and $\lambda_{r}$ are weighting
hyperparameters.

##  5\. Experiments

In this section, we introduce the experiments on ProtoAD under different
settings. We experiment over different real-world datasets with a variety of
anomaly types. In addition, to evaluate the model performance on specific data
characteristics, we also introduce a synthetic dataset with artificial
anomalies. Finally, we demonstrate the prototypes visually and analyze the
prototype properties w.r.t. a variety of parameter settings.

###  5.1. Experiment setup

####  5.1.1. Datasets

We experiment on one synthetic dataset and four common real-world benchmark
datasets in the time series anomaly detection domain. The dataset properties
are summarised in Table 1.

To understand the anomaly detection process and the learned prototypes, we
introduce a one-dimensional synthetic dataset sampled from a sine wave with
amplitude $1$ and period $100$ timestamps. A random noise $\epsilon\in[0,0.1]$
is added to every timestamp. In addition, we add a random factor
${\alpha}\in[0,1]$ every $100$ timestamps in the test set to simulate point
anomalies. We define a half period of the sine wave as the window length
(i.e., $L=50$), such that the model is supposed to learn the crests and
troughs as two types of prototypes.

The New York City Taxi (Taxi) dataset is a one-dimensional real-world dataset
with a clear periodical feature. It recorded the passenger counts over days in
2014. Extreme passenger counts on public holidays are considered anomalies.
Following (Cui et al., 2016) we aggregate the count numbers into $30$-minute
intervals. We take one day (i.e., $L=48$) as the window length.

SMAP (Soil Moisture Active Passive satellite) and MSL (Mars Science Laboratory
rover) are two multivariate telemetry datasets from NASA (Hundman et al.,
2018). The datasets contain both point and contextual anomalies. Domain
experts labeled the test sets. However, there are also anomaly data in the
training sets. The polluted training set can impact the purity of prototypes.
There is no common repeating pattern in the datasets. We set the window length
$L$ as $100$ for both datasets.

SMD (Server Machine Dataset) (Su et al., 2019) is a multivariate dataset
collected from servers of an Internet company. The dataset is already divided
into two equal-sized training and test sets by the provider and labeled
anomalies by domain experts based on incident reports. We only use the data
from one machine (machine-1-1) in our experiments. We set $L=100$ for SMD.

Table 1. Datasets | Length | Dimensionality | Anomaly rate (%)
---|---|---|---
Synthetic | $20000$ | $1$ | $1.0$
Taxi | $17520$ | $1$ | $4.4$
SMAP | $562800$ | $25$ | $13.1$
MSL | $132046$ | $55$ | $10.7$
SMD | $56958$ | $38$ | $9.5$

####  5.1.2. Evaluation metrics

We adopt the AUC score as the evaluation metric. Considering the essential
requirement of detecting both point and contextual anomalies, we only evaluate
on the window level. A data window is abnormal if it contains one or multiple
abnormal instance(s).

####  5.1.3. Competitors

To the best of our knowledge, this is the first work that engages time series
anomaly detection and prototype learning. The existing prototype learning
networks (Ming et al., 2019; Ni et al., 2021; Gee et al., 2019) commonly work
in a supervised manner, which requires labeled data for the training phase.
Therefore they are not directly relevant to our setting. We mainly compare our
method with the unsupervised anomaly detection approaches. Firstly, we compare
with the EncDecAD (Malhotra et al., 2016), which has a similar setting as
ours, but without the prototype layer. Thereby, we can determine whether the
prototype learning damages the original reconstruction-based anomaly
detection. Furthermore, we compare with one of the state-of-the-art
unsupervised time series anomaly detection OmniAnomaly (Su et al., 2019). We
follow most of the default hyperparameter settings in (Su et al., 2019) but
use the window length same as in our work for the sliding windowing.

####  5.1.4. Experimental details

In all experiments, we set ${\lambda}_{e}=0.025$, ${\lambda}_{d}=0.2$ and
${\lambda}_{r}=0.5$. During training, $25\%$ of the data is used for learning
the parameters $\mu$ and $\Sigma$ ($\sigma$ for univariate data). All models
are trained for $100$ epochs with batch size $20$, learning rate $0.0001$,
dropout rate $0.2$. We use the Adam optimizer (Kingma and Ba, 2014). All
experiments are conducted on a NVIDIA Quadro RTX 6000 24GB GPU. The
experimental results are averaged over three runs.

###  5.2. Evaluation results

####  5.2.1. Anomaly detection performance

Table 2. Anomaly detection performance | EncDecAD | ProtoAD | OmniAnomaly
---|---|---|---
Synthetic | $0.50$ | $0.54$ | $0.95$
Taxi | $0.53$ | $0.63$ | $0.52$
SMAP | $0.41$ | $0.40$ | $0.49$
MSL | $0.73$ | $0.73$ | $0.50$
SMD | $0.95$ | $0.95$ | $0.51$

Firstly we report the AUC score over different models in Table 2. For ProtoAD,
we take the number of prototypes $k=10$. There is no significant difference
between EncDecAD and ProtoAD, which indicates that the additional prototype
layer and corresponding learning process do not directly impact the anomaly
detection performance. ProtoAD even benefits from the prototype learning in
the Synthetic and Taxi datasets. OmniAnomaly shows worse AUC scores in
comparison with the other two models. Different from (Su et al., 2019), where
all possible thresholds over the predicted anomaly scores are traversed, and
only the threshold with the best $F1$ score is reported, the AUC score
reflects more general quality of the anomaly scores over multiple thresholds.

####  5.2.2. Parameter sensitivity

The hidden layer size $m$ and the number of prototypes $k$ are two major
hyperparameters in ProtoAD. In this section, we examine the performance
sensitivity to those two parameters. For each dataset, we try
$m\in[10,50,100,200,400,600,800]$ and $k\in[0,5,10,20,30,50]$, where $k=0$
reduces ProtoAD to EncDecAD. Heatmaps of the AUC scores on real-world datasets
are shown in Figure 2. The results indicate that the model is sensitive to
neither hidden size nor the number of prototypes, as far as the hidden size is
large enough to capture all information in the data windows. However, the
number of prototype $k$ should not be set too large. Otherwise, the model
tends to learn redundant prototypes (see Figure 4).

Figure 2. AUC score under different parameter setting

####  5.2.3. Latent space visualization

We investigate a visualization of the autoencoder hidden space in this section
to understand how time series data windows are embedded and how prototypes of
regular data are learned. We use UMAP (McInnes et al., 2018) to reduce the
high-dimensional latent representations into two dimensions. The result is
visualized in Figure 4. Here, we set $k=5$ for all datasets. The prototypes
shown in the plots are learned during the training phase. The plotted regular
and anomaly points are embedded from the test data.

In the synthetic data, the regular data lie in two regions. Four prototypes
are learned from the trough half (lower left) and one from the crest half
(upper right). Since the anomalies are always generated at the beginning of
the crest half, the anomalies lie nearer to the crest regular data embedding.

In the real-world datasets, especially the SMAP and MSL with polluted training
data, regular and abnormal data do not clearly show separated clusters, though
the learned prototypes represent the major blocks of dense regions showing
regular patterns. Specifically, the prototypes gather at the bottom right
corner for SMD, while no prototype is at the larger upper cluster. A possible
reason is that the high-dimensional server data contain many zero values. The
model can not summarize informative patterns in the training set, and slightly
different regular patterns in the test data are mapped into a different
region.

(a) Taxi regular

(b) Taxi anomaly

(c) Synthetic regular and anomaly

Figure 3. Prototype visualization (blue) with assigned regular (grey) and
anomaly sequences (red)

####  5.2.4. Prototype-based explanation

Finally, we map the prototypes learned in the latent space back to the human-
interpretable input space. Similar to (Chen et al., 2018; Ming et al., 2019),
we map the prototypes back to the input space using the nearest training data
embedding in the latent space to prevent unrealistic produced by the decoder.
Moreover, each neighbor can only be used once, so every prototype is unique.
We visualize the prototypes learned in the one-dimensional datasets Taxi and
Synthetic in Figure 3 with five prototypes ($P1$ to $P5$) for each dataset.

In 3(a) and 3(b), four similar prototypes ($P1$, $P2$, $P4$, $P5$) show a
increasing taxi usage pattern in the morning and turning down at night. $P3$
can be seen as a delayed version of the other four, which is a weekend
pattern. The light lines in the background are the regular (grey) and anomaly
(red) sequences with the smallest distance to the corresponding prototypes in
the latent space. Most of the regular patterns fit the assigned prototypes. A
considerable number of both regular and anomaly sequences have the smallest
latent space distance to $P3$. Some of them visually fit better with other
prototypes, while the distance comparison and prototype assignment do not
directly take place in the input space but in latent space. However, this is
effective for long and high-dimensional sequences. 3(b) depicts the
explanation of anomaly patterns, namely how different are the anomaly
sequences to their nearest prototype. 3(c) shows the three assigned regular
and anomaly sequences (if available) to each prototype. Since the point
anomalies are always generated at the beginning of the crest half, all
anomalies are assigned to crest prototype $P4$.

For the high dimensional datasets, we stay observing the prototypes in the
latent space and leave the searching for informative sub-input space
prototypes as future work. Similarly, we also plan to investigate the
reduction of redundant prototypes (e.g., $P1$, $P2$, $P3$, $P5$ in 3(c)).

Figure 4. The UMAP visualization of the ProtoAD latent space

####  5.2.5. Efficiency comparison

Training the autoencoder with an extra prototype layer does not bring much
training expense. We compare the epoch training time between EncDecAD ($k=0$)
and ProtoAD ($k\in[5,10,20,30,50]$) in Figure 5. As shown in the figure, there
is no significant increase in training time for ProtoAD. In the contrary, due
to the complex model structure, the epoch training time for OmniAnomaly is:
Synthetic $32s$, Taxi $39s$, SMD $225s$, MSL $888s$ and SMAP $3627s$.

Figure 5. Efficiency analysis

##  6\. Conclusion and discussion

In this paper, we explored using prototypes to explain the reconstruction-
based anomaly detection process. Specifically, we integrate the recent end-to-
end prototype learning into the LSTM-Autoencoder. We use the latent space
representations in the autoencoder, which are not directly used in the
conventional reconstruction-based anomaly detection models. In our empirical
evaluation, we figured out that adding a prototype learning step during the
training of the autoencoders does not damage the performance of the
autoencoder. The prototypes contribute to an intuitive understanding of the
regular pattern state.

Although the prototypes learned in the two one-dimensional datasets are
realistic and interpretable for humans, there are still two major problems to
be solved. Firstly, the selection of parameter $k$ is tricky. Pruning
techniques can be applied to reduce the redundancy in the prototypes.
Moreover, the prototypes of high-dimensional data can only be shown as regular
state patterns. However, it is not intuitive enough for humans to directly
figure out a small subset of dimensions of interest. In future work, we plan
to investigate learning prototypes in subspaces.

## References

  * (1)
  * Ahmad et al. (2017) Subutai Ahmad, Alexander Lavin, Scott Purdy, and Zuha Agha. 2017\.  Unsupervised real-time anomaly detection for streaming data.  _Neurocomputing_ 262 (2017), 134–147.
  * Ates et al. (2021) Emre Ates, Burak Aksar, Vitus J Leung, and Ayse K Coskun. 2021\.  Counterfactual Explanations for Multivariate Time Series. In _2021 International Conference on Applied Artificial Intelligence (ICAPAI)_. IEEE, 1–8.
  * Chen et al. (2018) Chaofan Chen, Oscar Li, Chaofan Tao, Alina Jade Barnett, Jonathan Su, and Cynthia Rudin. 2018\.  This looks like that: deep learning for interpretable image recognition.  _arXiv preprint arXiv:1806.10574_ (2018).
  * Choi et al. (2021) Kukjin Choi, Jihun Yi, Changhwa Park, and Sungroh Yoon. 2021\.  Deep Learning for Anomaly Detection in Time-Series Data: Review, Analysis, and Guidelines.  _IEEE Access_ (2021).
  * Cui et al. (2016) Yuwei Cui, Chetan Surpur, Subutai Ahmad, and Jeff Hawkins. 2016\.  A comparative study of HTM and other neural network models for online sequence learning with streaming data. In _2016 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 1530–1538.
  * Gama et al. (2014) João Gama, Indrė Žliobaitė, Albert Bifet, Mykola Pechenizkiy, and Abdelhamid Bouchachia. 2014\.  A survey on concept drift adaptation.  _ACM computing surveys (CSUR)_ 46, 4 (2014), 1–37.
  * Gee et al. (2019) Alan H Gee, Diego Garcia-Olano, Joydeep Ghosh, and David Paydarfar. 2019.  Explaining deep classification of time-series data with learned prototypes. In _CEUR workshop proceedings_ , Vol. 2429. NIH Public Access, 15\.
  * Hase et al. (2019) Peter Hase, Chaofan Chen, Oscar Li, and Cynthia Rudin. 2019\.  Interpretable image recognition with hierarchical prototypes. In _Proceedings of the AAAI Conference on Human Computation and Crowdsourcing_ , Vol. 7. 32–40.
  * He and Zhao (2019) Yangdong He and Jiabao Zhao. 2019.  Temporal convolutional networks for anomaly detection in time series. In _Journal of Physics: Conference Series_ , Vol. 1213. IOP Publishing, 042050\.
  * Hundman et al. (2018) Kyle Hundman, Valentino Constantinou, Christopher Laporte, Ian Colwell, and Tom Soderstrom. 2018\.  Detecting spacecraft anomalies using lstms and nonparametric dynamic thresholding. In _Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining_. 387–395.
  * Kidger et al. (2020) Patrick Kidger, James Morrill, and Terry Lyons. 2020\.  Generalised Interpretable Shapelets for Irregular Time Series.  _arXiv preprint arXiv:2005.13948_ (2020).
  * Kim et al. (2016) Been Kim, Rajiv Khanna, and Oluwasanmi O Koyejo. 2016\.  Examples are not enough, learn to criticize! criticism for interpretability.  _Advances in neural information processing systems_ 29 (2016).
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Li et al. (2020) Guozhong Li, Byron Koon Kau Choi, Jianliang Xu, Sourav S Bhowmick, Kwok-Pan Chun, and Grace LH Wong. 2020.  Efficient shapelet discovery for time series classification.  _IEEE Transactions on Knowledge and Data Engineering_ (2020).
  * Li et al. (2018) Oscar Li, Hao Liu, Chaofan Chen, and Cynthia Rudin. 2018\.  Deep learning for case-based reasoning through prototypes: A neural network that explains its predictions. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 32.
  * Malhotra et al. (2016) Pankaj Malhotra, Anusha Ramakrishnan, Gaurangi Anand, Lovekesh Vig, Puneet Agarwal, and Gautam Shroff. 2016.  LSTM-based encoder-decoder for multi-sensor anomaly detection.  _arXiv preprint arXiv:1607.00148_ (2016).
  * McInnes et al. (2018) Leland McInnes, John Healy, and James Melville. 2018\.  Umap: Uniform manifold approximation and projection for dimension reduction.  _arXiv preprint arXiv:1802.03426_ (2018).
  * Ming et al. (2019) Yao Ming, Panpan Xu, Huamin Qu, and Liu Ren. 2019\.  Interpretable and steerable sequence learning via prototypes. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 903–913.
  * Ni et al. (2021) Jingchao Ni, Zhengzhang Chen, Wei Cheng, Bo Zong, Dongjin Song, Yanchi Liu, Xuchao Zhang, and Haifeng Chen. 2021\.  Interpreting Convolutional Sequence Model by Learning Local Prototypes with Adaptation Regularization. In _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_. 1366–1375.
  * Ren et al. (2019) Hansheng Ren, Bixiong Xu, Yujing Wang, Chao Yi, Congrui Huang, Xiaoyu Kou, Tony Xing, Mao Yang, Jie Tong, and Qi Zhang. 2019.  Time-series anomaly detection service at microsoft. In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_. 3009–3017.
  * Schlegel et al. (2020) Udo Schlegel, Daniela Oelke, Daniel A Keim, and Mennatallah El-Assady. 2020.  An empirical study of explainable AI techniques on deep learning models for time series tasks.  _arXiv preprint arXiv:2012.04344_ (2020).
  * Su et al. (2019) Ya Su, Youjian Zhao, Chenhao Niu, Rong Liu, Wei Sun, and Dan Pei. 2019\.  Robust anomaly detection for multivariate time series through stochastic recurrent neural network. In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_. 2828–2837.
  * Ursu and Duchesne (2009) Eugen Ursu and Pierre Duchesne. 2009.  On modelling and diagnostic checking of vector periodic autoregressive time series models.  _Journal of Time Series Analysis_ 30, 1 (2009), 70–96.
  * Vinayavekhin et al. (2018) Phongtharin Vinayavekhin, Subhajit Chaudhury, Asim Munawar, Don Joven Agravante, Giovanni De Magistris, Daiki Kimura, and Ryuki Tachibana. 2018.  Focusing on what is relevant: Time-series learning and understanding using attention. In _2018 24th International Conference on Pattern Recognition (ICPR)_. IEEE, 2624–2629.
  * Wang et al. (2015) Bing Wang, Jie Sun, and Adilson E Motter. 2015.  Detecting structural breaks in seasonal time series by regularized optimization.  _arXiv preprint arXiv:1505.04305_ (2015).
  * Xu et al. (2018) Haowen Xu, Wenxiao Chen, Nengwen Zhao, Zeyan Li, Jiahao Bu, Zhihan Li, Ying Liu, Youjian Zhao, Dan Pei, Yang Feng, et al. 2018\.  Unsupervised anomaly detection via variational auto-encoder for seasonal kpis in web applications. In _Proceedings of the 2018 World Wide Web Conference_. 187–196.
  * Zhai et al. (2016) Shuangfei Zhai, Yu Cheng, Weining Lu, and Zhongfei Zhang. 2016\.  Deep structured energy based models for anomaly detection. In _International Conference on Machine Learning_. PMLR, 1100–1109.
  * Zhang et al. (2021) Zaixi Zhang, Qi Liu, Hao Wang, Chengqiang Lu, and Cheekong Lee. 2021.  ProtGNN: Towards Self-Explaining Graph Neural Networks.  _arXiv preprint arXiv:2112.00911_ (2021).
  * Zong et al. (2018) Bo Zong, Qi Song, Martin Renqiang Min, Wei Cheng, Cristian Lumezanu, Daeki Cho, and Haifeng Chen. 2018.  Deep autoencoding gaussian mixture model for unsupervised anomaly detection. In _International Conference on Learning Representations_.


# Chirped pulse control over the melting of superconductors

Maria Recasens  ICFO - Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, Av. Carl Friedrich Gauss 3, 08860
Castelldefels (Barcelona), Spain<EMAIL_ADDRESS>Valentin Kasper  ICFO
- Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain
<EMAIL_ADDRESS>Maciej Lewenstein  ICFO - Institut de Ciencies
Fotoniques, The Barcelona Institute of Science and Technology, Av. Carl
Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain  ICREA, Pg. Lluis
Companys 23, 08010 Barcelona, Spain  Allan S. Johnson  IMDEA Nanoscience,
Calle Faraday 9, 28023 Madrid, Spain<EMAIL_ADDRESS>
###### Abstract

Strong field terahertz pulses are increasingly used to excite and control
quantum materials at the ultrafast timescale. They have found widespread
application by enabling direct addressing of the superconducting gap or
Josephson resonances and are essential in Higgs spectroscopy. Large non-linear
optical signals can be induced by the strong coupling of the THz and
superconducting degrees of freedom. However, far less attention has been paid
to the strong bi-directional coupling between field and material this implies.
Here, we use the framework of the time-dependent Ginzburg-Landau equations to
study the full field and material evolution of a superconductor driven by
strong field terahertz pulses. We find that at high field strengths, the
backreaction of the superconductor induces large changes to the driving pulse,
which in turn leads to a runaway melting of the superconducting condensate.
This results in a surprisingly large sensitivity to the initial driving pulse
chirp, enabling these purely dynamical changes to result in order of magnitude
different levels of melting. We also find large-scale spectral shifting of the
driving pulse to occur in just a few hundred nanometers of propagation through
a superconductor. We attribute these effects to an inverse plasma redshift, in
which the driving field breaks Cooper pairs and decreases the free-electron
mobility, analogous to reducing the density of a plasma.

_Introduction._ Terahertz (THz) pulses have recently emerged as a powerful
tool for exciting and controlling quantum materials. The low photon energies
allow the THz radiation to couple the collective excitations efficiently in
various systems. Furthermore, with the advent of high-field pulses with peak
electric field strength of hundreds of kilovolts to megavolts per centimeter,
it has become possible to drive quantum materials far out of equilibrium Yang
_et al._ (2023). THz driving has been particularly successfully applied to
superconductors with various low-energy collective modes. Examples include the
use of THz driving to observe hidden stripe phases, parametrically amplify
superconducting fluctuations, and coherently drive Higgs modes, leading to new
forms of nonlinear spectroscopy, emergent phases, and control of
superconductivity Cavalleri (2018); Budden  _et al._ (2021); Shimano and Tsuji
(2020).

In many of these cases, it has been noted that the nonlinear response of the
superconductor to the THz radiation is exceptionally large, with cross
sections far greater than can be found in conventional materials Matsunaga
_et al._ (2014); Chu  _et al._ (2020). This is true both when the THz is
resonant to particular collective modes, but also works off-resonance when the
primary coupling is directly to the superconducting condensate Cea  _et al._
(2016); Robson  _et al._ (2017a). Especially, in the latter case, the light
field can drive large time-dependent supercurrents. Numerous studies and
models have been made to explain these large non-linear effects with focus
upon the accurate treatment of the superconducting response Cea  _et al._
(2016); Shimano and Tsuji (2020); Demsar (2020). However, at high-fields it
can be important to model the entire system of light and superconductor Robson
_et al._ (2017b); Robson and Biancalana (2018a), i.e., the back-reaction of
the superconductor on the THz field. This back-reaction can drastically modify
the response of the system, and the ability to drive and control
superconductors with light selectively Cea  _et al._ (2016). As larger field
strength THz pulses become available Rader  _et al._ (2022), this backreaction
can be expected to be even more significant.

Here we present work modeling the impact of the superconducting backreaction
on an intense THz driving field and how this strongly modifies the overall
system dynamics. We directly model and solve a one-dimensional system to study
the response of a superconducting condensate and the propagation of the THz
driving laser field. Particularly, we treat the superconducting response via
the phenomenological time-dependent Ginzburg Landau (TDGL) equations Robson
_et al._ (2017a) and the field propagation by the full Maxwell equations. We
focus on the effect of the initial chirp of the pulse in the interaction
Warren  _et al._ (1993). Our results demonstrate a large dependence of the
system on the initial pulse chirp at high THz field strengths, with realistic
chirp variations leading to differences in the superconducting melting up to
one order of magnitude. This dependence is attributed to a spectral shift in
the driving pulse occurring within a range of a few hundred nanometers. We
associate this spectral shift with the breaking of Cooper pairs caused by the
pulse, giving rise to an inverse plasma redshift effect Bloembergen (1973).
The substantial nonlinear response observed in the superconductor has the
potential to improve control and diagnostics of superconductivity.

##  I Model

A wide variety of models have been used previously to describe the dynamics of
superconductors out of equilibrium. These include microscopic models like the
time-dependent Gutzwiller variational approach or time-dependent Dynamical
Mean Field Theory, see Giannetti  _et al._ (2016) for a review. On the other
hand phenomenological models, such as the two-temperature model or the
Rothwarf-Taylor model Coslovich  _et al._ (2011), can return qualitative
understanding about the melting process but generally do not directly treat
the coupling between laser pulses and superconductors.

For our simulations, we use the time-dependent Ginzburg Landau theory coupled
to the Maxwell equations. The TDGL equation describes the dynamics of the
superconductor order parameter $\psi$ and the Maxwell equations the
propagation of the EM field. The TDGL equation was derived by Schmidt Schmid
(1966), and later Gor’kov and Eliashberg showed that it could also be derived
from the BCS theory in the case of gapless superconductors Gor’kov (1959). The
order parameter $\psi$ can be directly related to the Cooper pair density
$n_{c}$ via $n_{c}=|\psi|^{2}$ through BSC theory. Critically, the TDGL makes
concrete and quantitative predictions for a wide variety of dynamical
behaviors, such as the dynamics of vortices in type II superconductors Dorsey
(1992); Machida and Kaburaki (1993) or of fluctuations in the conductivity
above $T_{c}$ Schmidt (1968). TDGL has previously been applied to studying THz
driving of superconductors, where it was found that the TDGL could semi-
quantitatively describe the giant Kerr response in superconductors Robson  _et
al._ (2017a); Robson and Biancalana (2018a). We briefly note that an
alternative approach to treat the complete field+superconductor system would
be models based on the Boguliobov-de Gennes Hamiltonian Vaswani  _et al._
(2021); Luo  _et al._ (2023); Mootz  _et al._ (2020).

The TDGL equation Gropp (1996) is

$\displaystyle\frac{\hbar^{2}}{2m^{*}D}\frac{\partial\psi}{\partial
t}=\frac{\partial f}{\partial\psi^{*}}\,,$ (1)

where $D$ is the diffusion coefficient, $m^{*}$ is the effective mass of the
Cooper pair and $f$ is the free energy, where we already implicitly assume the
temporal gauge $\phi=0$. In the Ginzburg-Landau equation for
superconductivity, the free energy is given by

$\displaystyle
f=f_{n}+a|\psi|^{2}+\frac{b}{2}|\psi|^{4}+\frac{1}{2m^{*}}|(-i\hbar\nabla-e^{*}\textbf{A})\psi|^{2}\,,$
(2)

where $a=a_{0}\left(T-T_{c}\right),a_{0}$ and $b$ are constants of the
material, $e^{*}$ and $m^{*}$ are the effective charge and mass of the Cooper
pairs, and $\bm{A}$ is the vector potential that represents the external light
field. The parameters $a_{0}$ and $b$ can be determined through their
relationships to the experimentally measurable London penetration depth
$\lambda$ and coherence length $\xi$, which are defined as;

$\displaystyle\xi(T)$ $\displaystyle=\sqrt{\hbar^{2}/(2m^{*}|a(T)|)},$ (3a)
$\displaystyle\lambda(T)$
$\displaystyle=\sqrt{m^{*}/(4\mu_{0}e^{2}|\psi_{0}|^{2})}\,,$ (3b)
$\displaystyle|\psi_{0}|^{2}$ $\displaystyle=-a_{0}/b,$ (3c)
$\displaystyle\tau$ $\displaystyle=\xi^{2}/D$ (3d)

where $|\psi_{0}|^{2}$ is the solution of the TDGL equation in the case of no
spatial dependence and no electromagnetic field applied. The Ginzburg–Landau
parameter, $\kappa=\lambda/\xi$, can also be obtained through these
quantities.

By inserting the free energy of Eq. (2) in Eq. (1), we obtain

$\displaystyle\frac{\hbar^{2}}{2m^{*}D}\frac{\partial\psi}{\partial
t}=\frac{1}{2m^{*}}(-i\hbar\partial_{j}-2eA_{j})^{2}\psi-(a+b|\psi|^{2})\psi.$
(4)

Further, we rewrite the TDGL using the abbreviation introduced in Eq. (3)
leading to

$\displaystyle\tau\partial_{t}\psi=\xi^{2}\left(-i\partial_{j}-\frac{2e}{\hbar}A_{j}\right)^{2}\psi+\left(1-\frac{|\psi|^{2}}{|\psi_{0}|^{2}}\right)\psi.$
(5)

To describe the superconductor interaction with the EM pulse, it is also
necessary to include in the model the supercurrent equation, i.e., the
dissipation-less current,

$\displaystyle
J_{s}=\frac{e^{*}}{m^{*}}\text{Re}[\psi^{*}(-i\hbar\nabla-e^{*}\textbf{A})\psi]\,,$
(6)

where $\psi^{*}$ is the complex conjugate of $\psi$ and Re denotes the real
part Tinkham (2015).

Furthermore, to fully describe the coupled system, we also need to include the
response of normal-state materials. Here, we use the Drude model to consider
the remaining free electrons of the superconductor not bound into Cooper
pairs. The evolution of the electric field is, as a result, described by

$\displaystyle\partial_{t}E_{x}(y,t)$
$\displaystyle=\epsilon^{-1}_{0}\partial_{y}H_{z}(y,t)+\epsilon^{-1}_{0}[J_{d}(t)+J_{s}(y,t)],$
(7a) $\displaystyle\partial_{t}H_{z}(y,t)$
$\displaystyle=\mu^{-1}_{0}\partial_{y}E_{x}(y,t)\,,$ (7b)

where $\epsilon_{0}$ is the permittivity of free space and $J_{d}$ is the
Drude current, whose dynamics are described by

$\displaystyle\partial_{t}J_{d}(t)$ $\displaystyle=-\gamma
J_{d}(t)+\omega_{p}^{2}\epsilon_{0}E_{x}(y,t)\,,$ (8)

where $\gamma$ is a damping factor and $\omega_{p}$ the plasma frequency.
Thus, the entire system of equations we solve to describe the material
response and field evolution is given by Eq. (5) - (8). The TDGL theory is
valid within certain limits, specifically when the frequency of the time
variation of $\psi$ is greater than the gap frequency as discussed in Abrahams
and Tsuneto (1966); Tinkham (2015).

###  I.1 Geometry and setup

To investigate the back reaction of a superconductor on a THz field, we
simulate the case of a high-field THz pulse incident on bulk material.
Specifically, we focus on the effect of the initial chirp of the pulse, as
introducing a frequency chirp modifies the pulse without affecting the total
energy of the pulse. In particular, we model the scenario of a THz pulse
hitting a superconducting thin film as illustrated in Fig. 1, allowing us to
describe the system in one spatial dimension. For the explained geometry of
the simulation, the order parameter $\psi$ can hence be treated as only
dependent on $y$.

Figure 1: Simulation setup. We consider a thin, infinite two-dimensional
superconducting plate (green) with a chirped probing pulse (orange).

###  I.2 Pulse and Material

The dynamics of a superconductor driven by a THz pulse will depend on the
parameters describing the geometry, the incident pulse, and the material
properties. Concerning the incoming pulse, we mainly focus on the pulse chirp.
Tuning the chirp does not change the total input energy, which gives a unique
look into the dynamic response of the system.

Regarding the driving THz field, we start from a simple Gaussian pulse with an
electric field given by

$\displaystyle
E(t)=E_{0}e^{-(t-t_{0})^{2}/t_{g}^{2}}\cos{[\omega_{0}(t-t_{0})]}\,,$ (9)

where $E_{0}$ is the electric field strength, $t_{g}$ determines the width of
the envelope in the time domain, $\omega_{0}$ is the angular photon frequency,
and $t_{0}$ is the center of the pulse in the time domain.

We include the chirp by multiplying the Fourier transform of the pulse given
in Eq. (9) with

$\displaystyle e^{-i[a(2\pi f-\omega_{0})+b(2\pi f-\omega_{0})^{2}+d(2\pi
f-\omega_{0})^{3}+...]}$ (10)

and then Fourier transforming back to the time domain, where $a$ quantifies
the scalar phase shift and $b$ the linear chirp. This pure phase addition
stretches the pulse in the time domain but does not change the total energy or
the spectral content. Instead, it modifies the arrival time of the different
spectral components. To obtain a chirp of order $n-1$, we add a term
proportional to $(2\pi f-\omega_{0})^{n}$ in the exponent of Eq. (10).

Concerning the initial pulse, we explore a range of field strengths spanning
one order of magnitude, from
$0.9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ to
$9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ for the unchirped pulses.
These fields are high but technically feasible Rader  _et al._ (2022).
Further, we note using thinner samples allows us to imitate lower field
strengths in certain circumstances. We consider photon frequencies ranging
from the THz to the mid-infrared domain with $t_{g}=$132\text{\,}\mathrm{fs}$$
and $\omega_{0}=$4.8\text{\,}\mathrm{THz}$$. As for the frequency chirp, we
introduce a linear chirp, keeping the parameters $a=c=0$.

Motivated by experiment, we consider a $310\text{\,}\mathrm{nm}$ thick film of
lead, with a critical temperature of $T_{c}$ of $7.2\text{\,}\mathrm{K}$, see
Wesche (2017). The material choice determines the TDGL parameters. For a BCS
superconductors superconductor, $m^{*}$ and $e^{*}$ take straightforward
values related to the Cooper-pair structure, $m^{*}=2m_{e}$ and $e^{*}=2e$
Cabrera and Peskin (1989); Gor’kov (1959). For lead, the remaining parameters
take the following values $\xi=$82\text{\,}\mathrm{nm}$$, $\kappa=0.48$ Wesche
(2017) and
$\rho=$3.84\text{\times}{10}^{-10}\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}\mathrm{m}$$.
The simulations are initialized with the system in equilibrium at
$T=$0\text{\,}\mathrm{K}$$.

##  II Results

We calculate the propagation of the pulse and the dynamics of the material by
numerically integrating the system of equations Eq. (5) to Eq. (8). At the
boundaries of the material, we enforce $\psi=0$ and generate the propagating
pulse using a source far away from the material. To avoid reflections of the
EM pulse at the boundaries of the simulation, we added perfectly matched
layers far away from the source and the material. Finally, we discretize space
using finite differences and solve the time evolution with a 4th-order Runge-
Kutta solver. In the following subsections, we analyze our simulation results,
focusing on the influence of the chirp. We analyze the effect of the sign of
the chirp, the response of the superconducting order parameter concerning the
chirp, and the response of the electromagnetic field.

(a)

(b)

(c)

(d)

(e)

Figure 2: Pulse sequence and spatio-temporal distribution of order parameter
Interaction between the superconductor and the electromagnetic pulse for a
peak field strength of $8\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ for
the unchirped pulse. The red lines represent the boundaries of the
superconductor. In subfigure (a) we show the electric field of EM pulses with
different linear frequency chirps. Further, in subfigure (b) and (c) the order
parameter as a function of time and space normalized by the maximum value of
the initial thermal state, $|\psi|^{2}/|\psi_{\text{max}}|^{2}$, where
$|\psi_{\text{max}}|^{2}=$3.8517\text{\,}{\mathrm{nm}}^{-3}$$. In subfigures
(d) and (e) we show the time and space evolution of the electric field
amplitude squared, $E^{2}$ in
$\text{\,}{\mathrm{V}}^{2}\text{\,}{\mathrm{m}}^{-2}$. The explicit value of
the chirp is $b=-$20\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$$ in (b)
and (d), and $b=$20\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$$ in
subfigure (c) and (e).

(a)

(b)

Figure 3: Tuning order parameter with linear chirp. The maximum value of the
order parameter after the pulse has passed the superconductor, at
$t=$12\text{\,}\mathrm{ps}$$ for two different initial field strengths. The
value of the order parameter is normalized by the maximum value of the initial
state, $|\psi_{\text{max}}|^{2}/\text{max}(|\psi_{\text{therm}}|^{2})$ with
$\text{max}(|\psi_{\text{therm}}|^{2})=$3.8517\text{\,}{\mathrm{nm}}^{-3}$$.
Green vertical lines denote where the chirp is of magnitude
$5\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$. In subfigures (a) and (b),
we use a peak field strength of
$9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ and
$0.9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ respectively for $b=0$.

###  II.1 Positive chirp - negative chirp

We simulate the interaction between the superconductor and the pulse for two
different pulses, one with positive and one with negative chirp. Fig. 2a
illustrates the electric field of pulses with chirps of
$20\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$,
$0\text{\,}{\mathrm{fs}}^{2}$ and
$-20\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$. The high-frequency
components arrive first for positive chirps, while for negative chirps we
observe the opposite.

In Fig. 2, we study the spatial and temporal results of the described
simulations. Specifically, Fig. 2b and Fig. 2d show the results for the
negative linear chirp, whereas Fig. 2c and Fig. 2e depict the results for the
positive linear chirp. By comparing the left and right columns, we note that
the two scenarios show significant qualitative differences in the evolution of
the superconductor. Specifically, we observe a more significant suppression of
the order parameter for negative chirps, and the suppression occurs at earlier
times for the positive chirp case. In addition, the figures Fig. 2d and Fig.
2e show the evolution of the corresponding EM pulse inside the superconductor.
For the pulse with positive chirp, see Fig. 2e, we observe that low and high
frequencies are similarly suppressed inside the superconductor, whereas for
the negative chirp pulse, see Fig. 2d, high frequencies are more suppressed
than lower frequencies. These differences indicate tuning the response of the
material by the sign of the chirp.

###  II.2 Suppression of superconductivity using the chirp

To gain more insight into the chirp dependence of the dynamics, we
systematically vary the chirp and quantify the extinction of the order
parameter by studying the maximum value of $\psi_{\text{max}}$ inside the
superconductor at the end of the simulation, specifically $\psi_{\text{max}}$
at $12\text{\,}\mathrm{ps}$. The results of this study are given in Fig. 3
where we show $\psi_{\text{max}}$ as a function of the linear chirp of the
source pulse for two different field strengths one order of magnitude apart.

In both plots, for chirps with an absolute value larger than
$5\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$, we can see that negatively
chirped pulses suppress the order parameter more than positively chirped
pulses. Thus, negative chirped pulses more effectively destroy
superconductivity than positive chirped ones. As we discuss in the following
section, this marked difference arises from the dependence of the backreaction
of the superconductor on the sign of the chirp. This result holds for an
extensive range of field strengths, although the magnitude of the effect
changes significantly. In addition to this effect, there are clear
oscillations in $\psi_{\text{max}}$ with chirp, along with exceptionally high
suppression of the amplitude for chirps lower than
$5\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$ in the high field case.

###  II.3 Optical back-reaction

To understand the origin of these changes in suppression, we first examine the
transmission of the pulse as a function of chirp, as shown in Fig. 4, focusing
on the high field case where the effects are more dramatic. We immediately
note that the changes observed in Fig. 3 do not map directly to changes in
transmission; the asymmetry between positive and negative chirps is small, and
there are no oscillations. However, the absorption increases dramatically for
small chirp values at the highest pulse intensity. This behavior does not
follow a simple power law of the intensity. It thus suggests that for low-
chirp, high-intensity pulses, additional non-linear melting processes of the
condensate become significant. We have verified this assertion by performing
additional scans under zero chirp conditions in which we keep the total energy
in the pulse constant but change $t_{g}$, leading to a decrease in peak field,
which indeed confirms a nonlinear melting process occurs at high field
strengths, see Supplementary Fig. 6.

Figure 4: Transmission coefficient: For different chirps, we depict the total
transmitted energy after the material normalized by the total energy of the
source pulse for a peak field strength of
$9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$ at $b=0$. We observe a
clear minimal transmission coefficient close to zero chirp, and observe
asymmetric tuning of of the transmission for positive and negative chirps.

As changes in total absorption cannot explain the strongly asymmetric melting
process nor the oscillations, we next examine the backreaction of the
superconductor on the spectrum of the pulse. Fig. 5 displays the spectra of
the source pulse and the output pulse spectrum for a range of chirps. Solid
lines correspond to the spectra of pulses after going through the material,
and the dashed black line is the spectra of the source pulses. We see a clear
spectra shift towards lower frequencies for negative chirps lower than
$-5\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$. On the other hand, for
positive chirps higher than
$5\text{\times}{10}^{3}\text{\,}{\mathrm{fs}}^{2}$, we observe no shift and an
asymmetric narrowing of the spectra. This offers a partial explanation for the
behavior of the superconductor; low-frequencies couple more strongly to the
condensate Robson  _et al._ (2017a) and lead to a non-thermal suppression of
the superconducting condensate, and hence this giant redshift enhances the
non-thermal melting process. A redshift of this type has previously been
predicted for THz-driven superconductors, but only under the zero-chirp
condition where we see a significantly smaller redshift Robson and Biancalana
(2018a).

Figure 5: Spectrum of the laser pulse after interaction. Spectra of the pulse
after the material for a field strength of
$9\text{\,}\mathrm{MV}\text{\,}{\mathrm{cm}}^{-1}$. Red solid lines correspond
to the spectra of the pulses with negative chirps, purple solid lines to the
spectra of the Gaussian pulse without chirp, and blue solid lines to the
spectra of the positive chirped pulses. The dashed black line corresponds to
the input spectra of the source, which is the same for all chirps.

Finally, we discuss the origin of the oscillations in $\psi_{\text{max}}$ in
Fig. 3. The magnitude of these oscillations increases with increasing $E_{0}$
magnitude. However, their period remains constant, and they do not appear in
either the transmission or spectral shift. Changing the pulse central
frequency $\omega_{0}$ does not affect the oscillation period either.
Modifying the material thickness alters the period of the oscillations,
suggesting that they can be attributed to the finite-size effects of the
superconducting film.

###  II.4 Asymmetric spectral broadening and melting

Previous work Robson and Biancalana (2018b) related the melting of
superconductivity to a giant asymmetric third-order non-linear optical
response. Here, we offer an intuitive explanation of this effect and how this
leads to strong chirp sensitivity of the melting process. A driven material
with a time-varying refractive index modulates the phase of the driving field,
leading to new frequencies, i.e., self-phase modulation Boyd (2008).
Typically, the refractive index follows a simple rule
$n(t)=n_{0}+3\chi^{(3)}I(t)/(4n_{0}^{2}\epsilon_{0}c)$, where $n_{0}$ is the
static refractive index, $\chi^{(3)}$ is the third-order non-linear
susceptibility and $I(t)$ is the intensity of the driving pulse. Because
$n(t)$ then follows the temporal behavior of $I(t)$, there is a region where
$n(t)$ increases, then decreases. For a positive $\chi^{(3)}$, which describes
most materials, rising intensities lead to an increase in refractive index and
a redshift, while falling intensities lead to a blueshift Boyd (2008).

In a superconductor, however, the leading contributions to $\chi^{(3)}$ come
from the melting of the superconducting condensate Robson  _et al._ (2017a).
In this case, $n(t)$ does not follow $I(t)$ as the pulse can only melt the
superconductor, and the evolution of $n(t)$ is monotonic. This is conceptually
very close to the idea of the plasma-blue shift: in a gas ultrafast laser
ionization generates free carriers, leading to the formation of a plasma and a
blue shift in the driving laser frequency because of the associated time-
dependent decrease in the refractive index Bloembergen (1973); Giulietti and
Giulietti (2015). However, in a superconductor, the inverse effect is at play;
by breaking Cooper pairs (effectively melting a high-mobility plasma), we
dynamically reduce the effective free-carrier response, leading to a redshift.
Robson _et al._ Robson  _et al._ (2017b), derived an approximate expression
for $\chi^{(3)}$ in superconductors

$\displaystyle\chi^{(3)}=-\frac{128k_{b}^{2}(T-T_{c})^{2}e^{4}}{\epsilon_{0}\beta\omega_{0}^{6}\pi^{2}\hbar^{2}m_{e}^{2}},$
(11)

where $k_{b}$ is the Boltzmann constant. As can be seen, low frequencies
(small $\omega_{0}$) dramatically increase the Kerr-like response of the
superconductor. So, the red shift in the driving pulse, in turn, more
effectively melts the superconductor and dynamically reinforces the melting
process in a positive feedback loop, as alluded to above. This inverse plasma
redshift intuitively explains the giant redshift Robson and Biancalana (2018a)
and why the redshifted pulse is linked to more substantial melting.

Next, we consider the role of the chirp explicitly. That the sign of chirp
dramatically affects the degree of broadening induced by nonlinear processes
has long been known for conventional self-phase modulation Oberthaler and
Höpfel (1993), with the exact behavior depending on whether the newly
generated frequencies are additive or out-of-phase with existing spectral
components of the pulse. However, in the case of a melting superconductor,
there is an additional history effect, as the magnitude of $\chi^{(3)}$ is
directly related to the number of cooper pairs that remain in the sample, as
represented by the $(T-T_{c})^{2}$ factor in Eq. (11). A positively chirped
pulse thus red-shifts the low-frequency components effectively, leading to
very strong superconducting suppression in the leading edge of the pulse but
comparatively little red-shifting in the trailing edge. As such, only the
leading edge effectively melts the superconductor, as can be seen examining
Fig. 2c, leading to the most substantial absorption at these frequencies. In
contrast, the trailing high-frequency edge undergoes only a minor redshift.

In contrast, for a negative chirped pulse, the redshift is strongest for the
high-frequency components, and the overall melting process is more distributed
over the entire pulse duration, as seen in Fig. 2b. This leads to a more even
shift to low frequencies across the whole bandwidth and allows the total
spectral bandwidth to contribute to the suppression of the superconductivity.
Thus, this history effect explains why both the red-shift of negatively
chirped pulses and the narrowing of positively chirped pulses saturate so
effectively in Fig. 5, and why the effects depend intimately on the particular
nature of the melting of the superconducting condensate.

##  III Conclusion

With increasingly high-field terahertz pulses being used to manipulate
superconductors, the back-reaction of the superconductor upon the driving
pulse will become more and more significant. Here, we have presented a study
of the entire field and material evolution of a BCS-type superconductor driven
by high-field strength THz pulses through a clear and transparent approach
based on the TDGL theory. We find the dynamics of the superconductor to be
strongly dependent on the initial chirp of the driving pulse, leading to
variations up to one order of magnitude of the melting of the superconductor
and large-scale spectral shifting of the driving pulse within a few hundred
nanometers of propagation through a superconductor. We attribute these
dynamics to an inverse plasma redshift response of the condensate. Similarly
to the plasma-blueshift phenomena known in plasma physics, here the
destruction of Cooper pairs caused by the driving field effectively decreases
the free-carrier density/mobility, causing a time-dependent change of the
refractive index, which red-shifts the driving field, leading to stronger non-
thermal suppression of the condensate. In the limit of strong melting, this
shift depends upon the history of the pulse, so the input chirp can modulate
the overall degree of melting significantly.

This history effect and giant red-shift are rather unique features of
superconductors. Thus, the uncovered high dependence of the state of the
superconductor on purely dynamical variables, such as the chirp, indicates
that pulse shaping could be an interesting tool for the control and diagnosis
of superconductivity. It also suggests that superconductors could form an
interesting class of non-linear optic materials in the THz regime. Our
simulations show that macroscopic models of the superconductor-light
interaction can reveal unexpected effects in the dynamics of such systems and
may need to be explicitly treated to understand experiments performed at high
field strengths.

##  IV Acknowledgements.

This work was funded by the Spanish AIE (projects PID2022-137817NA-I00 and
EUR2022-134052).

ICFO group acknowledges support from: ERC AdG NOQIA; MCIN/AEI
(PGC2018-0910.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan
National FIDEUA PID2019-106901GB-I00, Plan National STAMEENA
PID2022-139099NB-I00 project funded by MCIN/AEI/10.13039/501100011033 and by
the “European Union NextGenerationEU/PRTR" (PRTR-C17.I1), FPI); QUANTERA MAQS
PCI2019-111828-2); QUANTERA DYNAMITE PCI2022-132919 (QuantERA II Programme co-
funded by European Union’s Horizon 2020 program under Grant Agreement No
101017733), Ministry for Digital Transformation and of Civil Service of the
Spanish Government through the QUANTUM ENIA project call - Quantum Spain
project, and by the European Union through the Recovery, Transformation and
Resilience Plan - NextGenerationEU within the framework of the Digital Spain
2026 Agenda; Fundació Cellex; Fundació Mir-Puig; Generalitat de Catalunya
(European Social Fund FEDER and CERCA program, AGAUR Grant No. 2021 SGR 01452,
QuantumCAT U16-011424, co-funded by ERDF Operational Program of Catalonia
2014-2020); Barcelona Supercomputing Center MareNostrum (FI-2023-1-0013); EU
Quantum Flagship (PASQuanS2.1, 101113690, funded by the European Union. Views
and opinions expressed are, however, those of the author(s) only and do not
necessarily reflect those of the European Union or the European Commission.
Neither the European Union nor the granting authority can be held responsible
for them); EU Horizon 2020 FET-OPEN OPTOlogic (Grant No 899794); EU Horizon
Europe Program (Grant Agreement 101080086 — NeQST), results incorporated in
this standard have received funding from the European Innovation Council and
SMEs Executive Agency under the European Union’s Horizon Europe programme).,
ICFO Internal “QuantumGaudi” project; European Union’s Horizon 2020 program
under the Marie Sklodowska-Curie grant agreement No 847648; “La Caixa” Junior
Leaders fellowships, La Caixa” Foundation (ID 100010434): CF/BQ/PR23/11980043.
Views and opinions expressed are, however, those of the author(s) only and do
not necessarily reflect those of the European Union, European Commission,
European Climate, Infrastructure and Environment Executive Agency (CINEA), or
any other granting authority. Neither the European Union nor any granting
authority can be held responsible for them.

ASJ acknowledges the support of the Ramón y Cajal Program (Grant
RYC2021-032392-I). IMDEA Nanociencia acknowledges support from the “Severo
Ochoa” Programme for Centers of Excellence in R&D (MICIN, CEX2020-001039-S).

## References

  * Yang  _et al._ (2023) C.-J. Yang, J. Li, M. Fiebig, and S. Pal, Nature Reviews Materials , 1 (2023).
  * Cavalleri (2018) A. Cavalleri, Contemporary Physics 59, 31 (2018).
  * Budden  _et al._ (2021) M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, F. Schlawin, D. Jaksch, and A. Cavalleri, Nature Physics 17, 611 (2021).
  * Shimano and Tsuji (2020) R. Shimano and N. Tsuji, Annual Review of Condensed Matter Physics 11, 103 (2020).
  * Matsunaga  _et al._ (2014) R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, K. Makise, Y. Uzawa, H. Terai, Z. Wang, H. Aoki, and R. Shimano, Science 345, 1145 (2014).
  * Chu  _et al._ (2020) H. Chu, M.-J. Kim, K. Katsumi, S. Kovalev, R. D. Dawson, L. Schwarz, N. Yoshikawa, G. Kim, D. Putzky, Z. Z. Li, H. Raffy, S. Germanskiy, J.-C. Deinert, N. Awari, I. Ilyakov, B. Green, M. Chen, M. Bawatna, G. Cristiani, G. Logvenov, Y. Gallais, A. V. Boris, B. Keimer, A. P. Schnyder, D. Manske, M. Gensch, Z. Wang, R. Shimano, and S. Kaiser, Nature Communications 11, 1793 (2020).
  * Cea  _et al._ (2016) T. Cea, C. Castellani, and L. Benfatto, Physical Review B 93, 180507 (2016).
  * Robson  _et al._ (2017a) C. W. Robson, K. A. Fraser, and F. Biancalana, Physical Review B 95, 214504 (2017a).
  * Demsar (2020) J. Demsar, Journal of Low Temperature Physics 201, 676 (2020).
  * Robson  _et al._ (2017b) C. W. Robson, K. A. Fraser, and F. Biancalana, Physical Review B 95, 214504 (2017b).
  * Robson and Biancalana (2018a) C. W. Robson and F. Biancalana, Materials Research Express 5, 046001 (2018a).
  * Rader  _et al._ (2022) C. Rader, Z. B. Zaccardi, S.-H. E. Ho, K. G. Harrell, P. K. Petersen, M. F. Nielson, H. Stephan, N. K. Green, D. J. H. Ludlow, M. J. Lutz, S. J. Smith, D. J. Michaelis, and J. A. Johnson, ACS Photonics 9, 3720 (2022), https://doi.org/10.1021/acsphotonics.2c01336 .
  * Warren  _et al._ (1993) W. S. Warren, H. Rabitz, and M. Dahleh, Science 259, 1581 (1993).
  * Bloembergen (1973) N. Bloembergen, Optics Communications 8, 285 (1973).
  * Giannetti  _et al._ (2016) C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D. Mihailovic, Advances in Physics 65, 58 (2016).
  * Coslovich  _et al._ (2011) G. Coslovich, C. Giannetti, F. Cilento, S. Dal Conte, G. Ferrini, P. Galinetto, M. Greven, H. Eisaki, M. Raichle, R. Liang, _et al._ , Physical Review B 83, 064519 (2011).
  * Schmid (1966) A. Schmid, Physik der Kondensierten Materie 5, 302 (1966).
  * Gor’kov (1959) L. P. Gor’kov, Sov. Phys. - JETP (Engl. Transl.); (United States) 9:6 (1959).
  * Dorsey (1992) A. T. Dorsey, Physical Review B 46, 8376 (1992).
  * Machida and Kaburaki (1993) M. Machida and H. Kaburaki, Physical Review Letters 71, 3206 (1993).
  * Schmidt (1968) H. Schmidt, Zeitschrift für Physik A Hadrons and nuclei 216, 336 (1968).
  * Vaswani  _et al._ (2021) C. Vaswani, J. H. Kang, M. Mootz, L. Luo, X. Yang, C. Sundahl, D. Cheng, C. Huang, R. H. J. Kim, Z. Liu, Y. G. Collantes, E. E. Hellstrom, I. E. Perakis, C. B. Eom, and J. Wang, Nature Communications 12, 258 (2021).
  * Luo  _et al._ (2023) L. Luo, M. Mootz, J.-H. Kang, C. Huang, K. Eom, J. Lee, C. Vaswani, Y. Collantes, E. Hellstrom, I. E. Perakis, _et al._ , Nature Physics 19, 201 (2023).
  * Mootz  _et al._ (2020) M. Mootz, J. Wang, and I. E. Perakis, Phys. Rev. B 102, 054517 (2020).
  * Gropp (1996) e. Gropp, Journal of Computational Physics 123, 254 (1996).
  * Tinkham (2015) M. Tinkham,  _Introduction to superconductivity_ , 2nd ed., Dover books on physics (Dover Publ, Mineola, NY, 2015).
  * Abrahams and Tsuneto (1966) E. Abrahams and T. Tsuneto, Physical Review 152, 416 (1966).
  * Wesche (2017) R. Wesche,  _High Temperature Superconductors - Handbook of Electronic and Photonic Materials_ (2017) publication Title: Springer Handbook of Electronic and Photonic Materials.
  * Cabrera and Peskin (1989) B. Cabrera and M. E. Peskin, Physical Review B 39, 6425 (1989).
  * Gor’kov (1959) L. P. Gor’kov, J. Exptl. Theoret. Phys. (U.S.S.R.) 36, 1364 (1959).
  * Robson and Biancalana (2018b) C. W. Robson and F. Biancalana, Materials Research Express 5, 046001 (2018b).
  * Boyd (2008) R. W. Boyd,  _Nonlinear Optics_ (Academic Press, 2008).
  * Giulietti and Giulietti (2015) A. Giulietti and D. Giulietti, Journal of Plasma Physics 81, 495810608 (2015).
  * Oberthaler and Höpfel (1993) M. Oberthaler and R. A. Höpfel, Applied Physics Letters 63, 1017 (1993), https://pubs.aip.org/aip/apl/article-pdf/63/8/1017/18497373/1017_1_online.pdf .

##  V Supplemental material

For low-chirp values, the pulse has high-intensity peaks, and our results thus
suggest this results in additional non-linear melting processes of the
condensate becoming significant. To verify this assertion, we performed
additional scans under zero chirp conditions. We keep the total energy in the
pulse constant but change $t_{g}$, decreasing the peak field. As we observe in
Fig. 6, lower values of $t_{g}$, with higher peak intensities, significantly
suppress the superconducting state, confirming the hypothesis.

Figure 6: Maximum value of the order parameter as a function of the pulse
duration Maximum value of the order parameter after the pulse has gone through
the superconductor, at $t=$12\text{\,}\mathrm{ps}$$, for incident unchirped
pulses with different values of $t_{g}$. The value is normalized by the
maximum value of the initial state,
$|\psi|^{2}_{\text{max}}/\text{max}(|\psi_{\text{therm}}|^{2})$, where
$\text{max}(|\psi_{\text{therm}}|^{2})=$3.8517\text{\,}{\mathrm{nm}}^{-3}$$.