# Fluctuating Nature of Light-Enhanced $d$-Wave Superconductivity:
A Time-Dependent Variational Non-Gaussian Exact Diagonalization Study

Yao Wang  Department of Physics and Astronomy, Clemson University, Clemson,
South Carolina 29631, USA  Tao Shi  CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  CAS Center for Excellence in Topological Quantum Computation,
University of Chinese Academy of Sciences, Beijing 100049, China  Cheng-Chien
Chen  Department of Physics, University of Alabama at Birmingham, Birmingham,
Alabama 35294, USA

###### Abstract

Engineering quantum phases using light is a novel route to designing
functional materials, where light-induced superconductivity is a successful
example. Although this phenomenon has been realized experimentally, especially
for the high-$T_{c}$ cuprates, the underlying mechanism remains mysterious.
Using the recently developed variational non-Gaussian exact diagonalization
method, we investigate a particular type of photoenhanced superconductivity by
suppressing a competing charge order in a strongly correlated electron-
electron and electron-phonon system. We find that the $d$-wave
superconductivity pairing correlation can be enhanced by a pulsed laser,
consistent with recent experiments based on gap characterizations. However, we
also find that the pairing correlation length is heavily suppressed by the
pump pulse, indicating that light-enhanced superconductivity may be of
fluctuating nature. Our findings also imply a general behavior of
nonequilibrium states with competing orders, beyond the description of a mean-
field framework.

##  I Introduction

Understanding, controlling, and designing functional quantum phases are major
goals and challenges in modern condensed matter physics Basov  _et al._
(2017); Wang  _et al._ (2018a). Among a few successful examples, light-induced
superconductivity, above its original transition temperature, has been a great
surprise and is believed to be a promising route to room-temperature
superconductors. Although this novel phenomenon has been realized
experimentally in various materials including cuprates, fullerides, and
organic salts Mankowsky  _et al._ (2014); Hu  _et al._ (2014); Kaiser  _et
al._ (2014); Mitrano  _et al._ (2016); Buzzi  _et al._ (2020), its mechanism
remains controversial. The underlying physics is partly attractive yet
mysterious for the enhanced $d$-wave superconductivity observed in pumped
charge-ordered high-$T_{c}$ cuprates La1.675Eu0.2Sr0.125CuO4 and La2-xBaxCuO4
near 1/8 doping Fausti  _et al._ (2011); Nicoletti  _et al._ (2014); Först
_et al._ (2014); Casandruc  _et al._ (2015); Khanna  _et al._ (2016);
Nicoletti  _et al._ (2018); Cremin  _et al._ (2019), partially due to their
original high critical temperatures at equilibrium. As illustrated in Fig.
1(a), the occurrence of the Cooper pairs above $T_{c}$ (reflected by the
Josephson plasma resonance in experiments), as a signature of light-induced
superconductivity, is observed when these materials are stimulated by a near-
infrared pulse laser.

Motivated by the experiments, various theories and numerical simulations have
been conducted to explain the observed light-induced phenomena. In the context
of conventional BCS superconductivity, simulations with both mean-field and
many-body models have demonstrated the feasibility to manipulate and enhance
local Cooper pairs (i.e. $s$-wave superconductivity) Sentef  _et al._ (2016,
2017); Coulthard  _et al._ (2017); Babadi  _et al._ (2017); Murakami  _et al._
(2017); Dasari and Eckstein (2018); Bittner  _et al._ (2019); Li and Eckstein
(2020); Tindall  _et al._ (2020a, b). The understanding becomes more
challenging for the unconventional $d$-wave superconductivity in cuprates, due
to the two-dimensional geometry and the strongly correlated nature. Insightful
theoretical perspectives have been proposed in the context of phenomenological
or steady-state theory, including the suppression of competing charge order
Patel and Eberlein (2016), Floquet engineering of the Fermi surface Kennes
_et al._ (2019), and (for the terahertz pump) parametric amplification Michael
_et al._ (2020). In contrast to the studies of $s$-wave superconductivity,
rigorous nonequilibrium simulations for nonlocal $d$-wave pairing instability
in microscopic models have been limited to undoped systems with truncated
phonon modes Wang  _et al._ (2018b), distinct from the conditions of existing
cuprate-based experiments. The extension to a doped system has been hindered
by the difficulty of treating both strong electronic correlation and electron-
phonon coupling in a quantum many-body simulation. In particular, the spatial
fluctuations of phonons and bosonic excitations in a 1/8-doped cuprates are
expected to be important due to the lack of a nesting momentum. Therefore, the
microscopic theory for light-induced superconductivity in a doped cuprate
remains an open question.

Figure 1:  Light-induced $d$-wave superconductivity experiment and theory. (a)
Sketch of existing light-induced or enhanced superconductivity experiments in
the charge-dominant phase in cuprates. (b) Summary of the simulation results
in this work. The green circles denote the enhancement of the pairing
correlation and the purple squares denote the suppression of the relative
correlation length, averaged between time $t=10-30\,t_{h}^{-1}$. The data are
extracted from Figs. 2 and 3. The left and right insets in (b) sketch real-
space distribution of pairing correlations, respectively, before and after
pump.

Moreover, recent experiments have revealed the presence of fluctuating
superconductivity above the transition temperature $T_{c}$ in overdoped
cuprates and FeSe Kondo  _et al._ (2015); He  _et al._ (2021); Faeth  _et al._
(2021); Xu  _et al._ (2021). In this regime, the superconductivity gap remains
open, but long-range order and zero resistance are absent due to the reduction
of correlation length. Recent ultrafast experiments also showed that the pump
light can destroy the coherence of Cooper pairs and fluctuate an existing
superconductor Boschini  _et al._ (2018); Yang  _et al._ (2018); Mootz  _et
al._ (2020). Thus, the resonance and gap in the transient optical conductivity
may cause a misassignment of the long-range superconductivity Chiriacò  _et
al._ (2018); Lemonik and Mitra (2019); Sun and Millis (2020). Both
observations raise the necessity to further investigate the coherence of
superconductivity induced by light in a quantum many-body model.

For this purpose, we study the photoinduced dynamics and superconductivity in
a light-driven Hubbard-Holstein model relevant to cuprates. To overcome the
numerical difficulties of simulating many-body dynamics with strong electronic
correlations and electron-phonon coupling, we construct on top of the recently
developed variational non-Gaussian exact diagonalization (NGSED) and develop
the time-dependent NGSED (see Sec. III). This hybrid method leverages the
merits of both the numerical many-body solver and the variational non-Gaussian
solver: The former is necessary to unbiasedly tackle strong electronic
correlations, while the latter avoids the phonon’s unbounded Hilbert-space
problem and has been demonstrated efficient in describing the ground and
excited states of systems with the Fröhlich-type electron-phonon coupling Shi
_et al._ (2018, 2020, 2019); Knörzer  _et al._ (2021). As an extension of the
equilibrium NGSED Wang  _et al._ (2019), this time-dependent method provides
an accurate description of far-from-equilibrium states through the Krylov-
subspace method and the Kählerization of the solvers. These advances of the
numerical method allow the simulation of light-induced dynamics in quantum
materials with both electronic correlations and electron-phonon couplings.

As summarized in Fig. 1(b), our results suggest that the $d$-wave pairing
correlation can be dramatically enhanced by a pulsed laser when the system is
charge dominant and close to a phase boundary, consistent with optical
experiments. However, we also find that the pairing correlation length is
heavily suppressed by the pump pulse. Therefore, light-enhanced $d$-wave
superconductivity may be of fluctuating nature. In addition to existing
ultrafast reflectivity measurements, our theoretical findings also predict
various observations verifiable by future transport and photoemission
experiments on pumped La1.675Eu0.2Sr0.125CuO4 and La2-xBaxCuO4.

The organization of this paper is as follows. We introduce our microscopic
model for the cuprate system and light-matter interaction in Sec. II. Next, we
briefly introduce the assumptions and framework of the time-dependent NGSED
method in Sec. III, while the detailed derivations and benchmarks are shown in
Appendix A and B. The main results about light-enhanced $d$-wave
superconductivity and the fluctuating nature are presented in Secs. IV and V.
We then discuss the frequency dependence of these observations in Sec. VI.
Finally, we conclude our paper and discuss relevant experimental predictions
in Sec. VII.

##  II Prototypical Models

Unconventional superconductivity is believed to emerge from the intertwined
orders of strongly correlated systems, where both spin and charge
instabilities exist Keimer  _et al._ (2015); Davis and Lee (2013). It was
widely believed that spin fluctuation, induced by electron correlation, is a
viable candidate to provide the pairing glue for $d$-wave superconductivity
Scalapino  _et al._ (1986); Gros  _et al._ (1987); Kotliar and Liu (1988);
Tsuei and Kirtley (2000). The minimal model to represent this strong
correlation is the single-band Hubbard model Zhang and Rice (1988). Based on
rigorous numerical simulations, many important experimental discoveries have
been reproduced using this model, such as antiferromagnetism Singh and Goswami
(2002), stripe phases Zheng  _et al._ (2017); Huang  _et al._ (2017, 2018);
Ponsioen  _et al._ (2019), strange metallicity Kokalj (2017); Huang  _et al._
(2019); Cha  _et al._ (2020), and superconductivity Zheng and Chan (2016); Ido
_et al._ (2018); Jiang and Devereaux (2019). However, increasing experimental
evidence reveals the significant role of phonons in high-$T_{c}$
superconductors Shen  _et al._ (2004); Lanzara  _et al._ (2001); Reznik  _et
al._ (2006); Devereaux and Hackl (2007); He  _et al._ (2018), in addition to
the strong electronic correlation. Together with the crucial lattice effects
observed in pump-probe experiments Mankowsky  _et al._ (2014); Hu  _et al._
(2014); Kaiser  _et al._ (2014), we believe that the minimal model to describe
the light-induced $d$-wave superconductivity must involve both interactions.

The Hubbard-Holstein model is the prototypical model for describing correlated
quantum materials with both electron-electron interaction and electron-phonon
coupling Hubbard (1963); Holstein (1959). Its Hamiltonian is written as

$\displaystyle{\mathcal{H}}$ $\displaystyle=$ $\displaystyle-
t_{h}\sum_{\langle\textbf{i},\textbf{j}\rangle,\sigma}\left(c^{\dagger}_{\textbf{i}\sigma}c_{\textbf{j}\sigma}+H.c.\right)+U\sum_{\textbf{i}}n_{\textbf{i}\uparrow}n_{\textbf{i}\downarrow}$
(1)
$\displaystyle+g\sum_{\textbf{i}\sigma}(a_{\textbf{i}}+a_{\textbf{i}}^{\dagger})n_{\textbf{i}\sigma}+\omega_{0}\sum_{\textbf{i}}a_{\textbf{i}}^{\dagger}a_{\textbf{i}}\,.$

Here, $c_{\textbf{i}\sigma}$ ($c_{\textbf{i}\sigma}^{\dagger}$) annihilates
(creates) an electron at site i with spin $\sigma$, and $a_{\textbf{i}}$
($a_{\textbf{i}}^{\dagger}$) annihilates (creates) a phonon at site i. To
reduce the number of model parameters, we restrict the hopping integral
$t_{h}$ to the nearest neighbor, electron-electron interaction to the on-site
Coulomb (Hubbard) repulsion $U$, electron-phonon coupling to the on-site
electrostatic coupling (Holstein) $g$, and phonon energy to the dispersionless
$\omega_{0}$. The electron-electron interaction and electron-phonon coupling
can be depicted by the dimensionless parameters $u=U/t_{h}$ and
$\lambda=g^{2}/\omega_{0}t_{h}$, respectively. Here, we set $\omega_{0}=t_{h}$
in accord with common choices Clay and Hardikar (2005); Nowadnick  _et al._
(2012, 2015); Johnston  _et al._ (2013); Wang  _et al._ (2016); Mendl  _et
al._ (2017). Throughout the paper, we focus on 12.5% hole doping simulated on
an $N=4\times 4$ cluster, corresponding to the La1.675Eu0.2Sr0.125CuO4 and
La1.875Ba0.125CuO4 experiments Fausti  _et al._ (2011); Nicoletti  _et al._
(2014); Först  _et al._ (2014); Casandruc  _et al._ (2015); Khanna  _et al._
(2016); Nicoletti  _et al._ (2018); Zhang  _et al._ (2018); Cremin  _et al._
(2019).

In an undoped system with a well-defined nesting momentum, one can restrict
the phonons to only the $\textbf{q}=(\pi,\pi)$ mode Wang  _et al._ (2016,
2018b). However, in a doped Hubbard-Holstein model without commensurability,
phonon modes at all momenta should be considered. As shown in Ref. Wang  _et
al._ , 2019, the phase diagram for a doped Hubbard-Holstein model is dominated
by a regime with strong charge susceptibility and a regime with strong spin
susceptibility, although both are not ordered like the half-filled case.
Unlike the striped order demonstrated in the Hubbard model at the
thermodynamic limit Zheng  _et al._ (2017); Huang  _et al._ (2017, 2018);
Ponsioen  _et al._ (2019), our system size does not support such a period-8
instability. Therefore, to increase the charge correlation and mimic the
charge-dominant cuprate system, we exploit relatively strong electron-phonon
couplings.

The light-driven physics is described (on the microscopic level) through the
Peierls substitution
$t_{h}c^{\dagger}_{\textbf{i}\sigma}c_{\textbf{j}\sigma}\rightarrow
t_{h}e^{i\int_{\textbf{r}_{\textbf{i}}}^{\textbf{r}_{\textbf{j}}}\mathbf{A}(t)\cdot
d\textbf{r}}c^{\dagger}_{\textbf{i}\sigma}c_{\textbf{j}\sigma}$, where the
vector potential $\mathbf{A}(t)$ of the external light pulse affects the many-
body Hamiltonian Eq. (1). In this paper, we simulate the pump pulse with an
oscillatory Gaussian vector potential:

$\mathbf{A}(t)=A_{0}\hat{e}_{\rm
pol}\exp\left[-\frac{(t-t_{0})^{2}}{2\sigma^{2}}\right]\cos(\Omega t)\,,$ (2)

and fix the polarization as diagonal $\hat{e}_{\rm
pol}=\frac{1}{\sqrt{2}}(\hat{x}+\hat{y})$ and the pump frequency as
$\Omega=4t_{h}$ (close to the 800-nm laser for $t_{h}=350$meV). In a strongly
correlated model like the Hubbard model or the extended Hubbard model, the
ultrafast pump pulse is coupled to electrons and may manipulate the delicate
balance between different competing orders Lu  _et al._ (2012); Al-Hassanieh
_et al._ (2013); Rincón  _et al._ (2014); Lu  _et al._ (2015). When a finite-
frequency phonon is involved, the photomanipulated competition of phases also
relies on the retardation effect of the phonons.

Throughout this paper, we restrict ourselves to the two-dimensional Hubbard-
Holstein model and transient equal-time correlation functions. The relation
between the in-plane Cooper pairs and the $c$-axis Josephson plasmon resonance
has been studied through semiclassical simulations Denny  _et al._ (2015);
Raines  _et al._ (2015); Okamoto  _et al._ (2017). Our simulation does not
consider specific experimental conditions including the aforementioned probe
schemes, material-specific matrix elements, and finite temperature,
nevertheless it helps explain the observed phenomena, in general. Therefore,
our correlation function analysis in Secs. IV-VI aims to address a matter-of-
principle question when the balance between charge-density wave (CDW) and
superconductivity, in a strongly correlated material, is altered by a pulsed
laser.

##  III Time-Dependent Variational Non-Gaussian Exact Diagonalization Method

The simulation of nonequilibrium quantum many-body systems requires either
Green’s function or wavefunction methods. The Green’s function methods,
constructed on the Keldysh formalism and represented by the dynamical mean-
field theory Freericks  _et al._ (2006); Aoki  _et al._ (2014), have achieved
great success in solving correlated materials and pump-probe spectroscopies in
the thermodynamic limit Werner and Millis (2007); Eckstein  _et al._ (2010).
However, when multiple instabilities compete in systems with strong electronic
correlations and electron-phonon couplings, the accuracy cannot be guaranteed
through perturbations which may lead to a biased solution. On the other hand,
wavefunction methods, such as exact diagonalization (ED) and density-matrix
renormalization group, have well-controlled numerical error, but are
restricted to small systems or low dimensions Dagotto (1994); White (1992).
This issue becomes even more severe when the electron-phonon coupling is non-
negligible due to the unbounded phonon Hilbert space Zhang  _et al._ (1998);
Brockt  _et al._ (2015).

To tackle the strong and dynamical electron-phonon coupling and overcome the
issue of phonon Hilbert space, we develop a time-dependent extension of the
NGSED method. The idea of this method is based on the following observations
for electrons and phonons, respectively. Electrons have a complicated form of
interactions (four-fermion terms) and intertwined instabilities, which thereby
have to be treated by an accurate solver. However, with the Pauli exclusion
principle, the electronic Hilbert-space dimension is relatively small, which
allows an exact solution for finite-size or low-dimensional systems. In
contrast, the phonon-phonon interaction (anharmonicity) is usually weak and
electron-phonon coupling has a linear form, with the complexity coming instead
from the unbounded local phonon Hilbert space. Therefore, if one can find an
entangler transformation involving a general form of entanglements between
electrons and phonons, the solution of the correlated Hamiltonian can be
mapped (after the transformation) to a factorizable wavefunction as a product
state consisting of electrons and phonons separately. Using this trick, we can
take advantage of the distinct properties mentioned above and solve the two
subsystems using different techniques.

In equilibrium, such an efficient entangler transformation has been found and
benchmarked, in terms of the generalized polaron transformation Shi  _et al._
(2018). This leads to the wavefunction ansatz Wang  _et al._ (2019)

$\displaystyle\big{|}\Psi_{\rm G}\big{\rangle}\\!=\\!U_{\rm plrn}|\psi_{\rm
ph}\rangle\\!\otimes\\!|\psi_{\rm
e}\rangle\\!=\\!e^{\frac{i}{\sqrt{N}}\sum_{\textbf{q}}\lambda_{\textbf{q}}p_{-\textbf{q}}\rho_{\textbf{q}}}|\psi_{\rm
ph}\rangle\\!\otimes\\!|\psi_{\rm e}\rangle\,,$ (3)

with the electron density operator
$\rho_{\textbf{q}}=\sum_{\textbf{i}\sigma}n_{\textbf{i}\sigma}e^{-i\textbf{q}\cdot\textbf{r}_{\textbf{i}}}$,
the phonon momentum
$p_{\textbf{q}}=i\sum_{\textbf{i}}(a_{\textbf{i}}^{\dagger}-a_{\textbf{i}})e^{-i\textbf{q}\cdot\textbf{r}_{\textbf{i}}}/\sqrt{N}$,
and the phonon displacement
$x_{\textbf{q}}=\sum_{\textbf{i}}(a_{\textbf{i}}+a_{\textbf{i}}^{\dagger})e^{-i\textbf{q}\cdot\textbf{r}_{\textbf{i}}}/\sqrt{N}$,
where $N$ is the number of lattice sites in the calculation. Here, the right-
hand side is a direct product of electron and phonon states: The electronic
wavefunction $|\psi_{\rm e}\rangle$ is treated as a full many-body state,
while the phonon wavefunction $|\psi_{\rm ph}\rangle$ is treated as coherent
Gaussian state (see Appendix. A for details). The entangler transformation
involves momentum-dependent variational parameters $\lambda_{\textbf{q}}$,
describing the polaronic dressing. Physically, a larger dressing
$\lambda_{\textbf{q}}\sim g/\omega_{0}$ accurately describes the phonon and
coupling energies, known as the Lang-Firsov transformation Lang and Firsov
(1962), while a smaller dressing allows a precise solution for the electron
energy. Thus, the variational parameters $\lambda_{\textbf{q}}$s are optimized
numerically as a balance between these two effects. Note that all
$\lambda_{\textbf{q}}$s are independent real numbers and entangle the phonon
momentum with electron density. This entangler transformation is demonstrated
to be sufficient in equilibrium Shi  _et al._ (2018); Knörzer  _et al._
(2021). For systems with both strong electron-phonon coupling and electronic
interactions, the electronic part of wavefunction $|\psi_{\rm e}\rangle$ can
be solved by ED and the above framework becomes the NGSED method. The
application of the NGSED to the equilibrium 2D Hubbard-Holstein model
successfully reveals the novel intermediate phases with superconducting
instability Wang  _et al._ (2019); Karakuzu  _et al._ (2017); Ohgoe and Imada
(2017), which is consistent with the recent DQMC-AFQMC simulations at the
thermodynamic limit Costa  _et al._ (2020).

In the real-time evolution, the electron density is driven by the pump laser,
the consequence of which is the varying force acting on the lattice and the
finite momentum of phonon. To characterize the electron-density-dependent
phonon momentum, we introduce an extra cubic coupling between the position and
the phonon density in the non-Gaussian transformation

$U_{\rm
NGS}(t)=e^{i\frac{1}{\sqrt{N}}\sum_{\textbf{q}}\lambda_{\textbf{q}}(t)(p_{-\textbf{q}}\cos\varphi_{\textbf{q}}(t)\,-x_{-\textbf{q}}\sin\varphi_{\textbf{q}}(t))\rho_{\textbf{q}}}$
(4)

to construct the variational ansatz

$\big{|}\Psi(t)\big{\rangle}=U_{\rm NGS}(t)|\psi_{\rm
ph}(t)\rangle\otimes|\psi_{\rm e}(t)\rangle\,.$ (5)

Here, the additional phase parameter $\varphi_{\textbf{q}}$ controls the ratio
of position and momentum displacements. The time dependence of
$\lambda_{\textbf{q}}$ and $\varphi_{\textbf{q}}$ allows dynamical fluctuation
of the polaronic dressing effect. At the same time, these two sets of
parameters, as a whole, disentangle the electrons and phonons with the
Fröhlich-type coupling and allow a relatively accurate description of the
transformed wavefunction via the direct-product state in the rightmost of Eq.
(5). The price one pays for the transformation is the complication of the
effective electronic Hamiltonian [see Eq. (19) in Appendix A], which is solved
by the time-dependent ED with well-controlled numerical errors [see
discussions in Appendix B]. As a minimal extension to the equilibrium
wavefunction ansatz Eq. (3), the structure of Eq. (4) involves an implicit
Kählerization, which guarantees the minimization of errors throughout the
evolution (see discussions in Appendix A). In each numerically small step of
the time evolution, i.e. $|\Psi(t)\rangle\rightarrow|\Psi(t+\delta t)\rangle$,
we simultaneously evolve both the variational parameters and the full
electronic wavefunction $|\psi_{e}(t)\rangle$, as explained in Appendixes A
and B, respectively.

Figure 2:  Light-driven dynamics of charge, spin, and pairing correlations.
(a). The square of local moment $\langle m_{z}^{2}\rangle$ for different model
parameters, subtracting off that for noninteracting electrons $\langle
m_{z}^{2}\rangle_{0}$. (b). The square of local moment $\langle
m_{z}^{2}\rangle$ (red) and the double occupation $\langle n_{d}\rangle$ as a
function of $u$ for a fixed $\lambda=4$, where $\langle m_{z}^{2}\rangle_{0}$
is indicated by a horizontal dotted line. The two parameter sets used on right
are marked as the dots in (a) and the dashed lines in (b). (c1-c2) The
evolution of charge correlation $N(\pi,\pi)$ after a linear pump of various
strengths $A_{0}$ ranging from 0.1 to 1 (denoted by different colors) with
$\lambda=4$ and (c1) $u=6.6$ and (c2) $u=7$, respectively. (d1-e2) The same as
(c1-c2) but for (d1-d2) the spin correlation $S(\pi,\pi)$ and (e1-e2) the
$d$-wave pairing susceptibility $P^{(d)}$. The pump pulse is sketched above.

##  IV Light-Enhanced Pairing Correlations

Figure 2(a) shows the square of local moment $\langle
m_{z}^{2}\rangle=\sum_{i}\langle(n_{i\uparrow}-n_{i\downarrow})^{2}\rangle/N$
calculated for different $u$ and $\lambda$ of the 12.5% doped Hubbard-Holstein
model, where $N$ is the system size. $\langle m_{z}^{2}\rangle$ provides an
indication of the overall spin moment and spin correlation. An abrupt change
of $\langle m_{z}^{2}\rangle$ with varying $u$ and/or $\lambda$ can signal a
phase transition. At the large-$u$ regime, the system is dominated by the spin
correlations inherent from the Hubbard model. In this regime, $d$-wave
superconductivity is identified in a quasi-1D system Jiang and Devereaux
(2019); Chung  _et al._ (2020), although inthe 2D thermodynamic limit its
existence is not fully established Maier  _et al._ (2005). In contrast, at the
large-$\lambda$ regime, the electrons are bound with phonons as bipolarons,
exhibiting a local moment lower than that of noninteracting electrons. This is
also reflected in the average double occupation $\langle
n_{d}\rangle=\sum_{i}\langle n_{i\uparrow}n_{i\downarrow}\rangle/N$ presented
in Fig. 2(b). There is an intermediate regime between these two limits, where
the charge correlation slightly dominates (reflected by the comparison to non-
interacting-electron local moment $\langle m_{z}^{2}\rangle_{0}$), while the
spin correlations remain finite. Based on recent studies using NGSED and QMC
methods, this intermediate regime persists at the thermodynamics limit and
reflects the realistic phases in a competing-order system Wang  _et al._
(2019); Costa  _et al._ (2020). We focus on this intermediate regime, which we
believe corresponds to the situation in cuprate experiments Fausti  _et al._
(2011); Nicoletti  _et al._ (2014); Först  _et al._ (2014); Casandruc  _et
al._ (2015); Khanna  _et al._ (2016); Nicoletti  _et al._ (2018); Cremin  _et
al._ (2019); Wang  _et al._ (2018b).

Figure 3:  Time-dependent correlation length and spatial distribution of
Cooper pairs. (a,b) Evolution of the correlation length $\xi$ of $d$-wave
Cooper pairs renormalized by its equilibrium value $\xi_{0}$, obtained for (a)
$u=6.6$ and (b) $u=7$, respectively. (c) The real-space distribution of
$\langle\Delta_{\textbf{r}}^{(d)\dagger}\Delta_{\bf 0}^{(d)}\rangle$ for
different r (considering translational invariance) evaluated at equilibrium,
$t=-3t_{h}^{-1}$, 0, and $3t_{h}^{-1}$, respectively. The upper and lower
present the dynamics induced by two different pump strengths $A_{0}=0.5$ and
$A_{0}=1$, respectively. The simulations are obtained for the $u=6.6$ and
$\lambda=4$ system. (d) The same as (c) but for the $u=7$ and $\lambda=4$
system.

We select two sets of model parameters inside this intermediate regime. Set 1
($u=6.6$ and $\lambda=4$) lies in the center of the intermediate regime, while
Set 2 ($u=7$ and $\lambda=4$) lies close to the boundary. We first focus on
Set 1 and examine the light-induced dynamics for the charge structure factor
$N(\textbf{q})=\langle\rho_{-\textbf{q}}\rho_{\textbf{q}}\rangle/N$ and spin
structure factor
$S(\textbf{q})=\langle\rho^{s}_{-\textbf{q}}\rho^{s}_{\textbf{q}}\rangle/N$,
where
$\rho^{s}_{\textbf{q}}=\sum_{\textbf{i}}(n_{\textbf{i}\uparrow}-n_{\textbf{i}\downarrow})e^{-i\textbf{q}\cdot\textbf{r}_{\textbf{i}}}$.
We focus on the momentum $\textbf{q}=(\pi,\pi)$, due to the important role of
spin fluctuations at this momentum on superconductivity; data at other momenta
are presented in the Appendix C. Figures 2 (c1) and (d1) show the pump
dynamics of $N({\pi,\pi})$ and $S({\pi,\pi})$ for Set 1. The same as the half-
filled case Wang  _et al._ (2018b), the pump pulse suppresses the charge
correlations and enhances spin fluctuations, which possibly serve as the
pairing glue, at short time ($t<0$). This leads to the rise of $d$-wave
superconductivity instability [see Fig. 2(e1)], characterized by the pairing
susceptibility

$\displaystyle
P^{(d)}=\frac{1}{N}\left\langle\Delta^{(d)\dagger}_{0}\Delta^{(d)}_{0}\right\rangle,$
(6)

where the $d$-wave pairing operator reads

$\displaystyle\Delta^{(d)}_{\textbf{k}}$ $\displaystyle=$
$\displaystyle\sum_{\textbf{k}^{\prime}}c_{\textbf{k}-\textbf{k}^{\prime}\downarrow}c_{\textbf{k}^{\prime}\uparrow}\left[\cos
k^{\prime}_{x}-\cos k^{\prime}_{y}\right]\,.$ (7)

As Set 1 is relatively far from the spin-dominant regime, spin fluctuation is
unstable, manifest as the drop of $S({\pi,\pi})$ for strong pump ($A_{0}>0.5$)
at longer timescales. As a consequence, the $d$-wave pairing susceptibility
stops increasing and starts to decrease for these strong pump conditions. This
nonmonotonic behavior reflects that the transient $d$-wave pairing emerges
from the dedicated balance among charge, spin, and electronic itineracy. While
the latter two can be induced by light and enhance the $d$-wave pairing, their
internal competition may reduce the enhancement if the pump strength keeps
increasing. Noticeably, when the nonequilibrium state finally melts into such
a “metallic” state, lattice distortions are released. In these cases, the
dynamics exhibit a period $\sim 2\pi/\omega_{0}$. This is reflected in the
corresponding Fourier spectra (see also the Appendixes).

The situation becomes different if one drives the system with Set 2 ($u=7$ and
$\lambda=4$), where the system resides in proximity to a phase boundary.
Because of the existing strong magnetic fluctuations, we find that the
$d$-wave pairing correlation is easier to enhance. As shown in Figs. 2(c2-e2),
the pump pulse increasingly enhances $S({\pi,\pi})$ and $P^{\rm(d)}$ after
suppressing the charge fluctuations. Although the enhancement of
$S({\pi,\pi})$ saturates at large pump strengths, spin fluctuations are more
robust against the increase of pump, due to proximity to a quantum phase
boundary. Therefore, they never drop below the equilibrium values. This change
of spin fluctuations causes the saturation of $d$-wave pairing correlations in
Fig. 2(e2). Similarly, this enhanced $d$-wave pairing concurs with a melting
of the competing CDW order, which leads to an oscillation of phonon energy.
Since the CDW order is already very weak near the phase boundary, the
intensity of the amplitude mode is no longer visible in the dynamics.

##  V Quantum Fluctuations and Correlation Length

The increase of pairing correlation reflects the formation of Cooper pairs
induced by the pump pulse. The experimental reflection of this correlation is
the opening of a gap (or Josephson plasma) characterized by ultrafast optical
reflectivity Fausti  _et al._ (2011); Nicoletti  _et al._ (2014); Först  _et
al._ (2014); Casandruc  _et al._ (2015); Khanna  _et al._ (2016); Nicoletti
_et al._ (2018); Cremin  _et al._ (2019). However, we emphasize that the
opening of a gap does not necessarily reflect the onset of superconductivity.
As recently shown in cuprates and FeSe experiments, above $T_{c}$ there exists
a fluctuating superconductivity phase that exhibits a superconducting gap but
also a finite resistance, reflecting preformed Cooper pairs Kondo  _et al._
(2015); He  _et al._ (2021); Faeth  _et al._ (2021); Xu  _et al._ (2021). This
is a signature of strong correlations in quantum material distinct from the
mean-field notions, due to strong thermal or quantum fluctuations. For the
nonequilibrium superconductivity in this paper, we focus on the zero-
temperature dynamics and fluctuations may arise from quantum instead of
thermal origins.

Figure 4:  Frequency dependence of the light-induced dynamics and
superconductivity. (a-c) Evolution of the (a) charge correlation $N(\pi,\pi)$,
(b) the $d$-wave pairing susceptibility $P^{(d)}$, and (c) the correlation
length $\xi$, at fixed pump strength $A_{0}=0.8$ but with different pump
frequencies $\Omega$. (d,e) The postpump (d) pairing correlation $P^{(d)}$ and
(e) correlation length $\xi$ averaged between time $t=10-30\,t_{h}^{-1}$ for
different pump frequencies, relative to their equilibrium values.

To better clarify the nature of this photoinduced many-body state, we further
consider the spatial fluctuation through the FFLO pairing correlation with
finite momentum
$P^{(d)}_{\textbf{k}}=\frac{1}{N}\left\langle\Delta^{(d)\dagger}_{\textbf{k}}\Delta^{(d)}_{\textbf{k}}\right\rangle$
Fulde and Ferrell (1964); Larkin and Ovchinnikov (1965). For $\textbf{k}=0$,
it recovers to the BCS pairing correlation, reflecting the total number of
Cooper pairs. The correlation length $\xi$ can be estimated through

$\xi=\lim_{|\textbf{k}|\rightarrow
0}\frac{1}{|\textbf{k}|}\sqrt{\frac{P^{(d)}_{0}}{P^{(d)}_{\textbf{k}}}-1}\,,$
(8)

based on the definition of correlation length satisfying
$\langle\Delta_{\textbf{i}}^{(d)\dagger}\Delta_{\textbf{j}}^{(d)}\rangle\sim
e^{-|\textbf{r}_{i}-\textbf{r}_{j}|/\xi}$. Figures 3(a) and (b) present the
evolution of $\xi$, estimated using the smallest momentum accessible in our
cluster $\textbf{k}=(\pi/2,0)$ and normalized by the equilibrium correlation
length $\xi_{0}$. For the $u=6.6$ system where charge order is robust [3 (a)],
the correlation length is dramatically suppressed for $A_{0}<0.5$, where the
$P^{(d)}$ is enhanced in Fig. 2(e1). When the pump strength increases to
$A_{0}\geq 0.5$, $\xi$ even drops to zero at the center of the pump,
indicating the destruction of local Cooper pairs and consistent with the drop
of $P^{(d)}$ in Fig. 2(e1). This overall drop of correlation length can be
visualized in the real-space distribution of
$\langle\Delta_{\textbf{r}}^{(d)\dagger}\Delta_{0}^{(d)}\rangle$ shown in Fig.
3(c): The Cooper pairs become strongly localized and spatially decoherent
after the pump.

What is more interesting is the system near the phase boundary ($u=7$), where
the BCS pairing correlation $P^{(d)}$ always increases after the pump, as
discussed in Fig. 2(e2). However, this enhancement of total Cooper pairs is
always accompanied by a drop of the correlation length, as shown in Fig. 3(b).
This means that the photoinduced Cooper pairing is more fluctuating than the
equilibrium one. Different from the $u=6.6$ system, here the decrease of $\xi$
saturates at a moderate value, about half of the equilibrium $\xi_{0}$. To
filter out the oscillatory dynamics, we extract the average enhancement of
pairing correlation $P_{d}(t)/P_{d}(-\infty)$ and the suppression of
$\xi(t)/\xi(-\infty)$ during the time window $10t_{h}^{-1}<t<30t_{h}^{-1}$,
long after the pump pulse disappears. As already shown in Fig. 1(b), the
relative changes (suppression or enhancement) in these two quantities are
comparable for all pump strengths. Considering that long-range correlations
asymptotically approach
$\langle\Delta_{\textbf{i}}^{(d)\dagger}\Delta_{\textbf{j}}^{(d)}\rangle\sim
P^{(d)}e^{-|\textbf{r}_{i}-\textbf{r}_{j}|/\xi}$, our simulation reflects an
enhancement of short-range superconductivity but a suppression of the long-
range one def .

The difference between a short-range enhancement and long-range suppression
can be visualized via the spatial distribution obtained by the Fourier
transform of $P^{(d)}_{\textbf{k}}$. As shown in Fig. 3(c), the local pairing
correlation increases and the range where the correlation is visible remains
finite after the pump. Both facts suggest that this type of short-range
superconductivity is visible in pump-probe spectroscopies, albeit with a
finite broadening. Thus, the fluctuating nature of light-enhanced
superconductivity does not conflict with existing optical experiments in
La1.675Eu0.2Sr0.125CuO4 and La1.875Ba0.125CuO4 Fausti  _et al._ (2011);
Nicoletti  _et al._ (2014); Först  _et al._ (2014); Casandruc  _et al._
(2015); Khanna  _et al._ (2016); Nicoletti  _et al._ (2018); Cremin  _et al._
(2019). Such a suppression of the coherence is associated with the nonlocal
nature of $d$-wave Cooper pairs, which is in contrast to the light-enhanced
strength and coherence of local pairing Sentef  _et al._ (2017); Coulthard
_et al._ (2017); Li and Eckstein (2020); Tindall  _et al._ (2020a). This
anticorrelation between the strength and coherence is also distinct from
previous light-induced competing orders Lu  _et al._ (2012); Al-Hassanieh  _et
al._ (2013); Rincón  _et al._ (2014); Lu  _et al._ (2015), possibly resulting
from the fact that superconductivity arises as a third instability emergent
from the competition.

##  VI Frequency Dependence

We further investigate the impact of the pump frequencies on superconductivity
and the competing orders. With a fixed pump amplitude $A_{0}=0.8$, Figs.
4(a-c) present the dynamics of various quantities discussed above for
different frequencies. As $\Omega$ increases from $t_{h}$ (350 meV) to
$10t_{h}$ (3.5 eV), the suppression of the charge correlations is weaker. This
can be attributed to the loss of resonances to the phonon dynamics, which
plays the crucial role in the melting of the charge order. Interestingly,
these suppressions of CDW do not always translate into a monotonic increase of
the $d$-wave superconductivity. As shown in Fig. 4(b), the low-frequency pump
with $\Omega=2t_{h}$ leads to a slightly weaker enhancement and a dramatic
decrease of the correlation length $\xi$. The transient correlation length is
suppressed to zero in the center of the pump [see Fig. 4(c)], indicating the
loss of superconductivity. We tentatively attribute this anomaly to the fact
that the parametric driving condition $\Omega\sim 2\omega_{0}$ is reached Knap
_et al._ (2016); Babadi  _et al._ (2017), leading to strong fluctuations of
phonons and overwhelming the $d$-wave pairing Macridin  _et al._ (2006); Mendl
_et al._ (2017).

Excepting the two frequencies which correspond to parametric driving, the
enhancement of the $d$-wave pairing susceptibility drops monotonically with
the pump frequency, consistent with the melting of the charge correlations. We
stress that such a phenomenon reflects that the dynamics are driven by light-
induced quantum fluctuations instead of a thermal effect, as the fluence
increases (instead of decreases) with the pump frequency. Although some
nonequilibrium effects are simply attributed to the sudden heating of the
electronic states caused by the pump pulse, this is not the situation we found
in the light-enhanced $d$-wave superconductivity. More rigorous discrimination
of the thermal effects relies on the calculation of pump-probe spectroscopies
and a fitting with the thermal distribution function Wang  _et al._ (2017);
Shvaika  _et al._ (2018); Wang  _et al._ (2018c); Matveev  _et al._ (2019);
Mootz  _et al._ (2020), which is beyond the scope of this paper.

Figures 4(d) and (e) summarize the frequency dependence of the pairing
susceptibility and the correlation length, averaged between time
$t=10-30\,t_{h}^{-1}$ after the pump, similar to the pump strength dependence
shown in Fig. 1. Except for the situation close to phonon parametric
resonance, we find that the maximal superconductivity enhancement is reached
at approximately $1.4$ eV and rapidly decreases for higher frequencies. As
this frequency is close to the Mott gap resonance (approximately $4t_{h}\sim
1.4\,$eV) for $U=8t_{h}$, this enhancement can be attributed to the generation
of spin fluctuations across the Mott gap. This observation is consistent with
the wavelength-dependence study in LBCO, where the 800-nm laser is found to
enhance superconductivity move efficiently than a 400-nm laser Casandruc  _et
al._ (2015). Although other mechanisms beyond the single-band model have been
proposed to address this wavelength dependence, such a consistency strengthens
the connection between our microscopic simulation and the real experiments.

##  VII Summary and Outlook

Altogether, our study provides a novel perspective of interpreting the light-
induced $d$-wave pairing in the context of strong correlation and electron-
phonon coupling. Although it does not rule out other interpretations, our
result reflects that the light-induced Cooper pairs may be spatially local.
Therefore, the Josephson plasma resonance and the optical gap may coexist with
a small Drude weight, which reflects the superfluid density and cannot be
directly resolved in ultrafast reflectivity experiments. Fortunately, this
mystery was recently answered by a state-of-the-art ultrafast transport
measurement, in another type of light-induced superconductor (K3C60) Budden
_et al._ (2021). Such a transport measurement directly provides the transient
resistivity of the material, which reflects the superconductivity long-range
order. Therefore, a similar transport measurement in La1.875Ba0.125CuO4 will
help to distinguish fluctuations from ordered $d$-wave superconductivity.

In addition to optical and transport properties, the single-particle
electronic structure measurable by photoemission is employed frequently to
characterize the $d$-wave superconductors. In the single-particle context, the
equilibrium fluctuating superconductivity translates into the (single-
particle) gap opening above the transition temperature $T_{c}$ He  _et al._
(2021); Faeth  _et al._ (2021); Xu  _et al._ (2021), while the distinction
from long-range ordering may hide in the quantitative spectral shape (e.g.,
quasiparticle width and temperature dependence of the Fermi-surface spectral
weight) and is still an ongoing research. Identifying this fluctuation would
be even harder, but also promising, for a nonequilibrium state measured by the
time-resolved ARPES. The findings of our simulations indicate that the light-
driven La1.875Ba0.125CuO4 should exhibit a (single-particle) superconducting
gap, but its quasiparticle peak should be damped by the fluctuations.
Complementary to regular angle-resolved photoemission spectroscopy (ARPES),
future developments on two-electron ARPES measurement Su and Zhang (2020);
Mahmood  _et al._ (2021), at ultrafast timescales, are promising directions to
distinguish long-range superconductivity from a fluctuating one.

Alternatively, a recent attempt has investigated the inverse process by
measuring the coherence of CDW states using x-ray scattering Wandel  _et al._
(2020). Light-induced CDW also has been recently observed in LaTe3 with two
competing CDW orders Kogar  _et al._ (2020). Here, we investigate specifically
light-driven cuprates, but appropriate modification Demler and Zhang (1995);
Demler  _et al._ (1998) of our study and microscopic model can provide a
general platform to examine light-driven fluctuations for intrinsically
competing orders, as pointed out in a recent phenomenological study Sun and
Millis (2020). The fluctuating nature of the enhanced quantum phase suggests
that the nonequilibrium states of competing-order systems are beyond the
description of a mean-field framework.

To tackle the dynamics of systems with strong electronic correlation and
electron-phonon coupling, we develop a new technique by generalizing the
variational non-Gaussian exact diagonalization framework out of equilibrium.
This method overcomes the issues of unbounded phonon Hilbert space, while
keeping the numerical accuracy through the exact diagonalization of the
(transformed) electronic problem. Although we restrict the application into
the Hubbard-Holstein model relevant for the cuprate $d$-wave
superconductivity, our derivation is rather general and can be applied to
dispersive phonons, extended electronic interactions, and more complicated
electron-phonon coupling. The formalism also can be applied to multiple phonon
branches as long as the phonon-phonon interaction can be ignored, but the
generalization to multiband electrons is restricted by the Hilbert-space size
issue in exact diagonalization. Future application of this method also
includes the evolution of pump-probe spectroscopies Freericks  _et al._
(2009); Lenarčič  _et al._ (2014); Shimizu and Yuge (2011); Wang  _et al._
(2017, 2018c) by tracking the time-dependent gauges and evaluating the off-
diagonal overlap of Gaussian states.

## Acknowledgement

The authors thank J.I. Cirac, E. Demler, J.-F. He, M. Mitrano, M. A. Sentef,
and I. E. Perakis for insightful discussions. Y.W. acknowledges support from
National Science Foundation (NSF) Grant No. DMR-2038011. C.-C.C. acknowledges
support from NSF Grant No. OIA-1738698 and OAC-2031563. The calculations were
performed on the Frontera computing system at the Texas Advanced Computing
Center. Frontera is made possible by NSF NSF Grant No. OAC-1818253.

##  Appendix A Detailed Derivation of the Time-Dependent NGSED Method

To obtain the equations of motion for variational parameters in Eq. (5), where
the phonon state is

$\displaystyle|\psi_{\rm
ph}\rangle=e^{-\frac{1}{2}R_{0}^{T}\sigma_{y}\Delta_{R}}e^{-i\frac{1}{4}\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}\xi_{\textbf{q}}R_{\textbf{q}}}|0\rangle=U_{\rm
GS}|0\rangle\,,$ (9)

we notice that the generalization from $U_{\rm plrn}$ to $U_{\rm NGS}$
involves a rotation of phonon operators in the phase space. Therefore, the
$U_{\rm NGS}$ can be expressed by a composite transformation, leading to an
equivalent form of Eq. (5):

$\displaystyle\big{|}\Psi(t)\big{\rangle}=U_{\rm rot}(t)U_{\rm
plrn}(t)|\psi_{\rm ph}(t)\rangle\otimes|\psi_{\rm e}(t)\rangle\,,$ (10)

where the rotation is reflected in the unitary (Gaussian) transformation

$\displaystyle U_{\rm
rot}(t)=e^{-i\frac{1}{4}\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}\eta_{\textbf{q}}R_{\textbf{q}}},$
(11)

and $\eta_{\textbf{q}}$ is a variational matrix. In the derivation of
variational equations of motion, it is usually convenient to employ the
linearized form $\bar{S}_{\textbf{q}}$ to parametrize the rotational
transformation, i.e. $U_{\rm rot}^{\dagger}R_{\textbf{q}}U_{\rm
rot}=\bar{S}_{\textbf{q}}R_{\textbf{q}}$. Thus, ignoring the redundant degrees
of freedom, which can be absorbed into $|\psi_{\rm ph}\rangle$, we obtain the
following relation between $\varphi_{\textbf{q}}$ and $\bar{S}_{\textbf{q}}$

$\displaystyle\bar{S}_{\textbf{q}}(t)=e^{i\sigma_{y}\eta_{\textbf{q}}}=\left(\begin{array}[]{cc}\cos\varphi_{\textbf{q}}(t)&-\sin\varphi_{\textbf{q}}(t)\\\
\sin\varphi_{\textbf{q}}(t)&\cos\varphi_{\textbf{q}}(t)\end{array}\right).$
(14)

In the above wavefunction prototype, $\Delta_{R}$ (vector), $\xi_{\textbf{q}}$
(matrix), $\lambda_{\textbf{q}}$ and $\varphi_{\textbf{q}}$ (scalar) are
variational parameters, which are allowed to evolve during the dynamics.

To evaluate the nonequilibrium dynamics within the hybrid wavefunction ansatz,
we project the Schrödinger equation
$i\partial_{t}\big{|}\Psi(t)\big{\rangle}={\mathcal{H}}(t)\big{|}\Psi(t)\big{\rangle}$
onto the tangential space in the variational manifold, which gives equations
of motion for both the variational parameters and the electronic wavefunction
$|\psi_{e}(t)\rangle$. In contrast to an energy minimization (imaginary-time
evolution), the algebra of this projection in a real-time evolution may be ill
defined. For the general time-dependent variational principle, the real and
imaginary parts in the projection of Schrödinger equation on the tangential
space

$i\langle\xi_{\alpha}|\partial_{t}\big{|}{\Psi}(t)\big{\rangle}=\langle\xi_{\alpha}|{\mathcal{H}}(t)\big{|}{\Psi}(t)\big{\rangle},$
(15)

with $|\xi_{\alpha}\rangle$ denoting any tangential vector, may give rise to
two sets of incompatible equations of motion, which correspond to the Dirac-
Frankel and McLachlan variational principles, respectively Hackl  _et al._
(2020). If and only if the variational manifold is Kähler manifold, the two
principles result in the same equation of motion. By designing the
wavefunction ansatz in Eqs. (5) and (4), we Kählerize the variational manifold
such that the two principles are compatible and the derivations below are
self-consistent. Physically, the Kählerity guarantees that the variational
principle Eq. (15) minimizes errors with respect to the realistic time-
dependent wavefunctions.

Taking advantage of the composite representation of $U_{\rm NGS}$, it is
convenient to define

$\displaystyle\big{|}\bar{\Psi}(t)\big{\rangle}=U_{\rm
rot}^{\dagger}\big{|}\Psi(t)\big{\rangle}=U_{\rm plrn}(t)|\psi_{\rm
ph}\rangle\otimes|\psi_{\rm e}\rangle\,,$ (16)

and solve the Schrödinger equation in a rotating frame
$i\partial_{t}\big{|}\bar{\Psi}(t)\big{\rangle}=\bar{{\mathcal{H}}}\big{|}\bar{\Psi}(t)\big{\rangle}$.
Thus, the rotated Hamiltonian becomes

$\displaystyle\bar{{\mathcal{H}}}$ $\displaystyle=$ $\displaystyle U_{\rm
rot}^{\dagger}{\mathcal{H}}U_{\rm rot}-iU_{\rm
rot}^{\dagger}\partial_{t}U_{\rm rot}$ (17) $\displaystyle=$
$\displaystyle\sum_{\textbf{k}\sigma}(\varepsilon_{\textbf{k}}+\mu)c^{\dagger}_{\textbf{k}\sigma}c_{\textbf{k}\sigma}+\frac{1}{4}\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}(\omega_{\textbf{q}}+\partial_{t}\varphi_{\textbf{q}})R_{\textbf{q}}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{\textbf{q}}g_{\textbf{q}}R_{\textbf{q}}^{\dagger}\bar{S}_{\textbf{q}}^{\dagger}e_{1}\rho_{\textbf{q}}+{\mathcal{H}}_{e-e}.$

Here, we denote the general electron-electron interaction terms as
${\mathcal{H}}_{e-e}$ and assume it to commute with the electronic density
operator $n_{\textbf{i}}$, which is satisfied in the Hubbard-Holstein model
Eq. (1). This assumption substantially simplifies the derivation and leads to
the following Eq. (19).

On the left-hand side of the rotated Schrödinger equation, we obtain

$\displaystyle i\partial_{t}\big{|}\bar{\Psi}\big{\rangle}\\!$
$\displaystyle=$ $\displaystyle\\!U_{\rm plrn}U_{\rm
GS}\Big{[}-i\frac{1}{4}\Delta_{R}^{T}\sigma_{y}\partial_{t}\Delta_{R}-i\frac{1}{2}R_{0}^{T}S_{0}^{T}\sigma_{y}\partial_{t}\Delta_{R}$
(18)
$\displaystyle-i\frac{1}{4}\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}S_{\textbf{q}}^{\dagger}\sigma_{y}\partial_{t}S_{\textbf{q}}R_{\textbf{q}}-\frac{1}{\sqrt{N}}\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}S_{\textbf{q}}^{\dagger}e_{2}\rho_{\textbf{q}}\partial_{t}\lambda_{\textbf{q}}$
$\displaystyle-\frac{1}{\sqrt{N}}\Delta_{R}^{T}e_{2}\rho_{0}\partial_{t}\lambda_{0}\Big{]}|0\rangle_{\rm
ph}\otimes|\psi_{e}\rangle$ $\displaystyle+U_{\rm plrn}|\psi_{\rm
ph}\rangle\otimes(i\partial_{t}\big{|}\psi_{e}\big{\rangle})\,.$

Here, $S_{\textbf{q}}=e^{i\sigma_{y}\xi_{\textbf{q}}}$ is the linearized
transformation matrix of the $U_{\rm GS}$ and the effective Hamiltonian
becomes

$\displaystyle{\mathcal{H}}_{\rm eff}$ $\displaystyle=$ $\displaystyle U_{\rm
plrn}^{\dagger}\bar{{\mathcal{H}}}U_{\rm plrn}$ (19) $\displaystyle=$
$\displaystyle-
t_{h}\sum_{\textbf{j}\sigma\alpha\delta_{\alpha}}e^{i\frac{1}{\sqrt{N}}\sum_{\textbf{q}}\lambda_{\textbf{q}}R_{\textbf{q}}^{\dagger}e_{2}e^{-i\textbf{q}\cdot\textbf{j}}\left(1-e^{-i\textbf{q}\cdot\delta_{\alpha}}\right)}c_{\textbf{j}+\delta_{\alpha},\sigma}^{\dagger}c_{\textbf{j}\sigma}$
$\displaystyle-\frac{1}{2N}\sum_{\textbf{q}}\sum_{\textbf{k},\textbf{k}^{\prime}\atop\sigma,\sigma^{\prime}}V_{\textbf{q}}\,c_{\textbf{k}+\textbf{q},\sigma}^{\dagger}c_{\textbf{k}\sigma}c_{\textbf{k}^{\prime}-\textbf{q},\sigma^{\prime}}^{\dagger}c_{\textbf{k}^{\prime}\sigma^{\prime}}+{\mathcal{H}}_{e-e}$
$\displaystyle-\frac{1}{\sqrt{N}}\\!\sum_{\textbf{q}}R_{\textbf{q}}^{\dagger}\\!\left[\lambda_{\textbf{q}}\\!\left(\omega_{0}\\!+\\!\partial_{t}\varphi_{\textbf{q}}\right)\\!-\\!g\bar{S}_{\textbf{q}}^{\dagger}\right]\\!e_{1}\rho_{\textbf{q}}$
$\displaystyle+\frac{1}{4}\sum_{\textbf{q}}(\omega_{0}\\!+\\!\partial_{t}\varphi_{\textbf{q}})R_{\textbf{q}}^{\dagger}R_{\textbf{q}}\,,$

for the product state $|\psi_{\rm
ph}\rangle\otimes\big{|}\psi_{e}\big{\rangle}$. Here, $\alpha$ takes the $x$
and $y$ directions, while $\delta_{\alpha}$ denotes the vector pointing to the
nearest neighbors along the corresponding directions.

After transforming the Hamiltonian to the basis of product states, the
electron-phonon dressing effects are reflected in both the interaction and the
kinetic energy. For the effective electronic interaction, the dressing of
phonons mediates an additional $V_{\textbf{q}}$ on top of the original
$H_{e-e}$, whose variational expression is

$\displaystyle
V_{\textbf{q}}=4g\cos\varphi_{\textbf{q}}\mathrm{Re}[\lambda_{\textbf{q}}]-2|\lambda_{\textbf{q}}|^{2}(\omega_{0}+\partial_{t}\varphi_{\textbf{q}}).$
(20)

(Note, that we employed the assumption that $H_{e-e}$ commutes with
$n_{\textbf{i}}$ and, therefore, $U_{\rm plrn}$.) For the kinetic energy, the
phonon-dressed effective hopping term can be reformulated in a closed form

$\displaystyle
t_{\alpha}=t_{h}e^{-\sum_{\textbf{q}}\frac{|\lambda_{\textbf{q}}|^{2}}{N}(1-\cos
q_{\alpha})e_{2}^{T}\Gamma_{\textbf{q}}e_{2}}\,,$ (21)

when taking the expectation with respect to the phonon Gaussian state
$|\psi\rangle_{\rm ph}$.

Now with both the time derivative and the effective Hamiltonian transformed
into the product-state basis, we can project Eqs. (18) and (19) onto various
tangential vectors sequentially. For the zeroth-order tangential vector, which
is proportional to the phonon vacuum state $|0\rangle_{\rm ph}$, we obtain the
equation of motion for the electronic state:

$\displaystyle i\partial_{t}\big{|}\psi_{e}\big{\rangle}$ $\displaystyle=$
$\displaystyle\Big{[}\langle\psi_{\rm ph}|{\mathcal{H}}_{\rm eff}(t)|\psi_{\rm
ph}\rangle+i\frac{1}{4}\Delta_{R}^{T}\sigma_{y}\partial_{t}\Delta_{R}$
$\displaystyle+\frac{1}{\sqrt{N}}\Delta_{R}^{T}e_{2}\rho_{0}\partial_{t}\lambda_{0}-\frac{1}{2}\sum_{\textbf{q}}\partial_{t}\varphi_{\textbf{q}}\Big{]}|\psi_{e}\big{\rangle}.$

The (phonon-traced) ${\mathcal{H}}_{\rm eff}$ together with the scalar terms
in the square brackets govern the time evolution of the electronic
wavefunction $|\psi_{e}\rangle$. In the NGSED framework, we treat
$|\psi_{e}\rangle$ as a full many-body state and evaluate this evolution
stepwisely by the Krylov-subspace method (see Appendix B).

Moreover, the equations of motion for (the variational parameters of)
$|\psi_{\rm ph}\rangle$ can be obtained by projecting Eqs. (18) and (19) onto
the first-order tangential vector ($R_{0}^{T}S_{0}^{T}|0\rangle_{\rm
ph}\otimes|\psi_{e}\rangle$) and the second-order tangential vector
($R_{\textbf{q}}^{\dagger}R_{\textbf{q}}|0\rangle_{\rm
ph}\otimes|\psi_{e}\rangle$). After some algebra, these two equations are

$\displaystyle\partial_{t}\Delta_{R}$ $\displaystyle=$ $\displaystyle
i\sigma_{y}\Big{\\{}\left(\omega_{0}+\partial_{t}\varphi_{0}\right)\Delta_{R}-2\sqrt{N}\rho_{\rm
f}\Big{[}\left(\lambda_{0}\omega_{0}+\lambda_{0}\partial_{t}\varphi_{0}\right.$
(23)
$\displaystyle\left.-g_{0}\bar{S}_{0}^{\dagger}\right)e_{1}-\partial_{t}\lambda_{0}e_{2}\Big{]}\Big{\\}},$

and

$\displaystyle-
iS^{\dagger}_{\textbf{q}}\sigma_{y}\partial_{t}S_{\textbf{q}}=S_{\textbf{q}}^{\dagger}\tilde{\Omega}_{\textbf{q}}(t)S_{\textbf{q}},$
(24)

respectively. Here, $\rho_{\rm f}=N_{e}/N$ is the average electron density,
and the renormalized phonon energy matrix is

$\displaystyle\tilde{\Omega}_{\textbf{q}}(t)$ $\displaystyle=$
$\displaystyle\frac{8|\lambda_{\textbf{q}}|^{2}}{N}\sum_{k\alpha}t_{\alpha}\big{[}1-\cos
q_{\alpha}\big{]}\langle n_{\textbf{k}}\rangle\cos k_{\alpha}e_{2}e_{2}^{T}$
(25) $\displaystyle+\omega_{\textbf{q}}+\partial_{t}\varphi_{\textbf{q}}.$

As we are interested in the equal-time measurements (see Sec. IV), the gauge
degrees of freedom are redundant. Therefore, we employ a gauge-invariant
covariance matrix $\Gamma_{\textbf{q}}=S_{\textbf{q}}S_{\textbf{q}}^{\dagger}$
and track the evolution of $\Gamma_{\textbf{q}}$ instead of $S_{\textbf{q}}$.
Equation (24) becomes

$\displaystyle\partial_{t}\Gamma_{\textbf{q}}=i\sigma_{y}\tilde{\Omega}_{\textbf{q}}(t)\Gamma_{\textbf{q}}-i\Gamma_{\textbf{q}}\tilde{\Omega}_{\textbf{q}}(t)\sigma_{y}.$
(26)

Figure 5:  Time evolution of the charge structure factor (blue curve), spin
structure factor (red curve), $d$-wave pairing susceptibility (green curve),
and $s$-wave pairing susceptibility (orange curve), for a pump pulse with
$A_{0}=0.4$ and frequency $\Omega=5t_{h}$. The solid and dotted lines denote
the results obtained by NGSED and the standard ED (with maximal phonon number
$M=10$), respectively. The inset guides the eye for the pump pulse.

Finally, the equations of motion of $\lambda_{\textbf{q}}$ and
$\varphi_{\textbf{q}}$, which appear in $U_{\rm NGS}$ and entangle phonons and
electrons, can be obtained by projecting to the third-order tangential vector
$a_{\textbf{q}}^{\dagger}\rho_{\textbf{q}}|0\rangle_{\rm
ph}\otimes|\psi_{e}\rangle$. These two equations are

$\displaystyle\partial_{t}\lambda_{\textbf{q}}=g_{\textbf{q}}\sin\varphi_{\textbf{q}}+2\lambda_{\textbf{q}}\frac{\Pi_{\textbf{q}}}{C_{\textbf{q}}}e_{1}^{T}\Gamma_{\textbf{q}}e_{2},$
(27)

and

$\displaystyle\partial_{t}\varphi_{\textbf{q}}=\frac{g_{\textbf{q}}\cos\varphi_{\textbf{q}}}{\lambda_{\textbf{q}}}-\omega_{\textbf{q}}+2\frac{\Pi_{\textbf{q}}}{C_{\textbf{q}}}e_{2}^{T}\Gamma_{\textbf{q}}e_{2},$
(28)

respectively, where the modulated electronic correlation is

$\displaystyle\Pi_{\textbf{q}}=\sum_{k\sigma\alpha}\\!t_{\alpha}[\cos
k_{\alpha}\\!-\\!\cos(k_{\alpha}\\!+\\!q_{\alpha})]\Big{\langle}\\!{\rho_{-\textbf{q}}}{c}_{\textbf{k}\sigma}^{\dagger}{c}_{\textbf{k}\\!+\\!\textbf{q},\sigma}\\!\Big{\rangle}\,,$
(29)

and the density correlation is
$C_{\textbf{q}}=\langle{\rho_{-\textbf{q}}}{\rho_{\textbf{q}}}\rangle$.

With the variational wavefunction Eq. (5), we can express all (equal-time)
observables using analytical expressions formed by the variational parameters
and the correlations of fermionic operated evaluated in the electronic
wavefunction $|\psi_{e}\rangle$ Wang  _et al._ (2019). The most important
observable in this manuscript is $d$-wave pairing correlation, which can be
explicitly written as

$\displaystyle P^{(d)}_{\textbf{k}}\\!$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\textbf{k}^{\prime}}\sum_{\textbf{k}_{1},\textbf{k}_{2}}\xi(\textbf{k}^{\prime}+\textbf{k},\textbf{k}_{1},\textbf{k}_{2})\left[\zeta_{x}(\textbf{k}^{\prime})\cos\left(\frac{k^{\prime}_{x}}{2}-k_{1x}\right)\cos\left(\frac{k^{\prime}_{x}}{2}-k_{2x}\right)+\zeta_{y}(\textbf{k}^{\prime})\cos\left(\frac{k^{\prime}_{y}}{2}-k_{1y}\right)\cos\left(\frac{k^{\prime}_{y}}{2}-k_{2y}\right)\right.$
(30)
$\displaystyle\left.-\zeta_{xy}(\textbf{k}^{\prime})\cos\left(\frac{k^{\prime}_{x}}{2}-k_{2x}\right)\cos\left(\frac{k^{\prime}_{y}}{2}-k_{1y}\right)-\zeta_{xy}(\textbf{k}^{\prime})\cos\left(\frac{k^{\prime}_{y}}{2}-k_{2y}\right)\cos\left(\frac{k^{\prime}_{x}}{2}-k_{1x}\right)\right]\,.$

This expression contains the pairing correlation in the electronic part of
wavefunction

$\displaystyle\xi(\textbf{k}^{\prime},\textbf{k}_{1},\textbf{k}_{2})=\langle\psi_{e}|c_{\textbf{k}_{2}\downarrow}^{\dagger}c_{\textbf{k}^{\prime}-\textbf{k}_{2}\uparrow}^{\dagger}c_{\textbf{k}^{\prime}-\textbf{k}_{1}\uparrow}c_{\textbf{k}_{1}\downarrow}|\psi_{e}\rangle$
(31)

and analytical dressing factors

$\displaystyle\zeta_{x}(\textbf{k}^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\textbf{r}}e^{-i\textbf{k}^{\prime}\cdot\textbf{r}}e^{-\sum_{\textbf{q}}\left[(1+\cos
q_{x})(1-e^{i\textbf{q}\cdot\textbf{r}})\right]\frac{2|\lambda_{\textbf{q}}|^{2}}{N}e_{2}^{T}\Gamma_{\textbf{q}}e_{2}}$
$\displaystyle\zeta_{y}(\textbf{k}^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\textbf{r}}e^{-i\textbf{k}^{\prime}\cdot\textbf{r}}e^{-\sum_{\textbf{q}}\left[(1+\cos
q_{y})(1-e^{i\textbf{q}\cdot\textbf{r}})\right]\frac{2|\lambda_{\textbf{q}}|^{2}}{N}e_{2}^{T}\Gamma_{\textbf{q}}e_{2}}$
$\displaystyle\zeta_{xy}(\textbf{k}^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\bar{\textbf{r}}}e^{-i\textbf{k}^{\prime}\cdot\bar{\textbf{r}}}\exp\bigg{\\{}-\sum_{\textbf{q}}\Big{[}(2+\cos
q_{x}+\cos q_{y})$ (32)
$\displaystyle-4e^{i\textbf{q}\cdot\bar{\textbf{r}}}\cos\frac{q_{x}}{2}\cos\frac{q_{y}}{2}\Big{]}\frac{|\lambda_{\textbf{q}}|^{2}}{N}e_{2}^{T}\Gamma_{\textbf{q}}e_{2}\bigg{\\}}\,.$

Here, the $\bar{\textbf{r}}$ in the last summation denotes the half-unit-cell-
shifted coordinates $\bar{\textbf{r}}=\textbf{r}+\hat{x}/2-\hat{y}/2$.

Thus, we have now extended the NGSED framework into nonequilibrium dynamics.
Each step of the time evolution is achieved by the evaluation of above coupled
differential equations, as well as large-scale Krylov-subspace method (for the
electronic wavefunction). To benchmark the accuracy, we compare the simulation
results with the ED in a 2D Hubbard-Holstein model with $g=t_{h}$ and
$U=8t_{h}$. The ED simulation is exact, except for the finite phonon
occupation, which is truncated at $M=10$ in our calculation. To allow the ED
simulation with acceptable computational complexity and accuracy, we choose a
small eight-site cluster and a relatively high phonon frequency
$\omega_{0}=5t_{h}$. For a medium pump strength ($A_{0}=0.4$) with frequency
$\Omega=5t_{h}$, the relative errors for two dominant correlations, i.e., the
spin structure factor $S(\pi,\pi)$ and the $d$-wave pairing susceptibility
$P^{(d)}$ are both below 1% even in the center of the pump pulse (see Fig. 5).
At the same time, the relative error for the $s$-wave pairing susceptibility
$P^{(s)}$ reaches approximately $5\%$. This relatively larger deviation mainly
originates from the fact that the variational wavefunction tends to primarily
capture the dominant correlations with limited parameters, leaving a slightly
larger deviation for other correlations. In addition, the truncation error for
the phonon Hilbert space in ED simulations is also more prominent out of
equilibrium, which also contributes to the numerical error shown in Fig. 5.

##  Appendix B Krylov-Subspace Method for Time-Evolution

With instantaneous variational parameters, the Hamiltonian can be reduced to
an effective electronic one

$\displaystyle{\mathcal{H}}_{\rm eff}(t)=\langle\phi_{\rm ph}(t)|U_{\rm
NGS}(t)^{\dagger}{\mathcal{H}}(t)U_{\rm NGS}(t)|\phi_{\rm ph}(t)\rangle.$ (33)

Then the single-step electronic state evolution becomes

$\displaystyle|\psi_{e}(t+\delta t)\rangle\approx e^{-i\mathcal{H}(t)\delta
t}|\psi_{e}(t)\rangle.$ (34)

Here, the large Hilbert-space dimension ($>10^{8}$) requires stable and
efficient evaluation of the wavefunction, which in this study is based on the
Krylov-subspace method. The Krylov subspace for an instantaneous wavefunction
$|\psi_{e}(t)\rangle$ and Hamiltonian ${\mathcal{H}}_{\rm eff}(t)$ is defined
as

$\displaystyle\mathcal{K}_{n}(t)=\textrm{span}\\{|\psi_{e}(t)\rangle,{\mathcal{H}}_{\rm
eff}(t)|\psi_{e}(t)\rangle,\cdots,{\mathcal{H}}_{\rm
eff}(t)^{n-1}|\psi_{e}(t)\rangle\\}.$

A widely used Krylov-subspace generator is the Lanczos algorithm, where the
wavefunction can be approximated by Manmana  _et al._ (2005); Balzer  _et al._
(2012); Hochbruck and Lubich (1997); Innerberger  _et al._ (2020)

$\displaystyle|\psi(t+\delta t)\rangle$ $\displaystyle\approx$
$\displaystyle\exp\left[-iU_{n}(t)T_{n}(t)U_{n}(t)^{\dagger}\delta
t\right]|\psi(t)\rangle.$ (35)

Here, $U_{n}(t)$ is the basis matrix of the Krylov subspace
$\mathcal{K}_{n}(t)$, and the $T_{n}(t)$ is the tridiagonal matrix generated
by the Lanczos algorithm. Since the projected matrix is with dimension $n$,
much smaller than the Hilbert-space dimension, the evaluation of Eq. (35) is
much cheaper. The error of evaluating the vector propagation is well
controlled by $\epsilon_{n}\leq 12e^{-{(\rho\,\delta
t)^{2}}/{16n}}\left({e\rho\,\delta t}/{4n}\right)^{n}$ given the Krylov
dimension $n\geq\rho\,\delta t/2$ and spectral radius
$\rho=|E_{\max}-E_{\min}|$ Hochbruck and Lubich (1997). Depending on the
simulated pump strengths, we adopt $n=50-100$ in this paper.

Figure 6:  Dynamics of $N(\textbf{q},t)$ and $S(\textbf{q},t)$ for other
momenta [left to right: $(\pi/2,0)$, $(\pi,0)$, $(\pi/2,\pi/2)$, and
$(\pi,\pi/2)$] not shown in the main text. The upper two rows are obtained for
parameter Set 1, while the lower two rows are for parameter Set 2.

##  Appendix C Spin and Charge Dynamics for Other Momenta

In Fig. 2 in the main text, we show the time evolution for the charge and spin
structure factors for only a single momentum $\textbf{q}=(\pi,\pi)$, which
plays the dominant role in $d$-wave superconductivity. In this section, we
also present the calculations for other momenta.

The upper two rows in Fig. 6 show the dynamics for $u=6.6$ and $\lambda=4$
(Set 1). In contrast to the suppression of $N(\pi,\pi)$, the charge structure
factors are always enhanced for other momenta, due to the light-melting of the
dominant order. The dynamics for spin structure factors are more complicated:
They are enhanced during the pump, but for large pump strengths ($A_{0}\geq
0.5$) the structure factors decrease rapidly after the pump pulse. As
discussed in the main text, this means that a strong pump further suppresses
the spin fluctuations and is, thereby, unfavorable for $d$-wave
superconductivity.

In contrast to the dynamics with parameter Set 1, the evolution with parameter
Set 2 (lower rows in Fig. 6) suggests unchanged spin correlations after the
pump. Together with the dominant $(\pi,\pi)$ momentum in the main text, these
spin fluctuations give rise to the $d$-wave superconductivity after melting
the charge-ordered state.

Using the momentum distribution of the spin and charge structure factors, we
can also obtain the evolution of the correlation length. As shown in Fig. 7,
the charge correlation length decreases rapidly after the pump pulse enters.
Together with the reduction of $N(\pi,\pi)$ shown in Fig. 2, it reflects the
light melting of the charge instability, otherwise dominant at equilibrium.
This trend is distinct from the light-driven dynamics of the $d$-wave pairing
correlation, whose intensity is enhanced but the correlation length is
reduced.

Figure 7:  Evolution of the correlation length $\xi$ of the charge structure
factor renormalized by its equilibrium value $\xi_{0}$, obtained for (upper)
$u=6.6$ and (lower) $u=7$, respectively.

##  Appendix D Frequency Distribution of Dynamics

Figure 8:  Fourier spectra for the pump-induced evolution of $N(\pi,\pi)$ and
$S(\pi,\pi)$ for (left) $u=6.6$, (right) $u=7$, and different pump strengths
shown in the main text.

Besides the increase and decrease of the correlations, the postpump dynamics
reflect some underlying physical quantities. Figure 8 presents the Fourier
spectra of the charge and spin dynamics shown in Fig. 2 in the main text. For
the model parameter Set 1, the dynamics have low-energy periodicity
approximately $2\pi/\omega_{0}$, reflecting the role of phonons in forming and
melting a charge-ordered state. For the model parameter Set 2 close to the
phase boundary, the phonon energy becomes highly renormalized Wang  _et al._
(2016). Therefore, the low-energy spectrum becomes complicated and dependent
on the pump strengths. This strength dependence is also reflected by the
dynamics of other momenta in Fig. 6.

## References

  * Basov  _et al._ (2017) D. Basov, R. Averitt, and D. Hsieh, Nat. Mater. 16, 1077 (2017).
  * Wang  _et al._ (2018a) Y. Wang, M. Claassen, C. D. Pemmaraju, C. Jia, B. Moritz, and T. P. Devereaux, Nat. Rev. Mater. 3, 312 (2018a).
  * Mankowsky  _et al._ (2014) R. Mankowsky, A. Subedi, M. Först, S. Mariager, M. Chollet, H. Lemke, J. Robinson, J. Glownia, M. Minitti, A. Frano, _et al._ , Nature 516, 71 (2014).
  * Hu  _et al._ (2014) W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M. C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, Nat. Mater. (2014).
  * Kaiser  _et al._ (2014) S. Kaiser, C. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Liu, M. Le Tacon, T. Loew, D. Haug, B. Keimer, _et al._ , Phys. Rev. B 89, 184516 (2014).
  * Mitrano  _et al._ (2016) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, D. Jaksch, and A. Cavalleri, Nature 530, 461 (2016).
  * Buzzi  _et al._ (2020) M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, _et al._ , Phys. Rev. X 10, 031028 (2020).
  * Fausti  _et al._ (2011) D. Fausti, R. Tobey, N. Dean, S. Kaiser, A. Dienst, M. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * Nicoletti  _et al._ (2014) D. Nicoletti, E. Casandruc, Y. Laplace, V. Khanna, C. R. Hunt, S. Kaiser, S. Dhesi, G. Gu, J. Hill, and A. Cavalleri, Phys. Rev. B 90, 100503 (2014).
  * Först  _et al._ (2014) M. Först, R. Tobey, H. Bromberger, S. Wilkins, V. Khanna, A. Caviglia, Y.-D. Chuang, W. Lee, W. Schlotter, J. Turner, _et al._ , Phys. Rev. Lett. 112, 157002 (2014).
  * Casandruc  _et al._ (2015) E. Casandruc, D. Nicoletti, S. Rajasekaran, Y. Laplace, V. Khanna, G. Gu, J. Hill, and A. Cavalleri, Phys. Rev. B 91, 174502 (2015).
  * Khanna  _et al._ (2016) V. Khanna, R. Mankowsky, M. Petrich, H. Bromberger, S. Cavill, E. Möhr-Vorobeva, D. Nicoletti, Y. Laplace, G. Gu, J. Hill, _et al._ , Phys. Rev. B 93, 224522 (2016).
  * Nicoletti  _et al._ (2018) D. Nicoletti, D. Fu, O. Mehio, S. Moore, A. Disa, G. Gu, and A. Cavalleri, Phys. Rev. Lett. 121, 267003 (2018).
  * Cremin  _et al._ (2019) K. A. Cremin, J. Zhang, C. C. Homes, G. D. Gu, Z. Sun, M. M. Fogler, A. J. Millis, D. N. Basov, and R. D. Averitt, Proc. Nat. Acad. Sci. 116, 19875 (2019).
  * Sentef  _et al._ (2016) M. A. Sentef, A. Kemper, A. Georges, and C. Kollath, Phys. Rev. B 93, 144506 (2016).
  * Sentef  _et al._ (2017) M. A. Sentef, A. Tokuno, A. Georges, and C. Kollath, Phys. Rev. Lett. 118, 087002 (2017).
  * Coulthard  _et al._ (2017) J. Coulthard, S. R. Clark, S. Al-Assam, A. Cavalleri, and D. Jaksch, Phys. Rev. B 96, 085104 (2017).
  * Babadi  _et al._ (2017) M. Babadi, M. Knap, I. Martin, G. Refael, and E. Demler, Phys. Rev. B 96, 014512 (2017).
  * Murakami  _et al._ (2017) Y. Murakami, N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. B 96, 045125 (2017).
  * Dasari and Eckstein (2018) N. Dasari and M. Eckstein, Phys. Rev. B 98, 235149 (2018).
  * Bittner  _et al._ (2019) N. Bittner, T. Tohyama, S. Kaiser, and D. Manske, J. Phys. Soc. Jpn 88, 044704 (2019).
  * Li and Eckstein (2020) J. Li and M. Eckstein, Phys. Rev. Lett. 125, 217402 (2020).
  * Tindall  _et al._ (2020a) J. Tindall, F. Schlawin, M. Buzzi, D. Nicoletti, J. Coulthard, H. Gao, A. Cavalleri, M. Sentef, and D. Jaksch, Phys. Rev. Lett. 125, 137001 (2020a).
  * Tindall  _et al._ (2020b) J. Tindall, F. Schlawin, M. A. Sentef, and D. Jaksch, arXiv:2011.04417 (2020b).
  * Patel and Eberlein (2016) A. A. Patel and A. Eberlein, Phys. Rev. B 93, 195139 (2016).
  * Kennes  _et al._ (2019) D. M. Kennes, M. Claassen, M. A. Sentef, and C. Karrasch, Phys. Rev. B 100, 075115 (2019).
  * Michael  _et al._ (2020) M. H. Michael, A. von Hoegen, M. Fechner, M. Först, A. Cavalleri, and E. Demler, Phys. Rev. B 102, 174505 (2020).
  * Wang  _et al._ (2018b) Y. Wang, C.-C. Chen, B. Moritz, and T. Devereaux, Phys. Rev. Lett. 120, 246402 (2018b).
  * Kondo  _et al._ (2015) T. Kondo, W. Malaeb, Y. Ishida, T. Sasagawa, H. Sakamoto, T. Takeuchi, T. Tohyama, and S. Shin, Nat. Commun. 6, 7699 (2015).
  * He  _et al._ (2021) Y. He, S.-D. Chen, Z.-X. Li, D. Zhao, D. Song, Y. Yoshida, H. Eisaki, T. Wu, X.-H. Chen, D.-H. Lu, _et al._ , Phys. Rev. X (in press) (2021).
  * Faeth  _et al._ (2021) B. D. Faeth, S. Yang, J. K. Kawasaki, J. N. Nelson, P. Mishra, L. Chen, D. G. Schlom, and K. M. Shen, Phys. Rev. X 11, 021054 (2021).
  * Xu  _et al._ (2021) Y. Xu, H. Rong, Q. Wang, D. Wu, Y. Hu, Y. Cai, Q. Gao, H. Yan, C. Li, C. Yin, _et al._ , Nat. Commun. 12, 2840 (2021).
  * Boschini  _et al._ (2018) F. Boschini, E. da Silva Neto, E. Razzoli, M. Zonno, S. Peli, R. Day, M. Michiardi, M. Schneider, B. Zwartsenberg, P. Nigge, _et al._ , Nat. Mater. 17, 416 (2018).
  * Yang  _et al._ (2018) X. Yang, C. Vaswani, C. Sundahl, M. Mootz, P. Gagel, L. Luo, J. Kang, P. Orth, I. Perakis, C. Eom, _et al._ , Nat. Mater. 17, 586 (2018).
  * Mootz  _et al._ (2020) M. Mootz, J. Wang, and I. E. Perakis, Phys. Rev. B 102, 054517 (2020).
  * Chiriacò  _et al._ (2018) G. Chiriacò, A. J. Millis, and I. L. Aleiner, Phys. Rev. B 98, 220510 (2018).
  * Lemonik and Mitra (2019) Y. Lemonik and A. Mitra, Phys. Rev. B 100, 094503 (2019).
  * Sun and Millis (2020) Z. Sun and A. J. Millis, Phys. Rev. X 10, 021028 (2020).
  * Shi  _et al._ (2018) T. Shi, E. Demler, and J. I. Cirac, Ann. Phys. 390, 245 (2018).
  * Shi  _et al._ (2020) T. Shi, E. Demler, and J. I. Cirac, Phys. Rev. Lett. 125, 180602 (2020).
  * Shi  _et al._ (2019) T. Shi, J. I. Cirac, and E. Demler, Phys. Rev. Research 2, 033379 (2019).
  * Knörzer  _et al._ (2021) J. Knörzer, T. Shi, E. Demler, and J. Cirac, arXiv:2108.13730 (2021).
  * Wang  _et al._ (2019) Y. Wang, I. Esterlis, T. Shi, J. I. Cirac, and E. Demler, Phys. Rev. Research 2, 043258 (2019).
  * Keimer  _et al._ (2015) B. Keimer, S. Kivelson, M. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Davis and Lee (2013) J. S. Davis and D.-H. Lee, Proc. Natl. Acad. Sci. U.S.A. 110, 17623 (2013).
  * Scalapino  _et al._ (1986) D. Scalapino, E. Loh Jr, and J. Hirsch, Phys. Rev. B 34, 8190 (1986).
  * Gros  _et al._ (1987) C. Gros, R. Joynt, and T. Rice, Z. Phys. B Condens. Matter 68, 425 (1987).
  * Kotliar and Liu (1988) G. Kotliar and J. Liu, Phys. Rev. B 38, 5142 (1988).
  * Tsuei and Kirtley (2000) C. Tsuei and J. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Zhang and Rice (1988) F. Zhang and T. Rice, Phys. Rev. B 37, 3759 (1988).
  * Singh and Goswami (2002) A. Singh and P. Goswami, Phys. Rev. B 66, 092402 (2002).
  * Zheng  _et al._ (2017) B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Huang  _et al._ (2017) E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Science 358, 1161 (2017).
  * Huang  _et al._ (2018) E. W. Huang, C. B. Mendl, H.-C. Jiang, B. Moritz, and T. P. Devereaux, npj Quantum Mater. 3, 22 (2018).
  * Ponsioen  _et al._ (2019) B. Ponsioen, S. S. Chung, and P. Corboz, Phys. Rev. B 100, 195141 (2019).
  * Kokalj (2017) J. Kokalj, Phys. Rev. B 95, 041110 (2017).
  * Huang  _et al._ (2019) E. W. Huang, R. Sheppard, B. Moritz, and T. P. Devereaux, Science 366, 987 (2019).
  * Cha  _et al._ (2020) P. Cha, A. A. Patel, E. Gull, and E.-A. Kim, Phys. Rev. Research 2, 033434 (2020).
  * Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).
  * Ido  _et al._ (2018) K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 97, 045138 (2018).
  * Jiang and Devereaux (2019) H.-C. Jiang and T. P. Devereaux, Science 365, 1424 (2019).
  * Shen  _et al._ (2004) K. Shen, F. Ronning, D. Lu, W. Lee, N. Ingle, W. Meevasana, F. Baumberger, A. Damascelli, N. Armitage, L. Miller, _et al._ , Phys. Rev. Lett. 93, 267002 (2004).
  * Lanzara  _et al._ (2001) A. Lanzara, P. Bogdanov, X. Zhou, S. Kellar, D. Feng, E. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, _et al._ , Nature 412, 510 (2001).
  * Reznik  _et al._ (2006) D. Reznik, L. Pintschovius, M. Ito, S. Iikubo, M. Sato, H. Goka, M. Fujita, K. Yamada, G. Gu, and J. Tranquada, Nature 440, 1170 (2006).
  * Devereaux and Hackl (2007) T. P. Devereaux and R. Hackl, Rev. Mod. Phys. 79, 175 (2007).
  * He  _et al._ (2018) Y. He, M. Hashimoto, D. Song, S.-D. Chen, J. He, I. Vishik, B. Moritz, D.-H. Lee, N. Nagaosa, J. Zaanen, _et al._ , Science 362, 62 (2018).
  * Hubbard (1963) J. Hubbard, Theoretical Physics 276, 238 (1963).
  * Holstein (1959) T. Holstein, Annals of physics 8, 325 (1959).
  * Clay and Hardikar (2005) R. Clay and R. Hardikar, Phys. Rev. Lett. 95, 096401 (2005).
  * Nowadnick  _et al._ (2012) E. Nowadnick, S. Johnston, B. Moritz, R. Scalettar, and T. Devereaux, Phys. Rev. Lett. 109, 246404 (2012).
  * Nowadnick  _et al._ (2015) E. Nowadnick, S. Johnston, B. Moritz, and T. Devereaux, Phys. Rev. B 91, 165127 (2015).
  * Johnston  _et al._ (2013) S. Johnston, E. Nowadnick, Y. Kung, B. Moritz, R. Scalettar, and T. Devereaux, Phys. Rev. B 87, 235133 (2013).
  * Wang  _et al._ (2016) Y. Wang, B. Moritz, C.-C. Chen, C. Jia, M. van Veenendaal, and T. P. Devereaux, Phys. Rev. Lett. 116, 086401 (2016).
  * Mendl  _et al._ (2017) C. B. Mendl, E. A. Nowadnick, E. W. Huang, S. Johnston, B. Moritz, and T. P. Devereaux, Phys. Rev. B 96, 205141 (2017).
  * Zhang  _et al._ (2018) S. Zhang, Z. Wang, L. Shi, T. Lin, M. Zhang, G. Gu, T. Dong, and N. Wang, Phys. Rev. B 98, 020506 (2018).
  * Lu  _et al._ (2012) H. Lu, S. Sota, H. Matsueda, J. Bonča, and T. Tohyama, Phys. Rev. Lett. 109, 197401 (2012).
  * Al-Hassanieh  _et al._ (2013) K. A. Al-Hassanieh, J. Rincón, E. Dagotto, and G. Alvarez, Phys. Rev. B 88, 045107 (2013).
  * Rincón  _et al._ (2014) J. Rincón, K. Al-Hassanieh, A. E. Feiguin, and E. Dagotto, Phys. Rev. B 90, 155112 (2014).
  * Lu  _et al._ (2015) H. Lu, C. Shao, J. Bonča, D. Manske, and T. Tohyama, Phys. Rev. B 91, 245117 (2015).
  * Denny  _et al._ (2015) S. Denny, S. Clark, Y. Laplace, A. Cavalleri, and D. Jaksch, Phys. Rev. Lett. 114, 137001 (2015).
  * Raines  _et al._ (2015) Z. M. Raines, V. Stanev, and V. M. Galitski, Phys. Rev. B 91, 184506 (2015).
  * Okamoto  _et al._ (2017) J.-i. Okamoto, W. Hu, A. Cavalleri, and L. Mathey, Phys. Rev. B 96, 144505 (2017).
  * Freericks  _et al._ (2006) J. Freericks, V. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
  * Aoki  _et al._ (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * Werner and Millis (2007) P. Werner and A. J. Millis, Phys. Rev. Lett. 99, 146404 (2007).
  * Eckstein  _et al._ (2010) M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. B 81, 115131 (2010).
  * Dagotto (1994) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Zhang  _et al._ (1998) C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev. Lett. 80, 2661 (1998).
  * Brockt  _et al._ (2015) C. Brockt, F. Dorfner, L. Vidmar, F. Heidrich-Meisner, and E. Jeckelmann, Phys. Rev. B 92, 241106 (2015).
  * Lang and Firsov (1962) I. G. Lang and Y. A. Firsov, Sov. Phys. JETP 16, 1301 (1962).
  * Karakuzu  _et al._ (2017) S. Karakuzu, L. F. Tocchio, S. Sorella, and F. Becca, Phys. Rev. B 96, 205145 (2017).
  * Ohgoe and Imada (2017) T. Ohgoe and M. Imada, Phys. Rev. Lett. 119, 197001 (2017).
  * Costa  _et al._ (2020) N. C. Costa, K. Seki, S. Yunoki, and S. Sorella, Commun. Phys. 3, 80 (2020).
  * Chung  _et al._ (2020) C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R. White, Phys. Rev. B 102, 041106 (2020).
  * Maier  _et al._ (2005) T. A. Maier, M. Jarrell, T. Schulthess, P. Kent, and J. White, Phys. Rev. Lett. 95, 237001 (2005).
  * Fulde and Ferrell (1964) P. Fulde and R. A. Ferrell, Physical Review 135, A550 (1964).
  * Larkin and Ovchinnikov (1965) A. Larkin and Y. N. Ovchinnikov, Soviet Physics-JETP 20, 762 (1965).
  * (99) It is debatable whether a short-range superconductivity correlation can be conceptually classified as a light-induced superconducting state. Apart from the abstract definitions, a practical statement relies on the experimental consequences. Here, the short-range correlation can be detected by the onset of superconductivity gap.
  * Knap  _et al._ (2016) M. Knap, M. Babadi, G. Refael, I. Martin, and E. Demler, Phys. Rev. B 94, 214504 (2016).
  * Macridin  _et al._ (2006) A. Macridin, B. Moritz, M. Jarrell, and T. Maier, Physical review letters 97, 056402 (2006).
  * Wang  _et al._ (2017) Y. Wang, M. Claassen, B. Moritz, and T. Devereaux, Phys. Rev. B 96, 235142 (2017).
  * Shvaika  _et al._ (2018) A. Shvaika, O. Matveev, T. Devereaux, and J. Freericks, Condensed Matter Physics 21, 33707 (2018).
  * Wang  _et al._ (2018c) Y. Wang, T. P. Devereaux, and C.-C. Chen, Phys. Rev. B 98, 245106 (2018c).
  * Matveev  _et al._ (2019) O. Matveev, A. Shvaika, T. Devereaux, and J. Freericks, Phys. Rev. Lett. 122, 247402 (2019).
  * Budden  _et al._ (2021) M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, _et al._ , Nature Physics 17, 611 (2021).
  * Su and Zhang (2020) Y. Su and C. Zhang, Phys. Rev. B 101, 205110 (2020).
  * Mahmood  _et al._ (2021) F. Mahmood, T. Devereaux, P. Abbamonte, and D. K. Morr, arXiv:2108.04260 (2021).
  * Wandel  _et al._ (2020) S. Wandel, F. Boschini, E. Neto, L. Shen, M. Na, S. Zohar, Y. Wang, G. Welch, M. Seaberg, J. Koralek, _et al._ , arXiv:2003.04224 (2020).
  * Kogar  _et al._ (2020) A. Kogar, A. Zong, P. E. Dolgirev, X. Shen, J. Straquadine, Y.-Q. Bie, X. Wang, T. Rohwer, I.-C. Tung, Y. Yang, _et al._ , Nat. Phys. 16, 159 (2020).
  * Demler and Zhang (1995) E. Demler and S.-C. Zhang, Phys. Rev. Lett. 75, 4126 (1995).
  * Demler  _et al._ (1998) E. Demler, H. Kohno, and S.-C. Zhang, Phys. Rev. B 58, 5719 (1998).
  * Freericks  _et al._ (2009) J. Freericks, H. Krishnamurthy, and T. Pruschke, Phys. Rev. Lett. 102, 136401 (2009).
  * Lenarčič  _et al._ (2014) Z. Lenarčič, D. Golež, J. Bonča, and P. Prelovšek, Phys. Rev. B 89, 125123 (2014).
  * Shimizu and Yuge (2011) A. Shimizu and T. Yuge, Journal of the Physical Society of Japan 80, 093706 (2011).
  * Hackl  _et al._ (2020) L. Hackl, T. Guaita, T. Shi, J. Haegeman, E. Demler, and J. I. Cirac, SciPost 9, 048 (2020).
  * Manmana  _et al._ (2005) S. R. Manmana, A. Muramatsu, and R. M. Noack, AIP Conf. Proc. 789, 269 (2005).
  * Balzer  _et al._ (2012) M. Balzer, N. Gdaniec, and M. Potthoff, J. Phys. Condens. Matter 24, 035603 (2012).
  * Hochbruck and Lubich (1997) M. Hochbruck and C. Lubich, SIAM J. Numer. Anal. 34, 1911 (1997).
  * Innerberger  _et al._ (2020) M. Innerberger, P. Worm, P. Prauhart, and A. Kauch, The European Physical Journal Plus 135, 1 (2020).