Signatures of ultrafast reversal of excitonic order in Ta$_2$NiSe$_5$

In the presence of electron-phonon coupling, an excitonic insulator harbors two degenerate ground states described by an Ising-type order parameter. Starting from a microscopic Hamiltonian, we derive the equations of motion for the Ising order parameter in the phonon coupled excitonic insulator Ta$_2$NiSe$_5$ and show that it can be controllably reversed on ultrashort timescales using appropriate laser pulse sequences. Using a combination of theory and time-resolved optical reflectivity measurements, we report evidence of such order parameter reversal in Ta$_2$NiSe$_5$ based on the anomalous behavior of its coherently excited order-parameter-coupled phonons. Our work expands the field of ultrafast order parameter control beyond spin and charge ordered materials.

In the presence of electron-phonon coupling, an excitonic insulator harbors two degenerate ground states described by an Ising-type order parameter. Starting from a microscopic Hamiltonian, we derive the equations of motion for the Ising order parameter in the phonon coupled excitonic insulator Ta2NiSe5 and show that it can be controllably reversed on ultrashort timescales using appropriate laser pulse sequences. Using a combination of theory and time-resolved optical reflectivity measurements, we report evidence of such order parameter reversal in Ta2NiSe5 based on the anomalous behavior of its coherently excited order-parameter-coupled phonons. Our work expands the field of ultrafast order parameter control beyond spin and charge ordered materials.
Exploring new pathways to optically switch Ising-type electronic order parameters is a major theme of current ultrafast science. In recent years, a variety of out-ofequilibrium protocols have been developed for rapidly switching ferromagnetic [1][2][3], ferrimagnetic [4][5][6], antiferromagnetic [7][8][9], and ferroelectric [10][11][12] order parameters. However, far less is understood about the mechanisms for switching more exotic order parameters that are not of magnetic and charge dipolar type.
A particularly interesting case is the excitonic insulator (EI), a strongly correlated electronic phase realized through condensation of bound electron-hole pairs [13]. The free energy landscapes of the complex electronic order parameter and the real lattice order parameter of an EI are typically characterized by a Mexican hat with continuous U(1) symmetry and a parabola, respectively [ Fig. 1(a)]. However, strong electron-phonon coupling (EPC) induces a tilting of both the lattice and electronic potentials [14][15][16][17][18], reducing the U(1) symmetry to a discrete Z 2 Ising-type symmetry. Like in magnetic or charge dipole ordered ferroic materials, this leads to two degenerate ground states characterized by order parameters of equal magnitude but opposite phase [ Fig. 1 There is currently no experimental method to switch nor directly measure the phase of an EI order parameter on ultrashort timescales. An alternative strategy is to measure the phase of the coupled structural order parameter. However, existing time-resolved x-ray and electron diffraction techniques are not phase sensitive. Optical phase-resolved second harmonic generation measurements [12,19,20] have been used to measure the phase of structural order parameters in noncentrosymmetric ferroic materials [21][22][23][24], but all known EI candidates are centrosymmetric [25,26]. The possible presence of 180 • EI domains further complicates such measurements because domain averaging would cause overall signal cancellation.
In this Letter, we demonstrate via theory and experiment a pathway to optically switch an EI order parameter and to probe this reversal through a coherent EI order-parameter-coupled phonon (OPCP). The time evolution of the coupled EI and structural order parameters following impulsive laser excitation are derived from modeling a prototypical system Ta 2 NiSe 5 , which we assume to harbor an EI phase, with an elementary spinless twoband Hamiltonian. Our simulations reveal that the EI order parameter is stably reversed above a critical laser fluence, identifiable indirectly via a saturation of the coherent OPCP amplitude. Using time-resolved coherent phonon spectroscopy measurements, we experimentally verify this scenario and also demonstrate how switching can be controlled through the relative timing between successive laser excitation pulses.
The quasi-one-dimensional (1D) direct band-gap semiconductor Ta 2 NiSe 5 is reported to undergo an EI transition at a critical temperature T c = 328 K [26,27], accompanied by a weak orthorhombic-to-monoclinic structural distortion due to EPC [16]. Impulsive laser excitation below T c has been shown to coherently excite at least five distinct Raman-active phonons with frequencies near 1, 2, 3, 3.7, and 4 THz. The 1, 2, and 4 THz modes are sensitive to the EI transition at T c [18,[28][29][30] and thus constitute the OPCPs, while the 3 and 3.7 THz modes are reportedly not coupled to the EI order parameter and thus serve as a control.
The low energy electronic structure of Ta 2 NiSe 5 consists of a conduction band with Ta 5d orbital character and a valence band with Ni 3d-Se 4p hybridized orbital character. The EI instability is well captured by a 1D spinless two-band Hamiltonian with EPC [31][32][33][34] Ni 3d an OPCP mode of energy ω 0 at site i. The conduction and valence band dispersions are given by k = ∆ 2 + 2J c sin (ka/2) 2 and µ k = − ∆ 2 − 2J v sin (ka/2) 2 , respectively, where 2J c and 2J v are their bandwidths and ∆ is the band gap [33,34]. An on-site interband electronelectron interaction term V drives the excitonic pairing, and an EPC term g couples the electronic and phononic subsystems. These microscopic parameters have been experimentally determined [33,34].
From Eq.(1), we derive the equations of motion for the EI and structural order parameters, defined as Φ i = c i † v i and X i = b i † + b i , respectively, in two steps. We first derive the exact expression for the nonequilibrium free energy functional of Φ and X in the Keldysh path integral framework, where we include the light excitation via Peierls substitution. Then we obtain the equations of motion as the saddle point of the free energy, ignoring population in the conduction bands, spatial fluctuations of the order parameter field, and higher-order contributions O(Φ 6 ) (see Supplemental Material [35]), Here D is an effective diffusion coefficient, q is the electron charge, A is the light vector potential, g is the renormalized EPC coefficient for the electronic channel, and m and U are the second-and fourth-order expansion coefficients of the electron-electron interaction term, respectively. These parameters are all functions of J c , J v , ∆, g, and V [35]. We further introduce two phenomenological constants γ e and γ ph to account for damping of electronic and structural modes due to thermal fluctuation, whose values can be experimentally determined [35].
We first present a qualitative picture of how order parameter reversal occurs in our model. In equilibrium the electronic potential V (Φ) = 1 2 mΦ 2 + 1 4 U Φ 4 + 2g XΦ has either a tilted parabolic (m > 0) or tilted Mexicanhat (m < 0) form capturing the absence or presence of exciton condensation. For homogeneous optical excitation, one can ignore spatial derivatives of Φ. For pulsed excitation with infrared light, whose frequency well exceeds the electronic Higgs-Goldstone [33,46] and OPCP mode frequencies, one can also average out the fast oscillations of the perturbation and retain only its Gaussian envelope. Under these conditions one can make the simplification ). Here σ is the temporal width of the Gaussian pulse, α is a positive constant scaling factor that can be calculated analytically [35], and F is the pump fluence, the only tunable parameter in our model. In the EI phase (m < 0), optical excitation therefore acts to instantaneously increase m(t). The subsequent reduction of the EI order parameter, which occurs on a timescale much shorter than 2π/ω 0 , results in a sudden shift in the lattice potential V (X) = 1 2 ω 0 2 X 2 + 2gω 0 Re(Φ)X due to EPC [35], launching coherent oscillations through displacive excitation. In the low fluence regime, where the phonon oscillation amplitude is small enough such that X does not change sign, the direction of tilt of both V (Φ) and V (X) remains unchanged and so no switching occurs. However, above a critical fluence F c , where the phonon oscillation amplitude becomes large enough to change the sign of X, the tilting of both potentials is reversed and the system can relax into the switched state.
Numerical simulations of our model using experimentally determined material parameters for Ta 2 NiSe 5 and σ = 100 fs were carried out in the EI phase [35]. Figure  2 displays simulation results for F slightly greater than F c . At the moment of excitation t = 0 [ Fig. 2(b)], there is an instantaneous change in V (Φ) from Mexican-hat to parabolic form caused by the light-induced enhancement of m(t). The EI order parameter evolves rapidly to the new potential minimum with overdamped dynamics and is quenched within the pulse duration. This leads to a rapid shift in V (X), launching coherent oscillations of the underdamped OPCP, which shakes the electronic potential via EPC at the phonon frequency. Once the pulsed excitation is over, V (Φ) recovers a Mexican-hat form. However, as X crosses zero within the first half period of oscillation [ Fig. 2(c)], the tilting of V (Φ) is re-t < 0 ps t > 1 ps t = 0.2 ps t = 0 ps  versed with respect to the pre-pumped (t < 0) case [ Fig.  2(a)], sending Φ toward the new minimum on the negative side. As Φ crosses zero, the tilt of V (X) is also reversed due to EPC, thus pushing X to the new minimum on the positive side. The system then continues to oscillate about the reversed minima at the phonon frequency until the OPCP is damped out [ Fig. 2(d)]. Note that this model treats light as a coherent drive without considering heating-and cooling-induced changes in m(t). However, accounting for the latter merely shifts F c [35].
In lieu of probing the phase of X and Φ, we propose that the reversal can be identified via the pump fluence dependence of the OPCP amplitude. We solve Eqs. (2) and (3) to obtain the time evolution of Φ and X and then define the OPCP amplitude by its peak height in the fast Fourier transform (FFT) of X. A conventional Raman-active phonon is coherently launched through either displacive excitation or impulsive stimulated Raman scattering (ISRS) [47,48] with an amplitude that is linearly proportional to F . For F < F c , the amplitude of an OPCP also scales linearly with F . In this regime the structural order parameter simply oscillates about the initial potential minimum and thus behaves like a conventional phonon. But once the initial displacement of V (X) is large enough to enable escape to the opposite minimum (F > F c ), the amplitude ceases to grow. Order parameter reversal is thus marked by a saturation in the amplitude versus fluence curve.
We argue that experimental evidence for this phenomenon already exists in published studies of Ta 2 NiSe 5 . Werdehausen et al. [28] performed ultrafast optical reflectivity measurements at T = 120 K using a pump photon energy of 1.55 eV and observed clear coherent oscillations of the 1 THz OPCP as well as the uncoupled 3 THz phonon (ISRS). The pump fluence dependence of these two mode amplitudes is reproduced in Fig. 3. The 3 THz mode scales linearly with F , consistent with its assignment as a conventional phonon. In contrast, the 1 THz mode scales linearly with F only at low fluences and then saturates above ∼0.4 mJ/cm 2 , consistent with an order parameter reversal. By overlaying our simulation results [35] atop these experimental curves, we find close agreement (Fig. 3).
Our theory also predicts that the exciton condensate should be transiently quenched [m(t) → 0] above a critical fluence F = F * where the condition f (t) = |m| is satisfied. Eqs. (2) and (3) do not constrain F * to coincide with F c and our simulation shows that F * is clearly lower than F c in Ta 2 NiSe 5 (Fig. 3). Recently Tang et al. [30] performed time-and angle-resolved photoemission spectroscopy (tr-ARPES) measurements on Ta 2 NiSe 5 at T = 30 K and tracked the dynamics of the charge gap, a measure of Φ, immediately after pumping with 1.77 eV light polarized perpendicular to the chain direction (equivalent geometry to Ref. [28]). They found that the instantaneous gap size decreases linearly with increasing pump fluence and saturates above 0.29 mJ/cm 2 , which was interpreted as the point where Φ transiently collapses. The fact that this fluence is lower than 0.4 mJ/cm 2 , and is expected to be even lower if the experiment were conducted at 120 K, is consistent with our theory.
The dynamical nature of the order parameter reversal process suggests that it can be controlled not merely by the total pump energy deposited, but also by its distribution in time. To show this, we consider a situation where the sample is pumped by two identical pulses separated by time δt, with individual fluences F < F c but 2F > F c . For δt → 0 the system is effectively pumped by a single pulse exceeding F c and so reversal occurs, while for δt → ∞, the system relaxes back to the initial ground state before the second pulse arrives and so no reversal occurs. To qualitatively understand the behavior at intermediate δt values where the system is still dynamically evolving when the second pulse arrives, we recall that in the single pulse case, switching occurs once the OPCP amplitude is large enough to change the sign of

Sim.
FIG. 3: Experimental pump fluence dependence of the 1 THz (red circles) and 3 THz (blue circles) coherent phonon amplitudes in Ta2NiSe5 reproduced from Ref. [28]. Simulation results for an OPCP (red line) and a conventional ISRS phonon (blue line) are overlaid and horizontally scaled (α ≈ 700) to match the experimental data. Vertical dashed lines mark the calculated Fc and F * . The non-monotonic behavior of the OPCP amplitude just above Fc arises from strong feedback between Φ and X immediately after excitation [35].
X. Therefore, in the two-pulse case, switching possibly occurs if coherent oscillations of the OPCP induced by the first pulse can be sufficiently amplified by the second time-delayed pulse.
Previous studies have shown that impulsively and displacively excited conventional Raman-active phonons can be coherently amplified (suppressed) by a second pump pulse when δt is an integer (half-integer) multiple of the phonon period [49][50][51][52][53][54], dubbed in-phase (IP) and outof-phase (OP) pumping respectively. Therefore we simulated the effects of both IP and OP pumping on the order parameters of Ta 2 NiSe 5 using the same material parameters as before. We chose to simulate the OPCP at ω 0 /2π = 2 THz rather than at 1 THz because recent tr-ARPES data show that the most pronounced modulations of the valence band maximum occur at 2 THz, suggesting strong coupling to Φ [18]. As shown in Fig.  4(a), pumping by a single pulse with F slightly less than F c causes a rapid but incomplete reduction of Φ, followed by a slower recovery back to its original value on a timescale set by the damping of the 2 THz phonon. For the two-pulse case, OP pumping of the 2 THz phonon similarly leads to partial suppression of Φ without reversal, but reversal is achieved with IP pumping. This phenomenon is again manifested through an unconventional behavior of the OPCP. As shown in Figs. 4(b) and (c), OP pumping leads to suppression of the 2 THz phonon amplitude relative to the single pump case, resembling a conventional phonon because the oscillation

Sim.
Sim. is around the initial potential minimum. But, in contrast to conventional behavior, IP pumping does not lead to further amplification once X is excited to the opposite minimum.
To verify this prediction, we performed transient optical reflectivity measurements on Ta 2 NiSe 5 single crystals [35] using two identical pump pulses (σ = 80 fs) with variable δt. The light was polarized perpendicular to the chain direction and the fluence of each pulse was tuned slightly below F c to match our simulations. We chose a pump photon energy of 1 eV to enhance the 2 THz oscillations [35]. Figure 4(d) shows the fractional change in reflectivity (∆R/R) versus time for both IP and OP pumping of the 2 THz phonon, as well as for pumping with only a single pulse. All three curves exhibit fast (∼1 ps) exponential decay following pump excitation, corresponding to the charge relaxation process. Oscillations from the beating of several coherently excited phonons are also clear. A FFT of the data shows the most pronounced peaks at 2, 3, and 3.8 THz [35]. A focus on the 2 THz mode reveals that OP pumping strongly suppresses its amplitude relative to the single-pump case whereas IP pumping does not amplify it [ Fig. 4(e)], in quantitative agreement with our simulations [Fig. 4(c)]. In contrast, strong amplification occurs for the 3 and 3.8 THz modes [35], consistent with their uncoupled nature.
Our field theory description of the EI order parameter goes beyond the phenomenological time-dependent Landau theory [55][56][57][58] in that it allows the order parameter to explore the tilted Mexican-hat potential in the complex plane and can be naturally linked to microscopic parameters of the underlying lattice model. While more details including extension beyond the mean-field limit (∇ = 0), temperature dependence with T > 0, diffusion perpendicular to the surface [57], time-dependent damping [55], and anharmonic phonon coupling [59] can be added to refine the simulations, our minimal microscopic theory already captures the most salient physics and experimental features. These ideas and dynamical protocols apply not only to excitonic insulators, but also to any system featuring a continuous-symmetry-breaking electronic order parameter induced by coupling to a structural order parameter, such as a charge ordered system coupled to a Peierls distortion or an orbital ordered system coupled to a Jahn-Teller distortion. Therefore the OPCP behavior revealed here may be a general fingerprint of electronic order parameter switching.
We thank Rick Averitt, Edoardo Baldini, Swati Chaudhary, Nuh Gedik, Xinwei Li, Tianwei Tang  In order to connect the dynamics of the OPCP and the exciton condensate with the microscopic description of Ta 2 NiSe 5 , we start from a commonly used two-band semiconductor Hamiltonian with spinless fermions [1][2][3][4] and inter-band interactions, The bands are formed by the quasi-one-dimensional lattice and i, k are the corresponding lattice sites and momentum. Each operator and parameter is defined in the main text. We also set = 1 for simplicity. The equations of motion for the complex exciton order parameter Φ i = c † i v i and the real lattice displacement X i = b † i + b i can be obtain by expressing Eq. (1) in terms of a nonequilibrium path integral (i.e. in the Keldysh framework) and introducing Φ as a dynamic, bosonic Hubbard-Stratonovich field. This allows one to integrate out the fermionic modes and derive the equations of motion for Φ, X via saddle point equations.
In a path integral approach, the partition function Z corresponding to the Hamiltonian in Eq. (1) is formally obtained from a field integral of the form where D represents the common field integral measure andc i , c i ,v i , v i are independent Grassmann fields, representing the fermion modes in conduction and valence bands, and b * i , b i , Φ * i , Φ i are complex fields, corresponding to the phonon and exciton condensate modes. The Keldysh action S is obtained in the canonical way [5] and reads as S = S f + S b with the fermion part and the boson part − V Φ * c,l,t Φ q,l,t + Φ * q,l,t Φ c,l,t .
Here, Ψ l,t = (c 1,l,t , c 2,l,t , v 1,l,t , v 2,l,t ) is the fermion spinor in Keldysh space and each field carries an index triplet (i, l, t), which labels Keldysh component i (i = 1, 2 for Grassmann fields and i = c, q for complex fields), lattice site l and time t. The matrices B, C, W are the Keldysh space Green's function for phonons, conduction band, and valence band. They are diagonal in frequency and momentum space and W −1 k,ω is identical with k → µ k . Also, arXiv:2012.12912v1 [cond-mat.str-el] 23 Dec 2020 with η → 0 + . The matrix M l,t describes the local coupling of the fermions to the exciton and phonon fields M l,t = 1 V Φ c,l,t + gX c,l,t + σ x V Φ q,l,t + gX q,l,t .
The fermion part S f is quadratic in Grassmann fields and integration according to Grassmann calculus yields the formal expression where the trace includes the sum over Keldysh indices and lattice sites and an integral over time. In order to eliminate the phonon field from the nonlinear action S f , one performs a polaron-type shift Φ c/q,l,t → Φ c/q,l,t − g V X c/q,l,t . For small exciton field amplitudes Φ, Eq. (8) can be expanded up to fourth order in the fields (n ≤ 2) and in powers of derivatives. The equations of motion for the exciton condensate and the displacement, which are represented by the classical fields Φ c,l,t , X c,l,t are obtained via the saddle-point equations (and their complex conjugates) δS δb q,l,t = δS δΦ q,l,t = 0.
This yields the equations of motion by further introducing the perturbation of light via Peierls substitution, which is argued to work better for an electronically localized system [4], and assuming the lattice constant d = 1: with the parametersD,m,Ũ , Z depending on integrals over Green's functions and therefore on the band structure of the material and the temperature T of the system. Assuming a band gap ∆ and k B T < ∆, i.e. negligible population in the conduction bands, we can write out these parameters The equations in the main text are obtained via m = m/Z, U =Ũ /Z, D =D/Z, g = g/(ZV ), and f = Dq 2 A 2 /Z αF exp (− 4 ln 2t 2 σ 2 ). Also, we can obtain real valued equations of motion by defining Φ = φ + iη. To characterize the dephasing of the phononic and electronic channels, we also add phenomenological decaying terms to both branches. It is straightforward to add a −2γ ph ∂ t X term to the structural dynamical equation and γ ph can be determined experimentally. For the complex electronic order parameter, on the other hand, we rewrite the order parameter dynamical equation as followed: Here γ e is dimensionless. It expresses the ratio of dephasing dynamics due to a non-zero temperature to the coherent dynamics of the order parameter. Both dynamics are generated by the same effective free energy functional. It determines the dephasing time of the electronic Higgs/Goldstone oscillations as validated in Section II.
Also, we ignore the spatial diffusion and the phonon frequency shift due to an order of magnitude estimate g ≈ ω 0 V . With all of the above rectifications we have: We then construct the initial conditions for the above equations, which guarantee that X is static and Φ remains real and static before the light excitation: η| t=0 = 0.
After establishing the equations and the initial conditions, we can numerically solve the differential equations and trace the dynamics of the complex electronic order parameter Φ and the real structural order parameter X [ Fig. S1], as well as the dynamical free energy landscapes. By taking the fast Fourier transform (FFT) of X in the time interval from 0 ps to 20 ps with different pumping fluence values, we obtain the order-parametercoupled phonon (OPCP) amplitude versus fluence curves [ Fig. S2]. To simulate the two-pulse pumping situation, one simply adds another f term to Eq. (17) and (18) that is identical to the first, except that this f term is centered at the time when the second pulse arrives at the sample. The initial conditions are the same. Here we apply a FFT to X in the time interval between the arrival of the second pulse and 20 ps thereafter. We thus obtain the OPCP amplitude versus fluence at different time delays [Fig. S5].

II. DETERMINATION OF THE MICROSCOPIC PARAMETERS
From our experiment we get ω 0 /(2π) = 2 THz. The chosen microscopic parameter values in Ref. [3,4] reproduce the equilibrium band structure qualitatively well, therefore we adopt these and set g = 2 THz, V = 60 THz, ∆ = 40 THz, J c = J v = 40 THz, m = −17 THz, U = 132 THz, D = 13.3 THz, and g = 0.83g. The Higgs mode frequency −2m = 34 THz qualitatively matches the gap size ∆ [6]. The corresponding fast oscillation is beyond the time resolution of our experimental setup, and hence cannot be resolved. The phonon dephasing time we measured is approximately 3 ps, thus γ ph ≈ 0.3 THz. This set of parameter choices is self-consistent but may not be unique. We also demonstrate that the specific choice of the above microscopic parameters does not change the main conclusion of this paper, i.e. the reversal of the EI order [ Fig. S2(b)]. We simulate our incident light as a σ = 100 fs Gaussian irradiating the sample at t = 0. The pump fluence F is thus the only tunable parameter.
There is uncertainty in the determination of the electron decay rate γ e . Recent theories have demonstrated that the electronic system can oscillate in an amplitude (Higgs) and phase (Goldstone) mode around the transient free energy minimum [3], but there is no experimental evidence of such modes so far. A large γ e describes the overdamped case where the electronic subsystem adiabatically evolves into the transient free energy minimum, while a small γ e captures the underdamped case where Φ explores a larger region of the Mexican-hat potential via rapid oscillations of the Higgs and the Goldstone modes upon light excitation [ Fig. S1(h)]. We demonstrate that the nature of the electronic decay, whether overdamped or underdamped, does not change our main finding, i.e. the reversal of the EI order and the concurrent anomalous phonon amplitude dependence. We simulated the aforementioned two cases using either γ e = 1, which typically characterizes an overdamped scenario with no oscillation Despite the distinction in γ e , the dynamics of Φ and X are comparable qualitatively at times longer than 0.25 ps after the light excitation. The clear reversal when the pump fluence is higher than the critical fluence is realized independent of the value of γ e . Further, alteration of the electron dephasing rate has a minimal affect on the reversal critical fluence [ Fig. S2(a)]. Since an investigation of the behavior of the Higgs/Goldstone mode is beyond the scope of this work, we only simulate with the overdamped case hereafter.

III. DIFFERENCE BETWEEN THE CALCULATED AND EXPERIMENTAL CRITICAL FLUENCE
In Section I, we defined f αF exp (− 4 ln 2t 2 σ 2 ). To get the value of the simulated critical fluence F c in real units, we explicitly write out the scaling factor α = 2 ω ph 2 c 0 σ , where D is the parameter defined in Section I, q is electron charge, d is the lattice constant, R is the reflectance, is the reduced Planck's constant, ω ph is the pump light frequency, c is the speed of light, 0 is the vacuum permittivity, σ is the pulse duration. Our model predicts that F c should be smaller when the pump polarization is parallel to the chain direction because d is smaller and R is larger in this geometry. This is consistent with the data reported in Ref. [6].
We obtained a simulated F c ∼ 6 mJ/cm 2 under the same reported experimental conditions [6], which needs to be scaled down by 16 times to match the real experimental F c . There are several possible factors that give rise to this discrepancy. First, our model is an elementary model with a simplified two-band structure. There are also uncertainties in the determination of the parameters J c , J v , ∆, g, and V . It is especially hard to determine g and V from experiment because they are not directly reflected in band structure measurements. Even though the chosen parameter values in Ref. [3,4] are claimed to reproduce the experimental data well in Ref. [7], we find the band dispersions displayed in several other ARPES papers differ subtly, giving rise to an uncertainty in the determination of the parameters [8][9][10][11]. In Fig. S2(b) we show several OPCP amplitude versus pumping fluence curves obtained using different sets of parameter choices. It is clear that F c changes with the value of the microscopic parameters, but the qualitative trend stays the same.
Second, F c is temperature-dependent. Increasing temperature makes the bandwidth and bandgap smaller and the Mexican-hat minimum shallower, both of which promote the order-parameter switch, thereby decreasing F c . Our simulation corresponds to T = 0 K case, and therefore gives an upper bound of F c . As shown in Ref. [6], a 100 K increase in the temperature can make F c several times smaller. As such, it is possible that F c at finite temperature is much smaller than 0 K. To give an accurate estimate, a temperature dependent model needs to be considered.
Third, transient heating effects are not accounted for in our model. For pump photon energies above the insulating gap, which is the case in our experiment, there will be finite absorption leading to transient heating of the electronic subsystem. This will contribute towards quenching the electronic potential in a similar manner to coherent electric field driving. Therefore neglecting transient heating effects will cause F c to be over-estimated. On the other hand, this approximation may be more realistic for the case of sub-gap pumping where absorption is suppressed.
Taken altogether, these three factors lead us to conclude that our microscopic model and the corresponding simulation qualitatively reproduce the experimental observation. We also emphasize that regardless of the choice of parameters, the anomalous behavior of an OPCP is well reproduced in all simulations, demonstrating our conclusions are robust against the specific value of the parameters.

IV. ULTRAFAST HEATING AND SUBSEQUENT COOLING
We note that ultrafast heating and subsequent cooling in both the electronic and structural channels upon pumping are not considered in our coherent-drive microscopic model. We demonstrate here that ignoring the heating and the subsequent cooling does not change the major conclusions.
We first consider the electronic channel. Upon pumping, the electrons are excited into the conduction bands and the effective electronic temperature increases dramatically. The ultrafast heating of the electrons is accompanied by the transient restoration of higher symmetry, i.e. Mexican hat becomes parabolic, producing a qualitatively similar effect to that imparted by our coherent driving model. In both cases, one can deterministically engineer the final state through pumping. In our model we introduce light perturbation via a coherent Peierls phase, implying an immediate relaxation to lower symmetry as soon as the pulse excitation is over. Although this theoretical treatment is extensively utilized, in reality hot electrons will thermalize with the lattice through EPC with a characteristic depopulation time of 1 ps in Ta 2 NiSe 5 as measured by time-resolved optical and electron spectroscopy [6][7][8][9][10]12]. This implies the higher symmetry exists beyond the time duration of the pulse. To take this depopulation time into account, we simulate the pulse as a 1 ps exponential decay convolved with a Gaussian. The dynamics of the order parameters are very similar to the dynamics in the coherent quench case, except for the fact that the order parameters reach their stable states after longer time (1-2 ps) due to thermalization. We show the OPCP fluence dependence with thermalization in Fig. S2(c). Compared with the coherent quench case, the switch to the counterpart state occurs at a slightly lower fluence but the trend is the same. Therefore, we conclude here that the heating and cooling of the electrons do not influence our major conclusion that the switch is achievable by increasing fluence and observable via the OPCP fluence dependence.
To estimate the lattice heating, we use the formula ∆T = (1−R)F Cρδ to calculate the lattice effective temperature increase, where R is the reflectance, ρ is the density, C is the specific heat capacity, and δ is the optical penetration depth for the pump photon energy (1 eV) [13,14]. Our experiments were performed at 80 K with a pump fluence of 0.5 mJ/cm 2 . Using these values, we obtain a temperature increase of 15 K. Therefore, the lattice temperature is far below T c after pumping and the lattice temperature induced change of band structure is negligible. Thus, the lattice heating can also be ignored.

V. NON-MONOTONIC BEHAVIOR OF THE OPCP AMPLITUDE VERSUS FLUENCE ABOVE THE CRITICAL FLUENCE
As depicted in Fig. S2, the OPCP amplitude as a function of pump fluence shows a non-monotonic behavior when F > F c . The amplitude of the "oscillation" and the fluence where they emerge are dependent on the specific values of the microscopic parameters. This behavior arises from the strong feedback between the electronic and structural order parameter dynamics immedi-ately after excitation. Because Φ always responds more rapidly than X, the subtle mismatch of the time when the two order parameters cross zero will influence the phonon amplitude.
We take the case with overdamped dynamics as an example. In such a case, Φ relaxes into the potential minimum instantaneously, while X takes a longer time to settle into the minimum depending on the phonon damping rate. When the pump fluence just surpasses F c , once the pump excitation is over and X starts to approach zero from the negative side, Re(Φ) does not cross zero but exhibits a partial regression back to the initial state [ Fig. S1(a) and (c)]. This incomplete return to the initial state in turn tilts the phonon potential in the opposite direction as X is evolving, thus exerting resistance to the phonon and decreasing the phonon amplitude. However, as the pump fluence further increases, X crosses zero more quickly and the aforementioned temporal mismatch between X and Φ crossing zero will be smaller, leading to a slight increase of the phonon amplitude. This explains the non-monotonic behavior that occurs just above F c in Fig S2(a). After careful examination of the dynamics of both channels at each pump fluence, we find that this temporal mismatch occurs twice as the pump fluence is increased, yielding the two "dips" above F c in the overdamped cases [ Fig. S2(a) and (b)]. The "dips" finally disappear after the pump fluence is high enough so that Φ directly crosses zero without returning partially back to the initial state. Thereafter X and Φ will not experience any mismatch and the phonon amplitude shows a smooth monotonic dependence on fluence.
The non-monotonic behavior is stronger for the underdamped case because Φ undergoes Higgs/Goldstone oscillations [ Fig. S2(a)]. The feedback between the two channels thus also lasts longer, creating more complicated dynamics. In addition to the mismatch between the time when X and Φ cross zero as discussed in the overdamped case, the order parameters can now also oscillate back-and-forth between the minima on either side of zero [ Fig. S1(e)-(h)]. This makes the final state more sensitive to pump fluence compared to the overdamped case. In other words, the order parameter is more susceptible to reversal upon small changes in pump fluence, leading to sharper modulations of the phonon amplitude.
We also note that in the overdamped case where we account for electronic heating and cooling, increasing pump fluence will induce a greater number of reversals back to the initial state within one phonon period (0.5 ps), and the oscillations in Fig. S2(c) actually stem from these back-and-forth reversals. In contrast, the reversal only occurs once in the overdamped coherent pumping case. This discrepancy reveals that the long-time dynamics of the coupled system is strongly influenced by the dynamics immediately after the excitation. While these differences lead to quantitatively different long time behaviors, they do not alter the main conclusion of this work, i.e. the observation of the first reversal and the non-monotonic behavior of the OPCP amplitude versus fluence above F c .

VI. ISRS/DECP PHONON SIMULATION
As mentioned in the main text, conventional Raman active phonon modes are launched through the impulsive stimulated Raman scattering (ISRS) or displacive excitation of coherent phonons (DECP) mechanisms. Several references have discussed and summarized the disparities between and the unification of these two mechanisms [15][16][17][18]. We simulated conventional ISRS/DECP phonons using the simplified formula from [19]: where in the displacive case, F (t) = Dθ(t) and in the impulsive case F (t) = F δ(t). Convolving with the for the DECP case and F (t) = F exp(− 4 ln(2)t 2 σ 2 ) for the ISRS case, where D and F are normalized fluences. The dynamics for both cases are generally the same except for a π/2 phase shift. As shown in Fig. S2(a), the amplitude of the ISRS/DECP-launched phonon has a linear fluence dependence.

VII. EXPERIMENT AND FITTING DETAILS
Single crystals of Ta 2 NiSe 5 were grown by chemical vapor transport reaction. First, a powder of Ta 2 NiSe 5 was synthesized by the solid-state reaction from a stoichiometric mixture of its elements. They were sealed in an evacuated quartz tube and heated at 900 ℃ for 5 hours. Next, the powder (around 2 g) and chunks of iodine (50 mg) were loaded in a sealed quartz tube, which was put in a two-zone furnace. The temperature for the growth sides were kept at 875 ℃ and 800 ℃ respectively for one month. The samples were cleaved along the (010) direction immediately before the experiment to obtain a fresh and smooth surface.
For double-pump experiments, the sample temperature was fixed at 80 K. In the experimental setup, a Ti:sapphire amplified laser operating at 1 kHz produces 800 nm pulses with 40 fs time duration. A small portion of the power is used as the probe. The remainder seeds an optical parametric amplifier (OPA) and generates near infrared light tuned to 1200 nm with a duration of 80 fs, which is used for the the pump pulse(s). The fluence of each pump was set to ∼0.5 mJ/cm 2 , which is around F c at 80 K, considering the temperature dependence of F c [6]. Our simulation further substantiates that the fluence is slightly below F c as shown later. The polarizations of both pulses were set to be perpendicular to the (100) axis  3 THz IP  1 THz IP  3.8 THz IP  2 THz IP  3.8 THz OP FIG. S3: Transient reflectivity in response to both a single pump (black) and two pumps (other colors). In the doublepump data, the delays between the two pumps are set to be either in-phase (IP) or out-of-phase (OP) with a certain phonon, as indicated in the legend. Each curve is vertically offset for clarity.
of the sample, since pumping with a parallel polarization was reported to generate phonons less efficiently [6]. Transient differential reflectivity curves with doubleand single-pump are shown in Fig. S3. All curves exhibit a clear beat pattern. The quick rise upon the arrival of the pump and the ensuing exponential decay (∼ 1 ps) characterize photocarrier generation and recombination respectively, in agreement with previous results [6,12]. The beat pattern indicates the coexistence of multiple coherent phonons. Three phonons centered at around 2 THz, 3 THz and 3.8 THz are identified after taking the fast Fourier transforms (FFT) of the transient reflectivity data, as depicted in Fig. S4. The reported 1 and 4 THz phonons are missing [9,10,12].
To test whether we are able to amplify or suppress the different observed phonons, we excite the sample with two pump pulses with the time delay between them tuned to be either in-phase (IP) and out-of-phase (OP) with each phonon respectively. As such, there are six configurations in total. Note that when the time delay is IP with the 2 THz phonon, it is nearly OP with the 3 THz phonon. Similarly, when the time delay is IP with the 3.8 THz phonon, it is nearly OP with the 2 THz phonon. Therefore, two configurations are redundant and we executed the double-pump experiment using 4 different fixed delays between the two pumps. The corresponding double-pump FFT spectra are displayed in Fig. S4, together with the single pump FFT spectrum. In Fig. S4(b), it is demonstrated that when the time delay between the two pump pulses is resonant with the 3 THz phonon, the 3 THz phonon amplitude is amplified by almost two times. Simultaneously, the 3.8 THz phonon is slightly enhanced because this time delay is also partially IP with its period, while the 2 THz phonon is suppressed due to the nearly OP time delay. Similar analysis can be used to interpret the 3.8 THz IP and OP pumping cases in Fig. S4(c) and S4(d). However, an anomalous behavior is observed in the 2 THz IP pumping configuration as shown in Fig. S4(a). Although the 3 THz phonon is suppressed due to the OP time delay, and the 3.8 THz phonon is amplified due to nearly IP pumping, the 2 THz phonon is not enhanced. The reported 1 THz OPCP [6] is missing in our singlepump spectrum. Fig. S4(e) displays the FFT spectrum of the 1 THz IP pumping. 1 THz IP pumping should be nearly IP with all the phonons, but only the 3 THz and 3.8 THz phonon are amplified. There is still no 1 THz phonon after pumping resonantly with it, and the enhancement of the 2 THz phonon is negligible, echoing the 2 THz IP double-pump results. Similar results were reproduced at different spots on two samples. This further demonstrates the 2 THz phonon is an OPCP unlike the 3 THz and 3.8 THz phonons.
We also conducted coherent phonon spectroscopy measurement on 5% S-doped Ta 2 NiSe 5 . The FFT spectrum 14 × T2T Hz (green, corresponds to the 1 THz IP pumping), 1 × T2T Hz (red, corresponds to the 2 THz IP pumping), 0.7 × T2T Hz (blue, corresponds to the 3 THz IP pumping), and 0.54 × T2T Hz (yellow, corresponds to the 3.8 THz IP pumping). The single pump results are displayed as a reference. Each curve is vertically offset for clarity and normalized by the peak value of the single pump FFT.
shows a much weaker 2 THz phonon than the undoped case, but similarly intense 3 THz and 3.8 THz phonons, as shown in Fig. S4(f). This observation is further evidence that the 2 THz phonon is coupled to the EI order while the 3 and 3.8 THz phonons are not, since S-doped Ta 2 NiSe 5 exhibits weaker EI order with a lower T c [14].
In principle we can fit the time traces with damped oscillations superposed atop an exponentially decaying background: where A denotes the electronic background amplitude due to photocarrier generation with a decay time τ 0 , and C characterizes the long heat escape time. B i , τ ph,i , ν i , φ i are the amplitude, lifetime, frequency and phase of the ith phonon respectively. Here i runs from 1 to 3 corresponding to 2, 3 and 3.8 THz phonons.
Equivalently, we can also fit the peaks in the FFT spectrum. A damped oscillation in the time domain transforms into a Lorentzian in the frequency domain. Thus, we can fit the FFT data with three Lorentzians with the same definition of the corresponding parameters as defined above:

VIII. DETAILED COMPARISON OF DOUBLE-PUMP EXPERIMENT AND SIMULATION
Here we present a detailed comparison between experiment and simulation of the complete pump time-delay dependence of the 2 THz phonon in the double-pump scheme. We zoom in on the 2 THz phonon FFT peak at the aforementioned five different time delays as shown in Fig. S5(a). Note that the time delays are expressed as multiples of the 2 THz phonon period (T 2T Hz ). We simulated the time evolution of X upon two-pulse excitation using the same time delays as in our double-pump experiment and obtained a FFT spectrum for each delay. The simulation results using a pump fluence F = 0.96F c per pulse reproduce the experimental data well [ Fig. S5(b)]. Simulation with pump fluence F > F c fails to reproduce the experimental data even qualitatively.
As a comparison, the 3 THz and 3.8 THz phonons exhibit enhancement with IP pumping and suppression with OP pumping, resembling the ISRS simulation very well [ Fig. S6] and thus demonstrating their uncoupled nature. Also note that the measured frequency of all three phonons after double-pumping redshifts compared with the single-pump case due to higher net pumping fluences [ Fig. S5 and Fig. S6]. This softening may come from carrier-excitation-induced lattice softening and phonon anharmonicity [20][21][22], which are ignored in our microscopic model and unrelated to the main results.