Universal Prethermal Dynamics in Gross-Neveu-Yukawa Criticality

TL;DR

This paper investigates the universal prethermal dynamics near a quantum critical point in the Gross-Neveu-Yukawa model, revealing two fixed points and unique fermionic signatures, with potential experimental implications.

Contribution

It identifies a new non-thermal fixed point in fermionic systems and characterizes the crossover between different prethermal scaling regimes.

Findings

Discovery of a non-thermal fixed point unique to fermionic fields.

Negative initial slip exponent for multiple fermion flavors.

Temporal crossover between two universal prethermal regimes.

Abstract

We study the prethermal dynamics of the Gross-Neveu-Yukawa quantum field theory, suddenly quenched in the vicinity of a quantum critical point. We find that the universal prethermal dynamics is controlled by two fixed points depending on the size of the quench. Besides the usual equilibrium chiral Ising fixed point for a shallow quench, a dynamical chiral Ising fixed point is identified for a deep quench. Intriguingly, the latter is a non-thermal fixed point without any equilibrium counterpart due to the participation of gapless fermionic fields. We also find that in the scaling regime controlled by the equilibrium fixed point, the initial slip exponent is rendered negative if there are enough flavors of fermions, thus providing a unique signature of fermionic prethermal dynamics. We then explore the temporal crossover between the universal scaling regimes governed by the two universality classes. Possible experimental realizations are also discussed.

Figures

Tables

Table 1. Table 1: Critical initial slip exponents of the ECIFP and the DCIFP for the shallow quench and deep quench, respectively. Note the different definitions of ϵ italic-ϵ \epsilon in these two cases. η 𝜂 \eta and ν 𝜈 \nu are also exhibited for comparison.

Prethermal universality class	$\theta$	$\eta \left(\equiv 2{\eta }_b\right)$	${\eta }_f$	${\nu }^{-1}$
Equilibrium chiral Ising ($ϵ=3-d$)	$\frac{3-37N+\sqrt{9+N\left(66+N\right)}}{24\left(3+N\right)}ϵ$	$\frac{N}{3+N}ϵ$	$\frac{1}{4\left(3+N\right)}ϵ$	$2-\frac{3+5N+\sqrt{9+N\left(66+N\right)}}{6+N}ϵ$
Dynamical chiral Ising ($ϵ=4-d$)	$\frac{ϵ}{12}$	$𝒪\left({ϵ}^2\right)$	$\frac{ϵ}{12}$	$2-\frac{ϵ}{3}$

Equations

$$\displaystyle H(r,g,u)\!=\!\int_{x}\!\Big{[}\frac{1}{2}\pi^{2}+\frac{1}{2}(\nabla\phi)^{2}+\frac{r}{2}\phi^{2}+\frac{u}{4!}\phi^{4}$$

$$\displaystyle\quad\quad\quad\quad\quad+\sum_{\alpha}i\bar{\psi}_{\alpha}\gamma\cdot\nabla\psi_{\alpha}+g\phi\sum_{\alpha}\bar{\psi}_{\alpha}\psi_{\alpha}\Big{]},$$

$$i{D}_K(k,t,t^{\prime })\approx \frac{{\omega }_{0k}}{2{\omega }_k^2}[\cos {\omega }_k(t-t^{\prime })-\cos {\omega }_k(t+t^{\prime })],$$

$$S_{\text{int}}$$

$$\delta S=⟨S_{\text{int}}{⟩}_{>}+\frac{i}{2}⟨S_{\text{int}}^2{⟩}_{>}-\frac{1}{6}⟨S_{\text{int}}^3{⟩}_{>}-\frac{i}{24}⟨S_{\text{int}}^4{⟩}_{>},$$

$$\frac{d{\bar{g}}_{c,q}^2}{dl}$$

$$\frac{d{\bar{u}}_{c,q}}{dl}$$

$$\frac{d{\bar{\Omega }}_0}{dl}$$

$$\displaystyle S_{s}=\frac{1}{2}\int_{k}\big{(}\omega_{0k}|\phi_{0q}({\bf k})|^{2}-\frac{|\dot{\phi}_{0q}({\bf k})|^{2}}{\omega_{0k}}\big{)},$$

$${\eta }_0=-\frac{{\bar{u}}_c}{32}+\frac{N}{4}{\bar{g}}_c{\bar{g}}_q.$$

$$D_R(\mathbf{k},t,t^{\prime })=\frac{t^{1+\theta }}{t^{\prime \eta +\theta }}\mathcal{F}(\mathbf{k}t),$$

$$\displaystyle H_{M_{0}}=\frac{1}{2}\int_{x}\Big{[}\pi^{2}+(\nabla\phi)^{2}+\Omega_{0}^{2}(\phi-M_{0})^{2}\Big{]}.$$

$$S_{s,M_0}=S_s+{\int }_x\frac{i{M}_0}{2\sqrt{2}}{\dot{\varphi }}_{0q},$$

$$M(t)$$

$$M(t)=M_0{t}^{\theta }\mathcal{M}(t^{{\mathcal{D}}_{{\Omega }_0}+\frac{\eta }{2}+\theta }{M}_0),$$

$$\frac{d\lambda }{dl}=ϵ\lambda -\frac{3}{8}{\lambda }^2,$$

$$S_s$$

$$S_b$$

$$L_E[\varphi ,{\psi }^{†},\psi ]$$

$$L_{E,b}[\varphi ]$$

$$L_{E,f}[{\psi }^{†},\psi ]$$

$$L[\varphi ,{\psi }^{†},\psi ]$$

$${\varphi }_{c/q}$$

$${\psi }_{c/q}$$

$${\psi }_{c/q}^{†}$$

$$S_0$$

$$L_{0,f}[{\Psi }^{†},\Psi ]$$

$$L_{0,b}[{\varphi }_{c/q}]$$

$$S_{\text{int}}$$

$$\displaystyle\hat{G}=-i\langle\Psi\Psi^{\dagger}\rangle=\left(\begin{array}[]{cccc}G_{R}&G_{K}\\ 0&G_{A}\end{array}\right),$$

$$\displaystyle\hat{D}=-i\langle\Phi\Phi^{\dagger}\rangle=\left(\begin{array}[]{cccc}D_{K}&D_{R}\\ D_{A}&0\end{array}\right).$$

$$\displaystyle W[\xi_{c/q}]=\ln\int D\mu_{f}\exp\Big{[}-\int d\tau L_{E,f}[\psi^{\dagger}_{i},\psi_{i}]+i\int dt\big{(}L_{0,f}[\Psi^{\dagger},\Psi]+(\xi_{c}^{\dagger}\psi_{c}+\xi_{q}^{\dagger}\psi_{q}+H.c.)\big{)}\Big{]},$$

$$({\partial }_{\tau }+\mathcal{H}){\psi }_i=0.$$

$${\psi }_i(\tau )=[P_{+}(k)e^{-k\tau }+P_{-}(k)e^{k\tau }]{\psi }_0,$$

$${\psi }_c(0)=\frac{1}{\sqrt{2}}(P_{+}(k)e^{-k\beta }+P_{-}(k)e^{k\beta }-1){\psi }_0,$$

$${\psi }_q(0)=\frac{1}{\sqrt{2}}(P_{+}(k)e^{-k\beta }+P_{-}(k)e^{k\beta }+1){\psi }_0,$$

$${\psi }_c^{†}(0)=\frac{1}{\sqrt{2}}{\psi }_0^{†}(P_{+}(k)e^{-k\beta }+P_{-}(k)e^{k\beta }+1),$$

$${\psi }_q^{†}(0)=\frac{1}{\sqrt{2}}{\psi }_0^{†}(P_{+}(k)e^{-k\beta }+P_{-}(k)e^{k\beta }-1).$$

Full text

*Introduction.—*Memory effects are ubiquitous phenomena in nature. The classic example is brain memory in living beings, but memory effects are also widespread in physics. In cosmology, the cosmic microwave background radiation can be regarded as a memory of the Big Bang. In condensed matter physics, memory effects often occur in relaxation dynamics Calabrese2005 ; tauberbook04 . For example, in classical critical systems, the short-time critical dynamics remembers the initial state and affects the critical relaxation process during a macroscopic initial stage janssen1988 . In this stage, the evolution of the system is called the critical initial slip and is characterized by an additional critical exponent, the initial slip exponent janssen1988 . These short-time critical dynamics have been widely exploited in determining the critical point and critical exponents in classical systems Li1995 ; Zheng1998 .

In the context of isolated quantum systems, a vibrant set of purely quantum memory phenomena have been studied. At long times, quantum chaotic systems are expected to effectively lose memory of their initial state, except for conserved quantities like the total energy. This is encoded in the celebrated eigenstate thermalization hypothesis which states that chaotic energy eigenstates look like equilibrium thermal states of the appropriate temperature deutsch1991 ; srednicki1994 ; rigol2008 ; vengalattor2011 ; Dziarmaga2010 . However, some long-lived prethermal states, which bring in additional universal initial state information into the dynamics, have been discovered wetterich2004 ; gring2012 ; langen2013 ; Eigen2018 ; mitra2018 ; Langen2016 ; ueda2018b . Various causes for these effects have been proposed, including proximity to integrability Kollar2011 ; Kastner2013 ; Marcuzzi2013 ; Smith2013 ; Bertini2015 ; Diehl2016 , existence of a dynamical phase transition after a global quench Calabrese2006 ; Calabrese2007 ; Eckstein2009 ; Sciolla2010 ; Fabrizio2010 ; demler2011 ; Tsuji2013a ; Tsuji2013b ; Sciolla2013 ; Heyl2013 ; Sondhi2013 ; Silva2015 , emergence of a non-thermal fixed point Berges2008 ; Gasenzer2011 ; Gasenzer2012 ; Berges2014 ; Berges2015 ; Erne2018 ; Oberthaler2018 ; Gasenzer2019 ; Fujimoto2019 ; Berges2019 , non-local initial entanglement effects ueda2015 ; ueda2018 , and so on. Among these, short-time quantum critical dynamics have recently received considerable attention Polkovnikov2013 ; Yin2014 ; Schmalian2014 ; Schmalian2015 ; mitra2015 ; Gambassi2015 ; mitra2016 ; marino2017 . As with its classical counterpart, a critical initial slip exponent can be defined to characterize the dynamical fingerprint of the initial state. However, unlike the classical case in which the short-time dynamics is completely controlled by a corresponding equilibrium critical point, short-time quantum critical dynamics in an isolated system can also be controlled by a dynamical fixed point, which is similar to a thermal fixed point for a purely bosonic field mitra2015 ; Gambassi2015 ; mitra2016 ; marino2017 .

In the context of equilibrium criticality, gapless Dirac fermions, which appear in various systems including graphene, the surface of topological insulators, and Weyl semimetals, lead to exotic quantum phase transitions sachdevbook . The best known example is the chiral Ising universality class nambu1961a ; nambu1961b ; gross1974 ; scherer2018 , in which the gapless Dirac fermion is coupled to a bosonic field via a Yukawa coupling. However, dynamical phase transitions in these systems have rarely been studied. Given ubiquity of Dirac systems in nature and the exotic universal physics associated to them, studies of the nonequilibrium behavior of these systems are urgently called for. Specifically, some important questions arise: To what extent is the prethermal dynamics affected by fermions? Can an associated dynamical fixed point always be cast as a thermal fixed point?

In this paper, we explore memory effects in Dirac systems. We study the prethermal dynamics of a Dirac system with $N$ flavors of two-component Dirac fermions nambu1961a ; nambu1961b ; gross1974 ; scherer2018 after a sudden quench to the critical point. Using the Keldysh renormalization group (RG) analysis kamenev2011 , we identify two non-trivial fixed points. The usual equilibrium chiral Ising fixed point (ECIFP) controls the near-equilibrium prethermal dynamics induced by a shallow quench, and a dynamical chiral Ising fixed point (DCIFP), which is a non-thermal fixed point, controls the prethermal dynamics after a deep quench. For both cases, we find that prethermal dynamics influenced by the initial condition is characterized by a critical initial slip exponent. Moreover, a negative initial slip exponent, induced by the fermionic degree of freedom, is found for shallow quenches. We also explore the temporal crossover between the two universal scaling regimes. Finally, the eventual long-time thermalization stage and possible experimental realizations are discussed.

*Quench protocol.—*The system consists of $N$ flavors of two-component Dirac fermions coupled to a real boson (which can be understood as the fluctuations of an order parameter) in $d$-spatial dimensions with Hamiltonian,

[TABLE]

where $\int_{x}\equiv\int d^{d}x$, $\phi$ and $\pi$ are the boson and its conjugate momentum, respectively. The fermion notation is as follows: $\psi_{\alpha}$ denotes Dirac fermions with flavor index $\alpha=1,...,N$, $\bar{\psi}\equiv\psi^{\dagger}\gamma^{0}$, $\gamma\cdot\nabla\equiv\sum_{i=1}^{d}\gamma^{i}\partial_{i}$, and $\gamma^{i}$ are the gamma matrices. $r$ is the mass of the boson, which controls the distance to the critical point, while $g$, $u>0$ refer to the Yukawa coupling and four-boson coupling, respectively. For $t<0$, the initial state is the ground state of $H(\Omega_{0}^{2},0,0)$. At $t=0$, the Hamiltonian is suddenly quenched to $H(r,g,u)$, and we are interested in the emergent universal behaviors after this global quench supp .

Without the coupling to fermions, $g=0$, it was shown that the nature of the universal prethermal dynamics depends on ratio of $\Omega_{0}$ and $\Lambda$, the UV energy cutoff: the system exhibits a dimensional crossover from an equilibrium “quantum” Wilson-Fisher fixed point for a shallow quench $\Omega_{0}\ll\Lambda$ to a dynamical “classical” fixed point for a deep quench $\Omega_{0}\gg\Lambda$ marino2017 . In lattice Dirac systems $\Lambda$ is set by the inverse lattice constant (in units where the Dirac velocity is one). It has also been shown that the “classical” fixed point has some of the features of the equilibrium transition with an effective temperature determined by $\Omega_{0}$ mitra2015 ; Gambassi2015 ; mitra2016 ; marino2017 .

When the boson-fermion coupling is turned on, $g>0$, the Dirac fermions can affect the prethermal dynamics. For a quench of the type described above, the initial state information is contained in the free boson Keldysh propagator which reads supp

[TABLE]

where $\omega_{0k}^{2}=k^{2}+\Omega_{0}^{2}$ and $\omega_{k}^{2}=k^{2}+r$, $k^{2}\equiv\sum_{i=1}^{d}k_{i}^{2}$. Note that the second term explicitly breaks time translation invariance due to the quench. If $\Omega_{0}\gg\Lambda$, $iD_{K}\approx\frac{\Omega_{0}}{2\omega_{k}^{2}}[\cos\omega_{k}(t-t^{\prime})-\cos\omega_{k}(t+t^{\prime})]$, one can compare it with equilibrium Keldysh propagator at low temperature, $iD_{K}\approx\frac{2T}{\omega_{k}^{2}}\cos\omega_{k}(t-t^{\prime})$, which shows that $\Omega_{0}$ plays a role of effective temperature, $T_{\text{eff}}=\frac{\Omega_{0}}{4}$.

Because $\Omega_{0}$ sets an effective temperature scale, one may expect a “quantum” to “classical” crossover similar to the bosonic system marino2017 . Moreover, given their inherently quantum nature, fermions are not expected to play a role in a classical phase transition Stephanov ; Hesselmann . Thus, one might expect that the fermions effectively decouple with a vanishing Yukawa coupling for a deep quench $\Omega_{0}\gg\Lambda$, but we will show that this is not the case.

*RG analysis.—*Our analysis is based on a renormalization group (RG) calculation. In the Keldysh formalism kamenev2011 ; supp , the interacting part of the action is given by

[TABLE]

where $\phi_{c/q}$ ($\psi_{c/q}$) refer to the classical and quantum fields kamenev2011 ; supp , respectively. $\Psi=(\psi_{c},\psi_{q})^{T}$, and $\tau$ is Pauli matrix acting on the classical/quantum fields. Summation over flavors of Dirac fields is implicit. Note that there is no symmetry that relates the classical and quantum fields. Thus, while they have the same bare values, $g_{c/q}$, $u_{c/q}$ are in principle independent coupling constants at lower energy scales.

The fast modes within the momentum shell, $k\in[\Lambda e^{-l},\Lambda]$, are integrated out to generate the RG equation janssen1988 ; mitra2015 . Here $l>0$ denotes the flow parameter. The interacting part Eq. (3) generates RG processes that renormalize the slow modes, namely, $S_{\text{eff}}=S_{0}+\delta S$,

[TABLE]

where $\langle A\rangle_{>}\!=\!\int D\phi_{>}D\bar{\Psi}_{>}D\Psi_{>}Ae^{iS_{0}}$ denotes the functional integration over fast modes, and the calculation is done to one-loop order.

In terms of dimensionless coupling constants, $\bar{\Omega}_{0}\equiv\Omega_{0}\Lambda^{-1}$, $\bar{g}_{c/q}^{2}\equiv K_{d}\Lambda^{d-3}(1+\bar{\Omega}_{0}^{2})^{\pm 1/2}g_{c/q}^{2}$ , $\bar{u}_{c/q}\equiv K_{d}\Lambda^{d-3}(1+\bar{\Omega}_{0}^{2})^{\pm 1/2}u_{c/q}$, where $K_{d}\equiv\frac{2\pi^{(d+1)/2}}{(2\pi)^{d}\Gamma(\frac{d+1}{2})}$, the RG equations supp are,

[TABLE]

where $D_{c,q}(\bar{\Omega}_{0})=3-d\pm\frac{\bar{\Omega}_{0}^{2}}{1+\bar{\Omega}_{0}^{2}}$ can be understood as the bare scaling dimensions of the coupling constants. For $\bar{\Omega}_{0}=0+\mathcal{O}(\epsilon)$, the RG Eqs (5), and (6) reduce to the corresponding equilibrium ones, giving rise to the usual equilibrium chiral Ising universality class (ECIFP, see e.g.,  scherer2018 ). Also, from Eq. (7) we know the quench parameter is a relevant variable with an unstable value $\bar{\Omega}_{0}=0$ and a stable one $\bar{\Omega}_{0}\rightarrow\infty$. As a result, there is a crossover, indicated by $D_{c}(0)=3-d$ and $D_{c}(\infty)=4-d$.

Shallow quench.— For a shallow quench, $\bar{\Omega}_{0}\ll 1$, it is expected that the universal prethermal dynamics is controlled by the ECIFP. Accordingly, the usual chiral Ising critical exponents are applicable to describe the prethermal dynamics. However, an additional initial slip exponent $\theta$ must be included to describe the effects induced by the initial condition.

To see this, we inspect the pre-quench action supp ,

[TABLE]

where $\int_{k}\equiv\int\frac{d^{d}k}{(2\pi)^{d}}$, and $\phi_{0q}$ and $\dot{\phi}_{0q}\equiv\partial_{t}\phi_{0q}$ are initial fields defined at $t=0$. Because the sudden quench explicitly breaks time translation invariance, the Feynman diagrams shown in Fig. 1 lead to a correction of scaling of initial fields. The diagrams give a contribution of $\delta S_{s}=(-\frac{\bar{u}_{c}}{16}+\frac{N}{2}\bar{g}_{c}\bar{g}_{q})S_{s}$, shifting the anomalous dimension of the initial fields,

[TABLE]

The sign difference between these terms arises from the minus sign associated with the fermion loop in Fig. 1 (b).

This shift manifests in behavior of the retarded Green’s function with initial fields $iD_{R}({\bf k},t,0)=\langle\phi_{c}({\bf k},t)\phi_{q0}(-{\bf k})\rangle$. According to the scaling form of boson fields, we have $D_{R}({\bf k},t,0)=t^{1+\theta}D_{R}({\bf k}t,1,0)$, where the critical initial slip exponent is $\theta\equiv-(\eta_{b}+\eta_{0})$, and $\eta_{b}=\frac{\eta}{2}$ is the anomalous dimension of the boson field. For the retarded Green’s function $D_{R}({\bf k},t,t^{\prime})$, when $t^{\prime}$ is very close to zero, we can approximate $\phi_{q}(t)\approx\sigma(t)\phi_{0q}$ janssen1988 , with $\sigma(t)=b^{\eta_{b}-\eta_{0}}\sigma(bt)$, then

[TABLE]

where $\mathcal{F}$ is a universal function.

The critical slip also manifests in the short-time scaling behavior of the order parameter. If the system is prepared with a finite order $M_{0}$, the order parameter shows early-time power-law time dependence set by the initial slip exponent. To set a finite order $\phi\sim M_{0}$, we modify the pre-quench Hamiltonian supp :

[TABLE]

Consequently, the pre-quench action is changed to

[TABLE]

and the order parameter in the presence of the initial magnetization is given by

[TABLE]

where $\int_{Y_{n}}=\int_{y_{1}}...\int_{y_{n}}$. Since the $n=0$ term vanishes,

[TABLE]

where $\mathcal{D}_{\Omega_{0}}$ is the canonical scaling dimension of $\phi_{c}$, and $\mathcal{M}$ is a universal function. Notice that Eq. (14) validates after a non-universal initial time scale $\Lambda^{-1}$.

The critical exponents including the initial slip exponent of ECIFP are shown in Table 1. It is interesting to note that, unlike the pure bosonic case, the initial slip exponent can be rendered negative if the number of fermion flavors $N$ is large enough (at one-loop, $N=1$ is already enough to make $\theta<0$), leading to a unique signature massless fermions in prethermal dynamics. The reason is that the large fluctuations of the fermion drive the critical point to be larger than its mean-field value.

*Deep quench.—*Next, for a deep quench, $\bar{\Omega}_{0}\gg 1$, the system is initially prepared far away from critical point and then suddenly quenched to the dynamical critical point. Unlike for a shallow quench, here $\epsilon=4-d$. The critical and quantum coupling constants are significantly different already at the level of their bare couplings: $D_{c}(\infty)=4-d$ while $D_{q}(\infty)=2-d$. In $d=4-\epsilon$ dimensions, the quantum coupling constants are highly irrelevant, i.e., $g_{q}^{*}=u_{q}^{*}=0$. So, the low-energy dynamics are largely controlled by the classical coupling constants. From Eq. (5), the flow equations reduce to

[TABLE]

for $\lambda=\bar{g}_{c}^{2},\bar{u}_{c}$.

Remarkably, Eq. (15) predicts a new fixed point with a $\bar{g}_{c}^{2}=\bar{u}_{c}=\frac{8\epsilon}{3}$. We call this new fixed point the dynamical chiral Ising fixed point (DCIFP) and observe that it is a non-thermal fixed point. For a conventional fixed point describing thermal criticality, the Yukawa coupling is irrelevant and the fermionic degrees of freedom decouple at low energies because their Matsubara frequencies do not include zero. For example, at finite temperatures, the classical thermal phase transition for the Dirac system, whose quantum phase transition at the zero temperature is the chiral Ising universality class, belongs to the conventional Ising universality class Stephanov ; Hesselmann . However, for the DCIFP, $g_{c}\neq 0$, indicating the participation of the fermions. It was argued that in the pure bosonic model the dynamical fixed point bears similarities to the equilibrium thermal one. Here, we find that the dynamical phase transition can be quite different from the thermal one, with the physics of the dynamical fixed point originating in the interplay between the quenched boson fields and the gapless fermion fields.

The critical exponents associated with the DCIFP are listed in Table 1. The ECIFP Yukawa coupling does not vanish, and it leads to a nonvanishing fermionic scaling dimension $\eta_{f}$. This anomalous dimension is observable, for example, in the energy $\varepsilon$ dependence of fermion density of states, $n(\varepsilon)\propto\varepsilon^{1+2\eta_{f}}$, for time-scales less than the thermalization time $t\ll t_{\text{th}}^{\prime}$. The thermalization time is given by supp $t^{\prime}_{\text{th}}\sim\Lambda^{-1}\epsilon^{-2}$. We emphasize that this non-equilibrium universality class does not have a classical counterpart and only exists in non-equilibrium states of matter.

Temporal crossover between two universality classes.—Above we have shown that depending on the quench variable $\bar{\Omega}_{0}$, the universal prethermal dynamics is controlled by the ECIFP or the DCIFP. Since $\bar{\Omega}_{0}$ is a relevant quantity, the prethermal dynamics near the ECIFP resulting from a shallow quench is unstable at low energies, and the dynamics will flow to a scaling regime governed by the DCIFP.

Note that the crossover RG flow can manifest itself as a temporal crossover. Because $t\sim\mu^{-1}$ with $\mu$ being the energy scale in the question, the running RG parameter is effectively set by $l=\log\Lambda t$. In this sense, the RG flow can be directly seen in time-resolved experiments. Figure 2 shows the temporal crossover between the two scaling regimes. The two plateaus of the Yukawa coupling $g_{c/q}$ in Fig. 2 (a) indicate the RG flow near the ECIFP and the DCIFP, respectively. Since the bare values of $g_{c/q}$ are the same, the flows are identical in the scaling regime controlled by the ECIFP. However, at longer times, a non-zero $\bar{\Omega}_{0}$ drives a difference between $g_{c/q}$ through the bare scaling dimensions $D_{c/q}$. The quantum Yukawa coupling $g_{q}$ flows to zero while the classical $g_{c}$ flows to another finite value corresponding to the DCIFP. Intriguingly, Fig. 2 (b) shows the temporal crossover of the initial slip exponent $\theta$, where $\theta<0$ for large enough $N$ in the regime governed by the ECIFP and $\theta>0$ in the regime governed by the DCIFP.

*Discussion.—*The observability of these non-equilibrium universal scaling behaviors relies on a long thermalization time scale. In the present case, the thermalization can be understood as the generation of the dissipative term, $\gamma\phi_{q}\dot{\phi}_{c}$, in the RG flow mitra2015 ; Gambassi2015 ; mitra2016 ; marino2017 ; mitra2011 . Accordingly, the RG equation of $\bar{\gamma}\equiv\gamma\Lambda^{-1}$ effectively sets the thermalization time scale. Initially, the prequench Hamiltonian describes a noninteracting system without dissipation, i.e., $\gamma=0$. The postquench interactions generate inelastic processes and render the dissipative coupling $\gamma$ nonzero. For a shallow quench supp , the thermalization time is t_{\text{th}}=\Lambda^{-1}\Big{(}1+N\frac{1-\epsilon}{16\alpha_{3}\epsilon^{2}}\Big{)}^{\frac{3+N}{3+(1-\epsilon)N}}, where $\alpha_{3}>0$, and $\epsilon=3-d$. One finds that the thermalization time scale is extremely long at large $N$, even when $\epsilon=1$. Indeed, the thermalization time scale is $t_{\text{th}}\propto(1+\frac{3}{16\alpha_{3}})^{N}$ for $d=2$ and $N\gg 1$. This feature is due to the integrability of zero-temperature Gross-Neveu-Yukawa model in the infinite-$N$ limit. On the other hand, for the deep quench, the thermalization time is t_{\text{th}}^{\prime}=\Lambda^{-1}\Big{(}1+\frac{9}{64\alpha_{1}\epsilon^{2}}\Big{)}, where $\alpha_{1}>0$ and $\epsilon=4-d$. For a deep quench, the thermalization time scale is only controlled by $\epsilon$ because the boson cannot efficiently mediate an interaction between different fermion flavors at the DCIFP.

The results obtained above provide several sharp experimental signatures, in particular, the negative initial slip exponent for a shallow quench and its sign reversal in the temporal crossover to a deep quench arise from the fermion fields. Recently, the generation and detection of nonequilibrium dynamics has been demonstrated in various systems. In particular, in cold atom systems, dynamical scaling has been observed in experiments Oberthaler2015 ; Eigen2018 ; Erne2018 ; Oberthaler2018 . Furthermore, the real-time dynamics has been measured in a tunable honeycomb lattice, which hosts gapless Dirac fermions Greif2015 ; Bloch2017 . Accordingly, the prethermal dynamics studied here appear within experimental reach. It would also be interesting if some aspect of this physics survived to strong coupling, where it might serve as a signature of emergent fermions in candidate spin liquid materials.

To summarize, we studied the prethermal dynamics of Dirac systems. Two fixed points were found to govern the quench dynamics, the usual equilibrium chiral Ising fixed point and a new dynamical chiral Ising fixed point. The latter fixed point is non-thermal due to a non-trivial coupling between the boson and fermion degrees of freedom. The initial slip exponent was calculated for both shallow and deep quenches, and the negative initial slip exponent in the shallow case and its sigh change in the temporal crossover to the deep case provide a sharp signature of prethermal dynamics in Dirac systems.

Acknowledgement: We thank H. Y. Xie and H. Yao for helpful discussions. This work is supported in part by the NSFC under grant 11825404 (SKJ and SY), the Simons Foundation via the It From Qubit Collaboration (SKJ and BS), and the Department of Energy award number de-sc0017905 (BS). S.Y. is also supported in part by China Postdoctoral Science Foundation (Grant No. 2017M620035).

I Supplemental Material

I.1 A. Propagators in Keldysh contour

To prepare the initial state away from the transition point in the disordered phase, we use the Keldysh contour as indicated by Fig. S1. The vertical line in Fig. S1 indicates the preparation of the pre-quench state. The action is $e^{-S_{s}+iS_{b}}$, where

[TABLE]

are the “surface” and “bulk” action. The vertical line is a constant time slice while horizontal lines stand for the bulk of spacetime. The Lagrangians are

[TABLE]

where $\int_{x}\equiv\int d^{d}x$, $d$ is the spatial dimension, and $\dot{\phi}\equiv\partial_{t}\phi$. In this paper, we mainly focus on $d=2$ dimensions. $\phi$ is a real boson serving as an order parameter, and $\psi_{\alpha}$ ($\psi_{\alpha}^{\dagger}$), $\alpha=1,...,N$, denotes the annihilation (creation) operator of two-component Dirac fermion for $\alpha$ flavor. The summation over $N$ flavors is assumed. The Dirac Hamiltonian is $\mathcal{H}=-i\sigma\cdot\nabla$, where $\sigma\cdot\nabla\equiv\sum_{i=1,2}\sigma_{i}\partial_{i}$ and $\sigma_{i}$ is the Pauli matrix. $\Omega_{0}^{2}$ and $r$denote the pre- and post-quench mass of the real boson, respectively. The quench from disordered phase towards the transition point obeys $\Omega_{0}^{2}\gg r$. In the bulk Lagrangian, $g$ and $u$ denote the Yukawa coupling and four-boson coupling, respectively.

To get the propagators in the Keldysh contour, we perform a Keldysh rotation, i.e.,

[TABLE]

In terms of classical/quantum fields, the bulk action is $S_{b}=S_{0}+S_{\text{int}}$, where we have separated the bulk action into a non-interacting part $S_{0}$ and an interacting part $S_{\text{int}}$. They are given by

[TABLE]

where $\int_{k}\equiv\int\frac{d^{d}k}{(2\pi)^{d}}$, $k\cdot\sigma\equiv\sum_{i=1,2}k_{i}\sigma^{i}$, $\omega_{k}\equiv\sqrt{k^{2}+r}$, and $\Psi=(\psi_{c},\psi_{q})^{T}$. Since there is no symmetry that relates the quantum and classical fields, we allow their couplings to evolve independently under RG. These independent coupling constants are $g_{c},g_{q}$ and $u_{c},u_{q}$. $\tau$ is a Pauli matrix acting on the classical/quantum space.

Now we are ready to derive the propagators in the presence of the surface action. Since in both surface and bulk noninteracting action, Eqs. (S3, S10), the fermion and boson are decoupled, we can derive their propagators independently. In terms of the classical/quantum fields, the propagators are given by

[TABLE]

In the following, we will first get the fermion propagator and then the boson propagator.

We introduce source field $\xi_{c/q}$ coupling to $\psi_{c/q}$, and integrate out all dynamical fields to get the generating functional as a function of $\xi_{c/q}$. The generating functional is defined by

[TABLE]

where $D\mu_{f}\equiv D\psi^{\dagger}_{i}D\psi_{i}D\psi^{\dagger}_{c}D\psi_{c}D\psi^{\dagger}_{q}D\psi_{q}$ is the functional measure of all dynamical fields. Integrating over $\psi^{\dagger}_{i}$, we obtain the equation of motion of the initial fields

[TABLE]

Since the equation of motion of fermion is a first order differential equation, only one boundary condition is needed. In order to simplify the notation, we define $\psi_{0}\equiv\psi_{i}(0)$. It is straightforward to get the solution of Eq. (S23), namely,

[TABLE]

where $P_{\pm}(k)\equiv\frac{1}{2}(1\pm\hat{k}\cdot\sigma)$, $\hat{k}\equiv(k_{1}/k,k_{2}/k)$, and $k=\sqrt{k_{1}^{2}+k_{2}^{2}}$. Combining the boundary condition $\psi_{-}(0)=-\psi_{i}(0),\psi_{+}(0)=\psi_{i}(\beta)$ and Eqs. (S8, S9), we have

[TABLE]

After integrating the initial fields, we arrive at

[TABLE]

where we manipulate to get

[TABLE]

Integrating over $\psi_{c}^{\dagger}$ and $\psi_{q}$ leads to the equation of motion of $\psi_{c}$ and $\psi_{q}^{\dagger}$ in the presence of source fields, namely,

[TABLE]

Recalling the boundary condition Eqs. (S25, S28), it is straightforward to get solutions of Eqs. (S31, S32),

[TABLE]

where $G_{R/A}$ are retarded/advanced Green function given in the following,

[TABLE]

Finally, integrating over $\psi_{c}$ and $\psi_{q}^{\dagger}$ is amount to substitute the solution into the action,

[TABLE]

Given the boundary condition Eqs. (S26, S28), we can integrate over $\psi_{0}$ to get the final result,

[TABLE]

from which we have the fermionic Keldysh propagator,

[TABLE]

Next, we will derive the bosonic propagator. We introduce a source field $j_{c/q}$ coupling to $\phi_{q/c}$, and integrate out all dynamical fields to get the generating functional,

[TABLE]

To integrate out $\phi_{i}$, we first solve the equation of motion for $\phi_{i}$,

[TABLE]

where $\omega_{0k}^{2}\equiv k^{2}+\Omega_{0}^{2}$ and we have used the boundary conditions $\phi_{i}(0)=\phi_{-}(0)$ and $\phi_{i}(\beta)=\phi_{+}(0)$. Then we plug the solution into the action to arrive at

[TABLE]

where $\phi_{0c}\equiv\frac{1}{\sqrt{2}}[\phi_{+}(0)+\phi_{-}(0)]$ and $\phi_{0q}=\frac{1}{\sqrt{2}}[\phi_{+}(0)-\phi_{-}(0)]$. For later consideration, we take the imaginary time derivative of the solution Eq. (S43), i.e.,

[TABLE]

and using $\dot{\phi}_{+}(0)=i\partial_{\tau}\phi_{i}(\beta)$, $\dot{\phi}_{-}(0)=i\partial_{\tau}\phi_{i}(0)$, we can get

[TABLE]

As the equation of motion is second-order in derivatives, there are two independent fields, so Eq. S44 can be equivalently expressed as

[TABLE]

Integrating over $\phi_{q}$, we obtain the equation of motion in the presence of source field $j_{c}$,

[TABLE]

It is straightforward to get the solution of Eq. (S48), i.e.,

[TABLE]

We plug the solution into the action and integrate out $\phi_{0c}$ and $\phi_{0q}$ to get

[TABLE]

where $K_{\pm}=\frac{1}{2}(\frac{\omega_{k}}{\omega_{0k}}\pm\frac{\omega_{0k}}{\omega_{k}})$, and in the calculation we have used Eq. (S46). From the generating functional, we can easily obtain the boson propagators,

[TABLE]

In the following, we summarize both fermion and boson propagators derived above,

[TABLE]

where $\omega_{k}\equiv\sqrt{k^{2}+r}$, $\omega_{0k}\equiv\sqrt{k^{2}+\Omega_{0}^{2}}$, and $K_{\pm}=\frac{1}{2}(\frac{\omega_{k}}{\omega_{0k}}\pm\frac{\omega_{0k}}{\omega_{k}})$. It is worth noting that the quench information is encoded in the bosonic Keldysh propagator, while all other propagators have time translation symmetry.

I.2 B. Propagator in the zero temperature limit

Since we are interested in the role of fermions in prethermalization dynamics near a critical point, we should focus on the zero-temperature limit. Since there is no symmetry to relate the classical and quantum fields, we allow the bosonic pre-quench masses to evolve independently. Namely, the pre-quench boson action Eq. (S44) can be generalized into

[TABLE]

where $\omega_{0c}\equiv\sqrt{k^{2}+\Omega_{0c}^{2}}$ and $\omega_{0q}\equiv\sqrt{k^{2}+\Omega_{0q}^{2}}$. It will modifies the bosonic Keldysh propagator $D_{K}(k,t,t^{\prime})$ by

[TABLE]

Because $\Omega_{0c}$ is the pre-quench a relevant scale, we take $1/\Omega_{0c}\rightarrow 0$ and define $\Omega_{0}\equiv\Omega_{0q}$. As we mentioned, Under such simplification, $K_{\pm}=\pm\frac{\omega_{0k}}{2\omega_{k}}$. And the propagators at zero-temperature limit are given by

[TABLE]

I.3 C. Renormalization group calculation

The fast modes within the momentum shell, $[\Lambda e^{-l},\Lambda]$, are integrated out to generate the RG equation. Here $\Lambda$ is the momentum cutoff and $l$ denotes the flow parameter which is a real number. To do that, we separate the fields into slow and fast modes, and the interacting part generates Feynman diagram that renormalize the slow modes, namely,

[TABLE]

where $\langle A\rangle_{>}=\int D\phi_{>}D\Psi^{\dagger}_{>}D\Psi_{>}Ae^{iS_{0}}$ denotes the functional integration over fast modes. In the following calculation, we will keep the calculation up to one-loop order.

We are ready to calculate one-loop RG equations. Feynman diagrams shown in Fig. S2 lead to the renormalization, $\delta S_{b}^{(2)}$ and $\delta S_{f}^{(2)}$, of two-point functions of both fermionic and bosonic fields. To illustrate the calculation, we show in the following the calculation of first two Feynman diagrams in Fig. S2:

[TABLE]

where

[TABLE]

where $\omega_{\Lambda}\equiv\sqrt{\Lambda^{2}+r}$, $\omega_{0\Lambda}\equiv\sqrt{\Lambda^{2}+\Omega_{0}^{2}}$, $K_{d}\equiv\frac{A[S^{d-1}]}{(2\pi)^{d}}$, and $A[S^{d-1}]$ is the area of unit sphere $S^{d-1}$. And

[TABLE]

where

[TABLE]

We will only keep time dependent contribution to $\int\phi_{q}\phi_{c}$ terms. Similarly, straightforward generalization of above calculations to other Feynman diagrams lead to the answer:

[TABLE]

and

[TABLE]

Next, Feynman diagrams shown in Fig. S3 lead to the renormalization, $\delta S^{(3)}$, of three-point vertices, namely,

[TABLE]

Finally, Feynman diagrams shown in Fig. S4 lead to the renormalization, $\delta S^{(4)}$, of four-point vertices, namely,

[TABLE]

Collecting the results from above calculations, the corrections from one-loop diagrams to the slow modes are given by,

[TABLE]

Now to get RG equations, we introduce the dimensionless coupling constants,

[TABLE]

According to the renormalized action, the anomalous dimensions are given by

[TABLE]

And the RG equations are given by

[TABLE]

where $D_{c}(\bar{\Omega}_{0})\equiv 3-d+\frac{\bar{\Omega}_{0}^{2}}{1+\bar{\Omega}_{0}^{2}}$, $D_{q}(\bar{\Omega}_{0})\equiv 3-d-\frac{\bar{\Omega}_{0}^{2}}{1+\bar{\Omega}_{0}^{2}}$, and

[TABLE]

Notice that $\bar{\Omega}_{0}$ is relevant with an unstable value $\bar{\Omega}_{0}=0$ and a stable one $\bar{\Omega}_{0}\rightarrow\infty$. To get some physical intuition, let’s first assume the system is prepared close to critical point, namely, $\bar{\Omega}_{0}=\bar{r}=0+O(1/N)$. One expects the RG equations will reproduce the equilibrium ones. Indeed, the RG equations are given by

[TABLE]

The similarity between $g_{c}$, $u_{c}$ and $g_{q}$, $u_{q}$ indicates only two of these coupling constants are needed, i.e.,

[TABLE]

with a stable fixed point

[TABLE]

and critical exponents listed in Table I. However, we should notice that $\bar{\Omega}_{0}$ is a relevant perturbation at this fixed point. From above lessons, it is straightforward to infer that the critical dynamics at short time after a soft quench is controlled by the Gross-Neveu fixed point with a new nonequilibrium exponent—the critical initial slip.

On the other hand, when the system is initially prepared far away from critical point and suddenly quenched to the critical point, i.e., $\bar{\Omega}\rightarrow\infty$ and $\bar{r}=0+O(\epsilon)$, where $\epsilon=4-d$. Now those critical and quantum coupling constants are very distinct from each other because the bare scalings, $D_{c}(\infty)=4-d$ while $D_{q}(\infty)=2-d$. The RG equations are given by

[TABLE]

which possesses a stable fixed point given by

[TABLE]

The corresponding critical exponents are listed in Table I.

I.4 D. Critical initial slip exponent

The surface action is given by

[TABLE]

The time-dependent contribution will leads to renormalization of boundary field $\phi_{0f/0c}$. The time-dependent contribution is given by

[TABLE]

which leads to anomalous dimension of initial fields, i.e., $\eta_{0}=-\frac{\bar{u}_{c}}{32(1+r)^{2}}+\frac{N}{4}\bar{g}_{c}\bar{g}_{q}$. Let’s consider the retarded Green’s function with initial fields $iD_{R}({\bf k},t,0)=\langle\phi_{c}({\bf k},t)\phi_{q0}(-{\bf k})\rangle$,

[TABLE]

Choosing $b=t^{-1}$, the retarded function has the scaling form

[TABLE]

Now we can introduce the cirtical initial slip $\theta\equiv-(\eta_{b}+\eta_{0})$ To investigate the retarded Green’s function $iD_{G}({\bf k},t,t^{\prime})$ with $t^{\prime}$ very close to zero, we can use $\phi_{q}(t)\approx\sigma(t)\phi_{0q}$, where $\sigma(t)=b^{\eta_{b}-\eta_{0}}\sigma(tb)$,

[TABLE]

Setting $b=t^{\prime-1}$, we have

[TABLE]

where $\mathcal{F}$ is a universal function.

A more experimental relevant quantity revealing the critical initial slip is the scaling form of order parameter. We change the boundary condition to set the initial order parameter to be $M_{0}$, namely, the boundary action is changed to

[TABLE]

At zero temperature limit,

[TABLE]

where we have used Eq. (S46) and take the limit $\Omega_{0c}\rightarrow\infty$ in the second line. It is easy to see the second term has a larger or equal scaling compared to the first term, so in order to get a scaling form, we can only keep the second term. The order parameter in the presence of the initial magnetization is given by

[TABLE]

It should be noticed that $n=0$ term vanishes. Then

[TABLE]

where $\mathcal{D}_{\Omega_{0}}$ is the bare scaling dimension of $\phi_{c}$, namely, $\mathcal{D}_{0}=\frac{d-1}{2}$ and $\mathcal{D}_{\infty}=\frac{d-2}{2}$. Above equation leads to the result presented in the main text, $M(t,M_{0})=M_{0}t^{\theta}\mathcal{M}(t^{\mathcal{D}_{\Omega_{0}}+\frac{\eta}{2}+\theta}M_{0})$.

I.5 E. Estimate of thermalization time

Thermalization is due to effectively irreversible processes. In an isolated system, although quantum time evolution is unitary, when we focus on one part of the system, say the low-energy modes in the Hilbert space, the other parts will effectively serve as a thermal bath. Then thermalization can be understood as the generation of a dissipative term, $\gamma\phi_{q}\dot{\phi}_{c}$, under the RG flow. The prequench Hamiltonian describes a noninteracting system without dissipation, so initially $\gamma=0$. When we quench the system, the interactions generate inelastic processes and drive the dissipative term $\gamma$ nonzero. Because the Keldysh boson propagator encodes the information of the quench protocol, the generation of a finite $\gamma$ is closely linked to the Keldysh boson propagator. The lowest order processes linking the dissipative term and Keldysh boson propagator come from the two-loop Feynman diagram shown in Fig. S5, which gives an RG equation for $\bar{\gamma}\equiv\gamma\Lambda^{-1}$,

[TABLE]

where $\alpha_{i}$ are some positive constants. To estimate the thermalization time, we plug the fixed point couplings into the above RG equation, and get the solution,

[TABLE]

Because $l=\log\Lambda t$, the thermal time scale is given by

[TABLE]

At the near-equilibrium chiral Ising fixed point, the thermalization time is

[TABLE]

where $\epsilon=3-d$. Apparently the thermalization time is extremely long at large $N$, even when $\epsilon=1$. Indeed, the thermalization time scale is exponentially long $t_{\text{th}}\propto(1+\frac{3}{16\alpha_{3}})^{N}$ for $\epsilon=1,N\gg 1$. On the other hand, at the dynamical chiral Ising fixed point, the thermalization time is

[TABLE]

where $\epsilon=4-d$. Distinct from the near-equilibrium chiral Ising fixed point, thermalization time scale is only controlled by small $\epsilon$.