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

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 [1,2]. 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 [3][4][5][6][7], 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 La 1.675 Eu 0.2 Sr 0.125 CuO 4 and La 2−x Ba x CuO 4 near 1/8 doping [8][9][10][11][12][13][14], 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) [15][16][17][18][19][20][21][22][23][24]. 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 [25], Floquet engineering of the Fermi surface [26], and (for the terahertz pump) parametric amplification [27]. 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 [28], 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 electronphonon 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.
Moreover, recent experiments have revealed the presence of fluctuating superconductivity above the transition temperature T c in overdoped cuprates and FeSe [29][30][31][32]. 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 [33][34][35]. Thus, the resonance and gap in the transient optical conductivity may cause a misassignment of the long-range superconductivity [36][37][38]. 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 electronphonon 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öhlichtype electron-phonon coupling [39][40][41][42]. As an extension of the equilibrium NGSED [43], this time-dependent method provides an accurate description of far-fromequilibrium states through the Krylov-subspace method and the Kählerization of the solvers. These advances of the numerical method allow the simulation of lightinduced 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 La 1.675 Eu 0.2 Sr 0.125 CuO 4 and La 2−x Ba x CuO 4 .
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 timedependent 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 [44,45]. It was widely believed that spin fluctuation, induced by electron correlation, is a viable candidate to provide the pairing glue for d-wave superconductivity [46][47][48][49]. The minimal model to represent this strong correlation is the single-band Hubbard model [50]. Based on rigorous numerical simulations, many important experimental discoveries have been reproduced using this model, such as antiferromagnetism [51], stripe phases [52][53][54][55], strange metallicity [56][57][58], and superconductivity [59][60][61]. However, increasing experimental evidence reveals the significant role of phonons in high-T c superconductors [62][63][64][65][66], in addition to the strong electronic correlation. Together with the crucial lattice effects observed in pump-probe experiments [3][4][5], 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 [67,68]. Its Hamiltonian is written as Here, c iσ (c † iσ ) annihilates (creates) an electron at site i with spin σ, and a i (a † i ) 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 ω 0 . The electron-electron interaction and electron-phonon coupling can be depicted by the dimensionless parameters u = U/t h and λ = g 2 /ω 0 t h , respectively. Here, we set ω 0 = t h in accord with common choices [69][70][71][72][73][74]. Throughout the paper, we focus on 12.5% hole doping simulated on an N = 4 × 4 cluster, corresponding to the La 1.675 Eu 0.2 Sr 0.125 CuO 4 and La 1.875 Ba 0.125 CuO 4 experiments [8][9][10][11][12][13][14]75].
In an undoped system with a well-defined nesting momentum, one can restrict the phonons to only the q = (π, π) mode [28,73].
However, in a doped Hubbard-Holstein model without commensurability, phonon modes at all momenta should be considered. As shown in Ref. 43, 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 [52][53][54][55], 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 where the vector potential 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: and fix the polarization as diagonalê pol = 1 √ 2 (x +ŷ) and the pump frequency as Ω = 4t h (close to the 800-nm laser for t h = 350meV). 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 [76][77][78][79]. 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 [80][81][82]. 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 [83,84], have achieved great success in solving correlated materials and pump-probe spectroscopies in the thermodynamic limit [85,86]. 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 [87,88]. This issue becomes even more severe when the electronphonon coupling is non-negligible due to the unbounded phonon Hilbert space [89,90].
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 lowdimensional systems. In contrast, the phonon-phonon interaction (anharmonicity) is usually weak and electronphonon coupling has a linear form, with the complexity coming instead from the unbounded local phonon Hilbert space. Therefore, if one can find an entangler transforma-tion 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 [39]. This leads to the wavefunction ansatz [43] with the electron density operator ρ q = iσ n iσ e −iq·r i , the phonon momentum 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 |ψ e is treated as a full many-body state, while the phonon wavefunction |ψ ph is treated as coherent Gaussian state (see Appendix. A for details). The entangler transformation involves momentum-dependent variational parameters λ q , describing the polaronic dressing. Physically, a larger dressing λ q ∼ g/ω 0 accurately describes the phonon and coupling energies, known as the Lang-Firsov transformation [91], while a smaller dressing allows a precise solution for the electron energy. Thus, the variational parameters λ q s are optimized numerically as a balance between these two effects. Note that all λ q s are independent real numbers and entangle the phonon momentum with electron density. This entangler transformation is demonstrated to be sufficient in equilibrium [39,42]. For systems with both strong electron-phonon coupling and electronic interactions, the electronic part of wavefunction |ψ e 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 [43,92,93], which is consistent with the recent DQMC-AFQMC simulations at the thermodynamic limit [94].
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 electrondensity-dependent phonon momentum, we introduce an extra cubic coupling between the position and the phonon density in the non-Gaussian transformation to construct the variational ansatz Here, the additional phase parameter ϕ q controls the ratio of position and momentum displacements. The time dependence of λ q and ϕ 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. (A9) 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. |Ψ(t) → |Ψ(t+δt) , we simultaneously evolve both the variational parameters and the full electronic wavefunction |ψ e (t) , as explained in Appendixes A and B, respectively.

IV. LIGHT-ENHANCED PAIRING CORRELATIONS
Figure 2(a) shows the square of local moment m 2 z = i (n i↑ −n i↓ ) 2 /N calculated for different u and λ of the 12.5% doped Hubbard-Holstein model, where N is the system size. m 2 z provides an indication of the overall spin moment and spin correlation. An abrupt change of m 2 z with varying u and/or λ 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 [61,95], although inthe 2D thermodynamic limit its existence is not fully established [96]. In contrast, at the large-λ 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 n d = i n i↑ n i↓ /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 m 2 z 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 [43,94]. We focus on this intermediate regime, which we believe corresponds to the situation in cuprate experiments [8][9][10][11][12][13][14]28].
We select two sets of model parameters inside this intermediate regime. Set 1 (u = 6.6 and λ = 4) lies in the center of the intermediate regime, while Set 2 (u = 7 and λ = 4) lies close to the boundary. We first focus on Set 1 and examine the light-induced dynamics for the charge structure factor N (q) = ρ −q ρ q /N and spin structure factor S(q) = ρ s −q ρ s q /N , where ρ s q = i (n i↑ − n i↓ )e −iq·r i . We focus on the momentum (c1-c2) The evolution of charge correlation N (π, π) after a linear pump of various strengths A0 ranging from 0.1 to 1 (denoted by different colors) with λ = 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(π, π) and (e1-e2) the d-wave pairing susceptibility P (d) . The pump pulse is sketched above. q = (π, π), 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 (π, π) and S(π, π) for Set 1. The same as the half-filled case [28], 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 where the d-wave pairing operator reads As Set 1 is relatively far from the spin-dominant regime, spin fluctuation is unstable, manifest as the drop of S(π, π) 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 ∼ 2π/ω 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 λ = 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(π, π) and P (d) after suppressing the charge fluctuations. Although the enhancement of S(π, π) 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 dwave 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 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 [8][9][10][11][12][13][14]. 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 [29][30][31][32]. This is a signature of strong correlations in quantum material distinct from the meanfield 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.
To better clarify the nature of this photoinduced manybody state, we further consider the spatial fluctuation through the FFLO pairing correlation with finite momentum P [97,98]. For k = 0, it recovers to the BCS pairing correlation, reflecting the total number of Cooper pairs. The correlation length ξ can be estimated through based on the definition of correlation length satisfying ∆ Figures 3(a) and (b) present the evolution of ξ, estimated using the smallest momentum accessible in our cluster k = (π/2, 0) and normalized by the equilibrium correlation length ξ 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 ≥ 0.5, ξ 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 ∆ 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 ξ saturates at a moderate value, about half of the equilibrium ξ 0 . To filter out the oscillatory dynamics, we extract the average enhancement of pairing correlation P d (t)/P d (−∞) and the suppression of ξ(t)/ξ(−∞) during the time window 10t −1 h < t < 30t −1 h , 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 ∆ ∼ P (d) e −|ri−rj |/ξ , our simulation reflects an enhancement of short-range superconductivity but a suppression of the long-range one [99].
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) 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 La 1.675 Eu 0.2 Sr 0.125 CuO 4 and La 1.875 Ba 0.125 CuO 4 [8][9][10][11][12][13][14]. Such a suppression of the coherence is associated with the nonlocal nature of d-wave Cooper pairs, which is in contrast to the lightenhanced strength and coherence of local pairing [16,17,22,23]. This anticorrelation between the strength and coherence is also distinct from previous light-induced competing orders [76][77][78][79], 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(ac) present the dynamics of various quantities discussed above for different frequencies.
As Ω 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 dwave superconductivity. As shown in Fig. 4(b), the lowfrequency pump with Ω = 2t h leads to a slightly weaker enhancement and a dramatic decrease of the correlation length ξ. 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 Ω ∼ 2ω 0 is reached [18,100], leading to strong fluctuations of phonons and overwhelming the d-wave pairing [74,101].
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 [35,[102][103][104][105], 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 −1 h 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 ∼ 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 [11]. 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 an-other type of light-induced superconductor (K 3 C 60 ) [106]. 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 La 1.875 Ba 0.125 CuO 4 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 [30][31][32], while the distinction from longrange 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 La 1.875 Ba 0.125 CuO 4 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 [107,108], 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 [109]. Light-induced CDW also has been recently observed in LaTe 3 with two competing CDW orders [110]. Here, we investigate specifically lightdriven cuprates, but appropriate modification [111,112] 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 [38]. 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 [102,104,[113][114][115]  To obtain the equations of motion for variational parameters in Eq. (5), where the phonon state is we notice that the generalization from U plrn to U NGS involves a rotation of phonon operators in the phase space. Therefore, the U NGS can be expressed by a composite transformation, leading to an equivalent form of Eq. (5): where the rotation is reflected in the unitary (Gaussian) transformation and η q is a variational matrix. In the derivation of variational equations of motion, it is usually convenient to employ the linearized formS q to parametrize the rotational transformation, i.e. U † rot R q U rot =S q R q . Thus, ignoring the redundant degrees of freedom, which can be absorbed into |ψ ph , we obtain the following relation between ϕ q andS q S q (t) = e iσyηq = cos ϕ q (t) − sin ϕ q (t) sin ϕ q (t) cos ϕ q (t) . (A4) In the above wavefunction prototype, ∆ R (vector), ξ q (matrix), λ q and ϕ 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∂ t Ψ(t) = H(t) Ψ(t) onto the tangential space in the variational manifold, which gives equations of motion for both the variational parameters and the electronic wavefunction |ψ e (t) . 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 with |ξ α 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 [116]. 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. (A5) minimizes errors with respect to the realistic time-dependent wavefunctions.
Taking advantage of the composite representation of U NGS , it is convenient to define and solve the Schrödinger equation in a rotating frame i∂ t Ψ (t) =H Ψ (t) . Thus, the rotated Hamiltonian becomes Here, we denote the general electron-electron interaction terms as H e−e and assume it to commute with the electronic density operator n i , which is satisfied in the Hubbard-Holstein model Eq. (1). This assumption substantially simplifies the derivation and leads to the following Eq. (A9).
On the left-hand side of the rotated Schrödinger equation, we obtain Here, S q = e iσyξq is the linearized transformation matrix of the U GS and the effective Hamiltonian becomes for the product state |ψ ph ⊗ ψ e . Here, α takes the x and y directions, while δ α 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 q on top of the original H e−e , whose variational expression is (Note, that we employed the assumption that H e−e commutes with n i and, therefore, U plrn .) For the kinetic energy, the phonon-dressed effective hopping term can be reformulated in a closed form when taking the expectation with respect to the phonon Gaussian state |ψ ph . Now with both the time derivative and the effective Hamiltonian transformed into the product-state basis, we can project Eqs. (A8) and (A9) onto various tangential vectors sequentially. For the zeroth-order tangential vector, which is proportional to the phonon vacuum state |0 ph , we obtain the equation of motion for the electronic state: The (phonon-traced) H eff together with the scalar terms in the square brackets govern the time evolution of the electronic wavefunction |ψ e . In the NGSED framework, we treat |ψ e 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) |ψ ph can be obtained by projecting Eqs. (A8) and (A9) onto the first-order tangential vector (R T 0 S T 0 |0 ph ⊗ |ψ e ) and the second-order tangential vector (R † q R q |0 ph ⊗ |ψ e ). After some algebra, these two equations are and respectively. Here, ρ f = N e /N is the average electron density, and the renormalized phonon energy matrix is 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 Γ q = S q S † q and track the evolution of Γ q instead of S q . Equation (A14) becomes Finally, the equations of motion of λ q and ϕ q , which appear in U NGS and entangle phonons and electrons, can be obtained by projecting to the third-order tangential vector a † q ρ q |0 ph ⊗ |ψ e . These two equations are and respectively, where the modulated electronic correlation and the density correlation is C q = ρ −q ρ q . 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 |ψ e [43]. The most important observable in this manuscript is d-wave pairing correlation, which can be explicitly written as This expression contains the pairing correlation in the electronic part of wavefunction −4e iq·r cos q x 2 cos q y 2 Here, ther in the last summation denotes the half-unitcell-shifted coordinatesr = r +x/2 −ŷ/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 ω 0 = 5t h . For a medium pump strength (A 0 = 0.4) with frequency Ω = 5t h , the relative errors for two dominant correlations, i.e., the spin structure factor S(π, π) 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 Then the single-step electronic state evolution becomes |ψ e (t + δt) ≈ e −iH(t)δt |ψ e (t) . (B2) 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 |ψ e (t) and Hamiltonian H eff (t) is defined as K n (t) = span{|ψ e (t) , H eff (t)|ψ e (t) , · · · , H eff (t) n−1 |ψ e (t) }.
A widely used Krylov-subspace generator is the Lanczos algorithm, where the wavefunction can be approximated by [117][118][119][120] |ψ(t + δt) ≈ exp −iU n (t)T n (t)U n (t) † δt |ψ(t) .(B3) Here, U n (t) is the basis matrix of the Krylov subspace 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. (B3) is much cheaper. The error of evaluating the vector propagation is well controlled by n ≤ 12e −(ρ δt) 2 /16n (eρ δt/4n) n given the Krylov dimension n ≥ ρ δt/2 and spectral radius ρ = |E max − E min | [119]. Depending on the simulated pump strengths, we adopt n = 50 − 100 in this paper.

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 q = (π, π), 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 λ = 4 (Set 1). In contrast to the suppression of N (π, π), 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 ≥ 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 (π, π) 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 (π, π) 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.

Appendix D: Frequency Distribution of Dynamics
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π/ω 0 , reflecting the role of phonons in forming and melting a chargeordered state. For the model parameter Set 2 close to the phase boundary, the phonon energy becomes highly renormalized [73]. 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.
FIG. 7. Evolution of the correlation length ξ of the charge structure factor renormalized by its equilibrium value ξ0, obtained for (upper) u = 6.6 and (lower) u = 7, respectively.