Time-resolved spectral density of interacting fermions following a quench to a superconducting critical point

TL;DR

This paper investigates the time evolution of fermions after a quench to a superconducting critical point, revealing universal dynamics and a singular spectral density correction linked to superfluid fluctuations.

Contribution

It demonstrates that fermion fluctuations follow universal model A dynamics near the critical point using a large N approximation.

Findings

Universal scaling form for fermion lifetime correction.

Singular behavior of spectral density due to incoherent Andreev reflections.

Analysis of dynamics near superconducting critical point.

Abstract

Results are presented for the time evolution of fermions initially in a non-zero temperature normal phase, following the switch on of an attractive interaction. The dynamics are studied in the disordered phase close to the critical point, where the superfluid fluctuations are large. The analysis is conducted within a two-particle irreducible, large $N$ approximation. The system is considered from the perspective of critical quenches where it is shown that the fluctuations follow universal model A dynamics. A signature of this universality is found in a singular correction to the fermion lifetime, given by a scaling form $t^{(3-d)/2}S_d(\varepsilon^2 t)$, where $d$ is the spatial dimension, $t$ is the time since the quench, and $\varepsilon$ is the fermion energy. The singular behavior of the spectral density is interpreted as arising due to incoherent Andreev reflections off superfluid fluctuations.

Figures

Equations

$$H_i=\sum _{k,\sigma =\uparrow ,\downarrow ,\tau =1\dots N}{ϵ}_k{c}_{k\sigma \tau }^{†}{c}_{k\sigma \tau }.$$

$$H_f=H_i+\frac{u}{N}\sum _q{\Delta }_q^{†}{\Delta }_q,$$

$${\Delta }_q=\sum _{k\tau }{c}_{k,\uparrow ,\tau }{c}_{-k+q,\downarrow ,\tau };{\Delta }_q^{†}=\sum _{k,\tau }{c}_{-k+q,\downarrow ,\tau }^{†}{c}_{k\uparrow \tau }^{†}.$$

$$G_R(1,2)$$

$$G_K(1,2)$$

$${G_R}^{-1}={g_R}^{-1}-{\Sigma }_R[G],$$

$${\Sigma }_R[G]\equiv \delta {\Gamma }^{\prime }/\delta G_A,$$

$${\Sigma }_R(1,2)=\frac{i}{N}\left[D_K(1,2)G_A(2,1)+D_R(1,2)G_K(2,1)\right]$$

$$D_R^{-1}$$

$$D_K$$

$$i{\Pi }^K(q,t,t^{\prime })$$

$$i{\Pi }^R(q,t,t^{\prime })$$

$$D_R(1,2)$$

$$D_K(1,2)$$

$${\Pi }_R(v_F|q|\ll T,\omega \ll T)$$

$$i{\Pi }_K(v_F|q|\ll T,\omega \ll T)$$

$$\Rightarrow i{\Pi }_K(q=0,t,t^{\prime })$$

$$\displaystyle\biggl{[}\partial_{t}+\gamma_{q}\biggr{]}D_{R}(q,t,t^{\prime})=-Z\delta(t-t^{\prime}),$$

$${\gamma }_q=T\left(l^2{q}^2+r\right).$$

$$D_R(q,t,t^{\prime })$$

$$\displaystyle iD_{K}(q,t,t^{\prime})=Z\frac{T}{2\gamma_{q}}\biggl{[}e^{-\gamma_{q}|t-t^{\prime}|}-e^{-\gamma_{q}(t+t^{\prime})}\biggr{]}.$$

$$F_K^0(x)=\frac{1-e^{-2x}}{2x},$$

$$\displaystyle iD_{K}(q,t,t^{\prime})=-T\biggl{[}D_{R}(q,t,t^{\prime})t^{\prime}F_{K}^{0}(\gamma_{q}t^{\prime})$$

$$\displaystyle+tF_{K}^{0}(\gamma_{q}t)D_{A}(q,t,t^{\prime})\biggr{]}.$$

$$2x{F}_K^0\left(x\right)-1\ne 0.$$

$$i{D}_K\left(q,t,t\right)=ZTt{F}_K^0\left(T{l}^2{q}^{2}t\right);v_{F}q\ll T.$$

$${\Sigma }_R^{\mathrm{WT}}\left(k,\omega ;t\right)\equiv \int d\tau e^{i\tau \omega }{\Sigma }_R\left(k,t+\frac{\tau }{2},t-\frac{\tau }{2}\right),$$

$$\frac{1}{\tau\!\left(k,t\right)}\equiv\textrm{Im}\biggl{[}\Sigma_{R}^{\rm WT}\left(k,\omega=\varepsilon_{k};t\right)\biggr{]},$$

$${\Sigma }_R\left(k;t_1,t_2\right)$$

$$\displaystyle+G_{A}\!\left(-k+q;t_{2},t_{1}\right)D_{K}\!\left(q;t_{1},t_{2}\right)\biggr{]}.$$

$${\Sigma }_R\left(k;t_1,t_2\right)$$

$$i\theta \left(t_1-t_2\right)\int \frac{d^{d}q}{(2\pi {)}^d}{e}^{i{\epsilon }_{q-k}(t_1-t_2)}i{D}_K\left(q,t,t\right).$$

$${\Sigma }_R^{WT}\left(k,\omega ;t\right)$$

$${\Sigma }^R\left(k,\omega ={\epsilon }_k,t\right)$$

$$\approx ZTt\int \frac{d^{d}q}{(2\pi {)}^d}\frac{F_K^0({\gamma }_{q}t)}{2{\epsilon }_k+\vec{q}\cdot {\vec{v}}_k+i\delta }$$

$${\Sigma }^R\left(k,\omega ={\epsilon }_k,t\right)$$

$$\displaystyle\hskip 10.0pt\int\!\!d\hat{n}\frac{F^{0}_{K}\!\left(tTl^{2}q^{2}\right)}{q\vec{v}_{k}\cdot\hat{n}+2\varepsilon_{k}+i\delta}\bigg{]},$$

$${\Sigma }^R\left(k,\omega ={\epsilon }_k,t\right)$$

$$\times {\int }_0^{\sqrt{Tt}}\frac{y^{d-1}dy}{(2\pi {)}^d}\int d\hat{n}\frac{F_K^0\left(y^2\right)}{y\cos {\theta }^{\prime }+2c{\epsilon }_k\sqrt{t/T}+i\delta }$$

Full text

I Introduction

Experiments involving pump-probe spectroscopy of solid-state systems Fausti et al. (2011); Smallwood et al. (2012, 2014); Beck et al. (2013); Mitrano et al. (2016) as well as dynamics of cold-atomic gases Regal et al. (2004); Zwierlein et al. (2004); Bloch et al. (2008); Endres et al. (2012); Navon et al. (2016) have opened up an entirely new temporal regime for probing correlated systems. A strong disturbance at an initial time, either by a pump beam of a few femto-second duration, or by explicit changes to the lattice and interaction parameters for atoms confined in an optical lattice, places the system in a highly nonequilibrium state. Following this, the dynamics can be probed with pico-second resolution for solid state systems, and milli-second resolution for cold atoms. At these short time-scales, the system is far from thermal equilibrium, and has memory of the initial pulse.

For generic systems, with only a finite number of conservation laws such as energy and particle number, relaxation to thermal equilibrium is expected to be fast Juchem et al. (2004); Eckstein et al. (2009); Santos and Rigol (2010); Tavora et al. (2014). However the relaxation may still possess rich dynamics when tuned near a critical point. As in equilibrium, observables will have a singular dependence on the detuning from the critical point. Such singular behavior in quench dynamics has already been identified for bosonic systems coupled to a bath Janssen et al. (1989); Huse (1989); Calabrese and Gambassi (2005); Gagel et al. (2014, 2015), and also for isolated bosonic systems Chandran et al. (2013); Chiocchetta et al. (2015); Maraga et al. (2015); Nicklas et al. (2015); Chiocchetta et al. (2016); Lemonik and Mitra (2016); Karl et al. (2017); Chiocchetta et al. (2017). In this paper, we complement this study with that for an isolated fermionic system.

We consider a gas of fermions in a lattice without disorder at finite temperature, where the fermions have an attractive interaction. This can be realized as a gas of cold atoms in an optical lattice. In equilibrium, there is a phase transition separating superfluid and normal (non-superfluid) phases. We study this phase transition as a quench process. In particular we study the dynamics when the system is initially in the normal phase, but where the interaction is suddenly increased so the system is close to the phase transition.

Our results may be understood as follows. The distance from the critical point is captured by the superconducting fluctuation $D(q,t)\equiv\langle\Delta^{\dagger}(q,t)\Delta(q,t)\rangle$, where $\Delta^{\dagger}(q,0)$ creates a Cooper pair at momentum $q$. In the initial system, deep in the normal phase, $D(q)$ is small and non-singular as a function of $q$. At the finite temperature equilibrium critical point, we expect $D(q,\infty)\sim q^{-2}$ as $q\rightarrow 0$. When a large interaction is turned on in the normal state, the system must interpolate between these two limits as a function of time. Since in a diffusive system information can only be communicated a distance $\propto\sqrt{t}$ in a time $t$ after the quench, we should expect a scale $\sim t^{-1/2}$ to function as the cutoff on the divergence of $D(q,t)$. Therefore the long wavelength fluctuations should show strong dependence on $t$.

We propose to detect these fluctuations through their effect on the spectral properties of the fermions, in the vein of fluctuation superconductivity. If we imagine a fixed background of superfluid fluctuations, the fermions can be thought of as Andreev reflecting off of this background. As the system is disordered the Andreev reflection is incoherent, and no gap is opened in the fermionic spectrum. However, the Andreev reflection is an energy conserving process for fermions at the chemical potential. This then contributes a decay channel for fermions at the Fermi surface which we can estimate by Fermi’s golden rule as $\int d^{d-1}{q}D(q,t)$ (where the $d-1$ dimensional integral corresponds to the Fermi surface of a $d$ dimensional system). In $d=2,3$ this integral is singular as $q\rightarrow 0,t\rightarrow\infty$, and we obtain a singular correction which goes like $\sqrt{t}$ in $d=2$ and $\log t$ in $d=3$. This behavior is summarized in Fig. 1.

Our analysis is in a regime complementary to studies such as Yuzbashyan et al. (2015); Liao and Foster (2015); Foster et al. (2014) where the initial state was already superconducting to begin with, and the dynamics of the superconducting order-parameter under an interaction quench was studied within a mean field approximation. It is also complementary to studies such as Sentef et al. (2016); Knap et al. (2016); Dehghani and Mitra (shed); Kennes et al. (2017) where the initial state was in the normal phase, and the external perturbation puts the system deep in the ordered phase where the dynamics were again studied in mean field. In contrast, we study the behavior when the system is always in the disordered phase and the mean field behavior is trivial. To controllably go beyond mean-field calculation we perform a $1/N$ expansion, where $N$ is an additional orbital degree of freedom for the fermions.

The paper is organized follows. In Section II we present the model, outline the approximations, and introduce an auxiliary or Hubbard-Stratonovich field that represents the Cooper pair fluctuations or Cooperons. In Section III the equations of motions are analyzed assuming the Cooperons are non-interacting, an assumption valid at short times. Further the effect of the fluctuations on the fermion lifetime is calculated. In Section IV, the longer time behavior is considered. This is done by mapping the dynamics to model-A Hohenberg and Halperin (1977), and solving the self-consistent equation for the self-energy. We conclude in section V. The mapping to model A dynamics is outlined in Appendices A and B, where App. B includes a perturbative estimate of the parameters of Model A. The implication of model A on Cooperon dynamics is relegated to App. C.

II Model

We study a quench where the initial Hamiltonian is that of free fermions,

[TABLE]

Above $k$ is the momentum, $\sigma=\uparrow,\downarrow$ denotes the spin, and $\tau$ is an orbital quantum number that takes $N$ values. We consider the initial state to be the ground state of $H_{i}$ at non-zero temperature $T$, and chemical potential $\mu$. The time-evolution from $t>0$ is in the presence of a weak pairing interaction $u$. We write the quartic pairing interaction in terms of pair operators $\Delta_{q}$ such that,

[TABLE]

The Hamiltonian above assumes contact interaction, so that only fermions with opposite spin quantum numbers scatter off of each other. In the superfluid phase $\langle\Delta_{q}\rangle\neq 0$. In this paper, since we are always in the normal phase, $\langle\Delta_{q}\rangle=0$.

In the strict $N\rightarrow\infty$ limit, the fermion number at each momentum $k$ is conserved. Thus the system is fully integrable, and fails to thermalize. We do not work in this limit and consider finite $1/N$ corrections to the behavior. In particular this allows us to study the back-reaction of the fluctuating Cooper pairs on the fermionic gas.

We use a two-particle irreducible ($2PI$) formalism Cornwall et al. (1974); Berges (2004) to obtain the equations of motions. The main ingredient is the sum of $2PI$ diagrams $\Gamma^{\prime}[G]$, which is a functional of the fermions Greens functions $G$,

[TABLE]

We here and generally suppress the spin and orbital indices, use numbers to indicate spacetime coordinates, and do not define the advanced function since it is the hermitian conjugate of the retarded part: $G_{A}(1,2)=G_{R}^{*}(2,1)$.

The Green’s function and interaction vertex correspond to the diagrams in Fig. 2, while at ${\cal O}(1/N)$, the Keldysh functional $\Gamma^{\prime}$ is the set of fermion loops shown in Fig. 3. Note the effect of the $1/N$ expansion is to select the Cooper interaction channel which is the most singular channel near the critical point.

The Green’s function is determined by the Dyson equation

[TABLE]

evaluated to $\mathcal{O}\left(1/N\right)$. The Cooperon can be understood as the correlator of an auxiliary or Hubbard-Stratonovich field $\phi$ conjugate to $\Delta$, used to decouple the fermionic quartic interaction in the Cooper channel, as outlined in Appendix A. In this language $D$ is defined by,

[TABLE]

We emphasize that Eqs. (4a-g), constitute a highly non-trivial set of coupled equations as $D$, $G$, and $\Sigma$ are defined self-consistently in terms of each other. In the rest of this paper we solve these equations in the nonequilibrium system following the quench. In Sec. III we do this in a short time approximation, which gives the essential qualitative behavior. In Sec. IV we remove the short time limit and consider the general behavior.

III Perturbative Regime

We begin by evaluating Eqs. (4) in the spirit of a short time approximation. We do this by replacing $G$ with it’s initial, noninteracting value $g$. The functions $D$ and $\Pi$ may then be straightforwardly obtained in terms of the dispersion $\epsilon_{k}$ and the initial occupation. At finite temperature $T$, and low frequency ($\omega/T\ll 1$) (see Appendix B) these can be estimated as,

[TABLE]

where $v_{F}$ is the Fermi velocity, $\nu$ the density of states and $a$, $b$ and $c$ system-dependent dimensionless constants. Then Eq. (4d) reduces to

[TABLE]

where we have $l=cv_{F}/T$ and $Z\sim T/\nu$, and $r$ is the distance from the critical point in units of $T$.

[TABLE]

The $D$ are overdamped as a consequence of the fermionic bath. When $r>0$ ($r<0$) $D_{R}$ decays (grows) with time indicating that the system is in the disordered phase (unstable to the ordered phase). The critical point $r=0$ separates the two regimes.

While $D_{R,A}$ are time translation invariant within the current approximation, $iD_{K}$ explicitly breaks time translation invariance. From Eqns. (4e), (7), and (9) it follows that

[TABLE]

As expected of a non-equilibrium system this violates the fluctuation dissipation theorem. We may quantify this by introducing a function

[TABLE]

Eq. (10) may be written as

[TABLE]

The violation the FDT (at the initial temperature T) is given by the fact that,

[TABLE]

Therefore this quantity measures the extent to which the Cooper pair fluctuations are out of equilibrium with the fermions.

At the critical point $r=0$, $iD_{K}\!\left(q,t,t\right)$ for $v_{F}q\ll T$ can be written in the scaling form

[TABLE]

This is a consequence of the fact that at the critical point the only length scale greater than $l$ is the one generated from $t$, $l\sqrt{Tt}$.

III.1 Fermion lifetime

We now show that the growing fluctuations may be detected through the spectral properties of the fermions, in particular the lifetime. There is some subtlety with defining the lifetime, as the system is not translationally invariant and so the response functions may not be decomposed in frequency space. However, we may take advantage of the fact the rate of change of $D_{K}$ is $\gamma_{q}$, while the typical energy of a thermal fermion is $T\gg\gamma_{q}$. Therefore it is reasonable to interpret the Wigner-Transform of the self energy,

[TABLE]

as being the self energy at $\omega$ near time $t$, and the quantity,

[TABLE]

as the fermion lifetime. In any case, the correct observable will be determined by the particular experimental protocol.

We now proceed to evaluate $\tau^{-1}$ within the present approximation. The self energy is determined from equation (4c) which in momentum space takes the form,

[TABLE]

We may now make several simplifications. First as $D_{K}/D_{R}\sim T/\gamma_{q}$ we may neglect the first term. Second, as the evaluation of $D_{K}(q;t_{1},t_{2})$ has shown that it varies on the scale of $\gamma_{q}$ which is much less then $T$, it is sufficient to replace $D_{K}\!\left(q,t_{1},t_{2}\right)$ with the equal time quantity $D_{K}(q;t,t)$, $t=(t_{1}+t_{2})/2$ evaluated at the average value of the two time coordinates. Third, continuing within the perturbative approximation, $G_{R}$ may be replaced with it’s non-interacting value.

Therefore we obtain the equation,

[TABLE]

Finally, we Fourier transform with respect to the time difference $t_{1}-t_{2}$,

[TABLE]

Setting $\omega=\varepsilon_{k}$, and using Eq. (14), we obtain,

[TABLE]

Focusing on the region where $ql<1$, we may set $\gamma_{q}=Tl^{2}q^{2}$. Going over to spherical coordinates this may be written as,

[TABLE]

where $\int d\hat{n}$ indicates the integral over the $d$-dimensional sphere. Now introducing new coordinates $y=ql\sqrt{Tt}$, produces

[TABLE]

[TABLE]

where $\theta^{\prime}$ is the angle between $\vec{v}_{k}$ and $\hat{n}$. Therefore the lifetime is given by

[TABLE]

We now evaluate this function, beginning in $d=2$. Recalling that $F^{0}_{K}(y)\rightarrow y^{-1}$ as $y\rightarrow\infty$ and $F^{0}_{K}(0)=1$, we see that in $d=2$ the integral is convergent and so we may neglect the condition that $q\ll l^{-1}$. Therefore $S^{0}_{2}(x)$ is not sensitive to how the integral is cutoff at $q\approx l^{-1}$, and is plotted in Fig. 5. The asymptotics may be extracted from the asymptotics of $F^{0}_{K}$,

[TABLE]

For completeness we give the asymptotic forms of Re$S_{2}^{0}$ which take the form.

[TABLE]

As a result, for $d=2$,

[TABLE]

The full curve is plotted in Fig. 5.

In $d=3$ the integral is logarithmically divergent. The upper limit of the integral over $y$ is $\sqrt{tT}$ whereas the divergence as $y\rightarrow 0$ is cutoff by the greater of $x$ and one. Thus, $S^{0}_{3}(x)\sim\log\left[\sqrt{Tt}/{\rm max}(1,x)\right]$, and

[TABLE]

For completeness we give asymptotic limits of Re$S_{3}^{0}$ which takes the form.

[TABLE]

Therefore for $d=3$, we have the dependence,

[TABLE]

The full behavior of the function $S^{0}_{3}$ is plotted in Fig. 6. Note that as this function is logarithmically dependent on the cutoff of the integral, a change in the precise form of the cutoff enforcing $q\ell\ll 1$ will shift the final result by an additive constant. This ambiguity would be fixed by comparing the asymptotic behavior of $\tau(\varepsilon)$, $\varepsilon\rightarrow\infty$ with a microscopic calculation of the lifetime at high energies.

The self-energy at $\varepsilon=0$ diverges as $t\rightarrow\infty$. Therefore at sufficiently long time the assumption that $G$ may be substituted with it’s non-interacting value is not valid. In the next section we lift this assumption.

IV Non-Perturbative Regime

We note that the $2PI$ formalism is not a short-time expansion, and that the Eqns. (4) are valid at all times. There are two shortcomings that must be remedied.

First, we have seen that the $D$ is a slow function of time, and that $\Sigma$ depends essentially on the equal time value $D_{K}(t,t)$. The estimation in Sec. III essentially estimates $D_{K}(t,t)$ by linear response. However as $D$ increases with time, it will eventually grow large enough to violate the linear response assumption. We remedy this by employing the theory of critical quenches in this section.

Second, we have that $G$ and $\Sigma$ must be self-consistent. We therefore solve the self-consistent version of Eq. (4). The results of this analysis are shown in Fig. 1.

IV.1 Propagator for interacting Cooperons

The equation of motion for $D$ given in Eq. (8), is equivalent to linear response. This is because the non-linear behavior comes from the dependence of $G$ on $\Sigma$, which in turns depends on $D$. Therefore once the perturbative approximation for $G$ fails it becomes necessary to consider the non-linear evolution of $D$.

Since we are only concerned with the long wavelength behavior of $D$ in the vicinity of the dynamical critical point, we do not need to solve for the full behavior of $D$. Instead we need only to identify the appropriate dynamical universality class. As the fluctuations of the bosonic field $\Delta$ are not conserved and are overdamped by the fermionic bath, the system belongs to the dynamical model-A transition Hohenberg and Halperin (1977). The standard manipulation mapping the original Hamiltonian to this model are relegated to Appendices A, B.

We quote the results for the transient dynamics of thermal aging in model A already discussed elsewhere Janssen et al. (1989); Calabrese and Gambassi (2005), and re-derived in Appendix C.

In $d=2$ there is no true dynamical critical point, as in equilibrium. Therefore our results are only valid in the perturbative regime in $d=2$. The behavior at intermediate to long time is expected to be described by a crossover to Kosterlitz-Thouless physics, where the amplitude of the Cooper fluctuations saturate but there are long range phase fluctuations. However, such a calculation is beyond the scope of this paper.

In $d=3$, the dynamics is characterized by three exponents $z,\eta,\theta$. Of these $z,\eta$ are already familiar in equilibrium $\phi^{4}$ theory and are the dynamical critical exponent and the scaling dimension of $\phi$ respectively.

The exponent $\theta$ is a non-equilibrium exponent known as the initial slip exponent Janssen et al. (1989). It is responsible for non-trivial aging dynamics, and is interpreted as the scaling dimension of a source field applied at short times after the quench. This is because such a source field will induce an initial order-parameter $M_{0}=\langle\phi_{c}\rangle$, to grow with time at short times after the quench as Janssen et al. (1989) $M_{0}\sim t^{\theta}$ even though the quench is still within the disordered phase. At long times, eventually the order-parameter will decay to zero. For the present calculation, the short time behavior is not directly relevant as we are interested in the regime when $|t-t^{\prime}|\ll t$. However the short-time exponent still affects the qualitative behavior of $D_{K}$.

The values of the Model A exponents $z$, $\eta$ and $\theta$ may be calculated using standard methods such as epsilon-expansions or large-N, where $N$ now controls the components of the bosonic field. We adopt the latter approach, and emphasize that this component $N$ is not the same as the fermion orbital index used to justify the form of $\Gamma^{\prime}$. The derivation of the exponent using large-N for the bosonic theory is equivalent to a Hartree-Fock approximation for model A (see App. C) giving $z=2$, $\eta=0$ and $\theta=\epsilon/4,\epsilon=4-d$. Other approximations will change the precise value of the exponent, but not the overall scaling form.

The results for the Cooperon dynamics from model A (see App. C) are as follows. When $t^{\prime}$ becomes comparable to $t$, we expect,

[TABLE]

In particular for equal times,

[TABLE]

The above form for the boson density $iD_{K}(q,t,t)$ is the same as that in Sec. III if one replaces the scaling function $F_{K}^{0}$ by $F_{K}$. The function $F_{K}$ has the asymptotic limits

[TABLE]

Note the leading asymptotic behavior of $F_{K}(x)$ as $x\rightarrow\infty$ is the same as that of $F_{K}^{0}(x)$.

We define a scaling function $S_{d}$, $d>2$ as the analogue of Eq. (23) for interacting bosons,

[TABLE]

This replacement makes only a small quantitative change to the final result in $d=3$. As the leading behavior at $F_{K}$ is unchanged, the derivation of the asymptotics Sec.  III may be followed precisely, leading to

[TABLE]

The only difference in the asymptotic behavior between ${\rm Im}S_{3}^{0}$ and ${\rm Im}S_{3}$ might be in the constants. However $\tau^{-1}$ only depends on ${\rm Im}S_{3}\left(\beta\right)$ where $\beta$ is a material dependent parameter, see Eq. (24). Further as ${\rm Im}S_{3}$ has a logarithmic dependence on the cutoff, changing the (material-dependent) cutoff shifts ${\rm Im}S_{3}(x)\rightarrow{\rm Im}S_{3}(x)+\gamma$. Thus the constants in Eq. (35) may be absorbed into constant $\beta,\gamma$.

Once these are fixed, ${\rm Im}S_{3}(x)$ and ${\rm Im}S_{3}^{0}(x)$ both crossover smoothly between the same asymptotics and therefore it is reasonable that they are qualitatively similar, see Fig. 7. As the final result does not appear significantly sensitive to the critical exponents we will not attempt to estimate the value of these exponents more accurately.

IV.2 Self-consistent solution

Although the intermediate regime is only controlled for the case of $d=3$, we include the self-consistent equation in both $d=2,3$ for completeness. Having considered the behavior of $D$ we may now directly solve the self consistent equation for the self-energy, Eq. (4). In Wigner coordinates, it is given by,

[TABLE]

Expanding the denominator in $q$ to first order, and assuming that the variation in $\Sigma^{R}$ with $q$ is negligible, gives

[TABLE]

Now the scaling for $D_{K}$ is

[TABLE]

To take advantage of this we rescale units by defining

[TABLE]

Giving

[TABLE]

The integral over $\hat{n}$ is an integral over unit vectors in $\mathbb{R}^{d}$. The condition that $v|q|\ll T$ is imposed by cutting off the integral at $y_{m}\sim\sqrt{tT}$. So, it goes to infinity as $t\rightarrow\infty$. This gives a self-consistent equation for $z$. The function $F_{K}(w)$ goes to a constant as $w\rightarrow 0$, decays like $1/w$ as $w\rightarrow\infty$. We introduce the function $S_{d}$ depending on dimension $d$ so that we can write,

[TABLE]

Since we need to solve the integral self-consistently we must understand how $S_{d}(z)$ behaves for $z$ in the upper half of the complex plane. This is greatly simplified since $S_{d}$ is analytic as a function of $z$ in the upper half complex plane, as the only singularity can come from the pole $y\hat{k}\cdot\hat{n}=-z$. The estimate of $S_{d}$ depends on the dimensions.

IV.2.1 $d=3$

In $d=3$, Eq. (43) leads to

[TABLE]

Recalling the position of the branch cut as $z\rightarrow i\delta$ we obtain

[TABLE]

so that as $|z|\rightarrow 0$ we get the leading behavior as $y_{m}\rightarrow\infty$,

[TABLE]

where $\zeta$ is an order one constant. The integral does not converge as $y_{m}$ goes to $\infty$. Therefore, $S_{3}$ does not depend only on the variable $z$ but also on $y_{m}$ and therefore exactly how the integral is cutoff at $q\approx q_{m}$. In particular by shifting $y_{m}$ to a new value $y^{\prime}_{m}$ changes $S_{3}\rightarrow S_{3}+i\log(y_{m}/y^{\prime}_{m})/4\pi$. Therefore the imaginary part of $S_{3}$ is ambiguous up to an overall additive constant.

To understand the large $z$ behavior, we split the integral into the regions $y\ll\zeta$ and $y\gg\zeta$, where $\zeta$ is some constant of order one. The small $y$ limit is

[TABLE]

And the large $y$ limit is

[TABLE]

As $z\rightarrow\infty$ this diverges logarithmically around $u=0$, therefore the integral is approximately $-\pi i\log\left(\zeta/z\right)$. Collecting the results we have that,

[TABLE]

The results are summarized in Fig. 6. Note the the substitution of $F_{K}$ for $F_{K}^{0}$, makes minimal difference in the calculation of $S_{3}$, see Fig. 7.

Returning to the self consistent equation

[TABLE]

If we assume that $z\approx z_{0}$, we obtain

[TABLE]

Plugging this back into the self consistent equation

[TABLE]

This implies the condition for validity of the perturbative solution is

[TABLE]

Substituting $z_{0}=2\varepsilon\sqrt{t/T}$, we see this condition is equivalent to,

[TABLE]

Therefore the short time dynamics is sufficient to explain the behavior of the tails of the distribution, which is reasonable as these saturate at short times.

Let us look for the self-consistent solution at $z_{0}=0$ and $\alpha(Tt)^{1/2}\gg 1$. Assuming $z\gg 1$ we get

[TABLE]

Bearing in mind the $y_{m}\propto\sqrt{t}$ we see that the above has a solution with $z\propto t^{1/2}$. Therefore the $\Sigma_{R}(0,0,t)$ saturates at a constant at long times, given by the equation

[TABLE]

To summarize, for $\varepsilon_{k}/T\gg\alpha$, the behavior is the same as in Sec. III, with saturation at $\log\varepsilon_{k}/T$. For smaller energies the logarithmic growth given earlier saturates at $Tt\sim\alpha^{-2}$. The general behavior is shown in Fig. 1.

The approximate behavior of $\tau^{-1}$ at $\varepsilon_{k}=0$ is shown in Fig. 8. Unfortunately calculating $S_{3}(z)$ over the upper half plane and solving Eq. (58) is numerically intensive. Instead we approximate

[TABLE]

which renders Eq. (58) analytically tractable. As this approximation has the same asymptotic limits as $S_{3}$ it should be sufficient for reproducing the qualitative shape of $\tau^{-1}$.

IV.2.2 $d=2$

We now analyze the self consistent equation in $d=2$.

[TABLE]

We estimate this integral as follows. First as $z\rightarrow 0+i\delta$ this goes to $i$ for some order one constant. Note the sign is determined by the branch cut and should be consistent with causality. On the other hand if $z\gg 1$ we split the integral at some $\zeta$ of order one:

[TABLE]

and the other half

[TABLE]

As $|z|\rightarrow\infty$ this is $\sim\log(2z/\zeta)/z$, which dominates the small $y$ contribution, and it’s effect is to renormalize the order one cutoff $\zeta$. So we may summarize the behavior as

[TABLE]

The real and imaginary parts of $S_{2}^{0}$ which is similar to $S_{2}$, are plotted in Fig. 5. We deal with the self consistent equation essentially as in $d=3$. The perturbative condition holds at large $z_{0}$,

[TABLE]

For $z_{0}=0$ this condition is always violated at the time scale $1/\alpha$. However, if we take $z_{0}=2\varepsilon\sqrt{t/T}\gg 1$, then using the asymptotics we estimate that

[TABLE]

The perturbative condition is only violated at an exponentially long time $t\propto\exp(\varepsilon^{2}/(\alpha T^{2})$ ).

We now seek a self-consistent solution when $\alpha Tt\gg 1$, but $z_{0}$ is small. We use the large $z$ asymptotics.

[TABLE]

Solving the quadratic equation treating $\log(z)$ as a constant we get,

[TABLE]

The choice of branch comes from matching the behavior as $\alpha\rightarrow 0$. In the regime of interest where $t\gg 1$, we can to good accuracy simply replace the $\log z$ on the RHS with $\log(\pi\alpha Tt)$.

Taking the $t\gg 1$ limit we obtain

[TABLE]

We see that $z\gg 1$ so the assumption of large $z$ is self-consistent.

Translating back to the self energy via $\Sigma_{R}=z\sqrt{T/t}$, we obtain

[TABLE]

The self energy apparently grows without bound at $z_{0}=0$ albeit extremely slowly. We interpret this unbounded growth as a symptom of the non-existence of the true critical point in $d=2$ and therefore the impossibility of a self-consistent treatment in this regime.

V Conclusions

In this paper we have analyzed the superfluid quench, wherein an attractive interaction is suddenly turned on in a normal fluid of fermions. This interaction enhances superfluid fluctuations. There are two regimes: a disordered phase at weak interaction strength where the fluctuations saturate at a finite value; and the ordered phase, for strong interaction strength, where the fluctuations grow exponentially, leading eventually to spontaneous symmetry breaking.

Between these two regimes is a dynamical critical point, where the fluctuations grow but order is not formed. We find that as with the usual equilibrium critical points, there is a notion of universality associated with this dynamical critical point. That is, once a small number of constants are fixed, the complete behavior of the superfluid fluctuations is determined by a function of the wavelength and time, with no further free parameters. The necessary parameters are the $r$ and $\ell$ given in Eq. (8).

Moreover, we find a signature of this universality in the lifetime of the fermions. The mechanism is essentially that the fermions near the Fermi energy scatter resonantly off of superfluid fluctuations. Thus the growing superfluid fluctuations lead to a singular feature in the fermion lifetime as a function of energy. We show that this singular feature inherits the universality of the dynamical critical point. In particular after fixing the Fermi velocity $v_{F}$ and normalized scattering rate $\alpha$, the energy and time dependence of the lifetime is completely determined.

The present work may be extended in several directions. One is the full development of the kinetic equation governing the fermion dynamics, to be published elsewhere. It would also be of interest to repeat this analysis for a disordered system to allow for comparison with pump-probe experiments. Lastly, extending this treatment to include other fermion symmetry breaking channels, such as magnetic orders, or charge-density waves, would be fruitful.

We note that the perturbative calculation in $d=3$ gives a logarithmic correction $\sim\log t$ which grows large with $t$. This suggests that a dynamical RG conducted around the critical dimension $d=3$ may be a fruitful alternative way to approach this problem.

In this paper we have consider the fermions to initially be at finite temperature before the quench. A natural problem would be to consider the quench starting with fermions at zero temperature. This problem is more delicate for at least two reasons. Firstly the superfluid phase transition always occurs at finite temperature, therefore to approach the critical regime one would have to consider the temperature that is dynamically generated by the self-heating of the fermions. Secondly before the temperature is generated, the fermions are controlled by quantum fluctuations, leading to complex prethermal dynamics Babadi et al. (2015). These difficulties aside, the problem appears deserving of future study.

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

Appendix A The $D$ propagator or Cooperon as correlators of Hubbard-Stratonovich fields

The final post-quench Hamiltonian is,

[TABLE]

We will highlight the meaning of $D$ in an imaginary time formalism as the generalization to real time Keldysh formalism is conceptually straightforward.

We may decouple the quartic interaction via a complex field $\phi_{q}$ for each momentum mode $q$,

[TABLE]

In this picture, the action is quadratic in the fermionic fields. After integrating out the fermions one may write the partition function $Z$ as,

[TABLE]

where $g^{-1}$ is the non-interacting fermionic Green’s function in $2\times 2$ Nambu space. On expanding the ${\rm Tr}\ln$, one obtains an action for the $\phi$ fields. Since the system is assumed to be in the normal phase, only even powers of the $\phi$ field enter the action. Thus, we obtain,

[TABLE]

where

[TABLE]

The Gaussian approximation involves keeping only quadratic terms in the $\phi$ fields. The coefficient of $\phi^{2}$ in the second term in the action is recognized as the polarization bubble $\Pi\equiv gg$. The equation of motion at Gaussian order is,

[TABLE]

where trace over the fermions gives an additional factor of $N$. In the next sub-section we show that Eq. (85) is equivalent to a classical Langevin equation for the Hubbard-Stratonivich fields when the fermions are at non-zero temperatures.

Appendix B Properties of the $\Pi$ and Relationship to Model-A

The fermionic distribution function before the quench is, $n_{\sigma}(k)=1/\left(e^{\frac{\xi_{k}}{T}}+1\right),\xi_{k}=\epsilon_{k}-\mu$. We measure all energies relative to the chemical potential.

The Keldysh component of the polarization bubble is found to be,

[TABLE]

while the retarded component is

[TABLE]

Since the $\Pi$ are time-translation invariant in this approximation, it is helpful to write them in frequency space,

[TABLE]

Now we use the fact that $1-2n(x)=\tanh(x/2)$ and using that $\coth(a)\coth(b)+1=\coth(a+b)(\coth(a)+\coth(b))$, one may show that fluctuation dissipation theorem (FDT) is obeyed,

[TABLE]

It should be emphasized that this FDT is simply inherited from the properties of the initial state. In a better approximation, the FDT will cease to hold as the system goes through the process of thermalization.

We are interested in the dynamics of the soft Cooperon mode, which evolves on a timescale much larger than $T^{-1}$. Therefore we expand $\Pi^{R}\!\left(q,\omega\right)$ in $\omega/T$, $q^{2}/T$. The constant term,

[TABLE]

is the usual Cooper logarithm, where $\nu$ is the density of states and $E_{F}$ is some bandwidth or Fermi energy. The coefficient of $\omega/T$ is

[TABLE]

which is purely imaginary in the absence of particle hole asymmetry.

For the coefficient of $q^{2}/T$, we expand the dispersion as $\epsilon_{k-q}=\epsilon_{k}-\vec{q}\cdot\vec{v}_{k}$ and obtain,

[TABLE]

where $\langle\cdot\rangle_{FS}$ is the average over the Fermi surface and the integral evaluates to the constant $28\zeta(3)/\pi^{2}\approx 3.41$.

With this expansion for $\Pi_{R}$, the FDT gives that $\Pi_{K}$ is given by

[TABLE]

and therefore in real time

[TABLE]

as discussed in the text. The fact that $\Pi_{K}$ is well approximated by a delta function is entirely a consequence of the fact that we are interested in timescales much longer than $T^{-1}$, because the Cooperon dynamics are governed by much longer timescales at the critical point.

Thus in summary, the above behavior for $\Pi$ together with how it affects the equation of motion of $D$ (Eq. (85)) show that the Cooperon obeys model-A dynamics close to the critical point.

Appendix C Interacting Cooperons in the Hartree-Fock Approximation

For the sake of completeness we outline how the results for interacting bosons used in the main text were obtained. We employ a Hartree-Fock approach, although the same scaling forms can be obtained with an $\epsilon$-expansion Janssen et al. (1989). The Hartree-Fock approximation for the bosons is justified as the $N\rightarrow\infty$ limit of a bosonic model where $N$ denotes the number of components of the boson field. This $N$ should not be confused with the orbital index of the fermions used in the main text.

The Hartree-Fock equations of motion are,

[TABLE]

where the mass obeys the equation of motion

[TABLE]

The overdamped dynamics of $D_{R}$ is entirely due to the underlying finite temperature Fermi sea which gives $u^{-1}-\Pi^{R}=r+i\omega$.

If we employ the Gaussian expression for $D_{K}(q,t,t)\rightarrow\frac{T}{q^{2}+r}\left[1-e^{-2(q^{2}+r)t}\right]$, we find that

[TABLE]

The above shows that scaling emerges only if we set $d=4$ in the above Gaussian result, showing that the upper critical dimension of the theory is $d=4$.

Thus with the ansatz,

[TABLE]

we obtain,

[TABLE]

For $D_{K}$ we have,

[TABLE]

For $q^{2}t^{\prime}\ll 1$, we obtain aging behavior,

[TABLE]

For equal times we may write,

[TABLE]

Note that $F(x=0)=0$ and $F(x=\infty)=1$.

In order to solve for $a$, we use that

[TABLE]

At criticality $r_{\rm eff}(\infty)=0$ and $iD_{K}(q,\infty,\infty)=T/q^{2}$. Using this,

[TABLE]

where $A_{d}=$ the surface area of a $d$-dimensional unit sphere.

The above may be recast as

[TABLE]

Thus we may write, defining $\epsilon=4-d$,

[TABLE]

The first integral above increases in time as $t^{\epsilon/2}$ unless

[TABLE]

Notice that in Eq. (104), to avoid infra-red singularity $2a<1$. Then,

[TABLE]

Since $F(0)=0$, we require $\epsilon/2<1$ or $d>2$ to make the $\int dxx^{-\epsilon/2}$ infra-red convergent. Substituting Eq. (110) in Eq. (109), we obtain

[TABLE]

Thus we have derived the quoted scaling forms, and also the initial slip exponent $a=\theta=\epsilon/4$.

It is also useful to note that for $t>t^{\prime}$ but general $q^{2}t,q^{2}t^{\prime}$, one obtains from Eq. (102),

[TABLE]

Defining

[TABLE]

[TABLE]