Non-equilibrium superconductivity in driven alkali-doped fullerides

TL;DR

This paper explores how periodic driving of interaction imbalances in alkali-doped fullerides can induce transient superconducting states that persist beyond the equilibrium critical temperature, shedding light on light-induced superconductivity.

Contribution

It demonstrates that dynamical modulation of interaction imbalance can create and sustain non-equilibrium superconducting states in fullerides, extending understanding of light-induced superconductivity.

Findings

Transient superconductivity can be induced above $T_c$ by periodic driving.

Driving can suppress or enhance superconductivity depending on frequency.

Interaction imbalance mechanism explains recent light-induced superconductivity experiments.

Abstract

We investigate the formation of non-equilibrium superconducting states in driven alkali-doped fullerides A$_3$C$_{60}$. Within a minimal three-orbital model for the superconductivity of these materials, it was recently demonstrated theoretically that an orbital-dependent imbalance of the interactions leads to an enhancement of superconductivity at equilibrium [M. Kim et al. Phys. Rev. B 94, 155152 (2016)]. We investigate the dynamical response to a time periodic modulation of this interaction imbalance, and show that it leads to the formation of a transient superconducting state which survives much beyond the equilibrium critical temperature $T_c$. For a specific range of modulation frequencies, we find that the driving reduces superconductivity when applied to a superconducting state below $T_c$, while still inducing a superconducting state when the initial temperature is larger than $T_c$. These findings reinforce the relevance of the interaction-imbalance mechanism as a possible explanation of the recent experimental observation of light-induced superconductivity in alkali-doped fullerenes.

Figures

Equations

$$\mathcal{H}=\sum _{\mathbf{k}\sigma }\sum _{a=1}^3ϵ(\mathbf{k})c_{\mathbf{k}a\sigma }^{†}{c}_{\mathbf{k}a\sigma }^{}+\sum _i{\mathcal{H}}_{loc}(i),$$

$$\begin{array}{rl}{\mathcal{H}}_{loc}& =\sum _{\alpha }{U}_{\alpha }{n}_{i\alpha \uparrow }{n}_{i\alpha \downarrow }\\ & +(U-2J_H)\sum _{\alpha \ne {\alpha }^{\prime }}{n}_{i\alpha \uparrow }{n}_{i{\alpha }^{\prime }\downarrow }\\ & +(U-3J_H)\sum _{\alpha <{\alpha }^{\prime }\sigma }{n}_{i\alpha \sigma }{n}_{i{\alpha }^{\prime }\sigma }\\ & +J_H\sum _{\alpha \ne {\alpha }^{\prime }}{c}_{i\alpha \uparrow }^{†}{c}_{i{\alpha }^{\prime }\downarrow }^{†}{c}_{i\alpha \downarrow }^{}{c}_{i{\alpha }^{\prime }\uparrow }^{}\\ & +J_H\sum _{\alpha \ne {\alpha }^{\prime }}{c}_{i\alpha \uparrow }^{†}{c}_{i\alpha \downarrow }^{†}{c}_{i{\alpha }^{\prime }\downarrow }^{}{c}_{i{\alpha }^{\prime }\uparrow }^{},\end{array}$$

$${\mathcal{H}}_{e-ph}\propto q(t{)}^2{n}_{\uparrow }{n}_{\downarrow },$$

$$\begin{array}{rl}{\mathcal{H}}_{e-ph}& =\sum _{\alpha }{C}_{\alpha }q(t{)}^2{n}_{\alpha \uparrow }{n}_{\alpha \downarrow }\\ & \equiv \sum _{\alpha }\Delta U_{\alpha }\left[1-\cos 2\Omega t\right]n_{\alpha \uparrow }{n}_{\alpha \downarrow },\end{array}$$

$$U_{x,y}(t)=U-r(t)\frac{\delta U}{2}\left(1-\cos 2\Omega t\right);U_z(t)=U,$$

$$|\Psi (t)⟩\simeq \prod _i{\mathcal{P}}_i(t)|{\Psi }_0(t)⟩,$$

$$\rho (t)=\mathcal{P}(t){\rho }_{\ast }(t)\mathcal{P}(t{)}^{†}$$

$$|\Psi (t)⟩\equiv \prod _i{\mathcal{P}}_i(t)|{\Psi }_0(t)⟩.$$

$$\delta \int ⟨\Psi (t)|i{\partial }_t-\mathcal{H}|\Psi (t)⟩=0.$$

$$⟨{\Psi }_0(t)|{\mathcal{P}}_i^{†}(t){\mathcal{P}}_i(t)|{\Psi }_0(t)⟩=1$$

$$⟨{\Psi }_0(t)|{\mathcal{P}}_i^{†}(t){\mathcal{P}}_i(t){\hat{\rho }}_i^{N,A}|{\Psi }_0(t)⟩=⟨{\Psi }_0(t)|{\hat{\rho }}_i^{N,S}|{\Psi }_0(t)⟩,$$

$${\left[{\hat{\rho }}_i^N\right]}_{\alpha ,\beta }=c_{i\alpha }^{†}{c}_{i\beta }^{}$$

$${\left[{\hat{\rho }}_i^A\right]}_{\alpha ,\beta }=c_{i\alpha }^{†}{c}_{i\beta }^{†}.$$

$${\mathcal{P}}_i=\sum _{\Gamma ,n}{\phi }_{i,\Gamma n}(t)|\Gamma ,i⟩⟨n,i|$$

$$⟨{\Psi }_0(t)|d_{i\alpha }^{†}{d}_{i\beta }^{}|{\Psi }_0(t)⟩={\delta }_{\alpha ,\beta }{n}_{i\alpha }^0(t)$$

$$⟨{\Psi }_0(t)|d_{i\alpha }^{†}{d}_{i\beta }^{†}|{\Psi }_0(t)⟩=0\forall \alpha ,\beta .$$

$${\phi }_{i,\Gamma n}(t)=\frac{{\Phi }_{i,\Gamma n}(t)}{\sqrt{P_{n,i}^0(t)}},$$

$$\begin{array}{rl}{P}_{n,i}^0(t)& =⟨{\Psi }_0(t)||n,i⟩⟨n,i||{\Psi }_0(t)⟩\\ & =\prod _{\alpha }{n}_{i,\alpha }^0(t{)}^{n_{\alpha }}(1-n_{i,\alpha }^0(t){)}^{1-n_{\alpha }}.\end{array}$$

$$\text{Tr}\left({\Phi }^{†}(t)\Phi (t)\right)=1$$

$$\text{Tr}\left({\Phi }^{†}(t)\Phi (t)d_{i\alpha }^{†}{d}_{i\beta }^{}\right)=⟨{\Psi }_0(t)|d_{i\alpha }^{†}{d}_{i\beta }^{}|{\Psi }_0(t)⟩={\delta }_{\alpha \beta }{n}_{\alpha }^0$$

$$\text{Tr}\left({\Phi }^{†}(t)\Phi (t)d_{i\alpha }^{†}{d}_{i\beta }^{†}\right)=⟨{\Psi }_0(t)|d_{i\alpha }^{†}{d}_{i\beta }^{†}|{\Psi }_0(t)⟩=0.$$

$$i{\partial }_t\hat{\Phi }(t)={\mathcal{H}}_{loc}(t)\hat{\Phi }(t)+⟨{\Psi }_0(t)|\frac{\delta {\mathcal{H}}_{\star }[\hat{\Phi }]}{\delta {\hat{\Phi }}^{†}(t)}|{\Psi }_0(t)⟩$$

$$i{\partial }_t|{\Psi }_0(t)⟩={\mathcal{H}}_{\star }[\hat{\Phi }]|{\Psi }_0(t)⟩.$$

$$c_{i\alpha }^{}\to \sum _{\beta }{R}_{\alpha \beta }[\Phi ]d_{i\beta }^{}+Q_{\alpha \beta }[\Phi ]d_{i\beta }^{†}.$$

$$R_{\alpha \beta }[\hat{\Phi }]=\frac{1}{\sqrt{n_{\beta }^0(t)(1-n_{\beta }^0(t))}}\text{Tr}\left({\hat{\Phi }}^{†}(t)c_{i\alpha }^{}\hat{\Phi }(t)d_{i\beta }^{†}\right)$$

$$Q_{\alpha \beta }[\hat{\Phi }]=\frac{1}{\sqrt{n_{\beta }^0(t)(1-n_{\beta }^0(t))}}\text{Tr}\left({\hat{\Phi }}^{†}(t)c_{i\alpha }^{}\hat{\Phi }(t)d_{i\beta }^{}\right)$$

$${\mathcal{H}}_{\star }=\sum _{\mathbf{k}}{\mathbf{\Psi }}_{\mathbf{k}}^{†}\left(\begin{array}{cc}{\mathbf{R}}^{†}{\hat{t}}_{\mathbf{k}}\mathbf{R}& {\mathbf{R}}^{†}{\hat{t}}_{\mathbf{k}}\mathbf{Q}\\ {\mathbf{Q}}^{†}{\hat{t}}_{\mathbf{k}}\mathbf{R}& {\mathbf{Q}}^{†}{\hat{t}}_{\mathbf{k}}\mathbf{Q}\end{array}\right){\mathbf{\Psi }}_{\mathbf{k}},$$

$${\mathcal{H}}_{loc}=\frac{U-3J_H}{2}\hat{N}\left(\hat{N}-1\right)-2J_H{\mathbf{S}}^2-\frac{J}{2}{\mathbf{L}}^2+\frac{5}{2}{J}_H\hat{N},$$

$$\text{Tr}\left({\hat{\Phi }}_k^{†}{\hat{\Phi }}_{k^{\prime }}^{}\right)={\delta }_{k,k^{\prime }}k,k^{\prime }=1,\dots ,N_{\Phi },$$

$$\hat{\Phi }(t)=\sum _{k=1}^{N_{\Phi }}{\alpha }_k(t){\hat{\Phi }}_k.$$

$${\left[{\hat{\rho }}_{\alpha \beta }^N\right]}_{k,k^{\prime }}=\text{Tr}\left({\hat{\Phi }}_k^{†}{\hat{\Phi }}_{k^{\prime }}{c}_{\alpha }^{†}{c}_{\beta }^{}\right)$$

$${\left[{\hat{\rho }}_{\alpha \beta }^S\right]}_{k,k^{\prime }}=\text{Tr}\left({\hat{\Phi }}_k^{†}{\hat{\Phi }}_{k^{\prime }}{c}_{\alpha }^{†}{c}_{\beta }^{†}\right)$$

$${\left[{\hat{r}}_{\alpha \beta }\right]}_{k,k^{\prime }}=\text{Tr}\left({\hat{\Phi }}_k^{†}{c}_{i\alpha }^{}{\hat{\Phi }}_{k^{\prime }}{d}_{i\beta }^{†}\right)$$

$${\left[{\hat{q}}_{\alpha \beta }\right]}_{k,k^{\prime }}=\text{Tr}\left({\hat{\Phi }}_k^{†}{c}_{i\alpha }^{}{\hat{\Phi }}_{k^{\prime }}{d}_{i\beta }^{}\right)$$

$${\left[{\hat{\mathcal{O}}}_i\right]}_{k,k^{\prime }}=\text{Tr}\left({\hat{\Phi }}_k^{†}{\mathcal{O}}_i{\hat{\Phi }}_{k^{\prime }}\right),$$

$$⟨\alpha (t)|\alpha (t)⟩=1$$

Full text

I Introduction

The optical stimulation of solids by means of strong light pulses has opened new routes for the investigation of collective phenomena in quantum materials [1]. A fascinating one consists in inducing superconductivity (SC) beyond the limits where it can be stabilized at equilibrium, which are set *e.g. *by temperature, external pressure or doping concentration. A series of experiments in different compounds revealed light-induced modifications of the electronic properties suggestive of the formation of a transient superconducting state extending above the equilibrium critical temperature $T_{c}$ [2, 3, 4, 5, 6]. Recently, the remarkable observation of a superconducting-like response above $T_{c}$ in the molecular compound K3C{}_{60}\leavevmode\nobreak\[6] enlarged this experimental panorama and raised new questions about the possible mechanisms leading to transient SC above $T_{c}$.

The possibility of a light-stimulated superconducting phase extending above the equilibrium critical temperature has been originally proposed by Eliashberg [7] who considered the quasiparticle redistribution induced by a laser excitation with frequencies below the equilibrium superconducting gap. The above experiments are far beyond this limit with excitation frequencies much larger than the superconducting gap, thus requiring the investigation of alternative mechanisms for light-induced superconductivity.

In the case of the K3C{}_{60}\leavevmode\nobreak\the laser frequencies for which the transient response is observed are close to the frequencies of four intra-molecular phonons, the $T_{1u}$ modes with frequencies in the mid-infrared range $60\lesssim\omega_{T_{1u}}\lesssim 180\leavevmode\nobreak\ \text{meV}$. The effect disappears for much larger excitation frequencies. This suggests that the observed effect might be related to the light-induced excitation of these phonon modes. [6, 8]

From the theoretical point of view, various mechanisms have been so far investigated such as the non-linear excitations of phononic modes [5, 9, 10] and their coupling to the electronic density [11, 12] or the effective slowing-down of the electronic motion [13, 14]. While all these mechanisms lead to an increase of the superconducting coupling which is expected to provide a source of transient SC in a broad class of superconductors, K3C{}_{60}\leavevmode\nobreak\appears as a peculiar case. Indeed, the absence of any enhancement or even the suppression of SC below $T_{c}$ reported in Ref. 6, together with the appearance of a transient response above $T_{c}$ is not fully understood within an effective SC coupling enhancement. Furthermore, SC in alkali-doped fullerides is strongly affected by the non-trivial interplay between pairing, electronic correlations and orbital degrees of freedom [15, 16, 17, 18, 19], requiring the investigation of the mechanisms for transient SC within a proper theoretical framework taking this interplay into account.

A first step in this direction has been recently taken by Kim et al. [8], working in the framework of the low-energy electronic description of fullerides based on the Jahn-Teller induced inversion of the effective Hund’s coupling [16], which provides one of the most successful description of the unconventional superconducting properties of these materials [17]. These authors demonstrated that a specific orbital-dependent perturbation of the on-site repulsive interactions does enhance SC at equilibrium. Such a perturbation was motivated by the possible modulations of the electronic interactions that result from the excitation of a local phononic mode, as already demonstrated for other correlated organic compounds [20, 21]. Furthermore, a first principle calculation [8] for K3C{}_{60}\leavevmode\nobreak\revealed that the favorable perturbation is indeed induced, under the assumption that light excites the $T_{1u}$ mode [6], as a result of the structural and electronic changes associated with this excitation.

This proposal relies entirely on equilibrium considerations however, and this raises the outstanding question of the relevance of this mechanism to the non-equilibrium response of the system and to the transient light-induced SC. In this work we address this question. We investigate the non-equilibrium dynamics induced by the time- and orbital- dependent modulation of the electron-electron repulsion. We show that this results in a transient superconducting state, which can be induced when the system is initially well above its equilibrium critical temperature. The properties of this transient state can be very dramatically different from the equilibrium expectation and depend on the frequency of the modulation. In particular, we uncover a regime of frequency in which the modulation leads to the reduction of SC below $T_{c}$ and to the creation of SC above $T_{c}$.

In the following we will first introduce the model and the non-equilibrium perturbation considered in this work together with the method used to describe the non-equilibrium dynamics. After presenting our results we will discuss possible implication for the description of the experimental observations.

II Model

The minimal description of strongly correlated superconductivity in alkali-doped fullerenes is given by the following three-bands model [16] arising from the $t_{1u}$ LUMO states of the C60 molecule half-filled with electrons donated by the alkali atoms

[TABLE]

where the local Hamiltonian $\mathcal{H}_{loc}(i)$ is of the Kanamori type [22] and takes into account intra- and inter- orbital electron electron repulsion, spin-flip and pair hopping terms, with a negative (inverted) Hund’s coupling $J_{H}$ resulting from the competition between the Hund’s coupling and the Jahn-Teller intramolecular interactions [16]. The explicit expression of the local Hamiltonian reads

[TABLE]

where $\alpha,\alpha^{\prime}$ and $\sigma$ indices indicates orbital and spin degrees of freedom respectively.

At equilibrium the interaction terms on each orbital are degenerate $U_{\alpha}=U$. The phonon excitation induces the modification of the local interactions energies $U_{\alpha}$ due to the coupling between the local electronic configurations and the coordinate of the displaced phononic mode along a given direction $q(t)=A\sin\Omega t$. In general, this is due to the fact that, for an odd parity mode, such as $T_{1u}$, at the lowest order the displaced phononic coordinate couples quadratically with the local double occupied states. [20, 21] For a single band case this leads to an additional term in the local Hamiltonian

[TABLE]

meaning the oscillation of the local interaction with a frequency $2\Omega$ around a renormalized value due to the fact that square of the mode displacement has a finite average $\langle{q(t)^{2}}\rangle\neq 0$. In the multi-band case, neglecting contributions coming from the coupling between different orbitals electronic configurations, Eq. (3) is generalized to

[TABLE]

where the coefficients $\Delta U_{\alpha}$ are specific properties of the phononic mode and of the direction of excitation, determined by the light pulse polarization. Here we take advantage of the first principle results of Ref. 8 showing that the displacement of the normal coordinate of the $T_{1u}$ mode along a given direction leads to the removal of the orbital degeneracy between the $U_{\alpha}$, leading to two orbitals with smaller interaction with respect to the third one. We insert this result in Eq. (4) by considering $\Delta U=-\left[\delta U,\delta U,0\right]$ with $\delta U>0$, so that the intra-orbital interaction terms in Eq. (2) become

[TABLE]

where $r(t)$ is a smooth ramping function, defined as $r(t)=1/2-3/4\cos\pi t/\tau+1/4\cos^{3}\pi t/\tau$ for $t<\tau$ and $r(t)=1$ for $t\geq\tau$, which phenomenologically takes into account the time $\tau$ during which the modulation of the $U$ is switched on.

We implemented and use the time-dependent Gutzwiller approximation (tdGA) [23, 24] extended to the multi-band superconducting case  [25, 26, 27]. The method is based on the variational ansatz for the time evolved state

[TABLE]

where $\left|\Psi_{0}(t)\right\rangle$ is an uncorrelated wavefunction describing the coherent quasiparticle dynamics and $\mathcal{P}(t)$ is a projector onto the local Hilbert spaces giving the weights of the local atomic multiplets. The dynamics of both quantities are determined via the time-dependent variational principle $\delta\int\langle\Psi(t)|i\partial_{t}-H|\Psi(t)\rangle=0$. At equilibrium the variational ansatz Eq. (6) is equivalent to the rotationally invariant slave bosons technique [28] which has been already successfully used to describe equilibrium strongly correlated SC in the present model [29].

The method is extended to the finite-temperature case by the introduction of a time dependent variational density matrix [30, 31]

[TABLE]

where the projector $\mathcal{P}(t)$ has the same definition as in (6) and $\rho_{*}(t)=\sum_{n}p_{n}|{\Psi_{n}(t)}\rangle\left\langle{\Psi_{n}(t)}\right|$ is the density matrix corresponding to a complete set of uncorrelated states $|{\Psi_{n}(t)}\rangle$ and a distribution $p_{n}$. The dynamical equations for (7) are obtained by applying the finite temperature generalization of the Dirac-Frenkel variational principle [32]. They are solved numerically, with an initial condition corresponding to the equilibrium thermal state. Details about the methods are reported in the Appendix.

III Results

In the following we will consider a semicircular density of states with a bandwidth $W=0.5$ eV and take $U=0.5\leavevmode\nobreak\ \text{eV}$ and $J_{H}=-0.02\leavevmode\nobreak\ \text{eV}$. In the inset of Fig. 1(c) we show that at equilibrium this corresponds to a superconductor on the weak correlation side of the superconducting dome determined by the electron-electron repulsion $U$ in the model Eq. (1). This is consistent with *e.g. *the pressure dependence of $T_{c}$ observed experimentally for K3C{}_{60}\leavevmode\nobreak\ [33]. We take the modulation frequency $\Omega$ as an adjustable parameter in a range reasonably including the typical frequencies of $T_{1u}$ modes and we fix $\delta U/U=0.1$. We choose a ramping time $\tau=0.9\leavevmode\nobreak\ \text{ps}$. The following results do not depend qualitatively on this choice.

In Fig. 1(a)-(b), we plot the dynamics of the orbital averaged amplitude of the order parameter $P_{sc}=\sum_{\alpha}P_{sc}^{\alpha}/3$, $P_{sc}^{\alpha}=|\langle{c^{\dagger}_{\alpha,\uparrow}c^{\dagger}_{\alpha,\downarrow}}\rangle|$, for two different driving frequencies and at two temperatures below and above the equilibrium critical temperature $T_{c}\simeq 42\leavevmode\nobreak\ \text{K}$ (T=0\leavevmode\nobreak\K in panel (a) and T=58\leavevmode\nobreak\K in panel (b)). A tiny symmetry breaking field is introduced for $T>T_{c}$ to allow SC to develop. 111The time needed for a finite order parameter to develop depends on this choise and, therefore, it is not related to any realistic time scale. On the other hand, the long time value of the order parameter is independent on this choice.

For the larger driving frequency ($\Omega=0.1875\leavevmode\nobreak\ \text{eV}$), we observe the increase of the superconducting order parameter at zero temperature and the formation of a finite order parameter for $T>T_{c}$. This establishes that the SC-enhancement mechanism based on the imbalance of $U$ does apply out of equilibrium. In particular, the formation of a finite order parameter above $T_{c}$ signals that the initial normal metal becomes an unstable state due to the increase of the critical temperature induced by the average perturbation. Both enhancement of the order parameter below $T_{c}$ and its formation for $T>T_{c}$ are expected from the equilibrium predictions of the average interaction imbalance and it is due to the energetic stabilization of the local multiplets with a singlet pair in the $x$ or $y$ orbitals [8].

We notice that, in principle, the continuous modulation of the interactions might lead to a continuos energy absorption inside the system which would eventually destroy the transient superconducting state. However, the persistence of the superconducting state as long as the interaction modulation is active suggests that no sizable energy absorption occours in this case. In the following we show that this strongly depends on the frequency of the interaction modulation.

The strong deviation from the above equilibrium expectations is observed when the driving frequency is lowered to $\Omega=0.15\leavevmode\nobreak\ \text{eV}$. Starting from the superconducting state, the order parameter undergoes a decrease instead of the expected increase. On the other hand, a finite order parameter is still established above the critical temperature, though its amplitude is smaller than the one established for $\Omega=0.1875\leavevmode\nobreak\ \text{eV}$.

In panel (c) we compare the transient order parameter as a function of the initial temperature to the equilibrium one. The former is extracted as a time average over a time interval $\Delta t=t_{max}-t_{min}$, $\overline{P}_{sc}=\int_{t_{min}}^{t_{max}}d\tau P_{sc}(\tau)/\Delta t$, for T=0\leavevmode\nobreak\K while we estimate it for $T>T_{c}$ from the value it takes at $t=3\leavevmode\nobreak\ \text{ps}$. For $\Omega=0.1875\leavevmode\nobreak\ \text{eV}$, the non-equilibrium perturbation enhances SC for all temperatures, up to about $T\lesssim\leavevmode\nobreak\ 100\leavevmode\nobreak\ \text{K}$, way above equilibrium $T_{c}$. On the contrary, the $\Omega=0.15\leavevmode\nobreak\ \text{eV}$ case shows a remarkable suppression of SC for $T<T_{c}$ and the formation of SC up to temperatures slightly below the previous case $T\lesssim\leavevmode\nobreak\ 90\leavevmode\nobreak\ \text{K}$.

The above results show that at the frequency $\Omega=0.15\leavevmode\nobreak\ \text{eV}$ the dynamical modulation of the inteaction leads to the suppression of the order parameter with respect to the value expected from the sole interaction imbalance. This strongly suggests that some energy absorption due to the continous modulation of the interaction occurs at this value of the driving frequency. In order to obtain insights into this, we study the orbital resolved dynamics of the zero-temperature superconductivity at three increasing frequencies (Fig. 2). We compare it to the time evolution obtained by an imbalance of $U$ equal to the average of the modulated case Eq. (5) switched on during the same ramp time $U_{x,y}(t)=U-r(t)\delta U/2$, hereafter called unmodulated dynamics.

The imbalance of $U$ lifts the orbital degeneracy between the components of the superconducting order parameter $P^{\alpha}_{sc}$. For the slowest driving frequency ($\Omega=0.0625\leavevmode\nobreak\ \text{eV}$ in panel (a)) the order parameter components display an oscillating behaviour with a fast component $\omega=2\Omega$ which reflects the periodical modulation of $U_{\alpha}$ and a slower component related to the amplitude of the superconducting gap. Only the slow component is retained in the case of the unmodulated dynamics (dashed lines).

In both cases a global enhancement of the order parameter is observed, though smaller with respect to what is expected from the equilibrium imbalanced case (see arrows in Fig. 2). This is understood from the fact that the average interaction imbalance is switched on in a finite time, wheares the equilibrium limit is expected to be recovered for an infinitely slow switching.

We notice that the modulated and unmodulated dynamics tend to separate at long times ($t\gtrsim 1.5\leavevmode\nobreak\ \text{ps}$) where a larger order parameter is established in the latter case. This suggests that the effect of the continous interaction modulation induces some energy absorption inside the system, leading to a slow decrease of the order parameter at long times. As already anticipated this effect is strongly dependent on the modulation frequency and it is almost absent at larger driving frequency ($\Omega=0.1875\leavevmode\nobreak\ \text{eV}$ in panel (c)), where the modulated and unmodulated dynamics become almost equivalent, the only difference being the fast oscillations of small amplitude in the modulated case.

A dramatic effect of the dynamical modulation occurs at intermediate frequencies. This is clearly seen on panel (b) of Fig.  2 for the driving frequency $\Omega=0.125\leavevmode\nobreak\ \text{eV}$ where the periodic modulation leads, in sharp contrast to the average interaction imbalance, to an almost complete suppression of the superconducting order already for times of the order of the ramping $\tau$.

The non-monotonic response of the transient state to the frequency of the interaction modulation is summarized in Fig. 3, where we report the existence of a frequency range $0.09\leavevmode\nobreak\ \text{eV}\lesssim\Omega\lesssim 0.16\leavevmode\nobreak\ \text{eV}$ in which SC is suppressed instead of enhanced. This frequency range is comprised between a lower frequency region in which the enhancement of SC is smaller than the case of the unmodulated dynamics and a higher frequency region in which the perturbation becomes completely anti-adiabatic with respect to the periodic change of $U$ and thus coincides with the unmodulated one.

We trace back the origin of the behaviour observed in the different regimes to the doublons excitations induced by the periodic modulation of $U$. Panels (b)-(d) of Fig. 3 display the dynamics of the double occupancies on each orbital $D_{\alpha}=\langle{n_{\alpha\uparrow}n_{\alpha\downarrow}}\rangle$ for the three frequencies representative of the different regimes. Due to the asymmetric value of the interaction the number of pairs is enhanced in the $(x,y)$ orbital and lowered on the $(z)$ one. However the dynamics is markedly different in the different regimes of frequency. At small and large frequencies (panel (b) and (d)), the dynamics closely follows the unmodulated one with superimposed $2\Omega$ oscillations indicating the creation of double occupancies on top of the orbitally imbalanced populations of non-equilibrium doublons. Such process becomes resonant in the intermediate frequency regime (panel (c)) where the strong amplification of the doublons oscillations is observed. As shown by the dynamics of the system’s internal energy $\mathcal{E}(t)=\langle{\mathcal{H}}\rangle$ (panel (e)), such a large amount of excitations leads to a sizable energy absorption which suppresses SC with respect to what is expected for the unmodulated dynamics (dashed line in panel (a)).

This shows that the energy absorption induced by the modulation frequency is responsible for the superconductivity suppression in the frequency region $0.09\lesssim\Omega\lesssim 0.16$ eV. We find that the origin of such behaviour is the resonance between the modulation frequency $2\Omega$ and the spectrum of the non-equilibrium excitations of the orbitally imbalanced populations of doublons induced by the asymmetric interaction (shaded area in panel (a)). Such spectrum is extracted from the frequency spectrum $\mathcal{F}_{D}$ of the dynamics of doublons following the sudden switch of a fixed $U-$imbalance, $U_{x,y}(t)=U-\theta(t)\delta U$ [35]. We compute $\mathcal{F}_{D}$ at fixed $\delta U$ and then integrate over a window equal to the amplitude of the interaction imbalance considered in the modulated dynamics $0<\delta U<0.05\leavevmode\nobreak\ \text{eV}$. The resulting spectrum has a broad three peaks structure which exactly matches the frequency region for which SC is suppressed.

An investigation of the structure of the spectrum revealed that the two side-bands mainly depend on the value of $J_{H}$ and disappear for $J_{H}=0$ indicating processes of inter-orbital origin. On the other hand, the central peak weakly depends on $J_{H}$ and decreases with $U$ (not shown) suggesting processes within the renormalized quasi-particle bandwidth as expected by the fact that the number of excited doublons is in a small quench regime ($\delta U/U=0.1$).

We finally observe that, as already anticipated in the discussion of Fig. 2, away from the resonance a small energy absorption is present for the slow driving case, whereas it is almost negligible for the fast one. This reflects in the mismatch between the modulated and the unmodulated dynamics for frequencies smaller than the resonance.

IV Discussion

The above results show that the dynamical modulation of the interaction that can be induced by the excitation of a molecular vibration may lead to significantly different effects as compared to equilibirum. This is due to energy absorption effects, that in the present case of interest is particularly evident in a range of driving frequency for which the modulation is resonant with the characteristic energies of the induced non-equilibrium excitations.

In such a region the effect of modulation competes with the effect of the interaction imbalance. The average imbalance favours the formation of a transient superconducting state with a larger order parameter extending much beyond the equilibrium critical temperature, while the time-dependent modulation induces energy absorption into the system. The latter effect may lead to a depletion of the initial superconducting state but, in spite of this, it does not preclude the SC order parameter to become finite above the equilibrium $T_{c}$ due to the former, even though with a smaller value with respect to the case in which the equivalence between the modulated and unmodulated dynamics is established. Therefore, the signature of a transient SC state above $T_{c}$ may survive also in the region where pairs are resonantly excited, as shown in Fig. 1 for $\Omega=0.15\leavevmode\nobreak\ \text{eV}$, thus realizing a non-trivial dynamical response for which SC is dynamically extended beyond $T_{c}$ and not enhanced for $T<T_{c}$.

In connection with the light-induced transient response in K3C60, we stress that our results are based on the assumption that the main effect of the light pulse is the excitation of the $T_{1u}$ phonon mode that in turn leads to the discussed interaction imbalance modulation. While this is not the only possible outcome of the light excitation and other types of excitation, *e.g. *of electronic origin [36], may play an important role, the above observations show that the discussed mechanism is a valid source for a transient superconductivity response above $T_{c}$ in alkali-doped fullerides.

Support to the considered mechanism comes from the facts that the superconducting-like transient response is observed for frequencies close to the phononic ones and for large enough laser fluence ($\gtrsim 1\leavevmode\nobreak\ \text{mJ cm}^{-2}$). Under the present assumption, the latter can be understood considering that a larger laser fluence translates into a larger displacement of the phonon mode and, therefore, to a larger amplitude of the interaction imbalance. As shown in Fig. 4 a larger interaction amplitude imbalance $\delta U/U=0.2$ leads to an enhancement of the discussed effects, correctly describing, on a qualitative level, the experimental dependence of the transient state on the laser fluence. In this respect, it should be noticed that the value of the interaction imbalance amplitude considered in this minimal model $\delta U/U=0.1$ might be an overestimation of the distorsion that can be induced in a more realistic description of fullerides [8].

An important aspect of our theoretical results is that, in a frequency regime, the dynamical response to the modulation can be different when the initial state is a superconducting state at $T<T_{c}$ and a normal state at $T>T_{c}$. With respect to the available experimental observations [6, 37] we mention that a similar effect has been observed in Ref. 6, where no enhancement of the transient superconducting gap below $T_{c}$ has been reported.

All the above elements emphasize the relevance of the present mechanism for the description of the transient response of stimulated fullerides. Despite this, additional evidence is needed in order to fully clarify the origin of such light-induced transient state. In this respect, experiments changing the alkali-atom, as e.g. the Mott insulating compound Cs3C60, might be useful to detect possible transient changes in the electronic correlations induced by the light pulse. At the same time an extensive frequency dependence investigation of the transient response below and above $T_{c}$ would be useful to better highlight the possible differences between the two cases, as described by our results.

In conclusion we have investigated the non-equilibrium dynamics of a minimal model describing SC in alkali-doped fullerides subject to the orbital asymmetric periodic modulation of the local interaction energies. This perturbation results from the assumption that the effect of the light-pulse is the excitation of a local vibrational mode and it is known to enhance SC at equilibrium [8]. We showed that this leads to the dynamical formation of a superconducting state extending beyond the equilibrium critical temperature $T_{c}$. Due to the correlated nature of such systems we showed that the transient state, while extending SC beyond $T_{c}$, may show no significant enhancement of the superconducting properties below $T_{c}$ where a SC suppression can even occour. This captures some non-trivial observations in the light-induced response of K3C{}_{60}\leavevmode\nobreak\which validates the mechanism as a possible source of light-induced SC in these systems.

V Acknowledgments.

We thank Minjae Kim, A. Cavalleri, M. Fabrizio, M. Capone and H. Strand for insightful discussions. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) ERC Grant Agreement nr. 319286 (Q-MAC).

Appendix: Gutzwiller approximation for multi-band superconductivity

In this section we give some details on the time dependent variational approach used to study the dynamics in the present model. For a datailed general formulation of the time-dependent Gutzwiller approximation we refer the reader to Ref. 24. In the following we will discuss the extension to the present multi-band superconducting case.

V.1 General Formulation

We consider the Hamiltonian defined in section II, which we divide into an hopping part $\mathcal{H}_{0}=\sum_{\langle i,j\rangle}\sum_{\alpha\beta}t_{i,j}^{\alpha,\beta}c^{\dagger}_{i\alpha}c^{\phantom{\dagger}}_{j\beta}$ and the local interaction $\mathcal{H}_{loc}$ defined in Eq. 2. Greek indeces include both orbital and spin degrees of freedom.

The dynamics is described starting from the following ansatz for the time-evolving wavefunction

[TABLE]

$\left|\Psi_{0}(t)\right\rangle$ is a single-particle wavefunction which is meant to describe the dynamics of delocalized quasiparticles and $\mathcal{P}_{i}(t)$ are operators acting on the local Hilbert space defined by the a of Fock states $|{\Gamma,i}\rangle$.

The variational dynamics is determined by the time-dependent variational principle

[TABLE]

An exact expression for the Lagrangian defining the above variational principle can be analytically obtained in the limit of infinite lattice coordination once the following constraints are imposed on the variational ansatz at each time $t$

[TABLE]

where $\hat{\rho}_{i}^{N,A}$ are the normal (N) and anomalous (A) components of the local single particle density matrix defined as

[TABLE]

A convenient representation of the local projectors is obtained in the so called mixed original-natural basis representation

[TABLE]

where $|{n,i}\rangle$ are the Fock states in the natural basis defined by a new set of creation and annhilation operators $d^{\dagger}_{i,\alpha}d^{\phantom{\dagger}}_{i,\alpha}$ for which the expectation values of the local single particle density matrix onto the uncorrelated wavefunction $|{\Psi_{0}(t)}\rangle$ is diagonal for the normal component and zero for the anomalous one

[TABLE]

The elements $\varphi_{i,\Gamma n}(t)$ are the set of local variational parameters defining the projectors $\mathcal{P}_{i}$ and they can be rewritten in the more convenient form

[TABLE]

where $P^{0}_{n,i}$ are the diagonal occupation probabilities onto the uncorrelated wavefunction of the local Fock states in the natural basis, expressed in terms of the variational density matrix $n_{i,\alpha}^{0}(t)$

[TABLE]

The matrices $\hat{\Phi}_{i}$, with elements $\Phi_{i,\Gamma n}$, contain the set of all the local variational parameters. In the present spatial homogeneous case all the $\hat{\Phi}_{i}$ matrices are equal and we neglect the site index $i$.

With the above definitions the constraints can be written in terms of the the local matrices $\hat{\Phi}$

[TABLE]

Plugging the above definitions into the expression of the Lagrangian $\mathcal{L}(t)=\langle{\Psi(t)}|{i\partial_{t}-\mathcal{H}}|{\Psi(t)}\rangle$ and imposing the condition $\delta\int_{0}^{t}d\tau\mathcal{L}(\tau)=0$ the following equations of motion are obtained [24]

[TABLE]

$\mathcal{H}_{loc}$ is the orginal time-dependent local Hamiltonian and $\mathcal{H}_{\star}[\hat{\Phi}]$ is a quadratic Hamiltonian which is obtained from the orginial tight binding Hamiltonian $\mathcal{H}_{0}$ upone the following transformation of the fermionic operators

[TABLE]

The transformation matrices $\mathbf{R}[\hat{\Phi}(t)]$ and $\mathbf{Q}[\hat{\Phi}(t)]$ entering in Eq. 24 depend on time through the local matrices $\hat{\Phi}(t)$ and their explicit expressions read

[TABLE]

With the above substitution the quasiparticle Hamiltonian $\mathcal{H}_{\star}$ acquires the general matrix form

[TABLE]

where $\boldsymbol{\Psi}_{\mathbf{k}}^{\dagger}=\left(d^{\dagger}_{\mathbf{k},1},\leavevmode\nobreak\ d^{\dagger}_{\mathbf{k},2},\leavevmode\nobreak\ \ldots,\leavevmode\nobreak\ d^{\dagger}_{\mathbf{k},N},\leavevmode\nobreak\ d^{\phantom{\dagger}}_{\mathbf{k},1},\leavevmode\nobreak\ d^{\phantom{\dagger}}_{\mathbf{k},2},\leavevmode\nobreak\ \ldots,\leavevmode\nobreak\ d^{\phantom{\dagger}}_{\mathbf{k},N}\right)$ and $\hat{t}_{\mathbf{k}}$ is the discrete Fourier transform of the hopping matrix elements in the Hamiltonian $\mathcal{H}_{0}.$

The coupled equations of motion 22-23 describe, respectively, the dynamics of the local degrees of freedom (Eq. 22) and of the delocalized quasiparticles (Eq. 23). The two dynamics are coupled in a mean-field way where each degree of freedom provides an effective field for the other, implying a mutual feedback between the localized and delocalized degrees of freedom. This is a big improvement of the present method with respect to standard mean-field (Hartree-Fock) approaches where only the delocalized quasiparticles are described. For this reason the method is able to capture correlations effects beyond mean-field, as the strongly correlated superconductivity discussed in the present case.

V.2 An explicit implementation

We now explicity illustrate the implementation of the above general formulation for the three orbitals model discussed in the main text in the presence of superconductivity.

The variational wavefunction depends on the choise of the variational matrix $\hat{\Phi}$ defined in the mixed basis of original and natural Fock states. In the three orbitals case the dimension of the local Hilber spase is 64, so that a total number of 64$\times$64=4096 variational parameters is contained in $\hat{\Phi}$. Such a large number of variational parameter is conveniently reduced exploiting the symmetries of the Hamiltonian as detailed described in Ref. 25.

The orginal Hamiltonian has a $U(1)\times SU(2)\times O(3)$ symmetry corresponding to charge conservation and invariance with respect to spin and orbital rotations respectively. This can be readly appreciated writing the local Hamiltonian in terms of the charge operators $N$ and the spin $\mathbf{S}$ and orbital $\mathbf{L}$ rotations generators [38]

[TABLE]

where, restoring the the separation between orbital and spin indeces, $\hat{N}=\sum_{a\sigma}c^{\dagger}_{a\sigma}c^{\phantom{\dagger}}_{a\sigma}$, $\mathbf{S}=\frac{1}{2}\sum_{a}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{a\sigma}\boldsymbol{\tau}_{\sigma\sigma^{\prime}}c^{\phantom{\dagger}}_{a\sigma}$ and $L_{a}=i\sum_{bc}\sum_{\sigma}\epsilon_{abc}c^{\dagger}_{b\sigma}c^{\phantom{\dagger}}_{c\sigma}$, $\boldsymbol{\tau}$ are Pauli matrices and $\epsilon_{abc}$ is the Levi-Civita tensor.

In order to describe superconductivity we break the $U(1)$ symmetry which corresponds to allow for non-zero $\hat{\Phi}$ matrix elements between states with different number of particles. In such a situation, the matrix $\mathbf{Q}$ defined in Eq. 26 is non zero, so that the quasiparticle Hamiltonian becomes an effective multi-band BCS Hamiltonian. The removal of the orbital degeneracy introduced by the imbalanced interaction terms explicitly breaks the $O(3)$ symmetry so that only the $SU(2)$ symmetry is left. We impose the $SU(2)$ symmetry onto the $\hat{\Phi}$ matrix following the method outlined in Ref. 25. This leads to $N_{\Phi}=429$ independent variational parameters.

The $\hat{\Phi}$ matrix can be expanded onto a basis of $N_{\Phi}$ matrices satifying

[TABLE]

so that all the information about the variational parameters can be stored in a time-dependent vector of complex number $|{\alpha(t)}\rangle$ of dimension $N_{\Phi}$ [25]

[TABLE]

Next we define the matrices of dimension $N_{\Phi}\times N_{\Phi}$ containg the possible combinations of traces over the matrix basis $\hat{\Phi}_{k}$

[TABLE]

where $\mathcal{O}_{i}$ represents any local many-body operators.

With the above definitions the constraints and the hopping renormalization matrices $\mathbf{R}$ and $\mathbf{Q}$ become

[TABLE]

Inserting the above definitions in the evolution of the $\hat{\Phi}$ matrix Eq. 22, we obtain the time evolution in terms of the $|{\alpha(t)}\rangle$ vector

[TABLE]

$\hat{\mathcal{H}}_{loc}(t)$ is $N_{\Phi}\times N_{\Phi}$ Hamiltonian defined as follow

[TABLE]

The matrices $\mathbf{\mathcal{D}}$ and $\mathbf{\mathcal{S}}$ are defined through the occupations onto the quasiparticle wavefunction $|{\Psi_{0}(t)}\rangle$

[TABLE]

[TABLE]

The dynamics described by Eq. 41 is coupled to the dynamics of the wavefunction $|{\Psi_{0}(t)}\rangle$ (Eq. 23) which can be fully expressed through the dynamics of the occupations $\Delta_{\mathbf{k}}^{\alpha,\beta}(t)=\langle{\Psi_{0}(t)}|{d^{\dagger}_{\mathbf{k}\alpha}d^{\phantom{\dagger}}_{\mathbf{k}\beta}}|{\Psi_{0}(t)}\rangle$ and $\Gamma_{\mathbf{k}}^{\alpha,\beta}(t)=\langle{\Psi_{0}(t)}|{d^{\dagger}_{\mathbf{k}\alpha}d^{\dagger}_{\mathbf{k}\beta}}|{\Psi_{0}(t)}\rangle$

[TABLE]

The system of coupled Eqs. 41, 45 and 46 fully describe the dynamics within the variational ansatz Eq. 8 in a lattice with infinite coordination number.

The number of coupled differential equations is further reduced by the fact that in the present case the hopping matrix is diagonal in both spin and orbital indeces and we consider only the possibility of intra-orbital pairing. This means that $\Delta_{\mathbf{k}}^{\alpha,\beta}=\delta_{\alpha,\beta}\Delta_{\mathbf{k}}^{\alpha}$ and $\Gamma_{\mathbf{k}}^{\alpha,\beta}=\delta_{a,b}(1-\delta_{\sigma\sigma^{\prime}})\Gamma_{\mathbf{k}}^{a\sigma,a\sigma^{\prime}}$, where in the last case we separated orbital $(a,b)$ and spin $(\sigma,\sigma^{\prime})$ degrees of freedom. Moreover, the matrix $\mathbf{R}$ is diagonal in both orbital and spin indeces, while the matrix $\mathbf{Q}$ is diagonal in the orbital index and couples only states with opposite spin.

The dynamics is unitary and preserves the normalization constraints Eq. 36 and the density constraint Eq. 37 and they need to be enforced only at equilibrium. On the contrary, we found that the constraint on the anomalous density is no more conserved by the unitary dynamics as the particle-hole symmetry is lifted by the imbalanced interaction. In that case we introduce a set of time-dependent Lagrange multipliers enforcing the constraints at each time step.

The equations of motions are solved using the explicit 4th-order Runge-Kutta method with a time discretization $\delta t=0.01$ starting from the variational estimation of the equilibrium ground state. The latter is obtained from the stationary limit of the equations of motion

[TABLE]

This corresponds to finding the ground state of a non-linear eigenvalue problem with appropriate Lagrange parameters enforcing the constraints at $t=0$. This can be done recursively at fixed value of the variational density matrix $n^{0}$. [25] Therefore, a full minimization of the obtained ground state energy with respect to $n^{0}$ gives the variational estimation of the ground state energy.

V.3 FINITE TEMPERATURE

The finite temperature extension closely follows the zero-temperature one where the variational ansatz for the time evolving ground state is replaced by the ansatz for the time dependent density matrix

[TABLE]

namely the Gutzwiller projection is done on each uncorrelated state $|{\Psi_{0,n}(t)}\rangle$ describing an uncorrelated density matrix $\rho_{\star}(t)=\sum_{n}p_{n}|{\Psi_{0,n}(t)}\rangle\left\langle{\Psi_{0,n}(t)}\right|$ through the distribution $p_{n}$. The Dirac-Frenkel extension of the time-dependent variational principle reads

[TABLE]

With the same definitions of projectors and constraints as in the previous section the equation of motions are equivalent to the zero-temperature case with the averages onto the uncorrelated wavefunction replaced by traces over the uncorrelated density matrix

[TABLE]

A practical solution of Eqs. 51-52 is obtained starting the initial estimation of the finite temperature expectation values for the quasiparticle occupations $\Delta_{\mathbf{k}}^{\alpha,\beta}(t)=\text{Tr}\left(\rho_{\star}(t)d^{\dagger}_{\mathbf{k}\alpha}d^{\phantom{\dagger}}_{\mathbf{k}\beta}\right)$ and $\Gamma_{\mathbf{k}}^{\alpha,\beta}(t)=\text{Tr}\left(\rho_{\star}(t)d^{\dagger}_{\mathbf{k}\alpha}d^{\dagger}_{\mathbf{k}\beta}\right)$ leading to the following equations of motions

[TABLE]

which can be written only in terms of combinations of $\Delta_{\mathbf{k}}^{\alpha,\beta}$ and $\Gamma_{\mathbf{k}}^{\alpha,\beta}$. Therefore the explicit time evolution for each $|{\Psi_{0,n}(t)}\rangle$ Eq. 52 is, in practice, not needed and the distribution $p_{n}$ is defined once for all by the initial variational estimation of the thermal state.

The initial thermal state is obtained by minimizaing the variational estimation of the free energy. [30] In doing so, for a simplification of the minimization procedure, we neglect the contribution coming from the entropy of the local degrees of freedom. This is a reasonable approximation since we focus on the weak coupling side of the superconducting dome far from the Mott transition, where the contribution of the entropy of the local degrees of freedom is small. Its inclusion would lead to a weak renormalization of the equilibrium transition temperature, but no qualitative difference in the dynamics is expected with respect to what reported in the main text.