Dynamical generation of superconducting order of different symmetries in hexagonal lattices

TL;DR

This paper investigates how different superconducting orders develop dynamically in hexagonal lattices after an interaction quench, revealing a transition from d-wave to s-wave dominance as quench strength increases, with implications for transient measurements.

Contribution

It demonstrates the dynamical evolution of superconducting symmetries in hexagonal lattices, including the effects of time-reversal symmetry and quench amplitude, highlighting a crossover from d-wave to s-wave dominance.

Findings

d-wave orders grow fastest for small quenches

TR symmetry breaking favors chiral d-wave superconductivity

Large quenches favor s-wave superconductivity

Abstract

The growth of superconducting order after an interaction quench in a hexagonal lattice is studied. The cases of both time-reversal (TR) preserving graphene, as well as the TR broken Haldane model are explored. Spin singlet superconducting order is studied where the $s$, $d+id$, and $d-id$ wave orders are the irreducible representations of the hexagonal lattice. For small quenches, the $d$-wave order parameter grows the fastest, a result also expected when the system is in thermal equilibrium. For the TR symmetry preserving case, the growth rate of the two $d$-wave orders is identical, while the TR-broken case prefers one of the chiral $d$-wave orders over the other, leading to a TR broken topological superconductor. As the interaction quench becomes larger, a smooth crossover is found where eventually the growth rate of the $s$-wave becomes the largest. Thus for large interaction quenches, the $s$-wave is preferred over the $d$-wave for both TR preserving and TR broken systems. This result is explained in terms of the high energy quasi-particles responsible for the dynamics as the interaction quench amplitude grows. The results are relevant for time-resolved measurements that can probe the symmetry of the superconducting fluctuations in a transient regime.

Figures

Equations

$$H_0(t)=-t_h\sum _{k\sigma ,\alpha =1,2,3}\left(\begin{array}{cc}{a}_{k\sigma }^{†}& b_{k\sigma }^{†}\end{array}\right)$$

$$\times \left(\begin{array}{cc}0& e^{i\vec{k}\cdot {\vec{a}}_{\alpha }+ia\vec{A}(t)\cdot {\vec{\delta }}_{\alpha }}=h_{ab}\\ e^{-i\vec{k}\cdot {\vec{a}}_{\alpha }-ia\vec{A}(t)\cdot {\vec{\delta }}_{\alpha }}=h_{ba}& 0\end{array}\right)$$

$$\times \left(\begin{array}{c}{a}_{k\sigma }\\ b_{k\sigma }\end{array}\right).$$

$${\vec{\delta }}_1=a\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right);{\vec{\delta }}_2=a\left(\frac{1}{2},-\frac{\sqrt{3}}{2}\right);{\vec{\delta }}_3=a\left(-1,0\right).$$

$${\vec{a}}_1={\vec{\delta }}_1-{\vec{\delta }}_3=a\left(3/2,\sqrt{3}/2\right),$$

$${\vec{a}}_2={\vec{\delta }}_2-{\vec{\delta }}_3=a\left(3/2,-\sqrt{3}/2\right),$$

$${\vec{a}}_3=0.$$

$$\displaystyle V=J\sum_{<ij>}\biggl{[}\vec{S}_{i}\cdot\vec{S}_{j}-\frac{1}{4}n_{i}n_{j}\biggr{]},$$

$$h_{ab}^m=\frac{1}{T}{\int }_0^{T}dt{e}^{-im\Omega t}{h}_{ab}(t)$$

$$=i^{-m}{J}_{-m}(A_{0}a)\sum _{j=1,2,3}{e}^{i\vec{k}\cdot {\vec{a}}_j+im{\alpha }_j},$$

$$h_{ba}^m=\frac{1}{T}{\int }_0^{T}dt{e}^{-im\Omega t}{h}_{ba}(t)$$

$$=(-i{)}^{-m}{J}_{-m}(A_{0}a)\sum _{j=1,2,3}{e}^{-i\vec{k}\cdot {\vec{a}}_j+im{\alpha }_j}.$$

$$H_{\mathrm{eff}}={\hat{H}}_{m=0}+\sum _{m\ne 0}\frac{{\hat{H}}_m{\hat{H}}_{-m}}{m\hslash \Omega }+\dots .$$

$$H_{\mathrm{eff}}=-t_h\sum _{k\sigma ,\alpha =1,2,3}\left(\begin{array}{cc}{a}_{k\sigma }^{†}& b_{k\sigma }^{†}\end{array}\right)$$

$$\times \left(\begin{array}{cc}{t}_h{\sum }_{m\ne 0}\frac{h_{ab}^m{h}_{ba}^{-m}}{\hslash m\Omega }& h_{ab}^{m=0}\\ h_{ba}^{m=0}& t_h{\sum }_{m\ne 0}\frac{h_{ba}^m{h}_{ab}^{-m}}{\hslash m\Omega }\end{array}\right)$$

$$\times \left(\begin{array}{c}{a}_{k\sigma }\\ b_{k\sigma }\end{array}\right).$$

$$\displaystyle H_{\rm eff}=-\sum_{k\sigma}\begin{pmatrix}a_{k\sigma}^{\dagger}&b_{k\sigma}^{\dagger}\end{pmatrix}\begin{pmatrix}2t_{2}\biggl{[}\sum_{ij=12,23,31}\sin(\vec{k}\cdot(\vec{a}_{i}-\vec{a}_{j}))\biggr{]}&t_{1}\sum_{j=1,2,3}e^{i\vec{k}\cdot\vec{a}_{j}}\\ t_{1}\sum_{j=1,2,3}e^{-i\vec{k}\cdot\vec{a}_{j}}&-2t_{2}\biggl{[}\sum_{ij=12,23,31}\sin(\vec{k}\cdot(\vec{a}_{i}-\vec{a}_{j}))\biggr{]}\end{pmatrix}\begin{pmatrix}a_{k\sigma}\\ b_{k\sigma}\end{pmatrix},$$

$$\displaystyle t_{1}=t_{h}J_{0}(A_{0}a),\,\,\,\,\,\,\,\,t_{2}=2t_{h}^{2}\sum_{m=1\ldots\infty}\frac{J_{m}^{2}(A_{0}a)}{\hbar m\Omega}\sin\biggl{[}\frac{2\pi m}{3}\biggr{]}.$$

$$V=-J\sum _{<ij>}{h}_{ij}^{†}{h}_{ij},$$

$$\displaystyle h_{ij}^{\dagger}=\frac{1}{\sqrt{2}}\biggl{(}a_{i\uparrow}^{\dagger}b_{j\downarrow}^{\dagger}-a_{i\downarrow}^{\dagger}b_{j\uparrow}^{\dagger}\biggr{)},$$

$$\displaystyle\phi_{k}={\rm Arg}\biggl{[}\sum_{\alpha}e^{i\vec{k}\cdot\vec{a}_{\alpha}}\biggr{]},$$

$$H_{\mathrm{eff}}=\sum _{k\sigma }\left(\begin{array}{cc}{a}_{k\sigma }^{†}& b_{k\sigma }^{†}\end{array}\right)$$

$$\times d_k\left(\begin{array}{cc}\cos {\theta }_k& e^{i{\varphi }_k}\sin {\theta }_k\\ e^{-i{\varphi }_k}\sin {\theta }_k& -\cos {\theta }_k\end{array}\right)\left(\begin{array}{c}{a}_{k\sigma }\\ b_{k\sigma }\end{array}\right).$$

$$|\tau =\pm {⟩}_k=\left(\begin{array}{c}{u}_{k\pm }\\ v_{k\pm }\end{array}\right),$$

$$|\tau =+{⟩}_k=\left(\begin{array}{c}\cos \frac{{\theta }_k}{2}\\ e^{-i{\varphi }_k}\sin \frac{{\theta }_k}{2}\end{array}\right);|\tau =-{⟩}_k=\left(\begin{array}{c}-\sin \frac{{\theta }_k}{2}\\ e^{-i{\varphi }_k}\cos \frac{{\theta }_k}{2}\end{array}\right).$$

$$\displaystyle\Delta_{\alpha}=J\biggl{\langle}b_{i+\alpha\downarrow}a_{i\uparrow}-b_{i+\alpha\uparrow}a_{i\downarrow}\biggr{\rangle},$$

$$H(t)=H_{\mathrm{eff}}$$

$$\displaystyle-\sum_{k,\alpha}\biggl{[}\Delta_{\alpha}(t)e^{i\vec{k}\cdot{\vec{a}}_{\alpha}}\biggl{(}a_{k\uparrow}^{\dagger}b^{\dagger}_{-k\downarrow}-a_{k\downarrow}^{\dagger}b^{\dagger}_{-k\uparrow}\biggr{)}+h.c.\biggr{]}.$$

$$\displaystyle\Delta_{\alpha}(t)=J\sum_{k}e^{-i\vec{k}\cdot{\vec{a}}_{\alpha}}\biggl{\langle}b_{-k\downarrow}(t)a_{k\uparrow}(t)-b_{-k\uparrow}(t)a_{k\downarrow}(t)\biggr{\rangle},$$

$${\Delta }_{\alpha }(t)=-iJ\sum _{k{k}^{\prime },\beta }{\int }_0^{t}d{t}^{\prime }{\Delta }_{\beta }(t^{\prime })e^{-i\vec{k}\cdot {\vec{a}}_{\alpha }+i\overrightarrow{k^{\prime }}\cdot {\vec{a}}_{\beta }}$$

$$\displaystyle\times\biggl{\langle}\biggl{[}\biggl{(}a_{k^{\prime}\uparrow}^{\dagger}(t^{\prime})b^{\dagger}_{-k^{\prime}\downarrow}(t^{\prime})-a_{k^{\prime}\downarrow}^{\dagger}(t^{\prime})b^{\dagger}_{-k^{\prime}\uparrow}(t^{\prime})\biggr{)},$$

$$\displaystyle\biggl{(}b_{-k\downarrow}(t)a_{k\uparrow}(t)-b_{-k\uparrow}(t)a_{k\downarrow}(t)\biggr{)}\biggr{]}\biggr{\rangle}.$$

$${\Delta }_{\alpha }(t)=J\sum _{\beta }{\int }_0^{t}d{t}^{\prime }{\Pi }_{\alpha \beta }^R(q=0,t,t^{\prime }){\Delta }_{\beta }(t^{\prime }).$$

$$\displaystyle\Pi^{R}_{\alpha\beta}(q,t,t^{\prime})=-2i\theta(t-t^{\prime})\sum_{k,\tau,\tau_{1}=\pm}e^{-i\epsilon_{k+q\tau}(t-t^{\prime})}e^{-i\epsilon_{-k\tau_{1}}(t-t^{\prime})}\biggl{\{}1-\rho_{k+q,\tau}-\rho_{-k,\tau_{1}}\biggr{\}}\biggl{[}e^{-i\vec{k}\cdot({\vec{a}}_{\alpha}-{\vec{a}}_{\beta})}|u_{k+q,\tau}|^{2}|v_{-k\tau_{1}}|^{2}$$

$$\displaystyle+e^{-i\vec{k}\cdot({\vec{a}}_{\alpha}+{\vec{a}}_{\beta})-i\vec{q}\cdot{\vec{a}}_{\beta}}u_{k+q\tau}v_{k+q\tau}^{*}v_{-k\tau_{1}}u^{*}_{-k\tau_{1}}\biggr{]},$$

$$1-{\rho }_{k,\tau }-{\rho }_{-k,{\tau }_1}={\delta }_{\tau ,{\tau }_1}\tanh \left(\frac{{ϵ}_{k,\tau }+\mu }{2T}\right)$$

$$\displaystyle+\delta_{\tau,-\tau_{1}}\frac{1}{2}\biggl{[}\tanh\left(\frac{\epsilon_{k,\tau}+\mu}{2T}\right)-\tanh\left(\frac{\epsilon_{k,\tau}-\mu}{2T}\right)\biggr{]}.$$

$$\displaystyle\delta=\frac{1}{2}\sum_{k}\biggl{[}\tanh\biggl{(}\frac{\epsilon_{k}+\mu}{2T}\biggr{)}-\tanh\biggl{(}\frac{\epsilon_{k}-\mu}{2T}\biggr{)}\biggr{]}.$$

$${\Pi }_{\alpha \beta }^R(q,\omega =0)=\sum _{k,\tau ,{\tau }_1=\pm }\frac{1}{{ϵ}_{k+q\tau }+{ϵ}_{-k{\tau }_1}}$$

Full text

I Introduction

The dynamics of how superconducting order develops, and the possibility of enhancing it under nonequilibrium conditions, has become an active area of research spanning solid-state systems Fausti et al. (2011); Graf et al. (2011); Smallwood et al. (2012, 2014); Beck et al. (2013); Mitrano et al. (2016) and cold-atomic gases Regal et al. (2004); Zwierlein et al. (2004); Bloch et al. (2008); Endres et al. (2012). These modern experiments have been complemented by intense theoretical activity attempting to understand dynamics of superconducting systems following an interaction quench Yin and Radzihovsky (2016); Yuzbashyan et al. (2015); Liao and Foster (2015); Foster et al. (2014) and also the possibility of driving a normal system superconducting via nonequilibrium lattice vibrations Sentef et al. (2016); Knap et al. (2016); Kennes et al. (2017); Sentef (2017).

Most theoretical studies of nonequilibrium superconductivity, with a few exceptions Peronaci et al. (2015); Sentef et al. (2017), assume that only one kind of order is relevant. However realistic systems allow for several competing superconducting orders corresponding to different irreducible representations of the crystal point group. In this paper we therefore ask the following question: is the symmetry of the superconducting order developing under nonequilibrium conditions, such as after a quench, the same or different from the one favored by the system in thermal equilibrium?

In this paper, we identify the dominant pairing symmetry following an interaction quench in a hexagonal lattice, both in the presence and absence of time-reversal (TR) symmetry. For the TR-preserving system we study doped graphene, while for the TR-broken system we consider doped graphene subjected to a high frequency circularly polarized laser. We assume that the system before the interaction quench is in the normal state, and then study dynamics of an infinitesimal superconducting fluctuation following the sudden switch on of an attractive interaction. The time-evolution of the initial superconducting fluctuation is studied within linear response, the goal being to simply identify the most unstable mode. Eventually heating effects arising due to the energy injected in the quench, and also due to the laser, could prevent true order from developing at the longest times. We do not address this long time behavior here. Thus our theoretical treatment assumes collisionless or prethermal dynamics. We are motivated by pump-probe experiments that clearly show that prethermal regimes exist, and that these transient regimes may be unstable to ordered phases different from that in thermal equilibrium Fausti et al. (2011); Mitrano et al. (2016).

Graphene subjected to a high frequency (as compared to the bandwidth) and circularly polarized laser has been shown to be equivalent to the Haldane model Oka and Aoki (2009); Kitagawa et al. (2011); Haldane (1988). In fact deviations from the Haldane model are found to be small even when accounting for how the quasi-energy bands are occupied, provided the laser is of sufficiently high frequency and weak amplitude Dehghani et al. (2015); Dehghani and Mitra (2015, 2016). Thus for this paper, we will model the TR-broken system as the Haldane model, and study its pairing susceptiblity following an interaction quench.

Studying the onset of superconducting instabilities in hexagonal lattices is also of experimental relevance due to the realization of such lattices using cold-atomic gases with tunable interactions Uehlinger et al. (2013); Jotzu et al. (2014). Moreover, recent experiments in bilayer graphene have shown that a strong non-linear coupling to phonons can be achieved Pomarico et al. (2017) further increasing the possibility of inducing superconductivity.

Superconductivity in graphene has a long history Castro Neto et al. (2009); Roy and Juričić (2014). It has been predicted that graphene, with sufficient doping so that one is near the van-Hove singularities, can realize superconducting order Nandkishore et al. (2012); Black-Schaffer and Honerkamp (2014). In particular this order can be either of the $d_{x^{2}-y^{2}}$ or $d_{xy}$ kind. However unlike $d_{x^{2}-y^{2}},d_{xy}$ orders that have nodes, the spontaneously TR symmetry breaking chiral $d_{x^{2}-y^{2}}\pm id_{xy}$ order opens up a gap everywhere on the Fermi surface. Thus, from purely energetic arguments, the preferred state in thermal equilibrium is one of the two chiral $d$ wave states Nandkishore et al. (2012); Black-Schaffer and Honerkamp (2014). The $s$-wave order, while an allowed symmetry of the lattice, is always associated with a lower critical temperature $T_{c}$ than the chiral $d$-wave order, so that the latter is the preferred phase in thermal equilibrium. In the rest of the paper, for notational convenience we will denote $d_{x^{2}-y^{2}}\pm id_{xy}$ by $d\pm id$.

The above equilibrium results for graphene will form a backdrop for comparing our results for the superconducting order following an interaction quench. Firstly we will show that for the Haldane model realized from graphene subjected to a high frequency laser, the broken TR symmetry naturally lifts the degeneracy between the two chiral modes, $d+id$ and $d-id$. Secondly, due to band flattening caused by the laser, one need not dope the Haldane model to the same degree as graphene to achieve the same $T_{c}$. Thirdly, and the main result of the paper, the interaction quench amplitude can influence which superconducting order parameter is preferred by the system. Thus although for weak interaction quenches, the symmetry of the order parameter is the same as in thermal equilibrium, for larger quench amplitudes, a smooth crossover to a phase with a different superconducting order is found.

Here we discuss some subtleties both in equilibrium and out of equilibrium related to studying superconductivity in two dimensional systems. In equilibrium, Mermin-Wagner theorem does not allow for true long range order but only quasi-long range order. Nevertheless, the quasi-long range order is associated with correlation lengths that are fairly long so that from the point of view of local measurements such as the spectral density, the system appears superconducting. For this reason, the study of superconducting instabilities in two dimensional systems has been a very active area of research. In this paper, in addition to two spatial dimensions, we are exploring superconductivity in a transient regime where the system has not fully thermalized. Time-resolved measurements Graf et al. (2011); Smallwood et al. (2012, 2014) are capable of probing the symmetry of the superconducting fluctuations, and we expect the results of the paper to be relevant for such experiments.

The paper is organized as follows. In Section II we describe the model, and outline the derivation of the pairing susceptibility and the equation of motion of the superconducting order parameter. In Section III we discuss the order parameter symmetries for graphene, while in Section IV we do the same for the Haldane model. We present our results for the time evolution of superconducting fluctuation after an interaction quench in both graphene and the Haldane model in Section V, and discuss the phase diagram in terms of the interaction quench amplitude and temperature of the initial state, for a particular choice of doping. Finally in Section VI we present our conclusions.

II Model

The system we will consider is a hexagonal lattice with nearest neighbor hopping, such as graphene, with and without a circularly polarized laser. In addition we will consider interactions as our goal is to explore superconductivity. Writing the full Hamiltonian as $H=H_{0}(t)+V$, where $H_{0}$ denotes the kinetic energy and $V$ the interactions, the kinetic part is

[TABLE]

Above $k$ is the crystal momentum, $\sigma$ denotes the electron spin, $t_{h}$ is the bare tunneling amplitude between the neighboring $A,B$ sites, $a_{k\sigma},b_{k\sigma}$ are electron annihilation operators for the $A,B$ sub-lattices respectively, $A(t)=A_{0}\left[\cos(\Omega t)\hat{x}-\sin(\Omega t)\hat{y}\right]$ denotes the circularly polarized laser of amplitude $A_{0}$ and frequency $\Omega$. $\vec{\delta}_{i}$ are the three nearest-neighbor vectors of the hexagonal lattice, which in terms of the lattice spacing $a$ are,

[TABLE]

Above we have implicitly performed the replacement $b_{k\sigma}\rightarrow b_{k\sigma}e^{i\vec{k}\cdot\vec{\delta}_{3}}$ in order to restore the periodicity of the unit-cell. Thus $\vec{a}_{i}$ are the translation vectors of the hexagonal lattice,

[TABLE]

The interactions correspond to density-density interactions on sites $i$,$j$, $Un_{i\uparrow}n_{j\downarrow}$. A convenient way to probe the consequences of large $U$ is to map the system to the $t$-$J$ model Baskaran (2002); Black-Schaffer and Doniach (2007)

[TABLE]

with $J=2t_{h}^{2}/U$. $\vec{S}_{i}$ and $n_{i}$ are respectively the spin and number operators on site $i$, and the interactions are only between nearest-neighbor sites. These operators act on a restricted Hilbert space that excludes doubly occupied sites.

Periodic drive and interactions will eventually lead to heating D’Alessio and Rigol (2014), however we are interested in very high frequency $\Omega>6t_{h}$ non-resonant driving, where the time scales for heating processes are exponentially long D’Alessio and Polkovnikov (2013); Abanin et al. (2015); Mori et al. (2016); Bukov et al. (2016). In addition, it has been argued that an effective Hamiltonian obtained from a high-frequency expansion may be used to capture the dynamics Kuwahara et al. (2016); Abanin et al. (2017). With this in mind, in the next sub-section, we derive an effective time-independent Hamiltonian by performing a high-frequency expansion. We retain the leading term that captures the effective TR symmetry breaking due to the circularly polarized laser. This leading term opens up a TR symmetry breaking gap at the Dirac points of graphene, mapping it to an interacting Haldane model.

II.1 Mapping to the interacting Haldane model

Let us Fourier transform the two off-diagonal elements of $H_{0}$ in Eq. (1),

[TABLE]

Above $\alpha_{1}=-\alpha_{2}=\frac{\pi}{3},\alpha_{3}=\pi$ and $J_{m}$ denotes the Bessel function. The amplitude of the laser will be given in terms of the dimensionless quantity $A_{0}a$, while the laser frequency and electron temperature will be given in units of the bare tunneling amplitude $t_{h}$.

The high-frequency expansion for a general Hamiltonian $\hat{H}(t)$ takes the form Eckardt and Anisimovas (2015),

[TABLE]

Keeping only the first two terms in the above series, the kinetic energy term becomes,

[TABLE]

Using the expressions for $h_{ab}^{m},h_{ba}^{m}$ derived above, we obtain,

[TABLE]

$H_{\rm eff}$ is the Haldane model for maximal flux of $\phi=\pi/2$ threading the plaquette Haldane (1988). The above high frequency expansion shows that in general a periodic drive generates longer ranged matrix elements such as the next to nearest-neighbor tunneling terms appearing on the diagonal of $H_{\rm eff}$. In addition, these tunneling matrix elements carry non-trivial phases that depend on the polarization of the drive. The generation of the term proportional to $\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}$ in $H_{\rm eff}$ indicates the broken TR symmetry due to the application of the circularly polarized laser. The high-frequency expansion is kept only up to the leading non-zero value as it is sufficient to lift the degeneracy at the Dirac point and open up a gap that corresponds to broken TR symmetry. Keeping higher order terms in the frequency expansion will not generate or break any further symmetries but will only renormalize the parameters of the Floquet Hamiltonian.

For the interaction it is sufficient to keep the $m=0$ term in the high-frequency limit. The interaction $V$ is given by Eq. (4), which we can rewrite as, Baskaran (2002); Black-Schaffer and Doniach (2007)

[TABLE]

where,

[TABLE]

is the nearest-neighbor spin-singlet creation operator. Recall that the elimination of on-site double occupancy in the $t-J$ model excludes on site singlet formation in favor of nearest-neighbor singlet formation. As discussed further, this nearest-neighbor spin singlet formation also gives more structure to the singlet order-parameter, allowing for $s$,$d+id$,$d-id$ symmetries depending on how the order-parameter varies between the three nearest-neighbor sites.

Writing

[TABLE]

it is convenient to parameterize the non-interacting part as,

[TABLE]

Using the pseudo-spin label $\tau=\pm$ to denote the eigenmodes of $H_{\rm eff}$ with energies $\pm|d_{k}|$, we use the standard practice of parameterizing the modes by the two angles $\phi_{k},\theta_{k}$ as follows,

[TABLE]

II.2 Quench dynamics

We now formulate the question as follows. Let us suppose that the interactions are initially zero so that the system is in the normal phase. An attractive interaction of magnitude $J$ is switched on at $t=0$. As the system evolves in time under the influence of this interaction quench, we would like to study the tendency of the system to become superconducting. In particular, we plan to study the time-evolution of superconducting order perturbatively in the interaction.

We consider a spin singlet superconducting order,

[TABLE]

where $\alpha=1,2,3$ denote the three nearest neighbor bonds to the site $i$. For a spatially uniform order parameter, we need not keep the site label $i$. The explicit dependence in $\alpha$ denotes how the order-parameter can have different phases depending on the orientation of the nearest-neighbor bonds. This piece of information is key in order to differentiate between $s$-wave ($\Delta_{\alpha}$ independent of $\alpha$) and a TR broken superconductor (where $\Delta_{\alpha}$ picks up different phases for clockwise and anti-clockwise rotations around a lattice site).

After a mean-field decoupling of the interaction, we obtain,

[TABLE]

where $\Delta_{\alpha}(t)$ is a self-consistently determined pairing field. The mean-field decoupling neglects fluctuation terms ($\Delta_{i}\Delta_{j}$), their role will be discussed later in the section.

While the interaction in Eq. (17) is written in the sub-lattice basis, it can also be written in the energy band basis where the kinetic energy is diagonal. In that basis there are two different pairing interaction terms, one where the electrons from the same band (intra-band) are paired, and one where electrons from different bands (inter-band) are paired. It is important to keep track of both the inter- and intra- band pairing amplitudes.

The equation of motion of the order parameter is,

[TABLE]

where the operators obey Heisenberg time-evolution for an effectively time-dependent Hamiltonian, $\hat{O}(t)=\tilde{T}e^{i\int_{0}^{t}dt^{\prime}H(t^{\prime})}\hat{O}(0)Te^{-i\int_{0}^{t}dt^{\prime}H(t^{\prime})}$. The time-dependence is due to the fluctuating pairing field $\Delta_{\alpha}(t)$ in the mean-field Hamiltonian Eq. (17).

Using Eq. (17), and to leading order in perturbation theory in the pairing amplitude $\Delta_{\alpha}$, one need only retain the commutator of the time-dependent pairing term with the observable. Thus one obtains,

[TABLE]

After performing the averages with respect to the initial state denoted by $\langle\rangle$, we obtain the equation of motion,

[TABLE]

Above $\Pi^{R}$ is the response function or pairing susceptibility of the free electron gas.

Since we assume a spatially uniform order-parameter, only momentum $q=0$ component of $\Pi^{R}$ appears above. For a general momentum $q$, $\Pi^{R}$ has the form,

[TABLE]

where $\rho_{k,\tau}$ is the occupation probability of the $k,\tau$ level in the absence of interactions.

We now briefly discuss the relation of Eq. (19) to what is usually done in equilibrium to determine the superconducting phase. In equilibrium, $t=\infty$, and going into frequency space, one solves $\Delta_{\alpha}=J\sum_{\beta}\Pi^{R}_{\alpha\beta}(q=0,\omega=0)\Delta_{\beta}$. Secondly, the response function $\Pi^{R}$ is evaluated for the full interacting electron gas in the normal phase. Thus the only assumption is that the pairing amplitude $\Delta_{\alpha}$ is small. Beyond this, an additional assumption can be made, and that is of weak $J/t_{h}$. This is the BCS approximation, where now $\Pi^{R}$ is evaluated for the free Fermi gas.

In our problem, namely a quench, since the average is always with respect to the initial state which is that of a free Fermi gas, the linear-response assumption is tied to perturbation theory in $J$ as well. Strictly speaking this assumption will break down when quartic terms ($\Delta_{i}\Delta_{j}$) or inelastic scattering between particles become important. The time scale for this from Fermi’s Golden rule, in units of the hopping $t_{h}$, is $t_{\rm in}\sim 1/(J^{2}T)$, where $T$ is the temperature of the electron gas, and $J,T$ are in units of the hopping $t_{h}$. At times longer than this time, thermalizing processes will become active.

Since finite temperatures are not detrimental to superconductivity in spatial dimensions $d>1$, including $d=2$ provided $T<T_{\rm BKT}$, where $T_{\rm BKT}$ is the Brezenskii Kosterlitz Thouless temperature, our conclusions are not completely invalidated in the inelastic scattering dominated regime. We discuss this issue in more detail in Section V.

Mean-field Yuzbashyan et al. (2015); Peronaci et al. (2015) and linear instability analyses Knap et al. (2016), despite their apparent simplicity, are relevant to experiments that can probe the short time regime where quasi-particles have formed even though the bulk system has not developed a traditional Meissner effect. In particular time-resolved ARPES can pick out the dispersion of the quasi-particles, and hence the symmetries of the order-parameter Graf et al. (2011); Smallwood et al. (2012, 2014).

In what follows we will assume that initially the free fermions were in equilibrium at temperature $T$ and chemical potential $\mu$. Thus for a free fermion dispersion $\epsilon_{k}$ which is also inversion symmetric, $\epsilon_{k,\tau}=\epsilon_{-k,\tau}$, we obtain,

[TABLE]

A given temperature $T$ and chemical potential $\mu$ denotes a doping $\delta$ away from half-filling,

[TABLE]

We will present results for a non-zero doping $\delta=0.1$. For zero doping, since the chemical potential is in the gap in equilibrium, the density of states at the Fermi energy vanishes and superconductivity does not occur in the BCS limit.

In equilibrium, in order to determine the critical temperature for some coupling $J$, one solves the gap equation, $J\Pi^{R}(q=0,\omega=0)=1$, where,

[TABLE]

In the next sections, we will directly solve Eq. (19) for some small, initial randomly chosen $\Delta_{\alpha}$. Since $\alpha$ takes three values, we denote the superconducting order as a vector $\vec{\Delta}$. Growing and decaying solutions in time will indicate whether the system is susceptible to pairing. In addition, how the three components $\Delta_{\alpha=1,2,3}$ of the order parameter grow will indicate the preferred symmetry of the superconducting order parameter.

It might seem that the only place the chemical potential appears is in the initial distribution function. However we will measure all energies with respect to this chemical potential. This choice is equivalent to multiplying the superconducting order-parameter by the phase $e^{-2i\mu t}$. Thus in Eq. (20) we will denote the energy of the upper band $\epsilon_{k+}=\epsilon_{k}+\mu$, and that of the lower band as $\epsilon_{k-}=-\epsilon_{k}+\mu$, where we have also used the fact that particle-hole symmetry in the problem makes the band dispersions symmetric about zero energy.

Although Eq. (19) was derived perturbatively in $J$, yet the solution of the integral equation is non-perturbative in $J$. To see this note that if initially at $t=0$ we had a small seed order parameter $\delta(t^{\prime})\Delta_{0}$, then the first order correction from Eq. (19) is $\Delta_{1}(t)=J\Pi^{R}(t,0)\Delta_{0}$. The second order correction is obtained from substituting this back in Eq. (19) to obtain $\Delta_{2}(t)=J^{2}\int dt_{1}\Pi^{R}(t,t_{1})\Pi^{R}(t_{1},0)\Delta_{0}$, and so on. Thus the solution for Eq. (19) can be recast as an integral equation $\Delta(t)=D^{-1}\Delta_{0}$, essentially the Dyson equation in real time, where $(1-J\Pi)=D$. This leads to a solution, which after a short transient $t<O(1)$ stabilizes to an exponentially growing or decaying solution, $\Delta\sim e^{(J-J_{c})t}$ where $J_{c}$ is a critical coupling. The appearance of $J$ in the argument of the exponential shows the non-perturbative role of $J$ in determining the superconducting phase.

From complex analysis, the relation $\Delta(t)=D^{-1}\Delta_{0}$ implies that locations of zeros of $D$ in the complex frequency plane will determine the growth or decay rate of the order parameter. Thus for later discussions we will find it convenient to interpret the results of the time evolution in terms of the location of the zeros of $D$ in the complex frequency (denoted by $z$) plane. The explicit form for $D(q,z)$ is,

[TABLE]

In section III we discuss the symmetries of the superconducting order for graphene while in section IV we discuss the same for the Haldane model. The results of the time-evolution are presented in Section V.

III Pairing symmetries of graphene

Since graphene has only nearest neighbor hopping, $t_{2}=0,t_{1}=t_{h}$, the modes are,

[TABLE]

$\phi_{k}$ is defined in Eq. (13).

Since under $k\rightarrow-k$, $H_{k}\rightarrow\tau_{x}H_{k}\tau_{x}$, then,

[TABLE]

This gives the following expression for the pair susceptibility,

[TABLE]

The dc component of the response function $\Pi^{R}(\omega=0)$ agrees with Ref. Black-Schaffer and Honerkamp, 2014, where the first line above corresponds to intra-band pairing, and the second line to inter-band pairing.

Let us now discuss the symmetries of the $\Pi^{R}$ matrix. For every $\vec{k}$, there are two others oriented by $2\pi/3,4\pi/3$ from it. We label the triad as,

[TABLE]

We label the nearest-neighbor vectors $\vec{\delta}_{\alpha}$ similarly,

[TABLE]

Thus,

[TABLE]

From above we see $\vec{k}_{m}\cdot\vec{\delta}_{j}=\vec{k}_{m\pm 1}\cdot\vec{\delta}_{j\mp 1}$. In addition if $j$ or $m$ change by $\pm 3$ or its multiples, the function comes back to itself. This implies the following relations,

[TABLE]

The above equalities simply reflect the $C_{3}$ symmetry of the hexagonal lattice. While the first two equalities above directly influence the energy eigenvalues and eigenvectors, the diagonal component of the pairing susceptibility $\Pi_{\alpha\alpha}$ is of the last form. Note that the vectors $\vec{\delta},\vec{a}$ are related by a constant shift (see Eq. (3)), thus the arguments for the symmetries of the $\Pi$ matrix hold irrespective of whether the $\Pi$ matrix is written in the $\vec{\delta}$ or $\vec{a}$ basis.

Thus the $C_{3}$ symmetry implies $\Pi_{11}=\Pi_{22}=\Pi_{33}=A$. Other components of $\Pi$ may be written as,

[TABLE]

The above implies $\Pi_{\alpha\beta}=\Pi_{\beta\alpha}$. This together with the fact that $L_{ji}=L_{j\pm 1,i\pm 1}$ gives $\Pi_{12}=\Pi_{23}=\Pi_{31}=B$ and $\Pi_{12}=\Pi_{13}=\Pi_{23}=B$.

In fact the above symmetries hold at any instant of time $t$, so that $\Pi^{R}$ has the general form,

[TABLE]

From the structure of Eq. (30) it follows that the eigenvalues and eigenvectors of the $\Pi^{R}$ matrix at any instant of time are

[TABLE]

While the eigenvalues depend explicitly on time, the eigenvectors do not. Thus all throughout the time-evolution, there are three mutually orthogonal directions that remain the same. These directions correspond to a non-degenerate $s$-wave solution and a doubly degenerate $d$-wave solution. Any random initial condition for the superconducting order parameter evolves independently along these three directions. We track the time-evolution, and the fastest growing order parameter determines the nature of the superconductor at steady state. Due to the degeneracy of the two $d$-wave modes, the growth rate of the two $d$-wave orders, or any linear combination of the two $d$-wave orders will be identical.

Here we should mention that in equilibrium, the $T_{c}$ of the superconducting phase and its symmetry is determined from $1-J\Pi^{R}(\omega=0)=0$. This treatment Black-Schaffer and Honerkamp (2014) gives the largest eigenvalue, and therefore the dominant instability to correspond to $d$-wave. Since this is doubly degenerate, the order parameter symmetry is not uniquely determined. Instead, energetic considerations indicate that the TR breaking combination d_{x^{2}-y^{2}}+id_{xy}=\frac{1}{\sqrt{3}}\biggl{[}1,e^{2\pi i/3},e^{4\pi i/3}\biggr{]} (or its complex conjugate $d_{x^{2}-y^{2}}-id_{xy}$) will be favored. This is because the chiral order parameter has no nodes.

For the case of the quench, such energy minimization considerations no longer hold. Instead if the time-evolution shows the $d$-wave to be the fastest growing order parameter, the precise order parameter could be any linear combination in the degenerate sub-space. If the fastest growing order parameter is the $s$-wave, then the order parameter symmetry is uniquely determined. For the Haldane model, as we discuss below, the breaking of TR lifts the degeneracy in the $d$-wave sub-space.

IV Pairing symmetries of the Haldane model

We now discuss the pairing susceptibility of the Haldane model. We again make an assumption that the fermions are in thermal equilibrium at temperature $T$ and chemical potential $\mu$ before the interaction quench. This gives,

[TABLE]

Like graphene, the Haldane model also has $C_{3}$ symmetry associated with invariance under rotations by $2\pi/3$. Thus following the arguments given for graphene in the previous section, we conclude, $\Pi_{11}=\Pi_{22}=\Pi_{33}=A$, $\Pi_{12}=\Pi_{23}=\Pi_{31}=B$, and $\Pi_{13}=\Pi_{21}=\Pi_{32}=C$. Thus we may write,

[TABLE]

Unlike graphene however, $\Pi_{\alpha\beta}\neq\Pi_{\beta\alpha}$, when $\alpha\neq\beta$. In the dc limit $\Pi^{R}(\omega=0)$ is Hermitian forcing $C(\omega=0)=B^{*}(\omega=0)$.

The eigenvalues and eigenvectors of the $\Pi(t)$ matrix are

[TABLE]

Graphene corresponds to the case where $\Pi_{\alpha\beta}=\Pi_{\beta\alpha}$ making $B=C$ and $\lambda_{d+id}=\lambda_{d-id}=A-B$. Thus we find that there are no degenerate eigenvalues for the Haldane model as broken TR prefers one of the chiral $d$-wave modes over the other.

Just as for graphene, here too the eigenvalues of $\Pi^{R}$ are time-dependent, while the mutually orthogonal directions stay static. We again study the time-evolution of an initial random but small vector $\vec{\Delta}$, project the evolution along the three mutually orthogonal directions corresponding to the eigenvectors of $\Pi^{R}$, and determine the phase from the fastest growing mode.

Before we present the results for the quench, we point out that in equilibrium, the Haldane model obtained from applying a high frequency circularly polarized laser to graphene, is more susceptible to pairing than graphene. This is because the laser flattens out the band as can be seen from the renormalization of the hopping amplitudes from $t_{h}\rightarrow t_{h}J_{0}(A_{0}a)$ in Eq. (10). As a consequence, the same $T_{c}$ can be obtained at much lower doping levels in the presence of the laser than in graphene. As an example, the pairing susceptibilities derived above give that graphene for $J=0.6t_{h}$ has a $T_{c}=0.01t_{h}$ at a doping of $\delta=0.11$. In contrast the Haldane model realized from a laser of amplitude $A_{0}a=1.0$ and frequency $\Omega=10t_{h}$ shows the same $T_{c}$ but at a much lower doping of $\delta\sim 0.01$.

V Results

We now present results for the solution of the evolution equation Eq. (19). We start with an initial small (precise number does not matter), random vector $\vec{\Delta}$ at $t=0$, and time-evolve it forward. We project the time-evolution along the three mutually orthogonal eigenvectors of the pairing susceptibility $\Pi^{R}$. For graphene these correspond to the $s$-wave order parameter and a doubly degenerate sub-space with $d$-wave symmetry. For the Haldane model, the degeneracy is lifted into two chiral solutions $d+id$, and $d-id$ respectively.

The typical time-evolution of $\vec{\Delta}$ for the Haldane model is shown in Fig. 1. The slopes are proportional to the growth rate, and for the parameters chosen, $d+id$ grows faster than $d-id$, followed by $s$-wave. The magnitude of the order-parameter at any given time depends on the initial condition, and can be rescaled away, and is therefore not of physical relevance. Only the growth rates convey the main physics. For this case we conclude that the preferred superconducting phase is $d+id$.

The time-evolution for doped graphene is similar to Fig. 1 except that the slopes for $d+id$ and $d-id$ are the same, reflecting the degenerate eigenvalues of $\Pi^{R}$. While Fig. 1 is for parameters where all three orders grow in time, the disordered phase is characterized by all three orders decaying exponentially in time. In addition, for $J$-values smaller than that shown in Fig. 1, we can have a situation, where only one order-parameter (typically $d$-wave for the doping levels discussed) grows, while the others decay in time. For all these scenarios, the fastest growing order-parameter determines the preferred phase of the system.

As the quench amplitude increases, the $d$-wave order parameter growth rate increases. However, after some critical value of $J$, the $s$-wave order parameter begins to grow faster, eventually out growing the $d$-wave order parameter. How the growth rate difference of the $d$-wave and $s$-wave order parameters vary with quench amplitude is shown in Fig. 2. Initially the $d$-wave becomes more unstable as the interaction parameter increases. This result is expected as in general the growth rate of any order parameter is zero at the critical point, and grows (decreases) faster as one moves further into (away) from the ordered phase. However, we find that after a critical $J$, the difference between the growth rates of the two order parameters rapidly approaches zero, with the $s$-wave growing faster than the $d$-wave beyond a critical quench amplitude.

The phase diagram determined by the fastest growing order parameter is shown in Fig. 3 for doped graphene, and Figures 4 and 5 for the Haldane model. The two figures for the Haldane model are obtained by applying a laser of two different strengths but same frequency to graphene. We still take care that the frequency of the laser is larger than the bandwidth of graphene so that the high-frequency expansion is valid.

For all these figures, the transition from the disordered (i.e, normal) to the $d$-wave ordered phase obtained from the time-evolution, coincides with an equilibrium calculation based on identifying the $T$ and $J$ values at which $1=J\Pi^{R}(\omega=0)$. Thus consistent with BCS theory, the boundary between the disordered and $d$-wave phase in Figures 3, 4, 5 are strongly dependent on the density of states at the Fermi energy. For this reason, the onset of superconductivity in Fig. 5 occurs at a smaller value of the interaction as compared to the other figures because the larger amplitude laser in Fig. 5 flattens the bands more, increasing the density of states at the Fermi energy.

For larger values of $J$ (i.e., as the quench amplitude increases), eventually the $s$-wave grows faster for graphene as well as the Haldane model. This result can be understood as follows. For smaller quench amplitudes, the dynamics is primarily from quasi-particle excitations in the vicinity of the Fermi surface. Thus the equilibrium phase predicting $d$-wave order is recovered. In contrast, for larger quench amplitudes, the dynamics is governed by highly excited quasi-particles. These quasi-particles cause an effective dephasing, leading to an averaging over the entire Brillouin zone (BZ). This averaging favors an $s$-wave rather than a $d$-wave because the latter order parameter changes sign in the BZ, so that its magnitude is effectively averaged out to zero by the dephasing.

Since at short times $t\sim O(1)$, quasi-particles everywhere in the BZ participate in the dynamics, the $s$-wave component of the order-parameter grows faster than the $d$-wave (see Fig. 1). If the quench amplitude is large, the initial impulse on the $s$-wave order-parameter is large enough to overtake the growth of the $d$-wave. For the particular case shown in Fig. 1, this initial impulse is not strong enough, and at long times the $d$-wave grows faster.

Note that the boundary between the $d$-wave and $s$-wave phases in Figures 3, 4 and 5 is rather flat as a function of temperature. This is simply reflecting the fact that since quasi-particles everywhere in the BZ are participating in the dynamics for such large quench amplitudes, these are not sensitive to the details of the initial distribution function. The location of the flat line between the $d$-wave and $s$-wave phases does however depend on the initial state through the doping level.

In all the three Figures 3, 4 and 5, we plot the phase diagram up to the temperatures where the $d$-wave and the $s$-wave growth rates are clearly different. At higher temperatures, all the order parameter growth rates are degenerate, and one has to account for effects beyond mean-field in order to lift their degeneracy, an analysis that is beyond the scope of the paper.

We make our observation regarding the favored superconducting order more formal by noting that the growth or decay rate of an initial order parameter fluctuation is determined by the location of the zeros of $D(q,z)=1-J\Pi^{R}(q,z)$ in the complex $z$ plane (see Eq. (24) and discussion above it). To locate the zero $z_{0}$, the real and imaginary parts of $z$ should be tuned together to satisfy the two conditions $1/J={\rm Re}[\Pi^{R}(z_{0})]=0,{\rm Im}[\Pi^{R}(z_{0})]=0$. The real part of the zero ${\rm Re}\left[z_{0}\right]$ determines the exponential growth or decay rate, while the imaginary part ${\rm Im}\left[z_{0}\right]$ denotes an overall oscillation.

In equilibrium, the critical temperature $T_{c}$ for a given $J$ is one where the zero approaches $z=0$ linearly in the deviation from the critical temperature or critical coupling. For an interaction quench deeper and deeper into the ordered phase, the pole shifts further out into the complex plane, indicating that the order parameter grows faster. This behavior can be seen in Fig. 6 where the growth rate is determined by the value of $z$ for which $1/J$ intersects with the real part of $\Pi^{R}(z)$. In addition, Fig. 6 clearly shows that as $J$ increases, first the $d$-wave grows faster. But for larger values of $J$, eventually the $s$-wave order grows fastest.

In Fig. 6, for simplicity we allow $z$ to be purely real, and so we assume that the zeros fall entirely on the imaginary axis, or equivalently the order-parameter has a purely exponentially growing or decaying component in time, and no oscillatory component. We look for solutions where $1/J={\rm Re}[\Pi^{R}(z)]$. Inspite of not doing the analysis in the full complex plane, the growth rate and critical couplings agree very well with the full time evolution. Since the full correct analysis is the real time evolution, this good agreement between the time-evolution and the zeros in the complex plane indicate that ${\rm Im}\left[z_{0}\right]\simeq 0$.

We now discuss the effect of inelastic scattering. This will have two effects. One is to thermalize the electron gas, giving rise to a temperature larger than the initial temperature $T$ of the electron gas. As long as the final temperature is less than $T_{\rm BKT}$, the system will show quasi-long range order at long times. The second effect of inelastic scattering is that it causes the breakdown of the mean-field approximation by causing scattering between superconducting fluctuations. For the parameters of Fig. 1, the Fermi Golden rule estimate for the inelastic scattering time (which is also the time for leaving the prethermal regime) is $t_{\rm in}=1/J^{2}T$ $\sim 30t_{h}^{-1}$. Fig. 1 shows that the physical quantity of interest, the growth rate or slope has stabilized well before this time. In addition the growth rate itself is order of the hopping amplitude, $\nu_{s,d\pm id}\sim t_{h}$, and is therefore much larger than the inelastic scattering rate ($t_{\rm in}^{-1}\sim 0.03t_{h}$). Thus for the parameters discussed in the paper, the transient regime is indeed dominated by the mean field estimate because the scattering rate between fluctuations is small as compared to the rate at which the order-parameter grows.

VI Conclusions

The growth of superconducting order under nonequilibrium conditions is an active area of research. In this paper we give an explicit example of how two competing superconducting states behave after an interaction quench. We show that tuning the interaction quench amplitude can favor one superconducting state over the other.

This result is relevant to experiments in cold-atom gases where interaction quenches into the superconducting phase can be performed. Recent experiments in pump-probe spectroscopy make these results also relevant to the solid-state. In fact a laser can modify the lattice parameters, and through it the effective interactions in the $t$-$J$ model used by us to study superconductivity. Thus even in the solid-state, a laser quench can effectively give rise to an interaction quench.

While our study was for a hexagonal lattice with and without TR symmetry, our result showing how the dominant phase can be tuned by nonequilibrium conditions, such as an external laser and interaction quenches, is rather general and applicable to a wide range of lattice models where different ordered phases compete.

Our study treated the effect of the TR symmetry breaking laser within a high frequency expansion. The dynamics in the presence of a resonant low-frequency laser is expected to be even richer, and is left for future studies. Another interesting and challenging direction of study is to go beyond mean-field by allowing for interactions between the superconducting fluctuations. Since a laser, by its spatial orientation at any given time, breaks the underlying lattice symmetries, the interactions in principle can couple different superconducting channels. This could affect the outcome of the symmetry of the dominant superconducting order-parameter at steady-state.

Acknowledgements: The authors thank Yonah Lemonik for many useful discussions and Daniel Yates for numerical help. This work was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0010821.