The theory of parametrically amplified electron-phonon superconductivity

TL;DR

This paper develops a comprehensive theoretical framework for light-induced superconductivity in driven electron-phonon systems, highlighting the role of lattice nonlinearities and parametric amplification in enhancing electron-phonon interactions.

Contribution

It introduces a Floquet generalization of Migdal-Eliashberg theory and derives quantum Floquet-Boltzmann equations to study nonequilibrium dynamics of driven superconductors.

Findings

Enhanced effective electron-phonon couplings via parametric amplification.

Identification of conditions leading to transient superconductivity.

Prediction of electronic Floquet band formation observable in experiments.

Abstract

The ultrafast optical manipulation of ordered phases in strongly correlated materials is a topic of significant theoretical, experimental, and technological interest. Inspired by a recent experiment on light-induced superconductivity in fullerenes [Mitrano et al., Nature 530, 2016], we develop a comprehensive theory of light-induced superconductivity in driven electron-phonon systems with lattice nonlinearities. In analogy with the operation of parametric amplifiers, we show how the interplay between the external drive and lattice nonlinearities lead to significantly enhanced effective electron-phonon couplings. We provide a detailed and unbiased study of the nonequilibrium dynamics of the driven system using the real-time Green's function technique. To this end, we develop a Floquet generalization of the Migdal-Eliashberg theory and derive a numerically tractable set of quantum Floquet-Boltzmann kinetic equations for the coupled electron-phonon system. We study the role of parametric phonon generation and electronic heating in destroying the transient superconducting state. Finally, we predict the transient formation of electronic Floquet bands in time- and angle-resolved photo-emission spectroscopy experiments as a consequence of the proposed mechanism.

Figures

Equations

$$\mathcal{L}[\phi ,\Psi ](t)=\sum _{\mathbf{k}}{\Psi }_{\mathbf{k}}^{†}\left(i{\partial }_t\mathbb{I}-{\xi }_{\mathbf{k}}{\hat{\sigma }}_3\right){\Psi }_{\mathbf{k}}^{}$$

$$-\frac{1}{2}\sum _{\mathbf{q}}\frac{1}{2{\omega }_{\mathbf{q}}}{\phi }_{\mathbf{q}}^{}\left({\partial }_t^2+{\omega }_{\mathbf{q}}^2\right){\phi }_{-\mathbf{q}}^{}-\sum _{j\in \mathrm{lattice}}{\mathcal{V}}^{\mathrm{ph}}({\phi }_j)$$

$$-\frac{1}{\sqrt{N}}\sum _{\mathbf{k},{\mathbf{k}}^{\prime }}{g}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{}{\phi }_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{}{\Psi }_{{\mathbf{k}}^{\prime }}^{†}{\hat{\sigma }}_3{\Psi }_{\mathbf{k}}^{}$$

$$+\frac{\Lambda }{2}|F(t){|}^2\sum _{j\in \mathrm{lattice}}{\phi }_j$$

$${\mathcal{V}}^{\mathrm{ph}}(\phi )=-\frac{{\kappa }_3}{3!}{\phi }^3-\frac{{\kappa }_4}{4!}{\phi }^4,$$

$$F(t)=F_{\mathrm{env}}(t)\cos ({\Omega }_{\mathrm{drv}}t),$$

$${\partial }_t^2⟨{\hat{\phi }}_0(t)⟩+{\omega }_0^2⟨{\hat{\phi }}_0(t)⟩=\frac{\Lambda {\omega }_0}{2}\sqrt{N}{F}_{\mathrm{env}}^2(t){\cos }^2({\Omega }_{\mathrm{drv}}t),$$

$$\begin{gathered} ⟨{\hat{\phi }}_0(t)⟩\approx \frac{\Lambda \sqrt{N}}{4{\omega }_0}{F}_{\mathrm{env}}^2(t)+ \\ \frac{\Lambda \sqrt{N}{\omega }_0}{4({\omega }_0^2-4{\Omega }_{\mathrm{drv}}^2)}{F}_{\mathrm{env}}^2(t)\cos (2{\Omega }_{\mathrm{drv}}t). \end{gathered}$$

$$\begin{gathered} -\frac{{\kappa }_3}{3!}\sum _{j\in \mathrm{lattice}}{\hat{\phi }}_j^3\to \\ -\frac{{\kappa }_3}{2}\left[{\phi }_0+{\phi }_1\cos (2{\Omega }_{\mathrm{drv}}t)\right]\sum _{\mathbf{q}\ne 0}{\hat{\phi }}_{\mathbf{q}}{\hat{\phi }}_{-\mathbf{q}}. \end{gathered}$$

$$\begin{gathered} -\frac{{\kappa }_4}{4!}\sum _{j\in \mathrm{lattice}}{\hat{\phi }}_j^4\to \\ -\frac{{\kappa }_4}{4}{\left[{\phi }_0+{\phi }_1\cos (2{\Omega }_{\mathrm{drv}}t)\right]}^2\sum _{\mathbf{q}\ne 0}{\hat{\phi }}_{\mathbf{q}}{\hat{\phi }}_{-\mathbf{q}}. \end{gathered}$$

$$H_{\mathrm{p},\pm \mathbf{q}}(t)$$

$$-{\kappa }_3{\phi }_1\cos (2{\Omega }_{\mathrm{drv}}t){\hat{\phi }}_{\mathbf{q}}{\hat{\phi }}_{-\mathbf{q}}$$

$$H_{\mathrm{ep},\pm \mathbf{q}}$$

$${\hat{Q}}_{+}$$

$${\hat{Q}}_{-}$$

$${\hat{P}}_{+}$$

$${\hat{P}}_{-}$$

$$H_{\mathrm{p},\mathbf{q}}$$

$$H_{\mathrm{ep},\mathbf{q}}$$

$${\Omega }_{\mathbf{q}}^2(t)$$

$$\alpha$$

$$\Delta S^{(2)}(t)=-\frac{|{\tilde{g}}_{\mathbf{q}}{|}^2}{\hslash }{\int }_{-\infty }^t\mathrm{d}{t}^{\prime }{e}^{-ϵ(t-t^{\prime })}{\mathcal{D}}_{QQ}^R(t,t^{\prime })\hat{\rho }(t)\hat{\rho }(t^{\prime }),$$

$$U(t)=\frac{|{\tilde{g}}_{\mathbf{q}}{|}^2}{\hslash }{\int }_{-\infty }^t\mathrm{d}{t}^{\prime }{e}^{-ϵ(t-t^{\prime })}{\mathcal{D}}_{QQ}^R(t,t^{\prime }).$$

$$\frac{\mathrm{d}\hat{Q}(t)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\hat{P}(t)}{\mathrm{d}t}$$

$$\hat{Q}(t)$$

$$\hat{P}(t)$$

$${\mathcal{D}}_{QQ}^R(t,t^{\prime })$$

$$=-\frac{\hslash }{M{\Omega }_{\mathrm{drv}}}\theta (t-t^{\prime }){\mathfrak{M}}_{QP}(t,t^{\prime }).$$

$$\underset{\alpha \to 0}{\lim }{\mathcal{D}}_{QQ}^R(t,t^{\prime })=\hslash \theta (t-t^{\prime })\frac{\sin {\omega }_{\mathbf{q}}(t-t^{\prime })}{M{\omega }_{\mathbf{q}}}.$$

$$\begin{gathered} \frac{U(t)}{U_{\mathrm{eq}}}=1-\frac{2\alpha {\omega }_{\mathbf{q}}^2\cos (2{\Omega }_{\mathrm{drv}}t)}{{\omega }_{\mathbf{q}}^2-4{\Omega }_{\mathrm{drv}}^2} \\ +\frac{2{\alpha }^2{\omega }_{\mathbf{q}}^2\left[{\omega }_{\mathbf{q}}^2-16{\Omega }_{\mathrm{drv}}^2+{\omega }_{\mathbf{q}}^2\cos (4{\Omega }_{\mathrm{drv}}t)\right]}{({\omega }_{\mathbf{q}}^2-16{\Omega }_{\mathrm{drv}}^2)({\omega }_{\mathbf{q}}^2-4{\Omega }_{\mathrm{drv}}^2)}+\mathcal{O}({\alpha }^4) \end{gathered}$$

$$⟨\frac{U(t)}{U_{\mathrm{eq}}}⟩$$

$$\mathrm{Var}\left[\frac{U(t)}{U_{\mathrm{eq}}}\right]$$

$$U(t)\approx U_0+U_1\cos (2{\Omega }_{\mathrm{drv}}t),$$

$$U_0$$

$$U_1$$

$$\left(1-U_0{F}_n\right){\Delta }_n+\frac{U_1}{2}{F}_n\left({\Delta }_{n-1}+{\Delta }_{n+1}\right)=0,$$

Full text

I Introduction

In the recent years, the rapid progress of the field of ultrafast pump-probe spectroscopy experiments has enabled an unprecedented exploration of many-body quantum dynamics in far-from-equilibrium states (for reviews, see Refs. Basov et al. (2011); Giannetti et al. (2016)). The application of strong ultrafast laser pulses can dramatically alter the equilibrium state, outreach the linear response regime, and enable the induction of novel ordered states and stimulation of phase transitions via transient nonequilibrium states.

One of the main motivations behind these experiments is to shed light on the interplay between competing orders in strongly-correlated superconductors, along with the tantalizing outlook of stimulating the superconducting transition at temperatures above the critical temperature $T_{c}$. An early experimental evidence proving the possibility of stimulating superconductivity via external fields is the Wyatt-Dayem effect Wyatt et al. (1966); Dayem and Wiegand (1967), where microwave radiation of superconducting micro-bridges in the MHz to GHz frequency range was found to increase $T_{c}$ by a few percents. This effect was explained theoretically by Eliashberg Eliashberg (1970) on the basis of the nonequilibrium shift of the quasiparticle occupation to high energies. Subsequent theoretical work Chang and Scalapino (1978) and experiments in double-barrier tunnel junctions and strips Pals and Dobben (1979); Blamire et al. (1991); Heslinga and Klapwijk (1993); Komissinski and Ovsyannikov (1996) found a much larger effect up to several times larger than the equilibrium $T_{c}$. An experimental proposal for investigating this effect using ultracold fermionic atoms has also been given Robertson and Galitski (2009).

Recently, Mitrano et al. Mitrano et al. (2016) have reported a large increase in carrier mobility and the opening of an optical gap upon stimulating the intercalated fullerene superconductor $K_{3}C_{60}$ with a femto-second mid-infrared light pulse in the frequency range $80$–$200$ meV ($19$–$48$ THz). These effects persist for several pico-seconds after pumping, and remarkably for initial temperatures up to $T_{i}\sim 100$ K, much higher than the equilibrium $T_{c}\sim 20$ K, providing a compelling evidence for a possibly light-induced superconducting state. The experimentally observed resonance with several $C_{60}$ vibrational modes suggests that the underlying mechanism for enhanced Cooper pairing in this experiment stems from lattice distortions and is distinct from the Wyatt-Dayem effect.

The application of a strong pump pulse alters the initial equilibrium state in various ways and is a complex function of the strength of the drive, strength of coupling to different degrees of freedom and energetic proximity to resonances. The accurate theoretical modeling of light-stimulated superconductivity in $K_{3}C_{60}$ is exacerbated by the structural complexity of $K_{3}C_{60}$, including the three-fold degeneracy of the conduction $t_{1u}$ bands and their coupling to eight $H_{g}$ intra-molecular Jahn-Teller phonons Gunnarsson and Fullerides (2004), strong electron-phonon coupling $\lambda\sim 0.5$–$1$, narrow conduction band $\omega_{\mathrm{ph}}/W_{\mathrm{el}}\sim 0.1$–$0.25$ ($\omega_{\mathrm{ph}}$ is the typical phonon energy scale and $W_{\mathrm{el}}$ is the conduction bandwidth), and strong Coulomb interaction $U_{c}/W_{\mathrm{el}}\sim 1.5$–$2.5$ Gunnarsson (1997). A reliable material-specific analysis must combine state-of-the-art ab initio modeling including nonlinear interactions and beyond-Migdal vertex corrections Grimaldi et al. (1995) into the framework of nonequilibrium field theory. We do not pursue this formidable goal here; rather, inspired by the experiment and embracing the Occam’s razor tradition, we explore a simplified model with fairly generic electron-phonon interaction which retains the essential features observed in experiments with light-stimulated superconductivity in $K_{3}C_{60}$.

At the simplest level, the pumping pulse with frequency $\Omega_{\mathrm{drv}}\sim 100$ meV strongly drives near-resonant infra-red (IR) active lattice vibrational mode, such as $T_{1u}$ modes in fullerenes. As a first approximation, one may average out fast oscillations at the scale of $\Omega_{\mathrm{drv}}^{-1}$. The presence of lattice anharmonicities and nonlinear coupling between vibrational modes results in the expansion and contraction of time-averaged lattice constants and electronic orbital configurations. The time-averaged electron-phonon coupling constants and electronic density of states are consequently renormalized. This approach is adopted in Ref. Mitrano et al. (2016) where an ab initio analysis in the static “frozen-phonon” approximation is performed and it is shown that time-averaged lattice deformations exhibit a favorable trend toward increasing $T_{c}$.

In this paper, we aim to show that the implications of a strong periodic drive and its interplay with lattice nonlinearities reaches beyond statically renormalized model parameters. In fact, we will show that the mechanism which yields the most striking enhancement of electron-phonon coupling is purely dynamical in nature and is not described by time-averaged Hamiltonians.

The phonon-mediated electron-electron attraction $U$ is usually understood using second-order perturbation theory: an electron distorts the lattice and the other electron is attracted to the lattice distortion, see Fig. 1 (left). In other words, this attractive potential is proportional to the retarded phonon response function. We will show that the enhancement of superconductivity in a driven nonlinear lattice is conceptually similar to the operation of a parametric amplifier circuit: the “nonlinear capacitor” is realized by the lattice nonlinearity, the “ac pump source” is realized by an excited lattice vibrational mode, the “input signal” is the phonon excitation caused by a momentum kick from an electron, and the “output signal” is the parametrically amplified phonon response observed by the other electron. In essence, lattice nonlinearities convert the coherent motion of the driven mode into a source of parametric drive for the phonon that couples to conduction electrons. When this drive is near parametric resonances, the retarded response will be significantly amplified, leading to a much stronger electron-electron attraction. Parametric driving also induces strong temporal oscillations in the effective electron-electron attraction, allowing it to visit very large values during each cycle. We will show that such temporal oscillations can significantly enhance $T_{c}$ even if the time-averaged attraction remains constant, see Figs. 1 (right) and 6.

A rigorous quantitative analysis of this simple mechanism and its consequences in a realistic electron-phonon model goes beyond the amplifier analogy as one must take into account several competing effects. Most importantly, the nearly resonant drive also results in parametric generation of high-energy phonons that dissipate their excess energy to electrons, leading to to higher scattering rates and heating. It is not a priori clear which subset of these phenomena prevails, even for short times, without resorting to an unbiased and rigorous framework. Ultimately, we find that without an external cooling mechanism, Cooper pairing may only be enhanced for a short time similar to the experiments, and the normal state takes over as high-energy phonon excitations equilibrate their energy with electrons.

Our goal in this paper is two-fold. First and foremost, we wish to present a transparent and physical analysis of the role of parametric resonances of the lattice in enhancing electron-phonon interactions and stimulating Cooper pair formation. The major part of this goal is achieved in the first part of the paper using perturbation theory, BCS theory, classical dynamics, along with a number of common-sense simplifications. Secondly, we aim to develop a rigorous theoretical formalism for analyzing the nonequilibrium dynamical nature of light-stimulated superconductivity experiments; a formalism that takes into account the detailed driven-dissipative evolution of phonons and heating of electrons while being flexible enough to include material-specific properties and paving the way for future investigations. To this end, we develop an extension of the Migdal-Eliashberg theory Migdal (1958); Eliashberg (1960) to periodically driven electron-phonon systems with lattice nonlinearities, and utilize it to substantiate the results of the first part as an immediate application.

The conventional Migdal-Eliashberg theory is a cornerstone of the modern theory of superconductivity, both for qualitative understandings and accurate ab initio calculations. The existing attempts at the real-time extension of the Migdal-Eliashberg theory are known to be intractably difficult to work with due to the complicated temporal structure of the equations Bennemann and Ketterson (2008); Sentef et al. (2016). Here, we combine ideas from effective actions, Floquet theory, dynamical mean-field theory, and quantum kinetic theory to develop a formalism that is well-suited for numerical and analytical studies of periodically driven systems. The Floquet quantum kinetic formalism trades fast drive-induced oscillations of nonequilibrium propagators with slowly-varying Floquet components, and memory convolution integrals with algebraic products along with derivative corrections Martinez (2003); Tsuji et al. (2008); Aoki et al. (2014); Genske and Rosch (2015). These controlled approximations effectively reduce the two-time Kadanoff-Baym integro-differential equations Kadanoff and Baym (1962) to (implicit) ordinary differential equations which are much easier to solve numerically. The extension of the quantum kinetic formalism to periodically driven systems has been considered before in Ref. Genske and Rosch (2015) in a different context and in the Boltzmann “quasiparticle” approximation. The latter is obtained by neglecting off-shell processes Danielewicz (1984); Berges (2004). We do not adopt this approximation here. As we pointed out earlier, suppression of electronic and phononic quasiparticle coherence is an important factor in the analysis of transient superconductivity. Hence, a detailed study of the changes in the spectral functions of Floquet quasiparticles will be an important ingredient of our theory.

The experimental observation of the light-induced superconducting state in $K_{3}C_{60}$ Mitrano et al. (2016) has inspired several theoretical works. Sentef et al. Sentef et al. (2016) have studied the transient dynamics of the superconducting gap following a change in the coupling constants. In an earlier work, we outlined the role of parametric driving in enhancing the electron-phonon coupling Knap et al. (2015) and analyzed the problem using a Floquet extension of the BCS theory. Komnik et al. Komnik and Thorwart (2016) have recently worked out a similar BCS framework using a more concise analytical approach. More recently, Kennes et al. Kennes et al. (2016) have suggested nonlinear electron-phonon couplings as another plausible source of enhancing $T_{c}$ in a highly-pumped state. Last but not least, Kim et al. Kim et al. (2016) have suggested light-induced changes in the screened Coulomb matrix elements as a factor for enhancing superconductivity in intercalated fullerenes. We would like to mention that none of these works, except for Ref. Knap et al. (2015), have studied the role of undesirable competing factors within their respective models. The electrons are always assumed to remain in the initial thermal state and heating is ignored. Given that superconductivity in fullerenes is mediated by high-frequency optical phonons, the issue of heating is a crucial aspect of the phenomenology even for short times. Another important goal of our paper to provide a first complete and unbiased analysis of the competition between processes that enhance and suppress Cooper pairing.

Finally, we would like to emphasize that the parametric amplification of electron-phonon coupling is not limited to enhancing Cooper pairing and is expected to find similar applications to other systems. For example, the same framework can be employed to study the recently observed enhancement of electron-phonon coupling in periodically distorted graphene Gierz et al. (2016) and driven opto-mechanical cavities Xu et al. (2015); Lemonde et al. (2016); Stehlik et al. (2016); Zeytinoglu et al. (2016).

I.1 Organization of the paper

This paper is organized as follows. We describe the model in Sec. II as the first step, and present its analysis in two separate stages. Before delving into the detailed formalism, we give a more intuitive account using perturbation theory, classical dynamics, and the BCS theory to demonstrate the idea of parametric amplification in nonlinear lattices and its implications in Sec. III. Many of the relevant details such as feedback to electrons, heating, and competing factors are left to the second stage.

Sec. IV and its multiple subsections are dedicated to developing the formalism of Floquet-Migdal-Eliashberg quantum kinetics. In particular, a pairing instability criterion is derived in Sec. V that generalizes the result of Scalapino, Schrieffer, and Wilkins Scalapino et al. (1966) to quasi-steady Floquet states. As a first application of the formalism, we study the stationary solutions of the driven-dissipative state of phonons while neglecting the heating of electrons. This allows us to gain insight about the parameter regimes of maximal electron-phonon coupling enhancement, and to study the role of individual factors in enhancing and suppressing Cooper pairing. We move on the fully non-equilibrium scenario in Sec. VI.2 where we discuss the dynamics of the coupled electron-phonon system and show that a window of transient superconducting instability can exist even if the heating of electrons are taken into account. Finally, we use our theory to make additional experimental predictions in Sec. VI.3, in particular, the dynamical formation of Floquet conduction bands which can be probed using time- and angle-resolved photo-emission spectroscopy (tr-ARPES).

Some of the technical details, in particular those pertaining to numerical methods, have been moved to the appendices. The appendices also include an extensive discussion of the role of electrons in generating phonon nonlinearities (see Sec. E). In particular, we show that the magnitude of electron-mediated phonon nonlinearities increases near parametric resonances and can make a significant contribution to intrinsic lattice nonlinearities.

Finally, Fig. 2 shows a summary of the non-equilibrium dynamics obtained from the Floquet-Migdal-Eliashberg formalism; refer to the figure caption for details.

II The model

We start with a general model for conduction electrons and a single phonon branch, along with an external drive that couples to the uniform lattice displacement, a local nonlinear lattice potential, and a linear electron-phonon coupling. The Lagrangian for this system is given as:

[TABLE]

Here, $\Psi^{\phantom{\dagger}}_{\mathbf{k}}=(c_{\mathbf{k}\uparrow},c_{-\mathbf{k}\downarrow}^{\dagger})^{T}$ is the Nambu spinor of the conduction electrons, $\xi_{\mathbf{k}}$ is the electron dispersion, $\varphi_{\mathbf{q}}\equiv b^{\dagger}_{\mathbf{q}}+b^{\phantom{\dagger}}_{-\mathbf{q}}=\sum_{j}e^{-i\mathbf{q}\cdot\mathbf{R}_{j}}\varphi_{j}/\sqrt{N}$ is the lattice displacement operator, $\omega_{\mathbf{q}}$ is the phonon dispersion, and ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ is the linear electron-phonon coupling constant. Furthermore, $\mathcal{V}^{\mathrm{ph}}(\varphi)$ is the local lattice anharmonic potential which, for low-amplitude deformations, can be modeled as:

[TABLE]

We assume $\kappa_{4}>0$ since the lattice potential generically softens for large deformations. The sign of $\kappa_{3}$ is inconsequential due to symmetries. We neglect Coulomb interaction to simplify the analysis. We will briefly comment on its effect later on and argue that it does not play a consequential role in the phenomenon that is the case here. Finally, $F(t)$ is the external classical drive that couples to the uniform $\mathbf{q}=0$ lattice displacement with strength $\Lambda/2$.

II.1 The origin of the drive term

The generic model we introduced in Eq. (II) is compatible with several scenarios suggested for modeling the role of the drive in pump-probe experiments of different materials. If $\varphi$ describes a polarizable (IR-active) phonon, $F(t)$ can be directly identified with the external electric field, in which case, the coupling strength $\Lambda$ will be proportional to the polarizability of $\varphi$. On the other hand, if $\varphi$ is a non-polarizable (Raman-active) phonon, even though the incident light does not directly influence it through dipole coupling, the classical drive term can still be obtained via nonlinear coupling to a driven “proxy” IR-active phonon. The leading order nonlinear IR/Raman coupling allowed by symmetries is the cubic $\propto\varphi_{\mathrm{IR}}^{2}\,\varphi$ interaction. In this case, we can identify $F(t)\sim\langle\varphi_{\mathrm{IR}}(t)\rangle$ as the coherent oscillations of the driven IR-active mode and $\Lambda$ as the strength of the cubic coupling to the Raman phonon $\varphi$. Regardless of the origin of the classical drive, we assume:

[TABLE]

where $\Omega_{\mathrm{drv}}$ is the principal frequency of the classical drive and $F_{\mathrm{env}}(t)$ is its slowly-varying envelope. Note that the classical drive couples to $\varphi$ in intensity, $|F(t)|^{2}$, such that the effective principal drive frequency is $2\Omega_{\mathrm{drv}}$.

If the drive term originates from nonlinear coupling to an IR-active phonon, $\Omega_{\mathrm{drv}}$ may no longer be identified with the frequency of the incident light after the pump pulse is ramped down; rather, the pump pulse coherently drives the proxy IR-active phonon out of its equilibrium position and subsequently, the coupled IR-active and Raman phonons oscillate together at a frequency predominantly determined by the IR-active mode. The proposed model still applies to this case with the appropriate choice of $\Omega_{\mathrm{drv}}$.

II.2 Different routes to parametric driving

For the purposes of this work, the necessary ingredient of the model is a mechanism to achieve parametric driving of the $\varphi$ phonon, i.e. a route for achieving the effective substitution $\omega_{\mathbf{q}}^{2}\rightarrow\omega_{\mathbf{q}}^{2}[1+2\alpha_{\mathbf{q}}\cos(2\Omega_{\mathrm{drv}}t)]$, where $\alpha_{\mathbf{q}}$ is the effective parametric driving amplitude. In the model proposed by Eq. (II), this is achieved from the interplay between the nonlinearities of the $\varphi$ phonon and its own coherent displacement, as further explained in the next section. There exists, however, a multitude of other physically realizable routes that all lead to parametric driving. This situation closely resembles the multitude of architectures proposed for building electronic parametric amplifiers over the years using elements such as variable capacitance diodes, nonlinear inductors, and Josephson junctions. The common theme remains the same: the interplay between pumping and nonlinear elements.

The strong pumping of a material with a complex crystal structure will induce coherent oscillations in a few primary modes. These oscillations trickle down to several other modes as a result of nonlinear couplings. Thus, every symmetry-allowed mode will be parametrically driven to a degree with strong enough pumping. With this understanding, the model proposed here is only one out of the numerous possible other routes to achieve parametric driving. For instance, the symmetry-allowed quartic coupling in fulleride superconductors $\sim\varphi_{\mathrm{IR}}^{2}\varphi^{2}$ directly translates the coherent motion of $\varphi_{\mathrm{IR}}$ to a parametric drive for $\varphi$. Here, $\varphi_{\mathrm{IR}}$ is one of the IR-active modes of $C_{60}$ such as $\mathrm{T}_{1u}(1-4)$, and $\varphi$ is a Raman-active mode such as $\mathrm{Hg}(7-8)$ that couples strongly to conduction electrons Mitrano et al. (2016); Varma et al. (1991). Even though achieving parametric driving is material-specific, it leads to the same qualitative physics. The present paper mainly deals with the universal consequences of parametric driving.

III Parametric amplification of phonon-mediated electron-electron attraction: a first look

Our goal in this section is to demonstrate the resonant amplification of the electron-phonon coupling in the presence of the drive. For the time being, we neglect the complex epiphenomena such as the nonequilibrium evolution of electrons, phonon dissipation and retardation, and the feedback between electrons and phonons. Instead, we resort to a perturbative treatment and elementary methods in order to elucidate the main ideas. We will revisit the problem again in later section and provide a comprehensive account using the nonequilibrium Migdal-Eliashberg theory. The latter treatment is naturally more cumbersome than the physical account given in this section. The present analysis serves a guideline to identify and interpret the results of the upcoming detailed analysis.

As a first step, we assume that the lattice nonlinearity $\mathcal{V}(\varphi)$ and the electron-phonon coupling $g_{\mathbf{k},\mathbf{k}^{\prime}}$ are both weak compared to the drive, such that we can study the coherent motion of the lattice in isolation. The classical equation of motion for $\langle\hat{\varphi}_{0}(t)\rangle$ (the $\mathbf{q}=0$ mode) is easily found as:

[TABLE]

The normalization constant $\sqrt{N}$ results from the definition of the Fourier operators given earlier, i.e. $\langle\hat{\varphi}_{0}\rangle=\sqrt{N}\,\langle\hat{\varphi}_{j}\rangle$ where $\langle\hat{\varphi}\rangle$ is the coordinate of an arbitrary single ion $j$. We assume that the temporal variation scale of $F_{\mathrm{env}}(t)$ is much longer than the drive period. Thus, for an adiabatically ramped up $F_{\mathrm{env}}(t)$, we find:

[TABLE]

Near the resonance $\Omega_{\mathrm{drv}}=\omega_{0}/2$, the oscillatory term dominates the dc term in amplitude. The precise value of the prefactors of the dc and ac terms are not important for the present discussion and in a more realistic setting, both get corrections from phonon damping, nonlinearities, etc. Quite generally though, we have $\langle\hat{\varphi}_{0}(t)\rangle\approx\sqrt{N}\,\varphi_{0}(t)+\sqrt{N}\,\varphi_{1}(t)\cos(2\Omega_{\mathrm{drv}}t)$ where $\varphi_{0}(t)$ and $\varphi_{1}(t)$ are slowly-varying functions of time. With this understanding, we drop the time labels from $\varphi_{0}$ and $\varphi_{1}$ hereafter and treat than as given quasi-steady parameters.

The local lattice nonlinearity terms couple the coherent uniform motion of the lattice to $\pm\mathbf{q}$ modes. For instance, the leading order correction resulting from the cubic nonlinearity $\sim\varphi^{3}$ is found by replacing one of the operators with $\langle\hat{\varphi}_{0}(t)\rangle$. Momentum conservation implies opposite momenta for the remaining two operators:

[TABLE]

Likewise, the leading order contribution from the quartic nonlinearity is found by replacing two of the operators with $\mathbf{q}=0$, which yields:

[TABLE]

The dc terms result in the renormalization of the phonon frequency, e.g. $\omega_{\mathbf{q}}^{2}\rightarrow\omega_{\mathbf{q}}^{2}-\kappa_{3}\varphi_{0}/2$ for the cubic nonlinearity, and $\omega_{\mathbf{q}}^{2}\rightarrow\omega_{\mathbf{q}}^{2}-\kappa_{4}(\varphi_{0}^{2}/4+\varphi_{1}/8)$. Such corrections are precisely the time-averaged renormalized lattice properties that we discussed earlier in the introduction and can enhance or suppress the effective electron-phonon coupling on their own account.

As we will show soon, the most intriguing effect is purely dynamical and stems from the ac term. We neglect the dc corrections for simplicity hereafter. At the present order in perturbation theory, only $\pm\mathbf{q}$ opposite momentum pairs couple and thus we may focus on a single momentum pair without the loss of generality. We also only consider the cubic nonlinearity. It will soon become apparent that both types of nonlinearity give rise to the same resonant amplification phenomenon. The Hamiltonian is given as $\hat{H}_{\pm\mathbf{q}}(t)=\hat{H}_{\mathrm{e}}+\hat{H}_{\mathrm{p,\pm\mathbf{q}}}(t)+\hat{H}_{\mathrm{ep,\pm\mathbf{q}}}$, where $\hat{H}_{e}=\sum_{\mathbf{k},\sigma}\xi_{\mathbf{k}}\,c^{\dagger}_{\mathbf{k},\sigma}c^{\phantom{\dagger}}_{\mathbf{k},\sigma}$ and:

[TABLE]

where $\hat{\pi}_{\mathbf{q}}=(b^{\phantom{\dagger}}_{-\mathbf{q}}-b^{\dagger}_{\mathbf{q}})/(2i)$ is the conjugate momentum to $\hat{\varphi}_{\mathbf{q}}$ 111Our definition of $\hat{\phi}_{\mathbf{q}}$ and $\hat{\pi}_{\mathbf{q}}$ deviates from the standard literature, resulting in unbalanced prefactors of the kinetic and potential terms in the harmonic oscillator Hamiltonian., and $g_{\mathbf{q}}$ is the linear electron-phonon coupling constant (we have assumed $g_{\mathbf{k},\mathbf{k}^{\prime}}\approx g_{\mathbf{k}-\mathbf{k}^{\prime}}$).

Finally, $\rho_{\mathbf{q}}=\sum_{\mathbf{k},\sigma}c^{\dagger}_{\mathbf{k}+\mathbf{q},\sigma}c^{\phantom{\dagger}}_{\mathbf{k},\sigma}$ is the electron charge density operator. The analysis can be further simplified by performing a canonical change of variables to standing wave phonon operators:

[TABLE]

and their corresponding conjugate momenta:

[TABLE]

Here, $M$ is the ion mass. It is readily verified that $[\hat{Q}_{+},\hat{P}_{+}]=[\hat{Q}_{-},\hat{P}_{-}]=i\hbar$ while all other commutators vanish. The Hamiltonian can be easily written in terms of $\hat{Q}_{\pm}$ and $\hat{P}_{\pm}$ operators:

[TABLE]

where $\tilde{g}_{\mathbf{q}}=g_{\mathbf{q}}\sqrt{2M\omega_{\mathbf{q}}/\hbar}$, $\hat{\rho}_{+}=(\hat{\rho}_{\mathbf{q}}+\hat{\rho}_{-\mathbf{q}})/\sqrt{2}$, $\hat{\rho}_{-}=i(\hat{\rho}_{\mathbf{q}}-\hat{\rho}_{-\mathbf{q}})/\sqrt{2}$, and:

[TABLE]

The Hamiltonian is the sum of two decoupled harmonic oscillators ($\alpha=\pm$) each with linear coupling to a standing electronic charge density wave. Since the $\alpha=\pm$ modes undergo a similar evolution, we will drop the $\pm$ index hereafter and focus on a one of the modes. In the absence of electron-phonon coupling, the problem reduces to a parametrically driven harmonic oscillator. Corrections arising from electron-phonon coupling can be be studied order by order using time-dependent perturbation theory in the weak-coupling limit $g_{\mathbf{q}}\nu(0)\ll 1$. This can be done, e.g. via a perturbative expansion of the unitary evolution operator $\hat{U}(t)=\mathcal{T}\exp\left[-i\hbar^{-1}\int_{0}^{t}\mathrm{d}t^{\prime}\,\hat{H}_{\pm\mathbf{q}}(t^{\prime})\right]$ in the powers of $\tilde{g}_{\mathbf{q}}$ Fujita et al. (2010). The leading order correction to the action is $\mathcal{O}(\tilde{g}_{\mathbf{q}}^{2})$ and is easily found as:

[TABLE]

where $\mathcal{D}^{R}_{QQ}(t,t^{\prime})=-i\theta(t-t^{\prime})\langle[\hat{Q}(t),\hat{Q}(t^{\prime})]\rangle$ is the retarded phonon correlator and $\epsilon$ is an infinitesimal. It is well-known that this correction implies an attractive interaction between the electrons in the long wave-length regime $|\xi_{\mathbf{k}\pm\mathbf{q}}-\xi_{\mathbf{k}}|\ll\omega_{\mathbf{q}}$. In this regime of interest, the phase winding of the electron charge density excitation is much slower than the phonon time scale. Thus, to simplify the discussion further, we simply neglect the relative phase winding of the electron charge density waves and set $\hat{\rho}(t^{\prime})\rightarrow\hat{\rho}(t)$ from the outset. Note that this coincides with the usual $\varepsilon_{\mathbf{k}}=\varepsilon_{\mathbf{k}^{\prime}}=\varepsilon_{F}$ approximation in the BCS treatment. This results in the following simple expression for the phonon-mediated attractive interaction, $\Delta S^{(2)}(t)=U(t)\,\hat{\rho}(t)\,\hat{\rho}(t)$ where:

[TABLE]

Since the Hamiltonian is time-dependent, $U(t)$ is expected to be time-dependent as well. In particular, for the periodic Hamiltonian given in Eq. (11a), $U(t)$ further admits a Fourier expansion $U(t)=\sum_{n=-\infty}^{+\infty}U_{n}\,e^{-2in\Omega_{\mathrm{drv}}t}$ in the harmonics of $2\Omega_{\mathrm{drv}}$. In plain words, $U(t)$ is proportional to the mean displacement of the oscillator up to time $t$ in response to a momentum boost at all prior times.

The retarded phonon propagator $\mathcal{D}_{QQ}^{R}(t,t^{\prime})$ is most easily calculated using Heisenberg equations for $\hat{Q}(t)$ and $\hat{P}(t)$:

[TABLE]

The Heisenberg equations exactly coincide with the classical equations of motion due to the Ehrenfest’s theorem. The formal solution of these equations can be expressed in terms of Mathieu functions. For given Heisenberg operators at time $t^{\prime}$, we find:

[TABLE]

The explicit expressions for the $\mathfrak{M}$-functions are given in Eq. (A) in terms of even and odd Mathieu functions. In the limit $\alpha\rightarrow 0$ (no drive), we have $\lim_{\alpha\rightarrow 0}\mathfrak{M}_{QQ}(t,t^{\prime})=\lim_{\alpha\rightarrow 0}\mathfrak{M}_{PP}(t,t^{\prime})=\cos[\omega_{\mathbf{q}}(t-t^{\prime})]$, and $(\omega_{\mathbf{q}}/\Omega_{\mathrm{drv}})\,\lim_{\alpha\rightarrow 0}\mathfrak{M}_{QP}(t,t^{\prime})=-(\Omega_{\mathrm{drv}}/\omega_{\mathbf{q}})\,\lim_{\alpha\rightarrow 0}\mathfrak{M}_{PQ}(t,t^{\prime})=-\sin[\omega_{\mathbf{q}}(t-t^{\prime})]$, thus reducing Eq. (III) the usual non-driven harmonic oscillator evolution. Furthermore, $\mathfrak{M}_{QQ}(t^{\prime},t^{\prime})=\mathfrak{M}_{PP}(t^{\prime},t^{\prime})=1$ and $\mathfrak{M}_{PQ}(t^{\prime},t^{\prime})=\mathfrak{M}_{QP}(t^{\prime},t^{\prime})=0$, as it is also required from the initial condition. The retarded correlator is immediately calculated using Eq. (III), giving an exceedingly simple result:

[TABLE]

This result presents several important features:

• •

The retarded phonon propagator $\mathcal{D}^{R}_{QQ}(t,t^{\prime})$, and consequently $U(t)$ (see Eq. 14), are fully determined by $\mathfrak{M}_{QP}(t,t^{\prime})$. The latter is in turn fully determined by the Heisenberg equations, and is independent of the initial wavefunction of phonons. In other words, for the Hamiltonian given in Eq. (11a), one obtains the same effective interaction $U(t)$ in all equilibrium and nonequilibrium phonon states. This is a direct consequence of the linearity of the harmonic oscillator evolution. As a corollary, this result immediately shows the peculiar cancellation of Bose enhancement factors between phonon absorption and emission processes in the conventional textbook diagrammatic calculation of $U$ in equilibrium Schrieffer (1983).

• •

$\mathcal{D}_{QQ}^{R}(t,t^{\prime})\propto\mathfrak{M}_{QP}(t,t^{\prime})$ admits a simple classical interpretation: it coincides with the displacement of a classical driven oscillator at time $t$ in response to a momentum jump at time $t^{\prime}$, i.e. $\mathcal{D}_{QQ}^{R}(t,t^{\prime})~{}=~{}\hbar\,\delta Q_{\mathrm{cl}}(t)/\delta P_{\mathrm{cl}}(t^{\prime})$. This result can be obtained independently and more directly using quantum phase-space methods Polkovnikov (2010).

• •

The origin of phonon-mediated attraction is purely quantum mechanical. Had $\hat{Q}$ been a classical operator, it would commute with itself at different times and we would obtain $\lim_{\hbar\rightarrow 0}\mathcal{D}^{R}_{QQ}(t,t^{\prime})=0$. Nonetheless, this analysis provides a classical analogy as a response to a momentum jump (see the previous remark).

The second remark implies that classical trajectories of a parametrically driven oscillator following a momentum jump encode the necessary and sufficient information to calculate $U(t)$. As a first example, let us consider the non-driven limit. In this case, we find $\lim_{\alpha\rightarrow 0}\mathfrak{M}_{QP}(t,t^{\prime})=-(\Omega_{\mathrm{drv}}/\omega_{\mathbf{q}})\,\sin\omega_{\mathbf{q}}(t-t^{\prime})$

[TABLE]

The above result indeed corresponds the $QP-$response of the classical harmonic oscillator up to a factor of $\hbar$. Plugging this result in Eq. (14), we find $U=-|\tilde{g}_{\mathbf{q}}|^{2}/(M\omega_{\mathbf{q}}^{2})=-2|g_{\mathbf{q}}|^{2}/(\hbar\omega_{\mathbf{q}})$ which is the usual time-independent equilibrium result Mahan (2013). Note that regularizing prefactor $e^{-\epsilon(t-t^{\prime})}$ is crucial for obtaining this result, without which the $t^{\prime}$-integration is ill-defined. In a more realistic model with finite phonon damping, regularization is unnecessary.

Fig. 3 shows the classical trajectories in response to a sudden momentum jump at $t=0$ in the presence of a finite parametric drive. The red circular orbit shows the response in the absence of the drive. It is noticed that the classical trajectory diverge on resonance $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}/2$ (middle panel), implying an infinitely enhanced response. Intuitively, the lattice becomes critically unstable in the presence of a resonant drive of finite amplitude such that smallest perturbation causes an infinitely large deformation. Just below the resonance (left panel), the trajectories are noticeably expanded in the phase space though remain bounded. This corresponds to a moderately enhanced $U(t)$. Finally, the response is attenuated just above the resonance (right panel), which is the expected asymmetric behavior near a parametric resonance.

The time-dependent effective attraction $U(t)$ can be obtained numerically using Eqs. (III) and (14). The results are shown in Fig. 4. According to the previously mentioned periodicity $U(t+\pi/\Omega_{\mathrm{drv}})=U(t)$, we have only shown one full period in $t$. The immediately noticeable feature is the large temporal variations of $U(t)$ near parametric resonances, which is a manifestation of the wild variations of classical trajectories in the phase-space in this regime. The perturbation series for $U(t)$ can be worked out using Eqs. (14), (III), and the expressions given in Appendix A. We quote the final result here:

[TABLE]

Note the parametric resonances at $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}/2$ in the first order term and $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}/4$ in the second order term, as well as the appearance of higher order harmonics of $2\Omega_{\mathrm{drv}}$.

We remark that the parametric resonances of the Mathieu oscillator, $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}/n$, $n\in\mathbb{N}$, do not necessarily imply a corresponding resonance in the effective attraction $U(t)$. As remarked after Eq. (14), $U(t)$ can be interpreted as the mean displacement of the oscillator in response to a momentum boost. The mean displacement can behave properly even for divergent trajectories. For example, the leading $n=1$ parametric resonance at $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}$ leaves $U(t)$ regular while leading to a divergent response at the same time. This can be easily noticed in the perturbation analysis: Eq. (A) shows a resonance at $\Omega_{\mathrm{drv}}=\omega_{\mathbf{q}}$ while $U(t)$ remains regular, see Eq. (19) and Fig. 5.

The time average and variance of $U(t)/U_{\mathrm{eq}}$ can be found readily using the above result:

[TABLE]

The leading order correction to $\langle U(t)\rangle$ is $\mathcal{O}(\alpha^{2})$. Fig. 5 shows the average and variance of $U(t)$ within a period. The dashed lines show the above perturbative results. As expected from the study of classical trajectories (see Fig. 3), $\langle U(t)\rangle$ is enhanced for frequencies below the resonance and is suppressed above the resonance. Perhaps more importantly, the temporal variance of $U(t)$ is significantly increased on either side of the resonance. This can be seen e.g. from Eq. (20b). In the next section, we discuss the important role of enhanced temporal variations of $U(t)$ in enhancing the superconducting transition temperature.

III.1 The superconducting transition temperature for a time-dependent effective interaction

We showed that the parametric drive of the lattice results in the enhancement of the effective attraction mediated between the electrons, $U(t)$. We further showed that $U(t)$ exhibits large oscillations in $t$ near parametric resonances (see Fig. 4), such that $U(t)$ takes on values that are significantly higher and lower than its equilibrium value. According to the BCS theory, $T_{c}[U]\sim\exp[-1/(\nu(0)U)]$ where $\nu(0)$ is the electronic density of states (EDOS) at the Fermi surface. It is tempting to naively propose a replacement $U\rightarrow U(t)$ in the BCS formula and propose $T_{c}^{\mathrm{drv}}\sim\langle T_{c}[U(t)]\rangle$. Since $T_{c}[U]$ is a convex function of $U$ in the weak-coupling limit, one would then conclude $\langle T_{c}[U(t)]\rangle>T_{c}[\langle U(t)\rangle]$, i.e. temporal variations of $U(t)$ can increase $T_{c}$ even if $\langle U(t)\rangle$ remains constant or even decrease. Of course, this argument lacks rigor and $T_{c}^{\mathrm{drv}}$ must be found within a proper Floquet extension of the BCS theory Knap et al. (2015). To this end, we assume:

[TABLE]

where $U_{0}$ and $U_{1}$, can be read from Eq. (19):

[TABLE]

The superconducting gap inherits the periodicity of $U(t)$ such that $\Delta(t)=\frac{1}{N}\sum_{\mathbf{k}}\langle\psi_{\mathbf{k},\uparrow}(t)\,\psi_{-\mathbf{k},\downarrow}(t)\rangle\rightarrow\sum_{n=-\infty}^{+\infty}\Delta_{n}\,e^{2ni\Omega_{\mathrm{drv}}t}$. At the onset of pairing, the Floquet BCS gap equation takes the following form Knap et al. (2015):

[TABLE]

where:

[TABLE]

Here, $T$ is the temperature and $\omega_{c}$ is the UV cutoff for $U$. For an Einstein phonon with frequency $\omega_{0}$, we expect $\omega_{c}\sim\omega_{0}$. In principle, $T_{c}$ must be found such that Eq. (23) admits a non-trivial solution for $\{\Delta_{n}\}$ and since it is a homogeneous equation in $\{\Delta_{n}\}$, it reduces to the vanishing determinant condition for an (infinitely large) matrix.

A closed form solution for $T_{c}$ seems to be out of reach in general and one must resort to numerical methods. We attempt to find an approximate analytic solution with additional assumptions $U_{0}\nu(0)\ll 1$ and $U_{1}\ll U_{0}$. Strictly speaking, the last assumption does not generally apply to our problem since $|U_{1}|\geq|U_{0}|$ near the resonances. Our main goal here is to demonstrate how temporal variations in $U(t)$ increases $T_{c}$ rather than presenting a rigorous analysis; the latter is the objective of the upcoming sections. Thus, the additional simplifying assumptions must be taken with this understanding in mind.

As a first step, we observe that $\Delta_{n+1}/\Delta_{n}\sim\mathcal{O}(U_{1}/U_{0})\ll 1$. Therefore, we may neglect $\Delta_{n}$ for $|n|\geq 2$ for small $U_{1}$. This reduces the infinite set of equations for $\{\Delta_{n}\}$ to only three equations for $\Delta_{0}$, $\Delta_{1}$, and $\Delta_{-1}$. Omitting $\Delta_{\pm 1}$ between the equations and assuming $\Delta_{0}\neq 0$, we find the following approximate pairing condition:

[TABLE]

Approximate expressions for $F_{n}$ can be found in the limit $T/\omega_{c},T/\Omega_{\mathrm{drv}}\ll 1$ using the Sommerfeld expansion technique. We quote the final result here:

[TABLE]

where $\gamma\simeq 0.577$ is the Euler-Mascheroni constant. The final result resembles the BCS formula for $T_{c}$, though, with $U_{0}$ replaced with an effective interaction $U_{\mathrm{eff}}$:

[TABLE]

We notice that $U_{\mathrm{eff}}>|U_{0}|$ in both cases. The last result provides a well-founded justification for our preliminary heuristic argument based on the convexity of $T_{c}[U]$.

Finally, we have plotted $T_{c}^{\mathrm{drv}}$ based on the approximate analytic results obtained in this section (Eqs. III.1, 27a, 27b) in Fig. 6. The left heat map plot shows the full result, when both dc and ac components of $U(t)$ are taken into account. It is noticed that $T_{c}$ is dramatically enhanced both below and above the resonance. The right plot shows the result when only the dc component, $U_{0}$, is kept. While $T_{c}^{\mathrm{drv}}$ is enhanced below the resonance, it is suppressed for $\Omega_{\mathrm{drv}}>\omega_{0}/2$. This result can be understood by appealing to Fig. 5. The time-averaged interaction $\langle U(t)\rangle=U_{0}$ is lower than $U_{\mathrm{eq}}$ for $\Omega_{\mathrm{drv}}>\omega_{0}/2$, and $T_{c}$ is decreased accordingly. In contrast, including the ac component brings in large oscillations $\propto U_{1}$ which offset the loss in $U_{0}$ by allowing $U(t)$ to exceed $U_{\mathrm{eq}}$ during a fraction of each cycle.

We will obtain a plot similar to Fig. 6 later using a Floquet extension of the Migdal-Eliashberg theory (see Fig. 10), and we will find that the approximate analytic picture provided here remains remarkably accurate. We conclude this section by noting that had we included higher-order harmonic corrections in $U(t)$, we would get additional parametric resonances beside to the main one. Those would subsequently lead to additional “tongues” in Fig. 6 around $\Omega_{\mathrm{drv}}=\omega_{0}/4,\omega_{0}/8,\ldots$ (see Fig. 10).

III.2 Higher-order nonlinearities, phonon damping, and parametric resonance

Our discussion in the last two sections was based on a cubic lattice nonlinearity. We showed that the nonlinearity leads to a constant renormalization of $\omega_{\mathbf{q}}$ and produces a $\sim\cos(2\Omega_{\mathrm{drv}}t)$ periodic correction to $\omega_{\mathbf{q}}$ at the leading order (see Eq. 6). The scenario remains similar for quartic nonlinearities, save for an additional $\sim\cos(4\Omega_{\mathrm{drv}}t)$ correction to $\omega_{\mathbf{q}}$. This can be seen by expanding the square brackets in Eq. (7). In fact, for a nonlinearity $\sim\hat{\varphi}^{n}$ ($n\geq 3$), we find $\hat{\varphi}^{n}\approx\langle\hat{\varphi}_{0}(t)\rangle^{n-2}\,\sum_{\mathbf{q}}\hat{\varphi}_{\mathbf{q}}\,\hat{\varphi}_{-\mathbf{q}}$. Subsequently, $\langle\hat{\varphi}_{0}(t)\rangle^{n-2}$ can be expanded in harmonics of $2\Omega_{\mathrm{drv}}$ such that \omega_{\mathbf{q}}^{2}\rightarrow\omega_{\mathbf{q}}^{2}\big{[}1+\sum_{n=0}^{n-2}\alpha_{n}\,\cos(2n\Omega_{\mathrm{drv}}t)\big{]} for some $\{\alpha_{n}\}$. This analysis can be naturally extended for any smooth nonlinear potential. Since all parametric drive terms are harmonics of $2\Omega_{\mathrm{drv}}$, the physics is expected to remain qualitatively similarly to the single harmonic case. For example, the case of double parametric drive terms is studied in Ref. Sofroniou and Bishop (2014), where it is shown if the drive frequencies are integer multiples of one another, the same instability “tongue” patterns are obtained. Finally, let us mention in passing that with a finite phonon damping and for a fixed driving strength, the infinite cascade of ideal parametric resonances at $\Omega_{\mathrm{drv}}=\omega_{0}/(2n)$ will be truncated above a certain order Gunderson et al. (1974). For instance, only two resonances are noticeable in Fig. 4 (black lines). The reason can be traced back to using a small phonon damping in the numerics.

III.3 Intermission

We derived a simple formula for the phonon-mediated electron-electron attraction $U(t)$ in the presence of parametric drive (Eq. 14). We showed that $U(t)$ is independent of the initial phonon wavefunction and is fully determined by the Hamiltonian and the drive. We related $U(t)$ to the momentum-jump response of a classical parametric oscillator and used this classical analogy to demonstrate the parametric amplification of $U(t)$ by appealing to the classical phase-space trajectories of the parametric oscillator. Finally, we calculated an analytical formula for $T_{c}$ using a Floquet generalization of the BCS theory Knap et al. (2015) and demonstrated that the oscillations $U(t)$, as well as its increased time-average, lead to enhanced $T_{c}$.

Parametric amplification, while enhancing $T_{c}$ at the first glance, leads to several undesirable consequences as well. It leads to parametric phonon generation, a well-known phenomenon in the context of early universe field theory Berges and Serreau (2003). The generated phonons heat up the electrons and decrease their coherence. The oscillatory effective electron-electron interaction $U(t)$ also generates electron-hole excitations on its own account Knap et al. (2015). An unbiased and consistent analysis of these competing effects is a challenging task and requires a rigorous and unified treatment. In the remainder of the paper, we will develop such a formalism and revisit the problem one more time in full detail.

IV The Floquet-Migdal-Eliashberg quantum kinetic theory

The major steps of the forthcoming calculation is summarized graphically in Fig. 7. Initially, the electron-phonon system is prepared in a normal-conducting equilibrium state with temperature $T_{i}>T_{c}$, where $T_{c}$ is the equilibrium critical temperature for Cooper pairing. The external drive is then ramped up and initiates the nonequilibrium quantum dynamics.

We study the coupled quantum dynamics of the lattice displacement, phonons, and electrons by deriving and numerically solving a set of quantum kinetic equations specifically tailored for investigating driven electron-phonon systems. We start this theoretical development from a two-particle irreducible effective active (2PI-EA) formulation of the Migdal-Eliashberg theory extended to nonlinear lattices in Sec. IV.1, followed by a Floquet generalization of the ensuing quantum kinetic equations in Secs. IV.3-IV.5.

Once the evolution of the system in the normal-conducting state is calculated, it is checked for the emergence of a Floquet superconducting instability during the evolution. This step in enabled by deriving a Floquet-Migdal-Eliashberg gap equation for the onset of pairing in Sec. V. A summary of the numerical methods is provided in Appendix F.

IV.1 Two-particle irreducible effective action (2PI-EA) and real-time evolution equations

We study the nonequilibrium dynamics of the driven electron-phonon system using the functional technique of 2PI effective actions (2PI-EA). A prominent feature of this approach is the guaranteed satisfaction of conservation laws and the absence of secular terms that arise in non-self-consistent perturbation theory Berges (2004). Both of these features are necessary for stable and physically meaningful description of nonequilibrium quantum dynamics. In the context of our problem, the effective action is a functional of the uniform coherent lattice displacement $\varphi$, the Nambu closed-time-path (CTP) electron propagator $\hat{\mathcal{G}}$, and the CTP phonon propagator $\mathcal{D}$. The lattice displacement field (a “classical” object) is coupled to the phonon propagator (a “quantum” object) via lattice nonlinearities and lead to the previously discussed parametric amplification effect.

We will work in the Migdal’s limit $\omega_{\mathrm{ph}}/W_{\mathrm{el}}\ll 1$ ($\omega_{\mathrm{ph}}$ is the typical optical phonon frequency and $W_{\mathrm{el}}$ is the conduction electron bandwidth) where electron-phonon vertex corrections are suppressed by the powers of $\omega_{\mathrm{ph}}/W_{\mathrm{el}}$ and may be controllably neglected Migdal (1958); Engelsberg and Schrieffer (1963). In the ideal Migdal limit $\omega_{\mathrm{ph}}/W_{\mathrm{el}}\rightarrow 0$, the 2PI-EA, $\mathsf{\Gamma}[\varphi,\mathcal{D},\mathcal{D}]$ truncated to the two-loop order in lattice nonlinearities is given as:

[TABLE]

The definition and symmetries of CTP propagators as well as the definition of various real-time components are given in Appendix B. The bare electron and phonon propagators are given as:

[TABLE]

and the classical action of the coherent lattice displacement is given as:

[TABLE]

We have introduced a small phenomenological Ohmic friction $\gamma_{0}\omega_{0}\partial_{t}$ for $\varphi$ to model the effect of classical dissipation due to coupling to other lattice modes. In principle, the friction must be accompanied by a stochastic Langevin noise for consistency. The presence of a strong driving term, however, dominates over the Langevin noise in practice and allows us to neglect the latter. Finally, the electron-hole bubble $\Pi_{\mathbf{q}}$ is defined as:

[TABLE]

and the local phonon propagator $\mathcal{D}_{\ell}$ appearing in the last three terms of $\eqref{eq:2PI}$ is defined as:

[TABLE]

We note that the external drive only appears in the classical action $S_{\mathrm{cl}}[\varphi]$ and couples to $\varphi_{\mathbf{q}=0}$. The latter directly couples to finite-$\mathbf{q}$ phonons via nonlinearities (see the fifth and sixth diagrams in Eq. 28), and to electrons if ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}}\neq 0$ (SSH-type phonons, see the forth diagram in Eq. 28).

The equations of motion are found by imposing the stationarity condition on $\mathsf{\Gamma}$ with respect to $\varphi$, $\mathcal{G}$, and $\mathcal{D}$. For the classical displacement, we find:

[TABLE]

where:

[TABLE]

Here, $\chi(t)$ is the phonon tadpole and $\eta(t)$ is the electron-mediated classical force on the uniform lattice displacement. The evolution of $\mathcal{D}$ is given by the following set of Kadanoff-Baym (KB) equations:

[TABLE]

where:

[TABLE]

is a local time-dependent potential acting on phonons. Aside from $\chi(t)$ which is a phonon self-interaction, the second and third terms are responsible for the parametric drive of phonons in connection to the analysis given in Sec. III.

Finally, the evolution of the CTP Nambu electron propagator $\mathcal{G}$ is given by the following set of KB equations:

[TABLE]

where the Nambu self-energy is given as:

[TABLE]

which is the usual Migdal-Eliashberg self-energy in real time. Explicit equations for the retarded (R), advanced (A), and Keldysh (K) components of $\hat{\mathcal{G}}$ and $\mathcal{D}$ can be worked out from Eqs. (33)-(33) and Eqs. (31)-(31) using Langreth rules Rammer (2007), respectively.

In principle, the solution of the coupled integro-differential equations derived in this section, while being a daunting task, yields a complete and unbiased analysis. Given that our goal in the present paper is to give a transparent account of the key mechanisms that play a role in enhancing or suppressing superconductivity, we find it rather beneficial to simplify the model to the greatest possible extent without sacrificing any qualitative physics.

IV.2 The trimmed-down model

In this section, we present and discuss several simplifying approximations which we adopt in the rest of the paper. These assumptions are presented in a single section, rather than incrementally, for better clarity.

Holstein-type electron-phonon coupling— Depending on the nature of the electron-phonon coupling, ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ may assume different dependencies on $\mathbf{k}$ and $\mathbf{k}^{\prime}$. For Holstein-type phonons, appropriate for describing longitudinal optical (LO) phonons, ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ only depends on the net momentum of the electron-hole excitation, i.e. $g_{\mathbf{k},\mathbf{k}^{\prime}}\sim g_{\mathbf{k}-\mathbf{k}^{\prime}}$. On the other hand, for Su-Shrieffer-Hieger (SSH) type phonons, ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ will depend on the individual momenta. We restrict our analysis to Holstein-type coupling here.

For a realistic description of an electron-phonon system, some aspects of the Coulomb interaction must be incorporated into the model, in particular, the renormalized electron dispersion $\xi_{\mathbf{k}}$ and screening of the electron-phonon coupling $g^{\mathrm{(scr)}}_{\mathbf{\mathbf{k}}-\mathbf{\mathbf{k}^{\prime}}}\approx{g}^{\phantom{\dagger}}_{\mathbf{k}-\mathbf{k}^{\prime}}/\varepsilon(\Omega,\mathbf{k}-\mathbf{k}^{\prime})$ Mahan (2013). Here, $\varepsilon$ is the dielectric function and $\Omega$ is the relevant energy scale of the dynamical screening process which is commonly set to zero. In our problem, the largest relevant frequency for dynamical screening is set by external drive and the optical phonon peak frequency, whichever is the largest. In the regime relevant to the experiments, both are an order of magnitude smaller than the typical plasma frequency $\omega_{p}$. For example, $\omega_{p}\approx 0.5\sim 2~{}{\rm eV}$ for fulleride superconductors whereas the typical frequency of optical phonons is $\omega_{\mathrm{ph}}\approx 10\sim 100~{}{\rm meV}$ Gunnarsson and Fullerides (2004). Therefore, we may safely use the static Thomas-Fermi (TF) dielectric function:

[TABLE]

The screened coupling in the long wavelength $|\mathbf{k}-\mathbf{k}^{\prime}|\ll q_{\mathrm{TF}}$ is almost perfectly screened. As a consequence, (1) the direct coupling of the uniform $\mathbf{q}=0$ lattice displacement to conduction electrons (described by the forth diagram in Eq. (28) and Eq. (30b)) is vanishingly small and can be neglected; (2) the time-dependent correction $g_{\mathbf{k},\mathbf{k}}\,\varphi(t)$ to the electron dispersion in Eqs. (33)-(33) vanishes as well.

Non-dispersive (Einstein) optical phonons— So far, we have assumed a general dispersion $\omega_{\mathbf{q}}$ for the optical phonon. We neglect the phonon dispersion hereafter, i.e. $\omega_{\mathbf{q}}\rightarrow\omega_{0}$, which an excellent approximation for a large class of materials, including fulleride superconductors Gunnarsson (1997).

Local approximation for $\Pi_{\mathbf{q}}$— In the absence of $\Pi_{\mathbf{q}}$ phonon self-energy correction appearing on the right hand side of Eqs. (31)-(31) and using Einstein phonons, the phonon propagator $\mathcal{D}_{\mathbf{q}}$ will have no $\mathbf{q}$-dependence. Therefore, the $\mathbf{q}$-dependence of $\mathcal{D}_{\mathbf{q}}$ is entirely induced by $\Pi_{\mathbf{q}}$. The $\mathbf{q}$-dependence of the latter is inherited from ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ and the electron dispersion $\xi_{\mathbf{k}}$ (See Eq. IV.1) and is non-universal. In an attempt to simplify the model, we propose a local approximation for $\Pi_{\mathbf{q}}$ in the spirit of dynamical mean field theory (DMFT):

[TABLE]

The local approximation for $\Pi_{\mathbf{q}}$ has an important practical advantage. It removes the $\mathbf{q}$-dependence from $\mathcal{D}_{\mathbf{q}}$ and allows us to study a single momentum mode.

Fermi-surface averaged (FSA) electron self-energy— The Cooper pairing process in the majority of conventional superconductors only involves electrons within in a thin shell about the Fermi surface. This observation is the basis of a widely used approximation where the Migdal-Eliashberg electron self-energy $\hat{\Sigma}_{\mathbf{k}}$ is replaced with its Fermi-surface averaged (FSA) approximation:

[TABLE]

This approximation, being akin to DMFT-type approximation, fully retains the dynamical structure of the self-energy while simplifying the momentum summations by removing the spatial structure of the self-energy. We remark that the FSA approximation is indeed an excellent approximation while adopting the previous two approximations: for a local Holstein-type electron-phonon coupling and Einstein phonons, $\hat{\Sigma}_{\mathbf{k}}$ naturally loses is $\mathbf{k}$-dependence (see Eq. 34).

Flat electronic density of states (EDOS)— We neglect the variations of the EDOS $\nu(\xi)$ and pin it to its value at the Fermi surface $\nu(0)$. This is an excellent approximation in three-dimensional systems.

Ideal Migdal’s limit— We work in the ideal Migdal’s limit $\omega_{\mathrm{ph}}/W_{\mathrm{el}}\rightarrow 0$.

IV.3 Floquet-Boltzmann quantum kinetic formalism

Despite the simplifications proposed to the model in Sec. IV.2, the solution of the KB equations and obtaining the two-time propagators remains a challenging numerical task. Provided that the external perturbing field varies on time and length scales longer than the intrinsic microscopic time and scale of the many-body system, the KB equations can be controllably reduced to one-time kinetic-type equation using the machinery of Wigner transforms and gradient expansion Kadanoff and Baym (1962); Rammer (2007). The case of driven systems is subtler though since the external drive $F(t)$ has both fast and slow components: even though the temporal variations of the amplitude $F_{\mathrm{env}}(t)$ occurs on long time scales (compared to the phonon period), the multiplicative oscillatory factor $\cos(\Omega_{\mathrm{drv}}t)$ varies on the same scale as the phonon frequency in the interesting nearly-resonant regime. In this section, we show that by introducing Floquet bands via short-time Fourier transforms, we can performs gradient expansion on the amplitude of Fourier components and derive Boltzmann-type quantum kinetic equations. This procedure is schematically shown in Fig. 8.

Floquet-Wigner (FW) transform— Let $\mathcal{A}(t_{1},t_{2})$ be an arbitrary two-time function, such as the $\hat{\mathcal{G}}(t_{1},t_{2})$ and $\mathcal{D}(t_{1},t_{2})$. The Wigner transform $\mathcal{A}(\omega,t)$ is formally defined as:

[TABLE]

Here, $t=(t_{1}+t_{2})/2$ is the center-of-mass (COM) time and $\tau=t_{1}-t_{2}$ is the relative time on which we perform a Fourier transform. The COM-time dependence vanishes at equilibrium. By continuity, we expect the presence of a slowly-varying external field to induce a similarly slowly-varying COM-dependence on $\mathcal{A}(\omega,t)$. This permits a controllable series expansion in successive COM-time derivatives of $\mathcal{A}(\omega,t)$ in the collision integrals Rammer (2007). In case of periodically driven systems, however, the fast oscillatory component of the external field can induce fast harmonics on the COM-dependence of all two-time quantities. To make a connection with Fig. 8, one must identify $\psi(t)$ as $\mathcal{A}(\omega,t)$ for a fixed value of $\omega$. In principle, we can resolve the $t$-dependence of $\mathcal{A}(\omega,t)$ into Fourier harmonics of the drive using short-time Fourier transforms:

[TABLE]

where formally $\mathcal{A}_{n}(\omega;t)=\int_{-\infty}^{+\infty}\mathrm{d}t^{\prime}\,W(t^{\prime}-t)\,e^{in\Omega t^{\prime}}\,\mathcal{A}(\omega,t^{\prime})$. We express the COM-time of the harmonics as a label, i.e. $\mathcal{A}_{n}(\omega;t)$, to emphasize on the quasi-static nature of the Fourier amplitudes. Here, $\Omega$ is the principal frequency of the driving term. In the present context, $\Omega=2\Omega_{\mathrm{drv}}$ since the external drive appears as $|F(t)|^{2}$ in the model. $W(t)$ is a window function normalized to unity and concentrated near $t=0$ with support in the time interval $\sim(-\tau_{W},+\tau_{W})$ where $\Omega^{-1}\ll\tau_{W}\ll|F_{\mathrm{env}}(t)|/|F^{\prime}_{\mathrm{env}}(t)|$. The shape of the window function is immaterial so long as this condition is satisfied Mallat (1999). We refer to the collection of $\{\mathcal{A}_{n}(\omega,t)\}$ as the Floquet-Wigner transform of $\mathcal{A}(t_{1},t_{2})$ and hereafter, we assume that the Floquet components admit a controlled series expansion in COM time derivatives.

Floquet-Groenewold-Moyal (FGM) product formula— The convolution integrals appearing on the right hand side of KB equations can be formally expressed as a series expansion using Groenewold-Moyal (GM) product formula Groenewold (1946); Moyal (1949):

[TABLE]

The left hand side represent the convolution integral of $\mathcal{A}$ and $\mathcal{B}$ followed by a Wigner transform. Expanding the exponentiated differential operator and truncating the series at finite orders yields an approximate expression for the convolution of $\mathcal{A}$ and $\mathcal{B}$ in terms of the time and frequency derivatives of their Wigner transforms. In particular, truncating the series at the first order yields the well-known gradient expansion formula which forms the basis of quantum kinetic equations Kadanoff and Baym (1962). This procedure can be readily generalized to the case of Floquet-Wigner transforms. To this end, we plug the Floquet-Wigner expansion of $\mathcal{A}$ and $\mathcal{B}$ (as given in Eq. 39) in Eq. (40) and take a short-time Fourier transform of both sides. The final result is:

[TABLE]

The COM time derivatives have been resolved into a part acting on the phasor and a part acting on the Fourier amplitudes, i.e. $\partial_{t}\rightarrow-in_{\mathcal{A}/\mathcal{B}}\Omega+\partial^{\mathcal{A}/\mathcal{B}}_{t}$. Performing the $t^{\prime}$-integral is trivial since by construction, we can neglect the COM time variations $\{\mathcal{A}_{n}(\omega;t\}$ and $\{\mathcal{B}_{n}(\omega;t\}$ within the support of the window function. Using the formal Taylor’s expansion formula $e^{\alpha x}f(x)=f(x+\alpha)$, we finally obtain:

[TABLE]

where we have defined:

[TABLE]

We refer to labels $n$ and $m$ in $\mathcal{A}_{n,m}(\omega;t)$ as Floquet band index and Floquet quasi-momentum, respectively, in analogy to the Bloch band theory. Eq. (42) will be referred to as Floquet-Groenewold-Moyal (FGM) product formula. Expanding the exponentiated differential operator to linear order, we obtain the FGM gradient expansion formula.

IV.4 Floquet-Boltzmann quantum kinetic equations for phonons

The formalism outlined in the previous section can be used to obtain quantum kinetic-like (“one-time”) equations for the nonequilibrium evolution of the lattice displacement and phonon propagators in the presence of a periodic drive with a slowly varying envelope.

Starting with Eq. (29), taking a short-time Fourier transform of the sides and neglecting the second-order derivatives, we find:

[TABLE]

Next, we consider the retarded component of the KB equation for $\mathcal{D}$. Taking a FW transform of the sides of Eqs. (31)-(31) and using the FGM gradient expansion formula in the collision convolution integrals, we find:

[TABLE]

where $\mathfrak{D}^{-1}(\omega)\equiv(\omega^{2}-\omega_{0}^{2})/(2\omega_{0})$ is the bare inverse phonon propagator in the frequency time. $V(t)$ was defined earlier in Eq. (131) and is the self-consistently determined potential that parametrically drives the phonons. The Fourier amplitudes of $V(t)$ are trivially found as $V_{n}(t)=-(\kappa_{4}/2)\sum_{n_{1}}\varphi_{n_{1}}(t)\,\varphi_{n-n_{1}}(t)-(\kappa_{4}/2)\,\chi_{n}(t)-\kappa_{3}\,\varphi_{n}(t)$. We have employed local approximation for $\Pi$. The phonon propagator thus can be thought of either as that of a single momentum-mode or as the momentum-summed (local) one. We do not need a separate evolution equation for the advanced propagator since it can be determined from the retarded propagator via the identity $\mathcal{D}^{A}_{n}(\omega;t)=[\mathcal{D}^{R}_{-n}(\omega;t)]^{*}$.

The kinetic equations for the Keldysh phonon propagator is found in a similar fashion. First, we use the Langreth rules to find an explicit KB equation for $\mathcal{D}^{K}$, followed by a FW transform and FGM gradient expansion of the collision integral convolutions. The final result is:

[TABLE]

Even though Eqs. (45)-(46) have a more complex presentation compared to the original KB equations, they are significantly simpler to work with in practice: convolution integrals have been reduced to discrete Floquet index summations, and the two-time structure has been reduced to the COM time and the relative frequency $\omega$ which does not appear in a convolution. Finally, we remark that the electron-hole bubble $\Pi_{\ell}$ acts as a dissipation source (“bath”) for phonons, and must be determined by solving the evolution equations for $\hat{\mathcal{G}}$ in a fully self-consistent treatment. We will argue later that it can be approximately calculated using bare electron propagators at the initial temperature as long as heating does not bring up the energy density of electrons to phonon energy scales. Approximate expressions for $\Pi_{\ell}$ have been provided in Sec. E.1 (see Eq. E.1). We show $\Pi_{\ell}$ acts as an quantum Ohmic bath for phonons, and gives rise to a Lamb shift of the Einstein oscillator.

IV.5 Floquet-Boltzmann quantum kinetic equations for electrons

We can obtain quantum kinetic equations for the nonequilibrium evolution of electrons in a similar to phonons. Before embarking on deriving these equations, we take a short detour to derive explicit expressions for the Migdal-Eliashberg and employ the approximations discussed in Sec. IV.2.

IV.5.1 Migdal-Eliashberg self-energy: the general case

Our goal in this section is derive explicit expressions for various real-time components of the Migdal-Eliashberg self-energy, starting from Eq. (34). To keep the notation uncluttered, we will work out the results for a general Wigner-transformed propagators as a first step. We employ the approximations discussed in Sec. IV.2 step by step. The simplified self-energy expressions for Floquet-Wigner-transformed propagators are readily found in the end as an special case.

As a starting point, we transform Eq. (34) to the Wigner representation and employ the decomposition of propagators defined in Eqs. (104a)-(104b) in terms of their Keldysh and spectral components. The calculation is elementary and the final result is:

[TABLE]

where the capped Nambu propagators are defined according to Eq. (106). We define the Eliashberg function 222In the Migdal-Eliashberg theory literature, it is customary to refer to az $F^{\rho}$ as $\alpha^{2}F$. We do not employ this cluttered notation in our treatment. $F^{\rho}_{\xi,\xi^{\prime}}(\nu)$, as well as a Keldysh Eliashberg function $iF^{K}_{\xi,\xi^{\prime}}(\nu)$ as:

[TABLE]

In the special case of Einstein phonons where $i\mathcal{D}^{K}_{\mathbf{q}}$ and $\rho_{\mathbf{q}}$ have no $\mathbf{q}$-dependence, we find:

[TABLE]

where the energy-resolved dimensionless coupling constant $\alpha^{2}(\xi,\xi^{\prime})$ is defined as:

[TABLE]

which in turn in the limit of flat EDOS and constant ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$ evaluates to $\alpha^{2}(\xi,\xi^{\prime})\rightarrow\nu(0)\,g^{2}\sim\mathrm{const}$. The Fermi-surface averaged (FSA) self-energy is readily found as:

[TABLE]

Since the FSA self-energy has no momentum dependence, the $\mathbf{k}$-dependence of the resulting electron propagators is induced from the bare electron dispersion $\xi_{\mathbf{k}}$. Hence, we have legitimately replaced the $\mathbf{k}^{\prime}$ momentum labels with $\xi^{\prime}$, and $\mathbf{k}^{\prime}$ momentum summations with EDOS-weighted $\xi^{\prime}$ integrals over the bandwidth $[-W_{\mathrm{el}}/2,+W_{\mathrm{el}}/2]$. At this point, we employ the remaining approximations discussed in Sec. IV.2, i.e. flat EDOS, infinitely large electronic bandwidth $W_{\mathrm{el}}$ compared to the phonon scale, constant ${g}^{\phantom{\dagger}}_{\mathbf{k},\mathbf{k}^{\prime}}$, and dispersionless phonons. In particular, the last two imply that the Eliashberg functions do not depend on $\xi^{\prime}$. The final result reads:

[TABLE]

where:

[TABLE]

are local electronic spectral and Keldysh functions as obtained by summing over all momentum states:

[TABLE]

Finally, we find it useful to parametrize the strength of the electron-phonon coupling in terms of the dimensionless mass enhancement factor of an ideal Einstein oscillator at equilibrium Bennemann and Ketterson (2008):

[TABLE]

IV.5.2 Migdal-Eliashberg self-energy: the normal state

The results of the previous section were worked out for a general Nambu electron propagator. We specialize the result to the normal non-paired state in this section. A Nambu functions are diagonal in the normal state and the Nambu structure of the self-energy and the ensuing KB equations can be simplified. Starting with the general ansatz,

[TABLE]

and using Lemma 1(b) from Appendix B, we find:

[TABLE]

where we have defined $\mathsf{A}_{\xi}(\omega,t)=-2\,\mathrm{Im}\,\mathcal{G}_{\xi}^{R}(\omega,t)$ as the normal spectral function. The kinetic energy variable $\xi$ appears as a convenient scalar surrogate for $\mathbf{k}$ after using FSA self-energies (see the discussion after Eq. 51). We have also used Lemma 1(d) to relate the time-reversed spectral function $\bar{\mathsf{A}}_{\xi}(\omega,t)\equiv\mathrm{Im}\,\bar{\mathcal{G}}_{\xi}^{R}(\omega,t)$ to $\mathsf{A}_{\xi}(\omega,t)$, i.e. $\bar{\mathsf{A}}_{\xi}(\omega,t)=\mathsf{A}_{\xi}(-\omega,t)$. We further define an unrestricted ansatz for $i\mathcal{G}_{\xi}^{K}$ in compliance with Lemma 1(c) and 1(e):

[TABLE]

where $\psi(\omega,t)$ is an odd real function of $\omega$ that encodes the statistics of electrons in the normal state. For example, in thermal equilibrium, the Kubo-Martin-Schwinger (KMS) boundary condition implies:

[TABLE]

We further define the local electron statistics as:

[TABLE]

In theory, the $\Sigma^{R/K}$ can be expressed as functionals of $\psi$, $i\mathcal{D}^{K}$, and $\rho$. An explicit formula for $\Sigma^{R/K}[\psi,i\mathcal{D}^{k},\rho]$ can be found using Eqs. (52)-(52), Eqs. (54a)-(54b), Eq. (58), and the following crucial lemma:

Lemma ($\mathcal{G}^{R}$ momentum summation formula): Assuming (1) infinite electronic bandwidth limit, and (2) a momentum-independent retarded self-energy as in Fermi surface averaging approximation and local approximation (DMFT), the following identity holds:

[TABLE]

The proof is given in Appendix C. An immediate corollary is:

[TABLE]

Combining Eqs. (52)-(52), Eqs. (54a)-(54b), Eq. (58), and Eq. (61), we find the sought after explicit self-energy functionals:

[TABLE]

where:

[TABLE]

It is noteworthy that $\Sigma^{R/K}$ only depends on the local electron statistics $\psi(\omega,t)$ and not the $\xi$-response $\psi_{\xi}(\omega,t)$. Finally, Lemma 3(a) and 3(b) from Appendix B and the above result imply $\bar{\Sigma}^{R/K}(\omega,t)=\Sigma^{R/K}(\omega,t)$, i.e. $\hat{\Sigma}$ is proportional to the identity matrix in the normal state. This result is strictly a consequence of the ideal Migdal limit. One can show $\Sigma^{R}-\bar{\Sigma}^{R}\propto 1/W_{\mathrm{el}}$ in a finite bandwidth model. This completes our discussion of the Migdal-Eliashberg self-energy in the normal state.

IV.5.3 Floquet-Boltzmann kinetic equation for electrons

An explicit Floquet-Boltzmann kinetic equation can be derived for electrons in the normal state using the result of the previous section and the KB equations. To this end, we write the Keldysh component of Eqs. (33)-(33) using FSA self-energy:

[TABLE]

Subtracting the sides of these equations from one another and performing a Wigner transformation, we find:

[TABLE]

Integrating both sides over $\xi$, using Eq. (61), and the definition of the local electron statistics (see Eq. 60), we find a simple evolution equation for $\psi(\omega,t)$:

[TABLE]

As an intermediate consistency check, at equilibrium where $\psi=\tanh(\beta\omega/2)$ and the GM product reduces to an algebraic product, the right hand side evaluates to $i\Sigma^{K}(\omega)+2i\mathrm{Im}[\Sigma^{R}]\tanh(\beta\omega/2)$ which vanishes in light of the KMS boundary condition. Thus, the thermal state remains stationary as expected.

For a periodically driven system with a slowly varying drive envelope, expanding the convolution integrals using first-order FGM product formula yields:

[TABLE]

The Floquet components of self-energy functionals are worked out easily from Eqs. (63)-(63b):

[TABLE]

where:

[TABLE]

The advanced self-energy is readily obtained using the symmetry relation $\Sigma^{A}_{n,m}(\omega;t)=\Sigma^{A}_{n}(\omega-m\Omega/2;t)=[\Sigma^{R}_{-n}(\omega-m\Omega/2;t)]^{*}=[\Sigma^{R}_{-n,m}(\omega;t)]^{*}$. Self-energies with finite Floquet quasi-momentum, e.g. $\Sigma^{R}_{n,m}(\omega;t)$ are found by shifting $\Sigma^{R}_{n,m}(\omega;t)\equiv\Sigma^{R}_{n}(\omega-m\Omega/2;t)$ (see Eq. 43).

In the fully self-consistent scheme, one must integrate Eq. (IV.5.3) together with the previously derived kinetic equations for the lattice displacement and phonon propagator self-consistently. We note that Eq. (IV.5.3) is an implicit integro-differential equation for $\partial_{t}\psi_{n,m}$ in disguise due to the presence of derivative terms $\partial_{t}\Sigma^{R/A}[\psi]$. We will discuss the numerical approach for solving this equation in Sec. F.

V Migdal-Eliashberg theory of Floquet superconducting instability

We derived a set of tractable evolution equations for the driven system in the normal state. In this section, we derive a criterion to identify the instability of the normal state toward forming a Floquet superconducting state. This criterion follows from a careful linear response analysis as follows: we introduce a small fictitious pairing potential (i.e. an off-diagonal self-energy term) to the time-dependent self-energy obtained in the normal state: $\hat{\Sigma}^{R}(\omega,t)\rightarrow\Sigma^{R}(\omega,t)\mathbb{I}+\hat{\phi}(\omega,t)$. Here, $\hat{\phi}(\omega,t)$ is off-diagonal in the Nambu space. The off-diagonal self-energy, in turn, induces an anomalous (off-diagonal) propagator $\delta\hat{\mathcal{F}}[\hat{\phi}]$ which in turn generates the pairing potential. The introduced pair potential may only persist if and only if $\hat{\Sigma}^{R}[\hat{\mathcal{G}}+\delta\hat{\mathcal{F}}[\hat{\phi}]]-\hat{\Sigma}^{R}[\hat{\mathcal{G}}]\equiv\hat{\phi}$. This procedure is shown diagrammatically as:

[TABLE]

Since the self-energy is a linear functional of $\hat{\varphi}$, satisfiability of the above equation for a non-trivial $\hat{\phi}$ requires the linear operator $\mathbb{I}-\mathcal{D}\mathcal{G}\mathcal{G}$ to have a non-trivial null space. This operator is precisely the inverse vertex operator that appears in the two-particle propagator $\sim\mathcal{G}\mathcal{G}(\mathbb{I}-\mathcal{D}\mathcal{G}\mathcal{G})^{-1}$. This, the pairing condition is formally equivalent to requiring a zero-energy pole in the two-particle propagator, the well-known Thouless criterion for spontaneous symmetry breaking Thouless (1960). The precise condition for a driven nonequilibrium system is complicated due to nonequilibrium propagators and Floquet bands, and requires a careful implementation of the outlined steps, which is the goal of the next sections. As a first step, we will work out the pairing instability criterion for an arbitrary normal state. The results will be used to find the pairing condition for quasi-steady Floquet states.

V.1 Pairing instability criterion for arbitrary nonequilibrium states

We start the analysis by revisiting the KB equations for the retarded Nambu propagator using the FSA self-energy and in the Wigner representation:

[TABLE]

As discussed in the earlier remarks, we assume the following ansatz at the onset of pairing:

[TABLE]

where $\hat{\phi}$ is the infinitesimal pairing potential, $\hat{\mathcal{G}}^{R/K}_{\xi}$ and $\Sigma^{R}\mathbb{I}$ denote the unperturbed Nambu propagators and self-energy in the normal state, respectively, $\psi_{\xi}$ is the electron statistics in the unperturbed normal state (see Eq. 58, and $\delta\hat{\mathcal{F}}^{R}_{\xi}[\phi]$ and $\delta\hat{\mathsf{A}}_{\xi}[\phi]$ denote the off-diagonal linear response of the retarded propagator and the spectral function, respectively. Using Lemma 1(b) of Appendix B, we have:

[TABLE]

Here, we have assumed that the pairs formed at the onset of transition have the same statistics $\psi_{\xi}$ as the normal electrons. Inserting the ansatz in Eqs. (73a)-(73b), summing the sides and keeping terms to linear order in $\hat{\phi}$, we find:

[TABLE]

Since $\delta\hat{\mathcal{F}}^{R}_{\xi}$ is fully off-diagonal, \big{\{}\hat{\sigma}_{z},\delta\hat{\mathcal{F}}^{R}_{\xi}\big{\}}=0. Integrating both sides of Eq. (76) over $\xi$, using Eq. (61) to replace the $\xi$-summed normal retarded propagators with the universal value of $-i\pi$, we find:

[TABLE]

Without the loss of generality, we have assumed $\hat{\phi}=\phi\,\hat{\sigma}_{x}$ and $\delta\hat{\mathcal{F}}^{R}[\phi]=\delta\mathcal{F}^{R}[\phi]\,\hat{\sigma}_{x}$ to distill the Nambu matrix structure on the last equation.

Solving the equation above for arbitrary $\Sigma(\omega,t)$ and $\phi(\omega,t)$ is a formidable task due to the intricate differential structure of the GM product formula. However, we will show later that it can be reduced to a simpler algebraic structure using the properties of Floquet-Boltzmann states. Here, we proceed with the general observation that $\delta\mathcal{F}^{R}[\phi]$, and subsequently $\delta\check{\mathsf{A}}[\phi]$ as given by Eq. (75), are computable linear functionals of $\phi$. Projecting out the off-diagonal component of the retarded Migdal-Eliashberg self-energy (Eq. 52), we find:

[TABLE]

This is the sought after self-consistency relation between the pairing potential and the induced anomalous response. Eq. (75) implies $\mathrm{Tr}\left[\hat{\sigma}_{x}\,\delta\check{\mathsf{A}}_{\xi}[\phi](\omega,t)\right]=+4\,\mathrm{Im}\,\delta\mathcal{F}^{R}_{\xi}[\phi](\omega,t)$. Furthermore, Lemma 1(d) of Appendix B implies that this quantity is a real odd function of $\omega$. These considerations allow us to simplify the pairing self-consistency condition, Eq. (78), to:

[TABLE]

where:

[TABLE]

Eq. (79) is a functional eigenvalue equation for $\phi(\omega,t)$. As mentioned earlier, the paired state is stable if and only if Eq. (79) admits a non-trivial solution for $\phi(\omega,t)$.

V.2 Pairing instability criterion for quasi-steady Floquet states

As mentioned earlier, solving Eq. (77) for arbitrary nonequilibrium states is a challenging task and requires resorting to numerical methods in general. This task is significantly simpler in special cases such as stationary states where all time derivatives vanish, or quasi-steady Floquet states where the Fourier amplitudes of all involved quantities are approximately stationary. In both cases, Eq. (77) can be cast into an algebraic equation and be solved either numerically or by perturbation. To study the case of quasi-steady Floquet states, we take a Fourier transform of the sides of Eq. (77) in $t$ and neglect the time derivatives of Fourier amplitudes in convolutions. Physically, the latter is justified if the pair formation rate is faster than the macroscopic time scale over which the quasi-stationary Floquet-Boltzmann state evolves. Replacing the GM products appearing in the right hand side of Eq. (77) with the leading order FGM product formula, we find:

[TABLE]

We note that $\mathsf{Q}^{n,m}_{n^{\prime},m^{\prime}}(\omega;t)$ only depends on the unperturbed retarded self-energy in the normal state. Diagrammatically, this step is equivalent to attaching the pair propagator to $\phi$ an in the right hand side of Eq. (72). This expression, together with Eq. (80a), yield the sought after explicit relation between the gap and the pairing potential:

[TABLE]

Taking a Fourier transform of the sides of Eq. (79) in $t$ and neglecting time derivatives as before, we find:

[TABLE]

where:

[TABLE]

Plugging $\phi_{n^{\prime},m^{\prime}}$ from Eq. (83) into Eq. (82) yields the final functional eigenvalue equation for the Floquet gap at the onset of pairing:

[TABLE]

It is worthwhile to take a moment and study this equation in some detail. We recall that $\mathsf{Q}^{n,m}_{n^{\prime},m^{\prime}}$ is derived from the normal-state self-energy and relates the anomalous response to the pairing potential (see Eq. 81b). Thus, this quantity brings in the physics of quasiparticle propagation in the normal state such as lifetime and residue. On the other hand, the phonon propagator enters through $K_{n,m}$ and therefore, it brings in the retarded phonon-mediated attraction between the quasiparticles (see Eqs. 84 and 71).

Finally, let us consider the static non-driven limit of Eq. (85) as a consistency check where all Floquet indices and summations can be dropped. In this limit, Eq. (81a) admits a simple algebraic solution:

[TABLE]

implying $\mathsf{Q}^{n,0}_{n^{\prime},m^{\prime}}(\omega)=-i\pi\,\delta_{n,0}\,\delta_{n^{\prime},0}\,\delta_{m^{\prime},0}\left[\omega-\Sigma^{R}(\omega)\right]^{-1}$. Plugging this into Eq. (85), we find:

[TABLE]

where $Z(\omega)=[\omega-\Sigma^{R}(\omega)]/\omega$ as it is usually defined in the context of Migdal-Eliashberg theory. This coincides with the result obtained earlier by by Scalapino, Schrieffer, and Wilkins Scalapino et al. (1966) for equilibrium systems. Our result is a proper generalization of the static Migdal-Eliashberg pairing criterion to arbitrary nonequilibrium states [Eq. (79)], and particularly to Floquet states [Eq. (85)].

In practice, we monitor the eigenvalue spectrum of the linear functional posed by Eq. (85) as the system evolves in the normal state. The normal state is deemed unstable as soon as the lowest eigenvalue crosses zero. The same procedure applies to equilibrium states, where one calculates the normal-state self-energy of electrons at different temperatures and decreases the temperature until the lowest eigenvalue crosses zero.

VI Revisiting the problem: results from the Floquet-Migdal-Eliashberg theory

The machinery of Floquet-Migdal-Eliashberg (FME) quantum kinetics allows us to fill out the missing details in the preliminary analysis given in Sec. III. In particular, we can study the role of competing factors such as parametric phonon generation and the heating of electrons in order to assess whether the mechanism laid out in Sec. III persists in transient dynamics.

At this point, we have developed all the necessary tools to solve the problem using the full FME formalism, Fig. 7. As outlined above, the electron-phonon system is initially prepared in an equilibrium normal state with temperature $T_{i}>T_{c}$, where $T_{c}$ is the critical superconducting transition temperature for the given system parameters. The drive is smoothly ramped up according to:

[TABLE]

where $I_{0}$ denotes the intensity of the drive. We restrict our numerical analysis to weak and intermediate couplings where the phononic and electronic quantities can be calculated iteratively as described below.

As a first step, the electrons are assumed to remain in the equilibrium state, effectively providing a fixed-temperature Ohmic quantum bath $\Pi^{(0)}_{\ell}(\omega)$ for the phonons. Explicit expressions for $\Pi^{(0)}_{\ell}$ are given in Appendix E.1). The Floquet-Boltzmann equations for the lattice displacement $\{\varphi_{n}(t)\}$ and phonon propagators $\{i\mathcal{D}^{K}_{n}(\omega;t),\rho_{n}(\omega;t)\}$ are then numerically integrated forward in time as described in Appendix F.1. Subsequently, the Floquet-Boltzmann equations for the energy distribution of electrons $\{\psi_{n}(\omega;t)\}$ are numerically solved as described in Appendix F.2. If deemed necessary, this two-pass iterative calculation is looped until a self-consistent nonequilibrium solution is obtained for both phonons and electrons. For our choice of parameters, we found additional iterations to be unnecessary by the virtue of the large separation of energy scales between phonons and electrons and weak coupling. Finally, we study the FME pairing instability condition throughout the evolution in order to assess whether the normal state exhibits the pairing instability at some time.

We remark that the system is assumed to evolve in the normal state throughout the simulation. The present formalism detects the instability toward Cooper pairing. Studying the full Floquet superconducting gap formation and its dynamics is a more challenging problem and is better suited to be studied via a Floquet extension of the time-dependent Landau-Ginzburg (TDGL) formalism.

We present the results in two stages. As a first step, in order to gain insight into parameter regimes of maximally enhanced superconductivity, we hold the electrons in thermal states with different temperatures (e.g. by coupling them to a large and infinitely efficient heat bath). We proceed by letting the driven phonons settle to a stationary Floquet driven-dissipative state and calculate the $T_{c}$ of electrons on its backdrop. This procedure is similar to calculating $T_{c}$ in equilibrium by studying the eigenvalue spectrum of Eq. (85), however, using the driven phonon propagators.

Finally, we study the full evolution of the coupled system by allowing the electrons evolve on par with phonons. This allows us to investigate the transient nature of the superconducting instability. We find that the heating of electrons destroys the instability at late times and stabilizes the normal state as seen in the experiments Mitrano et al. (2016).

Choice of physical parameters— The trimmed-down model is fully specified by a few physical parameters: electron-phonon coupling $g$, electronic density of states at the Fermi level $\nu(0)$, optical phonon frequency $\omega_{0}$, the phenomenological damping of the coherent lattice displacement $\gamma_{0}$, and the drive coupling strength $\Lambda$. We trade $g$ and $\nu(0)$ with the mass enhancement factor at equilibrium $\lambda_{\mathrm{eq}}$ and local phonon damping $\gamma_{\ell}$. These dimensionless quantities are defined in Eqs. (55) and (E.1) and we quote them again here for reference:

[TABLE]

We set $\lambda=0.5$ and $\gamma_{\ell}=0.2$ in the numerics, which correspond to typical values for fulleride superconductors Varma et al. (1991). We study cubic and quartic nonlinearities separately. The majority of the results are presented for a cubic nonlinearity. As we argued earlier in Sec. III.2, both types of nonlinearity lead to qualitatively similar phenomena. We set $\kappa_{3}=0.1\,\omega_{0}$, $\kappa_{4}=0$ for “cubic results”, and $\kappa_{3}=0$, $\kappa_{4}=0.1\,\omega_{0}$ for “quartic results”. These values are expected to reflect the typical intrinsic lattice nonlinearities.

Normalization constants— Quantities with the dimension of energy are presented in the units of $\Omega_{0}$, the renormalized phonon frequency at equilibrium defined in Eq. (153), an the time axes are scaled with respect to $\tau_{\mathrm{ph}}$, the renormalized period of phonons at equilibrium:

[TABLE]

Here, $\bar{\omega}_{L}$ is the effective Lamb shift of phonons as a matter of coupling to electrons which can be neglected in the weak coupling regime. $U_{0}$ is determined by self-consistently solving the set of equations given in Eq. (F.1) and represents the phonon frequency correction due to lattice nonlinearities. Finally, we parametrize the fully ramped-up drive amplitude by the dimensionless quantity $\mathcal{A}$ defined as:

[TABLE]

and set $\tau_{\mathrm{drv}}=5\,\tau_{\mathrm{ph}}$ in Eq. (88).

Notation— We often discuss period-averaged quantities along with their temporal variances, respectively defined as:

[TABLE]

where $O(t)$ is an arbitrary observable. Note that the effective drive period is $\pi/\Omega_{\mathrm{drv}}$ since the principal harmonic of all observables is $2\Omega_{\mathrm{drv}}$ (see the discussion after Eq. 39). If $O(t)$ is given as a Fourier series with slowly-varying amplitudes, i.e. $O(t)=\sum_{n=-\infty}^{\infty}O_{n}(t)\,e^{2ni\Omega_{\mathrm{drv}}t}$, then $\langle O(t)\rangle=O_{0}(t)$ and $\mathrm{Var}[O(t)]\equiv\sum_{n=1}^{\infty}|O_{n}(t)|^{2}$.

VI.1 Stage I: Driven phonons, thermal electrons

Fig. 9 shows a typical example of phonon propagators subject to nearly-resonant drive in the presence of lattice nonlinearities. Panel (a) shows a heat map of the period-averaged Keldysh phonon propagator $iD^{K}_{n=0}(\nu;t)$ at a function of relative frequency $\nu$ and COM time $t$. The most prominent feature is the red-shift of the phonon peak frequency. The red-shift is a direct consequence of the lattice nonlinearity: with a cubic nonlinearity $\sim\kappa_{3}\,\hat{\varphi}^{3}$, the drive shifts the equilibrium position of the lattice on average, producing a frequency renormalization $\Delta\omega^{2}_{\mathrm{ph}}\sim\kappa_{3}\langle\hat{\varphi}\rangle$. For a nearly-resonant drive, nonlinear effects dominate the value of $\langle\hat{\varphi}\rangle$ such that $\mathrm{sign}[\langle\hat{\varphi}\rangle]=-\mathrm{sign}(\kappa_{3})$. As a result, $\Delta\omega^{2}_{\mathrm{ph}}<0$ regardless of the sign of $\kappa_{3}$ for a strong nearly-resonant drive. In other words, the cubic nonlinearity always softens the lattice. This phenomenon resembles the usual physics of thermal expansion where the drive plays the role of heating.

Panels (b1-2) and (c1-2) show the higher Floquet components of the Keldysh phonon propagator, both of which show emergent features as the drive is ramped up. It is noticed that $|i\mathcal{D}^{K}_{2}|<|i\mathcal{D}^{K}_{1}|$ suggesting that the role of higher Floquet bands become increasingly smaller. Most strikingly, it is noticed that $n>0$ Floquet phonons do not admit a single coherent peak in contrast to the $n=0$ case. As a consequence, Eqs. (45)-(46) do not admit a reliable Boltzmann “quasiparticle” limit, justifying our usage of the more cumbersome quantum kinetic formalism.

One consequence of the drive is parametric generation of phonons Berges and Serreau (2003) and squeezing of lattice momentum fluctuation Knap et al. (2015). Panel (d) and (e) show the evolution of these quantities as a function of COM time $t$. Appendix D shows how these quantities can be calculated from $i\mathcal{D}^{K}(\omega,t)$. As expected, both of these observables significantly increase as the drive is ramped up. The emergence of squeezed states is a well-known hallmark of parametrically driven harmonic oscillators. Finally, panel (f) shows the instantaneous mass enhancement factor $\lambda(t)$ defined as:

[TABLE]

This quantity plays a similar role in the Migdal-Eliashberg theory as $\nu(0)U(t)$ in the BCS theory (the latter was defined in Sec. III). For instance, $T_{c}\approx\omega_{0}\,e^{-1/\lambda_{\mathrm{eq}}}$ in the equilibrium Migdal-Eliashberg theory Bennemann and Ketterson (2008). It is noticed that $\lambda(t)$ has a strong oscillatory component. The period-average of $\lambda(t)$ during one drive period and the lower and upper envelopes are shown as green, blue, and red solid lines.

The Floquet superconducting transition temperature $T_{c}^{\mathrm{Floq}}$ can be determined for each choice of $\Omega_{\mathrm{drv}}$ and $\mathcal{A}$ by calculating the eigenvalue spectrum of the FME gap functional (Eq. 85, $\mathrm{eigs}_{\Delta}$, for different electronic temperatures and locating the first zero-crossing of the lowest eigenvalue $\mathrm{min}~{}\mathrm{eigs}_{\Delta}$.

Fig. 10 shows the results separately for cubic and quartic nonlinearities. The top row corresponds to a calculation using full Floquet phonon propagators. The bottom row is obtained using only $n=0$ (period-averaged) phonon propagators. Strong driving near parametric resonances lead to the instability of the lattice due to nonlinearities. This stems of the our choice of $\mathcal{V}^{\mathrm{ph}}(\varphi)$ (see Eq. 2) which is only valid for low-amplitude deformations and becomes unphysical for large deformations. These unstable regions are hatched in the heat map plots and indeed coincide with the three first parametric resonances at $\Omega_{\mathrm{drv}}/\Omega_{0}\approx 1/2,1/4,1/8$.

The full Floquet result (top row) shows a dramatic enhancement of $T_{c}$, reaching beyond three times the equilibrium value near the resonances. Neglecting the ac components, only a moderate enhancement of $T_{c}$ is found, and only away from the resonances. In particular, $T_{c}$ is suppressed above the main resonance in the static approximation in contrast to the full Floquet result. This finding is strikingly similar to the analytic Floquet BCS analysis of Sec. III.1; see Fig. 6.

To shed light on this finding, we have plotted the density of excited phonons $\langle\Delta n_{\mathrm{ph}}(t)\rangle$ as well as the mean and variance of $\lambda(t)$ during a drive period as a function of $\Omega_{\mathrm{drv}}$ and $\mathcal{A}$ in Fig. 11. It is noticed that (1) $\langle\lambda(t)\rangle$ is enhanced and suppressed below and above the main resonance, respectively, similar to the analysis of Sec. III and as summarized in Fig. 5; (2) both $\langle\Delta n_{\mathrm{ph}}\rangle$ and $\mathrm{Var}[\lambda(t)]$ are increase significantly above the main resonance. Neglecting $n>0$ components of the phonon propagator and neglecting $\mathrm{Var}[\lambda(t)]$ go hand in hand. The high density of phonon excitations and the suppression of $\langle\lambda(t)\rangle$ above the resonance imply decreased quasiparticle lifetime and electron-electron attraction, respectively, both of which are unfavorable for Cooper pairing. This explains suppression of $T_{c}$ above the main resonance in the static approximation. Away from the resonances, the moderate enhancement of $\langle\lambda(t)\rangle$, which has its roots in the phonon frequency red-shift and is present in the static approximation, explains the moderate enhancement of $T_{c}$. Finally, we remark that including $n>0$ components brings in large $\mathrm{Var}[\lambda(t)]$. In the example shown in Fig. 9(f), the upper envelope of $\lambda(t)$ is nearly three times as large as $\lambda_{\mathrm{eq}}$. As we argued earlier in Sec. III.1, $T_{c}$ is a convex functional of the interaction parameter in the weak coupling limit such that temporal variation of interaction can increase $T_{c}$ even if the period-average remains fixed or even decrease.

VI.2 Stage II: Driven phonons, evolving electrons

The notion of superconducting transition temperature is only well-defined in thermal states. Once the electrons are allowed to evolve as a matter of coupling to phonons, a different diagnostic will be needed to assess the enhancement or suppression of the superconducting transition. Here, we attempt to model a realistic experimental scenario: we prepare the electron-phonon system in $T_{i}>T_{\mathrm{eq}}$, ramp up the drive, and calculate the ensuing nonequilibrium dynamics of phonons and electrons. Since $T_{i}>T_{c}$ in the beginning, the lowest eigenvalue of the FME gap functional begins as a positive value. Whether it remains positive throughout the evolution or crosses zero at some point is our diagnostic. This allows us to study the instability of the normal state toward Cooper pair formation but does not describe the physics of gap formation. The latter can be addressed with an extended formalism based on the the present developments.

Fig. 12 shows an example of the evolution of the Floquet components of the $\xi$-summed (local) energy statistics of electrons $\{\psi_{n}(\omega;t)\}$ defined in Eq. (60), along with their period-averaged scattering rate $\langle\Gamma(\omega,t)\rangle\equiv-2i\,\mathrm{Im}[\Sigma^{R}_{n=0}(\omega;t)]$. The prominent features are (1) the heating of electrons, manifested as the decreased slope of $\psi_{0}(\omega;t)$ at $\omega=0$ as the drive is ramped up, (2) emergence of electrons in $n>0$ Floquet bands, and (3) increased (decreased) spectral broadening (lifetime) of quasiparticles.

Fig. 13 shows the evolution of the lowest eigenvalue of the FME gap functional, $\mathrm{min\,eigs}_{\Delta}$, for $\Omega_{\mathrm{drv}}/\Omega_{0}=0.3,0.4,0.5,0.6$ and $\mathcal{A}=0.25,0.50,0.75$. The nonlinearity is cubic and the choice of physical parameters is as described in the introductory remarks of this section, implying $T^{\mathrm{eq}}_{c}\simeq 0.034\,\Omega_{0}$. The initial temperature chosen as $T_{i}=0.04\,\Omega_{0}\simeq 1.2\,T_{c}^{\mathrm{eq}}$. The dashed lines show the hypothetical case where the electrons are kept at $T_{i}$ (no heating). Red segments indicate where $\mathrm{min\,eigs}_{\Delta}<0$. The most favorable outcome occurs for lower frequency driving, e.g. $\Omega_{\mathrm{drv}}=0.3\,\Omega_{0}$, where the pairing instability persists for a long time. In all cases, heating of electrons tends to stabilize the normal state with long enough driving. This is most easily noticeable for $\Omega_{\mathrm{drv}}=0.4\,\Omega_{0}$ where the instability is confined to a short interval. For $\Omega_{\mathrm{drv}}=0.5\,\Omega_{0},0.6\,\Omega_{0}$, we find $\mathrm{min\,eigs}_{\Delta}>0$ at all times for all three drive strengths. The strong heating of electrons prohibits pairing even though in the absence of heating (dashed lines), pairing would have ensued.

The desirability of lower frequency driving for enhancing the pairing instability can be understood by appealing to the different nature the two competing effects, parametric amplification of the retarded response on one hand, and parametric phonon generation on the other hand. As we argued early on, the former is the main mechanism for enhancing $T_{c}$ and the latter is the main suppressant. Heating of electrons and the decreased coherence of quasiparticles are both consequences of the interaction with the generated high-energy phonons.

Parametric generation of phonons is an on-shell process. For low frequency driving, phonons may only be generated through accumulation of multiple energy quanta from the drive. These higher order processes, however, become increasingly less probable. In contrast, (1) the retarded response does not need to satisfy an on-shell energetic condition, and (2) the cascade of low-frequency parametric resonances at $\Omega_{0}/(2n)$ extend the range of parametric amplification to very low frequencies. For a fixed driving strength and finite damping, the infinite cascade of parametric resonance “tongues” will be truncated at a certain lower frequency. This is exemplified in Fig. 4 in which only two resonances are present, or in Fig. 10 where only three lattice instability tongues are found. Nevertheless, the presence of even a few higher order resonances enables the amplification of the retarded response for reasonably low frequency drives.

To substantiate these arguments with numerical results, we have plotted the evolution of $\langle\Delta n_{\mathrm{ph}}(t)\rangle$ and $\lambda(t)$ in Fig. 14, and electronic effective mass $\langle m^{*}(t)\rangle$ and damping $\langle\Gamma(t)\rangle$ in Fig. 15. As before, time-averaging is performed during one drive period. We notice that for $\Omega_{\mathrm{drv}}=0.3\,\Omega_{0}$ both $\langle\Delta n_{\mathrm{ph}}(t)\rangle$ and $\langle\Gamma(t)\rangle$ remain nearly two orders of magnitude smaller than the on-resonance drive $\Omega_{\mathrm{drv}}=0.5\,\Omega_{0}$. In contrast, $\langle m^{*}(t)\rangle$ and $\lambda(t)$ are at most a factor of four smaller. Thus, we indeed expect a negligible undesirable suppression while still benefiting from parametric amplification of $\lambda(t)$.

The energy distribution of the electrons at $t=30\,\tau_{\mathrm{ph}}$ is shown in Fig. 16. It is noticed that the distribution barely changes from the initial thermal state for $\Omega_{\mathrm{drv}}=0.3\,\Omega_{0}$, in agreement with the aforementioned arguments regarding suppressed parametric phonon generation below the resonance. In all cases, a decent Fermi-Dirac fit can be obtained. For a on-resonant driving frequency $\Omega_{\mathrm{drv}}=0.5\,\Omega_{0}$, the effective temperature reaches $T_{f}^{\mathrm{eff}}\approx 0.22\,\Omega_{\mathrm{drv}}$ (see the figure caption). Even in such cases, the energy density of electrons remains low enough, obviating performing a self-consistency feedback loop to phonons (see Fig. 7).

Finally, we note that the light-induces superconducting state is expected to persist beyond the predictions of the present analysis in the experiments. On the one hand, the formation of a superconducting gap leads to increased coherence of electrons and suppress scatterings. On the other hand, the bulk electrons and phonons that remain unaffected by the pump pulse act as a low temperature heat bath for the driven subsystem and keep it from excessive heating. Exploring these aspects of the problem is left for future works and is further discussed in Sec. VII.

VI.3 Predictions for time-resolved angle-resolved photo-emission spectroscopy experiments (tr-ARPES)

Up to this point, we have developed the theory of parametric amplification of electron-phonon coupling within a model that is general and material-agnostic. Whether or not, and how, the mechanism is realized in a specific material depends on a number of factors such as the phonon spectra, strength of nonlinearities, and the selection rules that dictate and presence or absence of the required nonlinear phonon couplings. These questions can be investigated either by performing ab initio calculation or through further experimental scrutiny. The interaction of electrons with periodically driven phonons will necessarily induce a certain degree of periodicity in electronic observables. This can be see, for example, in Fig. 12 where the energy statistics develop Floquet components as the drive is ramped up. Here, we will show that the signal measured in tr-ARPES experiments will help reveal the formation of Floquet bands.

The tr-ARPES signal can be theoretically calculated from the lesser electron Green’s function Sentef et al. (2013) as follows:

[TABLE]

where $\sigma_{\mathrm{pr}}$ is the temporal resolution of the probe field which generically satisfies $\sigma_{\mathrm{pr}}\gg\Omega_{\mathrm{drv}}^{-1}$. The Gaussian window functions thus simply serve as picking up the “period-averaged” lesser Green’s function, which coincides with $n=0$ Floquet component of $\mathcal{G}$. Thus,

[TABLE]

Up to transitory effects, the lesser Green’s function can be calculated as $\mathcal{G}^{<}=\mathcal{G}^{R}\star\Sigma^{<}\star\mathcal{G}^{A}$ Danielewicz (1984). To leading order, we may approximate convolutions with zeroth-order FGM product formula. Each of the required ingredients for calculating $\mathcal{G}^{<}$ can be obtained using the results already available to us. As a first step, we solve the Dyson’s equation for $\mathcal{G}^{R}_{n;\xi}(\omega;t)$:

[TABLE]

Having calculated $\Sigma^{R}_{n,m}(\omega;t)$ previously from solving the Floquet-Boltzmann for electrons, $\mathcal{G}^{R}_{\xi;n,m}(\omega;t)$ is obtained by truncating the above equations in Floquet bands and solving it as a proper linear system. The advanced Green’s function is found immediately using the identity $\mathcal{G}^{A}_{\xi;n,m}(\omega;t)=\mathcal{G}^{R*}_{\xi;-n,m}(\omega;t)$. By definition, the lesser electron self-energy $\Sigma^{<}$ is related to $\Sigma^{R/A/K}$ as $\Sigma^{<}=(\Sigma^{K}-\Sigma^{R}+\Sigma^{A})/2$. Taking a FW transform, we find:

[TABLE]

where we have used the identity $\Sigma^{A}_{n,m}(\omega;t)=\Sigma^{R*}_{-n,m}(\omega;t)$. With the knowledge of $\Sigma^{<}_{n,m}$, $\mathcal{G}_{\xi;n,m}^{R}$, and $\mathcal{G}_{\xi;n,m}^{A}$, $\mathcal{G}^{<}_{\xi;n=0,m=0}$ is calculating by employing the zeroth-order FGM formula twice. The final result is:

[TABLE]

where we have dropped the common $(\omega;t)$ arguments for brevity. This procedure is formed for a range of $\xi$, and Eq. (95) is used to find the intensity of the tr-ARPES signal.

Fig. 17 shows an example of the predicted tr-ARPES signal for $\Omega_{\mathrm{drv}}=0.4\,\Omega_{0}$ and $\mathcal{A}=0.75$ as a function of $\omega$ and $\xi$. We have traded the momentum transfer with $\xi\equiv\varepsilon_{\mathbf{k}}-\varepsilon_{F}$, the kinetic energy measures from the Fermi surface. The period-averaged quasiparticle dispersion $E_{\xi}$ is obtained by locating the main coherent peak of the period-averaged retarded propagator by solving $\mathrm{Re}\,[\mathcal{G}^{R}_{\xi;n=0}(E_{\xi};t)]^{-1}=0$ and is shown as white dashed lines. Before the drive is ramped up, the signal matches what is expected from a coupled electron / optical phonon system at equilibrium Bennemann and Ketterson (2008): filled states for $\omega<0$, decreased quasiparticle coherence at $\omega=\pm\Omega_{0}$, and a larger effective mass for $|\omega|<\Omega_{0}$. The effective mass is obtained as $m/\langle m^{*}(t)\rangle=\mathrm{d}E_{\xi}(\omega=0)/\mathrm{d}\xi|_{\xi=0}$ and is shown in the second column of Fig. 15 for the same drive parameters as Fig. 17.

An intriguing feature of Fig. 17 is the progressive formation of electronic Floquet bands as the drive is ramped up. The dynamical formation of Floquet bands in driven systems has been observed before experimentally in tr-ARPES spectroscopy of topological insulators Wang et al. (2013). The frequency spacing between the emerging Floquet bands is set by $2\Omega_{\mathrm{drv}}$. Therefore, the ARPES experiment along with the available spectroscopic measurements of the phonon spectra will inform about the origin of the persistent drive after the pump pulse is gone. We remark that the higher Floquet bands might be challenging to observe from noisy measurements due to the small weight of these extra features. For example, at $t=15\,\tau_{\mathrm{ph}}$ and for the strong drive parameters used in Fig. 17, the intensity of the first Floquet band is nearly four orders of magnitude smaller than the main quasiparticle peak.

VII Conclusions and Outlook

In this paper, we studied the parametric resonances of driven nonlinear lattices and discussed its role in enhancing the effective phonon-mediated electron-electron attraction. We presented the analysis in two stages. First, we gave a qualitative and intuitive account using perturbation theory, classical dynamics, and the Floquet BCS theory in Sec. III in order to elucidate the mechanism of parametric amplification of phonon-mediated Cooper pairing. Next, we developed a quantum kinetic formalism based on an extension of the Migdal-Eliashberg theory to driven systems and nonlinear lattices in Sec. IV.3 and revisited the problem one more time and in full detail in Sec. VI. The numerically tractable quantum kinetic formalism allowed us to study both the intricate transient and long-time dynamics of the system following the pump pulse. In particular, we investigated the role of parametric phonon generation and subsequent heating of electrons in destroying the transient superconducting instability. Finally, we predicted the transient formation of electronic Floquet bands as an experimentally observable consequence of parametrically driven phonons in Sec. VI.3. This prediction can be tested in time-resolved ARPES experiments and can be used to better understand material-specific mechanisms of parametric drive generation.

This work can be extended in several directions. So far, we have studied the evolution of the system in the normal-conducting state and treated Cooper pairing as an instability. An important extension of this work is to take into account dynamical symmetry breaking and the formation of the Floquet superconducting gap. This can be done most naturally by generalizing the Floquet-Boltzmann kinetic equation of electrons to symmetry broken states and deriving a time-dependent Ginzburg-Landau theory for the slowly-varying Floquet components of the gap $\{\Delta_{n}(\mathbf{x},t)\}$. This extension allows us to address a broad range of largely unexplored theoretical questions, such as the scaling behavior of the coherently-driven system in the critical regime, and the nature of Kibble-Zurek defects Kibble (1980); Zurek (1996) formed as a result of non-adiabatic preparation of the ordered state. Furthermore, extension to gapped states allows us to calculate the nonequilibrium optical conductivity and make a more direct connection to pump-probe experiments Mitrano et al. (2016). A related problem is the question of the lifetime of the transient superconducting state. We find superconductivity as a transient phenomenon as shown in Fig. 13. It arises when electron-phonon interaction is already enhanced and before electrons have been heated too much. To give a more detailed analysis of the duration of the transient regime, we need to allow for opening of the quasiparticle gap which we expect to make transient superconductivity last longer.

The role of light-induced changes in the screened Coulomb interaction has been recently highlighted in the phenomenology of the light-induced superconductivity in $K_{3}C_{60}$ in Refs. Kim et al. (2016); Mitrano et al. (2016). Furthermore, the shortcomings of Migdal-Eliashberg theory for providing an accurate description of fullerene superconductors and necessity of beyond-Migdal vertex corrections have been indicated in Ref. Cappelluti et al. (2001). Therefore, it is desirable to extend the present formalism to include both Coulomb interaction and beyond-Migdal vertex corrections Grimaldi et al. (1995) and to study their role to the extent relevant to the mechanism discussed in this paper. In equilibrium, the effects of retarded Coulomb interaction can be incorporated in the Migdal-Eliashberg theory using the Morel-Anderson (MA) pseudo-potential Morel and Anderson (1962); Schrieffer (1983); Bauer et al. (2012). A nonequilibrium extension of this result is lacking and must be worked out. A naive application of the equilibrium result suggests that the MA pseudo-potential $\mu^{*}=\mu_{c}/[1+\mu_{c}\log(\varepsilon_{F}/\omega_{\mathrm{ph}})]$ directly decreases the mass enhancement factor, i.e. $\lambda(t)\rightarrow\lambda(t)-\mu^{*}$. Here, $\mu_{c}=\nu(0)U_{c}$ and $U_{c}$ is the typical screened Coulomb interaction between conduction electrons. In this paper, we showed that parametric driving enhances Cooper pairing by increasing $\langle\lambda(t)\rangle$ and its temporal variations. Since the Coulomb interaction does not directly play a role in the parametric resonance of the lattice, we expect our conclusions to remain valid. Moreover, Ref. Kim et al. (2016) suggests that $\mu_{c}$ effectively decreases in the pumped system, in which case, the parametric amplification of $\lambda(t)$ and decreased $\mu^{*}$ both work toward enhancing Cooper pair formation in $K_{3}C_{60}$. The role of dynamical vertex corrections and the status of Migdal’s theorem, in particular in the presence of the external drive, is less clear and must be carefully reassessed via real-time techniques in the spirit of the analysis provided for dynamical electron-mediated nonlinearities in Appendix E.

We note that photo-induced enhancement of superconductivity has also been observed in high-Tc cuprates Mankowsky et al. (2014) along with several theoretical proposals for explaining these experiments Raines et al. (2015); Okamoto et al. (2016); Höppner et al. (2015); Coulthard et al. (2016); Sentef et al. (2016); Patel and Eberlein (2016). Cuprate superconductors are considerably more complicated than conventional electron-phonon superconductors that we considered in this paper. Superconductivity in these materials is likely to be of non-phononic origin and there are several competing orders. However, we expect that the ideas explored here may be relevant for light-enhanced superconductivity in these materials as well. For example, periodic lattice modulation changes the strength of magnetic exchange interactions and may lead to parametrically amplified electron-paramagnon coupling. Paramagnons are expected to play the role of phonons in unconventional superconductors.

Last but not least, another intriguing future research direction which is also of much technological interest, is to extend the present analysis to open driven-dissipative systems along with accurate material-specific ab initio calculations. The transient light-induced superconducting state can be enhanced further or even stabilized by continuous pumping of the lattice and simultaneous cooling. Such a hybrid “pumped-and-cooled” device may operate more efficiently compared to the usual refrigerated superconductor depending on the highest achievable effective critical temperature and the pump absorption power of the material.

VIII Acknowledgments

We thank A. Cavalleri, A. Georges, V. Galitski, C. Kollath, A. Millis, and B. Halperin for useful discussions. MB and GR acknowledge support from the Institute for Quantum Information and Matter (IQIM), an NSF Physics Frontiers Center with support of the Gordon and Betty Moore Foundation. MK acknowledges support from the Technical University of Munich - Institute for Advanced Study, funded by the German Excellence Initiative and the European Union FP7 under grant agreement 291763, and from the DFG grant No. KN 1254/1-1. IM acknowledges support from the Materials Sciences and Engineering Division, Basic Energy Sciences, Office of Science, US Dept. of Energy. ED acknowledges support from Harvard-MIT CUA, NSF Grant No. DMR-1308435, AFOSR Quantum Simulation MURI, and AFOSR MURI Photonic Quantum Matter.

Appendix A The parametrically driven harmonic oscillator

We studied the problem of parametrically driven harmonic oscillator formally in Sec. III in the context of the resonant amplification of phonon response. Some of the technical details were left out and we present them here.

We presented the solution of the Heisenberg equation in terms of four special functions $\mathfrak{M}_{\alpha\beta}(t,t^{\prime})$, $\alpha,\beta=P,Q$ (see Eq. III). These functions can be expressed in terms of even and odd Mathieu functions and their derivatives as follows:

[TABLE]

where $\mathrm{s}(z)\equiv\mathrm{Se}(\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2},-\alpha\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2},z)$ and $\mathrm{c}(z)\equiv\mathrm{Ce}(\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2},-\alpha\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2},z)$ denote the odd and even Mathieu functions with characteristic value $\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2}$ and parameter $-\alpha\omega_{\mathbf{q}}^{2}/\Omega_{\mathrm{drv}}^{2}$, respectively, and the prime sign denotes derivatives with respect to $z$. We showed that $\mathfrak{M}_{QP}(t,t^{\prime})$ is of particular interest and determines the retarded phonon response $\mathcal{D}^{R}_{\mathsf{Q}}(t,t^{\prime})$ (see Eq. III). Here, we present a series expansion of this function in terms of the parameter $\alpha$ (see Eq. III), i.e. $\mathfrak{M}_{QQ}(t,t^{\prime})=\sum_{n=0}^{\infty}a^{n}\,\mathfrak{M}^{(n)}_{QQ}(t,t^{\prime})$. The first two terms in the series are given as:

[TABLE]

The higher order terms are increasingly more complex but can be easily worked out using a computer algebra system.

Appendix B Definition and properties of the CTP Green’s functions

In this appendix, we briefly review the definition of CTP Green’s functions, their various real-time components, and their symmetries.

The CTP Nambu electron propagator is defined as:

[TABLE]

where $\Psi^{\phantom{\dagger}}_{\mathbf{k}}=(c_{\mathbf{k}\uparrow},c_{-\mathbf{k}\downarrow}^{\dagger})^{T}$, $\Psi^{\dagger}_{\mathbf{k}}=(c_{\mathbf{k}\uparrow}^{\dagger},c_{-\mathbf{k}\downarrow})$, and $\hat{\mathcal{G}}_{\mathbf{k}}(t_{1},t_{2})$ is a $2\times 2$ matrix in the Nambu space. Here, $\mathcal{C}=\mathcal{C}^{+}\cup\mathcal{C}^{-}$ denotes the round-trip Keldysh contour where $\mathcal{C}^{+}=[t_{0},+\infty)$ and $\mathcal{C}^{-}=(+\infty,t_{0}]$, and $T_{\mathcal{C}}$ is the fermionic (anti-symmetric) time ordering operator on $\mathcal{C}$. Similarly, the real phonon propagator is defined as:

[TABLE]

where $\varphi_{\mathbf{q}}=b_{\mathbf{q}}^{\dagger}+b_{-\mathbf{q}}^{\phantom{\dagger}}$ is the Fourier transform of the lattice displacement operator, and $T_{\mathcal{C}}$ is the bosonic (symmetric) time ordering operator in $\mathcal{C}$. The lesser ($<$) and greater ($>$) real-time Green’s functions are defined as specific orderings of the two contour times where $t_{1}<_{\mathcal{C}}t_{2}$ and $t_{1}>_{\mathcal{C}}t_{2}$, respectively:

[TABLE]

The retarded ($R$), advanced ($A$), and Keldysh ($K$) propagators are defined as $\mathcal{A}^{R}(t_{1},t_{2})=\theta(t_{1}-t_{2})[\mathcal{A}^{>}(t_{1},t_{2})-\mathcal{A}^{<}(t_{1},t_{2})]$, $\mathcal{A}^{A}(t_{1},t_{2})=-\theta(t_{2}-t_{1})[\mathcal{A}^{>}(t_{1},t_{2})-\mathcal{A}^{<}(t_{1},t_{2})]$, and $\mathcal{A}^{K}(t_{1},t_{2})=\mathcal{A}^{>}(t_{1},t_{2})+\mathcal{A}^{<}(t_{1},t_{2})$, respectively, where $\mathcal{A}$ is either $\hat{\mathcal{G}}$ or $\mathcal{D}$. We define spectral/statical decomposition of lesser/greater electron and phonon Green’s functions as follows:

[TABLE]

These definitions can be thought of as definitions of electron and phonon spectral functions:

[TABLE]

Similar definitions apply to Green’s functions in (Floquet-)Wigner representation, and for momentum-summed Green’s functions. For all Nambu matrix quantities such as $\hat{\mathsf{A}}$, $i\hat{\mathcal{G}}^{K}$, etc., we define capped Nambu matrices as:

[TABLE]

We finish this appendix by listing a number of useful symmetry relations in the Wigner representation.

Lemma 1 (symmetries of Nambu functions): We define time reversal symmetric (TRS) states as being invariant under operation $(\mathbf{k},\uparrow)\leftrightarrow(-\mathbf{k},\downarrow)$. The following identities hold for a TRS state:

[TABLE]

(proof) The proofs are elementary and readily follow from the definitions. The last two identities are less trivial and require a careful examination of the matrix elements of $\hat{\mathcal{G}}_{\mathbf{k}}(\omega,t)$.

Lemma 2 (symmetries of Eliashberg functions): We define an inversion symmetric (IS) state as being invariant under operation $\mathbf{q}\leftrightarrow-\mathbf{q}$. The following identities hold exactly for inversion symmetric states:

[TABLE]

(proof) The proofs are elementary and follow from the definition of Eliashberg functions (Eqs. 48-48) and phonon propagators.

Lemma 3 (symmetries of the Nambu self-energy): The following identities holds for TRS and IS states:

[TABLE]

(proof) Both identities are easily established by calculating $\hat{\Sigma}^{R/K}(-\omega,T)$ using Eqs. (52)-(52), changing integration variables $\omega^{\prime},\nu\rightarrow-\omega^{\prime},-\nu$ and using Lemma 1 and 2 identities to change the sign of the frequencies that appear in the electron and phonon propagators.

Appendix C Proof of $\mathcal{G}^{R}$ momentum summation formula

In this section, we give a proof for $\mathcal{G}^{R}$ momentum summation formula (Eq. 61) using perturbation theory. One of the assumptions of the lemma is the independence of $\Sigma^{R}$ from the momentum variable $\mathbf{k}$. As a result, $\mathcal{G}^{R}$ depends on $\mathbf{k}$ only via the electronic dispersion $\xi_{\mathbf{k}}$. Therefore, we may trade the momentum variable in $\mathcal{G}^{R}$ with $\xi$ without loss of generality. The Dyson series for $\mathcal{G}^{R}$ is:

[TABLE]

where:

[TABLE]

is the non-interacting retarded Green’s function. Let us consider the second term in the series:

[TABLE]

Since $\partial_{t}\mathcal{G}^{R}_{\xi,0}=0$, if in addition we had $\partial_{t}\Sigma^{R}(\omega,t)=0$, we would simply get $[\mathcal{G}^{R}_{0,\xi}]^{2}\,\Sigma_{R}$. Expanding the differential operators in the exponents, it is easily noticed that every $t$-derivative of $\Sigma^{R}$ is accompanied either by $\partial_{\omega}\mathcal{G}^{R}_{0,\xi}=-[\mathcal{G}^{R}_{0,\xi}]^{2}$, or by $\mathcal{G}^{R}_{0,\xi}\partial_{\omega}\Sigma^{R}$. Therefore, derivative corrections due to $t$-dependence of $\Sigma^{R}$ are accompanied by at least one extra power of $\mathcal{G}^{R}_{0,\xi}$. Thus,

[TABLE]

This result is easily generalized to the $n$th term in the Dyson series:

[TABLE]

With this observation, let us integrate the sides of Eq. (110) over $\xi$, considering only the first $n$ terms in the series. The integral over the first term is trivial:

[TABLE]

Using Eq. (113), it is easily shown that the integral over the second term vanishes:

[TABLE]

This result is due to the fact that every term in the expansion of the Groenewold-Moyal series has at least a second order pole. This the same result holds for all higher order terms in the Dyson series. Assuming that the order of limit and integrations can be interchanged, and that the Dyson series converges, we find that the only non-trivial contribution stems from the non-interacting Green’s function. This proves the sought after result:

[TABLE]

Appendix D Calculating phonon excitation density and phonon squeezing from $i\mathcal{D}^{K}$

We formulated the problem in Sec. IV in terms of the real phonon propagator $\mathcal{D}_{\mathbf{q}}(t_{1},t_{2})\equiv-i\left\langle T_{\mathcal{C}}\left[\varphi_{\mathbf{q}}(t_{1})\,\varphi_{-\mathbf{q}}(t_{2})\right]\right\rangle$. While this formulation is convenient and compact, it does not immediately yield useful physical observables such a phonon number $n_{\mathbf{q}}(t)\equiv\langle b_{\mathbf{q}}^{\dagger}(t)b^{\phantom{\dagger}}_{\mathbf{q}}(t)\rangle$ or the anomalous correlations $\kappa_{\mathbf{q}}(t)\equiv\mathrm{Re}[\langle b_{\mathbf{q}}^{\phantom{\dagger}}(t)b^{\phantom{\dagger}}_{-\mathbf{q}}(t)\rangle]$. Here, we show that both quantities can be readily calculated from the Keldysh phonon correlator in the Wigner representation, $\mathcal{D}^{K}_{\mathbf{q}}(\omega,t)$, by performing appropriate frequency integrations. This is enabled by the observation that the interaction and drive terms in the Hamiltonian both commute with $\hat{\varphi}_{\mathbf{q}}$. The only non-commuting term is the lattice kinetic energy. Thus, the Heisenberg equation for $\hat{\varphi}_{\mathbf{q}}$ takes the following simple form:

[TABLE]

We assume the $\mathbf{q}\leftrightarrow-\mathbf{q}$ symmetry in this section and set $\hbar=1$. The last equation allows us obtain $\pi\pi$ correlators by calculating appropriate time derivatives of $\mathcal{D}^{K}$. Defining $\mathcal{P}_{\mathbf{q}}(t_{1},t_{2})\equiv-i\left\langle T_{\mathcal{C}}\left[\pi_{\mathbf{q}}(t_{1})\,\pi_{-\mathbf{q}}(t_{2})\right]\right\rangle$, Eq. 118 immediately implies:

[TABLE]

At equal times, the Keldysh $\varphi\varphi$ and $\pi\pi$ correlators evaluate to a combination of our sought-after observables, $n_{\mathbf{q}}$ and $\kappa_{\mathbf{q}}$:

[TABLE]

which together with Eq. (119) yields:

[TABLE]

In the Wigner representation, $\partial_{t_{1}}\partial_{t_{2}}\rightarrow\nu^{2}$ and we find:

[TABLE]

The anomalous phonon density can be related to phonon squeezing with additional considerations. First, we observe $\mathrm{Im}[\langle b_{\mathbf{q}}^{\phantom{\dagger}}(t)b^{\phantom{\dagger}}_{-\mathbf{q}}(t)\rangle]=0$ in our problem since nonlinearities produce modulation of $\varphi_{\mathbf{q}}\,\varphi_{-\mathbf{q}}$ (as opposed to $\pi_{\mathbf{q}}\,\pi_{-\mathbf{q}}$; see Sec. III). Assuming that a low temperature state is maintained at all times and weak electron-phonon coupling, the phonon state can be approximated as product of two-modes squeezed states of $\pm\mathbf{q}$ on the top of a coherent state for $\mathbf{q}=0$, i.e. $|\Phi(t)\rangle\sim\{\prod_{\mathbf{q}}\exp[\xi_{\mathbf{q}}(t)(b^{\dagger}_{\mathbf{q}}\,b^{\dagger}_{-\mathbf{q}}-b^{\phantom{\dagger}}_{-\mathbf{q}}\,b^{\phantom{\dagger}}_{\mathbf{q}})]\}\,\exp[\sqrt{N}\varphi_{0}(t)(b^{\dagger}_{0}-b^{\phantom{\dagger}}_{0})/2]\,|0\rangle$ throughout the evolution. Here, $\xi_{\mathbf{q}}(t)$ is the momentum-squeezing strength and $\varphi_{0}(t)$ is the coherent displacement. For $\mathbf{q}\neq 0$, this ansatz provides the following relation between the squeezing parameter $\xi_{\mathbf{q}}(t)$ and the anomalous phonon density:

[TABLE]

For weak nonlinearities, the squeezing is also weak $|\xi_{\mathbf{q}}(t)|\ll 1$ and we find $\xi_{\mathbf{q}}(t)=\kappa_{\mathbf{q}}(t)+\mathcal{O}(\kappa_{\mathbf{q}}^{3})$. Thus, the anomalous phonon density directly yields the squeezing parameter. We have shown the period-averaged $\langle\cosh[\xi_{\mathbf{q}}(t)]\rangle-1\approx 2\langle\kappa_{\mathbf{q}}^{2}(t)\rangle$ in Fig. 9(d). It is noticed that squeezing significantly increases as the drive is ramped-up, consistent with the physics of the parametrically driven harmonic oscillator.

Appendix E Electron-mediated phonon dissipation and nonlinearities

The evolution equation for phonon propagators was derived in Sec. IV.1 as well as their counterparts in the Floquet-Boltzmann kinetic approximation in Sec. IV.4. So long as the evolution of phonons is concerned, electrons play the role of a quantum bath through memory convolution integrals $\Pi_{\mathbf{q}}\star\mathcal{D}_{\mathbf{q}}$ and $\mathcal{D}_{\mathbf{q}}\star\Pi_{\mathbf{q}}$ appearing in Eqs. (31) and (31), respectively.

In this section, we derive approximate expressions for $\Pi_{\mathbf{q}}$ assuming that the electrons remain in the initial low-temperature degenerate regime. Meanwhile, we also study the contribution of electrons to phonon nonlinearities. Both objectives can be achieved by integrating out the electrons from the Lagrangian $\mathcal{L}[\varphi,\Psi]$ and obtaining a phonon-only effective action $S_{\mathrm{eff}}[\varphi]$. Expanding the effective action in the electron-phonon coupling, we obtain the bath term at the second order. Higher order terms give the electronic contribution to lattice nonlinearities. Since these corrections have a strong dynamical nature, it is conceivable that they could become large when certain resonance conditions are met; indeed we find this to be case. In other words, even though the intrinsic lattice nonlinearities might be small, coupling to conduction electrons effectively produces large nonlinearities in the presence of a near-resonant drive.

We start our discussion with the electron-phonon Lagrangian:

[TABLE]

where $\mathcal{L}_{0}[\varphi]=-\sum_{\mathbf{q}}(2\omega_{\mathbf{q}})^{-1}\,{\varphi}^{\phantom{\dagger}}_{\mathbf{q}}\left(\partial_{t}^{2}+\omega_{\mathbf{q}}^{2}\right){\varphi}^{\phantom{\dagger}}_{-\mathbf{q}}/2+(\Lambda/2)\,|F(t)|^{2}\,\sqrt{N}\,{\varphi}^{\phantom{\dagger}}_{\mathbf{q}=0}$ is the quadratic part, including the external drive. It is most convenient to perform the calculations in the real time formalism in order to avoid tedious analytical continuation procedure required in the Matsubara formalism. Integrating out the electrons, we find:

[TABLE]

Here, $\hat{\mathcal{G}}_{0,\mathbf{k}}^{-1}(t,t^{\prime})=(i\partial_{t}-\xi_{\mathbf{k}}\,\hat{\sigma}_{3})\delta_{\mathcal{C}}(t,t^{\prime})$ and the trace implies momentum summation, contour time integration, and Nambu space summation. Expanding the second term in powers of $g$, we find:

[TABLE]

The first order correction $S_{1}[\varphi]$ vanishes for Holstein-type screened electron-phonon couplings (see the discussion after Eq. 35). The sum of higher order vertices can be diagrammatically represented as:

[TABLE]

It is convenient to make the forward/backward contour time indices explicit and perform a Keldysh rotation of $\varphi^{\pm}$ fields into symmetric (“classical”) and anti-symmetric (“quantum”) components, $\overline{\varphi}^{\alpha}=U_{\alpha\beta}\,\varphi_{\beta}$:

[TABLE]

The bare action in Keldysh representation reads:

[TABLE]

Likewise, the higher order terms in the Keldysh representation read:

[TABLE]

where:

[TABLE]

and $\Gamma_{\alpha\mu\nu}=\sum_{\beta}U_{\alpha\beta}U_{\mu\beta}U_{\nu\beta}\sigma^{z}_{\beta\beta}$ is a vertex in the Keldysh space. Finally, $\overline{\hat{\mathcal{G}}}_{0}$ is the bare propagator in the Keldysh space:

[TABLE]

Note that each matrix element additionally carries a $2\times 2$ Nambu structure. We restrict our analysis to the normal state hereafter, in which case, the Nambu structure is immaterial. The Nambu space traces reduce to a multiplicative factor of $2$ (= total spin degeneracy) for each electron loop. In the following sections, we briefly study the first few vertices in succession.

E.1 The second order correction: Landau damping

A direct calculation using Eq. (131) gives the matrix elements of $\Pi_{\mathbf{q}}^{\alpha_{1}\alpha_{2}}(t_{1},t_{2})\equiv-V_{\mathbf{q},-\mathbf{q}}^{\alpha_{1}\alpha_{2}}(t_{1},t_{2})$:

[TABLE]

where:

[TABLE]

The bare propagators in equilibrium are functions of $t_{1}-t_{2}$ and admit the following standard Fourier representation:

[TABLE]

Calculating $\Pi^{R/A/K}_{\mathbf{q}}$ in equilibrium is standard and yields the well-known Lindhard function Mahan (2013):

[TABLE]

Combining $S_{2}[\varphi]$ with the bare action $S_{0}[\varphi]$ yields the full quadratic part of the effective phonon-only action $S_{\mathrm{quad}}[\varphi]$:

[TABLE]

The quantity $\Gamma^{(\mathrm{ph})}_{\mathbf{q}}(\omega)\equiv-2\mathrm{Im}[\Pi_{\mathbf{q}}^{R}(\omega)]$ is of particular interest and represents the spectrum of the dissipative bath the electronic degrees of freedom provide for phonons. We observe that $\Gamma_{\mathbf{q}=0}(\omega)=0$ for finite $\omega$, which is an expected consequence of momentum conservation. Thus, the uniform lattice displacement $\langle\overline{\varphi}^{c}(t)\rangle$ experiences no friction from electrons. Employing the simplifications introduced in Sec. IV.2, i.e. flat EDOS, local approximation for $\Pi$, and $\mathbf{k},\mathbf{q}$-independent e-p coupling, we can calculate the local bath spectrum $\Gamma_{\ell}^{(\mathrm{ph})}$analytically at $T=0$:

[TABLE]

Given that the energy density of electrons is much lower than Einstein phonon energy scale in the present context, it is an excellent approximation to use the last equation even for nonequilibrium electronic states. The last result is akin to the well-known Allen’s formula Allen (1987) which is often used to infer electron-phonon coupling from the phonon linewidth broadening. Finally, one can calculate $\mathrm{Re}\,\Pi^{(\ell)}(\omega;T=0)$ within the same approximations to find:

[TABLE]

In the Migdal limit, $\omega\ll W_{\mathrm{el}}$ for all $\omega$ of interest. To leading order in $\omega/W_{\mathrm{el}}$, we find $\mathrm{Re}\,\Pi^{R}_{\ell}(\omega)\approx-4|g|^{2}\nu(0)^{2}\ln(2)\,W_{\mathrm{el}}$ which corresponds to a constant Lamb shift.

In summary, within the validity limits of the simplified model of Sec. IV.2, the dissipative effect of electrons on the dynamics of optical phonons can be modeled as a local quantum Ohmic bath:

[TABLE]

where $\omega_{L}\approx-4|g|^{2}\nu(0)^{2}\,\ln(2)\,W$ is the Lamb shift and $\gamma_{\ell}=2\pi\,\omega_{0}\,\nu(0)\,\lambda$ is the dimensionless friction constant expressed in terms of the mass enhancement factor $\lambda$ (see Eq. 55). Furthermore, the above expressions remain valid as long as the electrons approximately remain in a quantum degenerate state.

E.2 The third-order correction: electron-mediated cubic nonlinearity

The cubic vertex $V^{\alpha_{1},\alpha_{2},\alpha_{3}}_{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}}$ has a complicated spatial and temporal structure due the nonlocality of electrons. Here, we rather focus on calculating the retarded phonon self-energy correction that arises from this cubic vertex rather than a general analysis. Recalling that the lattice has a large coherent uniform displacement in our problem, we find that the leading self-energy correction is obtained by contracting one of the legs (the third leg without the loss of generality) with the classical displacement $\varphi_{0}$. The resulting self-energy correction, $\Pi_{\Delta,\mathbf{q}}(t_{1},t_{2})$, has the following diagrammatic representation:

[TABLE]

Integration over $\tau$, the time argument of $\varphi_{0}(\tau)$, is implied. Contracting the third leg with $\varphi_{0}$ sets $\alpha_{3}$ to 1, i.e. to the “classical” Keldysh index. The retarded phonon self-energy is obtained by further choosing $\alpha_{1}=2$ and $\alpha_{2}=1$ (e.g. see Eq. 133). Performing the intermediate Keldysh space traces in Eq. (131), we find:

[TABLE]

Notice that the electron lines, starting form the $\mathbf{k}+\mathbf{q}$ line and traversing counter clockwise, assume the following Keldysh space labels: KAA, RKA, RRK. It can be shown that the same structure applies to higher order single-electron-loop vertices: with $N$ fermion propagators in a loop and $N-2$ classical field contractions, the retarded self-energy comprises $N$ terms, and the electron propagators in each term have Keldysh space labels $[R\ldots R]K[A\ldots A]$ in a counter clockwise fashion. The index subsets $[R\ldots R]$ and $[A\ldots A]$ comprise $N-1$ indices, and either subset can be empty (for example, see the next section for the quartic vertex).

The overall symmetry factor can be worked out as follows: $1/3!$ from the definition of $S_{3}[\varphi]$, 3 choices for the classical leg, 2 choices for attaching one of the two remaining legs to the left external point, and a factor of $i^{2}$ from the two phonon propagators, amounting to $2\times 3\times i^{2}/3!=-1$. Thus, we obtain:

[TABLE]

Note that ${\varphi}^{\phantom{\dagger}}_{\mathbf{q}=0}(\tau)=(1/\sqrt{N})\sum_{j}\varphi_{j}(\tau)=\sqrt{N}\varphi_{0}(\tau)$, where $\varphi_{j}(\tau)=\varphi_{0}(\tau)$ is the uniform ionic displacement at site $j$. Also, the factor $2/\sqrt{2}$ arises from the definition of the “classical” component, i.e. $\overline{\varphi}^{c}=(\varphi^{+}_{0}+\varphi^{-}_{0})\sqrt{2}=2\varphi_{0}/\sqrt{2}$. To proceed, we assume $\varphi_{0}(\tau)=A\,e^{i\Omega\tau}$ and a thermal state for electrons. Taking a Wigner transform $(t_{1},t_{2})\rightarrow(\omega,t)$, we find:

[TABLE]

Performing the time integrals and subsequently the frequency integrals over $\omega_{1}$ and $\omega_{2}$ is a lengthy calculation and we quote the final result:

[TABLE]

We recall that the uniform displacement $\varphi_{0}(\tau)$ approximately takes the form $\varphi_{0}(\tau)\approx\varphi_{0}+\varphi_{1}\,\cos(2\Omega_{\mathrm{drv}}t)=\varphi_{0}+(\varphi_{1}/2)\,e^{2i\Omega_{\mathrm{drv}}t}+(\varphi_{1}/2)\,e^{-2i\Omega_{\mathrm{drv}}t}$ for a slowly ramped-up drive (see Eq. 5). Accordingly, the complete cubic self-energy correction is the sum of three terms obtained from replacing $(A,\Omega)\rightarrow(\varphi_{0},0)$, $(\varphi_{1}/2,+2\Omega_{\mathrm{drv}})$ and $(\varphi_{1}/2,-2\Omega_{\mathrm{drv}})$ in Eq. (144).

Since $\Pi_{\Delta,\mathbf{q}}^{R}(\omega,t)\propto g_{\mathbf{k},\mathbf{k}}$, it vanishes for screened Holstein-type electron-phonon couplings. Had $g_{\mathbf{k},\mathbf{k}}$ been finite, however, $\Pi_{\Delta,\mathbf{q}}^{R}(\omega,t)$ would show a divergent behavior for $\mathbf{q}\approx 0$ and $\omega\simeq\Omega_{\mathrm{drv}}\simeq\omega_{\mathbf{q}}$. In any event, $\Pi_{\Delta,\mathbf{q}}^{R}(\omega,t)$ remains $\mathcal{O}(g^{3})$ and non-divergent for $\Omega_{\mathrm{drv}}\sim\omega_{\mathbf{q}}/2$. Therefore, the electronic contribution to cubic lattice nonlinearity and its corresponding phonon self-energy corrections are negligible. We show in the next section that the situation is very different for the forth order correction.

E.3 The forth-order correction: electron-mediated quartic nonlinearity

We can similarly calculate the contribution of the quartic vertex (the last diagram in Eq. 127) to the phonon self-energy. In this case, two diagrams with different topologies comprise the leading order contribution to the quartic self-energy correction:

[TABLE]

We attach $A_{1}e^{i\Omega_{1}\tau_{1}}$ and $A_{2}e^{i\Omega_{1}\tau_{2}}$ to two of the external legs. This can be done in $3\times 2=6$ and $2\times 2=4$ different ways for (a) and (b) topologies, respectively, this:

[TABLE]

Performing the intermediate Keldysh space summations in Eq. (131), we find:

[TABLE]

where $\mathcal{I}_{4}=\{KAAA,RKAA,RRKA,RRRK\}$ denotes the set of Keldysh space labels of the four electron propagators. After a lengthy but straightforward calculation, we find the contribution of the first diagram to be:

[TABLE]

Similarly, the contribution of the second diagram is found as:

[TABLE]

The $\omega_{4}$-integration is easily performed since for each choice of Keldysh space labels $(a_{1}\ldots a_{4})\in\mathcal{I}_{4}$, one of the electrons is on shell (Keldysh) and fixes the value of $\omega_{4}$. Assuming a coherent displacement like $\varphi_{0}(\tau)=\varphi_{0}+\varphi_{1}\,\cos(2\Omega_{\mathrm{drv}}\tau)$, the resulting self-energy contributions will have 3 contributions: a constant (dc) contribution, a contribution $\propto\cos(2\Omega_{\mathrm{drv}}t)$, and a contribution $\propto\cos(4\Omega_{\mathrm{drv}}t)$. Each contribution can be found by making appropriate choices for $(A_{1},\Omega_{1})$ and $(A_{2},\Omega_{2})$. The most interesting contribution is the one $\propto\cos(4\Omega_{\mathrm{drv}}t)$ which is found by substituting $(A_{1},A_{2};\Omega_{1},\Omega_{2})\rightarrow(\varphi_{1}/2,\varphi_{1}/2;\pm 2\Omega_{\mathrm{drv}},\pm\Omega_{\mathrm{drv}})$ in Eqs. (147)-(148) and summing up the four contributions. We quote the final result from this lengthy calculation:

[TABLE]

where:

[TABLE]

As in the cubic vertex case, this contribution also vanishes for a perfectly screened Holstein-type electron-phonon coupling since $\Pi^{R,\mathrm{ac}}_{\square,\mathbf{q}}(\omega,t)\propto|g_{\mathbf{k},\mathbf{k}}|^{2}$. In a more realistic model, $g_{\mathbf{k},\mathbf{k}}$ is generically non-vanishing.

The $\mathbf{k}$-integral in Eq. (149) can be calculated in the limit $\mathbf{q}\approx 0$ and assuming a constant electronic density of states and zero temperature. In the vicinity of the parametric resonance $\Omega_{\mathrm{drv}}\sim\omega_{\mathbf{q}}/2$, we find:

[TABLE]

The logarithmic divergence could be anticipated from $D_{\mathbf{k},\mathbf{q}}\propto(\omega_{\mathbf{q}}^{2}-4\Omega_{\mathrm{drv}}^{2})$ in the limit $\omega=\omega_{\mathbf{q}}$ and $\mathbf{q}\approx 0$.

The above finding has a consequential implication: the electronic contribution to the lattice nonlinearity, even though is $\sim\mathcal{O}(g^{2})$ and small in general, in the presence of coherent lattice oscillations leads to a self-energy correction that diverges logarithmically in the vicinity of $\Omega_{\mathrm{drv}}\sim\omega_{\mathbf{q}}/2$. Thus, even if purely ionic contributions to the lattice nonlinearity is small, large nonlinearities will be dynamically generated as a matter of coupling to electrons. Also, note that $\Pi^{R,\mathrm{ac}}_{\square,\mathbf{q}\approx 0}(\omega_{\mathbf{q}},t)\propto\cos(2\omega_{\mathbf{q}}t)$ which is precisely the COM time-dependence required for giving rise to parametric amplification of the lattice response as discussed in Sec. III.

Appendix F Summary of numerical methods

In this appendix, we provide a summary of numerical methods for solving the quantum Floquet-Boltzmann kinetic equations for the lattice and electronic degrees of freedom. In reality, the two systems are coupled and must propagated forward in time self-consistently. The perturbative framework adopted in this work (when physically permissible), allows us to study the two systems in iterations: the dynamics of the lattice is worked out assuming unperturbed equilibrium electron propagators, the nonequilibrium correction to electron propagators electrons are worked out on the backdrop of the driven lattice, and so on; see Fig. 7. This iterative procedure is expected to converge to the self-consistent solution of the fully coupled system in the weak-coupling limit.

F.1 Solving the Quantum Floquet-Boltzmann kinetic equation for lattice displacement and phonon propagators

The quantum Floquet-Boltzmann kinetic equations for the lattice displacement and phonon propagators were worked out in Sec. IV.4. The final result is the coupled system of equations given in Eq. (44), Eqs. (45)-(45), and Eqs. (46)-(46). Coupling to electrons only appear in the bath term, $\Pi^{(\ell);R/A/K}_{n,m}(\omega;t)$, which we assume is given to us in this section. We take a further simplifying step and neglect the COM time dependence of $\Pi^{(\ell)}$ which in indeed the case if the bath is approximately calculated using equilibrium electron propagators (see Sec. E.1). The following analysis can be easily generalized for time-dependent baths, e.g. as required for the next iterations if one were to follow the perturbative decoupling recipe mentioned above.

The major difficulty in time-stepping Eqs. (45)-(45) and Eqs. (46)-(46) using ODE solvers is three-fold:

Time-derivatives appears on both sides of equations, and $\partial_{t}$ of different Floquet components $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\omega;t)$ are coupled due to the lattice nonlinearity and the bath. In other words, $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\omega;t)$ is only implicitly given by Eqs. (45)-(45) and Eqs. (46)-(46).

• 2.

We have two sets of evolution equations for the retarded and Keldysh propagators: one obtained from the forward KB equation (Eqs. 45 and 46), and another from the backward KB equation (Eqs. 45 and 46). In general, these two are complementary. For instance, the direct numerical solution of two-time propagators requires the forward and backward equations to step the propagators forward in the first and second times, respectively, e.g. see Ref. Semkat et al. (1999). In the kinetic approximation, however, only the COM time is stepped forward while the relative time is transformed to the frequency domain and is carried as a label. In theory, one may choose to work with either of the forward or backward equations for time stepping as both are correct to $\mathcal{O}(\partial_{t})$. However, the mixing of large non-gradient terms and small gradient terms leads to undesirable numerical errors.

• 3.

We found the system to be marginally stiff, requiring a robust ODE solver with adaptive time-stepping and local error control. This leads to unavoidably long run times. The application of stiff solvers is challenging as the Jacobian of the system is dense and is difficult to calculate.

Let us note that we do not need to calculate $\mathcal{D}^{A}$ as a separate quantity since the identity $\mathcal{D}^{R}(t_{1},t_{2})=\mathcal{D}^{A}(t_{2},t_{1})$ implies $\mathcal{D}^{A}_{n}(\omega;t)=\mathcal{D}^{R}_{n}(-\omega;t)$. Furthermore, the exact identities $\mathcal{D}^{R}_{-n}(\omega)=[\mathcal{D}^{R}_{n}(-\omega)]^{*}$, $i\mathcal{D}^{K}_{-n}(\omega)=[i\mathcal{D}^{K}_{n}(-\omega)]^{*}$, and $i\mathcal{D}^{K}_{n}(-\omega)=i\mathcal{D}^{K}_{n}(\omega)$ allow us to restrict the numerical calculation to non-negative Floquet indices.

Setting up the linear system and calculating the explicit $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\omega;t)$— The second issue mentioned above can be circumvented using anti-symmetric and symmetric combinations of Eqs. (45)-(45) and Eqs. (46)-(46), respectively. The issue of implicitness, however, remains challenging. In particular, the $\omega$-integral appearing in $\chi_{n}(t)$ and the appearance of $\partial_{t}U_{n}(t)$ in the kinetic equations implies that neither Floquet indinces, nor $\omega$ are “good” numbers. In other words, the kinetic equations of the lattice displacement and phonon propagators pose a dense linear system for $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\omega;t)$ and $\partial_{t}\varphi_{n}(t)$. In order to find $\partial_{t}\varphi_{n}(t)$ and $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\omega;t)$ explicitly, at each time $t$, we carefully index $\partial_{t}\mathcal{D}^{R/K}_{n,m}(\pm\omega;t)$ and $\partial_{t}\varphi_{n}(t)$ for all $(\omega,n,m)$, cast the coupled kinetic equations into a linear system and solve it via LU decomposition.

We perform the calculations on a regular frequency grid $\omega\in[-\omega_{M},\omega_{M}]$ where $\omega_{M}$ is a high frequency cutoff. The grid spacing is chosen as rational fraction of $\Omega/2$ close to $0.1\gamma_{\ell}$ in order to ensure that $\omega\pm n\Omega/2$ belongs to the grid. This allows us to identify a large fraction of unknown time-derivatives and matrix elements with one another and greatly reduce the dimension of the linear system. The $\omega$-derivatives are calculated using the 5-point finite difference approximation, and the $\omega$-integral appearing in Eq. (44c) is approximated using the trapezoid rule. We choose the Floquet cutoff $n_{D}=2$, and the frequency cutoff $\omega_{M}=2\omega_{0}+5\gamma_{\ell}+(n_{D}+1)\Omega$. This choice ensures that all involved propagators remain small and negligible for $|\omega|>\omega_{M}$. We carefully checked that increasing $\omega_{M}$ and $n_{D}$ had a negligible and controllably small effect on the results. For an $\omega$-grid with $\sim 500$ points, one needs to solve a linear system of size $\sim 6000\times 6000$ for each calculation of the explicit time derivatives.

Initial thermal state and renormalized phonon frequency— The lattice is in a thermal equilibrium state at the bath temperature before the drive ramped up. To find conditions describing the equilibrium state, we set the external drive and time derivatives to zero in the described evolution equations, assume $\varphi_{n}(t)\rightarrow\varphi_{0}\,\delta_{n,0}$, $U_{n}(t)\rightarrow U_{0}\,\delta_{n,0}$, $\chi_{n}(t)\rightarrow\chi_{0}\,\delta_{n,0}$, and $\mathcal{D}_{n}^{R/A/K}(\omega)\rightarrow\delta_{n,0}\,\mathcal{D}^{R/A/K}_{0}(\omega)$. This leads to the following set of coupled equations:

[TABLE]

The first and last two equations must be solved self-consistently, leading to a renormalized phonon frequency:

[TABLE]

Numerical time stepping— Provided that $\{\varphi_{n}(t)\}$, $\{U_{n}(t)\}$, $\{\chi_{n}(t)\}$, and $\{i\mathcal{D}^{R,K}_{n}(\omega;t)\}$ are known for all $\omega$ on a regular grid, we obtain the explicit time-derivatives of these quantities using by solving a linear system as described earlier. We can then invoke an explicit ODE solver to perform time-stepping. Here, we integrated the ODE using the adaptive Runge-Kutta-Fehlberg(4,5) method with local relative error tolerance of $10^{-6}$.

F.2 Solving the Quantum Floquet-Boltzmann kinetic equation for electrons

The quantum Floquet-Boltzmann kinetic equation for $\{\psi_{n}(\omega;t)\}$ (see Eq. 58) was derived in Sec. IV.5.3. Similar to the kinetic equation for the phonons, this kinetic equation is also a formidably dense implicit integral equation for $\{\partial_{t}\psi_{n}(\omega;t)\}$ in which all frequencies and Floquet indices are coupled and defies the immediate application of an explicit ODE solver. In this section, we describe a numerical strategy for solving this equation.

Preliminaries— As a first step, we use the exact identities $\psi^{*}_{n,m}(\omega)=\psi_{-n,m}(\omega)$ and $\Sigma^{A}_{n,m}(\omega)=[\Sigma^{R}_{-n,m}(\omega)]^{*}$ to cast Eq. (IV.5.3) into a more useful form:

[TABLE]

We have dropped the common $(\omega;t)$ argument from all quantities for brevity. Summation over repeated indices is implied everywhere in this section. The numerical integration of this equation is complicated by the fact that the self-energies are functionals of $\psi$, so that $\partial_{t}\Sigma^{R}$ terms implicitly involve $\partial_{t}\psi$. This functional dependence can be made explicit using Eqs. (70)-(70b):

[TABLE]

where:

[TABLE]

The following useful identities can be established using the symmetries of $\psi$, $F^{K/\rho}$, and the properties of Kramers-Kronig transforms:

[TABLE]

As a result, we only need to calculate each quantity only for $\omega>0$. Also, save for $\Sigma^{R}_{n}$, all other quantities can be calculated for $n\geq 0$.

The frequency grid— We proceed by generating a grid $X_{\omega}$ in the interval $[0,\omega_{c}]$. Here, $\omega_{c}$ is an appropriate cutoff $\omega_{c}\gg 1/\beta,\omega_{0},\Omega$. We generate the grid $X_{\omega}$ such that for all $\omega\in X_{\omega}$, if $\omega_{m}\equiv\omega+m\Omega/2<\omega_{c}$, then $\omega_{m}\in X_{\omega}$. We call such a grid $X_{\omega}$ as a Floquet-closed grid. We will shortly see that a Floquet-closed grid leads to a significant reduction in computational complexity by allowing us to reuse previously calculated integrals. In practice, it is necessary to generate a non-uniform grid that emphasizes on the $\omega\lesssim 1/\beta$ region. To this end, we create two uniform grids, $X^{\mathrm{th}}_{\omega}\in[0,c/\beta]$, and $X^{\mathrm{>}}_{\omega}\in[c/\beta,\omega_{c}]$, and concatenate them. Crucially, we choose the grid spacings $\delta\omega^{\mathrm{th}}$ and $\delta\omega^{>}$ such that both are integer multiples of $\Omega/(2N)$ for some $N$. Once we have this basic two-scale grid, we pool together $|X^{\mathrm{>}}_{\omega}\cup X^{\mathrm{th}}_{\omega}+m\Omega/2|$ for $|m|<m_{c}$ and keep the unique points to find $X_{\omega}$.

Calculating the required matrix elements— We assume that $\{\psi_{n}(\omega;t)\}$ are known for $\omega\in X_{\omega}$ and $0\leq n\leq N_{\psi}$ for some cutoff $N_{\psi}\geq N_{D}$. It is trivial to calculate $i\Sigma^{K}$ numerically based on Eq. (155) using a quadrature formula. To find $\Sigma^{R}$ and $\partial_{t}\Sigma^{R}$, we first calculate $\mathsf{F}[i\mathcal{D}^{K}]_{n}$ and $\partial_{t}\,\mathsf{F}[i\mathcal{D}^{K}]_{n}$, both of which are trivial. To calculate the contribution from $\{\rho_{n}\}$, we calculate the following quantities:

[TABLE]

for $\omega_{j},\omega_{k}\in X_{\omega}$. The $\tilde{\mathsf{K}}$ terms stem from $\int_{\omega_{c}}^{\infty}\mathrm{d}\omega^{\prime}$ assuming $\psi_{n}(\omega^{\prime})\approx\delta_{n,0}$ for $\omega^{\prime}>\omega_{c}$. The proper $\omega^{\prime}$-integrals must be approximated with quadratures much finer that $X_{\omega}$ grid spacing, and this is necessary since the integrands can vary on shorter scales than the grid spacing of a practically sized $X_{\omega}$. The improper $\omega^{\prime}$-integrals can be calculated using Möbius transformation and then using standard proper quadratures. Calculating $\mathsf{K}_{n,jk}[\rho](t)$ and $\mathsf{K}_{n,jk}[\partial_{t}\rho](t)$ is quite expensive as the integrand is given by a Kramers-Kronig integral and must be obtained numerically for every integration point. Having calculated these quantities, we may compose the full expression for $\Sigma^{R}_{n}(\omega_{j};t)$ and $\partial_{t}\Sigma^{R}_{n}(\omega_{j};t)$ approximately as:

[TABLE]

Here, $N_{\omega}$ is the number of grid points in $X_{\omega}$. We have also used linear interpolation for the values of $\psi$ between consecutive grid points. If $\omega_{c}\gg 1/\beta$, the high energy tail of Fermi distribution indeed remains intact (i.e. we assume $\omega_{c}$ is large enough so that no particles will be excited to energies above $\omega_{c}$). The important point about using the Floquet closed $X_{\omega}$ is that once we calculate $\mathsf{K}_{n}(\omega_{j},\omega_{k})$ for $\omega_{j}\in X_{\omega}$, we immediately get $\mathsf{K}_{n}(\omega_{j}+m\Omega/2,\omega_{k})$ for all $m$ using a combination of shifts, inversions, and conjugation (see Eq. 157). In other words, we do not need to perform the expensive calculation of $\mathsf{K}$ for all $m$-shifted $\omega$ points.

Setting up the linear system and time-stepping— Plugging Eq. (F.2) expression into Eq. (154), we find an explicit linear system for $\partial_{t}\psi_{n}(\omega_{j};t)$ for $\omega_{j}\in X_{\omega}$. The $\omega$-derivatives are found using 5-point finite difference approximation on the $X_{\omega}$ grid. This system can be mapped to a matrix equation by (1) indexing $\{\mathrm{Re}\,\psi_{n}(\omega_{j};t),\mathrm{Im}\,\psi_{n}(\omega_{j};t)\}$ for all $n$ and $\omega_{j}\in X_{\omega}$, and (2) setting up a mapping from $\{\mathrm{Re}\,\psi_{n}(\omega_{j}+m\Omega/2;t),\mathrm{Im}\,\psi_{n}(\omega_{j}+m\Omega/2;t)\}$ to the corresponding indexed values for all $m$ using the symmetries of $\psi$. Having a recipe to calculate $\partial_{t}\psi_{n}(\omega;t)$, we proceed and integrate the ODE using the adaptive Runge-Kutta-Fehlberg(4,5) method with local relative error tolerance of $10^{-6}$.

F.3 Numerical analysis of the spectrum of Floquet-Migdal-Eliashberg gap functional

Calculating the spectrum of the Floquet-Migdal-Eliashberg (FME) gap functional, given in Eq. (85), involves three steps: (1) calculating $\Sigma^{R/K}_{n,m}(\omega;t)$ in the normal state, (2) calculating the Floquet matrix elements of the anomalous response $\mathsf{Q}^{n,m}_{n^{\prime},m^{\prime}}(\omega;t)$, and (3) calculating the spectrum of the FME gap functional.

The first step is identical to the procedure described in Sec. F.2. In the equilibrium-electron approximation, we use $\psi_{n,m}(\omega;t)\rightarrow\delta_{n,0}\,\tanh[\beta(\omega-m\Omega/2)/2]$ in calculating the self-energies rather than using $\psi_{n,m}(\omega;t)$ found from solving the Floquet-Boltzmann equation for electrons.

The second step involves inverting the coefficient matrix of $\delta\mathcal{F}^{R}_{n,m}(\omega;t)$ which can be read from Eq. (81a). To this end, we truncate the intermediate $n^{\prime}$ Floquet band index summation to $|n^{\prime}|\leq N_{\phi}$ and Floquet quasi-momentum indices to $|m|\leq N_{m}$. The truncated system of equations is then carefully mapped to a proper linear system, $\sum_{n^{\prime}=-N_{\phi}}^{+N_{\phi}}\sum_{m^{\prime}=-N_{m}}^{+N_{m}}\mathcal{C}^{n,m}_{n^{\prime},m^{\prime}}(\omega;t)\,\delta\mathcal{F}^{R}_{n,m}(\omega;t)=-2\pi i\,\phi_{n,m}(\omega;t)$. The Floquet matrix elements of the anomalous response is readily found by inverting $\mathcal{C}^{n,m}_{n^{\prime},m^{\prime}}$ in the space of paired Floquet indices $(n,m)$:

[TABLE]

In practice, we found the final results to be accurate to $10^{-4}$ with the choice $N_{\phi}=N_{m}=N_{D}+2$ where $N_{\Sigma}$ is the previously chosen Floquet band cutoff in calculating the retarded self-energy.

The last step is slightly more involved. The overall strategy is to formally interpolate $\Delta_{n}(\omega;t)$ over a finite grid $G_{\omega}$, plug the interpolation formula in Eq. (85), read off the coefficients of $\Delta_{n}(\omega\in G_{\omega};t)$, and calculate its spectrum. Even though a brute-force discretization is equally applicable in principle, the uniform grid must be very dense in order to obtain accurate results, leading to calculating the spectrum of intractably large matrices. The interpolation procedure allows us to obtain accurate results using much coarser grids.

Setting up the grid— We generate $G_{\omega}$ by concatenating three grids $G_{\omega}=G^{(1)}_{\omega}\cup G^{(2)}_{\omega}\cup G^{(3)}_{\omega}$ where $G^{(1)}_{\omega}$ is a uniform grid for $\omega\in[0,10/\beta)$ where $\beta^{-1}\sim 0.05\,\Omega_{0}$ is the typical effective temperature of electrons, $G^{(2)}_{\omega}$ in another uniform grid for $\omega\in[10/\beta,\omega_{c})$ where $\omega_{c}\sim 10\,\Omega_{0}$ is a typical scale beyond which variations of $\Delta_{n}(\omega;t)$ becomes negligible, and finally, $G^{(3)}_{\omega}$ is a log-scaled grid for $\omega\in[\omega_{c},\infty)$. In practice, we found allocating 100 points for each sub-grid produced results accurate to $10^{-4}$.

Setting up the coefficient matrix— We approximate $\Delta_{n}(\omega;t)$ over $G_{\omega}$ using a linear interpolant:

[TABLE]

where $j_{\omega}$ is the nearest grid point to the left of $\omega$. Plugging this ansatz into Eq. (85), we find:

[TABLE]

where:

[TABLE]

The end points, $k=|\mathcal{G}_{\omega}|$ and $k=1$, only get contributions from the first and second integrals, respectively. The $\omega^{\prime}$-integrals are performed via an adaptive Gauss-Kronrod quadrature and is refined until a tolerance of $10^{-6}$ is achieved. The kernel $K_{n,m}(\omega,\omega^{\prime};t)$ is given by Eq. (84) and each evaluation requires performing a numerical Kramers-Kronig transform. Calculating the matrix elements $\mathsf{K}^{\Delta}_{n,m}(\omega_{j},\omega_{k};t)$ is the most computationally expensive part of this section. Finally, we decompose $\Delta_{n}(\omega;t)$ into real and imaginary parts and use the relations $\mathrm{Re}[\Delta_{n}(\omega;t)]=\mathrm{Re}[\Delta_{-n}(\omega;t)]$, $\mathrm{Im}[\Delta_{n}(\omega;t)]=-\mathrm{Im}[\Delta_{-n}(\omega;t)]$ to cast Eq. (167) as a matrix equation. The coefficient matrix acts on space of bundled labels $(n,j,o)$ where $n$ is the Floquet index, $j$ is the grid point index, and $o=0,1$ indicates real and imaginary component of $\Delta_{n}$.

Finally, we impose cutoffs $N_{\Delta}$ and $N_{m}$ over the Floquet band index of $\Delta$ and the internal $m^{\prime}$ quasi-momentum summation. The Floquet cutoff for $K$ is $N_{K}=N_{D}+N_{\psi}$ where $N_{D}$ and $N_{\psi}$ are the previously chosen Floquet cutoffs for phonon propagators and electron energy statistics, respectively. We found $N_{\Delta}=N_{m}=N_{K}+2$ to produce results accurate to $10^{-4}$. Assuming $N_{\psi}=N_{D}=2$ and $|G_{\omega}|=300$, the final coefficient matrix is has a dimension $3900\times 3900$ and its spectrum can be easily found numerically.