Twisting Anderson pseudospins with light:Quench dynamics in THz-pumped BCS superconductors

Twisting Anderson pseudospins with light:
Quench dynamics in THz-pumped BCS superconductors

Yang-Zhi Chou Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA Department of Physics and Center for Theoretical Quantum Matter, University of Colorado Boulder, Boulder, Colorado 80309, USA    Yunxiang Liao Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA    Matthew S. Foster Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA Rice Center for Quantum Materials, Rice University, Houston, Texas 77005, USA
July 15, 2019

We study the preparation (pump) and the detection (probe) of far-from-equilibrium BCS superconductor dynamics in THz pump-probe experiments. In a recent experiment [R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Phys. Rev. Lett. 111, 057002 (2013)], an intense monocycle THz pulse with center frequency was injected into a superconductor with BCS gap ; the subsequent post-pump evolution was detected via the optical conductivity. It was argued that nonlinear coupling of the pump to the Anderson pseudospins of the superconductor induces coherent dynamics of the Higgs (amplitude) mode . We validate this picture in a two-dimensional BCS model with a combination of exact numerics and the Lax reduction method, and we compute the nonequilibrium phase diagram as a function of the pump intensity. The main effect of the pump is to scramble the orientations of Anderson pseudospins along the Fermi surface by twisting them in the -plane. We show that more intense pump pulses can induce a far-from-equilibrium phase of gapless superconductivity (“phase I”), originally predicted in the context of interaction quenches in ultracold atoms. We show that the THz pump method can reach phase I at much lower energy densities than an interaction quench, and we demonstrate that Lax reduction (tied to the integrability of the BCS Hamiltonian) provides a general quantitative tool for computing coherent BCS dynamics. We also calculate the Mattis-Bardeen optical conductivity for the nonequilibrium states discussed here.


I Introduction

A series of recent optical pump-probe experiments has reignited interest in far-from-equilibrium superconductivity Fausti et al. (2011); Matsunaga and Shimano (2012); Matsunaga et al. (2013); Mansart et al. (2013); Matsunaga et al. (2014); Hu et al. (2014); Mankowsky et al. (2015); Matsunaga and Shimano (2015); Mitrano et al. (2016); Nicoletti and Cavalleri (2016). Most of the experiments fall into two broad classes. In the first class, mid-infrared radiation ( 10 THz) is injected with the goal of enhancing pairing, either by destabilizing competing orders Fausti et al. (2011) or by exciting optical phonons that participate in pairing Hu et al. (2014); Mankowsky et al. (2015); Mitrano et al. (2016). The essential idea of the latter is that by modulating the “pairing glue,” one might induce transient high(er) temperature superconductivity than is possible in equilibrium Knap et al. (2016); Kennes et al. (2017).

A complication of any high frequency excitation with ( is the BCS gap) is that the radiation can break Cooper pairs into hot quasiparticles, and these can serve as an efficient mechanism for rapid dissipation and thermalization. The second class of experiments Matsunaga and Shimano (2012); Matsunaga et al. (2013, 2014) used lower frequency radiation ( 1 THz), which might mitigate heating. In Ref. Matsunaga et al. (2013), a near monocycle pulse with center frequency was injected into a low-temperature BCS superconductor. Because most of the spectral weight lies below the optical gap edge , a weak pulse would not be expected to couple strongly to the sample. However, an intense pulse as used in Ref. Matsunaga et al. (2013) couples nonlinearly Tsuji and Aoki (2015) to the Cooper pairs (“Anderson pseudospins” Anderson (1958)) of the superconductor, and it was argued that this leads to a coherent excitation of the Higgs amplitude mode Volkov and Kogan (1974); Littlewood and Varma (1981); Barankov et al. (2004); Yuzbashyan et al. (2005a, b); Barankov and Levitov (2006a); Yuzbashyan et al. (2006); Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006); Gurarie (2009); Papenkort et al. (2007, 2008); Krull et al. (2014); Tsuji and Aoki (2015); Yuzbashyan et al. (2015); Cea et al. (2016); Varma (2002); Pekker and Varma (2015). The advent of ultrafast pump-probe detection can allow the observation of coherent many-body quantum dynamics in the solid state within a finite temporal window (“prethermalization plateau” Kinoshita et al. (2006); Rigol et al. (2007)). In the experiments Matsunaga and Shimano (2012); Matsunaga et al. (2013, 2014) detection is performed over a window of about 10 picoseconds (ps), well before thermalization occurs (likely due to acoustic phonons on a timescale of 100 ps Demsar et al. (2003)).

Figure 1: The Anderson pseudospin textures and the Lax roots of the spectral polynomial immediately after the application of an intense THz pulse. (a)-(c) The spin textures after light exposure, corresponding to the roots shown in (d)-(f), respectively. The pair of roots away from the real axis is called isolated and these roots encode the key properties of the BCS state. The system is a quarter-filled square lattice tight-binding model with 24-by-24 sites (note that our subsequent quantitative results are obtained from much larger systems). is the ground state order parameter of the BCS Hamiltonian; is the tight-binding parameter. The pump energies [ defined in Eq. (6)] are (a),(d) (very weak deformation from the BCS ground state), (b),(e), (intermediate strength deformation), and (c),(f) (strong deformation). The isolated roots for the deformed spin textures encode the asymptotic value of in the “pre-thermalization plateau.” (The chemical potential is absorbed into the dispersion in these plots.) More intense pulses push the isolated roots towards the real axis. For strong pulses, the isolated roots merge with the real axis and the system enters a phase of dynamical gapless superconductivity Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006) with ; see the “quench phase diagram” in Fig. 2.

In this paper, we model pump-probe experiments as typified by Ref. Matsunaga et al. (2013), using a combination of numerics and the Lax reduction method Yuzbashyan et al. (2006). We view the THz pump excitation as the preparation of an initial condition to the subsequent free (unperturbed) BCS time evolution, i.e. as a type of quantum quench Greiner et al. (2002); Barankov et al. (2004); Barankov and Levitov (2006b); Yuzbashyan et al. (2005a, b, c, 2006); Barankov and Levitov (2006a); Calabrese and Cardy (2006); Gurarie (2009). We show that the main effect of the pump is to scramble the orientations of Anderson pseudospins along the Fermi surface by twisting them in the -plane. We also show that more intense pump pulses can induce a far-from-equilibrium phase of gapless superconductivity (“phase I”), originally predicted in the context of interaction quenches in ultracold atomic systems Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006). An intense THz pump can reach phase I at much lower energy densities than an interaction quench.

Most existing theoretical studies of coherent far-from-equilibrium BCS superfluid dynamics focus on interaction quenches of isolated systems, as could be realized in ultracold atom experiments Bloch et al. (2008); Polkovnikov et al. (2011); Eisert et al. (2015). Due mainly to the integrability of the BCS model Yuzbashyan et al. (2005a, b, 2006), many asymptotically exact results are already known. The most important is the identification of three different nonequilibrium phases (dubbed “I, II, and III” in Foster et al. (2013); Yuzbashyan et al. (2015), reviewed below), predicted to occur in quenches of an -wave superconductor Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006), superfluid Foster et al. (2013, 2014); Liao and Foster (2015), BCS-BEC condensate Yuzbashyan et al. (2005c); Gurarie (2009); Yuzbashyan et al. (2015), and spin-orbit-coupled fermion condensate Dong et al. (2015); Dzero et al. (2015).

Ultrafast pump-probe spectra in BCS superconductors have been studied numerically using the density matrix formalism Papenkort et al. (2007, 2008); Krull et al. (2014). These studies revealed asymptotic post-pump dynamics identical to a small “phase II” interaction quench, in which approaches a nonequilibrium value as (neglecting dissipative processes that would ultimately induce thermalization). The approach to involves a characteristic damped oscillation at frequency Volkov and Kogan (1974); Yuzbashyan et al. (2006), also seen in the experiment Matsunaga et al. (2013). The calculations in Refs. Papenkort et al. (2007, 2008); Krull et al. (2014) assumed coupling to the electromagnetic field due to the finite (but large) photon wavelength m. In the present paper, we classify the dynamics according to the interaction quench phase diagram Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006); Foster et al. (2013); Yuzbashyan et al. (2015), and we assume that the strongest coupling to the pump excitation is due to non-quadratic band curvature Tsuji and Aoki (2015). A related experiment Matsunaga et al. (2014) shows third harmonic generation in the steady-state response; this was studied theoretically in Tsuji and Aoki (2015); Cea et al. (2016).

In this work, we demonstrate that the Lax reduction method Yuzbashyan et al. (2006) (tied to the integrability of the BCS Hamiltonian) provides a general quantitative tool for computing coherent nonlinear BCS dynamics. We also calculate the Mattis-Bardeen Mattis and Bardeen (1958); Tinkham (1996) optical conductivity for the nonequilibrium superconducting states predicted here.

The rest of this paper is organized as follows: In Sec. II, we review the interaction quench dynamics of -wave BCS superconductors. We then highlight our main findings for a THz pump probe, including the dynamical phase diagram given by Fig. 2 and predictions for the optical conductivity. In Sec. III, we provide details of the light-superconductor interaction that twists Anderson pseudospins during the application of the pump. In Sec. IV, we calculate the optical conductivity for nonequilibrium superconductors. In Sec. V, we introduce the “R-ratio” to quantify the relative energy injected by the pump, and we compare the threshold energies for reaching phase I in THz pump versus interaction quenches. We summarize in Sec. VI.

Ii Quantum quench via pump-probe: Main results

In this section, we review the dynamical phases of nonequilibrium BCS superconductors induced by interaction quenches. We then summarize our main results for a THz pump probe. Details of our calculations appear in Secs. III and IV.

ii.1 BCS model, review of interaction-quench dynamics, and relation to the “THz pump quench”

The BCS Hamiltonian can be expressed in terms of Anderson pseudospins Anderson (1958),


where , is the single particle dispersion, is the chemical potential, and is the coupling strength. The spin operators are


and . In the thermodynamic limit, Eq. (1) is equivalent to the mean field Hamiltonian,


The instantaneous order parameter is self-consistently determined


The equation of motion for an Anderson pseudospin is


The vector potential will be used to encode the pump electric field Tsuji and Aoki (2015).

Previous studies of quench dynamics in superconductors focused on interaction quenches, in which the ground state of Eq. (1) with an initial coupling strength is time-evolved according to the post quench Hamiltonian, given by Eq. (1), with “final” coupling strength (and at all times). The long-time behavior of the order parameter following the quench can be predicted with the Lax reduction method Yuzbashyan et al. (2006), which exploits the integrability of .

Interaction quenches from an initial BCS state can produce three possible nonequilibrium phases Yuzbashyan et al. (2006); Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006); an interaction-quench dynamical phase diagram can be constructed in the - plane, where () denotes the ground state BCS gap of the pre-quench (post-quench) Hamiltonian, and corresponds to equilibrium (no quench) Foster et al. (2013); Yuzbashyan et al. (2015). Small quenches with reside in phase II, wherein asymptotes to a non-zero constant value as . Note that for a quench. Large quenches from strong to weak pairing () can induce phase I, in which Yuzbashyan and Dzero (2006); Barankov and Levitov (2006b); this is a phase of gapless (fluctuating) superconductivity. Large quenches from weak to strong pairing () can induce phase III, in which exhibits persistent oscillations Barankov et al. (2004); Barankov and Levitov (2006b). Phase III is a self-generated Floquet phase without external periodic driving Foster et al. (2014); Liao and Foster (2015), and is connected to the instability of the normal state.

In the THz pump-probe experiments of Refs. Matsunaga et al. (2013); Matsunaga and Shimano (2015), the injected pump pulse can be viewed as preparing an initial condition for subsequent free evolution in the absence of the field. In this picture, the interaction strength is not modified by the radiation. Instead, for a system at zero temperature, the pump pulse deforms the BCS ground state into a new pure state over the duration of the pulse [wherein ]. The deformed state subsequently evolves with the original Hamiltonian and . In analogy with the interaction quench, we identify the original Hamiltonian as the post-quench Hamiltonian, i.e., the ground state gap is denoted in the following.

Figure 2: The nonequilibrium phase diagram of an -wave BCS superconductor subject to an ultrashort ( monocycle) THz pulse with center frequency (“pump quench”). We consider Eq. (1) on a 2D square lattice with a linear size (10001000 sites) at quarter filling. Values of are extracted from the isolated Lax roots corresponding to the spin configuration immediately after cessation of the pump pulse. We verified that these results agree with extracted from Eq. (4), using the numerical integration of Eq. (5) to large times. We plot as a function of the peak relative field intensity [Eq. (6)] with different values of . Here denotes the ground state BCS gap, and is the hopping strength. The polarization angle for all cases, i.e., the radiation is polarized along the -axis. The results weakly depend on . More intense pulses produce larger deformations of the pseudospins along the Fermi surface (see Fig. 1), leading to a suppression of the asymptotic BCS gap . For small values of , finite size effects become significant, and it becomes difficult to resolve the isolated roots via the Lax method. For these cases the results were obtained via numerical integration of Eq. (5). Phase I () can be achieved with (yellow shaded region) in all cases. The phase boundary for is close to .
Figure 3: The polarization angle dependence of for the pump quench. The system is the same as in Fig. 2. We plot as a function of pump energy () with two polarization angles, . Here we set . Radiation polarized along the -axis () generates nonequilibrium dynamics most efficiently, due to the density of states accumulation along the parallel Fermi line segments with . For , we did not observe phase I dynamics. To determine whether the latter is possible, it would be necessary to retain higher order terms in Eqs. (10) and (11).
Figure 4: The real time evolution of the superconducting order parameter. The system is a square lattice with linear size at quarter filling. The pump-pulse polarization angle is set to 0, and . The dashed line corresponds to pump energy (phase II); the solid line corresponds to pump energy (phase I). The red regions indicate the evolution during the application of the pump pulse. The blue regions mark the subsequent free time evolution with the time-independent BCS Hamiltonian. Inset: The phase I dynamics of . The order parameter decays exponentially in time Yuzbashyan and Dzero (2006).
Figure 5: The real part of the Mattis-Bardeen conductivity in phase II. The result depends on , but also on the full nonequilibrium distribution function that determines the occupation of quasiparticle states in the asymptotic steady state. For this plot and Fig. 6, we use the distribution function for an interaction quench whose analytical form is exactly known. We expect that results for the pump quench are qualitatively identical. The interaction quench is specified by the quench parameter , where () is the interaction strength of the pre-quench (post-quench) Hamiltonian. Since the Mattis-Bardeen formulae Mattis and Bardeen (1958) assume the dirty limit ( is the elastic scattering time), we consider a 2D particle-hole symmetric superconductor with a uniform density of states. The dynamical phase II-I boundary is located at . The real part of conductivity shows a cusp at for all phase II quenches. Although our system is specified at all times by a nonequilibrium BCS pure state, the optical conductivity for a quench “looks” thermal. In particular, as approaches the phase II-I boundary , the real part of the optical conductivity behaves the same way as an equilibrium superconductor at temperature approaching from below Tinkham (1996). The steady-state optical conductivity for phase I is indistinguishable from a normal metal.
Figure 6: The prefactor of the imaginary part of the Mattis-Bardeen conductivity, corresponding to the phase II quenches shown in Fig. 5. The horizontal axis specifies the quench parameter . The imaginary part of the conductivity can be fitted by for . The value in equilibrium with is . is the critical value of the quench parameter that separates phase II and phase I.

ii.2 Pump-pulse quench: Results and phase diagram

Light couples via the vector potential to the -component of the Anderson pseudospin -field [Eq. (3b)]. We neglect the small photon momentum in the following, and encode the uniform electric field of the THz pump via . The strength of the pump excitation can be measured by the peak field intensity (energy)


where denotes the lattice spacing and is the hopping strength that sets the electronic bandwidth. For pump energies , the dispersion can be expanded in powers of . In a clean single band superconductor, linear coupling to does not induce time evolution from a pure BCS state Mahan (1981); Dai and Lee (2017); the same is true for optical subgap excitation at frequency in a dirty superconductor. An intense pump pulse can couple through the band curvature at order as long as the dispersion is not purely parabolic [see Eqs. (8) and (9)].

We consider a two-dimensional (2D) square lattice model at quarter-filling qua (). Application of an intense, monocycle pump pulse with center frequency ( denotes the ground state gap) generates dynamics that deform the Anderson pseudospin texture of the initial state. Within the self-consistent mean field framework given by Eq. (3), the light-matter interaction does not break Cooper pairs (remove Anderson pseudospins), and the system resides in a BCS product state at all times (albeit one with time-evolving coherence factors). In a real experiment such as Ref. Matsunaga et al. (2013), the superconductor is disordered and quasiparticles are always excited by the linear coupling to , since some spectral weight of the ultrashort pump pulse extends above . We assume that the effects of these quasiparticles on the coherent evolution of can be neglected on the timescale of the experiment, which is typically 10 ps Matsunaga and Shimano (2012); Matsunaga et al. (2013, 2014); Matsunaga and Shimano (2015).

In the presence of the pump-pulse light, we use the fourth order Runge-Kutta method to numerically integrate the spin equations of motion, given by Eq. (5). After cessation of the pump, we extract the spin configuration and construct the corresponding spectral polynomial via the Lax vector Yuzbashyan et al. (2005a, b) (also see Appendix A). The spectral polynomial is a conserved quantity under time evolution with the [unperturbed, ] BCS Hamiltonian; its roots (“Lax roots”) are therefore also conserved. The isolated Lax roots occur in complex conjugate pairs, and encode the nature of the nonequilibrium state, and this is a robust, topological classification scheme Yuzbashyan et al. (2006). In particular, interaction quenches that produce phases I, II, and III (described above) respectively correspond to spectral polynomials with zero, one, or two pairs of isolated roots Yuzbashyan and Dzero (2006); Barankov and Levitov (2006b). This scheme also applies to topological -wave superfluids Foster et al. (2013) and strongly paired BCS-BEC -wave fermionic condensates Yuzbashyan et al. (2015). The isolated roots exactly determine in the limit. For phases II and I, , where is a nonzero constant (zero) in phase II (I). In phase II, is equal to the modulus of the imaginary part of the single isolated root pair.

For the pump-pulse quenches, we found that the spectral polynomial evaluated in terms of the post-pump Anderson pseudospin texture exhibits one or zero pairs of isolated roots, depending upon the degree of deformation from the initial BCS ground state. These results are consistent with the real time numerical simulations. The spin textures and the roots of the spectral polynomial are demonstrated for small system sizes in Fig 1. A quench phase diagram is constructed using much larger systems in Fig. 2. In phase II, only weakly depends on the value of (see Fig 2); here is the order parameter of the ground state, and is the tight-binding parameter of the square lattice. This result suggests that Eq. (9) is a reasonable approximation for capturing the nonlinear light-superconductor coupling. In Fig. 3, the dependence of on the polarization angle of the pump-pulse electric field is shown.

We have confirmed the isolated root predictions for phase II numerically up to with fourth order Runge-Kutta dynamics. ( is the linear size of the square lattice.) For sufficiently large pump energies, we find , consistent with phase I dynamics. The time evolution of the order parameter amplitudes in phases II and I is shown in Fig. 4.

ii.3 Optical conductivity in phases I and II

Besides the mechanism of the pump pulse quench, we have investigated the optical conductivity in phases II and I. The steady-state optical conductivity in the prethermalization plateau can be determined by the absorption or reflection of the probe pulse. Real low-temperature superconductors typically reside in the dirty limit, wherein . Here is the lifetime due to elastic impurity scattering. We reformulate the Mattis-Bardeen formulae Mattis and Bardeen (1958) for the real and imaginary parts of the optical conductivity in dirty superconductors, adapting them to steady-state nonequilibrium superconductivity.

For phases II and I, the generalized Mattis-Bardeen formulae depend on , but also on the full nonequilibrium quasiparticle distribution function . Here we employ the exact result for previously determined Yuzbashyan and Dzero (2006) for interaction quenches. The optical conductivity for the pump-probe quench will be qualitatively the same within each dynamical phase.

For a small phase II quench, the order parameter approaches to in the long time limit. The real part of the optical conductivity exhibits a cusp at frequency , but spectral weight extends down to zero (filling the gap) frequency for any nonzero quench. This is very similar to an equilibrium superconductor at finite temperature , despite the fact that our system is described by a pure nonequilibrium BCS state at all times. In fact, we find that the steady-state Mattis-Bardeen formulae for phase II are identical to the equilibrium finite-temperature results with replaced by and the thermal quasiparticle distribution replaced by the nonequilibrium one . The real part of the Mattis-Bardeen optical conductivity for phase II interaction quenches is presented in Fig. 5. The imaginary part of the optical conductivity exhibits a power law decay. The prefactor of the imaginary part depends on . Results are shown in Fig. 6. These Mattis-Bardeen results for nonequlibrium superconductors in phase II are consistent with the experimental observations via optical measurements Matsunaga et al. (2013); Matsunaga and Shimano (2015).

For phase I (strong deformation), the optical conductivity becomes indistinguishable from the normal metal. Namely, and . The approach to the phase II-I boundary from within phase II is qualitatively similar to an equilibrium superconductor with , see Figs. 5 and 6. Although vanishes, it is known that phase I is not equivalent to a zero temperature metallic state or any thermal ensemble; instead, it is a pure state with anomalous coherences, e.g., . The order parameter vanishes due to dephasing between the different momenta contributing to Eq. (4). We also show that the superfluid density (defined via the London penetration depth) vanishes in phase I; see Appendix D. Our results for the Mattis-Bardeen conductivity and superfluid density indicate that a different type of measurement (e.g., momentum-resolved) is required to reveal the anomalous coherences that persist in phase I.

Iii Pump: Twisting Anderson pseudospins

In the THz pump-probe experiments Matsunaga et al. (2013); Matsunaga and Shimano (2015), a pump pulse is injected into a superconductor (with ground state BCS gap 0.5 THz) at the beginning. If we consider a system at zero temperature, then the superconducting state is deformed from the ground state by the pulse. After the light exposure, the new state evolves under the original BCS Hamiltonian, given by Eq. (3), with . In this section, the effect of the pump pulse is discussed in detail. The coupling between light to the superconducting order parameter can only arise in the nonlinear order Varma (2002); Pekker and Varma (2015). We also discuss the connection between and the twisting of the Anderson pseudospin texture by the pump pulse.

iii.1 Light-superconductor coupling

The electric field of the pump pulse in the self-consistent Hamiltonian given by Eq. (3) is given by . We choose a Gaussian envelope,


for , and zero outside of this window. Here, is the unit polarization vector and is the peak amplitude of the vector potential. The detailed shape of the pulse only modifies the results quantitatively. We select a particular pulse profile to simulate the experimental setup in Ref. Matsunaga et al., 2013. Different from the previous numerical Bogoliubov approaches Papenkort et al. (2007, 2008); Krull et al. (2014), we assume that the linearly polarized THz pump pulse excitation acts on the system uniformly due to the large wavelength of THz light ( m). If the initial state is a parity-symmetric BCS ground state, then we can rewrite the Anderson pseudospin Hamiltonian in Eq. (3a) using a different field Tsuji and Aoki (2015),


where we have expressed in terms of a symmetric combination that only admits even powers of . The linear coupling vanishes, as expected for a clean single-band superconductor Varma (2002); Pekker and Varma (2015). For a 2D dispersion with , we expand to the order,


For a pure isotropic quadratic dispersion, the effect of the light can be viewed as a time-dependent chemical potential shift that does not alter the amplitude of the order parameter. It is therefore crucial to adopt a dispersion with energy-dependent curvature. Throughout this paper we study a square lattice nearest-neighbor tight-binding model at quarter filling qua (). The dispersion relation is where is the nearest-neighbor hopping strength and is the lattice constant. Equation (9) becomes


where is the relative angle between the polarization axis of and the direction, and defined by Eq. (6) is the pump energy that characterizes the strength of light-superconductor coupling.

The function encodes the anisotropy of the dispersion curvature, which depends on the lattice structure. The effect of such coupling induces a momentum-dependent -component magnetic field for the Anderson pseudospins. The Anderson spins with different momenta precess differently under the light exposure. In all cases, the time-evolving order parameter reaches a steady value and , where is the ground state BCS gap. The results are insensitive to (See Fig. 2). This suggests that the quadratic coupling to captures the mechanism of the pump-pulse quench. The values of depend strongly on the polarization angle. In Fig. 3, polarization angles and give quantitatively different but qualitatively similar results. The angular dependence is not generic but is determined by the underlying lattice structure. We do not expect the same angular dependence in the THz pump-probe experiments Matsunaga et al. (2013); Matsunaga and Shimano (2015).

iii.2 Twisted spin configurations and

After the pump pulse, the real time dynamics of the superconducting order parameters in Fig. 4 exhibit behavior consistent with interaction quenches Barankov et al. (2004); Yuzbashyan et al. (2005a, b); Barankov and Levitov (2006b); Yuzbashyan et al. (2006); Yuzbashyan and Dzero (2006) in phases II () and phase I (). In Appendix A we construct the Lax vector and spectral polynomial that can be employed to determine the Lax roots of the post–pulse-quench Anderson pseudospin texture. We find two possible situations, zero (phase I) or one pairs (phase II) of isolated roots. For the phase II, the values of the isolated roots , where and correspond to the asymptotic values of the nonequilibrium chemical potential and amplitude of the order parameter in the long time limit.

In Fig. 1, the Anderson pseudospin textures with different pump energies are plotted for a small square lattice system at quarter-filling. The main effect of the pump pulse is to twist the Anderson pseudospins near the Fermi surface in the plane. We can mimic this twisting with a simple construction that leads to , discussed in Appendix. A. We have also confirmed that the isolated root predictions with the numerical simulations of the real-time pseudospin dynamics. Due to finite size effects, it is difficult to locate the phase boundary between phase II and phase I from the isolated roots alone. We numerically simulate systems with linear sizes , , and . The pump-pulse quench phase diagram is given in Fig. 2.

Some numerical evidence for vanishing in THz pumped self-consistent mean field dynamics was previously reported in Ref. Papenkort et al. (2009), although this was not interpreted in terms of phase I. Moreover, the nature of the excitation in that work was different, since the ultrashort pump pulse assumed in Papenkort et al. (2009) had center frequency greater than the ground state optical gap, . Our pump pulse in Eq. (7) has center frequency roughly equal to , and most of the spectral weight falls below (similar to the experiment Matsunaga et al. (2013)). Our idea is that by exciting the “Higgs dynamics” entirely through the band nonlinearity, we are driving Anderson pseudospins in a way that would minimize linear absorption in a real, non–mean-field superconductor and maximize coherence. Understanding pair-breaking and decoherence by incorporating integrability-breaking terms would be necessary to prove this argument, which is an important avenue for future work. Our identification of the pump-quench induced phase II-I boundary with the disappearance of the isolated Lax root pair unifies the interaction and pump-quenches within the Lax classification Yuzbashyan et al. (2006) of BCS dynamics.

Iv Probe: Optical conductivity

In this section we extend the Mattis-Bardeen formulae Mattis and Bardeen (1958) for the real and imaginary parts of the optical conductivity of a dirty superconductor to the phase-II and phase-I nonequilibrium steady states.

iv.1 Brief review: Mattis-Bardeen formula

We start with a brief review of the Mattis-Bardeen formula Mattis and Bardeen (1958) for a superconductor at zero temperature. In the presence of elastic impurity scattering, momentum conservation is violated but energy conservation is preserved. Instead of deriving the explicit expression of the current-current correlation function in the presence of disorder, Mattis and Bardeen computed the ratio of the ac conductivity of a superconductor to the constant conductivity of a normal metal. The zero-temperature Mattis-Bardeen formula for the ac conductivity is


where is the conductivity for the normal metal, and are the ground state coherence factors for single particle energy , , and is an infinitesimal positive number. Equation (12) assumes a constant density of states. This generalization can be understood as follows. In the absence of the translational invariance, one can replace the labels and in Eq. (43) by independent energy labels and . The matrix elements of the current operators are randomized, and absorbed into the normal state conductivity . This level of approximation neglects quantum interference corrections.

The real part of Eq. (12) shows a gap for , an essential signature of zero temperature superconductivity Tinkham (1996). In order to evaluate the imaginary part of Eq. (12), one has to regularize the formula by subtracting the normal metal contribution at zero frequency Mattis and Bardeen (1958). The result can be approximated by for ; is proportional to the superfluid density. Despite the simplifying assumptions made to obtain the Mattis-Bardeen formula, the finite temperature version is consistent with experimental observations for BCS superconductors Tinkham (1996).

iv.2 Weak quench: Phase II

The long-time steady state of phase II can be expressed in terms of time-dependent coherence factors. These satisfy the asymptotic Bogoliubov-de Gennes equation,


where , and and are the long time asymptotic values of the nonequilibrium chemical potential and the order parameter, respectively.

The state is a pure BCS product wavefunction at all times,


The time-dependent coherence factors are


where is the nonequilibrium quasiparticle distribution function (see Appendix C) and . The time-independent coherence factors are


On can compute the two-time paramagnetic current polarization function for any BCS state of the form in Eq. (14); the result is given by Eq. (B). The current polarization function corresponding to Eq. (14) is a function separately of and . For the optical conductivity in the long-time asymptotic state, we only retain terms that can be expressed as functions of ; terms that oscillate with average out. The Mattis-Bardeen version of the ac conductivity then follows,




and is the nonequilibrium quasiparticle distribution function. The energy () integrations each run over the real line; the chemical potential can be set equal to zero for the constant (disorder-averaged) density of states.

The Mattis-Bardeen conductivity consists of two components. contains “ground state Cooper pair” and “excited state Cooper pair” contributions proportional to and , respectively. The real part of is zero for , while contributes to all frequencies. In fact, Eq. (17) turns out to be identical to the form obtained in thermal equilibrium at temperature Mattis and Bardeen (1958) if we replace and the thermal quasiparticle distribution function by the nonequilibrium one ; see Eqs. (49)–(52).

The real and imaginary parts of the ac conductivity for asymptotic phase II states are summarized in Figs. 5 and 6. The detailed calculations are discussed in Appendix C. Since the results depend upon both and the full nonequilibrium distribution function , we employ the exact analytical result for the latter that applies to instantaneous interaction quenches. We expect that the results are qualitatively unchanged for the THz pump quench. For interaction strength quenches, it is useful to introduce the quench parameter to locate the position in the phase diagram. This is defined via


where () is the interaction strength of the pre-quench (post-quench) Hamiltonian. For a 2D particle-hole symmetric superconductor with a uniform density of states, is the critical point that separates phase I () and phase II Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006). The real part of the ac conductivity is very similar to that of a finite temperature superconductor with optical gap replaced by . With a quench parameter , a peak arises at zero frequency, identical to the equilibrium situation with . When the quench strength approaches , the ac conductivity converges to the normal metal result. The imaginary part of the conductivity exhibits behavior for . We fit the imaginary part to


In Fig. 6, the prefactor as a function of the quench parameter is plotted. For , (which is the equilibrium result for ). deviates from as approaches ; the exact behavior depends upon the explicit form of the nonequilibrium quasiparticle distribution function .

iv.3 Strong quench: Phase I

For strong quenches, the order parameter vanishes to zero dynamically. In Eq. (16), reduces to . Moreover, the time-dependent coherence factors in Eqs. (15a) and (15b) become


The Mattis-Bardeen conductivity for phase I is


where we have used the same short-hand notations as in Eq. (17).

Equation (IV.3) depends only on the nonequilibrium distribution function . For an interaction quench, this is known explicitly and takes the same form as in phase II Yuzbashyan and Dzero (2006). One can show that the real part of Eq. (IV.3) is equal to and the imaginary part vanishes, so long as decays faster than in the limit . This is true for all phase I interaction quenches, including the quench to zero pairing strength. In the latter case the result can be obtained analytically. Since the pump-pulse quench induces phase I mainly by twisting pseudospins near the Fermi surface, and does so at much lower energy densities than an interaction quench (demonstrated in Sec. V, below), we unfortunately conclude that the optical conductivity will not distinguish a thermalized normal metallic state from the quantum coherent phase I of gapless superconductivity.

iv.4 Relation to experiments

In the pump-probe experiments Matsunaga et al. (2013); Matsunaga and Shimano (2015), the value of the superconducting order parameter was extracted from the kink in the real part of the measured optical conductivity. Our results (Fig. 5) confirmed that this interpretation is consistent for the coherent nonequilibrium steady states discussed here. We only compute the long time asymptotic behavior in this work, but the full two-time paramagnetic current polarization function necessary to monitor linear response at all times can be determined using the Mattis-Bardeen version of Eq. (B). We have neglected inelastic processes such as electron–acoustic-phonon scattering that are ultimately responsible for thermalization.

V Internal energy of nonequilibrium states

Figure 7: The -ratio for the pump-pulse quench. We consider the square lattice with linear system size at quarter filling. The red dashed line marks . The phase II-I boundary is at () for (). In particular, the ratio at the phase boundary () for . All the phase I states show . The polarization angle .
Figure 8: The -ratio for the interaction quench. The blue circles (green squares) correspond to a quarter-filled square lattice with (). is the strength of the nearest-neighbor hopping. The blue line (green line) corresponds to a 2D particle-hole symmetric model with (). In the latter, the single particle energies are uniformly distributed over to . The red dashed line marks . The phase II-I boundary is at () for the quarter-filled square lattice (2D particle-hole symmetric model). All the data plotted here belong to phase II, and there is a sizable range with . All the phase I states show in the interaction quench.

So far, we have only considered the post-pump dynamics of an integrable model for nonequilibrium superconductivity, i.e., time evolution according to the reduced BCS Hamiltonian in Eq. (1). Real superconductors contain many integrability-breaking perturbations such as inelastic electron–optical-phonon and electron–acoustic-phonon scattering, as well as Cooper pair breaking due to residual quasiparticle interactions. Thermalization to an equilibrium state always occurs on sufficiently long time scales; the novelty of the ultrafast “pump quench” in Matsunaga et al. (2013) is the ability to probe many-body dynamics on time scales shorter than this (i.e., within the pre-thermalization plateau). Of recent interest Basko et al. (2006); Nandkishore and Huse (2015); Schreiber et al. (2015) are two different thermalization schemes: (a) the system thermalizes as a generic (nonintegrable), but closed system, and (b) the system equilibrates with the external heat bath.

Here we pose the following interesting thermalization scenario. Suppose the pump-quench–induced nonequilibrium superconductor thermalizes to an equilibrium state before the external bath acts to absorb the excess heat. Since the pump dumps a large amount of energy into the system, the thermalized state can be either an equilibrium superconductor with , or a normal metal. We can distinguish these alternatives by comparing the injected internal energy of the nonequilibrium superconducting pure state to the internal energy of the equilibrium superconductor at the critical temperature, relative to the ground state (zero temperature) energy. We define the “-ratio” as follows,


where denotes the expectation of the post-quench Hamiltonian in the pumped state, is the ground state energy of the post-quench Hamiltonian, and is the thermal average of the post-quench Hamiltonian (internal energy density) at the critical temperature. The post-quench Hamiltonian corresponds to Eq. (1) with . If the nonequilibrium state carries more energy than the energy of a critical superconductor (), the system will thermalize to a normal state. Otherwise, it will thermalize to a finite temperature superconductor ().

For the pump pulse quench, we plot the -ratio as a function of the pump energy in Fig. 7. We found that all phase I states show . This indicates that thermalization in a closed system will not produce a superconducting equilibrium state. Some of the phase II states close to the phase boundary also show . The region decreases for smaller . Interestingly, the phase II-I boundary shows for the smallest studied .

In order to better understand the internal energy dependence of the nonequilibrium phases, we also compute the -ratio for two types of interaction quenches: (a) 2D square lattice with quarter filling and (b) 2D uniform particle-hole symmetric model with a finite bandwidth . Here, denotes the expectation value of the post-quench Hamiltonian in the pre-quench ground state. The results are shown in Fig. 8. Although interaction quenches (a) and (b) give nearly identical values of throughout most of phase II, both exceed far from the phase II-phase I boundary. In other words, the injected internal energy of interaction quenches approaching phase I is much larger than for the pump pulse quench. It means that the internal energy alone is not a good indicator for the dynamical boundaries between the nonequilibrium phases. This is not surprising, because the integrable nonequilibrium BCS dynamics are constrained by an infinite number of conserved quantities.

The pump-pulse quench reaches phase I at much lower energy densities than the interaction quench. This is because the former strongly scrambles the Anderson pseudospin texture for states near the Fermi energy, but barely modifies those away from it; see Fig. 1. On the other hand, the interaction quenches affect all pseudospins regardless of the single particle energy. The interaction quench injects more energy than minimally required to generate nontrivial nonequilibrium dynamics. Due to the lower internal energy, one should expect that the window for observing nonequilibrium dynamics is larger in the pump pulse versus interaction quench protocols. In addition, the minimal required -ratio might encode information about the phase II-I boundary. It is possible that at the phase II-I boundary for the optimized initial state(s) without any “redundant” internal energy.

Vi Discussion and conclusion

In this paper, inspired by the experiments in Refs. Matsunaga et al. (2013); Matsunaga and Shimano (2015) we have studied quantum quench dynamics following the exposure of a BCS superconductor to an intense, approximately monocycle THz pulse. We connected the “pump-pulse quench” protocol to the interaction quench protocol using the Lax reduction method Yuzbashyan et al. (2005a, b). We constructed the quench phase diagram in Fig. 2. Ultrafast spectroscopy with intense THz sources opens a new vista for solid state physics: the study of coherent, many-body quantum dynamics.

Due to the large wavelength ( m) of the THz light, we assumed that the pump-pulse light couples to the entire system uniformly in space. Such coupling can be described in the Anderson pseudospin language. We focused on the second order (in vector potential amplitude) coupling, which is a curvature-dependent -component effective magnetic field. We simulated the pump-pulse quench on a square lattice model at quarter-filling. We showed the possibility of accessing a gapless phase of fluctuating superconductivity (phase I) by increasing the pump energy.

To connect the pump-pulse quench to the existing experiments Matsunaga and Shimano (2012); Matsunaga et al. (2013); Matsunaga and Shimano (2015), we estimate the parameters including bandwidth, superconducting gap, and pump energy. The bandwidth of NbN is roughly 9.53eV Mattheiss (1972). One can estimate the tight binding parameter eV. The optical superconducting gap in NbN is meV. This gives . In Ref. Matsunaga et al., 2013, the reported pump energies vary from for generating nonequilibrium superconductivity in phase II. The phase diagram in Fig. 2 shows weak dependence upon the value of . Our predicted phase boundary for the quarter-filled 2D square lattice model is with the polarization angle . We expect that phase I can be generated with the comparable magnitude of pump energies in the THz pump-probe experiments.

Optical conductivity is the quantity measured in Refs. Matsunaga and Shimano (2012); Matsunaga et al. (2013); Matsunaga and Shimano (2015). We computed the Mattis-Bardeen optical conductivity for dirty superconductors in phases II and I. The phase II result is the same as the prediction for an equilibrium superconductor with , except that it depends on the asymptotic nonequilibrium order parameter and the nonequilibrium quasiparticle distribution function. These results back up the gap-extraction procedures employed in the existing experiments Matsunaga et al. (2013); Matsunaga and Shimano (2015). The phase I optical conductivity cannot be distinguished from a normal metal. In Appendix D, we show that the superfluid density (defined via the London penetration depth) also vanishes in phase I.

We have assumed integrable post-pump dynamics throughout this work. Integrability-breaking processes induce thermalization and limit the time window (pre-thermalization plateau) for observing quench dynamics. Moreover, phonons Matsunaga et al. (2013); Krull et al. (2014); Kemper et al. (2015); Murakami et al. (2016); Sentef et al. (2016) and other coexisting collective modes Akbari et al. (2013); Moor et al. (2014); Fu et al. (2014); Cea and Benfatto (2014); Dzero et al. (2015); Krull et al. (2016); Sentef et al. (2017) might modify the evolution in a superconductor. A detailed study of thermalization and decoherence out of the far-from-equilibrium superconductor states discussed here remains a challenging project for future work. Based on the results in Sec. V, the pump-pulse quench is a much more efficient protocol for accessing nonlinear nonequilibrium dynamics, compared to interaction quenches. It would be interesting to try to reach other dynamical phases (e.g., the quench-generated Floquet phase III) by engineering different ultrafast excitation schemes. One possibility would be to design a sequence of pump pulses of variable duration and delays. Realizing the pump-pulse quench protocol in an ultracold fermion condensate might also yield more understanding for THz pump-probe experiments.


We thank Maxim Dzero, Junichiro Kono, Ryo Shimano, and Emil Yuzbashyan for helpful discussions. This research was supported by the Welch Foundation grant No. C-1809 and by NSF CAREER grant no. DMR-1552327 (Y.-Z.C., Y.L., and M.S.F.). This work was also supported by NSF grant no. DMR-1001240, through the Simons Investigator award of the Simons Foundation, and by a Center for Theoretical Quantum Matter postdoctoral fellowship at University of Colorado Boulder (Y.-Z.C.).

Appendix A Lax construction and spectral polynomial; twisted pseudospins

For a classical system consisted of pseudospins, we define the Lax vector components Yuzbashyan et al. (2005a),


where is the strength of the BCS interaction in Eq. (1) and denotes a complex parameter. The Lax norm is defined as


We impose the canonical Poisson bracket relation for the pseudospins,


Based on the above condition, one can show that the Lax vector components satisfy


We define the spectral polynomial


is a polynomial of degree in the parameter ; it is conserved under evolution by the BCS Hamiltonian in Eq. (1) with . The spectral polynomial for the BCS ground state can be written as


There are real roots and 2 complex roots in . The two complex roots (isolated roots) encode the ground state BCS gap ,


For interaction quenches, it has been shown that most of the roots of the spectral polynomial lie along the real axis in the thermodynamic limit. The complex roots (isolated roots) determine the asymptotic long time behavior following the quench. The procedure used to extract this behavior is called Lax reduction Yuzbashyan et al. (2006). Interaction quenches from a BCS initial state exhibit 0, 1, or 2 pairs of isolated roots, associated to exactly solvable 0, 1, or 2 collective-pseudospin problems Yuzbashyan et al. (2006); Barankov and Levitov (2006b); Yuzbashyan and Dzero (2006). These respectively correspond to phases I, II, and III discussed in Sec. II.1. In phase II, is encoded in the single isolated root pair , where is the effective nonequilibrium chemical potential.

As demonstrated by the Anderson pseudospin textures exhibited in Fig. 1, the main effect of the THz pump pulse is to twist pseudospins along the Fermi surface in the pseudospin plane. This “scrambling” induces an asymptotic which is strictly less than the ground state gap , as shown in the dynamical phase diagram Fig. 2.

We can mimic the pump-induced twisting with the following simple construction. We consider a situation such that the - and -components of the pseudospin texture are deformed relative to the BCS ground state configuration. We write down the Lax vector as a sum over equal energy sectors,


where indicates the distinct energy level, is the spin with energy , and is the net polarization of the spins with energy .

The Lax vector for the twisted configuration is defined via


where characterizes the degree of twisting in the sector of single energy . For simplicity we take (independent of ). Then the spectral polynomial for the twisted texture can be written as


where is the same degree polynomial that appears in the ground state [Eq. (32)]. The isolated roots give . The asymptotic value of the order parameter is therefore . This suggests that one can achieve phase I (effective 0 spin problem with ) by increasing the degree of the twisting.

Appendix B Linear response theory and Kubo formula

We obtain the general expression for the Kubo conductivity associated to a generic, time-evolving pure BCS product state in this section. In the continuum limit, the paramagnetic current operator is expressed as


where is the effective mass of electron. Linear response theory gives




is the retarded paramagnetic current polarization function. is the vector potential for a monochromatic source with frequency . The expectation is performed in terms of the initial BCS ground state, while the Heisenberg picture current operators incorporate the effects of the pump pulse as well as the subsequent free time-evolution under the BCS Hamiltonian.

In the ground state, is translational invariant in time. The retarded polarization function in the frequency-momentum space is Mahan (1981)


The optical conductivity is


where the second term is the diamagnetic contribution. Since , the real part of the optical conductivity is exactly zero for a clean (single band) superconductor. Low-temperature superconductors typically reside in the dirty limit due to the presence of quenched disorder. Relaxing momentum conservation in Eq. (43) leads to the Mattis-Bardeen formulae Mattis and Bardeen (1958).

For the BCS wavefunction in Eq. (14) with generic time-dependent coherence factors, the polarization function can be evaluated. The result is