Transport and spectral signatures of transient fluctuating superfluids in the absence of long-range order

TL;DR

This paper investigates the non-equilibrium dynamics of transient fluctuating superfluids without long-range order, deriving a quantum kinetic framework to analyze optical conductivity and density of states evolution after a quench.

Contribution

It introduces a quantum kinetic equation for fluctuating superfluids in non-equilibrium, accounting for conservation laws and fluctuation effects, with applications to optical conductivity and density of states.

Findings

Enhanced low-frequency density of states due to Andreev scattering

Density of states growth near the critical point after a quench

Derived kinetic equations linking fluctuations to transport properties

Abstract

Results are presented for the quench dynamics of a clean and interacting electron system, where the quench involves varying the strength of the attractive interaction along arbitrary quench trajectories. The initial state before the quench is assumed to be a normal electron gas, and the dynamics is studied in a regime where long-range order is absent, but nonequilibrium superconducting fluctuations need to be accounted for. A quantum kinetic equation using a two-particle irreducible formalism is derived. Conservation of energy, particle-number, and momentum emerge naturally, with the conserved currents depending on both the electron Green's functions and the Green's functions for the superconducting fluctuations. The quantum kinetic equation is employed to derive a kinetic equation for the current, and the transient optical conductivity relevant to pump-probe spectroscopy is studied. The general structure of the kinetic equation for the current is also justified by a phenomenological approach that assumes model-F in the Halperin-Hohenberg classification, corresponding to a non-conserved order-parameter coupled to a conserved density. Linear response conductivity and the diffusion coefficients in thermal equilibrium are also derived, and connections with Aslamazov-Larkin fluctuation corrections are highlighted. Results are also presented for the time-evolution of the local density of states. It is shown that Andreev scattering processes result in an enhanced density of states at low frequencies. For a quench trajectory corresponding to a sudden quench to the critical point, the density of states is shown to grow in a manner where the time after the quench plays the role of the inverse detuning from the critical point.

Equations

$$H(t)$$

$$H_0$$

$$H_{\mathrm{int}}$$

$$\displaystyle\partial_{t}j^{i}(t)-E^{i}(t)=-\frac{1}{\tau_{r}}\biggl{[}A(t)j^{i}(t)$$

$$\displaystyle-\alpha\int_{0}^{t}dt^{\prime}\biggl{\{}B(t,t^{\prime})j^{i}(t^{\prime})+C(t,t^{\prime})E^{i}(t^{\prime})\biggr{\}}\biggr{]}.$$

$$\mathcal{F}\left[\Delta ,\rho \right]\equiv \int d^{d}r[e\Phi \delta \rho +\frac{1}{2C}(\delta \rho {)}^2+r(t)|\Delta {|}^2$$

$$+\frac{1}{2M}|(i\vec{\nabla }-2e\vec{A})\Delta {|}^2+u|\Delta {|}^4+{\gamma }_0(\delta \rho )|\Delta {|}^2].$$

$$\left(\frac{\partial \Delta }{\partial t}-\left\{\Delta ,\rho \right\}\frac{\partial \mathcal{F}}{\partial \rho }\right)$$

$$\displaystyle\left(\frac{\partial\Delta}{\partial t}+2i\Delta\left[e\Phi+\frac{\delta\rho}{C}+\gamma_{0}|\Delta|^{2}\right]\right)=-\Gamma_{\Delta}\biggl{[}r(t)+\gamma_{0}\delta\rho$$

$$\displaystyle+2u|\Delta|^{2}+\frac{1}{2M}(i\nabla^{i}-2eA^{i})^{2}\biggr{]}\Delta+\xi_{\Delta},$$

$$\frac{\partial \delta \rho }{\partial t}$$

$${\vec{j}}_{\mathrm{dis}}$$

$$j_{\mathrm{dis}}^i={\Gamma }_{j_{\mathrm{dis}}}{E}^i.$$

$$\left\{\rho ,\Delta \right\}\frac{\partial \mathcal{F}}{\partial \Delta }$$

$$\displaystyle-\Delta^{*}(i\nabla^{i}-2eA^{i})^{2}\Delta\biggr{]}=\vec{\nabla}\cdot\vec{j}_{\rm sc},$$

$${\vec{j}}_{\mathrm{sc}}$$

$$\frac{\partial \delta \rho }{\partial t}=\vec{\nabla }\cdot {\vec{j}}_{\mathrm{sc}}+\vec{\nabla }\cdot {\vec{j}}_{\mathrm{dis}}.$$

$${\vec{j}}_{\mathrm{sc}}$$

$$=\frac{2}{M}\int \frac{d^{d}q}{(2\pi {)}^d}\vec{q}⟨{\Delta }_q{\Delta }_q^{\ast }⟩$$

$$=2\int \frac{d^{d}q}{(2\pi {)}^d}\frac{\vec{q}}{M}F(q,t),$$

$$\propto \int d^{d}q\frac{\vec{q}}{M}{\partial }_{t}F(q,t),$$

$${\partial }_t{j}_{\mathrm{dis}}^i-E^i(t)=-{\tau }_r^{-1}A(t)j_{\mathrm{dis}}^i.$$

$$\frac{\partial {\Delta }_q}{\partial t}$$

$$+{\xi }_{\Delta }.$$

$$\delta {\rho }_q(t)=c{j}^i(t)q^i,$$

$${\lambda }_q(t)={\Gamma }_{\Delta }\left[r(t)+\frac{q^2}{2M}\right]={\Gamma }_{\Delta }{ϵ}_q(t),$$

$${\Delta }_q(t)={\int }_0^{t}d{t}^{\prime }{\xi }_{\Delta }(t^{\prime })e^{-{\int }_{t^{\prime }}^{t}d{t}^{\prime \prime }{\lambda }_q(t^{\prime \prime })}.$$

$$F(q,E=0,t)=⟨{\Delta }_q{\Delta }_q^{\ast }⟩=2{\Gamma }_{\Delta }T{\int }_0^{t}d{t}^{\prime }{e}^{-2{\int }_{t^{\prime }}^{t}d{t}^{\prime \prime }{\lambda }_q(t^{\prime \prime })}.$$

$$\delta {\lambda }_q(s)=\frac{{\Gamma }_{\Delta }}{M}(2e)\vec{q}\cdot {\int }_s^{t}d{t}^{\prime }\vec{E}(t^{\prime })+{\Gamma }_{\Delta }{\gamma }_0{q}^i{j}^i(s),$$

$$\displaystyle F(q,t)\simeq 2\Gamma_{\Delta}T\int_{0}^{t}dt^{\prime}e^{-2\int_{t^{\prime}}^{t}dt^{\prime\prime}\lambda_{q}(t^{\prime\prime})}\biggl{[}1-2\int_{t^{\prime}}^{t}du\delta\lambda_{q}(u)\biggr{]}.$$

$$F(q,t)=F(q,E=0,t)+\delta F(q,t),$$

$$\delta F(q,t)=-2{\int }_0^{t}du\delta {\lambda }_q(u)e^{-2{\int }_u^{t}d{t}^{\prime \prime }{\lambda }_q(t^{\prime \prime })}F(q,E=0,u).$$

$$\delta F(q,t)=-2F(q){\int }_0^{t}du\delta {\lambda }_q(u)e^{-2(t-u){\lambda }_q}.$$

$$\delta F_q={\Gamma }_{\Delta }^{-1}\frac{T}{{ϵ}_q^3}\frac{\vec{q}}{M}\cdot e\vec{E}.$$

$$j^i=2e\sum _q\frac{q^i}{M}\delta F(q)={\Gamma }_{\Delta }^{-1}2e^2\sum _q\frac{T}{{ϵ}_q^3}\frac{q^i{q}^j}{M^2}{E}^j,$$

Full text

I Introduction

Recent years have seen impressive advances in ultrafast measurements of strongly correlated systems Smallwood et al. (2012); Zhang and Averitt (2014); Smallwood et al. (2014); Wall et al. (2018); Harter et al. (2018); McIver et al. (shed); Mitrano et al. (2018); Tengdin et al. (2018). In these experiments, a pump field strongly perturbs the system, while a weak probe field studies the eventual time-evolution over time scales that can range from femto-seconds to nano-seconds. Thus one can study the full nonequilibrium dynamics of the system from short times, to its eventual thermalization at long times. Moreover, the access to probe fields ranging from x-rays to mid-infrared allows one to probe dynamics from short length and time scales to longer length and time scales associated with collective modes.

Among the various pump-probe studies, notable examples are experiments that show the appearance of highly conducting, superconducting-like states, that can persist from a few to hundreds of pico-seconds Fausti et al. (2011); Mitrano et al. (2016); Pomarico et al. (2017); Cremin et al. (shed). The microscopic origin of this physics is not fully understood. Explanations range from destabilization by the pump laser of a competing order in favor of superconductivity Okamoto et al. (2016); Rajasekaran et al. (2018), to selective pumping of phonon modes that can engineer the effective electronic degrees of freedom so as to enhance the attractive Hubbard-U, and hence $T_{c}$ Knap et al. (2016); Sentef et al. (2016); Sentef (2017); Kennes et al. (2017a); Murakami et al. (2017). There are also theoretical studies that argue that the metastable highly conducting state observed in experiments may not have a superconducting origin Chiriacò et al. (2018).

The resulting superconducting-like states are transient in nature, where traditional measurements of superconductivity, such as the Meissner effect Tinkham (1996), cannot be performed. On the other hand, time-resolved transport and angle-resolved photo-emission spectroscopy (tr-ARPES) are more relevant probes for such transient states. Thus, a microscopic treatment that makes predictions for transient conductivity and spectral features of a highly nonequilibrium system, accounting for dynamics of collective modes, is needed. This is the goal of the paper.

In what follows, we model the pump field as a quantum quench Calabrese and Cardy (2006); Mitra (2018) in that it perturbs a microscopic parameter of the system, such as the strength of the attractive Hubbard-U. We assume that, as in the experiments Mitrano et al. (2016), the electronic system is initially in the normal phase. The quench involves tuning the magnitude of the attractive Hubbard-U along different trajectories. How superconducting fluctuations develop in time, and how they affect transient transport and spectral properties are studied. Denoting the superconducting order-parameter as $\Delta$, we assume that all throughout the dynamics, the system never develops true long range order, $\langle\Delta(t)\rangle=0$. Thus our approach is complementary to other studies where the starting point is usually a fully ordered superconducting state, and the mean-field dynamics of the order-parameter $\langle\Delta(t)\rangle$ is probed Yuzbashyan et al. (2006); Barankov and Levitov (2006); Peronaci et al. (2015); Yuzbashyan et al. (2015); Chou et al. (2017); Scaramazza et al. (2019). In our case, the non-trivial dynamics appears in the fluctuations $\langle\Delta(t)\Delta(t^{\prime})\rangle$ which we study using a two particle irreducible (2PI) formalism. We justify the selection of diagrams using a $1/N$ approach where $N$ denotes an orbital degree of freedom. The resulting choice of diagrams is equivalent to a self-consistent random-phase-approximation (RPA) in the particle-particle channel. The physical meaning of the RPA is that the superconducting fluctuations are weakly interacting with each other.

While in this paper we study a clean system, a complementary study of the transient optical conductivity of a disordered system appears in Ref. Lemonik and Mitra, 2018a, where a Kubo formalism approach was employed, and the quench dynamics arising from fluctuation corrections of the Aslamazov-Larkin (AL) Azlamazov and Larkin (1968a, b) and Maki-Thompson (MT) Maki (1968a, b); Thompson (1970) type was explored. We avoid a Kubo-formalism approach in this paper because conserving approximations are harder to make, as compared to a quantum kinetic equation approach. This is an observation well known from other studies for transport in thermal equilibrium performed for example in the particle-hole channel Catelani and Aleiner (2005). Since we are in a disorder-free system, our quantum kinetic equations conserve momentum. Thus, to obtain any non-trivial current dynamics, we have to break Galilean invariance. We do this through the underlying lattice dispersion, while neglecting Umpklapp processes.

Complementary studies of transient conductivity involve phenomenological time-dependent Ginzburg-Landau (TDGL) theories Kennes and Millis (2017); Iwazaki et al. (shed), and microscopic approaches based on a Bardeen-Cooper-Schrieffer (BCS) mean-field treatment of the order-parameter dynamics Kennes et al. (2017b). A microscopic mean field approach has also been used to study transient spectral densities Xu et al. (2019), while nonequilibrium dynamical mean field theory has been used to study transient spectral properties of a BCS superconductor Stahl and Eckstein (shed).

The outline of the paper is as follows. In Section II we present the model. In this section we also outline a qualitative analysis for the transient conductivity employing a classical Model-F Hohenberg and Halperin (1977) scenario for the fluctuation dynamics. In Section III we derive the quantum kinetic equations using a two particle irreducible (2PI) formalism, and outline how the transient dynamics respects all the conservation laws (particle, energy, momentum) of the system. In Section IV we simplify the 2PI quantum kinetic equation to that of an effective classical equation for the decay of the current employing a separation of time-scales. In the process we justify the model-F scenario discussed in Section II. In Section V we assume thermal equilibrium, and derive the conductivity and diffusion constants. We discuss these in the context of Aslamazov-Larkin (AL) and Maki-Thompson (MT) corrections Larkin and Varlamov (2002). In Section VI we present results for the transient conductivity for some representative quench profiles. In section VII we present results for the transient local density of states, and finally in Section VIII we present our conclusions. Many intermediate steps of the derivations are relegated to the appendices.

II Model

The Hamiltonian describes fermions with spin and an additional orbital degree on a regular lattice in $d$ spatial dimensions. Fermions are created by the operator $c^{\dagger}_{\sigma n}(r_{i})$, where $\sigma=\pm 1$ labels spin, $n=1\ldots N$ labels the orbital index and $r_{i}$ gives the lattice site. These have the Fourier transform $\tilde{c}_{\sigma n}(k)=\sum_{r_{i}}e^{i\vec{k}\cdot\vec{r}_{i}}c_{\sigma n}(r_{i})$. The dynamics of the fermions is governed by the Hamiltonian:

[TABLE]

where $\epsilon_{k}$ is the dispersion, and $\mu$ is the chemical potential. We will assume for simplicity that at this chemical potential there is only a single Fermi surface. We explicitly allow for the interaction $U(t)$ to vary with time.

In addition, to probe current dynamics, we introduce an electric field $\vec{E}=-\partial_{t}\vec{A}$ by minimal coupling $\vec{k}\rightarrow\vec{k}-q\vec{A}$. Accounting for the charge of the electron $q=-e$, this appears in the Hamiltonian via $\vec{k}\rightarrow\vec{k}+e\vec{A}$, or equivalently, $\vec{\nabla}/i\rightarrow\vec{\nabla}/i+e\vec{A}$. In the microscopic derivation we will set $\hbar=1$.

In thermal equilibrium at a temperature $T$, and for spatial dimension $d>2$, the above system becomes superconducting below a critical temperature $T<T_{C}$. At $d=2$, there is a Berezenskii-Kosterlitz-Thouless (BKT) transition at $T<T_{\rm BKT}$ where the system shows only quasi-long range order. For $T>T_{\rm BKT}$, although the system has no long range order, yet superconducting fluctuations play a key role in transport. The lower the spatial dimensions, the more important the role of fluctuations. In particular for spatial dimension $d=2$, well known results for the optical conductivity exist for a strongly disordered system Larkin and Varlamov (2002). While our derivation is valid in any spatial dimension, we will present results for spatial dimensions $d=2$ where fluctuation effects are most pronounced. The experiments also study either a layered superconducting system, or three dimensional systems where the pump field is effectively a surface perturbation. This makes the spatial dimension $d=2$ also experimentally relevant.

Before we go into details, we note that our microscopic treatment leads to the following expression (reported in Ref. Lemonik and Mitra, 2018b) for the dynamics of the current $j^{i}$, generated by an applied electric field $E^{i}$,

[TABLE]

Above $\tau_{r},A,B,C$ are kinetic coefficients. We have set $e=1,\rho/\bar{m}=1$, where $\rho$ is the electron density, and $\bar{m}$ is the electron effective mass. When $B=C=0$, the above equation simply implies a time-dependent Drude scattering rate $\tau_{r}^{-1}A(t)$, which in our model arises because electrons scatter off superconducting fluctuations whose density is changing with time due to the quench. The appearance of memory terms $B,C$ is non-trivial, and before deriving them, we will justify the appearance of these memory terms through a phenomenological approach in the next subsection.

II.1 Model F dynamics

In a finite temperature phase transition the dynamics of the fluctuations may be understood by treating them as classical stochastic fields. This will be seen by direct calculation later in the paper. However for now, we seek to write a Langevin equation for the (classical) field $\Delta$. In the theory of dynamical critical phenomena, we must account not only for the fluctuating mode, $\Delta$, but also for any conserved quantities which are coupled to $\Delta$. The pairing field $\Delta$ couples directly to the density $\rho$, as given by the commutator $\left[\Delta,\rho\right]=2\Delta$. This can be translated into classical terms by the usual prescription $\left[\cdot,\cdot\right]/i\hbar\rightarrow\{\cdot,\cdot\}$ where $\{\cdot,\cdot\}$ is the classical Poisson bracket.

In the framework of dynamical critical phenomena we specify the static free energy (see Ref. Hohenberg and Halperin, 1977, Eq. (6.3))

[TABLE]

Here $\delta\rho$ is density minus the equilibrium density, $\Phi$ is the scalar potential, $C$ is the compressibility, $r(t)$ is the time dependent detuning from the critical point, $\vec{A}(t)$ is an external vector potential, and $u$, $\gamma_{0}$ are phenomenological parameters. Note that model E Hohenberg and Halperin (1977) has $\gamma_{0}=0$, while model F has $\gamma_{0}\neq 0$.

The field $\Delta$ then has the equation of motion,

[TABLE]

Using Eq. (5), we obtain,

[TABLE]

where ${\xi}_{\Delta}$ is a Gaussian white noise obeying $\langle\xi_{\Delta}(r,t)\xi_{\Delta^{*}}(r^{\prime},t^{\prime})\rangle=2T\Gamma_{\Delta}\delta(r-r^{\prime})\delta(t-t^{\prime})$, and $\Gamma_{\Delta}$ is a dimensionless phenomenological damping.

We will denote the total current by a dissipative component $j_{\rm dis}$, and a superfluid component $j_{\rm sc}$, $j^{i}=j^{i}_{\rm sc}+j^{i}_{\rm dis}$. The equation of motion for the density $\rho$ is given by,

[TABLE]

where the right hand side (rhs) can be interpreted as a definition of $j_{\rm dis}$. Under assumptions that the current may be written as a gradient of a scalar, and that it reaches a well defined steady state in the presence of a dc electric field,

[TABLE]

Above, $\vec{\xi}_{j_{\rm dis}}$ is a Gaussian white noise, whose strength is now controlled by $T\Gamma_{j_{\rm dis}}$, where $\Gamma_{j_{\rm dis}}$ is another phenomenological dissipation. Using the fact that $\partial_{t}\vec{A}=-\vec{E}$, the steady-state dissipative current for a spatially homogeneous system is (setting $e=1$)

[TABLE]

Now using,

[TABLE]

we find that the definition of the superfluid current naturally emerges,

[TABLE]

Substituting Eq. (11) in Eq. (8), we have,

[TABLE]

We now relate these equations to our equations of motion Eq. (4). The expectation value for the supercurrent is

[TABLE]

where $F(q)$ is the equal time $\Delta_{q}\Delta_{q}^{*}$ correlator at momentum $q$. Now we show that the integral or memory term proportional to $\alpha$ in Eq. (4) comes from a term of the form

[TABLE]

and that this term is precisely $\partial_{t}j_{\rm sc}$. The remainder can therefore be identified with the dissipative current giving an equation

[TABLE]

This should be compared to the expectation value for the dissipative current Eq. (10). (Note that the damping coefficient $A$ is not to be confused with the vector potential $A^{i}$). Thus in the dc limit of the kinetic equation, we identify $\Gamma_{j_{\rm dis}}=\tau_{r}A^{-1}$. Note that the model F dynamics are less general in the sense that they assume this dc limit, and therefore that all perturbations are slow compared with $\Gamma_{j_{\rm dis}}$. The kinetic equation does not make this assumption.

Now let us evaluate $\partial_{t}j_{\rm sc}$. We simplify the generalized time-dependent Ginzburg-Landau theory in Eq. (7) by going into Fourier space, and dropping all non-linear terms,

[TABLE]

Above we have encoded the effect of the electric field in two places, one in the minimal coupling $\vec{q}\rightarrow\vec{q}+2e\vec{A}$, and second in the change in the electron density at momentum $q$, where by symmetry,

[TABLE]

where $c$ is a phenomenological constant. Going forward, we absorb this constant into a redefinition of $\gamma_{0}$ appearing in the combination $\gamma_{0}\delta\rho_{q}(t)$.

Defining,

[TABLE]

in the absence of an electric field, the solution is

[TABLE]

Above we have adopted boundary conditions where the superconducting order-parameter and fluctuations are zero at $t=0$. This is because we are interested in quenches where there are initially $t\leq 0$ no attractive interactions, with these being switched on at some arbitrary rate from $t>0$.

The fluctuations involve averaging over noise, giving,

[TABLE]

In the presence of an electric field, and to leading order in it, $\lambda_{q}$ changes to $\lambda_{q}+\delta\lambda_{q}$, with

[TABLE]

where the first term is due to the order-parameter being charged, and the second term is due to the coupling of the order-parameter to the normal electrons whose density is perturbed by the electric field. We have also used that the electric field is $\vec{E}=-\partial_{t}\vec{A}$, so that $\vec{A}(t)=-\int_{s}^{t}dt^{\prime}\vec{E}(t^{\prime})$.

Expanding Eq. (21) in $\delta\lambda_{q}$, which is equivalent to expanding in the electric field,

[TABLE]

Splitting the integral in the exponent, $\int_{t^{\prime}}^{t}dt^{\prime\prime}=\int_{u}^{t}dt^{\prime\prime}+\int_{t^{\prime}}^{u}dt^{\prime\prime}$, and using Eq. (21), we arrive at the following expression for the superconducting fluctuations at momentum $q$, to leading non-zero order in the applied electric field, for an arbitrary time-dependent detuning $r(t)$,

[TABLE]

II.1.1 Aslamazov-Larkin (AL)

We now briefly show how one may recover the familiar AL conductivity in thermal equilibrium. For this it suffices to revert to model E by setting $\gamma_{0}=0$. For a static electric field, and a system in thermal equilibrium, all couplings are time-independent. Thus Eq. (23) becomes,

[TABLE]

Using Eq. (22) for $\gamma_{0}=0$, $\delta\lambda_{q}(u)=\frac{\Gamma_{\Delta}}{M}(2e)(t-u)\vec{q}\cdot\vec{E}$. Taking the long time $t\rightarrow\infty$ limit, and employing the equilibrium expression for $F(q)=T/\epsilon_{q}$, where $\epsilon_{q}=r+q^{2}/2M$, the change in the density of superconducting fluctuations due to the applied electric field is,

[TABLE]

The current is

[TABLE]

leading to the AL conductivity in $d$ dimensions,

[TABLE]

In $d=2$, this reduces to

[TABLE]

A microscopic treatment involving electrons in a disordered potential shows that Larkin and Varlamov (2002) $\Gamma_{\Delta}^{-1}=\pi/8$, giving the well known expression, $\sigma_{\rm AL}^{d=2}=\frac{e^{2}}{16}\frac{T}{r}$. In our disorder-free model, $\Gamma_{\Delta}$ takes a different value.

II.1.2 Kinetic equation for the current

We now return to the derivation of the dynamics of the current in model F. For the current dynamics we need to differentiate Eq. (23) with time,

[TABLE]

The first term on the rhs, being symmetric in momentum space, does not contribute to the current. The second term on the rhs simply provides a time-dependent Drude scattering rate. It is the last term, namely the memory term proportional to the electric field, that is unique to having long lived superconducting fluctuations. We focus only on this term, and using Eq. (22), write it as,

[TABLE]

Above, we have dropped the $E=0$ label in $F$. In the second term, we will find it convenient to replace the time integrals as follows $\int_{0}^{t}du\int_{u}^{t}dt^{\prime}\rightarrow\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}du$.

Comparing equations (4) and (30), we identify,

[TABLE]

The second memory term may be identified as,

[TABLE]

Comparing Eq. (31) and Eq. (32), we obtain,

[TABLE]

The microscopic treatment that follows not only justifies model F, it also provides an explicit calculation for the parameters $\tau_{r}^{-1},A(t),\alpha,\Gamma_{\Delta},\Gamma_{j_{\rm dis}},M,r,\gamma_{0}$.

III $2$PI Equations of Motion

III.1 Properties of $2$PI formalism

We briefly recapitulate the $2$PI formalism here. A detailed explanation can be found for example in Refs Cornwall et al., 1974; Ivanov et al., 1999; Arrizabalaga et al., 2005; Berges, 2004. The $2$PI formalism begins with the Keldysh action $S$ for the Hamiltonian $H$. This is written in terms of Grassmann fields $\psi_{\sigma,n,\pm}(r,t)$,$\bar{\psi}_{\sigma,n,\pm}(r,t)$ where $\sigma$ labels the spin, $n$ labels the orbital quantum number, and $\pm$ label the two branches of the Keldysh contour Kamenev (2011). In terms of these, the Keldysh action is given as

[TABLE]

where $H\left[\bar{\psi}_{\pm}(t),\psi_{\pm}(t)\right]$ indicates the substitution of $\psi_{\pm}(t)$ for $c(t)$ in the Hamiltonian, Eq. (3), in the obvious way. Now we consider a classical field $h$ that couples to a general bilinear of the Grassmann field, and is therefore given by $h_{\sigma\,n\,p}^{\sigma^{\prime}n^{\prime}p^{\prime}}(r,t,r^{\prime},t^{\prime})$, where $p$,$p^{\prime}$ run over $\pm$. In order to simplify notation, we combine the five indices, spin, orbital, Keldysh, position and time into a single vector index, so that the Grassman field becomes a vector $\psi_{i}$ and the source field a matrix $h_{ij}$.

From the action and the source field, we form a generating functional $\mathcal{I}[h]$:

[TABLE]

The first derivative of $\mathcal{I}$ is precisely the Green’s function in the presence of the source field $h$

[TABLE]

In order to work with the physical Green’s function, rather than the unphysical source field $h$, we perform a Lagrange inversion. First we invert Eq. (36) which gives $G$ as a function of $h$, to implicitly define $h$ as a function $G$. Then we construct the Lagrange transformed functional $\mathcal{W}$,

[TABLE]

This has the property that,

[TABLE]

In particular in the physical situation when $h=0$ we have that $\delta\mathcal{W}/\delta G=0$. Thus given the functional $\mathcal{W}$, we have reduced the problem of calculating $G$ to a minimization problem.

The functional $\mathcal{W}$ can be calculated in terms of a perturbation series in $U$, giving

[TABLE]

here the Tr operates on the entire combined index space, and $g^{-1}$ is the bare electron Green’s function. The functional $\Gamma_{2}[G]$ is constructed as follows. First the Feynman rules for the Hamiltonian $H$ are written down. In the present case, this is a single fermionic line and a four-fermion interaction. Then all two-particle irreducible bubble diagrams are drawn. A diagram is two particle irreducible (2PI) if it cannot be disconnected by deleting up to two fermionic lines. The functional $\Gamma_{2}[G]$ is given by the sum of all the diagrams. Lastly the diagrams are interpreted as in usual diagrammatic perturbation theory except that the fermionic line is not the bare Green’s function $g$, but instead the full Green’s function $G$.

The end result is the minimization condition

[TABLE]

This is the Dyson equation for the Green’s function $G$ where the self energy $\Sigma[G]\equiv 2i\delta\Gamma_{2}/\delta G$ is a self-consistent function of $G$. As no approximations have been made Eq. (41) is a restatement of the original many-body problem, and clearly cannot be solved. The advantage of the formalism is that we may replace $\Gamma_{2}$ with an approximate functional $\Gamma_{2}^{\prime}$ and still be guaranteed to preserve conservation laws and causality, which is not the case if one directly approximates $\Sigma$.

We assume that all symmetries are maintained by the solution to Eq. (41). Therefore the Green’s function can be written in terms of the original spin, orbital, Keldysh, space, and time indices as,

[TABLE]

so that only the Keldysh indices $p_{1}$, $p_{2}$ will be explicitly noted.

III.2 RPA approximation to the 2PI generating functional

The approximation we make is the random phase approximation (RPA) in the particle-particle channel, or, equivalently the $N\rightarrow\infty$ where $N$ is the number of orbital degrees of freedom. That is we approximate $\Gamma_{2}$ by $\Gamma^{\prime}_{2}$ where the latter include the set of all closed ladder diagrams,

[TABLE]

shown diagrammatically in Fig. 1, which gives the equation,

[TABLE]

Above we have used that $\Sigma(1,2)=2i\delta\Gamma_{2}/\delta G(2,1)$. The equations may be equivalently derived from the functional for two Green’s functions $G$ and $D$ given by

[TABLE]

This may be interpreted as the functional for electrons interacting with a fluctuating pairing field whose Green’s function is $D$. The $\Gamma^{\prime}_{2}$ is the minimal diagram which has interaction between the fermions and the fluctuating pairing field.

III.3 RPA equations of motion

We now proceed to analyze the RPA equations of motion. As is customary, we perform a unitary rotation of the Keldysh space, defining $G_{K,A,R}$ by Kamenev (2011)

[TABLE]

the last equality holding by causality. We likewise define $\Sigma_{K,A,R},D_{K,A,R},\Pi_{K,A,R}$. In this basis we may write the equations of motion as (setting $e=1$)

[TABLE]

and where $x$ and $y$ stand for combined space, and time indices. The function $G_{A}$ is given by $G_{A}=G_{R}^{\dagger}$, where $O(x,y)^{\dagger}=O(y,x)^{*}$, and likewise for $D_{A},\Sigma_{A},\Pi_{A}$. We note that we also have the relationship $G_{K}^{\dagger}=-G_{K}$

We now convert the definition of $G_{K}$ into a more useful form by evaluating $G_{R}^{-1}\cdot G_{K}-G_{K}\cdot G_{A}^{-1}$. On the one hand, substituting in the definition of $G_{R,A}^{-1}$ this is,

[TABLE]

on the other hand substituting in $G_{K}=G_{R}\cdot\Sigma_{K}\cdot G_{A}$ we obtain

[TABLE]

Therefore we obtain the fundamental equation

[TABLE]

Setting the two times equal and using the form of $g_{R}$ we obtain,

[TABLE]

Above $y$ includes both space and time coordinates, and we have given the collision term $S$ on the rhs for unequal times as we will need it to prove different conservation laws in the next section.

III.4 Conservation laws

We emphasize that the true Green’s function are found by minimizing with respect to the full 2PI potential $\Gamma_{2}$, whereas the functions $G$,$D$ defined above are only approximations, found by minimizing with respect to the approximate potential $\Gamma^{\prime}_{2}$. Therefore it is not apriori clear that there are any conserved quantities that correspond to the conserved quantities of the true Hamiltonian. However, the existence of such quantities is guaranteed by the 2PI formalism. Here we show that the naive expression for the conserved density

[TABLE]

is correct. Setting $r_{1}=r_{2}$ in Eq. (58),

[TABLE]

Inspecting the second term on the left hand side (lhs) we see that it is a total derivative and the equation can be written as

[TABLE]

where the velocity $\vec{v}$ generalized for an arbitrary dispersion $\epsilon$, is defined as

[TABLE]

This has the form of a continuity equation for the conserved quantity $n$ because

[TABLE]

We demonstrate this in Appendix A.1. Thus the particle number $n(r,t)=iG_{K}(rt;rt)$ is a conserved quantity.

Conservation of momentum follows similarly by taking the definition of the momentum

[TABLE]

Evaluating the time-derivative of $p^{\alpha}$, and using Eq. (58), one obtains,

[TABLE]

Above, the electric field is $E^{\alpha}=-\partial_{t}A^{\alpha}$ and the magnetic field is $B^{\alpha}=\frac{1}{2}\epsilon^{\alpha\beta\gamma}\left(\nabla^{\beta}A^{\gamma}-\nabla^{\gamma}A^{\beta}\right)$, and $T^{\alpha\beta}$ is the momentum tensor. Above, we have used the relation

[TABLE]

which is proved in Appendix A.2. We also note that since $S(rt,rt)=0$, $A^{\alpha}(r,t)S(rt,rt)=0$. Thus, the momentum tensor has a contribution from the full electron Green’s function $G$, and also from the full superconducting fluctuations $D$.

Finally conservation of energy is obtained by operating on the kinetic equation $\frac{i}{2}\left(\partial_{t_{1}}-\partial_{t_{2}}\right)$ and setting, $r_{1}=r_{2}$, $t_{1}=t_{2}$. This yields,

[TABLE]

The necessary relation for $S$ needed to derive the above is

[TABLE]

The proof of the above is in Appendix A.3. The above shows that the energy is not conserved due to Joule heating $\vec{E}\cdot\vec{j}$, from an applied electric field, and also when the parameters of the Hamiltonian are explicitly time-dependent $\partial_{t}U(t)\neq 0$.

IV Reduction to a classical kinetic equation

In this section we project the full quantum kinetic equation onto the dynamics of the current. This can be done by noting some important separation of energy-scales associated with critical slowing down of the superconducting fluctuations.

IV.1 Critical slowing down

Although the RPA approximation produces a closed set of equations of motion, Eq. (54), these are a set of non-linear coupled integral-differential equations, and therefore cannot be solved without further reduction. The key to reducing the equations further is the phenomenon of critical slowing down Hohenberg and Halperin (1977). This may be demonstrated by considering the particle-particle polarization in equilibrium. As the equilibrium system is time and space-translation invariant, we consider the Fourier transform

[TABLE]

Assuming for the moment that $G_{R}$ has a Fermi-liquid form

[TABLE]

where $Z$ is the quasiparticle residue, $\tau^{-1}$ is the decay rate and is $\propto{\rm max}(\omega^{2},T^{2})$, while $G_{\rm inc}$ is a smooth incoherent part whose effect is negligible at low energies. We also assume quasi-equilibrium where $G_{K}$ is related to $G_{R}$ by the fluctuation dissipation theorem $G_{K}(k,\omega)=2{\rm Im}\left[G_{R}(k,\omega)\right]\tanh(\omega/2T)$. We expect that $\Pi_{R}(r,t;r,t^{\prime})$ decays exponentially on the scale $T^{-1}$. Thus from the perspective of any dynamics that is slower than $T^{-1}$, the polarization is effectively short ranged in time. Thus, in terms of the Fourier expansion, we only need the leading behavior in $\omega/T$. Note that the short time dynamics showing the onset of quasi-equilibrium of the electron distribution function to a temperature $T$ can also be studied by the numerical time-evolution of the RPA Dyson equations Dasari and Eckstein (2018).

Now the Green’s function $D_{R}$ is precisely the linear response to a superconducting fluctuation. Therefore, if $U<U_{c}$, where $U_{c}$ denotes the critical value of the coupling, the system is in the disordered phase and $D_{R}$ should decay exponentially with time. If $U>U_{c}$ the system is super-critical and superconducting fluctuations should grow exponentially in time. In terms of the Fourier transform $D_{R}(q,\omega)$, this means that there is a pole in the complex $\omega$ plane which crosses the real axis at precisely the phase transition, $U_{c}(T)$. Therefore we may write

[TABLE]

where the function $\lambda(q=0,U_{c}(T))=0$. There is a critical regime in $q$ and $U-U_{c}$ where $\lambda(q,U)\ll T$. We further define a thermal wavelength $q_{T}$ by the expression $\lambda(q_{T},U_{c})=T$, which we estimate as $q_{T}\sim T/v_{F}$, where $v_{F}$ is the Fermi velocity.

On the assumption, to be explicitly shown later, that the transport is dominated by processes that are much slower than $T^{-1}$, we therefore replace

[TABLE]

where we have explicitly included the possibility that $U(t)$ is time dependent. Therefore the equations of motion for $D$ are,

[TABLE]

where $r(t)$ is an effective detuning arising from the the time-dependence of the interaction $U(t)$. Once $D_{R}$ is known then, $D_{K}$ follows from Eq. (54e), and the fact that $\Pi_{K}(t,t^{\prime})=4i\nu\delta(t-t^{\prime})$.

We consider two cases below. One where $\lambda_{q}(t)=\theta(t)\lambda_{q}+\theta(-t)r_{i}$, with $r_{i}/T=\mathcal{O}(1)$. This represents a rapid quench, where the detuning changes from an initial value to a final value at a rate that is faster than or of the order of the temperature. The second situation we consider is an arbitrary trajectory for the detuning $r(t)$.

The solution for bosonic propagators for the former, namely the rapid quench, at times $t,t^{\prime}\geq 0$ then become,

[TABLE]

Note the function $F(q,t)$ changes with rate $\lambda_{q}$ which for small $q$ is $\ll T$. Also note that at $t=t^{\prime}=0$, the density of superconducting fluctuations as measured by $D_{K}(t=t^{\prime}=0)$ or $F(t=0)$ is zero, consistent with an initial condition where the initial detuning is large and positive, and hence far from the critical point.

For an arbitrary trajectory $r(t)$ and hence $\lambda_{q}(t)$, the solutions for the bosonic propagators, generalize as follows,

[TABLE]

We show in Appendix B that the above equations for $D_{R,A,K}$ imply an effective model for a bosonic field obeying Langevin dynamics, where the Langevin noise is delta-correlated with a strength proportional to the temperature $T$.

IV.2 Linear response to electric field

The previous derivation of the kinetic equations Eq. (58) is correct for arbitrary electric fields strengths, and therefore includes non-linear effects. We now expand the Dyson equation to linear order in $E(t)$ in order to evaluate the linear response. Note that since $g_{R}(k,t,t^{\prime})=e^{-i\int_{t^{\prime}}^{t}ds\epsilon(k^{i}+A^{i}(s))}$, it is the combination, $g_{R}(\vec{k}-\vec{A}(t))=e^{-i\int_{t^{\prime}}^{t}ds\epsilon(k^{i}+A^{i}(s)-A^{i}(t))}$ which is gauge-invariant. The leading change in the electron Green’s functions at the gauge invariant momentum is,

[TABLE]

where $n_{k}=\int d^{d}(r_{1}-r_{2})e^{-ik(r_{1}-r_{2})}iG_{K}(r_{1},r_{2})$, and $\delta n$ is assumed to be of order $E$.

We now write the kinetic equation (58) after Fourier transforming with respect to the spatial coordinates $r_{1}-r_{2}$. The change in density at the gauge-invariant momentum is

[TABLE]

where $I_{C}$ is the rhs of the kinetic equation (58) evaluated at the gauge-invariant momentum $\vec{k}-\vec{A}(t)$. Substituting for $I_{C}$, and shifting the internal momentum $q$ by $q-2A(t)$, we obtain

[TABLE]

where the $\delta$ in Eq. (89) indicates the first order variation with respect to the electric field. As expected the conservation laws continue to hold after this approximation. Above, the quantity $\Pi^{\prime}(k,q,t,t^{\prime})$ is defined as,

[TABLE]

Note that all the momenta appear in the gauge invariant combination $\vec{k}-\vec{A}(t)$. Above, the last line highlights the relation between $\Pi^{\prime}$, and the polarization $\Pi$.

The fluctuations are modified by the electric field through their dependence on $\Pi_{R}$:

[TABLE]

In what follows, since we are interested in the response of the current, we multiply the linearized kinetic equation (89) by $v_{k}$ and integrate over $k$.

IV.3 Large fluctuations limit

We make one further simplification, which is based on the fact that in the critical regime $D_{K}\gg D_{R}$ by a factor $T/\lambda_{q}$. Thus we will only keep the terms that are highest order in $D_{K}$.

IV.4 Projections of equation to finite number of modes

The kinetic equation is an integral-differential equation in which the unknowns $\delta n_{k}(t)$ appear linearly. Conceptually therefore it may be solved by standard methods. However, even though it is linear, it is non-local in time and not time translation invariant. This makes direct analytical and numerical solution difficult.

We therefore make an additional simplification that instead of considering the full space of solutions, we instead project entirely onto the current mode. This means that we would like to fix,

[TABLE]

where $\rho=-\sum_{k}n^{\rm eq}_{k}$ is the density of fermions, $\bar{m}$ is the effective mass and $J^{i}(t)$ is some function to be determined. The $?$ in Eq. (95) is to question the correctness of this equation. This is because conservation of momentum $P^{i}=\sum_{k}k^{i}n_{k}$ causes Eq. (95) to fail qualitatively. As $\partial_{t}P^{i}=\rho E^{i}$, an electric field pulse, say for simplicity taken to be a delta-function in time, will generate a net momentum. Moreover, by conservation of momentum, the initial perturbation

[TABLE]

with $V^{i}$ arbitrary, will never decay. Therefore instead of Eq. (95) we decompose the occupation number as

[TABLE]

The parameter $\gamma$ gives the amount of current that is carried by the momentum mode, which does not relax. In the limit $\gamma\rightarrow 0$ there is no overlap between the modes, the momentum mode carries no current, and is thus neglected. On the other hand, when $\gamma\rightarrow 1$, as in a Galilean invariant system, the current and momentum are proportional and there is no relaxing current.

Multiplying the kinetic equation by $v_{k}^{i}-\gamma k^{i}/\bar{m}$ and $k^{i}$, and summing over $k$ we obtain,

[TABLE]

Above, in the first equation we have used that $\sum_{k}(v_{k}-\gamma k/\bar{m})\delta n_{k}=(1-\gamma)J^{i}_{r}$ as follows from Eq. (100). The second equation above just follows from conservation of momentum. We also do not write the explicit dependence of the momentum labels on the vector potential as it is understood that it is always the gauge-invariant combination that appears.

As the momentum mode has trivial dependence due to momentum conservation, we will simply drop it and set $\gamma=0$ in Eq. (101). In this limit,

[TABLE]

These are now a closed set of equations in terms of the single unknown function $J_{r}(t)$ which we will denote simply as $J(t)$. The remaining step is to evaluate the various terms. We note that $D_{K}(t,t^{\prime})$ is generally larger than $D_{R}$ by the factor $T/\lambda_{q}$. Thus in what follows we will only retain the first and third terms on the right hand side of Eq. (101).

IV.5 Time dependent kinetic coefficients

We begin by evaluating the first term, proportional to $\delta\Pi^{i}_{R}(q,t,t^{\prime})$ in Eq. (101). There are two contributions, one from varying $n_{k}$, and we denote it by $K_{J}$. The second is from varying $g_{R}$, and we denote it by $K_{E}$

[TABLE]

In Appendix C.1 we show that the term coming from varying $n_{k}$ is

[TABLE]

To estimate the magnitude of the above term, we neglect factors of order one, write $\rho\sim\nu k^{2}_{F}/m$, and obtain that this term is $\sim(q/k_{F})^{2}/\nu$.

For the term coming from varying $g_{R}$, we show in Appendix C.2 that,

[TABLE]

The resulting term has the form $\propto(q/k_{F})^{2}F_{q}(t)\rho E(t)/\nu T\bar{m}$, and is thus a parametrically small correction to the drift term $\rho E/\bar{m}$. Therefore we neglect it for the remainder.

We now turn to the third term $\sim\delta D_{K}\Pi^{\prime}_{A}$ in Eq. (101). We begin by considering the change in $\delta D_{K}$. This in turn depends on $\delta D_{R}$, which from Eq. (77) is given by

[TABLE]

There are two contributions to $\delta\lambda_{q}$. One is from varying $\delta n$ via

[TABLE]

Writing the above as a small $q$ expansion, we define $M$ such that,

[TABLE]

where we find that, on expanding,

[TABLE]

where $d$ is the spatial dimension. The above equation also shows, $M\sim T/v_{F}^{2}$

Thus, the change in $\lambda_{q}$ coming from the change in the electron distribution $\delta n_{k}$, may be evaluated as follows, (below we also use $\delta n_{k}=-\delta n_{-k}$),

[TABLE]

Using the relation between $\delta n_{k}$ and the current in Eq. (103), and the expression for $M$ in Eq. (110), we find,

[TABLE]

The second reason for the change in $\lambda_{q}$ is due to the direct coupling to the electric field, where defining $\delta A^{i}(s)=A^{i}(s)-A^{i}(t)$,

[TABLE]

Thus the total change in $\lambda_{q}$ is

[TABLE]

$\delta\lambda_{q}$ changes both $D_{R}$ and $D_{K}$, where the change in the former is

[TABLE]

Note that $\delta D_{K}=(\delta D_{R})\Pi_{K}D_{A}+D_{R}(\delta\Pi_{K})D_{A}+D_{R}\Pi_{K}(\delta D_{A})$. Since the variation $\delta\Pi_{K}$ produces a term that is smaller by a factor of $F_{q}$, we neglect it. Using the zeroth order calculation that $i\Pi_{K}\propto\nu\delta(t-t^{\prime})$, we have,

[TABLE]

Substituting in for the result for $\delta D_{R}$, using, $\delta D_{R}(t,t^{\prime})=\delta D_{A}(t^{\prime},t)^{*}$, we obtain (see Appendix C.3 for details),

[TABLE]

The coefficient $\alpha$ is given by,

[TABLE]

The above shows that $\alpha q^{2}\sim(q/k_{F})^{2}/\nu$, which is comparable to the local terms. Moreover the sign of $\alpha$ for generic band structures is such as to oppose the local in time term coming from $iD_{K}\delta\Pi_{A}$.

The function $\kappa(\epsilon)$ may be Taylor expanded to give $\kappa(\epsilon)\approx\frac{\nu}{\bar{m}}(1+b\epsilon/E_{F}+\cdots)$, where the parameter $b$ is a material dependent parameter of order $1$. The sign of $b$ is not fixed in general. However for a ’simple’ band structure, the parameter $b$ is negative. Thus,

[TABLE]

For generic quench profiles, the only change that is needed, is the replacement

[TABLE]

In what follows we set $\rho/\bar{m}=1$. In these units, the ratio $J/E$ gives a Drude scattering time, and also equals the conductivity. The kinetic equation for the current for general quench trajectories becomes Eq. (4), and for convenience we rewrite it below,

[TABLE]

We define the dimensionless quantity,

[TABLE]

The memory terms are,

[TABLE]

We now simplify the expressions for $A,B,C$ using the derived expressions for the time-dependent superconducting fluctuations.

[TABLE]

Similarly the memory term is,

[TABLE]

For the case of the sudden quench, since $r(t)=\theta(t)r+r_{i}$, the above expressions for $A,B,C$ simplify considerably. For $r_{i}/T=\mathcal{O}(1)$ so that at $t=0,F(q,t=0)=0$, and for a quench to the critical point $r=0$, equations (125), (127) and (124) give,

[TABLE]

Above we have regularized the integrals such that a short time cutoff of $T^{-1}$ has been introduced. It is also helpful to study the system at non-zero detuning. For a sudden quench to a distance $r$ from the critical point, equations  (125),  (127) and (124) give,

[TABLE]

IV.6 Charge Diffusion

We now adapt the kinetic equation to the case where the current is being driven by a density gradient. From Eq. (58), Fourier transforming with respect to the difference in position coordinates $r_{1}-r_{2}$, and performing a gradient expansion with respect to the center of mass coordinate $r=\frac{r_{1}+r_{2}}{2}$,

[TABLE]

where $S(k)=\int d^{d}(r_{1}-r_{2})e^{-ik(r_{1}-r_{2})}S(r_{1}t,r_{2}t)$. A spatial density gradient leads to a current $J^{i}=\sum_{k}v_{k}^{i}n_{k}$. Multiplying Eq. (134) with $v_{k}$, taking a sum on $k$, and using the fact that $v_{k}S(k)$ in the presence of the current can be written as a contribution coming from small changes to the superconducting propagator, and small changes to the polarization bubble as summarized on the rhs of Eq. (89). The rhs can be simplified as before, leading to the rhs of Eq. (121). This leads to,

[TABLE]

The coefficients $A,B$ are given in equations  (125), (127) for a general quench profile, and in equations  (128), (129) for a rapid quench to the critical point. Although there is no external electric field, the density gradient induces a current, and the electric field on the rhs is a linear response to this current.

V Conductivity and Diffusion coefficient in thermal equilibrium

In thermal equilibrium, one may simply take the long time limit of the expressions in equations (128), (129), and  (130). Then, at zero detuning $r=0$,

[TABLE]

It is also useful to analyze these coefficients when the system has equilibrated at a non-zero detuning $r$ away from the critical point. In this case, for $d=2$ we have,

[TABLE]

Above due to the separation of time-scales discussed in the previous section, $|r|/T\ll 1$. We note that in Fourier space,

[TABLE]

Above $\gamma$ is Euler Gamma.

We now discuss the linear response conductivity in thermal equilibrium. In this case, all the coefficients $A,B,C$ are time-translation invariant. Fourier transforming Eq. (121), we obtain,

[TABLE]

Writing,

[TABLE]

where the first is a dissipative current arising due to Drude scattering, while the second is a current arising due to the superconducting fluctuations. Then, taking the dc limit of Eq. (140),

[TABLE]

while the current from the superconducting fluctuations (neglecting $B$ which gives small correction to the Drude scattering rate) is,

[TABLE]

The above implies that the fluctuation conductivity $\sigma_{\rm fl}$ gives a correction to the dissipative conductivity $\tau_{r}$ by an amount $=\tau_{r}^{-1}\alpha\tilde{C}_{\rm eq}(r,\omega=0)$. Noting that $\tau_{r}\sim T^{-1}$, and using Eq. (139), this implies a fluctuation conductivity correction that goes as $\propto\alpha T/r$. While this is qualitatively the same as the fluctuation AL conductivity Azlamazov and Larkin (1968a, b), we note that for the ultra clean case considered here, material dependent parameters such as $\alpha$ which determine the strength of Galilean invariance breaking lattice effects, are unavoidable, and give a non-universal material dependent prefactor.

The MT correction is also discussed in the context of a disordered system Thompson (1970); Maki (1968a, b); Larkin and Varlamov (2002), and is well defined only as long as the mean free path is shorter than the inelastic scattering time Randeria and Varlamov (1994); Stepanov and Skvortsov (2018). Thus the MT conductivity does not emerge naturally in the clean limit we consider here.

Now we discuss the diffusion constant $D$ defined as $J^{i}=-D\nabla^{i}n$. Due to time-translation invariance in equilibrium, we write Eq. (135) in Fourier space. We also note that the electric field generated as a response to the current, is given by $J=\tau_{r}E$, up to fluctuation corrections. This gives the diffusion constant,

[TABLE]

Close to the critical point, $C$ dominates over $B$. Moreover, Taylor expanding in $\alpha$, the dc limit of the diffusion constant is,

[TABLE]

Thus the fluctuation correction for the diffusion constant has qualitatively the same form as that for the conductivity in being $\propto\alpha T/r$, where we have used that $\tau_{r}\sim T^{-1}$.

VI Solving the kinetic equation for the conductivity

The solution of the kinetic equation for the current was presented in much detail in Ref. Lemonik and Mitra, 2018b. For completeness we summarize some of the findings. Two kinds of quench trajectories were studied. One was a rapid quench from deep in the disordered phase to a distance $r\geq 0$ away from the critical point. The second was a smooth quench trajectory which started from the disordered phase, approached the critical point $r=0$, and returned back to the disordered phase. Regardless of the details of the trajectory, the current showed slow dynamics because of slowly relaxing superconducting fluctuations. In the kinetic equation Eq. (121), this physics is encoded in the memory terms ($B,C$).

In addition, the dynamics of the current for a critical quench to $r=0$, was found to show universal behavior, with a power-law aging in the conductivity $\sigma(t,t^{\prime})$, a result also found for a disordered system Lemonik and Mitra (2018a). Moreover, the dynamics at non-zero detuning $r$ was shown to obey scaling collapse. However it should be noted that, for the clean system studied here, the exponents entering the scaling behavior were non-universal in that they depend on Lemonik and Mitra (2018b) $T\tau_{r}$.

For a smooth quench it was shown that transiently enhanced superconducting fluctuations create transient a low resistance current carrying channel, whose signature is a suppression of the Drude scattering rate at low frequencies. Note that the Drude scattering rate in the absence of time-translation invariance is defined as,

[TABLE]

with

[TABLE]

In this section we revisit the smooth quench, defined by the trajectory

[TABLE]

Above, $r/T$ starts out being $1$, smoothly approaches $-\epsilon$ at a time $Tt_{*}=30$, and then smoothly returns back to its initial value of $r/T=1$. We consider two cases, one where $\epsilon=0$ (see Fig. 2), and one where the quench is super-critical in that $\epsilon=0.05$. For the latter case, for a finite time, the parameters of the Hamiltonian correspond to that of an ordered phase since $r<0$ (see Fig. 3). Despite this, for such a transient regime, the order-parameter is still zero, although superconducting fluctuations are significantly more enhanced than if $r\geq 0$. In particular, critical slowing down prevents true long range order to develop over time scales over which the microscopic parameters are varying.

Note that in the numerical simulation, we apply an electric field pulse which is a delta-function in time. In particular, for an electric field pulse of unit strength centered at $t^{\prime}$, the conductivity equals the current, $\sigma(t,t^{\prime})=J(t)$. Fourier transforming with respect to $t-t^{\prime}$, we can obtain the conductivity for arbitrary $\omega$. However, in an actual experiment, the finite temporal width of the electric field pulse places a limit on the lowest frequency accessible. Nevertheless, the physics of a suppressed Drude scattering at low frequencies is visible over a sufficiently broad range of frequencies, as shown in figures 2 and 3, for this to be a feasible experimental observation.

Top panels of Fig. 2 and Fig. 3 show how the superconducting fluctuations, at several different wavelengths, evolve for a critical and a super-critical quench respectively. The lower panels of the same figures show how the corresponding Drude scattering rate, for different frequencies, evolve in time. When the density of superconducting fluctuations peak, the Drude scattering rate at low frequencies $\omega<T$ dips, with the effect being more enhanced for the super-critical quench. Moreover the Drude scattering rate is strongly dispersive in that different frequency components of the Drude scattering rate peak at different times.

Recall that when the system returns to the normal phase, since we are in the clean limit, the steady-state Drude scattering rate is zero. A slow power-law approach to thermal equilibrium is also visible in the long time tails of the Drude scattering rate that are found to persist long after the detuning has returned back to its initial value in the normal phase. The slow dynamics highlighted above provide signatures in time-resolved optical conductivity where, although the system lives for too short a time for true long range order to develop, the transient dynamics can still show clear signatures of superconducting fluctuations.

VII Spectral properties

We now present results for the time-evolution of the electron spectral properties. Results for the electron life-time obtained from the imaginary part of the electron self-energy, appear Ref. Lemonik and Mitra, 2017, where it was shown that the behavior is non-Fermi liquid like due to enhanced Andreev scattering of the electrons at the Fermi energy. In this section we discuss how this scattering affects the local density of states.

For simplicity we consider a rapid quench from deep in the normal phase where $F(q)=0$ initially, to a final arbitrarily small detuning $r\geq 0$ away from the critical point. For this quench, the superconducting fluctuations obey the equation,

[TABLE]

where the solution of the above is identical to Eq. (82) with $\lambda_{q}/T=r+v^{2}q^{2}/T^{2}$. We have assumed that $F$ depends only isotropically on $q$, and follow a convention where $r$ is dimensionless. We define the dimensionless variable $x=v^{2}q^{2}/T^{2}$. The electron self energy using Eq. (54f) and in the limit of $D_{K}\gg D_{R}$, and in two spatial dimensions is,

[TABLE]

Going into Wigner coordinates which involves Fourier transforming with respect to relative coordinates $t_{1}-t_{2}$, while keeping the mean time $t=(t_{1}+t_{2})/2$ dependence in $D_{K}$ Lemonik and Mitra (2017), and performing an angular integral, we obtain,

[TABLE]

We start in the perturbative limit by setting $\Sigma_{A}=i0^{+}$, $\Sigma_{R}=-i0^{+}$. First let us study the equilibrium case. Here $F=1/(x+r)$. Thus we have,

[TABLE]

The density of states is,

[TABLE]

To evaluate this integral, we need to locate the pole in $\varepsilon$ complex plane: $\omega-\varepsilon^{*}(\omega)-\Sigma^{R}(\varepsilon^{*}(\omega),\omega)=0$. Then the value of the integral follows from the residue

[TABLE]

Thus, the perturbative change in the density of states is obtained from Taylor expanding,

[TABLE]

Using Eq. (157), one obtains the following asymptotic behavior,

[TABLE]

Note that this shift is perturbative for all $\omega$, that is $\delta\nu/\nu\ll 1$, when $\text{Gi}/4\pi r\ll 1$.

Now we turn to the critical quench, where $r$ is $\mathcal{O}(1)$ for $t<0$ and $r=0$ for $t>0$. In this case, the solution of Eq. (149) gives,

[TABLE]

The above implies,

[TABLE]

In Appendix D we show that $h$ has the following asymptotic form,

[TABLE]

Proceeding to the change in the density of states we have

[TABLE]

Note that the asymptotic behaviors of the critical and equilibrium cases agree if we set $r=\frac{1}{2Tt}$. This follows from the fact that the fluctuations $F(x,t)$ behave in a comparable manner in these two cases.

VII.1 Self-Consistent Regime (Results for $\omega=0$)

We now turn to the full solution that does not require $\delta\nu\ll\nu$. We will first work out the equilibrium case. The critical quench case will then follow from the substitution of $r=\frac{1}{2Tt}$.

By adding $\pm\omega+\varepsilon$, the self-consistent equation for the self-energy in equations  (151), (152) may be written as,

[TABLE]

where

[TABLE]

The above self-consistent equations can be solved in certain limits (see Appendix E) giving the following results for the density of states at zero frequency,

[TABLE]

The critical quench case then follows by the substitution $r=\frac{1}{2Tt}$,

[TABLE]

Thus we find that there is no true steady state of the zero frequency density of states due to the solution continuing to evolve logarithmically in time. This time-evolution, although slow, happens because there is no non-zero temperature critical point in $d=2$. To obtain the correct thermalized steady-state, vortices need to be accounted for in our treatment. However, for the transient regime relevant for experiments, where no superconducting gap has developed yet, and the dynamics is governed by weakly interacting superconducting fluctuations, our results are still valid.

VIII Conclusions

We have studied quench dynamics of an interacting electron system along general quench trajectories using a 2PI approach. To leading non-trivial order in $1/N$, our treatment reduced to RPA in the particle-particle channel. While the quantum kinetic equations are exact, leading to proper definitions of the conserved densities and currents, progress was made in solving them by exploiting a separation of time-scales. This could be seen by noting that thermalization of the electron distribution function sets in on a time-scale of $T^{-1}$, while the collective modes relax at a rate given by the detuning to the superconducting critical point, which can be arbitrarily slow. This is the phenomena of critical slowing down, and we derived an effective classical equation for the decay of the current in this regime. The dynamics of the current was also justified using a phenomenological mapping to model F in the Halperin-Hohenberg classification scheme. Going forward, we envision that such mappings may be helpful in other contexts when understanding time-resolved experiments that aim to probe collective modes.

Results were also presented for the local density of states which was found to be enhanced at low frequencies due to Andreev scattering processes. For a critical quench, the self-consistent equations showed an equivalence between the time after the quench and the inverse detuning.

Future directions involve studying the coarsening regime, in particular how true long range order develops. The reverse quench where the initial state is ordered, also needs to be explored, going beyond mean-field. Experiments involving ultra-fast manipulation of spin and charge order also require generalizing our study to quench dynamics in the particle-hole channel.

Acknowledgments: This work was supported by the US National Science Foundation Grant NSF-DMR 1607059.

Appendix A Proving conservation laws

A.1 Proof of $S(rt,rt)=0$

Using Eq. (54),

[TABLE]

The above terms can be rearranged to give,

[TABLE]

Using Eq. (54e), the second and fourth terms above are rewritten to give,

[TABLE]

Using Eq. (54d) and writing, $\Pi_{A,R}D_{A,R}=U^{-1}D_{A,R}-1$, or equivalently $D_{A,R}\Pi_{A,R}=D_{A,R}U^{-1}-1$,

[TABLE]

A.2 Proof of Eq. (67)

It is convenient to write $S$ for unequal positions and times, and express the self-energy $\Sigma$ in terms of the $G,D$ propagators. Doing this we obtain,

[TABLE]

In order to prove momentum conservation we need to study the action of the spatial gradients on $S$,

[TABLE]

where we have used that \left(\partial_{r_{1}}-\partial_{r_{2}}\right)\left[G_{a}(1^{\prime},1)G_{b}(1^{\prime},2)+G_{b}(1^{\prime},1)G_{a}(1^{\prime},2)\right]\biggl{|}_{1=2}=0. Now using that

[TABLE]

and performing a similar analysis for second last term in Eq. (182), we obtain,

[TABLE]

A.3 Proof of Eq. (72)

The manipulations involved are,

[TABLE]

Appendix B Equivalence to Langevin dynamics

Using the notation of Ref. Kamenev, 2011, the bosonic quantum ($\Delta_{q}$) and classical ($\Delta_{c}$) fields defined on the Keldysh contour, have the following correlators,

[TABLE]

above $1,2$ denotes both spatial and temporal indices. If, $D_{K}=D_{R}\cdot\Pi_{K}\cdot D_{A}$, then the above correlators can equivalently be written as a path integral over the bosonic fields as follows,

[TABLE]

According to Eq. (77), and the line below,

[TABLE]

One may introduce Kamenev (2011) an auxiliary field $\xi$ that decouples the $|\Delta_{q}|^{2}$ term in the action, (for notational convenience we only highlight the time coordinate)

[TABLE]

Writing the action in terms of this auxiliary field,

[TABLE]

Requiring $\delta S/\delta\Delta_{q}^{*}=0$ gives,

[TABLE]

where $\langle\xi^{*}(t)\xi(t^{\prime})\rangle=4\nu\delta(t-t^{\prime})$. Giving the explicit expression for $D_{R}$, and restoring the momentum label,

[TABLE]

It is convenient to rescale $\Delta_{c}\rightarrow\sqrt{\frac{\nu}{T}}\Delta_{c}$, then, the Langevin equation is,

[TABLE]

implying that $\Delta_{c}$ obeys Langevin dynamics where the noise is delta-correlated with a strength equal to the temperature $T$.

Appendix C Derivation of kinetic coefficients

In this section we derive the kinetic coefficients assuming a rapid quench profile $\lambda_{q}(t)=\theta(t)\lambda_{q}+\theta(-t)r_{i}$, where $r_{i}/T=O(1)$. This condition simply ensures that the density of superconducting fluctuations is zero at $t=0$. The generalization to general quench profiles is formally straightforward, and summarized in the main text.

C.1 Derivation of Eq. (105)

The term obtained from varying $n_{k}$ is,

[TABLE]

Substituting Eq. (81) for $D_{K}$, and Eq. (103) for $\delta n_{k}$, we see that the above can be rewritten as

[TABLE]

Following the discussion in the main text, the term $Q^{ij}_{J}$ decays exponentially with rate $T$, whereas $J$ and $F$ change at a slower rate. Therefore $Q_{J}$ appears like a delta function when integrated against $F(t)J(t)$, and we may write

[TABLE]

In the second line above we assume $t^{-1}\gg T$ so we are not too close to the quench. In the third line we take $q\ll T/v_{F}$, and we also assume $\lambda_{q}/T\ll 1$.

C.2 Derivation of Eq. (106)

Using Eq. (81) for $D_{K}$, and accounting for the change in $g_{R}$ due to the electric field in Eq. (86), we find,

[TABLE]

Assume that the frequency of the electric field $\omega\ll T$ - in this case we may approximate the integral

[TABLE]

Following the discussion in the main text, the term $Q^{ij}_{E}$ decays exponentially with rate $T$, whereas $J$ and $F$ change at a slower rate. Therefore $Q_{E}$ appears like a delta function when integrated against $F(t)J(t)$, and we may write, as we did for $Q^{ij}_{J}$,

[TABLE]

C.3 Derivation of Eq (117)

Substituting Eq. (114) into Eq. (115), and assuming, $t>t^{\prime}$, we obtain,

[TABLE]

Using the definition of $D_{K}=D_{R}\cdot\Pi_{K}\cdot D_{A}$ with $\Pi_{K}(t,t^{\prime})=4\nu\delta(t-t^{\prime})$, the above becomes,

[TABLE]

Then we use Eq. (81) to write,

[TABLE]

where in the last line we used the fact that $t-t^{\prime}$ is of order $T^{-1}$ as $iD_{K}$ appears along with $\Pi_{A}$. We discuss this last term further below and show it to be parametrically smaller than the local terms calculated earlier.

We insert the second part into the collision integral giving,

[TABLE]

Thus this term is parametrically smaller than the local term already calculated.

Now we consider the full term

[TABLE]

Above, in the second last line, we used the expression for $\delta\lambda_{q}$ in Eq. (114).

Appendix D Derivation of Eq. (167)

Let us now discuss the limiting behaviors of $h(z_{0})$. If $z_{0}\ll 1$, the integral is restricted to the region $\bar{x}\ll 1$ and we may approximate $\left(1-e^{-\bar{x}}\right)/\bar{x}\approx 1$. Therefore in this limit we have $\text{Re}\left[h\right]\approx 2z_{+}^{0}$. If $|z_{+}^{0}|\gg 1$ then we split the integral at $\bar{x}=1$. The lower part

[TABLE]

Whereas in the part from 1 to $\infty$ the exponential may be neglected leaving the leading part

[TABLE]

Thus we have Eq. (167).

Appendix E Derivation of Eq. (175)

We record here some facts about $f(z)$ as a function in the complex plane. It may be written as

[TABLE]

It is analytic in the upper and lower half planes separately with a branch cut along the real axis. The points $z=\pm i$ which appear singular are in fact smooth. One may also convince one self that $|f(z)|\leq\pi$ for all $z$, reaching this limit only when $z=0\pm i\delta$. It is also an injective function of $z$. In addition, we have that as $z\rightarrow\infty$,

[TABLE]

and that $f$ is pure imaginary on the imaginary axis. Causality requires that $z_{\pm}$ have no zeros or branch cuts in the upper half complex $\omega$ plane, but there will be analytic structure as a function of $\varepsilon$.

Let us first test the validity of the perturbative approximation. To zeroth order we set $z_{\pm}=z_{\pm}^{0}$. To first order we substitute this into the rhs of equation (170), (171), and obtain,

[TABLE]

and to second order

[TABLE]

Comparing the first and second correction we see that the latter is smaller when

[TABLE]

Considering the various limits, we see that the rhs is always at least $O(1)$ so the perturbative treatment is always valid if $\gamma\ll 1$. Making the very coarse approximation $|f|\sim\min(1,|z^{-1}|)$, we get roughly

[TABLE]

This may still be satisfied even if $\gamma\geq 1$. If $\omega/T\gg 1/\text{Gi}$, the perturbative limit will hold for all $\varepsilon$. With the weaker condition $\omega^{2}/T^{2}\gg\text{Gi}$, the perturbative limit holds except in the narrow region $|\varepsilon\pm\omega|/T\ll\sqrt{\epsilon}$, which gives only a small contribution to $\delta\nu(\omega)$.

Now let us consider the behavior in the simplest case when $z_{\pm}^{0}=0$. In this case we have a self-consistent solution when $z_{\pm}=\mp iy$, $y$ real, $y>0$, since this is equivalent to finding a solution to

[TABLE]

which is possible since $f$ is imaginary on the imaginary axis. In the limiting case of $\gamma\gg 1$, we get $y^{2}\sim\gamma\log y^{2}$ which gives to logarithmic accuracy,

[TABLE]

In the limit $\gamma\ll 1$, using the expansion of $f(z)$ for small argument, we obtain $z_{\pm}\approx\mp i\gamma\pi+\mathcal{O}(\gamma^{2})$.

Now let us advance to the case $z_{\pm}^{0}=c\in\mathbb{R}$. Recall $\omega=0$. This has a self-consistent solution when $z_{+}=z=z_{-}^{*}$,

[TABLE]

since $f(z)^{*}=f(z^{*})$. The pertubative condition in this case is

[TABLE]

Thus the perturbative regime begins when $c\gg y(\gamma)$. In this limit, we obtain the usual perturbative result by setting $z\approx c$.

In the opposite limit, let us assume $|z|\gg 1,|c|$. On these assumptions, writing $z=-i|z|e^{-i\phi}$, and comparing real and imaginary parts of

[TABLE]

we obtain,

[TABLE]

On the assumption $|z|\gg|c|$ the first line above gives $|z|=y(\gamma)$. The second line can be rearranged to give,

[TABLE]

Thus $\phi$ goes from $\pi/2$ to $-\pi/2$ as $c$ goes from $-\infty$ to $\infty$. We may split this into two limits. When $c\ll\gamma/y(\gamma)=y(\gamma)/\log y(\gamma)$, we have that $\phi\ll 1$ and so the $\cos$ may be set to $1$ giving $\phi=-cy(\gamma)/2\gamma$, and thus

[TABLE]

In the other limit, $cy(\gamma)/\gamma\gg 1$ is large so $\phi$ is close to $\pm\pi/2$. Let us take $c>0$, then

[TABLE]

so that $z$ is given by,

[TABLE]

This is almost identical to the perturbative calculation with a slight correction to the real part.

We can collect all of these limits,

[TABLE]

where in the last case, we have used the perturbative expansion in Eq. (211).

We are now in a position to estimate the density of states for $\omega=0$, $\gamma\gg 1$, and $c=\varepsilon/T\sqrt{r}$. Starting from

[TABLE]

Combining the positive and negative integrals over $\varepsilon$ or $c$,

[TABLE]

Next use the fact that at $\omega=0$, $z(c)=-z(-c)^{*}$, which follows from the above calculation. This can also be seen from noting that $f(-z)=-f(z)$. Then, using $z_{+}(c)=c+i\delta+\gamma f(z_{-}(c))=c+i\delta+\gamma f(z_{+}^{*}(c))$. Taking $c\rightarrow-c$ and conjugating, this becomes, $z_{+}^{*}(-c)=-c-i\delta+\gamma f(z_{+}(-c))=-c-i\delta-\gamma f(-z_{+}(-c))$. Thus $z_{+}(-c)^{*}=-z_{+}(c)$ is a solution of the above equation as this simultaneously requires $f(-z_{+}(-c))=f(z_{+}^{*}(c))$. Using this gives,

[TABLE]

Now we split this integral into three regions, according to the asymptotics above,

[TABLE]

where,

[TABLE]

[TABLE]

Here we are treating $\log\log\gamma$ as an $\mathcal{O}(1)$ constant, which has been the order of our precision throughout.

[TABLE]

Thus we see that all terms contribute $\propto 1/\log{\gamma}$ and therefore we have that

[TABLE]