Nonequilibrium steady states and transient dynamics of conventional superconductors under phonon driving

TL;DR

This study investigates how phonon driving affects superconductivity in the Holstein model, revealing that despite potential enhancement of interactions, external driving generally suppresses superconducting order due to heating effects, especially at parametric resonance.

Contribution

It provides a comprehensive nonequilibrium analysis of phonon-driven superconductors using dynamical mean-field theory, highlighting the suppression of superconductivity under external modulation.

Findings

Parametric phonon driving can enhance effective attraction via Floquet side bands.

Superconducting order is suppressed or melted under phonon driving.

Heating effects dominate transient dynamics, weakening superconductivity.

Abstract

We perform a systematic analysis of the influence of phonon driving on the superconducting Holstein model coupled to heat baths by studying both the transient dynamics and the nonequilibrium steady state (NESS) in the weak and strong electron-phonon coupling regimes. Our study is based on the nonequilibrium dynamical mean-field theory, and for the NESS we present a Floquet formulation adapted to electron-phonon systems. The analysis of the phonon propagator suggests that the effective attractive interaction can be strongly enhanced in a parametric resonant regime because of the Floquet side bands of phonons. While this may be expected to enhance the superconductivity (SC), our fully self-consistent calculations, which include the effects of heating and nonthermal distributions, show that the parametric phonon driving generically results in a suppression or complete melting of the SC order. In the strong coupling regime, the NESS always shows a suppression of the SC gap, the SC order parameter and the superfluid density as a result of the driving, and this tendency is most prominent at the parametric resonance. Using the real-time nonequilibrium DMFT formalism, we also study the dynamics towards the NESS, which shows that the heating effect dominates the transient dynamics, and SC is weakened by the external modulations, in particular at the parametric resonance. In the weak coupling regime, we find that the SC fluctuations above the transition temperature are generally weakened under the driving. The strongest suppression occurs again around the parametric resonances because of the efficient energy absorption.

Equations

$$H_{\mathrm{tot}}(t)=H_{\mathrm{sys}}(t)+H_{\mathrm{mix}}+H_{\mathrm{bath}},$$

$$H_{\mathrm{sys},0}=$$

$$+{\omega }_0\sum _i{a}_i^{†}{a}_i+g\sum _i(a_i^{†}+a_i)(n_i-1),$$

$$H_{\mathrm{ph}}(t)={\omega }_{0P}(t)\sum _i\frac{P_i^2}{4}+{\omega }_{0X}(t)\sum _i\frac{X_i^2}{4}.$$

$$H_{\mathrm{ph}}(t)={\omega }_0\sum _i\frac{P_i^2}{4}+{\omega }_{0X}(t)\sum _i\frac{X_i^2}{4}$$

$$\text{with }{\omega }_{0X}(t)={\omega }_0+\Delta {\omega }_0\cos (\Omega t),$$

$$H_{\mathrm{ph}}(t)={\omega }_0(t)\sum _i\left(\frac{P_i^2}{4}+\frac{X_i^2}{4}\right)$$

$$\text{with }{\omega }_0(t)={\omega }_0+\Delta {\omega }_0\cos (\Omega t).$$

$$G(t,t^{\prime })$$

$$D(t,t^{\prime })$$

$$\check{G}(t,t^{\prime })\equiv \left[\begin{array}{cc}{G}^{11}(t,t^{\prime })& G^{12}(t,t^{\prime })\\ G^{21}(t,t^{\prime })& G^{22}(t,t^{\prime })\end{array}\right].$$

$$\underline{G}$$

$${\int }_{-\infty }^{\infty }d{t}_1[{\underline{G}}_0^{-1}(t,t_1)-\underline{\Sigma }(t,t_1)]\cdot \underline{G}(t_1,t^{\prime })=\underline{I}\delta (t-t^{\prime }),$$

$${\int }_{-\infty }^{\infty }d{t}_1[{\underline{D}}_0^{-1}(t,t_1)-\underline{\Pi }(t,t_1)]\cdot \underline{D}(t_1,t^{\prime })=\underline{I}\delta (t-t^{\prime }).$$

$$\Gamma (\omega )\equiv -\mathrm{Im}{\Sigma }_{\mathrm{bath}}^R(\omega ).$$

$${\Sigma }_{\mathrm{tot}}={\Sigma }_{\mathrm{bath}}+{\Sigma }_{\mathrm{int}}.$$

$$\Gamma (\omega )={\gamma }_{\text{el}}\sqrt{1-{\left(\frac{\omega }{W_{\mathrm{ebath}}}\right)}^2},$$

$$H_{\mathrm{mix},\mathrm{ph}}$$

$$H_{\mathrm{bath},\mathrm{ph}}$$

$${\Pi }_{i,\mathrm{bath}}(t,t^{\prime })$$

$$D_{0,p}(t,t^{\prime })$$

$$-i[{\theta }_c(t^{\prime },t)+f_b({\omega }_p)]e^{-i(t^{\prime }-t){\omega }_p},$$

$${\Pi }_{\mathrm{tot}}={\Pi }_{\mathrm{bath}}+{\Pi }_{\mathrm{int}}.$$

$${\Pi }_{\mathrm{bath}}^R(\omega )$$

$$B_{\mathrm{bath}}(\omega )$$

$$=\pi \sum _p{\mathcal{V}}_p^2[\delta (\omega -{\omega }_p)-\delta (\omega +{\omega }_p)].$$

$${\Pi }_{\mathrm{bath}}^R(\omega )$$

$${\Pi }_{\mathrm{bath}}^A(\omega )$$

$${\Pi }_{\mathrm{bath}}^K(\omega )$$

$$\displaystyle=-2i\coth\Bigl{(}\frac{\beta\omega}{2}\Bigl{)}B_{\rm bath}(\omega).$$

$$B_{\mathrm{bath}}(\omega )={\mathcal{V}}^2\left\{\frac{\gamma }{(\omega -{\omega }_D{)}^2+{\gamma }^2}-\frac{\gamma }{(\omega +{\omega }_D{)}^2+{\gamma }^2}\right\}.$$

$$\begin{array}{rl}{G}^{\mathrm{R},\mathrm{A},\mathrm{K}}(t_{\mathrm{r}};t_{\mathrm{av}})=& G^{\mathrm{R},\mathrm{A},\mathrm{K}}(t_{\mathrm{r}};t_{\mathrm{av}}+\mathcal{T}),\\ D^{\mathrm{R},\mathrm{A},\mathrm{K}}(t_{\mathrm{r}};t_{\mathrm{av}})=& D^{\mathrm{R},\mathrm{A},\mathrm{K}}(t_{\mathrm{r}};t_{\mathrm{av}}+\mathcal{T}),\end{array}$$

$${\mathbf{F}}_{mn}(\omega )\equiv$$

$$\frac{1}{\mathcal{T}}{\int }_0^{\mathcal{T}}d{t}_{\mathrm{av}}{\int }_{-\infty }^{\infty }d{t}_{\mathrm{r}}{e}^{i(\omega +\frac{m+n}{2}\Omega )t_{\mathrm{r}}}{e}^{i(m-n)\Omega t_{\mathrm{av}}}F(t_{\mathrm{r}};t_{\mathrm{av}}).$$

$$[{\underline{\mathbf{G}}}_0^{-1}(\omega )-\underline{\mathbf{\Sigma }}(\omega )]\cdot \underline{\mathbf{G}}(\omega )=\underline{\mathbf{I}},$$

$$[{\underline{\mathbf{D}}}_0^{-1}(\omega )-\underline{\mathbf{\Pi }}(\omega )]\cdot \underline{\mathbf{D}}(\omega )=\underline{\mathbf{I}}.$$

$$H_{\mathrm{ph}}(t)={\omega }_{0P}(t)\frac{P^2}{4}+{\omega }_{0X}(t)\frac{X^2}{4}.$$

Full text

I Introduction

The prospect of nonequilibrium exploration and control of material properties has caught the interest and attention of broad segments of the condensed matter community. Giannetti2016 In particular, the resonant excitation of mid-infrared phonon modes has opened a new pathway for manipulating properties such as metallicity,Cavalleri2007 magnetismForst2011b and superconductivity (SC), Cavalleri2011 ; Kaiser2014 ; Hu2014 ; Mankowsky2014 ; Mitrano2016 ; Cantaluppi2017 and for modifying lattice structures in ways which may not be achievable in equilibrium. Especially, the light-induced SC-like behavior in the optical conductivity for high $T_{c}$ cuprates and doped fullerides Cavalleri2011 ; Kaiser2014 ; Hu2014 ; Mankowsky2014 ; Mitrano2016 has stimulated several theoretical investigations.Sentef2015 ; Knap2016 ; Komnik2016 ; Kennes2017 ; Kim2016 ; Okamoto2016 ; Sentef2017 ; Babadi2017 ; Mazza2017 ; Fernandes2017 So far, two possible scenarios for the enhancement of SC have been discussed in the literature based on the enhanced pairing interaction out of equilibrium: (i) Under time-periodic driving, the time-averaged Hamiltonian can be modified in such a way that superconductivity is favored, Sentef2015 ; Kim2016 ; Mazza2017 and (ii) the interaction can be effectively enhanced in the driven state by dynamical effects. Knap2016 ; Komnik2016 ; Kennes2017 ; Sentef2017 Both effects can be induced by the resonant excitation of mid-infrared phonon modes, and previous theoretical investigations have tried to clarify their influence on SC.

Concerning the effect of time-periodic driving on the time-averaged Hamiltonian, it has been pointed out that the resonantly excited mid-infrared phonon can lead to favorable parameter changes,Sentef2015 ; Kim2016 ; Mazza2017 such as an orbital imbalance of the Coulomb interaction. As for the effects of the dynamical part of the Hamiltonian, it has been argued that the parametric driving of the Raman modes through nonlinear couplings with the mid-infrared phonon modesKnap2016 ; Komnik2016 or the excited mid-infrared phonon through a nonlinear electron-phonon couplingKennes2017 ; Sentef2017 effectively pushes the system into the stronger-coupling regime, leading to an enhancement of the SC. Still, many important questions remain to be answered. In particular, when the enhancement of the SC is attributed to a change of the time-averaged Hamiltonian or the effective interaction, it is often assumed that the system remains in thermal equilibrium. However, in general, a periodically driven and interacting many-body system is expected to heat up (either by resonant excitations, or on a longer timescale by multi-photon absorption processes) and eventually reach an infinite temperature state.Rigol2014 ; Ponte2015 This would imply a melting of long-range order after a long time, even if the pairing interactions are enhanced by the driving. In a real system, this detrimental heating effect may potentially be circumvented in two ways: Either the heating processes are slow enough so that an enhancement of a long-range order can be observed at least transiently, or the heating is eventually balanced by energy dissipation to the environment, thus leading to the so-called nonequilibrium steady state (NESS). The transient dynamics has been studied by looking at the SC instabilities of the normal phase, in the form of negative eigenvalues of the transient pairing susceptibility.Knap2016 ; Babadi2017 While such negative susceptibilities can occur for short times, it is an open question on which time scale the long-range order could evolve in such a situation, and whether a potential increase can be fast enough to overcome the heating. Even less is known about the nature of the NESSs, where excitation and dissipation are balanced. For an understanding of the light-enhanced SC, it is therefore essential to go beyond the study of effective interactions and Hamiltonians, and to consider a formalism that can take account of the transient or steady-state modifications of the electron and phonon distribution functions, and of the SC order or fluctuations, in the nonequilibrium state.

In this paper, we focus on the dynamical aspect of a phonon-driven system and discuss the influence of the parametric phonon driving and the resulting heating and nonthermal distribution on the SC. Specifically, we study the time-periodic nonequilibrium steady state (NESS) and the nonequilibrium transient dynamics of the Holstein model, a prototypical model of electron-phonon systems, using the nonequilibrium dynamical mean-field theory (DMFT) combined with the Migdal approximation.Aoki2013 ; Murakami2015 ; Murakami2016 ; Murakami2016b To investigate the time-periodic NESS of the electron-phonon (el-ph) system with externally driven phonons, we extend the Floquet DMFT formalism, which has been previously only applied to purely electronic systems,Schmidt2002 ; Joura2008 ; Tsuji2008 ; Tsuji2009 ; Lee2014 ; Mikami2016 to the Holstein case.

The Floquet DMFT is based on the Floquet Green’s function formalism, which can naturally take into account the effect of many-body interactions, periodic driving by external fields and a dissipative coupling to the environment on equal footing. The Floquet DMFT for el-ph systems developed here is therefore not limited to the SC order, and should be useful to study a wide range of nonequilibrium situations that involve phonon excitations triggered by “nonlinear phononics”.Forst2011a

Our systematic analysis, which covers the weak and strong coupling regimes, reveals that the parametric phonon driving generally suppresses the SC (see Fig. 8 and Fig. 15), even though in some parameter regimes the retarded attractive interaction is increased and an enhancement of the SC is naively expected (see Fig. 4(a)(d)). In particular in the parametrically resonant regime, where the system absorbs energy most efficiently independent of the coupling strength, the SC order and fluctuations are strongly suppressed.

The paper is organized as follows. Section II describes the setup for the Holstein model under external phonon driving and attached to thermal baths. The framework used to study the problem is explained in Sec. III and the nonequilibrium Green’s functions are introduced in Sec. III.1. In Sec. III.2, we discuss the heat bath directly coupled to the phonon degrees of freedom, and in Sec. III.3 we introduce the Floquet Green’s function and derive its expression for the case of parametric phonon driving. In Sec. III.4, we explain the Floquet DMFT for the Holstein model and provide expressions for relevant physical quantities. Important remarks about the usage and limitations of the Floquet DMFT for the parametric phonon driving problem can be found in Sec. III.5. The results of our systematic analysis are shown in Sec. IV. In Sec. IV.1, we focus on the modification of the phonon propagator in the driven state in the weak coupling limit, and show that an enhancement of the attractive interaction is realized in some parameter regimes, which naively suggests an enhancement of SC. Sections IV.2 and IV.3 are dedicated to the strong coupling regime. In Sec. IV.2, we study the NESS under the parametric phonon driving using the Floquet DMFT and in Sec. IV.3 we discuss the transient dynamics toward the NESS. In Sec. IV.4, we focus on the weak coupling regime, and study how the SC fluctuations above $T_{c}$ react to the phonon modulation. A summary and conclusions are provided in Sec. V

II Model

In this paper, we consider the case where the system coupled to heat baths is periodically driven by an external field. Such a situation is expressed by the Hamiltonian

[TABLE]

see Fig. 1(a). Here $H_{\rm sys}(t)$ is the Hamiltonian for the system, which includes terms representing periodic excitations. $H_{\rm bath}$ is the Hamiltonian for heat baths, while $H_{\rm mix}$ represents the coupling between the system and the baths. When the system is continuously excited by periodic fields, it is expected that it finally reaches a time-periodic nonequilibrium steady state (NESS) in which the energy injection from the field and the energy flow into the baths are balanced. This situation is different from isolated systems under periodic driving.Rigol2014 ; Ponte2015 ; Kuwahara2016 ; Mori2016

The system Hamiltonian $H_{\rm sys}(t)$ consists of the unperturbed part $H_{\rm sys,0}$ and the external periodic field, $H_{\rm ext}(t)$. Here, we focus on a simple but relevant example of an electron-phonon coupled system, the Holstein model, with the Hamiltonian

[TABLE]

where $c_{i}^{\dagger}$ is the creation operator for an electron with spin $\sigma$ at site $i$, $v$ is the electron hopping, $\mu$ is the electron chemical potential, $n_{i}=c_{i,\uparrow}^{\dagger}c_{i,\uparrow}+c_{i,\downarrow}^{\dagger}c_{i,\downarrow}$, $\omega_{0}$ is the bare phonon frequency, $a_{i}^{\dagger}$ is the creation operator for the Einstein phonon, and $g$ is the el-ph coupling. The retarded attractive interaction between electrons is mediated by the phonons, which leads to an s-wave superconducting state at low enough temperatures. The SC order parameter is defined as $\phi=\frac{1}{N}\sum_{i}\langle c_{i\downarrow}c_{i\uparrow}\rangle$ which is assumed to be real without loss of generality. The strength of the phonon-mediated attractive interaction is characterized by the dimensionless el-ph coupling, $\lambda=-\rho_{0}(0)g^{2}D^{R}(\omega=0)$, where $\rho_{0}(\omega)$ is the density of states for the free electrons, and $D$ is the retarded phonon propagator.Bauer2011

Recent developments in THz or mid-infrared laser techniques make it possible to selectively and strongly excite a certain phonon mode. The infrared active phonons can couple to the electron sector through further excitation of another phonon mode (Raman mode) that is linearly coupled to the electron sectorSubedi2014 ; Knap2016 or through a nonlinear-coupling mechanism.Mitrano2016 ; Singla2015 ; Kennes2017 Our present interest is in the former case. Considering the mid-infrared phonon oscillations as an external driving force to the Raman mode through the nonlinear coupling, we study the following type of excitations for the phonon sector:

[TABLE]

Here $H_{\rm ph}(t)$ consists of the time-dependent part ($H_{\rm ext}(t)$) and the time-independent part ($\omega_{0}\sum_{i}a^{\dagger}_{i}a_{i}$ in $H_{\rm sys,0}$). $\omega_{0P}(t)$ and $\omega_{0X}(t)$ are periodic in time with a period $\mathscr{T}\equiv\frac{2\pi}{\Omega}$, $X_{i}=a_{i}+a_{i}^{\dagger}$ and $P_{i}=(a_{i}-a_{i}^{\dagger})/i$. This type of modulation of the phonon degrees of freedom has been discussed as a potential origin of the enhancement of superconductivity.Knap2016 ; Komnik2016 ; Babadi2017

To be specific, we focus on two types of periodic phonon excitations,

[TABLE]

and

[TABLE]

The first type, Eq. (4), is included in excitations through some of the nonlinear couplings between the mid-infrared mode and the Raman mode. Its potential effect on superconductivity has been discussed in Ref. Knap2016, . Potential effects of the second type of excitations, Eq. (5), have been discussed in Ref. Komnik2016, . As we will show below, it turns out that the type 1 and type 2 excitations give qualitatively the same results concerning the enhancement of the attractive interaction and the suppression of SC in the NESS and the transient dynamics. Besides these two types of excitations, one may also want to consider the case where $\frac{1}{\mathscr{T}}\int^{\mathscr{T}}_{0}\omega_{0X}(t)dt\neq\omega_{0}$ as in Ref. Knap2016, . If $\frac{1}{\mathscr{T}}\int^{\mathscr{T}}_{0}\omega_{0X}(t)dt<\omega_{0}$, the system goes to the stronger-coupling regime and the SC should be enhanced. However, since this effect is rather trivial we do not discuss it in this paper.

As for the bath, we use free electron and free phonon baths, which are introduced in detail in the following section.

III Formalism

III.1 Nonequilibrium Green’s functions

The transient time evolution toward the NESS can be described in the nonequilibrium Green’s function method formulated on the Kadanoff-Baym (KB) contour $\mathcal{C}_{\rm KB}$,Aoki2013 where we define the contour-ordered Green’s functions,

[TABLE]

for the electrons and phonons, respectively. Here $c^{\dagger}$ is the creation operator for an electron in the system, $X$ represents the phonon distortion and $\mathcal{T}_{\mathcal{C}_{\rm KB}}$ is the contour-ordering operator on the KB contour. All the operators are in the Heisenberg picture, and $\langle\cdots\rangle$ indicates the average with respect to the initial ensemble. In this formalism, the system is in equilibrium at the initial time $t_{0}$ and is then driven by the periodic external field. In Secs. IV.3,IV.4 we use this formalism to study the transient dynamics. The DMFT formalism for the real-time dynamics, both in electron and electron-phonon-coupled systems, has been discussed in previous worksAoki2013 ; Werner2013 ; Murakami2015 ; Murakami2016 and will not be presented in detail here.

The KB formalism can be inefficient for studying the NESS, which the system reaches in the long-time limit. Since the initial correlations are wiped out through the coupling to the thermal bath, they should be irrelevant to the NESS.Aoki2013 Hence, we can simplify the problem by assuming $t_{0}=-\infty$ and neglecting the left-mixing, right-mixing and Matsubara components in the KB formalism (which reduces to the Keldysh formalism), and focus directly on the NESS. 111We assume that at $t_{0}=-\infty$ the system is free and decoupled from the baths, and that the interactions and the coupling to the baths are adiabatically turned on. In the Keldysh formalism, we focus on the contour $\mathcal{C}_{K}=\mathcal{C}_{1}\cup\mathcal{C}_{2}$ [see Fig. 1(b)], and hence only need to consider

[TABLE]

Here $G^{ij}(t,t^{\prime})$ indicates $t\in\mathcal{C}_{i}$ and $t^{\prime}\in\mathcal{C}_{j}$.

The physical representation for the Green’s function is defined as

[TABLE]

where $\check{\sigma}_{3}$ is the third component of the Pauli matrices. This is the so-called Larkin-Ovchinnikov form. We can apply the same transformation to the phonon Green’s function $D$.

If there is no interaction or coupling to external baths, one can evaluate the bare Green’s function $G_{0}$ for electrons ($D_{0}$ for phonons), even in the presence of an external field. The correlations coming from the interaction and baths can be taken into account through the self-energy $\Sigma$ for electrons ($\Pi$ for phonons). Since the initial correlation can be neglected, the Dyson equation in the Keldysh formalism becomes simple:

[TABLE]

Here $\cdot$ indicates the matrix product, $\underline{I}$ is the $2\times 2$ identity matrix, $\int dt_{1}\underline{G}_{0}^{-1}(t,t_{1})\underline{G}_{0}(t_{1},t^{\prime})=\underline{I}\delta(t-t^{\prime})$ (the same holds for $\underline{D}_{0}^{-1}$). Precisely speaking, Eq. (9b) is true as long as the average phonon distortion $\langle X(t)\rangle$ remains zero, which is the case considered in this paper. When $\langle X(t)\rangle$ is finite, the phonon propagator has an additional contribution $-i\langle X(t)\rangle\langle X(t^{\prime})\rangle$.

III.2 Thermal baths

In this paper, we consider two types of thermal baths. One is attached to the electron degrees of freedom and the other is attached to the phonon degrees of freedom.

For the bath attached to the electrons we employ the so-called Büttiker model, which consists of multiple sets of free electrons in equilibrium connected to each site of the system. The model has been used in the previous Floquet DMFT studies for electron systems.Tsuji2009 ; Mikami2016 ; Aoki2013 Here we briefly review the model. Due to the noninteracting nature of the bath, it can be traced out exactly to yield the self-energy correction $\Sigma_{{\rm bath}}$, which is characterized by

[TABLE]

All of the physical components of $\Sigma_{{\rm bath}}$ can be obtained through the fluctuation-dissipation theorem for $\Sigma_{\rm bath}(\omega)$.Aoki2013 The total self-energy is the sum of the contributions from the bath and the interactions in $H_{\rm sys}$,

[TABLE]

In the following, we use a heat bath with a finite band width,

[TABLE]

to avoid the divergence of the bath self-energy in the time representation at $t=t^{\prime}$. We note that this bath tends to suppress the SC order, since it reduces the lifetime of quasiparticles.

As for the heat bath attached to the phonon degrees of freedom, we consider a free-boson bath as in the Caldeira-Leggett model,Caldeira1981 which is expressed as

[TABLE]

Here $b^{\dagger}_{i,p}$ is the creation operator for a bath phonon with an index $p$, which is attached to the system phonon at site $i$. We can trace out the bath phonons in the path-integral formalism, which yields an additional self-energy correction to the phonon Green’s function in the system,

[TABLE]

where $D_{0,p}$ is the bare bath phonon propagator in equilibrium, $f_{b}$ is the boson distribution function with the bath temperature $T$, and $\theta_{c}(t,t^{\prime})$ is the Heaviside step function along the Keldysh contour. We note that the total phonon self-energy includes the contribution from the bath and that from the interaction in $H_{\rm sys}$,

[TABLE]

Since $\Pi_{{\rm bath}}(t,t^{\prime})$ is time-translation invariant, we can consider its Fourier components. The retarded part is

[TABLE]

and the bath information is contained in its spectrum

[TABLE]

When $B_{\rm bath}(\omega)$ is given, the self-energy from the bath can be expressed as

[TABLE]

In the following, we consider a bath spectrum of the form

[TABLE]

Here we assume $\gamma=\omega_{D}$, which yields an almost linear structure in $B_{\rm bath}(\omega)$ around $\omega\in[-\omega_{D},\omega_{D}]$. This can be regarded as an Ohmic bath with a soft cut-off. Because of the non-vanishing real part of the bath self-energy, the bath coupling softens the phonon frequency, leading to a stronger attractive interaction. We put $\gamma_{\rm ph}=\mathcal{V}^{2}/\omega_{D}$ in the following.

III.3 Floquet Green’s functions

We apply the Keldysh formalism introduced in Sec. III.1 to a time-periodic NESS realized under an external periodic field, where we assume that relevant physical observables become periodic in time. We can take advantage of this to simplify the Keldysh formalism. For example, the Green’s function becomes periodic with respect to the average time $t_{\rm av}$,

[TABLE]

where $t_{\rm r}=t-t^{\prime}$ and $t_{\rm av}=\frac{t+t^{\prime}}{2}$. The self-energies have the same property in the NESS. Hence, the full $t_{\rm av}$-dependence of $G^{\rm R,A,K}(t_{\rm r};t_{\rm av})$ and $D^{\rm R,A,K}(t_{\rm r};t_{\rm av})$ contains redundant information. We can eliminate this redundancy by considering the Floquet representation of these functions in the following manner. Let us assume that a function $F(t_{\rm r};t_{\rm av})$ has the periodicity of Eq. (21). The Floquet representation (the Floquet matrix form) of this function is defined as Aoki2013

[TABLE]

This representation has some important properties. At each $\omega$, ${\bf F}_{mn}(\omega)$ can be regarded as a matrix whose indices are $m$ and $n$ (the Floquet matrix). In this representation, the convolution of two functions satisfying Eq. (21) in the two-time representation becomes the product of the Floquet matrices at each $\omega$. Aoki2013 This simplifies the Dyson equation, Eq. (9), to

[TABLE]

Here $\underline{\bf G}(\omega)$ indicates the Larkin-Ovchinnikov form with each physical component expressed in the Floquet representation, and $\cdot$ is the matrix product with respect to the Floquet representation and the Larkin-Ovchinnikov form. Hence the Dyson equation, which has been originally a Volterra equation in the two-time representation, becomes a matrix inversion problem in the Floquet representation. In addition, from the definition it follows that $F_{m,n}(\omega)=F_{m+l,n+l}(\omega-l\Omega)$. We deal with this redundancy by using the reduced Brillouin-zone scheme restricting the frequency to $\omega\in[-\frac{\Omega}{2},\frac{\Omega}{2})$.

When an external periodic field is applied to the electrons, we need to incorporate this effect into the bare electron propagator, which acquires a time-periodic structure as in Eq. (21). The corresponding expressions have been discussed in the previous works for the single-bandTsuji2008 and multi-bandMikami2016 systems, so we do not repeat this here. Instead, we focus on the direct phonon modulation originating from a THz or mid-infrared laser pump, Eq. (3). In this case, we need to first derive the expression for the bare phonon propagator in the presence of such a modulation, in order to understand the effect of direct phonon excitations and to construct the perturbation theory on top of it.

The free phonon system under the periodic driving is described by the Hamiltonian

[TABLE]

The operator $X$ in the Heisenberg representation satisfies the equation

[TABLE]

From this, we find that

[TABLE]

satisfies

[TABLE]

Hence, in the Dyson equation for the full phonon Green’s function, Eq. (9b), we can use $\underline{D}_{0}^{-1}(t,t^{\prime})=\underline{I}D^{-1}_{0}(t,t^{\prime})$. In the Floquet representation, $D^{-1}_{0}(t,t^{\prime})$ becomes

[TABLE]

where

[TABLE]

We note that the inverse of the Floquet matrix of the retarded part, $[{\bf D}^{R}_{0}(\omega)]^{-1}$, is equal to $[{\bf D}_{0}(\omega+i0^{+})]^{-1}$. In the numerical implementation, we solve the Dyson equation for the phonon Green’s function using Eq. (28) and Eq. (23b).

Now we focus on more specific situations described by the driving protocols of type 1 [Eq. (4)] and type 2 [Eq. (5)]. For the excitation of type 1 [Eq. (4)], the inverse of the bare phonon Green’s function is obtained by substituting

[TABLE]

into Eq. (28). Since the equation of $D^{R}_{0}(t,t^{\prime})$ turns out to be the Mathieu’s differential equation (see also Sec. III.5), we can express it as a linear combination of its two Floquet solutions.

For the excitation of type 2 [Eq. (5)], the inverse of the bare phonon Green’s function is obtained from

[TABLE]

where we put $r=m-n$ and ${}_{2}F_{1}(\alpha,\beta;\gamma;z)$ is the hypergeometric function. One can also easily express $D^{R}_{0}$ in the time domain as

[TABLE]

where $S(t)\equiv t\omega_{0}+\frac{\Delta\omega_{0}}{\Omega}\sin(\Omega t)$.

III.4 Dynamical mean-field theory

In order to evaluate the full Green’s functions, we employ the dynamical mean-field theory (DMFT), Metzner1989 ; Georges1992 ; Georges1996 ; Aoki2013 which becomes exact in the limit of infinite spatial dimensions. Georges1996 The idea of DMFT is to map the lattice problem to an effective impurity problem. In the present case of the Holstein model with thermal baths, the form of the effective impurity model in the path-integral formalism is

[TABLE]

where

[TABLE]

$\int_{\mathcal{C}_{K}}$ denotes an integral on the Keldysh contour, $\Psi^{\dagger}(t)\equiv[c_{\uparrow}^{\dagger}(t),c_{\downarrow}(t)]$ a Nambu spinor, and a hat symbol a $2\times 2$ matrix in the Nambu form. In this expression, the thermal baths and $P$ (phonon momentum) in the impurity problem have been integrated out.Murakami2016 The hybridization function $\hat{\Delta}$ is determined such that the impurity Green’s function, $\hat{G}_{\text{imp}}(t,t^{\prime})=-i\langle T_{\mathcal{C}_{K}}\Psi(t)\Psi^{\dagger}(t^{\prime})\rangle$, and the impurity self-energy $\hat{\Sigma}_{\rm imp,int}$ are identical to the local lattice Green’s function for the electrons, $\hat{G}_{\rm loc}(t,t^{\prime})=-i\langle T_{\mathcal{C}_{K}}\Psi_{i}(t)\Psi_{i}^{\dagger}(t^{\prime})\rangle$, and the local lattice self-energy, respectively.

Since we are considering a time-periodic NESS, the effective impurity model also has the same time periodicity, where the hybridization function, $\hat{\Delta}$, is assumed to show the periodic behavior of Eq. (21). Hence the Dyson equation involved in the solution of the effective impurity model can be dealt with as an inverse problem of the Floquet matrix. Precisely speaking, the matrix to invert involves three indices, i.e., those of the Floquet, Larkin-Ovchinnikov, and Nambu forms.

III.4.1 Migdal approximation

In order to solve the effective impurity model, we employ the self-consistent Migdal approximation, Murakami2015 ; Murakami2016 ; Murakami2016b ; Bauer2011 ; Freericks1994 ; Capone2003 ; Hague2008 ; Leeuwen2015 ; Pavlyukh2016 which is justified when the phonon frequency $\omega_{0}$ is small enough compared to the electron bandwidth.Bauer2011 ; Freericks1994 ; Capone2003 ; Hague2008 Here, the electron self-energy ($\hat{\Sigma}$) and phonon self-energy ($\Pi$) are given by

[TABLE]

where we do not explicitly write the Hartree term. Remember that in DMFT $\hat{\Sigma}_{\rm int}(t,t^{\prime})$ is identified with the lattice self-energy of the electrons at each momentum. The Hartree term is proportional to $\langle X(t)\rangle$, and, in the half-filled case, which we focus on in this paper, it is exactly zero. Away from half-filling, $\langle X(t)\rangle$ can be finite and time dependent. Still, since we are considering homogeneous excitations, this term can be regarded as a time-dependent chemical potential and does not affect the dynamics of relevant physical quantities without couplings to the Büttiker-type bath. The self-consistent Migdal approximation is a conserving approximation, which allows one to trace the energy flow from the system to the thermal baths.

III.4.2 Observables

We summarize the physical observables that are used to discuss the results in the paper. Let us denote the Heisenberg representation as $\mathcal{O}(t)=\mathcal{U}_{\text{tot}}(t_{0},t)\mathcal{O}\mathcal{U}_{\text{tot}}(t,t_{0})$, where $\mathcal{U}_{\rm tot}(t,t_{0})$ is the unitary evolution operator from time $t_{0}$ to $t$.

The change of the system energy can be expressed as

[TABLE]

The first term represents the work done by the external field, which we express as $W(t)$. The second term represents the energy dissipation from the system to the bath, which we express as $-I(t)$. If the system is in a time-periodic NESS, the average energy injected from the field and dissipated to the baths should be the same,

[TABLE]

where the overline indicates the time average over one period, $\frac{1}{\mathscr{T}}\int_{0}^{\mathscr{T}}dt$. We also note that the energy dissipated from the system should be equal to the energy injected to the bath. This leads to

[TABLE]

see also Appendix A for more details.

In the following, we give the explicit expression for each term in Eq. (37) using the self-energies and Green’s functions.

(1) The work done by the field is

[TABLE]

where $\langle X_{i}^{2}(t)\rangle$ can be obtained from the lesser part of the phonon Green’s function.

(2) The dissipated energy is

[TABLE]

where

[TABLE]

and

[TABLE]

Here the symbol $*$ indicates the convolution along the Keldysh contour $\mathcal{C}_{K}$ and the symbol $\circledast$ represents a multiplication of the 2$\times$2 Nambu matrices and the convolution i.e., $[\hat{A}\circledast\hat{B}]_{\alpha,\alpha^{\prime}}(t,t^{\prime})=\sum_{\gamma}\int_{\mathcal{C}}d\bar{t}A_{\alpha,\gamma}(t,\bar{t})B_{\gamma,\alpha^{\prime}}(\bar{t},t^{\prime})$. ${\rm Tr}$ refers to the trace of Nambu matrices.

The system observables (that $H_{\rm sys}$ involves) can be expressed in the following manner, in analogy to the Kadanoff-Baym formalism.Murakami2015

1. Kinetic energy

[TABLE]

2. Interaction energy

[TABLE]

3. Phonon energy

[TABLE]

Here $D_{\rm PP}(t,t^{\prime})=-i\langle T_{c}P(t)P(t^{\prime})\rangle$. This can be obtained by taking derivatives of $D(t,t^{\prime})$.Murakami2015 The Larkin-Ovchinnikov form of this function becomes

[TABLE]

We note that the derivative in the two-time representation becomes a multiplication of factors depending on $\omega,m$ and $n$ in the Floquet representation , which can be used for the numerical implementation.

Another important quantity is the optical conductivity, i.e. $\sigma_{\alpha\alpha}(t,t^{\prime})\equiv\delta\langle j_{\alpha}(t)\rangle/\delta E_{\alpha}(t^{\prime})$, where $\langle j_{\alpha}(t)\rangle$ is the current for the direction $\alpha$, $E_{\alpha}(t^{\prime})$ is the probe field, and $\delta$ indicates the functional derivative.Aoki2013 This can be expressed in terms of the susceptibility $\chi_{\alpha\alpha}(t,t^{\prime})\equiv\delta\langle j_{\alpha}(t)\rangle/\delta A_{\alpha}(t^{\prime})$ as

[TABLE]

$\chi_{\alpha\alpha}(t,t^{\prime})$ consists of the diamagnetic term $\chi^{\rm dia}_{\alpha\alpha}(t,t^{\prime})$ and the paramagnetic term $\chi^{\rm para}_{\alpha\alpha}(t,t^{\prime})$.Aoki2013 In principle, the paramagnetic term includes the vertex correction, but in the present case, since the phonon driving does not break the parity symmetry (${\bf k}\leftrightarrow{\bf-k}$), its contribution vanishes. In the end, we obtain the following expression in the Nambu form,

[TABLE]

where $\epsilon$ denotes the energy of noninteracting electrons, $\Phi(\epsilon)=\frac{1}{\Omega_{\rm vol}}\sum_{{\bf k},\sigma}\left(\frac{\partial\epsilon_{\bf k}}{\partial k_{\alpha}}\right)^{2}\delta(\epsilon-\epsilon_{\bf k})$, $\frac{d\Phi(\epsilon)}{d\epsilon}=\frac{1}{\Omega_{\rm vol}}\sum_{\bf k}\frac{\partial^{2}\epsilon_{\bf k}}{\partial k_{\alpha}^{2}}\delta(\epsilon-\epsilon_{\bf k})$, and $\Omega_{\rm vol}$ is the volume of the system. In particular, for the Bethe lattice, $\Phi(\epsilon)=(N/3d)[(W/2)^{2}-\epsilon^{2}]\rho_{0}(\epsilon)$ is used, where $N$ is the system size and 2$d$ the coordination number.Bethe_trans

In the NESS, $\sigma(t,t^{\prime})$ and $\chi(t,t^{\prime})$ become time periodic, and their Floquet expressions are connected through

[TABLE]

To obtain an intuitive picture, let us imagine that $E_{\alpha}(t)=\delta E_{0}e^{-i\nu t}$ is applied as a probe field. The induced current is

[TABLE]

We note that in this expression we unfold the restricted zone scheme for $\boldsymbol{\sigma}^{ij}_{m0}(\nu)$ by using $\boldsymbol{\sigma}_{m0}(\nu)=\boldsymbol{\sigma}_{m+l,l}(\nu-l\Omega)$. Equation (51) means that the response includes $e^{-i(\nu+m\Omega)t}$ components in addition to $e^{-i\nu t}$. Previous studies Tsuji2009 ; Mikami2016 have mainly focused on the diagonal component ($m=0$), while in the present work, we will also discuss the off-diagonal components.

III.5 Breakdown of the Floquet Green’s function method

It is important to remark that there are instabilities in the system under the excitation Eq. (4). Without the heat bath and the el-ph coupling, the equation for the phonon Green’s function [Eq. (27) combined with Eq. (4)] becomes the so-called Mathieu’s differential equation. It is known that the solution of this equation is unstable around $\frac{n}{2}\Omega=\omega_{0}$, with $n$ a natural number . In this unstable regime, the amplitude of oscillations in the solution of the equation monotonically increases, and there is no well-defined Floquet form of the Green’s function since we cannot define the Fourier components.

The heat bath and the el-ph coupling introduce a damping effect through the self-energy corrections. Even though the damping effect can reduce the unstable regime,Zerbe1994 ; Zerbe1995 the Green’s function still exhibits an instability when the strength of the external field is large enough. In such cases the Floquet Green’s function cannot be defined and the Floquet formalism becomes inapplicable. Physically, this means that the NESS, if any, is unstable against infinitesimally small perturbations to the phonon part.

In Fig. 2, we demonstrate that even with the phonon bath the behavior of the retarded part of the phonon Green’s function changes qualitatively as a function of driving amplitude. The correct results are obtained by solving the Dyson equation in the time domain (differential-integral Volterra equation), Eq. (9b), which produces either a stable or unstable solution depending on the excitation parameter. For comparison, we solve the Dyson equation in the Floquet form, Eq. (23b), and carry out an inverse transformation to the two-time representation. These two ways of solving the Dyson equation give the same solution, when the solution is stable. However, we have found an unphysical solution for the unstable regime in the Floquet form that is different from the transient solution of the Volterra equation in the time domain. In particular, the unphysical solution violates the initial condition (the sum rule in the $\omega$ space) $\partial_{t}D_{0}^{R}(t,t^{\prime})|_{t^{\prime}=t}=-2\omega_{0}$, see also Sec. IV.1.

Therefore, when using the Floquet form of the Green’s function, one has to be careful about the unstable regime, where the system cannot reach a physically stable NESS. The excitation Eq. (4) at the parametric resonance $\frac{\Omega}{2}=\omega_{0}$ is potentially in such a regime. However, there is no such problem in the calculation of the transient time evolution on the Kadanoff-Baym contour.

IV Results

In this paper, we investigate the effect of parametric phonon driving on superconductivity (SC) both in the weak coupling ($\lambda\ll 1$, the red region in Fig. 3) and strong coupling ($\lambda\sim 1$, the blue region in Fig. 3) regimes. In this paper, we use the semi-elliptic density of states, $\rho_{0}(\omega)=\frac{1}{2\pi v^{2}_{\ast}}\sqrt{4v^{2}_{\ast}-\epsilon^{2}}$, i.e. the Bethe lattice. We set $v_{\ast}=1$, i.e. the electron bandwidth $W$ is 4, and we focus on the half-filled case. We choose the bare phonon frequency $\omega_{0}=0.4$, which is much smaller than the electron bandwidth. If we take the bandwidth as $0.4$ eV, the bare phonon frequency is $40$ meV, corresponding to a period of about $100$ fs. The gap size of the SC spectrum around $\lambda\sim 1$ is approximately $10$ meV, i.e. this period is about $400$ fs. Our unit temperature corresponds to $0.1$ eV or 1160 K.

IV.1 Effective attractive interaction in the weak coupling limit : Bare phonon propagator

In this section, we discuss how the effective attractive interaction mediated by phonons, which is responsible for the SC, is modulated under the phonon driving in the weak coupling limit. Our discussion is based on the Migdal-Eliashberg theory, where the retarded attractive interaction mediated by the phonons is represented by the phonon propagator, $g^{2}D^{R}(\omega)$, and $g^{2}D^{R}(0)$ is conventionally taken as a measure for the strength of the attractive interaction. In the weak coupling limit, the renormalization of the phonon (the self-energy correction) is small and $g^{2}D^{R}_{0}(\omega)$ serves as a measure of the effective attractive interaction. With this idea in mind, the time-averaged $\bar{D}^{R}_{0}(0)=\frac{1}{\mathscr{T}}\int_{0}^{\mathscr{T}}dt_{\rm av}\int^{\infty}_{-\infty}dt_{\rm r}D^{R}_{0}(t_{\rm r};t_{\rm av})$ has been analyzed for the type 2 excitation in Ref. Komnik2016,

The effect of the type 1 excitation, on the other hand, has been discussed in Ref. Knap2016, using the polaron picture and the Lang-Firsov transformation. The polaron approximation works well when the phonon timescale is comparable or faster than the electronic timescale. Murakami2013a Within this framework, the hopping of polarons is renormalized from the bare electron hopping parameter by the Frank-Condon factor, $\exp(-g^{2}/\omega_{0}^{2})$, and the effective attractive interaction is $-\lambda_{0}=\frac{2g^{2}}{\omega_{0}}$. However, in many materials the phonon frequencies are much smaller than the bandwidth, in which case the Migdal-Eliashberg theory is applicable. We note that these two different theories can lead to different predictions concerning the behavior of the SC order. The isotope effect serves as an example. By the substitution of isotopes, usually $\omega_{0}$ decreases as $\propto\frac{1}{\sqrt{M}}$ ($M$ is the mass of the atom), while $\lambda_{0}$ does not change. In the polaron picture, the system moves to the stronger coupling regime, since the Frank-Condon factor decreases. The convex shape of the superconducting phase boundary in the weak-coupling regime (see Fig. 3) thus implies an increase of the transition temperature. On the other hand, the BCS theory, to which the Migdal-Eliashberg theory reduces in the weak coupling regime, predicts the opposite behavior. Here, the transition temperature is proportional to $\omega_{0}\exp(-\frac{1}{\lambda_{0}\rho_{0}(0)})$, which decreases if $\omega_{0}$ is reduced by the substitution of isotopes. Hence for realistic phonon frequencies we should discuss the effects of parametric phonon excitations within the Migdal-Eliashberg theory.

In Fig. 4, we show the numerical results for the time-averaged phonon propagator, $\bar{D}^{R}_{0}(\omega)=\frac{1}{\mathscr{T}}\int_{0}^{\mathscr{T}}dt_{\rm av}\int^{\infty}_{-\infty}dt_{\rm r}D^{R}_{0}(t_{\rm r};t_{\rm av})e^{i\omega t_{\rm r}}$, which corresponds to the diagonal components in the Floquet representation, for both types of excitations. In the evaluation, we formally invert the Floquet matrix of $D^{\rm R-1}_{0}$. Therefore, the results also include unphysical ones as we explained in Sec. III.5, but we can detect the unstable regime from the Mathieu’s equation. In Fig. 4(a)(d), we plot the change in the effective attractive interaction expressed as $({\rm Re}\bar{D}^{R}_{0}(0)-{\rm Re}D_{\rm 0,eq}(0))/{\rm Re}D_{\rm 0,eq}(0)$ in the plane of $\Delta\omega_{0}$ and $\Omega$. It turns out that both excitations result in a similar behavior. There is a strong enhancement of the attractive interaction $\bar{D}^{R}_{0}(0)$ at $\Omega\lesssim\omega_{0}$. On the other hand, for $\Omega\gtrsim\omega_{0}$, the attractive interaction is strongly suppressed and it can even become repulsive as is analytically shown in Ref. Komnik2016, .

This behavior is related to the position of the Floquet sidebands of the phonons in the phonon spectrum, $\bar{B}_{0}(\omega)\equiv-\frac{1}{\pi}{\rm Im}\bar{D}^{R}_{0}(\omega)$. At $\Omega\lesssim\omega_{0}$, the first replica peak ($\sim\omega_{0}-\Omega$) is just above $\omega=0$. At $\Omega\gtrsim\omega_{0}$, the first side peak moves to $\omega<0$ and there emerges a negative peak in the spectrum at $\omega>0$, which is the first sideband of the phonon peak at $-\omega_{0}$, in the time-averaged spectrum $\bar{B}_{0}(\omega)$, see Fig. 4 (c)(f). Since $\bar{D}^{R}_{0}(\omega)=\int d\omega^{\prime}\frac{\bar{B}_{0}(\omega^{\prime})}{\omega-\omega^{\prime}+i0^{+}}$, in the former case ($\Omega\lesssim\omega_{0}$) the sideband acts as an additional phonon coupled to electrons with a positive frequency, which leads to an additional attractive interaction. In the latter case, the sideband resembles a negative frequency phonon attached to electrons, which yields an effective repulsive interaction. There are analogous singularities around $\Omega=\omega/n$ related to the zero crossing of the $n$-th Floquet sideband of the phonons. Since the spectral weight of such high-order sidebands is small, a substantial effect only appears when the excitation strength $\Delta\omega_{0}$ is large enough. Even though the strong enhancement of the attractive interaction is restricted to some energy window smaller than $\omega_{0}$, see Fig. 4(b)(e), in the weak coupling limit, the effect of the energy window is smaller than that of the strength of the attractive interaction as can be seen from the BCS expression for the transition temperature. The former is a prefactor, while the latter is in the exponential. Hence it is natural to expect that this may lead to an enhancement of SC assuming that the electron distribution is the same as in equilibrium.

Here we note that the unstable regime in the case of the type 1 excitation is located around $\Omega=\frac{2\omega_{0}}{m}$. The boundary of the unstable region is depicted by colored dotted lines in Fig. 4(a). Therefore, when $m=2n$, the unstable regime overlaps with the $\Omega$-regime in which $\bar{D}^{R}_{0}(0)$ exhibits drastic changes. On the other hand, when $m=2n-1$ the existence of the unstable regime does not leave any trace in $\bar{D}^{R}_{0}(0)$. It turns out, instead, that the effect of the instability emerges in the finite frequency part of $\bar{D}^{R}_{0}(\omega)$. As can be seen in the result for $\Omega=0.7$ and $\Omega=0.75$ in Fig. 4 (b)(c), there is a large change in $\bar{D}^{R}_{0}(\omega)$ around $\omega\sim\omega_{0}$. The results for $\Omega=0.75$ are unphysical as discussed in Sec. III.5, and the spectrum shows dominant negative weight at $\omega>0$, which leads to the violation of the sum rule, $2\omega_{0}=\int^{\infty}_{-\infty}\omega\bar{B}(\omega)d\omega$. Based on the behavior of $\bar{D}^{R}_{0}(0)$, we thus do not expect any enhancement of SC around the parametric resonant regime at $m=2n-1$ for the type 1 excitation, and the analysis of $\bar{D}^{R}_{0}(0)$ results in a picture which differs from the enhancement of SC at the parametric resonance discussed in Ref. Knap2016, . For the type 2 excitation, there is no unstable solution, and we do not observe any large change in $\bar{D}^{R}_{0}(\omega)$ near $\omega=\omega_{0}$, see Fig. 4 (e)(f).

Even though the increased attractive interaction $g^{2}\bar{D}^{R}_{0}(0)$ seems a plausible scenario for the enhancement of SC, this naive expectation ignores the following issues: (i) The system gains energy from the external field, which results in heating and nonthermal distribution functions. (ii) When the electron-phonon coupling is not small, the phonon spectrum is renormalized and broadened, and the behavior of $D$ might be different from $D_{0}$. (iii) $\bar{D}^{R}_{0}(0)$ is an average of $D^{R}_{0}(t_{\rm r};t_{\rm av})$, and hence includes information of infinite past and infinite future times. Therefore, it is not clear to what extent $\bar{D}^{R}_{0}(0)$ is relevant to the transient dynamics of the superconducting order.

In the following sections, using the nonequilibrium DMFT and the Holstein model, we provide a comprehensive study of the effects of the parametric phonon driving on SC which takes the above issues into account. We also note that it turns out that the type 1 and type 2 excitations give the same results on a qualitative level. Therefore, in the main text, we only show the results for the type 1 excitation, and present the results for the type 2 excitation in Appendix B.

IV.2 Nonequilibrium steady states under phonon driving

In this section, we study the effects of the phonon driving on SC in the strong electron-phonon coupling regime by studying the nonequilibrium steady states. In the following we use, as a representative set of parameters for the strong coupling regime, $g=0.41,\omega_{0}=0.4,\beta=120$ and the bath parameters $\gamma_{\rm ph}=0.02,\omega_{D}=0.6,\gamma_{\text{el}}=0.005,W_{\rm ebath}=2.0$. This gives $\lambda=1.08$, the renormalized phonon frequency $\omega_{\rm r}\sim 0.2$, and the inverse transition temperature $\beta_{c}=55$. $\gamma_{\text{el}}$ is chosen to be much smaller than the SC gap in order not to break the SC order, while $\gamma_{\rm ph}$ is chosen so that the damping of phonons (the phonon peak width in the spectrum) mainly originates from the electron-phonon interaction. We have also considered stronger $g$, different temperatures and both stronger and weaker $\gamma_{\rm ph}$ and $\gamma_{\rm el}$, and confirmed that the behavior discussed below does not change qualitatively.

First, in Fig. 5, we show $\text{Re}\bar{D}^{R}(\omega)$, which includes the effects of $\Omega$ and $\Delta\omega_{0}$ as well as the renormalization through the electron-phonon coupling. In equilibrium, because of the strong el-ph coupling, the sign change in $\text{Re}\bar{D}^{R}(\omega)$ around the renormalized phonon frequency, $\omega_{r}$, is smeared out, see $\beta=50$ (above $T_{c}$) in Fig. 5(a). Below $T_{c}$ [$\beta=120$ in Fig. 5(a)], even though $|\bar{D}^{R}(0)|$ remains almost unchanged, a strong renormalization occurs for $\omega\neq 0$. This is caused by the opening of the SC gap, within which the scattering between electrons and phonons is suppressed. The peak in $|\text{Re}\bar{D}^{R}(\omega)|$ around $\omega=0.08$ corresponds to this energy scale. This effect is known as the phonon anomaly in strongly coupled el-ph superconductors,Axe1973 ; Allen1997 ; Weber2008 ; Murakami2016 see also Fig. 6(c). In the NESS, as in the case of the bare phonon propagator, the sign of the change in the effective attractive interaction $g^{2}\bar{D}^{R}(0)$ depends on the excitation frequency $\Omega$, despite the strong renormalization and the absence of sharp structures. Namely, for $\Omega\lesssim\omega_{\rm r}$, the strength $|g^{2}\bar{D}^{R}(0)|$ increases [see Fig. 5(a)], while for $\Omega\gtrsim\omega_{\rm r}$, it decreases [see Fig. 5(b)]. However, this enhancement of $D$ turns out not to lead to an enhancement of the superconductivity.

Figure 6 illustrates how the (time-averaged) nonequilibrium spectrum and the distribution functions for the electrons and phonons change as we increase the amplitude of the excitation. Here we choose $\Omega=0.125$ for which an enhancement of the effective attractive interaction has been observed [Fig. 5(a)]. The nonequilibrium spectrum is defined as $\bar{A}(\omega)\equiv-\frac{1}{\pi}{\rm Im}\bar{G}^{R}(\omega)$ and the distribution function is $\bar{f}_{\text{el}}(\omega)\equiv\bar{N}(\omega)/\bar{A}(\omega)$ with $\bar{N}(\omega)=\frac{1}{2\pi}{\rm Im}\bar{G}^{<}(\omega)$. Here $\bar{G}^{R(<)}(\omega)=\frac{1}{\mathscr{T}}\int_{0}^{\mathscr{T}}dt_{\rm av}\int^{\infty}_{-\infty}dt_{\rm r}G^{R(<)}(t_{\rm r};t_{\rm av})e^{i\omega t_{\rm r}}$. As we increase the excitation amplitude $\Delta\omega_{0}$, the size of the SC gap decreases and eventually closes completely, see Fig. 6(a). At the same time, the slope of the nonthermal distribution function for the electrons decreases around $\omega=0$ and high energy states become more occupied, which roughly resembles the results for equilibrium states with higher temperatures, see Fig. 6(b).

One interesting feature in the nonequilibrium SC before melting is that there emerges a dip in $\bar{A}(\omega)$ and a hump in $\bar{f}_{\text{el}}(\omega)$ at $\omega=\Omega/2$. These features are only pronounced in the NESS with SC and disappear in the normal state. They can be explained as follows as a consequence of the two-band structure in the SC phase. The hump in $\bar{f}_{\text{el}}(\omega)$ results from excitations from the lower band of the Bogoliubov quasiparticles ($-E_{\rm k}$) to the upper band ($E_{\rm k}$), which become efficient at $E_{\rm k}=\frac{\Omega}{2}$, where $E_{\rm k}$ is the energy of the quasiparticles. Although Bogoliubov excitations can recombine or relax within the bands due to the scattering with phonons and the coupling to the environment, the continuous creation of particles at the particular energy $E_{\rm k}=\frac{\Omega}{2}$ can lead to an enhanced distribution in the steady state at the corresponding energy. The dip in the spectral function can be explained as a hybridization effect between the first Floquet sideband of the upper (lower) Bogoliubov band with the lower (upper) Bogoliubov band, which occurs at $E_{\rm k}=\frac{\Omega}{2}$. The hybridization opens a gap at the crossing point, which is not fully developed in our case because of the finite correlations, resulting in a dip in the spectrum.

The phonon spectra, which are defined in the same manner as the electron spectra using $D$, show a hardening of the renormalized mode as we increase the field strength [Fig. 6(c)], again consistent with the behavior seen in equilibrium with increased temperature.Murakami2015 At the same time, the phonon distribution function defined as $\bar{f}_{\rm ph}(\omega)\equiv-\frac{1}{2\pi}{\rm Im}\bar{D}^{<}(\omega)/\bar{B}(\omega)$ shows an enhancement of the phonon occupation [see the inset in Fig. 6(c)]. When we further increase the excitation strength, the Floquet sidebands of the phonons become more prominent [see $\Delta\omega=0.08$ data in Fig. 6(c)], in which case the spectrum does not resemble that at elevated temperatures.

Next, let us discuss how the nonthermal distribution is different from a thermal one. In Fig. 7, we show the comparison between the nonthermal distribution and the thermal distribution with an effective temperature derived from $\beta_{\text{eff}}=-4\partial_{\omega}\bar{f}_{\text{el}}(\omega)|_{\omega=0}$. We can see that this thermal distribution always underestimates the nonequilibrium distribution for $\omega>0$. Therefore, it is expected that superconductivity will be weaker than what is expected from this effective temperature, which is indeed the case, as shown below.

In Fig. 8, we summarize our representative Floquet DMFT results in the strong el-ph coupling regime for quantities time-averaged over one period. The SC order parameter is generally reduced from the equilibrium value (in fact, despite an extensive search, we could not find any parameter set where the SC order is enhanced in the strong-coupling regime, $\lambda\sim 1$). In particular, around $\Omega=0.4$, which is almost twice the renormalized phonon frequency, the superconducting order vanishes even for phonon modulations with small amplitude $\Delta\omega_{0}$. This is consistent with the amount of energy injected by the excitation, and the effective temperature estimated from the nonthermal distribution function $\partial_{\omega}\bar{f}_{\text{el}}(\omega)|_{\omega=0}$, see Fig. 8(b). We note that the peak position in Fig. 8(b) slightly shifts to higher energy as we increase the excitation strength. This is because the injected energy changes the renormalization of the phonons, which affects the value of the parametric resonance ($\sim 2\omega_{\rm r}$). We also note that the SC order parameter is always lower than the equilibrium value at $T_{\text{eff}}$. This is consistent with the fact that the equilibrium distribution function at $T_{\text{eff}}$ underestimates the occupancy of the states at $\omega>0$.

Let us briefly comment on the effect of stronger couplings to heat baths, since one expects in this case more efficient energy absorption. We have checked that stronger couplings to the phonon bath lead to a broader phonon peak in the spectrum, which results in a stronger suppression of SC and higher effective temperatures at $\Omega\lesssim\omega_{r}$. Stronger couplings to the electron bath generally lead to a lower effective temperature and a less dramatic suppression of SC. However, the equilibrium value of the SC order becomes smaller and in any case we only find a suppression of SC as far as we have investigated.

So far we have focused on the time-averaged values, whereas in the NESS all the quantities are oscillating around the averaged value during one period. In Fig. 9, we plot the trajectories of $E_{\rm kin}(t)$ and $E_{\rm int}(t)=E_{\rm nX}(t)+E_{\rm ph}(t)$ for the NESS over one period, and compare them to the averaged $\bar{E}_{\rm kin}$ and $\bar{E}_{\rm int}$ of the NESS (violet lines) and the equilibrium $E_{\rm kin}$ and $E_{\rm int}$ for different temperatures (black lines). First, we can see that the results for the time averaged values of $E_{\rm kin}$ and $E_{\rm int}$ are in general not on top of the equilibrium curve, and that the deviation strongly depends on the excitation frequency. This indicates that it is difficult to fully describe the properties of NESSs through an analogy to equilibrium states with elevated temperatures. We also notice that the oscillations around the averaged values are large, and that they become largest at the parametric resonance, $\Omega=0.4\simeq 2\omega_{r}$. Even at the times where the Hamiltonian becomes temporarily identical to the unperturbed one ($t=\frac{\pi}{2\Omega},\frac{3\pi}{2\Omega}$) the energies can strongly deviate from the thermal line. This indicates that the dynamics which occurs during one period cannot be captured using an adiabatic picture either in these cases.

In Fig. 10, we show the optical conductivity in the NESS. In equilibrium, for $T<T_{c}$, the real part shows a reduction of the spectral weight below the SC gap and the imaginary part exhibits a $1/\omega$ dependence, whose coefficient corresponds to the superfluid density. In the NESS, the optical conductivity $\sigma(t,t^{\prime})$, like any correlation function, can be written in the Floquet matrix form Eq.(22), see the discussion around Eq. (51). In the diagonal component in the NESSs, the weight of the real part inside the SC gap increases with increasing field strength, and becomes close to the result at elevated temperatures [see Fig. 10(a)]. Meanwhile, the $1/\omega$ component in the imaginary part is suppressed, which indicates the suppression of the superfluid density and is consistent with the closure of the gap and the decrease of the SC order parameter [see Fig. 10(b)]. In the off-diagonal component, one can also observe a $1/\omega$ around $\omega=0$ behavior both in the real and imaginary parts in the SC phase [see Fig. 10(c)(d)]. (We note that $\sigma_{m,n}(\omega)^{*}=\sigma_{-m,-n}(-\omega)$.) This is caused from the imperfect cancellation between $\chi^{\rm para}$ and $\chi^{\rm dia}$, and corresponds to the following asymptotic form for $|t-t^{\prime}|\rightarrow\infty$ of the conductivity,

[TABLE]

We note that $n_{S,m}$ can be a complex number and $n_{S,m}^{*}=n_{S,-m}$. This indicates that if a delta-function electric-field pulse is applied to the system, the induced current continues to flow forever (supercurrent) oscillating with a period of $2\pi/\Omega$, and the phase of $n_{S,m}$ determines the phase of the oscillations in the suppercurrent. The $1/\omega$ component in $\sigma_{m0}(\omega)$ corresponds to $n_{S,m}$. For example, without the phonon driving $n_{S,0}\approx 1.18$, while for $\Delta\omega_{0}=0.03$ $n_{S,0}\approx 0.62$ and $n_{S,1}\approx 0.05-i0.01$.

IV.3 Transient dynamics towards the NESS

In order to study the transient dynamics towards the NESS, we consider a situation where the system at time $t=0$ is in equilibrium (attached to the bath) and at $t>0$ is exposed to a phonon driving of type 1 with

[TABLE]

In order to follow the transient dynamics, we use the nonequilibrium DMFT formulated on the Kadanoff-Baym contour.Aoki2013 ; Murakami2015 ; Murakami2016 For this, we first need to evaluate the bare phonon propagator subject to the periodic driving by solving Eq. (27). Apart from that, the calculation is identical to the scheme explained in Refs. Murakami2015, ; Murakami2016, .

In Fig. 11, we show the transient dynamics of the SC order parameter, $\phi_{\rm SC}(t)$, for various excitation conditions. First we note that when $\omega_{0\rm X}$ becomes small with $g$ and $T(<T_{c})$ fixed, $\phi_{\rm SC}$ becomes larger because the system moves to the stronger coupling regime. This leads to the expectation that in the adiabatic regime, $\phi_{\rm SC}$ becomes larger with decreasing $\omega_{0\rm X}$. In other words, for $\Delta\omega_{0}>0$ (solid lines in Fig. 11) $\phi_{\rm SC}$ should initially decrease, while for $\Delta\omega_{0}<0$ (dashed lines) it should increase. However, in the initial dynamics for non-adiabatic excitations, the solid (dashed) lines show a positive (negative) hump, which is opposite to the above expectation. The hump becomes smaller with decreasing $\Omega$, and the dynamics becomes more consistent with the expectations from the adiabatic picture. We also note that, during the first period, one can see a significant transient enhancement of the order parameter in some cases, see Fig. 11. There the time-averaged $\phi_{\rm SC}(t)$ for one period can become larger than the equilibrium value if we carefully choose the phase of the excitation, i.e. the sign of $\Delta\omega_{0}$.

As time evolves further, $\phi_{\rm SC}(t)$ is however strongly suppressed from the initial value. One can see that the fastest suppression occurs around $\Omega=0.4\sim 2\omega_{\rm r}$, which is at the parametric resonance. This is consistent with the analysis of the NESS with the Floquet DMFT, where we have observed the strongest energy injection and the strongest increase of the effective temperature in this driving regime. The results in Fig. 11 show that even in the transient dynamics, the parametric resonance is associated with a rapid destruction of the SC order.

In Fig. 12, we show how the energies approach the NESS described by the Floquet DMFT. These results demonstrate that it takes several cycles for the system to reach the NESS.

To gain more insights into the dynamics of the order parameter, we compare it with the linear-response result. In the linear-response regime, the variation of the order parameter can be expressed as

[TABLE]

with $\chi_{B_{\bf 0},X^{2}}(t)\equiv-i\theta(t)\langle[B_{\bf 0}(t),X^{2}(0)]\rangle$, $B_{\bf 0}=\sum_{i}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}+H.c.$, and $F(t)=\Delta\omega_{0}\sin(\Omega t)$. The linear-response function can be numerically evaluated by considering, instead of the periodic driving, $H_{\rm ex}(t)=d_{f}\delta(t)\sum_{i}X_{i}^{2}$ and taking $d_{f}$ small enough. The resulting response function is shown in Fig. 13. The oscillation frequency in $\chi_{B_{\bf 0},X^{2}}(t)$ is $\omega\sim 0.4\sim 2\omega_{\rm r}$. 222With the phonon bath it turns out that the phonon renormalization and the effective phonon-phonon interaction discussed in Ref. Murakami2016, become weaker than without. Hence the oscillation frequency in $\chi_{B_{\bf 0},X^{2}}$ remains close to $2\omega_{\rm r}$ in the present case..

In Fig. 14, we compare the full dynamics (solid lines) and the linear-response component (dashed lines), which is obtained by taking the convolution, Eq. (54). The difference corresponds to the contributions from the non-linear components. The linear-response theory captures the initial behavior of the order parameter including the temporal enhancement or suppression. The initial hump whose direction is opposite to the adiabatic expectation is also captured. This hump originates from the positive value of $\chi_{B_{\bf 0},X^{2}}(t)$ just after $t=0$. Since the first derivative of $\chi_{B_{\bf 0},X^{2}}(t)$ at $t=0$ is zero, as can be analytically shown, this corresponds to the positive value of $\partial^{2}_{t}\chi_{B_{\bf 0},X^{2}}(t)|_{t=0}=8ig\omega_{0}\sum_{i}\langle(-c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}+c_{i\downarrow}c_{i\uparrow})X_{i}\rangle$. If we decrease $\Omega$ with $\Delta\omega_{0}$ fixed (as we approach the adiabatic regime), the contribution from the initial positive hump in $\chi_{B_{\bf 0},X^{2}}(t)$ decreases, and the hump in $\phi_{\rm SC}$ becomes smaller. At longer times, $\phi_{\rm SC}(t)$ starts to deviate from the linear-response result, which generally overestimates $\phi_{\rm SC}(t)$.

This indicates that the general effect of the non-linear component is to reduce the superconducting order. It is also illustrative to compare the results to the expectation in the adiabatic limit, in which the temperature is fixed to the bath temperature and the phonon potential is changed (see the dash-dotted curves in Fig. 14). (Note that in an isolated system, the adiabatic limit would be defined by constant entropy instead of constant temperature.) The comparison confirms that the dynamics approaches the adiabatic expectation as we decrease the excitation frequency.

IV.4 Superconducting fluctuations above $T_{c}$ under phonon driving

So far we have focused on the superconducting phase in the strong-coupling regime and showed that the phonon driving generally weakens the superconductivity and that this tendency is particularly pronounced in the parametric resonance regime. Now we move on to the weak-coupling regime ($\lambda\ll 1$), where the BCS theory has usually been applied. In this regime the transition temperature becomes very small and a direct investigation of the superconducting state is numerically difficult. Hence we study the time evolution of the superconducting fluctuations in the normal state, in systems with and without phonon driving. Specifically, we evaluate the dynamical pair susceptibility under the excitation described by Eq. (53). For this, as in the equilibrium case,Murakami2016 we add $H^{\prime}_{\rm ex}(t)=F_{\rm ex}(t)B_{\bf 0}$ on top of the phonon modulation, where $B_{\bf 0}=\sum_{i}(c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}+c_{i\downarrow}c_{i\uparrow})$ and $F_{\rm ex}(t)=d_{\rm f}\delta(t)$. We take $d_{\rm f}$ small enough and compute the time evolution of $B_{\bf 0}$ to obtain the dynamical pair susceptibility,

[TABLE]

In the normal state in equilibrium, this susceptibility decays as $t-t^{\prime}$ increases. The decay time increases as we approach the boundary to the superconducting state and diverges at the boundary. The interesting question thus is how this decay depends on the phonon driving in the weak-coupling regime.

In Fig. 15, we show the results for $g=0.3,\omega_{0}=0.4,\beta=320$ without heat baths. This condition corresponds to $\lambda\sim 0.23$ and $\omega_{\rm r}=0.32$, and the inverse transition temperature is $\beta_{c}\sim 1000$. By changing the temperature (down to $\beta=640$), $g$ and the couplings to the heat baths, we have confirmed that the behavior described below is the generic one in the weak-coupling regime. In Fig. 15(a)(c), we show the results around $\Omega=\omega_{\rm r}$. We remind the reader that from the analysis of $D_{0}$ in Sec. IV.1, an strong enhancement of $D^{R}(0)$ is expected for $\Omega\lesssim\omega_{\rm r}$. However, what we observe is that generally in the presence of periodic phonon modulations the lifetime of the fluctuations becomes shorter than without phonon modulations and that at $\Omega=\omega_{\rm r}$ the lifetime is particularly strongly reduced. This tendency becomes clearer with increasing $\Delta\omega_{0}$. In Fig. 16(a), we show the evolution of the energy absorption around the parametric resonance, $\Omega=\omega_{\rm r}$. It turns out that at the parametric resonance the energy absorption is also very efficient. These results suggest that even though there might be some enhancement of the attractive interaction, the energy absorption and resulting disturbance of long-range correlations is the dominant effect in a wide parameter range. In particular, at the parametric resonance, the maximized energy absorption rate leads to a particularly rapid heating of the system, and efficient destruction of the SC fluctuations. We also note that the decay of the fluctuations accelerates as $t$ increases, as is evidenced by the concave form of $\chi^{R}_{\rm pair}(t,0)$ in the log scale. This can also be interpreted as a heating effect. The more energy the system gains as $t$ increases, the faster the decay of the fluctuation becomes.

This dominant heating effect and quick destruction of the SC fluctuation can also be observed in another resonant regime $\Omega\sim 2\omega_{\rm r}$, where in the polaron picture an enhancement of the superconductivity is expected,Knap2016 see Fig. 15(b)(d). Here, $\Omega=0.65\sim 2\omega_{\rm r}$ shows the faster decay than the neighboring driving frequencies. The efficient absorption of energy resulting from this parametric phonon driving is shown in Fig. 16(b).

Now, we point out that there is a regime where the phonon driving slightly enhances the SC fluctuations. In Fig. 15(e)(f), we show the SC fluctuations for a smaller excitation frequency than those above, where less heating is expected, and for various excitation strengths. In all cases, $\chi^{R}_{\rm pair}(t,0)$ stays around the equilibrium curve, but exhibits a shift above or below the equilibrium value depending on the sign of $\Delta\omega_{0}$. For small amplitudes, the change of $\chi^{R}_{\rm pair}(t,0)$ with respect to the equilibrium value $\chi^{R}_{\rm pair,eq}(t)$ is linear in $\Delta\omega_{0}$. In order to identify the component linear in $\Delta\omega_{0}$, we rescale the difference of $\chi^{R}_{\rm pair}(t,0)-\chi^{R}_{\rm pair,eq}(t)$ at the smallest amplitude $\Delta\omega_{0}=0.02$ to larger amplitudes, and depict the results with dashed lines in Fig. 15(e). The linear components consistently explain the general behavior. The enhancement of the SC fluctuations for $\Delta\omega_{0}<0$ can be explained as follows. At $t\in[0,\pi/\Omega]\sim[0,20]$, the phonon frequency is temporarily reduced and the system is temporarily in a stronger coupling regime compared to equilibrium (this is not the effect of the Floquet sidebands). Hence the decay of the SC fluctuations is slower in this time interval. During $t\in[\pi/\Omega,2\pi/\Omega]$ the situation is opposite and we observe a faster decay. However, because of the enhancement during $t\in[0,\pi/\Omega]$, $\chi^{R}_{\rm pair}(t,0)$ stays above the equilibrium curve. The opposite is the case when $\Delta\omega_{0}$ is positive.

The non-linear component can result in a slight but systematic increase of the SC fluctuations, in particular at $\Delta\omega_{0}=\pm 0.06$, see the difference between the solid and dotted lines in Fig. 15(e). We note that both the heating and the creation of the Floquet sidebands of the phonons, which leads to the enhancement of the attractive interaction, start from the second order in the field strength. Therefore, this observation may indicate that the enhancement of the attractive interaction slightly dominates the heating effect. However this enhancement is very subtle and as we further increase the field strength this enhancement becomes weak and eventually the SC fluctuations decay faster than in the equilibrium case, see Fig. 15(f).

Finally we note that even though we have also studied a parameter set for $\omega_{0}>W$ in the weak coupling regime with various excitation conditions, we could not observe any enhancement of the SC fluctuation. 333As far as the system is in the weak-coupling regime the Migdal approximation may still work.

V Conclusions

We have systematically investigated the effects of the parametric phonon driving on the superconductivity (SC) using the simple Holstein model and the nonequilibrium dynamical mean-field theory (DMFT). In order to study the time-periodic nonequilibrium steady state, we extended the Floquet DMFT formalism, which has previously been applied to purely electronic systems, to this electron-phonon model. We have also studied the transient dynamics under periodic phonon driving with the nonequilibrium DMFT formulated on the Kadanoff-Baym contour. This simulation method, which does not involve a gradient approximation, can describe dynamics on arbitrarily fast timescales.

In the strong electron-phonon coupling regime, we have studied both the nonequilibrium steady states and the transient dynamics towards these states. Even in the presence of the renormalization from the electron-phonon coupling, the effective attractive interaction characterized by $\bar{D}^{R}(\omega=0)$ shows an enhancement in some parameter regime. However, taking into account heating and the nonthermal nature of the driven state, the ultimate effect of the driving turns out to be a suppression of the SC gap, the SC order parameter or the superfluid density both in the time-periodic steady states and the transient dynamics. The strongest suppression of SC is observed in the parametric resonant regime, $\Omega\sim 2\omega_{r}$.

In the weak-coupling regime, we have studied the SC fluctuations above $T_{c}$. We have found that, generally, they are suppressed by the phonon driving, in the sense that the pairing susceptibility decays more rapidly in the driven state. The decay becomes particularly fast in the parametric resonant regimes, which can be attributed to a particularly efficient energy absorption.

From these analyses, we conclude that the previously predicted enhancement of the superconductivity in the parametric resonant regime is not a general result. In the simple Holstein model with a phonon energy scale smaller than that of electrons, such effects, if any, are overwhelmed by the absorption of energy in a wide parameter range, and particularly strongly so near the parametric resonances. Even though our analysis focuses on a specific model, it demonstrates the importance of a careful treatment of heating effects and/or nonthermal distributions when discussing the potential effects of phonon driving. We note that in addition to the coupling of light to phonons, one could also take into account the direct coupling of light to electrons, by means of a Peierls substitution. This can be expected to lead to further heating effects.

As for the relation to the K3C60 experiment mentioned in the introduction,Mitrano2016 from our analysis it seems likely that the apparent enhancement of the SC originates from the static change of the HamiltonianKim2016 ; Mazza2017 or some other dynamical effectsKennes2017 ; Sentef2017 ; Nava2017 , which are not captured in the present simple model. Since alkali doped fullerides are multi-orbital systems and an effective negative Hund’s coupling from the Jahn-Teller screening plays an important role,Capone2009 ; Nomura2015 ; Karim2016 ; Hoshino2016 we may need to go beyond the Holstein model to fully understand the mechanism.

Finally, we remark that the diagrammatic method for steady states and the transient dynamics used here can be easily extended to other scenarios for the enhancement of SC, such as non-linear couplings to electrons,Kennes2017 ; Sentef2017 which is an interesting direction for future studies. The investigation of nonequilibrium states with an unconventional pairing such as the d-wave superconductivity in cuprates is also an interesting future work. In particular it is important to understand whether or not a particular mechanism for the enhancement of SC survives even if heating effects are taken into account.

Acknowledgements.

The authors wish to thank D. Gole$\check{\rm z}$, H. Strand and J. Okamoto for fruitful discussions. YM is supported by the Swiss National Science Foundation through NCCR Marvel. PW acknowledges support from FP7 ERC Starting Grant No. 278023. NT is supported by JSPS KAKENHI (Grant No. 16K17729). ME acknowledges support by the Deutsche Forschungsgemeinschaft within the Sonderforschungsbereich 925 (projects B4).

Appendix A The expression for the energy dissipation in the time-periodic NESS

The explicit expression for Eq. (39) is

[TABLE]

Using the bath’s fluctuation-dissipation relation, we obtain

[TABLE]

with $F_{f}(\omega)=\tanh(\frac{\beta\omega}{2})$ and

[TABLE]

with $F_{b}(\omega)=\coth(\frac{\beta\omega}{2})$.

These expressions imply that the energy dissipation $\bar{I}$ is related to the violation of the system’s fluctuation-dissipation relation in the time-periodic NESS. This is the Floquet generalization of the so-called Harada-Sasa relation.Harada2005 We have confirmed that this expression is indeed satisfied within the Floquet DMFT implemented with the self-consistent Migdal approximation.

Appendix B Results for the type 2 excitation

In this section, we show the results for the the type 2 excitation [Eq. (5)], whose effect has been discussed in Ref. Komnik2016, . It turns out that the general behavior is very similar to the results for the type 1 excitation [Eq. (4)] shown in the main text. The ultimate effect of the type 2 excitation is to generally suppress the superconductivity both in the strong- and weak-coupling regimes. This tendency becomes strong at $\Omega\simeq 2\omega_{r}$ and/or at $\Omega\simeq\omega_{r}$ as in the case of the type 1 excitation.

B.1 Nonequilibrium steady states

Here we show the results for the nonequilibrium steady states for the type 2 excitation in the strong electron-phonon coupling regime as in Sec. IV.2. Even in this regime, the type 2 excitation shows an increase of the effective attractive interaction characterized by $|g^{2}\bar{D}^{R}(0)|$ at $\Omega\lesssim\omega_{\rm r}$, and a decrease at $\Omega\gtrsim\omega_{\rm r}$. In Fig. 17, we show the time-averaged spectra and distribution functions for the electrons and phonons for $\Omega=0.125$, where we observed an increase of the attractive interaction. The general features are the same as in the case of the type 1 excitation (Fig. 6). The SC gap is suppressed with increasing driving amplitude and is completely wiped out at $\Delta\omega_{0}\geq 0.05$, see Fig. 17(a). The slope of the electron distribution functions [Fig. 17(b)] at $\omega=0$ decreases compared to the equilibrium case, which resembles the effect of heating. Precisely speaking, the equilibrium distribution functions at the effective temperature underestimate the number of the occupied states above the Fermi level. The characteristic features in the NESS with SC are the kink in $\bar{A}(\omega)$ and the hump in $\bar{f}(\omega)$ at $\omega=\Omega/2$. The renormalized phonon frequency increases as we increase the strength of the excitation, which is naturally expected from the equilibrium state with elevated temperatures, see Fig. 17(c). When we further increase the excitation strength, the Floquet sidebands of the phonons become more prominent [see $\Delta\omega=0.08$ data in Fig. 17(c)], in which case the spectrum does not resemble that of the equilibrium state with elevated temperatures.

The summary of the results for the type 2 excitation is shown in Fig. 18. Generally, in the presence of the phonon driving the SC order is weakened, and this tendency becomes strong around $\Omega=2\omega_{r}$ as in the type 1 case, see Fig. 18 (a). Again the order parameter evaluated from the effective temperature overestimates its size in the NESS, but its general behavior can be qualitatively reproduced. The (time-averaged) total energy and the effective temperature shows a peak around $\Omega=2\omega_{r}$, see Fig. 18(b). The peak position is shifted to higher frequency as the excitation strength is increased, which is attributed to the hardening of the renormalized phonon frequency as in the type 1 excitation. Quantitatively, under the type 2 excitation the SC order is more robust than under the type 1 excitation, which is consistent with the lower energy and the lower effective temperatures (compare Figs. 8 and 18).

B.2 Transient dynamics towards the NESS

In Fig. 19, we show the transient dynamics towards the NESS under the type 2 parametric phonon driving. It corresponds to Fig. 11 in the main part. One again finds that (i) the decease of the order parameter is fastest around $\Omega=0.4\simeq 2\omega_{r}$, (ii) the initial hump goes in the opposite direction from the adiabatic excitation, and becomes less prominent as we approach the adiabatic limit, and (iii) depending on the phase of the excitation, there can be an enhancement of the SC order parameter during the first cycle. Quantitatively speaking, the speed of the melting of SC is slightly slower than in the type 1 case at the same set of the parameters. This is consistent with the results for the NESS, where the order parameter is more robust compared to the type 1 and the effective temperature is lower than the type 1 case.

B.3 Superconducting fluctuations above $T_{c}$

Here we show the SC fluctuations above $T_{c}$ in the weak coupling regime and under the type 2 driving. Figure 20 shows the results corresponding to Fig. 15 in the main text. First we can again see that generally the decay of the fluctuations in the driven system is faster than in the absence of phonon driving, which indicates that the driving does not enhance the superconductivity, see Fig. 20 (a-d). As in the type 1 case, around $\Omega=\omega_{r}$, the decay becomes particularly fastest at $\Omega=\omega_{r}$, even though at $\Delta\Omega=-0.06$ the tendency is not so clear. Near $\Omega=2\omega_{r}$, there is again a point where the decay becomes particularly fastest. We also note the concave form of the $\chi_{\rm pair}(t,0)$ curve in the log scale, which is explained by the continuous heating of the system.

In Fig. 20 (e)(f), we show the results for an excitation frequency small compared to $\omega_{r}$. With the weak to moderate excitation strength, depending on the phase of $\Delta\omega_{0}$, $\chi_{\rm pair}(t,0)$ is shifted above or below the equilibrium curve, which is consistently explained by the linear component in $\Delta\omega_{0}$ (the dashed lines). The non-linear component slightly but systematically shifts the SC fluctuations above the linear component, which indicates that there is more than the heating effect. However, this slight enhancement is washed out when we further increase the excitation strength. These features are the same as in the type 1 case.