S-JEPA: TOWARDS SEAMLESS CROSS-DATASET TRANSFER THROUGH DYNAMIC SPATIAL
ATTENTION
Pierre Guetschel1, Thomas Moreau2, Michael Tangermann1
1Donders Institute for Brain, Cognition and Behaviour, Radboud University,
Nijmegen, Netherlands
2Université Paris-Saclay, Inria, CEA, Palaiseau 91120, France
E-mail<EMAIL_ADDRESS>
ABSTRACT: Motivated by the challenge of seamless cross-dataset transfer in EEG
signal processing, this article presents an exploratory study on the use of
Joint Embedding Predictive Architectures (JEPAs). In recent years, self-
supervised learning has emerged as a promising approach for transfer learning
in various domains. However, its application to EEG signals remains largely
unexplored. In this article, we introduce Signal-JEPA for representing EEG
recordings which includes a novel domain-specific spatial block masking
strategy and three novel architectures for downstream classification. The
study is conducted on a 54 subjects dataset and the downstream performance of
the models is evaluated on three different BCI paradigms: motor imagery, ERP
and SSVEP. Our study provides preliminary evidence for the potential of JEPAs
in EEG signal encoding. Notably, our results highlight the importance of
spatial filtering for accurate downstream classification and reveal an
influence of the length of the pre-training examples but not of the mask size
on the downstream performance.

Figure 1:  S-JEPA training procedure. The framework takes as input EEG
recordings with $C$ channels and $T$ time samples, and binary masks of length
$L$. First, the Local encoder independently transforms $t$ windows from each
channel into $C\times t=L$ embedding vectors, called tokens, of dimensionality
$d$. Then, the tokens are marked according to their originating channel and
temporal position, and are flattened into a sequence of length $L$.
Subsequently, only the unmasked tokens are passed to the Contextual encoder,
while the full tokens sequence is given to the Contextual target encoder to
generate training targets. Finally, the Predictor attempts to reconstruct the
masked tokens and its predictions are compared with the targets using an L1
loss. During the optimisation, the parameters of the Contextual target encoder
are not trained via gradient backpropagation but follow those of the
Contextual encoder by Exponential Moving Average (EMA). Figure inspired from
[1].

INTRODUCTION

Electroencephalography (EEG) allows capturing neural activity directly from
the scalp, offering a high temporal resolution signal to investigate brain
functions. The interpretation of EEG signals with machine learning methods
opens the door to build brain-computer interfaces (BCIs). Despite the
potential of BCIs, their practical application is hindered by the intensive
requirement for calibration data, which is both time-consuming and demanding
for participants. Recently, transfer learning has been explored to mitigate
the constraints imposed by data scarcity and calibration demands. Self-
supervised learning (SSL) presents itself as a pivotal strategy to tackle
transfer learning, enabling models to learn rich representations from
unlabeled data that can be used to efficiently solve downstream tasks. A
particular approach to SSL consists in strategically masking parts of the
input data and training a model to predict these masked elements. A key
advantage of masking-based SSL methods is that they can be applied to
virtually any type of data.

Masking Strategies proved a determining factor for the success of SSL methods
in various domains, including image [2, 3, 4], speech [5, 2], text [2], and
video [1] processing. For instance, random masking strategies, where the
regions to be masked are sparsely selected, typically deliver inferior results
compared to block masking strategies, where larger continuous regions are
masked, requiring the model to gain a deeper understanding of the data
distribution [2].

Masked Autoencoders (MAEs) stand as the quintessential entry point to masking-
based SSL [3]. These models aim to reconstruct the masked sections of the
input directly, and the training objective is computed by comparing the
original input sections to the reconstructed ones. Unfortunately,
reconstructing the input and comparing elements in the original space is not
without challenges. When the original input space has a high dimensionality,
the reconstruction can be computationally expensive, and necessitate the use
of domain-specific constraints to produce valid signals. The original signals
can also be noisy, increasing the difficulty of encoding meaningful parts of
the signal. Moreover, the reconstruction’s difficulty can vary significantly
across different areas of the input data: with images, reconstructing a
monochromatic, non-structured sky is less difficult than reconstructing a
structured object like a hand. Such structural disparities are one of the
reasons, why conventional metrics like mean-square error often fail to assess
the reconstruction quality. These challenges question the scalability and
adaptability of MAEs as a universal SSL approach.

Joint-Embedding Predictive Architectures (JEPAs) offer a promising alternative
to address the limitations associated with direct reconstruction. JEPA-like
methods avoid reconstructing the input in its original space and focus on
predicting latent representations, or embeddings, of the data [4]. This
approach confers two major benefits: first, it is computationally efficient,
especially with high-dimensional input spaces as the embeddings can reduce the
dimensionality; second, the metric’s selection in the embedding space is less
critical, as the embeddings are learned adaptively to the chosen metric.
However, as the "ideal" embedding vectors are unknown, the reconstruction
objective is undefined a priori. This challenge is addressed by constructing
target embedding vectors during the training by using a bootstrapping
procedure which will be further explained in the S-JEPA section. The potential
of JEPA-like frameworks has been highlighted by their promising results with
images [2, 4], speech [5, 2], text [2], and videos [1].

Applications to the EEG Domain of masking-based SSL techniques have started to
emerge. Valesco and colleagues used a masked autoencoder approach, with random
masking over the spatial and temporal dimensions [6]. Chien and colleagues
experimented with MAEs with block masking over the temporal dimension [7].
Kostas et al. and Foumani et al. also used block masking over the temporal
dimension but coupled it with JEPA-like training strategies [8, 9].

Motivation. Despite the advancements in applying SSL to EEG data, the
exploration of block masking strategies over EEG channels remains uncharted
territory. Such an approach holds the potential for developing robust channel
attention mechanisms, and thus could facilitate dynamic spatial filtering.
This capability could prove instrumental in adapting to recordings with
varying channel sets, thereby facilitating cross-dataset transfer learning or
tackling corrupted channels. This paper seeks to bridge this gap by
investigating the implications of a channel-based block masking strategy on
SSL efficacy in EEG signal processing.

The main application of models trained with SSL is the following fine-tuning
on the actual task of interest, the so-called downstream task. While most
research within the EEG domain has focused on fine-tuning for sleep stage
classification tasks, their application to a BCI context remains largely
untapped, with only two studies by Kostas et al. and Velasco et al. exploring
the impact of SSL on BCI tasks [8, 6]. As BCI systems suffer from data
scarcity there is a consistent goal to minimize the amount of calibration data
required before each online session, specifically as the calibration phase
requires sustained attention from the participant and thus is tiring. This
highlights a significant opportunity to explore the effectiveness of SSL
models across various BCI paradigms, including but not limited to motor
imagery protocols, thereby contributing to the broader understanding and
application of SSL in enhancing BCI performance.

Furthermore, the application of pre-trained SSL models for solving downstream
tasks often involves the addition of a linear layer atop the embedding
dimension. This practice, however, may not be optimal in high-dimensional
embedding spaces. Through a comparative analysis of six different strategies
for leveraging pre-trained architectures in downstream tasks, this work aims
to support our understanding and application of SSL in EEG data processing.

Research Questions and Plan. This manuscript is an explorative study
investigating what approaches should be adopted for training SSL algorithms on
EEG signals, and what domain-specific considerations warrant attention. We
approach this through three research questions: 1) What constitutes the most
efficacious masking strategy for SSL when applied to EEG data? 2) How does the
temporal length of examples used influence the SSL training process? 3) What
fine-tuning strategies lead to the best downstream performance?

To answer these, we first propose a novel masking strategy as part of the
Signal-based Joint-Embedding Predictive Architecture (S-JEPA) framework in the
S-JEPA section. Then, we introduce fine-tuning strategies tailored for S-JEPA
in the Spatial Block Masking section. The datasets used for pre-training and
for downstream evaluation are presented in the Downstream classification
architectures section. Finally, the Downstream classification architectures
and Downstream classification scores sections will report and critically
discuss the outcomes and implications of the experiments we conducted.

S-JEPA FRAMEWORK

The Signal-based Joint-Embedding Predictive Architecture (S-JEPA) framework is
illustrated in Figure 1. It is used to pre-train models. Its architecture is
inspired by BENDR and MAEEG, introduced in the pioneering studies by
respectively Kostas et. al. [8] and Chien et. al. [7]. A key modification is
introduced in the design of the local encoder to enable our novel masking
strategy. This section details the architecture’s components, the spatial
masking strategy, and the training process.

Local encoder. The local encoder is implemented as a convolutional neural
network (CNN) with five layers, each formed of a convolution with a Gaussian
error linear unit (GELU) non-linearity. Contrary to the encoders in BENDR and
MAEEG, which accept multi-channel input windows, our encoder processes windows
from a singular channel. Each window is encoded into a 64-dimensional
embedding vector, hereafter referred to as a token. The windows are 1.19 s
long, with a stride of 1.0 s. The first convolutional kernel covers 0.25 s
while the following ones simply combine feature vectors in pairs, i.e., both
the kernel temporal lengths and strides are 2.

Contextual encoder. This encoder consists of a transformer architecture with
eight layers, as introduced by Vaswani and colleagues [10]. It processes the
unordered sequence of tokens generated by the local encoder, necessitating the
addition of position-encoding information to localize them temporally and
spatially. The temporal positioning of each token is defined using a cosine
encoding [10] over the first 34 dimensions, whereas spatial positioning is
achieved through trainable embeddings for each channel, initialized using
cosine encoding based on their three-dimensional coordinates. The contextual
encoder receives tokens that only contain local information; its role is to
establish relationships between them.

(a) Cz as center.

(b) C6 as center.

(c) Oz as center.

Figure 2:  Visualisation of the spatial block masking strategy for three
example mask centres (red electrodes). The dark to light green spheres
represent masks of diameters 40 %, 60 % and 80 % of the head size, as used in
our experiments. Assuming a top-down view upon the scalp, the depth of the
electrodes is denoted by their intensity (black: close, grey: distant). For a
given mask, all electrodes within the corresponding sphere are hidden from the
contextual encoder and must be predicted by the predictor.

Spatial Block Masking strategy. Unique to our methodology is the independent
encoding of each channel by the local encoder. While it avoids the learning of
spatial filters at early stages, it paves the way for innovative spatial
dimension masking strategies. Literature in both, image processing [2] and EEG
signal analysis [9] suggests that block masking yields superior results
compared to random masking. Motivated by these findings, our novel approach
extends block masking to the spatial domain of EEG channels.

Given the irregular distribution of EEG channels, the concept of a contiguous
block of tokens does not translate directly from its application in images or
temporal signals, where the pixels and time samples are regularly spaced. Thus
our approach masks all channels within a predetermined radius of a randomly
chosen central channel, as illustrated in Figure 2. In this work, we compare
three mask sizes with diameters approx. 40 %, 60 % and 80 % of the head size.
This strategy inherently introduces variability in the number of masked
tokens.

S-JEPA pre-training. With the operational principles of both the local and
contextual encoders established, we introduce two ancillary components
exclusively utilized during the training phase. The first, termed the
Contextual target encoder, is a non-trainable duplicate of the contextual
encoder which serves to generate the training targets. Its parameters are
updated via exponential moving average (EMA). The second component, the
Predictor, is a transformer decoder architecture with four layers as
delineated by Vaswani et al. [10]. The comprehensive training methodology is
illustrated in Figure 1.

The models are trained until no improvement of the validation loss is observed
for 10 epochs (i.e., complete pass through the entire dataset), at which point
the best model is saved for subsequent fine-tuning.

DOWNSTREAM EVALUATION

Upon completion of the network’s pre-training through the SSL task, which
holds no intrinsic value beyond training purposes, we proceed to evaluate its
efficacy on practical downstream classification tasks which, in our case, are
BCI tasks. This step is crucial for determining the real-world applicability
of the pre-trained model.

Downstream classification architectures. Using pre-trained models for
downstream classification tasks necessitates altering their architecture to
allow predicting class probabilities. The most widely adopted modification,
linear probing, consists of adding a linear classification layer directly
above the embedding space [4]. However, the individual tokenization of each
channel in our approach leads to a high-dimensional latent space, close to the
dimensionality of the raw input examples, which would make linear probing
inefficient. In response to this challenge, we enrich the architecture with
two layers instead which are Spatial aggregation and Fully-Connected, as
explained in Figure 3.

The integration of these layers is explored in three distinct configurations.
(3(a)) The Contextual downstream architecture places both layers after the
contextual encoder as depicted in 3(a). (3(b)) The Post-local downstream
architecture discards the pre-trained contextual encoder and adds the novel
layers atop the local encoder as shown in 3(b). (3(c)) The Pre-local
downstream architecture also discards the pre-trained contextual encoder but
then places the spatial averaging layer before the local encoder as
illustrated in 3(c). This third alternative allows the network to perform a
spatial EEG filtering step, as commonly present in BCI architectures.

(a) Contextual downstream architecture.

(b) Post-local downstream architecture.

(c) Pre-local downstream architecture.

Figure 3:  Downstream classification architectures. In each of the three
alternative alterations of the pre-trained networks, two new layers are added.
1) Spatial aggregation is a convolutional layer that realizes weighted
combinations of the elements in the channels dimension into $V\ll C$ "virtual"
channels. 2) Fully-connected is a linear layer that predicts $c$ class
probabilities.

Fine-tuning. We examine two distinct fine-tuning strategies for these
downstream architectures. The first method, indicated by the prefix new-,
involves exclusively training the newly introduced layers, keeping the pre-
trained components frozen. The alternative strategy, denoted with the prefix
full-, consists in fine-tuning the entire network. This second strategy starts
with a warm-up phase of 10 epochs where only the newly added layers undergo
training, preventing the deterioration of the pre-trained layers’ performance
due to irrelevant feedback [11], before the previously existing layers are
included into the training.

For both strategies, the model is fine-tuned until no improvement of the
validation loss is observed for 50 epochs, at which point the best model is
restored for testing. Additionally, it should be noted that the temporal
length of examples used during the fine-tuning phase is determined by the
requirements of the downstream task, which is independent of the length of
examples utilized during the SSL pre-training phase.

DATASETS

For the exploratory investigation in this work, we used the dataset introduced
by Lee and colleagues [12], subsequently referred to as the lee2019 dataset.
It contains EEG recordings from 54 subjects, each undergoing three distinct
BCI paradigms: steady-state visual evoked potentials (SSVEP) with four
classes, visual event-related potentials (ERP), and left vs. right hand motor
imagery (MI). The recordings are from 62 EEG channels, the spatial
distribution of which is detailed in Figure 2. The dataset was loaded from the
MOABB framework [13].

We used the first 40 subjects to pre-train any model. The subsequent 7
subjects were used for the validation during this pre-training phase. The
remaining 7 subjects were reserved for the downstream performance evaluation.

Pre-training data. SSL methods do not necessitate labels. As such, the
training and validation examples are slices of the continuous recordings taken
at a fixed interval of 16.9 seconds. This study compared three example lengths
approximately distributed on a logarithmic scale, namely 1, 4, and 16 seconds.
All allocated subjects and all paradigms were used collectively during the
pre-training phase which yielded a total of 36,576 training examples and 6,528
validation examples. No artefact rejection method was applied.

Downstream evaluation data. For each subject allocated for downstream
evaluation and each paradigm, the fine-tuning performance of the different
pre-trained models is assessed using a 5-folds within-subject stratified
cross-validation procedure. The examples used are 4.19 seconds long for MI and
SSVEP, and 1.19 seconds long for ERP. These lengths are chosen to respect the
durations defined in the original dataset [12] while maintaining compatibility
with our tokenization process.

RESULTS

Experimental details. For pre-training, our setup compares all combinations of
signal durations (1s, 4s, and 16s) and mask sizes (40 %, 60 %, and 80 % of
head size), along with a no-pre-training baseline.

For downstream performance evaluation, we assess each pre-trained model using
the three downstream architectures—contextual, post-local, and pre-local—and
the two fine-tuning approaches—full and new. The no-pre-training is only
assessed with the full fine-tuning approach. A given pairing of pre-training
and fine-tuning configurations is referred to as a pipeline.

Figure 4:  Pre-training curves of the different configurations tested. The
solid and dashed lines indicate the loss on the training and validation sets.
While the validation loss was tested once per epoch only, the training loss
was logged after every optimisation step. The train loss on individual
optimisation steps is visible in the background, corresponding epoch-wise
averages are outlined in white. A star marks the lowest validation loss per
curve, the early stopping time point and consequently the checkpoint from
which any fine-tuning started.

Pre-training dynamics. The pre-training phase, see Figure 4, provides insights
into the training process under the S-JEPA framework. It reveals that training
curves under the 16s condition are significantly smoother than the other
configurations. Additionally, the early stopping mechanism concludes the
training of the 16s-80% configuration prematurely at 12 epochs due to an early
trough in the loss curve. Similarly, the 16s-40% setup also encounters an
early trough yet manages to recover and complete a longer training. The
longest training durations are observed in the 16s-40% and 16s-60%
configurations, enduring for 74 and 58 epochs, respectively, translating to
approximately 12 and 10 hours of training.

Figure 5:  Global downstream classification ranking of all the combinations of
pre-training configurations and fine-tuning schemes. Each of the three test
datasets has 7 subjects and 5 folds per subject, which makes a total of 105
folds. In the legend, the combinations are ordered according to their average
rank over all folds. The vertical span of a coloured "pixel" in the plot
represents the number of folds in which this configuration has obtained the
rank indicated by the x-axis.

Pipelines ranking on downstream performances. The comprehensive ranking of all
tested pipelines is detailed in Figure 5, offering a comparative overview of
the performance of the different combinations over the experimental protocols.
Notably, the top-performing pipelines in the downstream tasks are the 16s-60%
and 16s-40% models, especially when paired with any pre-local fine-tuning
strategy. They notably occupy the first rank in two-thirds of the cases. In
contrast, lower-performing pipelines feature new-contextual or full-contextual
fine-tuning, particularly when combined with 16s and no-pre-training
configurations.

Paradigm-wise downstream performances. on the individual paradigms is reported
in Figure 6. For the subsequent analysis, we discern several key observations:
1) The pipelines obtaining the best score on the ERP, SSVEP and MI tasks are
respectively 16s-40%-full-pre-local with a 97% AUC, 16s-60%-new-pre-local with
a 94% accuracy, and 16s-40%-new-pre-local with a 65% accuracy. 2) Pipelines
combining 16s or no-pre-training with contextual architectures frequently
result at chance level on average. 3) Pipelines combining 1s or 4s with
contextual architectures also perform at chance level on the SSVEP paradigm
but above chance level on the MI and ERP ones. 4) Most pipelines manage to
achieve respectable scores on the ERP dataset. 5) Only a select subset of
pipelines excel on the SSVEP dataset. 6) The MI paradigm scores exhibit
notable variability, as indicated by the considerable standard deviation.

Figure 6:  Downstream classification scores of all the pipelines on the three
test datasets. The height of the coloured bars corresponds to the average
classification score over all the test subjects and cross-validation folds,
while the thin black bars correspond to their standard deviation.

DISCUSSION

In the light of observations made in the Downstream classification
architectures section, we aim to answer the research questions posed, draw
conclusions, and provide guidelines for future research.

Examples’ length largely influences downstream performance. In particular, the
16s pipelines consistently show the best performance, highlighting the
advantage of a longer context during pre-training. Conversely, when
considering solely the results using the contextual downstream architecture,
the 16s pipelines often perform at chance level. On the other hand, the 1s and
4s pipelines yield better results, which may be attributed to their better
alignment of the attention mechanism’s training with the short signal windows
of the downstream tasks.

Mask radius’ impact on downstream performance uncertain, as our results do not
demonstrate a clear trend between the masks compared. It is possible that the
range of masks we compared is not optimal or that other factors are
influencing the results. Future work should consider comparing our spatial
masking strategy with temporal masking to better understand the relative
strengths and weaknesses of these approaches.

The best fine-tuning strategy implements spatial filtering The pre-local
architecture emerges as the best for downstream classification. This
architecture’s approach to linearly combine channels before computing features
enables effective spatial filtering, thereby enhancing the signal-to-noise
ratio. This finding underscores the critical role of spatial filtering in
boosting model performance by leveraging the inherent spatial properties of
EEG data.

State-of-the-art comparison. According to Chevallier and colleagues’
manuscript in preparation [14], which benchmarks numerous decoding algorithms
across all the datasets available in MOABB [13], the current state-of-the-art
(SOTA) performances for within-session classification on the Lee2019 dataset
are: ERP at $98.41\pm 2.03\%$, SSVEP at $89.44\pm 13.84\%$, and MI at
$84.74\pm 13.19\%$. Notably, all pipelines establishing the SOTA utilize
Riemannian geometry. Our approach matches the SOTA for ERP, enhances it for
SSVEP, but falls short on MI. A critical difference in our evaluation
methodologies should be noted: our downstream evaluation only focuses on the
last 7 subjects, unlike their analysis on all 54 subjects. Specifically, Lee
and colleagues have identified 6 out of these 7 subjects as hard to classify
on the MI task [12]. We believe this exceptionally high rate of challenging
subjects might explain our low MI performance.

Choice of Dataset. The need for large datasets is paramount when training
transformers, potentially explaining the underperformance of the contextual
downstream strategy. Although this exploratory study on the Lee2019 dataset
provides valuable insights, future research should pivot towards larger
datasets to fully harness the capabilities of contextual architectures.

Conclusion. This exploratory work introduces a novel masking strategy and
three fine-tuning approaches, positioning our method competitively within the
realm of BCI tasks. We achieve SOTA performance on two out of three evaluated
downstream tasks. Our findings suggest that long pre-training windows favor
the local features encoder, while short windows benefit the contextual
encoder. Therefore, future research should aim at successfully training both
the local and contextual encoders. However, no influence of the mask radius on
the downstream performance was found. Finally, the best downstream
architecture includes a spatial filtering step and discards the contextual
encoder.

ACKNOWLEDGEMENTS

This work is in part supported by the Donders Center for Cognition (DCC) and
is part of the project Dutch Brain Interface Initiative (DBI2) with project
number 024.005.022 of the research programme Gravitation which is (partly)
financed by the Dutch Research Council (NWO).

REFERENCES

## References

  * [1] Adrien Bardes et al. “Revisiting Feature Prediction for Learning Visual Representations from Video”
  * [2] Alexei Baevski, Arun Babu, Wei-Ning Hsu and Michael Auli “Efficient Self-supervised Learning with Contextualized Target Representations for Vision, Speech and Language” arXiv, 2022 DOI: 10.48550/arXiv.2212.07525
  * [3] Kaiming He et al. “Masked Autoencoders Are Scalable Vision Learners” In _2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ New Orleans, LA, USA: IEEE, 2022, pp. 15979–15988 DOI: 10.1109/CVPR52688.2022.01553
  * [4] Mahmoud Assran et al. “Self-Supervised Learning from Images with a Joint-Embedding Predictive Architecture” arXiv, 2023 arXiv:2301.08243 [cs, eess]
  * [5] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed and Michael Auli “Wav2vec 2.0: A Framework for Self-Supervised Learning of Speech Representations” In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 12449–12460
  * [6] S Pérez Velasco et al. “Evaluación del Impacto del Aprendizaje Auto-Supervisado en la Precisión de Interfaces Cerebro-Ordenador basadas en Imaginación Motora”, 2023
  * [7] Hsiang-Yun Sherry Chien, Hanlin Goh, Christopher M. Sandino and Joseph Y. Cheng “MAEEG: Masked Auto-encoder for EEG Representation Learning” arXiv, 2022 DOI: 10.48550/arXiv.2211.02625
  * [8] Demetres Kostas, Stéphane Aroca-Ouellette and Frank Rudzicz “BENDR: Using Transformers and a Contrastive Self-Supervised Learning Task to Learn From Massive Amounts of EEG Data” In _Frontiers in Human Neuroscience_ 15, 2021, pp. 653659 DOI: 10.3389/fnhum.2021.653659
  * [9] Navid Mohammadi Foumani et al. “EEG2Rep: Enhancing Self-supervised EEG Representation Through Informative Masked Inputs” arXiv, 2024 arXiv:2402.17772 [cs, eess]
  * [10] Ashish Vaswani et al. “Attention Is All You Need” In _Advances in Neural Information Processing Systems_ 30 Curran Associates, Inc., 2017
  * [11] Ananya Kumar et al. “Fine-Tuning Can Distort Pretrained Features and Underperform Out-of-Distribution” In _ICLR_ arXiv, 2022 DOI: 10.48550/arXiv.2202.10054
  * [12] Min-Ho Lee et al. “EEG Dataset and OpenBMI Toolbox for Three BCI Paradigms: An Investigation into BCI Illiteracy” In _GigaScience_ 8.5, 2019, pp. giz002 DOI: 10.1093/gigascience/giz002
  * [13] Bruno Aristimunha et al. “Mother of All BCI Benchmarks”, Zenodo, 2023 DOI: 10.5281/ZENODO.10034223
  * [14] Sylvain Chevallier et al. “The Largest EEG-based BCI Reproducibility Study for Open Science: The MOABB Benchmark”, 2024


# Photo-induced phase transition and associated time scales in the excitonic
insulator Ta2NiSe5

Tanusree Saha  Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>Laboratory of Quantum Optics, University of Nova
Gorica, 5001 Nova Gorica, Slovenia.  Denis Golež  Center for Computational
Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Department of Theoretical Physics, Institute Jožef Stefan, Jamova 39, SI-1001
Ljubljana, Slovenia.  Department of Complex Matter, Institute Jožef Stefan,
Jamova 39, SI-1001 Ljubljana, Slovenia.  Giovanni De Ninno  Laboratory of
Quantum Optics, University of Nova Gorica, 5001 Nova Gorica, Slovenia.
Elettra-Sincrotrone Trieste, Area Science Park, 34149 Trieste, Italy.  Jernej
Mravlje  Department of Theoretical Physics, Institute Jožef Stefan, Jamova 39,
SI-1001 Ljubljana, Slovenia.  Yuta Murakami  Department of Physics, Tokyo
Institute of Technology, Meguro, Tokyo 152-8551, Japan  Barbara Ressel
Laboratory of Quantum Optics, University of Nova Gorica, 5001 Nova Gorica,
Slovenia.  Matija Stupar  Laboratory of Quantum Optics, University of Nova
Gorica, 5001 Nova Gorica, Slovenia.  Primož Rebernik Ribič ∗ Laboratory of
Quantum Optics, University of Nova Gorica, 5001 Nova Gorica, Slovenia.
Elettra-Sincrotrone Trieste, Area Science Park, 34149 Trieste, Italy.

###### Abstract

We investigate the non-equilibrium electronic structure and characteristic
time scales in a candidate excitonic insulator, Ta2NiSe5, using time- and
angle-resolved photoemission spectroscopy with a temporal resolution of 50 fs.
Following a strong photoexcitation, the band gap closes transiently within 100
fs, i.e., on a time scale faster than the typical lattice vibrational period.
Furthermore, we find that the characteristic time associated with the rise of
the photoemission intensity above the Fermi energy decreases with increasing
excitation strength, while the relaxation time of the electron population
towards equilibrium shows an opposite behaviour. We argue that these
experimental observations can be consistently explained by an excitonic origin
of the band gap in the material. The excitonic picture is supported by
microscopic calculations based on the non-equilibrium Green’s function
formalism for an interacting two-band system. We interpret the speedup of the
rise time with fluence in terms of an enhanced scattering probability between
photo-excited electrons and excitons, leading to an initially faster decay of
the order parameter. We show that the inclusion of electron-phonon coupling at
a semi-classical level changes only the quantitative aspects of the proposed
dynamics, while the qualitative features remain the same. The experimental
observations and microscopic calculations allow us to develop a simple and
intuitive phenomenological model that captures the main dynamics after
photoexcitation in Ta2NiSe5.

##  I Introduction

Excitonic insulators were originally proposed as materials where a spontaneous
condensation of excitons takes place in the vicinity of the semiconductor-
semimetal transitionMott (1961); Knox (1963). Such a condensation is
associated to the violation of separate charge conservation in the conduction
and valence bands (internal $U(1)$ gauge symmetry). While this phase has been
unambiguously identified for heterostructures in the quantized Hall regime
Butov  _et al._ (2002); Eisenstein and MacDonald (2004); Wang  _et al._
(2019), only a few single-material systems have been suggested to exhibit a
ground-state excitonic insulator phase. Among them, Ta2NiSe5Wakisaka  _et al._
(2009, 2012) and 1T-TiSe2Cercellier  _et al._ (2007) are the most prominent
ones. However, in real materials, there are always effects that may break the
$U(1)$ symmetryMazza  _et al._ (2020). In addition, the presence of electron-
phonon coupling Kaneko  _et al._ (2013); Murakami  _et al._ (2017); Subedi
(2020) can lead to a structurally ordered phase with properties similar to
that of an excitonic insulator. Indeed, Ta2NiSe5 shows a structural phase
transition from an orthorhombic to a monoclinic lattice at $328$K Kaneko  _et
al._ (2013); Di Salvo  _et al._ (1986). This raises an important long-standing
question: what is the dominant origin of the order in Ta2NiSe5 ? A number of
experiments have been performed, employing different methods such as Raman
spectroscopy Kim  _et al._ (2020a, b); Volkov  _et al._ (2020), optical
techniques Larkin  _et al._ (2017, 2018), and angle-resolved photoemission
spectroscopy (ARPES) Seki  _et al._ (2014); Watson  _et al._ (2020), to answer
this question. However, the results were interpreted in either the excitonic
or phononic picture and no consensus has been reached so far.

Tracking the system as it evolves through a nonequilibirum phase transition
triggered by a short optical pulse can provide further clues on the origin of
the symmetry breaking in Ta2NiSe5. However, also in this case, the situation
is controversial, as several groups performing time-resolved ARPES (trARPES)
experiments have reported rather contradicting results Mor  _et al._ (2017);
Okazaki  _et al._ (2018a); Baldini  _et al._ (2020); Tang  _et al._ (2020);
Suzuki  _et al._ (2020). Some studies Mor  _et al._ (2017); Baldini  _et al._
(2020) found a short-lived increase of the gap, which was interpreted as a
transient enhancement of the excitonic order Mor  _et al._ (2017), stimulating
further theoretical studies Murakami  _et al._ (2017); Tanaka  _et al._
(2018); Tanabe  _et al._ (2018); Fujiuchi  _et al._ (2019). On the other hand,
several recent experimental works have observed a decrease of the gap Mor  _et
al._ (2017); Okazaki  _et al._ (2018a); Tang  _et al._ (2020), albeit marginal
in recent high-resolution data Baldini  _et al._ (2020), that is naturally
expected as the excited electrons can scatter with condensed excitons Golež
_et al._ (2016). Differences have also been observed in the measured
characteristic time scales of the photoemission signal after photoexcitaion.
Ref. Okazaki  _et al._ , 2018a reports a characteristic time scale of around
100 fs, which is faster than the typical lattice oscillation in Ta2NiSe5. On
the other hand, the authors of Ref. Baldini  _et al._ , 2020 argue that this
time scale is comparable to the lattice distortion period (characteristic time
above 250 fs) and consider this as an evidence of a lattice-driven transition.
More generally, the behaviour of the time-scales vs. the excitation fluence is
also controversial. In particular, this concerns the apparent disagreement
between the dynamical slowing down, observed in both experimental Zong  _et
al._ (2019a); Tomeljak  _et al._ (2009) as well as some theoretical studies
focused on symmetry-broken states Tsuji  _et al._ (2013, 2014); Babadi  _et
al._ (2015); Golež  _et al._ (2016); Babadi  _et al._ (2017); Picano and
Eckstein (2020); Dolgirev  _et al._ (2020), and the measured rise time of the
photoemission intensity above the Fermi energy, which becomes shorter for
increasing photo-excitation strengths Okazaki  _et al._ (2018a). Summing up
the above contrasting results, the question on the dominant driving force of
the order in Ta2NiSe5 still remains open.

In this paper, we study the non-equilibrium electron dynamics of Ta2NiSe5 in a
high photoexcitation regime, which was not explored in previous works. Using
trARPES with a temporal resolution of 50 fs, we show that a strong
photoexcitation induces a phase transition from an insulating to a semi-
metallic phase and that the melting of the band gap occurs on a time scale
faster than the period of a phonon oscillation. Next, we show that the
dependence of different characteristic time scales on the photo-excitation
strength shows contrasting behaviours. On one hand, the rise time of the
photo-emission intensity above the Fermi energy becomes shorter with
increasing excitation strength. On the other hand, the relaxation time of the
electron population towards equilibrium becomes longer with increasing
fluence, reaching saturation above a critical value.

We corroborate these experimental findings by carrying out numerical
simulations accounting for the non-equilibrium response of a two-band tight-
binding system, with interband Coulomb repulsion and an excitonic ground
state. The calculations confirm closing of the band gap for high excitation
strengths and the speed up of the photoemission signal in the conduction band
with increasing fluence. They also allow us to analyse in detail the dynamics
of the order parameter, which can be separated into two parts: a) a fast
initial decay that speeds up with increasing fluence, b) a longer time scale
that becomes slower on approaching the critical fluence, as shown in previous
theoretical works Golež  _et al._ (2019); Tsuji  _et al._ (2013); Babadi  _et
al._ (2015). We provide a simple explanation for these different trends,
showing their compatibility. In addition, we analyse the effect of electron-
phonon coupling and show that it modifies only the quantitative features of
these dynamics, while the qualitative aspects remain unchanged. Finally, we
introduce a simple phenomenological model that captures the above dynamics and
provides an intuitive picture of the physics after photoexcitation in an
excitonic insulator. Our combined experimental/theoretical study supports a
pre-eminently excitonic character of the order parameter in Ta2NiSe5.

##  II Experimental setup

The dynamics of photo-excited electrons in Ta2NiSe5 was probed by means of
trARPES. The measurements were performed at the CITIUS high-harmonic
generation (HHG) light source Grazioli  _et al._ (2014). The system is driven
by a mode-locked Ti:Sapphire laser delivering 800-nm pulses, with a duration
of 40 fs at a repetition rate of 5 kHz. The driving laser was split in two
beams: the major part of the intensity was used to generate extreme-
ultraviolet (EUV) probe pulses through HHG, with Ar as the generating medium,
and the remaining part was used as the pump. The intensity of the pump pulses
on the sample was controlled with a variable attenuator - in all experimental
plots, the fluence refers to the peak energy density (in mJ/cm2), determined
from the expression $2E_{p}/(\pi w^{2})$, where $E_{p}$ is the energy per
pulse and $w$ is the beam waist at the sample position. Both pump and probe
pulses were p-polarized, with an in-plane component parallel to the atomic
chains of Ta2NiSe5 ($\mathbf{E}_{in}\parallel\mathbf{a}$) and an out-of-plane
component perpendicular to the sample surface
($\mathbf{E}_{out}\parallel\mathbf{b}$). The photon energy of the probe was
selected by a monochromator grating with off-plane geometry, which preserved
the pulse duration Poletto  _et al._ (2009). Because Ni 3d and Se 4p electron
orbitals have a high photo-ionisation cross-section at photon energies close
to 26 eV, a probe energy of 26.35 eV (harmonic 17 of the seed laser) was used
in the experiments. In order to preserve the ultrafast response, the energy
resolution of the source was limited to about 150 meV. This allowed us to
achieve a temporal resolution of around 50 fs. The energy and momenta of the
photo-emitted electrons were measured by a R3000 hemispherical analyser from
VG Scienta. A closed-cycle Helium cryostat was used to control the sample
temperature. The samples were cleaved in situ at a pressure lower than 6x10-9
mbar to generate a clean mirror-like surface, prior to each measurement. The
experiments were performed at a pressure lower than 1x10-10 mbar.

Figure 1: (a) Static (no pump) ARPES spectrum of Ta2NiSe5 at 100 K, measured
using a photon energy of 26.35 eV (left-most panel), and trARPES snapshots at
$\Delta t$ = 150 fs after photoexcitation for different pump fluences. The
Fermi level is marked by the horizontal dashed line; (b) momentum-integrated
EDCs at different delays ($\Delta t$), showing the temporal evolution of the
gap dynamics for different pump fluences. The vertical dashed lines (at
$\Delta t$ = -450 fs) indicate the considered region for momentum integration.

##  III Experimental results

We investigated the evolution of the topmost valence band as a function of the
incident pump fluence. This valence band is formed mainly by hybridized Ni
3$d$ and Se 4$p$ orbitals. As discussed in previous ARPES studiesWakisaka  _et
al._ (2009); Seki  _et al._ (2014), the band flatness is a possible indication
of an excitonic phase in the material, which sets in below Tc = 328 K. The
equilibrium (no pump) ARPES spectrum of Ta2NiSe5 around $\Gamma$ (Brillouin
zone center)-point at a temperature of 100 K is shown in Fig. 1(a) (left-most
panel). The Fermi level EF, which is represented by the horizontal dashed
line, is determined with respect to the valence band spectrum of a clean gold
sample and is about 0.2 eV above the topmost valence band. After
photoexcitation of the system with an infrared pump pulse ($h\nu$ = 1.55 eV),
the spectral weight of the valence band decreases and is transferred to
unoccupied energy levels above EF. In Fig. 1(a), one can observe that, at a
pump fluence of 1.2 mJ/cm2, the amount of spectral weight transfer from the
topmost valence band to the conduction band is relatively small. Increasing
the pump fluence results in an increase of the spectral weight transfer. While
this increase is modest at 2.4 mJ/cm2, the spectral weight transfer is
considerably enhanced for a pump fluence of 3.6 mJ/cm2.

These observations are made more evident in Fig. 1(b), where we show the
momentum integrated energy distribution curves (EDCs) before the arrival of
the pump pulse (at $\Delta$t = -450 fs) and at different delays after
photoexcitation, as a function of fluence. The integrated angular range for
all the EDCs is indicated by the black vertical lines in the ARPES spectrum
shown in the leftmost panel of Fig. 1(a). Comparing the EDCs before and after
photoexcitation, we observe that for the lowest fluence of 1.2 mJ/cm2, the
dynamics of the valence band shows a depopulation and broadening, accompanied
by a small energy shift towards the Fermi level at subsequent delay times. At
an intermediate pump fluence of 2.4 mJ/cm2, which is above the value reported
in previous studies, the spectral weight transfer to the conduction band is
increased, while the energy shift of the valence band towards the Fermi level
remains relatively small. Finally, at a pump fluence of 3.6 mJ/cm2, the energy
shift of the valence band towards EF is pronounced and there is a significant
increase of the photoemission intensity up to several hundreds of meV above
the Fermi level. These observations suggest that a strong pump fluence can
induce a phase transition in which the energy gap closes and the system is
driven from an insulating to a semi-metallic state Okazaki  _et al._ (2018b);
Tang  _et al._ (2020).

The key information that can be extracted from the above results is the time
scale after photoexcitation, within which the band gap closes. This provides
details on the origin of the band gap Hellmann  _et al._ (2012), i.e., whether
the dominant interaction governing the band dynamics is electronic or
phononic. Figure 1(b) shows that, at a pump fluence of 3.6 mJ/cm2, the band
gap closes within $\sim$100 fs, a time interval that is shorter compared to
the relevant vibrational time, which is about 500 fs for the 2 THz phonon
mode) Mor  _et al._ (2018); Subedi (2020); Kim  _et al._ (2020a, b); Volkov
_et al._ (2020). Even faster response is observed at the highest pump fluence
used in this study; see Fig. 3. This observation suggests that the induced
phase transition has a strong electronic character. Notice that the gap
closing is seen only for the highest fluence in Fig. 1, while for the two
lower fluences, only marginal energy shifts are observed. This suggest the
existence of a critical fluence and might also help explaining the relatively
small energy shifts observed in some of the previous trARPES experimentsMor
_et al._ (2017); Baldini  _et al._ (2020). We will discuss this issue further
in the theoretical analysis.

Figure 2: (a) trARPES spectrum around the $\Gamma$ point for the integrated
angular range indicated by the vertical dashed lines in Fig. 1(a). (b)
Temporal evolution of the integrated intensity (filled blue circles) in the
energy interval [0.0, 0.2] eV, shown by the horizontal dashed lines in (a).
The red curve represents the fit (see text for details). The measurements were
acquired at a pump fluence of 3.6 mJ/cm2. Inset: Fourier analysis of the
oscillatory part of the photoemission signal.

.

We note here that the above findings are not fully conclusive. For example,
one could argue that the relevant time scale for a pure structural phase
transition is a quarter of the phononic oscillation only Ligges  _et al._
(2018), which is comparable to the observed time scale in Fig. 1. In addition,
the structural and electronic dynamics are coupled: the structural monoclinic
displacement and the excitonic order act in a similar way on the electronic
dispersion Mazza  _et al._ (2020); Subedi (2020). Coherent oscillations in the
photoemission spectrum with the relevant phonon frequency are therefore
expected. Indeed, our data show a modulation of the integrated photoemission
signal in the energy range from 0 to 0.2 eV, see Fig. 2b). However, the fact
that the most intense phonon peak is located at 2 THz suggest that the lattice
remains in its low-temperature phase even under relatively high fluences used
in our experiment Werdehausen  _et al._ (2018); Mor  _et al._ (2018); Andrich
_et al._ (2020); Kim  _et al._ (2020a). This indicates that the dominant order
is not structural, as argued earlier in Ref. Tang  _et al._ , 2020. An
excitonic origin of the band gap is further supported below, where we discuss
the fluence dependence of the characteristic time scales observed in the
photoemission signal.

Figure 2(a) shows the angle-integrated photoemission spectral intensity as a
function of electron binding energy and time delay between the pump and probe
pulses for a pump fluence of 3.6 mJ/cm2. The temporal evolution of the
integrated intensity in the energy interval [0, 0.2] eV is plotted in Fig.
2(b). We observed two characteristic time scales: a fast rise time,
$\tau_{\mathrm{r}}$, of the angle- and energy-integrated photoemission signal
above Fermi energy and a slow relaxation time, $\tau_{\mathrm{d}}$, of the
electron population towards equilibrium. The intensity profile is fitted using
the following function:

$I(t)=H(t)[Ae^{(-t/\tau_{d})}+B],$ (1)

where H(t) describes the finite rise-time of the intensity which is modelled
by an error function, given by,

$H(t)=1+\mathrm{erf}(t/\tau_{c});$ (2)

A and $\tau_{d}$ are the amplitude and the decay time of the electronic
response, respectively, and B is a constant. Here, $\tau_{c}$ is the
convolution of the natural electron-intensity rise time $\tau_{r}$ and the
pump-probe cross-correlation time, which represents the finite temporal
resolution of our system.

The energy-integrated, time-resolved photoemission intensity at different pump
fluences is shown in Fig. 3(a). The rise time shows obvious shortening with
increased excitation strengths, as emphasized by the dashed line marking the
peak position. In Fig. 3(b), we plot the rise time $\tau_{r}$, determined from
fits shown in Fig. 3(a), as a function of pump fluence. The rise time has a
maximum value of about 120 fs for the lowest pump fluence and decreases down
to around 60 fs at the highest pump fluence used in the experiment. A fluence-
dependent rise time further supports an excitonic origin of the order
parameter in Ta2NiSe5 and will be interpreted in terms of an enhanced
scattering probability between photo-excited electrons and excitons in the
theoretical section. Although some degree of fluence dependence of the rise
time might also occur in the purely phononic scenario, the phonon modes are
expected to soften upon photoexcitation Murray  _et al._ (2007), which would
result in a trend opposite to the one observed in the experiment.

We also examined the fluence dependence of the electron relaxation time
$\tau_{d}$ towards equilibrium, see Fig. 3(c). Contrary to the behaviour of
$\tau_{r}$, the relaxation time $\tau_{d}$ shows a monotonic increase of up to
the critical fluence, after which saturation is observed. This trend is
consistent with previous experiments Okazaki  _et al._ (2018a) and we
attribute it to the scattering of phonons with excitons: at high excitation
fluences a substantial population of excited vibrational modes can transfer
energy back to the perturbed excitonic condensate and suppress its
reestablishment Hedayat  _et al._ (2019).

Figure 3: (a) Energy-integrated photoemission intensity for several pump
fluences, where the smooth curves are fits to the experimental data. The
shortening of the rise time for increasing pump fluence is indicated by a
black dashed line. (b) Rise time $\tau_{\mathrm{r}}$ as a function of pump
fluence; (c) Dependence of the hot electron relaxation time
$\tau_{\mathrm{d}}$ on pump fluence.

##  IV Theoretical analysis

We now turn to the theoretical discussion, which will help us understand the
early stages of the photoinduced dynamics in Ta2NiSe5. In the first
subsection, we will introduce a model of an excitonic insulator. We will then
discuss the theoretical results and compare them with the experimental
findings shown in the previous sections. In the last subsection we will
develop a simple phenomenological framework that captures the main dynamics
after photoionization in an excitonic insulator.

###  IV.1 Model and method

We consider a two-band tight-binding model in one dimension, which is
appropriate for Ta2NiSe5, due to the strong anisotropy along the atomic chains
Seki  _et al._ (2011, 2014). For simplicity, we consider spinless fermions and
couple the electrons in the two bands by a local interband repulsion. The
total Hamiltonian is given by

$\begin{split}H=H_{\text{kin}}+H_{\text{int}}+H_{\text{dip}}+H_{\text{ph}}.\end{split}$
(3)

The kinetic energy reads

$\displaystyle H_{\text{kin}}=\sum_{\begin{subarray}{c}k\\\
\alpha\in\\{0,1\\}\end{subarray}}h_{k,\alpha}c_{k,\alpha}^{\dagger}c_{k,\alpha}$
(4)

where $h_{k,\alpha}=-2J_{\alpha}\cos(k)-\mu+(-1)^{\alpha}\Delta\epsilon$ is
the single-particle dispersion relation, $\Delta\epsilon$ the band-level
splitting, and $J_{0}=-J_{1}$ the hopping integral; $c_{k,\alpha}^{\dagger}$
and $c_{k,\alpha}$ are, respectively, the creation and annihilation operators
of an electron with momentum $k$ in the orbital $\alpha$. The chemical
potential $\mu$ is chosen such that the system is at half filling. The
interaction part includes the local interband density-density interaction with
a strength $V$:

$\begin{split}H_{\text{int}}=V\sum_{i}n_{i,0}n_{i,1},\end{split}$ (5)

where $n_{i,\alpha}=c_{i,\alpha}^{\dagger}c_{i,\alpha}$ is the electronic
density at site $i$ and band $\alpha$. The Hamiltonian conserves the charge in
individual bands, but at sufficiently low temperatures the associated symmetry
is spontaneously broken and the corresponding order parameter $\phi=\langle
c_{0}^{\dagger}c_{1}\rangle$ acquires a finite value. In the previous
expression, the operator index refers to the considered orbital, while the
particle momentum is averaged. The order parameter $\phi$ is a complex
quantity and represents the degree of spontaneous hybridization between the
conducting and the valence bands (i.e., the number of condensed excitons).

We simulate the coupling of the electromagnetic field to the electrons by the
term

$\begin{split}H_{\text{dip}}=E(t)\sum_{i}c_{i,1}^{\dagger}c_{i,0}+\text{h.c.},\end{split}$
(6)

where the index $i$ runs over all sites and $E(t)$ is the time-dependent
amplitude of the electric field, which is responsible for the dipolar
excitation between the conduction and valence bands.

In previous studies, the time evolution of the excitonic order was described
within the time-dependent Hartree-Fock theory Murakami  _et al._ (2017);
Fujiuchi  _et al._ (2019), which is however only applicable for low photo-
doping and short times, due to the lack of scattering and heating effects. In
this work, we consider a strong photo-doping regime, where the order is
completely destroyedGolež  _et al._ (2016). For this purpose, we employ the
second Born approximation, describing the non-equilibrium dynamics within the
Keldysh formalism. This approach is non-perturbative in the strength of the
electric field and properly describes electron-electron scattering, as well as
the associated electron thermalization Balzer  _et al._ (2010); Schlünzen  _et
al._ (2019); Tuovinen  _et al._ (2019, 2020).

The phonon term in the Hamiltonian is given by

$\begin{split}H_{\text{ph}}=g\sum_{i}(b_{i}+b_{i}^{\dagger})(c_{i,1}^{\dagger}c_{i,0}+c_{i,0}^{\dagger}c_{i,1})+\omega_{0}\sum_{i}b_{i}^{\dagger}b_{i}.\end{split}$
(7)

Here $b_{i}^{\dagger}$ and $b_{i}$ are, respectively, the phonon creation and
annihilation operators and $g$ is the electron-phonon interaction strength. We
consider the case of a dispersionless phonon with frequency $\omega_{0}$. We
choose the form of the electron-phonon coupling in a way that is relevant for
the mode involved in the orthorhombic to monoclinic transition in
Ta2NiSe5Kaneko  _et al._ (2013), namely, a nonvanishing phonon distortion
leads to a hybridization of the two bands, and thus to a finite value of
$\phi$. We treat the electron-phonon interaction on the mean-field level
leading to a self-consistent classical equation of motion, see Appendix A for
details.

Figure 4: Snapshots of the simulated photoemission spectrum Ik($\omega$,tp):
a) at equilibrium; b) at a fixed pump-probe delay tp=70 fs and $E_{0}=0.4$,
i.e., below the value corresponding to the dynamical phase transition; c) at
tp=70 fs and $E_{0}=0.6$, above the critical value corresponding to gap
closure. The Fermi level corresponds to $\omega=0$. The model parameters are
$V/J_{0}=3.0$, $\lambda/J_{0}=0.0$ and the equilibrium state is prepared at
$T=100$ K. The white dashed lines represent the valence band dispersion in
equilibrium with an artificially suppressed order parameter. d) Dependence of
the order parameter $|\phi^{F}|=|\phi(t\rightarrow\infty)|$ (red curve) and
the polarization $\Delta n$ (black curve) on the excitation strength $E_{0}$
in the extrapolated long-time limit .

The electron-phonon interaction includes two free parameters, namely the
phonon frequency $\omega_{0}$ and the interaction strength $g$. We fix the
phonon frequency to $\omega_{0}$=0.013 eV Larkin  _et al._ (2017); Mor  _et
al._ (2018); Werdehausen  _et al._ (2018). As both the electron-electron $V$
and the electron-phonon $g$ interaction can contribute to the gap opening, we
have adjusted these two parameters so that they produce the same single-
particle gap. At a mean-field level, the gap size due to the excitonic and
lattice effects is governed by the effective interaction $V^{*}=V+2\lambda,$
where $\lambda=2g^{2}/(\omega_{0})$. We have checked that the effective
interaction $V^{*}$ is a good measure of the gap size in the framework of the
second Born approximation and all subsequent calculations are done with fixed
$V^{*}/J_{0}=3.0$. We also fixed the effective interaction $V^{*}$ and the
crystal field $\Delta\epsilon$, such that the system is in the BCS parameter
regime, where the normal (disordered) phase is semimetallic.

###  IV.2 Theoretical results and comparison with experiment

We have adjusted the theoretical calculations to the equilibrium ARPES
spectrum and set the hopping parameter to $J_{0}=0.26eV$ Mor  _et al._ (2017),
which means that the inverse hopping corresponds to $\approx 2.5$ fs. For
easier comparison with the experimental part, we report energies in units of
eV and times in fs. The system is excited with a short laser pulse that is
parameterized as

$E(t)=E_{0}\sin[\omega(t-t_{0})]e^{-4.6(t-t_{0})^{2}/t_{0}^{2}},$ (8)

where $t_{0}=2\pi/\omega$ and $\hbar\omega$=1.5 eV, consistent with the
experimental value. The pulse strength E0 fixes the amount of photo-doping
(polarization). With these choices, the pulse duration is about $t\approx 5$
fs.

We start our analysis by simulating the time-dependent photoemission spectrum,
see Fig. 4, which can be directly compared to the experimental data in Fig. 1.
Its expression is $I_{k}(\omega,t_{p})=-\text{Im}[\sum_{\alpha=0}^{1}\int
dt^{\prime}S^{2}(t^{\prime})\exp^{\mathrm{i}\omega(t^{\prime}-t)}G^{<}_{k,\alpha\alpha}(t^{\prime},t),$
where $G^{<}_{k,\alpha\alpha}$ is the lesser component of the Green’s function
(see Appendix A) and $S(t)=\exp(-t^{2}/2\delta^{2})$ is the envelope of the
probe pulse Freericks  _et al._ (2009) with a duration $\delta=40$ fs, in
agreement with the experimental setting. At equilibrium, see Fig. 4(a), only
the states below the Fermi level are occupied and their dispersion shows a
two-peak structure. This is characteristic of the so-called BCS regime, in
which the excitonic correlation length is larger than the lattice spacing, and
the corresponding disordered phase is semi-metallic. The excitonic
hybridization between the conduction and valence bands opens up a gap at the
degeneracy points Zenker  _et al._ (2012); Kaneko  _et al._ . After
photoexcitation, the electrons are promoted from the valence to the conduction
band. The valence band broadens and moves towards the Fermi level, see Fig.
4(b). Above the critical value $E_{0}=0.51$, the gap is suppressed, see Fig.
4(c), and one observes a semi-metallic dispersion, which looks consistent with
the experimental results corresponding to the largest excitation strength in
Fig. 1.

The reduction of the gap is driven by the suppression of excitonic order. The
value of the order parameter after a time much longer than the duration of the
pump envelope is displayed in Fig. 4(d). One observes a progressive
suppression of the order parameter with fluence. Notice that only on
approaching the critical fluence, the associated suppression of the gap
becomes sizeable, which helps interpreting recent experiments Mor  _et al._
(2017); Baldini  _et al._ (2020) that observed marginal effects on the
spectra. Figure 4(d) also shows the photo-induced polarization

$\Delta
n(t)=\sum_{\alpha=0}^{1}\int_{0}^{\infty}d\omega[G^{<}_{\text{loc},\alpha\alpha}(\omega,t)-G^{<}_{\text{loc},\alpha\alpha}(\omega,0)]$
(9)

in the long-time limit. Increasing fluences leads to higher occupancy of the
conduction band, which corresponds to higher $\Delta n$ values.

We now turn to time dependence and start with the photoemission spectra.
Figure 5(a) shows snapshots of the momentum integrated photoemission spectrum
$I_{\text{loc}}(\omega,t_{p})=\frac{1}{N_{k}}\sum_{k}I_{k}(\omega,t_{p})$,
where $N_{k}$ (=128) is the number of momentum points we used in the
simulation, at different delay times tp. The data indicate a reduction of the
gap size, and a substantial population of the conduction band. The evolution
of the electronic population above the Fermi level shows a fast (order of 100
fs) relaxation. The number of photo-excited electrons, $\Delta n(t)$,
increases even after the photoexcitation, Fig. 5(b). We interpret this result
by considering the scattering of electrons having excess kinetic energy with
excitons. The resulting exciton dissociation generates two electrons and one
hole. This leads to an increase in the number of conducting electrons, as it
was previously observed in Ref. Golež  _et al._ , 2016. This anomalous
increase in the number of conducting electrons was also observed in a
different excitonic insulator candidate, namely 1T-TiSe2, by analysing the
growth of the photo-induced Drude weight in the time-resolved optical
conductivity Mathias  _et al._ (2016).

For a direct comparison with the experiments reported in the previous section,
we have integrated the momentum-integrated photoemission spectrum over the
energy range [0-0.2] eV, as $\langle
I_{\text{loc}}\rangle(t_{p})=\int_{0}^{0.2}d\omega
I_{\text{loc}}(\omega,t_{p})$. In Fig. 5(c), we plot $\langle
I_{\text{loc}}\rangle(t_{p})$, normalized to its maximal value, for several
excitation strengths. For weak excitation (i.e., $E_{0}\leq 0.2$), $\langle
I_{\text{loc}}\rangle$ shows a small overall change and is rather oscillating
in time (data not shown), with a frequency corresponding to the amplitude mode
excitation, i.e., $\approx$ 8 fs, well below our experimental resolution. For
excitation strengths closer to the critical value, the redistribution of the
spectral weight is larger and the oscillations are suppressed. Interestingly,
at large fluences the dynamics of the intensity build-up shows a clear speed-
up, consistent with experimental findings.

In order to address this aspect more quantitatively, we fitted the curves
shown in Fig. 5(c) with a single exponential
$I(t)=I_{0}(1.0-\exp(-t/\tilde{\tau}_{r}))$, and extracted the rise time
$\tilde{\tau}_{r}$. The latter is shown in Fig. 5(d), as a function of
polarization $\Delta n$, for different strengths of the electron-phonon
coupling $\lambda$; the vertical dashed lines indicate critical fluences. Let
us first focus on the curve corresponding to $\lambda=0$, i.e., no electron-
phonon coupling; we will turn to the other two cases later in the text.

We note that both the magnitude of $\tilde{\tau}_{r}$ and the functional
dependence with the excitation strength (polarization) are in good agreement
with measurements, see Fig. 3(b). The fluence dependence of the rise time is
relatively weak for low fluences and becomes stronger at larger ones, showing
a very weak trend in the region around the critical fluence. The observed
discrepancy in the magnitude of experimental and theoretical $\tau_{r}$ values
can be explained by the strong dependence of the latter on the microscopic
parameters, such as the Coulomb interaction $V$, the electron-phonon
interaction $g$, and the electronic hopping $J_{0}$. The general trend of a
faster rise time with increasing fluence is surprising due to the expected
slowdown of the dynamics close to the critical point Tsuji  _et al._ (2013);
Golež  _et al._ (2016); Dolgirev  _et al._ (2020); Zong  _et al._ (2019b).
Before explaining this result, we contrast it to the dynamics of the order
parameter.

Figure 5:  (a) Time evolution of the momentum averaged component of the
photoemission spectra, after an excitation with a strength of $E_{0}=0.4$. The
dashed line corresponds to the equilibrium spectrum (i.e., no pump). (b) Time
evolution of the polarization $\Delta$n after excitation with $E_{0}=0.4$ and
$E_{0}=0.6$. (c) Time evolution of the photoemission integral over the energy
range [0, 0.2] eV, for different excitation strengths. (d) Rise time
$\tau_{I}$ of the energy integrated photoemission signal versus the excitation
strength, for the purely electronic theory ($\lambda/J_{0}=0.0$) and in the
presence of a finite electron-phonon coupling ($\lambda/J_{0}=0.03$ and
$\lambda/J_{0}=0.1$). The vertical lines mark the critical excitation
strength, above which the order parameter, in the long-time limit, goes to
zero. (e) Time evolution of the full-width at half maximum $\sigma$ for the
dominant peak in the photoemission spectra (below the Fermi energy), at
different excitation strengths. (f) Rise time $\tau_{\sigma}$ of the full-
width at half maximum $\sigma$ for the purely electronic theory
($\lambda/J_{0}=0.0$) and in the presence of a finite electron-phonon coupling
($\lambda/J_{0}=0.03$ and $\lambda/J_{0}=0.1$).

As the order parameter is a complex number
$\phi=|\phi|\exp[\mathrm{i}\theta]$, we will separately follow its absolute
value $|\phi|$ and phase $\theta$. Fig. 6 shows the time dependence of
$|\phi|$ and $\theta$, as predicted by the model. After photo-excitation,
$|\phi|$ diminishes with time until it levels off to a final, fluence-
dependent value. For a large enough fluence, after a certain time, $|\phi|$
vanishes, which indicates the closing of the gap (dynamical critical point).
With increasing fluence, the overall time scale becomes longer on approaching
the critical fluence, see Fig. 6(a). Intriguingly, the behaviour of the order
parameter around the critical fluence, which was remarked also in different
contexts Tsuji  _et al._ (2013); Golež  _et al._ (2016); Dolgirev  _et al._
(2020); Zong  _et al._ (2019b), is in contrast with that manifested by the
rise time of the electron population in the conduction band, see Fig. 5(d).
This raises the following question: Can this slowing down also be observed in
the photoemission spectra? It turns out that the momentum-integrated spectra
do show such a behaviour. We extract the full width at half maximum, $\sigma$,
of the photoemission peak below the Fermi level, see Fig. 5(e). The time
dependence of this quantity does become slower with fluence, see fitted rise
times in Fig. 5(f). The fact that $\sigma$ follows the dynamics of the order
parameter can be understood by recalling that we are considering momentum
integrated spectra: the suppression of the order parameter shows up in the
modification of the dispersion that brings the upper edge of the valence band
and lower edge of the conduction band closer; this leads to an increased width
of the two (momentum-integrated) bands. Due to the limited energy resolution
of our experiment, this theoretical prediction could not be directly verified
with our present setup. We propose it as the subject of a future experimental
study. One could interpret this behaviour as a manifestation of the non-
thermal criticality, as discussed in previous theoretical studies Tsuji  _et
al._ (2013); Golež  _et al._ (2016); Dolgirev  _et al._ (2020).

While we have exposed the discrepancy between the behaviour of the signal’s
rise time and the order parameter vs. excitation strength, we would like to
understand the origin of this dichotomy. If we look at early stages after the
excitation ( $t<$ 50 fs) in Fig. 6(a), the absolute value of $d|\phi(t)|/dt$
increases with excitation strength, showing a speedup similar to the one
observed for the rise time $\tilde{\tau}_{r}$ (there is no critical slowing
down). In the later stages, this dynamics is altered and we observe critical
slowing down of the order parameter with increasing fluence. However, we note
that the absolute change of the order parameter in the early stages dominates
over that in the later stage when the critical slowing down sets in. The
short-time behaviour of the order parameter is consistent with a simple
picture, where the condensed excitons act as a heat bath for relaxation of the
photoexcited electrons. The energy transfer from electrons to the excitons
leads to the decay of excitons and the process is faster if more electrons are
excited; in this case, the order parameter also decays faster.

Based on the above considerations, we believe that the rise-time of the
photoemission signal is related to the early-stages of the dynamics of the
order parameter, where no critical slowing down is observed. The increase of
the intensity in the momentum and energy-integrated spectrum can originate
from two processes; i) relaxation of the excited electrons to lower energies
due to scattering and ii) shift of the conduction band bottom into the
integrated energy range. In i), if scattering of the high-energy electrons
with condensed excitons is the dominant relaxation process, this relaxation
stops when the excitons are completely consumed. Therefore, under strong
photoexcitation, $|\phi|$ initially decays faster, and as the exciton bath is
depleted, the relaxation also stops rapidly. In ii), because the order
parameter is proportional to the band gap, the increase of the signal in the
integrated region is directly proportional to the early-stages of the
evolution of $|\phi|$. These two processes are therefore cooperating to
shorten $\tilde{\tau}_{r}$ at high fluences.

We now discuss how the electron-phonon interaction modifies the above picture.
In Fig. 6, we compare the behaviour for the purely electronic case
($U/J_{0}$=3, $\lambda/J_{0}$=0) with the case including a weak electron-
phonon interaction ($U/J_{0}$=2.94, $\lambda/J_{0}=0.03$). As a first remark,
we note that a nonzero $\lambda$ value makes the order parameter $|\phi|$ more
robust against photo-excitation, with the difference compared to the
$\lambda=0$ case being more visible for larger excitation strengths, see Fig.
6a. The importance of the electron-lattice interaction is even more dramatic
for the phase of the order parameter $\theta$, see Fig. 6b. The latter is
either freely oscillating in the case $\lambda=0$, or confined in the case of
electron-phonon coupling with $\lambda/J_{0}=0.03.$ The deconfinement of the
phase in the electron-lattice coupled situation for the strongest excitation,
namely $E_{0}=0.53,$ coincides with an almost complete destruction of the
order. The robustness of the order parameter in the electron-phonon coupled
systems can be understood as a combination of two effects: i) slower response
of the lattice contribution to the gap closing; ii) opening of the finite mass
of the phase mode of the excitonic order parameter, which prevents dephasing
Murakami  _et al._ (2017). The latter effect only takes place for couplings
that break the $U(1)$ gauge symmetry and make the phase mode massive Zenker
_et al._ (2014); Murakami  _et al._ (2020); Guseǐnov and Keldysh (1973);
Littlewood and Zhu (1996). As the electron-phonon interaction is increased,
$|\phi|$ becomes more robust after photo-excitation, and the phase more
confined.

Figure 6: Time evolution of the absolute value of the order parameter $|\phi|$
(a) and its phase $\theta$ (b), for different excitation strengths $E_{0}$.
The dashed lines represent the evolution predicted by a purely electronic
theory ($U/J_{0}$=3, $\lambda=0$), while the solid lines are the predictions
including the electron-lattice coupling ($U/J_{0}$=2.94,
$\lambda/J_{0}=0.03$). The vertical dashed-dotted line in a) separates the
short and long time regions. For clarity, we have omitted the plot for
$E_{0}=0.6$ in the purely electronic case in a).

Coming now to the effect of the electron-phonon interaction on the dynamics of
the rise time $\tilde{\tau}_{r}$, shown in Fig. 5(d), we considered the cases
of different electron-phonon coupling strengths, namely $U/J_{0}$=2.94,
$\lambda/J_{0}=0.03$ and $U/J_{0}$=2.8, $\lambda/J_{0}=0.1$. The first
observation is that the critical excitation strength increases with $\lambda$:
$\Delta n=0.057$ for $\lambda/J_{0}=0.03$ and $\Delta n=0.07$ for
$\lambda/J_{0}=0.1$, to be compared with $\Delta n=0.053$ for the case
$\lambda=0$. This theoretical finding, together with the experimental
observation that the gap closes only after a strong photo-excitation,
indicates that electron-phonon coupling might play an important role in
Ta2NiSe5, and provides a possible explanation for the observed small energy
shifts in the trARPES spectra in Refs. Mor  _et al._ , 2017; Baldini  _et al._
, 2020 . Also in the cases with $\lambda\neq 0$, the rise time shows a speed-
up with increasing photo-excitation and the shape of the curve is similar to
that found experimentally. Note however that in the curve corresponding to the
largest $\lambda$ value, the relaxation time does not show any significant
speedup with increased excitation strength below the critical one, in contrast
to the experimental data. This gives an upper limit of the theoretical
estimate of the electron-phonon coupling in the material. Similar conclusions
are obtained for the time scale of the broadening $\tau_{\sigma},$ see Fig.
5(f). The main signature of the increased electron-phonon coupling is an
increase of the critical excitation strength.

###  IV.3 Phenomenological description

As the microscopic calculations presented in previous sections are rather
involved, we will introduce a phenomenological model, which provides a simple
explanation of the initial speedup of order parameter dynamics with increasing
fluence, followed by a slowdown when approaching the critical excitation
strength. The basic assumption of the model is that the condensed excitons act
as a heat bath for photoexcited electrons. We have confirmed this assumption
by explicit calculations of photo-excitation in the semi-metallic and
semiconducting states without excitonic order and found that the relaxation is
much weaker and slower as compared to the relaxation in the case of an
excitonic insulator (data not shown).

The model considers, besides the order parameter $\phi$ (for brevity, we will
use the notation $\phi$ instead of $|\phi|$ in the following), also excited
electrons that are characterized by a density $\mathcal{N}$. The relaxation of
electrons is assumed to dominantly take place through coupling with excitons
and vice-versa - the relaxation of the excitonic order parameter occurs
through coupling with $\mathcal{N}$. The probability of the relaxation process
is proportional to the number of excitations and to the number of excitons,
i.e. to $\mathcal{N}$ and $\phi$. Hence, the two equations that govern the
behaviour of the two coupled variables can be written as

$d\phi/dt=-c_{\phi}\phi\mathcal{N},$ (10)
$d\mathcal{N}/dt=-c_{\mathcal{N}}\phi\mathcal{N},$ (11)

where the $c_{\phi}$ is the exciton dissociation rate due to scattering with
high-energy electrons and $c_{\mathcal{N}}$ is the decay rate of excited
electrons. Let us first discuss the extreme cases of very weak and very high
fluences. In the former case, the value of the order parameter remains close
to the initial one, $\phi_{I}$. Conversely, $\mathcal{N}$ is exponentially
damped and its relaxation time is given by $1/(c_{\phi}\phi_{I})$, i.e., it
becomes short for large values of the order parameter $\phi$. In the latter
case, $\mathcal{N}$ and $\phi$ exchange their roles and now $\phi$ is
exponentially damped and its relaxation time is given by
$1/(c_{\mathcal{N}}\mathcal{N}_{I})$, i.e., it becomes short for large values
of excitation $\mathcal{N}$.

Away from those two extremal regimes, one can integrate the two equations.
Because the two time derivatives are proportional to each other,
$d\mathcal{N}/dt=(c_{\mathcal{N}}/c_{\mathcal{\phi}})d\phi/dt$ (in this
approximation, the relaxation of an electron implies dissociation of an
exciton), $\phi$ and $\mathcal{N}$ are linearly related

$\mathcal{\phi}=c_{\phi}/c_{\mathcal{N}}\mathcal{N}+\mathrm{const.}$ (12)

This enables rewriting of the problem into a simple differential equation that
can be integrated. For weak fluences (below the critical one), the long-time
value of $\mathcal{N}$ vanishes, but the order parameter reaches a
nonvanishing value $\phi_{F}$, hence one takes
$\mathcal{N}=(\phi-\phi_{F})c_{\mathcal{N}}/c_{\phi}$. This leads to a
solution (denoting the initial condition as $\phi(t=0)=\phi_{I}$)

$\phi(t)=\frac{\phi_{I}\phi_{F}e^{t/\tau_{\phi}}}{\phi_{F}-\phi_{I}+\phi_{I}e^{t/\tau_{\phi}}},$
(13)

where $1/\tau_{\phi}=c_{\mathcal{N}}\phi_{F}$. At long times, $\phi$ thus
decays towards $\phi_{F}$ with a decay time proportional to $1/\phi_{F}$,
which directly implies a critical slowing down on approaching the critical
fluence. On the other hand, at short times, the decay is faster than that
inferred from continuation of exponential decay to short times. More
precisely, the initial slope $d\phi/dt=-c_{\mathcal{N}}\phi_{I}\Delta\phi$ is
proportional to the change of the order parameter
$\Delta\phi=(\phi_{I}-\phi_{F})$; i.e., the speed of the initial decay is
higher for increasing fluence, which explains the short-time behaviour of the
order parameter and the rise time in numerical simulations.

Figure 7: (a) Absolute value of the order parameter $\phi$ (full) fitted with
a heuristic function (dashed) and compared also to an exponential decay
towards the long-time value $\phi_{F}$ (dashed-dotted). (b) Logarithmic plot
of $\phi-\phi_{F}$. (c) Parameterless fits to data.

The fit of Eq. 13 to the numerical data is shown in Fig. 7(a). One sees that
at short times $\phi$ indeed decays faster, consistent with the heuristic
function, and yet the slope of the approach towards $\phi_{F}$ at long-times,
shown in Fig. 7(b), becomes smaller; the order parameter dynamics slows down
as the long-time value $\phi_{F}$ becomes small (note that the the slope is
reproduced only qualitatively). Fitted $c_{\mathcal{N}}$ depend on each case
very weakly and a single value can be used for all fluences. This is
demonstrated in Fig. 7(c) where we took $c_{\mathcal{N}}$=3.7 for all the
fluences. Note that $\phi_{F}$ is independently determined from the long-time
value, and there are thus no free parameters.

Finally, it is interesting to consider also fluences above the critical one,
where $\phi_{F}=0$ and there is a remaining density of excitations
$\mathcal{N}_{F}$ in the limit of long times. In this case $\phi$ decays with
a rate $1/\tau_{N}\propto\mathcal{N}_{F}$, that is, the decay becomes faster
at higher fluences, which is indeed what one sees numerically, cf. Fig. 6.
Moreover, one can estimate $\mathcal{N}_{F}$ for fluences above the critical
one from the behaviour for fluences below the critical one and find that the
long-time tails indeed decay proportionally to $\mathcal{N}_{F}$.

The key observations of the phenomenological analysis are thus: (i) the
critical slowing down of the dynamics of $\phi$ occurs on approaching the
critical fluence and is described naturally by the lack of excitons acting as
an efficient bath, (ii) this slowing down of the dynamics occurs in the limit
of long times only and is related to the fact that the time scale is given by
the inverse value of the order parameter (which itself becomes small), (iii)
the dominant change of the order parameters happens at short times, where the
derivative is independent of time for fluences below the critical one, (iv)
for fluences above the critical one, the evolution becomes faster also
initially.

##  V Conclusion

Summing up, our experimental results demonstrate that a phase transition in
Ta2NiSe5 from an insulating to a semi-metallic state can be induced under
strong photoexcitation, which belongs to a fluence regime that has not been
explored before(Mor  _et al._ , 2018; Okazaki  _et al._ , 2018a). We observe
that the closing of the band gap occurs on a time scale faster than the
typical phonon oscillation in the material. Furthermore, we find that the
(fast) rise time of photo-emission intensity above the Fermi energy and the
(slow) relaxation time of the electron population towards equilibrium show
contrasting behaviours as a function of the excitation strength: while the
rise time becomes shorter with increasing fluence, the relaxation time towards
the equilibrium becomes longer, reaching saturation above the critical
fluence.

The observed speed up of the photoemission signal is confirmed by the results
provided by a theoretical two-band tight-binding model, which allowed us to
interpret the experimental trends in terms of scattering of the photoexcited
electrons with the excitonic condensate. We found that the inclusion of a
phonon mode does not modify the physics qualitatively, however, it makes the
model system more robust towards photoexcitation. This might explain the
apparent discrepancies between different trARPES studies reported in the
literature Mor  _et al._ (2017); Okazaki  _et al._ (2018a); Baldini  _et al._
(2020); Tang  _et al._ (2020); Suzuki  _et al._ (2020). We note, however, that
this aspect needs to be considered in the future in a more realistic setting
with genuine symmetry of the involved phonon mode and the excitonic order
parameter Subedi (2020).

Finally, the above results point to a predominately electronic origin of the
order parameter in Ta2NiSe5 and allow us to develop simple phenomenological
model, which describes the dynamics after photoexcitation in an excitonic
insulator.

###### Acknowledgements.

D. G. and J. M. acknowledge the support of the Slovenian Research Agency
(ARRS) under Programs No. J1-2455, J1-1696, J1-2458, and P1-0044. J. M.
acknowledges insightful discussions with Alaska Subedi. The Flatiron Institute
is a division of the Simons Foundation. Y.M. acknowledges the support by a
Grant-in-Aid for Scientific Research from JSPS, KAKENHI Grant Nos. JP19K23425,
JP20K14412, JP20H05265, and JST CREST Grant No. JPMJCR1901. This work has
received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 654360 NFFA-Europe. We thank A.
Ciavardini for valuable help during experiments.

##  Appendix A Application of the time-dependent second Born approximation to
the excitonic insulator

To allow for the excitonic symmetry breaking, it is useful to introduce the
spinor

$\Psi_{k,\alpha}\equiv\begin{pmatrix}c_{k,0}\\\ c_{k,1}\\\ \end{pmatrix}$ (14)

and the corresponding 2$\times$2 Green’s function is given by

$G_{k}(t,t^{\prime})=-\mathrm{i}\langle\mathcal{T_{C}}\Psi_{k}(t)\Psi^{\dagger}_{k}(t^{\prime})\rangle$
(15)

and it is determined from a solution of the Dyson’s equation. For latter usage
we will introduce the single-particle density matrix
$\rho_{k,\alpha\beta}=\langle\Psi^{\dagger}_{k,\beta}\Psi_{k,\alpha}\rangle,$
where the Greek indices indicate the orbital space. The local component of the
density matrix is given by
$\rho_{\text{loc},\alpha\beta}=\sum_{k}\rho_{k,\alpha\beta}.$ The occupation
of the valence (conducting) band is therefore given by
$n_{0}=\rho_{\text{loc},00}~{}(n_{1}=\rho_{\text{loc},11}).$ Similarly, the
excitonic order parameter is given by the off-diagonal component of the
density matrix $\phi=\rho_{\text{loc},01}.$ The details of the theoretical
solution are presented in the following.

The employed approximation is determined by the self-energy and we will treat
the Coulomb interaction within the 2nd Born approximation, which includes the
inelastic scattering and captures the heating effects. For the Coulomb
interaction, the Hartree term is given by

$\Sigma^{H}_{k,ij}(t)=\delta_{ij}V\rho_{\text{loc},\bar{i}\bar{i}}(t),$ (16)

where the overline marks the oposite band and the Fock term is given by

$\Sigma^{F}_{k,ij}(t)=V\rho_{\text{loc},ij}(t).$ (17)

The coupling with phonons will be described on the Hartree-Fock level with the
instantaneous term

$\Sigma^{ph}_{k,ij}(t)=gX(t)\delta_{|i-j|=1},$ (18)

and the phonon motion is described by a system of differential equations for
the phononic distortion $X$ and the corresponding momentum $P$

$\begin{split}&\partial_{t}P(t)=\omega_{0}P(t)\\\
&\partial_{t}X(t)=-\omega_{0}X(t)-2g(\rho_{\text{loc},01}+\rho_{\text{loc},10}).\end{split}$
(19)

Figure 8: The interacting self-energy diagrams corresponding to the 2nd Born
approximation.

As the Hartree and Fock self-energies are instantaneous, it is useful to
redefine the single particle Hamiltonian as
$h_{k}^{HF}=h_{k}+\Sigma^{H}_{k}+\Sigma_{k}^{F}$ and introduce the
renormalized free propagator as

$G_{k}^{HF}=(\mathrm{i}\partial_{t}-h_{k}-\Sigma^{H}-\Sigma^{F}_{k}-\Sigma^{ph}_{k})^{-1}$
(20)

In the 2nd Born approximation, the correlation self-energy is given by the
expansion in the interaction strength which is then resumed using the Dyson
equation, see Fig. 8 for the diagrammatic representation of the self-energy.
The analytical expression for the self-energy is given by

$\begin{split}\Sigma^{2B}_{q,ij}(t,t^{\prime})=&-\mathrm{i}V^{2}\sum_{k}G_{k,\bar{j}\bar{i}}(t^{\prime},t)\\\
&\left[\chi_{k+q,i\bar{j}}(t,t^{\prime})-\chi_{k+q,\bar{i}\bar{j}}(t,t^{\prime})\right],\end{split}$
(21)

where we have introduced the charge susceptibility bubble

$\chi_{q,ij}(t,t^{\prime})=\mathrm{i}\sum_{k}G_{q-k,i\bar{j}}(t,t^{\prime})G_{k,\bar{i}j}(t,t’).$
(22)

The Greens function is then obtained from the solution of the Dyson equation

$G_{k}=G_{k}^{HF}+G_{k}^{HF}\ast\Sigma_{k}^{2B}\ast G_{k},$ (23)

where $\ast$ denotes the convolution in time and the matrix multiplication in
the orbital space. We solve this problem numerically using the library NESSi
employing an efficient MPI parallelization over the momentum points Schüler
_et al._ (2020).

## References

  * Mott (1961) N. F. Mott, Philosophical Magazine 6, 287 (1961), http://dx.doi.org/10.1080/14786436108243318 .
  * Knox (1963) R. Knox, Academic, New York (1963).
  * Butov  _et al._ (2002) L. Butov, A. Gossard, and D. Chemla, Nature 418, 751 (2002).
  * Eisenstein and MacDonald (2004) J. Eisenstein and A. H. MacDonald, Nature 432, 691 (2004).
  * Wang  _et al._ (2019) Z. Wang, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, Nature 574, 76 (2019).
  * Wakisaka  _et al._ (2009) Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, Physical review letters 103, 026402 (2009).
  * Wakisaka  _et al._ (2012) Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, N. L. Saini, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, J Supercond Nov Magn 25, 1231 (2012).
  * Cercellier  _et al._ (2007) H. Cercellier, C. Monney, F. Clerc, C. Battaglia, L. Despont, M. G. Garnier, H. Beck, P. Aebi, L. Patthey, H. Berger, and L. Forró, Phys. Rev. Lett. 99, 146403 (2007).
  * Mazza  _et al._ (2020) G. Mazza, M. Rösner, L. Windgätter, S. Latini, H. Hübener, A. J. Millis, A. Rubio, and A. Georges, Phys. Rev. Lett. 124, 197601 (2020).
  * Kaneko  _et al._ (2013) T. Kaneko, T. Toriyama, T. Konishi, and Y. Ohta, Phys. Rev. B 87, 035121 (2013).
  * Murakami  _et al._ (2017) Y. Murakami, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. Lett. 119, 247601 (2017).
  * Subedi (2020) A. Subedi, Phys. Rev. Materials 4, 083601 (2020).
  * Di Salvo  _et al._ (1986) F. J. Di Salvo, C. H. Chen, R. M. Fleming, J. V. Waszczak, R. G. Dunn, S. A. Sunshine, and J. A. Ibers, Journal of The Less-Common Metals 116, 51 (1986).
  * Kim  _et al._ (2020a) K. Kim, H. Kim, J. Kim, C. Kwon, J. S. Kim, and B. Kim, arXiv preprint arXiv:2007.08212 (2020a).
  * Kim  _et al._ (2020b) M.-J. Kim, A. Schulz, T. Takayama, M. Isobe, H. Takagi, and S. Kaiser, Phys. Rev. Research 2, 042039 (2020b).
  * Volkov  _et al._ (2020) P. A. Volkov, M. Ye, H. Lohani, I. Feldman, A. Kanigel, K. Haule, and G. Blumberg, arXiv preprint arXiv:2007.07344 (2020).
  * Larkin  _et al._ (2017) T. I. Larkin, A. N. Yaresko, D. Pröpper, K. A. Kikoin, Y. F. Lu, T. Takayama, Y.-L. Mathis, A. W. Rost, H. Takagi, B. Keimer, and A. V. Boris, Phys. Rev. B 95, 195144 (2017).
  * Larkin  _et al._ (2018) T. I. Larkin, R. D. Dawson, M. Höppner, T. Takayama, M. Isobe, Y.-L. Mathis, H. Takagi, B. Keimer, and A. V. Boris, Phys. Rev. B 98, 125113 (2018).
  * Seki  _et al._ (2014) K. Seki, Y. Wakisaka, T. Kaneko, T. Toriyama, T. Konishi, T. Sudayama, N. Saini, M. Arita, H. Namatame, M. Taniguchi, _et al._ , Physical Review B 90, 155116 (2014).
  * Watson  _et al._ (2020) M. D. Watson, I. Marković, E. A. Morales, P. L. Fèvre, M. Merz, A. A. Haghighirad, and P. D. C. King, Phys. Rev. Research 2, 013236 (2020).
  * Mor  _et al._ (2017) S. Mor, M. Herzog, D. Golež, P. Werner, M. Eckstein, N. Katayama, M. Nohara, H. Takagi, T. Mizokawa, C. Monney, and J. Stähler, Phys. Rev. Lett. 119, 086401 (2017).
  * Okazaki  _et al._ (2018a) K. Okazaki, Y. Ogawa, T. Suzuki, T. Yamamoto, T. Someya, S. Michimae, M. Watanabe, Y. Lu, M. Nohara, H. Takagi, _et al._ , Nature communications 9, 4322 (2018a).
  * Baldini  _et al._ (2020) E. Baldini, A. Zong, D. Choi, C. Lee, M. H. Michael, L. Windgaetter, I. I. Mazin, S. Latini, D. Azoury, B. Lv, _et al._ , arXiv preprint arXiv:2007.02909 (2020).
  * Tang  _et al._ (2020) T. Tang, H. Wang, S. Duan, Y. Yang, C. Huang, Y. Guo, D. Qian, and W. Zhang, Phys. Rev. B 101, 235148 (2020).
  * Suzuki  _et al._ (2020) T. Suzuki, Y. Shinohara, Y. Lu, M. Watanabe, J. Xu, K. L. Ishikawa, H. Takagi, M. Nohara, N. Katayama, H. Sawa, _et al._ , arXiv preprint arXiv:2002.10037 (2020).
  * Tanaka  _et al._ (2018) Y. Tanaka, M. Daira, and K. Yonemitsu, Physical Review B 97, 115105 (2018).
  * Tanabe  _et al._ (2018) T. Tanabe, K. Sugimoto, and Y. Ohta, Physical Review B 98, 235127 (2018).
  * Fujiuchi  _et al._ (2019) R. Fujiuchi, T. Kaneko, Y. Ohta, and S. Yunoki, Phys. Rev. B 100, 045121 (2019).
  * Golež  _et al._ (2016) D. Golež, P. Werner, and M. Eckstein, Phys. Rev. B 94, 035121 (2016).
  * Zong  _et al._ (2019a) A. Zong, P. E. Dolgirev, A. Kogar, E. Ergeçen, M. B. Yilmaz, Y.-Q. Bie, T. Rohwer, I.-C. Tung, J. Straquadine, X. Wang, Y. Yang, X. Shen, R. Li, J. Yang, S. Park, M. C. Hoffmann, B. K. Ofori-Okai, M. E. Kozina, H. Wen, X. Wang, I. R. Fisher, P. Jarillo-Herrero, and N. Gedik, Phys. Rev. Lett. 123, 097601 (2019a).
  * Tomeljak  _et al._ (2009) A. Tomeljak, H. Schäfer, D. Städter, M. Beyer, K. Biljakovic, and J. Demsar, Phys. Rev. Lett. 102, 066404 (2009).
  * Tsuji  _et al._ (2013) N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. 110, 136404 (2013).
  * Tsuji  _et al._ (2014) N. Tsuji, P. Barmettler, H. Aoki, and P. Werner, Phys. Rev. B 90, 075117 (2014).
  * Babadi  _et al._ (2015) M. Babadi, E. Demler, and M. Knap, Phys. Rev. X 5, 041005 (2015).
  * Babadi  _et al._ (2017) M. Babadi, M. Knap, I. Martin, G. Refael, and E. Demler, Phys. Rev. B 96, 014512 (2017).
  * Picano and Eckstein (2020) A. Picano and M. Eckstein, arXiv preprint arXiv:2009.04961 (2020).
  * Dolgirev  _et al._ (2020) P. E. Dolgirev, M. H. Michael, A. Zong, N. Gedik, and E. Demler, Phys. Rev. B 101, 174306 (2020).
  * Golež  _et al._ (2019) D. Golež, L. Boehnke, M. Eckstein, and P. Werner, Physical Review B 100, 041111 (2019).
  * Grazioli  _et al._ (2014) C. Grazioli, C. Callegari, A. Ciavardini, M. Coreno, F. Frassetto, D. Gauthier, D. Golob, R. Ivanov, A. Kivimäki, B. Mahieu, _et al._ , Review of Scientific Instruments 85, 023104 (2014).
  * Poletto  _et al._ (2009) L. Poletto, P. Villoresi, F. Frassetto, F. Calegari, F. Ferrari, M. Lucchini, G. Sansone, and M. Nisoli, Review of Scientific Instruments 80, 123109 (2009).
  * Okazaki  _et al._ (2018b) K. Okazaki, Y. Ogawa, T. Suzuki, T. Yamamoto, T. Someya, S. Michimae, M. Watanabe, Y. Lu, M. Nohara, H. Takagi, _et al._ , Nature communications 9, 1 (2018b).
  * Hellmann  _et al._ (2012) S. Hellmann, T. Rohwer, M. Kalläne, K. Hanff, C. Sohrt, A. Stange, A. Carr, M. Murnane, H. Kapteyn, L. Kipp, _et al._ , Nature communications 3, 1 (2012).
  * Mor  _et al._ (2018) S. Mor, M. Herzog, J. Noack, N. Katayama, M. Nohara, H. Takagi, A. Trunschke, T. Mizokawa, C. Monney, and J. Stähler, Phys. Rev. B 97, 115154 (2018).
  * Ligges  _et al._ (2018) M. Ligges, I. Avigo, D. Golež, H. U. R. Strand, Y. Beyazit, K. Hanff, F. Diekmann, L. Stojchevska, M. Kalläne, P. Zhou, K. Rossnagel, M. Eckstein, P. Werner, and U. Bovensiepen, Phys. Rev. Lett. 120, 166401 (2018).
  * Werdehausen  _et al._ (2018) D. Werdehausen, T. Takayama, M. Höppner, G. Albrecht, A. W. Rost, Y. Lu, D. Manske, H. Takagi, and S. Kaiser, Science advances 4, eaap8652 (2018).
  * Andrich  _et al._ (2020) P. Andrich, H. M. Bretscher, P. Telang, A. Singh, L. Harnaga, A. K. Sood, and A. Rao, arXiv preprint arXiv:2007.03368 (2020).
  * Murray  _et al._ (2007) E. D. Murray, S. Fahy, D. Prendergast, T. Ogitsu, D. M. Fritz, and D. A. Reis, Phys. Rev. B 75, 184301 (2007).
  * Hedayat  _et al._ (2019) H. Hedayat, C. J. Sayers, D. Bugini, C. Dallera, D. Wolverson, T. Batten, S. Karbassi, S. Friedemann, G. Cerullo, J. van Wezel, S. R. Clark, E. Carpene, and E. Da Como, Phys. Rev. Research 1, 023029 (2019).
  * Seki  _et al._ (2011) K. Seki, R. Eder, and Y. Ohta, Phys. Rev. B 84, 245106 (2011).
  * Balzer  _et al._ (2010) K. Balzer, S. Bauch, and M. Bonitz, Phys. Rev. A 82, 033427 (2010).
  * Schlünzen  _et al._ (2019) N. Schlünzen, S. Hermanns, M. Scharnke, and M. Bonitz, Journal of Physics: Condensed Matter 32, 103001 (2019).
  * Tuovinen  _et al._ (2019) R. Tuovinen, F. Covito, and M. A. Sentef, The Journal of chemical physics 151, 174110 (2019).
  * Tuovinen  _et al._ (2020) R. Tuovinen, D. Golež, M. Eckstein, and M. A. Sentef, Phys. Rev. B 102, 115157 (2020).
  * Freericks  _et al._ (2009) J. K. Freericks, H. R. Krishnamurthy, and T. Pruschke, Phys. Rev. Lett. 102, 136401 (2009).
  * Zenker  _et al._ (2012) B. Zenker, D. Ihle, F. X. Bronold, and H. Fehske, Phys. Rev. B 85, 121102 (2012).
  * (56) T. Kaneko, S. Ejima, H. Fehske, and Y. Ohta, “Excitonic bcs-bec crossover in double-layer systems,” in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2013)_, https://journals.jps.jp/doi/pdf/10.7566/JPSCP.3.017006 .
  * Mathias  _et al._ (2016) S. Mathias, S. Eich, J. Urbancic, S. Michael, A. Carr, S. Emmerich, A. Stange, T. Popmintchev, T. Rohwer, M. Wiesenmayer, _et al._ , Nature communications 7, 1 (2016).
  * Zong  _et al._ (2019b) A. Zong, A. Kogar, Y.-Q. Bie, T. Rohwer, C. Lee, E. Baldini, E. Ergeçen, M. B. Yilmaz, B. Freelon, E. J. Sie, _et al._ , Nature Physics 15, 27 (2019b).
  * Zenker  _et al._ (2014) B. Zenker, H. Fehske, and H. Beck, Phys. Rev. B 90, 195118 (2014).
  * Murakami  _et al._ (2020) Y. Murakami, D. Golež, T. Kaneko, A. Koga, A. J. Millis, and P. Werner, Phys. Rev. B 101, 195118 (2020).
  * Guseǐnov and Keldysh (1973) R. Guseǐnov and L. Keldysh, JETP 36, 1193 (1973).
  * Littlewood and Zhu (1996) P. Littlewood and X. Zhu, Physica Scripta 1996, 56 (1996).
  * Schüler  _et al._ (2020) M. Schüler, D. Golež, Y. Murakami, N. Bittner, A. Herrmann, H. U. Strand, P. Werner, and M. Eckstein, Computer Physics Communications 257, 107484 (2020).