Thermalization and chaos in QED$_{3}$

TL;DR

This paper investigates the thermalization process and quantum chaos in a 2+1 dimensional QED with many fermion flavors, revealing exponential growth in fermionic anti-commutators and calculating the Lyapunov exponent.

Contribution

It provides the first leading-order calculation of fermionic anti-commutator growth and quantum Lyapunov exponent in this gauge theory, linking chaos to conformal field theory behavior.

Findings

Anti-commutators grow exponentially, indicating chaos.

Calculated quantum Lyapunov exponent for the theory.

Discussed implications for locality and gauge invariance.

Abstract

We study the real time dynamics of $N_F$ flavors of fermions coupled to a $U(1)$ gauge field in $2+1$ dimensions to leading order in a $1/N_F$ expansion. For large enough $N_{F}$, this is an interacting conformal field theory and describes the low energy properties of the Dirac spin liquid. We focus on thermalization and the onset of many-body quantum chaos which can be diagnosed from the growth of initally anti-commuting fermion field operators. We compute such anti-commutators in this gauge theory to leading order in $1/N_F$. We find that the anti-commutator grows exponentially in time and compute the quantum Lyapunov exponent. We briefly comment on chaos, locality, and gauge invariance.

Figures

Equations

$$|\delta \mathbf{x}(t)|=e^{{\lambda }_{L}t}|\delta \mathbf{x}(0)|.$$

$$\mathcal{F}(t,\mathbf{x})=⟨\{W(t,\mathbf{x}),V(0)\}\{W(t),V(0){\}}^{†}{⟩}_{\beta }$$

$$\mathcal{F}(t,\mathbf{x})\sim ϵe^{{\lambda }_L\left(t-\frac{|\mathbf{x}|}{v_B}\right)}$$

$$\mathcal{L}=-\frac{1}{4e^2{N}_f}{F}^2+\sum _{a=1}^{N_F}i{\bar{\psi }}_a(x){\gamma }^{\mu }{D}_{\mu }{\psi }_a(x)$$

$$\mathcal{F}(t,\bm{x})=\frac{1}{N^{2}_{f}}\sum^{N_{F}}_{a,b=1}\text{Tr}\Big{[}\sqrt{\hat{\rho}}\{\psi_{a}(t,\bm{x}),\bar{\psi}_{b}(0)\}\sqrt{\hat{\rho}}\{\psi_{a}(t,\bm{x}),\bar{\psi}_{b}(0)\}^{\dagger}\Big{]}.$$

$${\Gamma }_{\mathbf{k}}\sim \frac{T}{N_F}$$

$${\lambda }_L=\frac{2\pi T}{N_F}\mathcal{C}$$

$$k^2=(k^0{)}^2-{\mathbf{k}}^2$$

$$K^2=k_n^2+{\mathbf{k}}^2$$

$${\int }_x$$

$${\int }_k$$

$${\int }_X$$

$$\displaystyle\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{K}$$

$$\mathcal{Z}=\int \mathcal{D}A\mathcal{D}\psi \mathcal{D}\bar{\psi }\text{exp}\left(i{S}_M\right)$$

$${\mathcal{L}}_E=-\frac{1}{4e^2{N}_F}{F}^2+\sum _{a=1}^{N_F}{\bar{\psi }}_a(X){\tilde{\gamma }}_{\mu }{D}_{\mu }{\psi }_a(X)$$

$${\mathcal{Z}}_E=\int \mathcal{D}A\mathcal{D}\psi \mathcal{D}\bar{\psi }\text{exp}\left(-{\int }_X{\mathcal{L}}_E\right)$$

$${\tilde{G}}^{ab}(t,\mathbf{x})$$

$$G_E^{ab}(\tau ,\mathbf{x})$$

$$G_R^{ab}(t,\mathbf{x})$$

$$G_W(/k)$$

$$G_E(/K)$$

$$G_R(/k)$$

$$/\rho (k)=-2\Im [G_R(/k)]$$

$$/\rho (k^0,\mathbf{k})=\frac{\pi /k}{E_{\mathbf{k}}}[\delta (k^0-E_{\mathbf{k}})-\delta (k^0+E_{\mathbf{k}})]$$

$$G_W(/k)$$

$$G_R(/k)$$

$${\mathcal{D}}_W^{\mu \nu }(t,\mathbf{x})$$

$${\mathcal{D}}_E^{\mu \nu }(\tau ,\mathbf{x})$$

$${\mathcal{D}}_R^{\mu \nu }(t,\mathbf{x})$$

$${\mathcal{D}}_W^{\mu \nu }(k)$$

$${\mathcal{D}}_E^{\mu \nu }(K)$$

$${\mathcal{D}}_R^{\mu \nu }(k)$$

$${\mathcal{A}}^{\mu \nu }=-2\Im {\mathcal{D}}_R^{\mu \nu }(k)]$$

$${ϵ}_T^{\mu }(k)=\frac{1}{|\mathbf{k}|}\left(0,-k^y,k^x\right)$$

$$i{\mathcal{D}}^{\mu \nu }(k^0,\mathbf{k})$$

$$P^{\mu\nu}_{T}(k)=\Big{(}\delta^{ij}-\frac{k^{i}k^{j}}{|{\bm{k}}|^{2}}\Big{)}g^{\mu}_{i}g^{\nu}_{j}$$

Full text

.

1 Introduction

Remarkably, closed quantum systems exhibiting chaotic dynamics can act as their own heat bath, leading to an effective thermal state at late time [1][2][3][4]. One of the great challenges in condensed matter physics is to understand this process of thermalization in closed many-body systems in the vicinity of interacting quantum critical points. These systems have fast dynamics and relax on a timescale set by the temperature, $\tau_{\varphi}\sim\frac{\hbar}{k_{B}T}$. At early times, the timescales of thermalization can be probed via the time dependence of simple local correlation functions. At intermediate times, a process called “scrambling” spreads information about the initial state or perturbation throughout the entire system [5][6][7][8]. Scrambling after local relaxation has occurred is not easily accessible in local correlation functions, as it must be described by objects that probe degrees of freedom of the entire system.

In quantum systems with many degrees of freedom, the rate of scrambling can be characterized by a quantum Lyapunov exponent which controls the onset of quantum chaos [9][5][10]. In classical systems, one aspect of chaos is the divergence of initially nearby phase space trajectories. This gives a schematic definition for a classical Lyapunov exponent $\lambda_{L}$ as

[TABLE]

An analogous notion in quantum systems is obtained from the behavior of an out-of-time-order correlation function (OTOC) [9][5][10]. This object has an elegant description in AdS/CFT as the correlation between two entangled black holes connected by a wormhole [5][11]. In fermionic systems, this object is contained in the regulated squared anti-commutator of two fermionic operators. We denote this object as $\mathcal{F}(t,\bm{x})$ which is defined as follows

[TABLE]

where $V$ and $W$ are two fermionic operators separated in space-time.

The squared anti-commutator can be understood as characterizing the effect that one local measurements has on another. For spatially local systems, causality implies that initially separated $V$ and $W$ must have zero anti-commutator (since we are dealing with fermions). However, the squared anti-commutator grows exponentially with time in systems with many local degrees of freedom, and it saturates to its late time value after the “scrambling” time $t_{\ast}$. During this process, the anti-commutator behaves like

[TABLE]

where $\lambda_{L}$ is identified as a Lyapunov exponent in analogy with classical chaos, and $v_{B}$ is called the butterfly velocity which acts as an effective speed limit. We emphasize that this is expected to hold when there are a large number of local degrees of freedom. More generally, it has been argued recently that the exponential form above is modified by quantum fluctuations in a non-trivial way. Here we focus for simplicity on the spatially averaged object, $\int d^{2}x\mathcal{F}(\bm{x},t)$, which diagnoses scrambling within the local Hilbert space.

Quantum Lyapunov exponents have by now received intense scrutiny. They can be measured experimentally [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22], are related to operator growth [23] [24] [25] [26] [27] [28] [29] [30], and can be computed either numerically or analytically in many model systems [10] [23] [24] [25] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41][42] [43] [44] [45] [46] [47] [48] [49] [50].

In particular, chaos exponents have been computed for several interesting field theories including scalar matrix theory, the O(N) model, the critical fermi surface, the diffusive metal, and graphene [51]. In this paper, we study many-body chaos in $\text{QED}_{3}$ with $N_{f}$ flavors of massless fermions in the large $N_{f}$ limit. This theory is described by the Lagrangian

[TABLE]

where we define the covariant derivative $D_{\mu}=\partial_{\mu}-\frac{iA_{\mu}}{\sqrt{N_{F}}}$ and the Lorentzian Dirac matrices are given by $\gamma^{\mu}=(\sigma^{z},i\sigma^{x},i\sigma^{y})$ and we work in the $+--$ signature. $\text{QED}_{3}$ appears throughout condensed matter as the proposed low energy theory for several lattice systems such as the the Kagome antiferromagnet for $N_{f}=2$. It also relevant for a variety of quantum critical points [52] [53]. While these lattice models also contain monopole operators which can drastically modify the long-wavelength physics, for large enough $N_{f}$ these become irrelevant and the theory flows to a strongly coupled renormalization group fixed point. It is known that the large $N_{f}$ expansion breaks down below some critical value of $N_{f}$, below which the theory is confining [54][55][56][57].

Our motivations for studying this model are the following. First, we wanted to understand possible subtleties in the quantum Lyapunov exponent story associated with gauge invariance. In particular, the anti-commutator by itself is not gauge invariant, although the squared anti-commutator is. The anti-commutator can be made gauge invariant by attaching a Wilson line, so this non-locality might conceivably modify the chaos story. Second, it is interesting to understand how to detect failures of the large $N_{f}$ expansion. One idea is that the chaos bound of Ref. [58] might be in tension with the large $N_{f}$ result for $\lambda_{L}$, indicating a transition. Of course, the leading $1/N_{f}$ term is not reliable away from large $N_{f}$, but it might still be suggestive of an issue. Third, we wanted to compare chaos in conventional and unconventional quantum critical points, for example, how does chaos in the O(N) model compare to chaos in $\text{QED}_{3}$?

To analyze chaos in $\text{QED}_{3}$, we choose to study the index averaged squared anti-commutator

[TABLE]

We restrict attention to the spatial average of $\mathcal{F}$ to avoid issues associated with spatial propagation and focus on chaos within the large local Hilbert space. However, our results could be easily extended to the compute the butterfly velocity at large $N_{f}$.

We briefly summarize our main results. We numerically calculate the photon polarization bubble in real time and analyze analytically it in the limit of high temperature and the limit of low frequency. We use it to numerically calculate the fermion decay rate which scales as

[TABLE]

and vanishes as ${\bm{k}}\rightarrow 0$. The Lyapunov exponent takes the form

[TABLE]

with $\mathcal{C}\approx 0.65$. This form is analogous to the form of $\lambda_{L}$ in the O(N) model with here the value of $\mathcal{C}$ is slightly higher. Because $\mathcal{C}<1$, the chaos bound is satisfied for all values of $N_{F}$ despite the known breakdown of the large $N_{F}$ expansion. Also, chaos is not significantly stronger in this unconventional quantum critical point as compared to the usual O(N) critical theory, at least at large $N$.

Notation guide: Throughout we will use lower case letters to denote Lorentz vectors $k^{\mu}=(k^{0},{\bm{k}})$ where the inner product is defined as

[TABLE]

For Euclidean vectors we will use capital letter $K_{\mu}=(k_{n},{\bm{k}})$ vectors and the Euclidean inner product is defined as

[TABLE]

Integrals in Minkowski space are defined as

[TABLE]

Integrals and sums in Euclidean space are defined as are defined

[TABLE]

For Feynman slash notation $\not{k}\equiv\gamma^{\mu}k_{\mu}$, we will use the shorthand $\not{k}^{0}=\gamma^{0}k_{0}$, $\not{{\bm{k}}}=\gamma^{i}k_{i}$ and $\not{k}=\gamma^{0}k_{0}+\gamma^{i}k_{i}=\gamma^{0}k_{0}-\bm{\gamma}\cdot{\bm{k}}$ where $\bm{\gamma}=(i\sigma^{x},i\sigma^{y})$. In Euclidean time $\not{K}\equiv\tilde{\gamma}^{\mu}K_{\mu}$ $\not{K}=\tilde{\gamma}^{0}ik_{n}+\tilde{\gamma}^{i}{\bm{k}}_{i}$ To go from to $\not{K}$ back to $\not{k}$ we have $i\not{K}\rightarrow\not{k}$. Further conventions are specified in the appendices.

2 Preliminaries

In this section, we set up the formalism required to compute $\mathcal{F}(t,\bm{x})$ in a $1/N_{f}$ expansion The generating function corresponding to Eqn. 1.4 is

[TABLE]

In the limit $e\rightarrow\infty$ , we can neglect the Maxwell term.

It will also be useful to consider the theory in Euclidean time. Setting $it\rightarrow\tau$ and $A_{0}\rightarrow iA_{\tau}$ we obtain the following Euclidean action

[TABLE]

where the Euclidean and Minkowski Dirac matrices are related by $\tilde{\gamma}^{0}\equiv\gamma^{0}$ and $\tilde{\gamma}^{i}\equiv-i\gamma^{i}$ and $\tilde{\gamma}^{i}=\tilde{\gamma}_{i}$. We will use repeated lower indices and capital letters when using the Euclidean metric This gives us the generating functional

[TABLE]

where $\int_{X}\equiv\int^{\beta}_{0}\mathrm{d}\tau\int_{\bm{x}}$.

To compute $F(t,\bm{x})$, we will need the Wightman, Euclidean, and Euclidean propagators for the fermions and the gauge field. For fermions these are defined as

[TABLE]

In momentum space these propagators can be written in the spectral representation

[TABLE]

where we have defined the spectral function as

[TABLE]

For free fermions the spectral function is

[TABLE]

which gives us the retarded and Wightman propagators

[TABLE]

Similarly, for the gauge field we define the propagators as

[TABLE]

In the spectral representation we can write these as

[TABLE]

where we define the spectral function as

[TABLE]

The gauge field propagator $\mathcal{D}^{\mu\nu}$ is $O(1)$ in $N_{f}$ and hence a non-perturbative object. Since the self energy dominates at low energies and $\mathcal{D}^{\mu\nu}$ completely in terms of the self energy.

At finite temperature, the choice of frame breaks Lorentz invariance. As a result, longitudinal and transverse become distinct. In $2+1d$, the photon has one longitudinal and one transverse polarization with associated polarization vectors $\epsilon^{\mu}_{L}(k)$ $\epsilon^{\mu}_{T}(k)$. A thermal mass can appear in the longitudinal component of the longitudinal part without violating gauge invariance. To make calculation simpler we will use Coulomb gauge throughout with gauge fixing condition $\nabla\cdot\bm{A}=0$. In this gauge, we define the longitudinal polarization vector $\epsilon^{\mu}_{L}(k)=u^{\mu}=(1,0,0)$, where $u^{\mu}$ is the rest frame of the medium. For ${\bm{k}}$ lying along the $x$ axis we can write the transverse polarization vector as following

[TABLE]

In this gauge, we can write the bare propagator in the following form

[TABLE]

where $P^{\mu\nu}_{T}(k)$ is a transverse projection tensor defined as

[TABLE]

We can also write $P^{\mu\nu}_{T}(k)$ as a tensor product of the transverse polarization vectors

[TABLE]

For the self energy, current conservation requires

[TABLE]

This allows us to write it with the following tensor structure.

[TABLE]

where longitudinal projector is defined as

[TABLE]

Since the tensor 2.17 has only two independent components, we can calculate everything in terms of the spatial components of the self energy $\Pi^{xx}(k^{0},{\bm{k}})$, $\Pi^{yy}(k^{0},{\bm{k}})$ and write $\Pi^{00}(k^{0},{\bm{k}})$ in terms of these functions.

We start with the following expression for $\mathrm{\Pi}_{\mu\nu}(K)$ in Euclidean time

[TABLE]

For For Euclidean gamma matrices we have

[TABLE]

which gives us

[TABLE]

We can simplify the $K^{\mu}P^{\nu}+K^{\nu}P^{\mu}$ terms by shifting $P\rightarrow K-P$ in half of this term to get

[TABLE]

We leave the rest of the calculation of the polarization bubble to appendix A.

We write the full photon propagator in terms of $\Pi^{00}(k^{0},{\bm{k}})$ and $\Pi^{yy}(k^{0},{\bm{k}})$ the as

[TABLE]

Since the self energy is $O(1)$, it will always dominate over the bare term. Thus we drop the bare term and write

[TABLE]

At finite temperature, a thermal mass $m_{thermal}=\frac{\sqrt{e}T\log 2}{\pi}$ will appear in the longitudinal part of the propagator $\Pi^{00}$. This does not violate gauge invariance and is due to the Debye screening of the gauge field.

Next we consider the first correction to the fermion self energy. Since the fermion mass is $1/N_{f}$ suppressed, to leading order the physical effects are controlled by the imaginary. This introduce a damping rate $\Gamma_{\textbf{k}}$ leading to a finite lifetime to the fermions. In terms of the photon and free fermion spectral functions, we can write the self energy as

[TABLE]

Doing the Matsubara sum gives us

[TABLE]

Analytically continuing back to real time and inserting the free fermion spectral function gives us

[TABLE]

where $p_{s}=(sE_{\bm{p}},\bm{p})$ and $E_{\bm{p}}=|\bm{p}|$. Defining $\Gamma_{{\bm{k}}}$ as the imaginary part of the pole in the propagator, $\Gamma_{{\bm{k}}}$ can be written in terms of the self energy as

[TABLE]

Taking the imaginary part of Eq. 2 we have

[TABLE]

Expanding the photon propagator into its longitudinal and transverse components, this becomes

[TABLE]

We can evaluate Eq. 2 in Coulomb gauge by choosing the longitudinal polarization vector $\epsilon^{\mu}_{L}=u^{\mu}$. The complete expression is given in B. We also note that this expression is consistent with the finite temperature optical theorem. Implementing $\Gamma_{{\bm{k}}}$ as a shift in the pole, we can write the retarded propagator as

[TABLE]

Note that the fermion 2-point function is not gauge invariant, so this decay rate is not manifestly gauge invariant and is defined in the Coulomb gauge which we use. It is consistently used below in the $1/N_{F}$ expansion to compute gauge-invariant observables, and it presumably manifests in other gauge-invariant response functions.

3 Diagrammatic expansion for the squared fermion anti-commutator

We now consider the diagrammatic expansion of the squared anti-commutator of two fermion operators by expanding $\mathcal{F}(t,\bm{x})$ in the interaction on both time folds and keeping only diagrams to order $1/N_{F}$ which are uncrossed. The rules are very similar to those summarized in [59] except all of the propagators and hence rungs have a matrix structure. In summary:

1. Vertex insertions are restricted to lie on real time folds

2. Horizontal lines correspond to retarded or advanced photon or fermions propagators including self energy corrections

3. Vertical lines correspond to full photon propagators or bare fermion Wightman propagators

4. Each rung or pair of rails is a matrix in spinor space and can be written in terms of a tensor product basis of gamma matrices. A rung or pair of rails is added to the ladder by matrix multiplication.

We now consider the two rungs at this order. The first rung involves a single photon Wightman function sandwiched between two gamma matrices. The whole rung can be giving us the following expression

[TABLE]

where we define $\mathcal{Q}_{B}(\omega)=(2\sinh(\frac{\beta\omega}{2}))^{-1}$ for bosons and $\mathcal{Q}_{F}(\omega)=(2\cosh(\frac{\beta\omega}{2}))^{-1}$ for fermions. We can write the longitudinal and transverse contributions to $\mathcal{R}_{1}(k-k^{\prime})$ in terms of a tensor products of polarization vectors. In $2+1d$, we can express the transverse projector $P^{\mu\nu}_{T}(k)$ as the tensor product of a single transverse polarization vector

[TABLE]

where we define $\epsilon^{\nu}_{T}({\bm{k}})$ lying along the $x$ axis as

[TABLE]

Defining $\bar{k}=k-k^{\prime}$ as the relative momentum, we can write a general expression for $\mathcal{R}_{1}(\bar{k})$

[TABLE]

where we have separated the longitudinal and transverse parts of the rung. Using the definitions of the longitudinal and transverse polarization vectors this becomes

[TABLE]

The second rung is a box structure consisting of two fermion Wightman functions and two retarded photon propagators. We can write this as

[TABLE]

where $p^{\mu}=(p^{0},\textbf{p})$ Inserting the free fermion spectral functions, we get

[TABLE]

defining the absolute Lorentz momentum $K=\frac{k+k^{\prime}}{2}$ corresponding to the Lorentz vector $K^{\mu}=(K^{0},\bm{K})$ and the relative momentum $\bar{k}=k-k^{\prime}$ corresponding to the Lorentz vector $\bar{k}^{\mu}=(\bar{k}^{0},\bar{{\bm{k}}})$. Shifting $p\rightarrow p+K$, we define $E_{\pm}=E_{\bar{{\bm{k}}}/2\pm\bm{p}}$ given by

[TABLE]

where $\theta$ is the angle between $\bar{{\bm{k}}}$ and $\bm{p}$ This gives us

[TABLE]

performing the $p^{0}$ integral in Eq. 3, we get

[TABLE]

where $\tilde{p}^{\mu}=(\bar{k}^{0}/2-s^{\prime}E_{-},\bm{p})$. In three dimensions $\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}=g^{\mu\nu}\gamma^{\rho}-g^{\mu\rho}\gamma^{\nu}+g^{\nu\rho}\gamma^{\mu}$ which allows us to write Eqn. 3 in terms of the polarization vectors as follows

[TABLE]

Since the external momentum ${\bm{q}}$ is $O\left(\frac{1}{N_{F}}\right)$, the momenta of the two photon propagators are approximately anti-parallel. Therefore, we only need to consider terms in the sum over $\lambda$ where $\lambda=\lambda^{\prime}$.

We can now write the whole rung in the tensor product form

[TABLE]

where we define the coefficients of $\gamma_{\mu}\otimes\gamma_{\nu}$ as

[TABLE]

The expressions for the matrix elements defined in 3 are contained in appendix.

We can do the integral over $|\bm{p}|$ where the zeroes of the delta function are given by

[TABLE]

where $\theta_{\bm{p}}$ is the angle between the loop momenta $\bm{p}$ and the relative momenta $\bar{{\bm{k}}}$.

For the root to be positive, we have condition $(k^{0})^{2}<|{\bm{k}}|^{2}\cos^{2}\theta$ or $(k^{0})^{2}>|{\bm{k}}|^{2}$, which occurs for $(s,s^{\prime})=(1,1)$ and $(s,s^{\prime})=(\pm 1,\mp 1)$.

4 Bethe Salpeter equation and the Lyapunov exponent

Our diagrammatic expansion of the squared fermion anti-commutator will give us an integral equation we will solve to obtain the Lyapunov exponent and butterfly velocity. This equation is reminiscent of the original Bethe-Salpeter equations first used to study bound states in quantum electrodynamics in [60][61][62].

In matrix form, the equation for the ladder is basis dependent. However it, turns out that at long times, we can take advantage of the matrix structure of this equation and turn it into two scalar equations for components of a “wave-function” To illustrate this, we start define the Lorentz vectors $q^{\mu}=(\nu,{\bm{q}})$ $k^{\mu}=(k^{0},{\bm{k}})$ and the following on-shell projectors

[TABLE]

where $i=1,2$ refers to the particle index and $\bm{\gamma}\cdot{\bm{k}}=-\gamma^{i}k_{i}$. These have the following properties

[TABLE]

The eigenvectors of $\Lambda^{(i)}_{\pm}({\bm{k}})$ are

[TABLE]

From Eqns. 4.3a and 4.3b we see $u_{-}({\bm{k}})=u_{+}(-{\bm{k}})$. Acting with the projectors on $u_{\pm}({\bm{k}})$ we get

[TABLE]

From this, we conclude that $u_{\pm}({\bm{k}})$ are eigenvectors of $\Lambda_{\pm}({\bm{k}})$ with eigenvalue one and $u_{\pm}({\bm{k}})$ are an eigenvectors of $\Lambda_{\mp}({\bm{k}})$ with eigenvalue zero. We define $\Lambda^{(1)}_{s}({\bm{k}})\Lambda^{(2)}_{s^{\prime}}({\bm{k}})$ as the tensor product of two projection operators. Eq. 4.2 gives us the additional four properties

[TABLE]

The eigenvectors of $u^{(1,2)}_{\pm}$. $\Lambda^{(1)}_{s}({\bm{k}})\Lambda^{(2)}_{s^{\prime}}({\bm{k}})$ can be written as tensor products of the form

[TABLE]

We can use the free fermion spectral function to write the retarded Green’s function in the form $\not{k}g_{R}(k)$.

[TABLE]

Consider the product $g^{\ast}_{R}(q-k)g_{R}(k)$. Writing it explicitly in terms of momentum and frequency arguments we have

[TABLE]

Evaluating the $k^{0}$ integral we have four poles: $k^{0}=\nu-sE_{{\bm{k}}-{\bm{q}}}+i\epsilon$ and $k^{0}=s^{\prime}E_{{\bm{k}}}-i\epsilon$ which give residues proportional to $\sim(\nu+sE_{{\bm{k}}-{\bm{q}}}-s^{\prime}E_{{\bm{k}}})^{-1}$. To get the longtime behavior, we keep only the contribution from the most singular terms $\sim(\nu\pm E_{{\bm{k}}-{\bm{q}}}\mp E_{{\bm{k}}})^{-1}$ where we close the contour clockwise on the lower half plane

In the end we have

[TABLE]

Note that for ${\bm{q}}=0$, these terms become proportional to $\nu^{-1}$. Assuming $\Gamma_{{\bm{q}}-{\bm{k}}}\approx\Gamma_{{\bm{q}}}$ we can shift the poles in 4.9 to get

[TABLE]

For small ${\bm{q}}$, we can express Eqn. 4.9 in terms of projection operators as follows

[TABLE]

where we used using $(\gamma^{0})^{2}=\mathbb{1}$. Consider the first term in the ladder, which is the tensor product of two retarded Green’s functions

[TABLE]

In momentum space we can define the object $f(q,k)$ as

[TABLE]

And write the first term as

[TABLE]

We can now write the infinite ladder as

[TABLE]

Where $q^{\mu}=(\nu,{\bm{q}})$, $k^{\mu}=(k^{0},{\bm{k}})$ and $k^{\prime\mu}=(k^{\prime 0},{\bm{k}}^{\prime})$, and $\mathcal{R}(k,k^{\prime})$ is the sum of the two rungs. It is important to note that we have computed the two rungs in a different basis, and we would need to use the Fierz identities to explicitly write $\mathcal{R}(k,k^{\prime})$ in the correct matrix form.

Assuming exponential growth at long times, the ladder will remain invariant under the addition of one rung. For ${\bm{q}}$ small, we can make the approximation $E_{{\bm{k}}}-E_{{\bm{k}}-{\bm{q}}}=\delta E_{{\bm{q}}}\approx{\bm{q}}\cdot v_{{\bm{k}}}$ where $v_{{\bm{k}}}=\nabla_{{\bm{k}}}E_{{\bm{k}}}$ is the group velocity. Dropping the constant term and inserting Eqn. 4.11 into 4.15, we get

[TABLE]

Acting on both sides of Eqn. 4.16 with $(\Lambda^{(1)}_{s}({\bm{k}})\gamma^{0}\Lambda^{(2)}_{s}({\bm{k}})\gamma^{0})$, we see from Eqn. 4.5 that Eqn. 4.16 is just multiplied by one or zero. So instead of considering the matrix $f(q,k)$, we can consider an eigenvalue equation for a vector $v(q,k)$ which is an eigenvector of the projection operators. We can now propose the following on-shell form of $v(q,k)$

[TABLE]

Since $v_{\pm}(\nu,{\bm{q}},{\bm{k}})$ are eigenstates of $\Lambda^{(1)}_{\pm}({\bm{k}})\Lambda^{(2)}_{\pm}({\bm{k}})$ for small $q$ we can write $v^{(i)}_{\pm}(\nu,{\bm{q}},{\bm{k}})$ in the form

[TABLE]

where $\psi_{\pm}(\nu,{\bm{q}},{\bm{k}})$ is a scalar function. To write equations for $\psi_{\pm}(\nu,{\bm{q}},{\bm{k}})$ we take the inner product of both sides, defining a scalar kernel written in terms of both rungs and as

[TABLE]

Because inner product is basis independent, it does not matter that the rungs are written in different bases.

[TABLE]

Eqns. 4.20 now looks like two coupled effective Dirac equations for a scalar wave function. We can now define a kernel matrix in the $ss^{\prime}$ indices as

[TABLE]

and $\Psi(\nu,{\bm{q}},{\bm{k}})=(\psi_{+}(\nu,{\bm{q}},{\bm{k}}),\psi_{-}(\nu,{\bm{q}},{\bm{k}})^{T}$

[TABLE]

In the time domain, this becomes

[TABLE]

5 Numerics

We define the Lyapunov exponent to be the largest eigenvalue of $\hat{K}$ at ${\bm{q}}=0$ with eigenvector $\psi_{1}$. We assume $\psi_{1}$ is rotationally symmetric so the eigenvalue equation can be simplified by projecting both sides onto states with angular momentum $l=0$. The resulting eigenvalue equation involves a one dimensional integral transform, involving only the magnitude of the momentum. Below we summarize the computational procedure, with further details in Appendix E.

First, we compute $\Pi^{xx}$ and $\Pi^{yy}$ numerically on a fine grid in $(\omega,|{\bm{k}}|)$ space and use these to construct the full photon spectral function matrix. We then use these to calculate the fermion decay rate and the matrix components of the rung functions. Since these quantitites effectively vanish when the norms of their arguments greatly exceed T, we use a simple linearly spaced grid with a hard momentum cutoff of order $10$ T . The temperature is the only scale in the problem, so the entire numerical computation is set up in terms of the scaled variables $\mathfrak{k}=k/T$ and $\mathfrak{w}=\omega/T$ .

We can write $\lambda_{L}$ in the form

[TABLE]

Numerics give $\mathcal{C}\approx 0.65$. $\lambda_{L}$ takes the same form for the O(N) model with $\mathcal{C}=0.51$. For nonzero ${\bm{q}}$ we can also calculate a butterfly velocity $v_{B}$ which we expect to be $v_{B}\sim 1$. We do not give a precise number as computing $v_{B}$ is more numerically involved and the spatial propagation is expected to be significantly modified by higher order in $\frac{1}{N_{F}}$ corrections [63].

6 Discussion

In this work we studied thermalization and chaos in $\text{QED}_{3}$ at finite temperature in a large $N_{F}$ expansion. We calculated the photon polarization tensor in real time numerically and also obtained exact expressions in the limits of high temperature and high frequency. We numerically computed the fermion decay rate and obtained the Lyapunov exponent from a two-body Dirac-like equation for the squared anti-commutator.

While the object we have computed, $\mathcal{F}(t,\bm{x})$, is gauge invariant, the fermion anti-commutator by itself is not. Furthermore, because of their anti-commutation relations, fermions are inherently non-local, thus the squared anti-commutator we are computing is actually a non-local object even though the growth is causal.

We can compare the growth of the squared fermion anti-commutator with the growth of local operators. The most natural local operator in $\text{QED}_{3}$ is the fermion current. This is defined as

[TABLE]

We can consider the squared commutator of two current operators

[TABLE]

Another object we could have considered is the squared commutator for the gauge field. Due to current conservation, we know that the photon and current correlation functions are related up to contact terms which vanish on shell.

We expect composite operators will grow as fast as elementary operators [64]. However, it is not guaranteed that these grow at the same rate.

We can compare $\mathcal{C}(t,\bm{x})$ to $\mathcal{F}(t,\bm{x})$ by looking at the structure of the ladder diagrams. The technical details are left to appendix D.1). We see that all of the rungs in $\mathcal{F}(t,\bm{x})$ are also contained in $\mathcal{C}(t,\bm{x})$ as well as additional structures with crossed legs and particle arrows oriented in opposite directions. The internal structure of this ladder contains the same terms in the squared fermion anti-commutator as well as other pieces. If we expand $\mathcal{C}(t,\bm{x})$ further, we obtain diagrams containing particle flows in opposing directions as shown in Fig. 8. In principle, these can also occur in the ladder expansion for $\mathcal{F}(t,\bm{x})$. However, these terms are suppressed by $\frac{1}{N_{F}}$, so they do not contribute to chaos to leading order. Thus, we argue that the diagrams in $\mathcal{F}(t,\bm{x})$ are all contained in $\mathcal{C}(t,\bm{x})$, so chaos must grow at least as fast as the rate obtained from $\mathcal{F}(t,\bm{x})$.

In the future it would be interesting to see if one could adapt a kinetic method such as described in [45] to simplify the calculation. It would also be interesting to improve numerics so that the Lyapunov exponent can be calculated to higher accuracy and get an estimate for the butterfly velocity. Our work also serves as a starting point to compute other dynamical quantities in finite temperature $\text{QED}_{3}$ such as the charge diffusion coefficient in systems with disorder. While the decay rate we computed is not gauge invariant, it is the first step in the computation of physical response functions.

7 Acknowledgements

We thank G. Cheng, D. Chowdhury, M. Pate, A. Patel, S. Sachdev, and S. Whitsitt for valuable disussions. We espcially thank G. Cheng for a valuable check on our calculations. JS is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1144152. BGS acknowledges support from the Simons Foundation as part of the It From Qubit Collaboration.

Appendix A Real time Polarization bubble

In this appendix we will describe the calculation of the components photon polarization tensor

We start with the expression for the spatial components in Euclidean time

[TABLE]

where $K$ and $P$ are Euclidean vectors. Note that the numerator does not involve the summed-over Matsubara frequencies. The Matsubara sums can be performed to give

[TABLE]

and

[TABLE]

The first term, after subtracting the zero temperature value, gives

[TABLE]

where $\bm{k}=(|{\bm{k}}|,0)$ and $\bm{p}=|\bm{p}|(\cos\theta,\sin\theta)$. The numerators of the second term are

[TABLE]

and

[TABLE]

A.1 Tensor structure

The full real-time polarization can be expanded in two tensor structures. While the transverse component is spatially transverse, the longitudinal component is spacetime transverse. This tensor structure ensures that the full polarization bubble obeys the Ward identity required by current conservation.

At finite temperature there is a new vector available set by the rest frame of the medium. The tensor structures are

[TABLE]

and

[TABLE]

From these two structures we can define

[TABLE]

These obey

[TABLE]

[TABLE]

and

[TABLE]

To give an invariant definition of $P^{\mu\nu}_{T}$, let $u^{\mu}$ denote the velocity of the medium, which in the medium’s rest frame is $u=(1,0,0)$. Let the transverse $u$ projector be

[TABLE]

and define

[TABLE]

In the medium’s rest frame, $k=(k^{0},\bm{k})$, then

[TABLE]

The transverse projector to $\bar{k}$ is

[TABLE]

which in the medium’s rest frame is

[TABLE]

[TABLE]

and

[TABLE]

Hence we may identify $P^{\mu\nu}_{T}$ with $u^{\mu}u^{\nu}+P^{\mu\nu}_{T}(\bar{k})$ (with $\tilde{P}^{\mu\nu}$ with $\tilde{P}^{\mu\nu}(k)$). We write the full polarization as

[TABLE]

The non-vanishing terms are

[TABLE]

[TABLE]

and

[TABLE]

For kinematics where $\bm{k}=(|{\bm{k}}|,0)$, these expressions simplify to

[TABLE]

[TABLE]

and

[TABLE]

Hence $\Pi^{xx}_{R}$ and $\Pi^{yy}_{R}$ are sufficient to determine the two functions $\Pi$ and $\tilde{\Pi}$.

A.2 Imaginary part, full expressions

The first contribution to both $\Pi^{xx}_{R}$ and $\Pi^{yy}_{R}$ is the constant term,

[TABLE]

where $\Pi_{T=0}$ is the zero temperature contribution. Because this contribution is pure real, it drops out of the imaginary part.

Thus the full imaginary part of $\Pi^{ij}_{R}$ can be obtained from

[TABLE]

where the numerator is

[TABLE]

Using kinematics where $\bm{k}=(|{\bm{k}}|,0)$ and $\bm{p}=(|\bm{p}|\cos\theta,|\bm{p}|\sin\theta)$, we find

[TABLE]

and

[TABLE]

We can check that $\Im\Pi^{ij}_{R}$ is an odd function of frequency. For positive $k^{0}$, only the $(1,1)$, $(1,-1)$, and $(-1,1)$ terms contribute. Furthermore, the $(1,-1)$ and $(-1,1)$ terms give the same result. The $(1,1)$ contribution is

[TABLE]

The $(1,-1)$ contribution (counted twice) is

[TABLE]

The delta function condition (with $s=1$) is

[TABLE]

To solve it, let us square both sides,

[TABLE]

The result is

[TABLE]

or

[TABLE]

We need either $k^{0}>|{\bm{k}}|$ or $k^{0}<|{\bm{k}}|\cos\theta$ to have a positive solution for $|\bm{p}|$. The former is realized for $ss^{\prime}=1$ while the latter is realized for $ss^{\prime}=-1$. To complete the evaluation of the delta function, we also need

[TABLE]

The $ss^{\prime}=1$ contribution is

[TABLE]

The $ss^{\prime}=-1$ contribution is

[TABLE]

The real part is then obtained from a subtracted Kramers-Kronig relation. For a function $\chi(\omega)=\chi_{1}+i\chi_{2}$ obeying the assumptions of Kramers-Kronig, the real part in terms of the imaginary part is

[TABLE]

A.3 Low temperature limit

Suppose $\omega\gg T$ and $|k|=0$. Only the $ss^{\prime}=1$ term contributes to $\Im\Pi^{ij}_{R}$. The contribution is

[TABLE]

with $L^{xx}\rightarrow(k^{0})^{2}-4|\bm{p}|^{2}\cos^{2}\theta$ and $L^{yy}\rightarrow(k^{0})^{2}-4|\bm{p}|^{2}\sin^{2}\theta$. Because of the simple dependence on $\theta$, we conclude the $\Im\Pi^{xx}_{R}=\Im\Pi^{yy}_{R}$ in this limit. (Physically, since $|{\bm{k}}|=0$, there is no difference between longitudinal and transverse.)

We find

[TABLE]

A.4 High temperature limit

In the high temperature limit, $k/T\ll 1$, we neglect $|{\bm{k}}|$ and $k^{0}$ relative to $T$ and the loop momentum whenever the dependence on $k$ is non-singular. The numerators limit to

[TABLE]

and

[TABLE]

The $xx$ polarization is

[TABLE]

Both $ss^{\prime}=1$ terms inside the parenthesis contribute the same value (in the small $k/T$ limit, both are real), as do both $ss^{\prime}=-1$ terms. These are

[TABLE]

and

[TABLE]

where $z=(k^{0}+i\epsilon)/|k|$. Further processing the $ss^{\prime}=-1$ term gives

[TABLE]

Then we see that the constant term, the $ss^{\prime}=1$ terms, and the constant part of the $ss^{\prime}=-1$ terms combine to give zero. Thus the $xx$ polarization is

[TABLE]

The only difference between $xx$ and $yy$ is the replacement $\cos^{2}\theta\rightarrow\sin^{2}\theta$ when going from $L^{xx}$ to $L^{yy}$. Hence it follows that

[TABLE]

The two integrals we must calculate are

[TABLE]

and

[TABLE]

The sum of $I_{1}$ and $I_{2}$ is

[TABLE]

where we converted the expression to a contour integral around the unit circle. The integrand has poles at $w=w_{\pm}=z\pm\sqrt{z^{2}-1}$. The $w=w_{-}$ pole is inside the unit circle for $\Im z>0$. The integral is thus

[TABLE]

The difference of $I_{1}$ and $I_{2}$ is

[TABLE]

where again we converted to a contour integral. The integrand now has poles at $w=0$ (double pole) and at $w=w_{\pm}=z\pm\sqrt{z^{2}-1}$. The integral is thus

[TABLE]

Simplifying gives

[TABLE]

or using $w_{-}^{2}+1=2zw_{-}$,

[TABLE]

Hence

[TABLE]

or

[TABLE]

Similarly,

[TABLE]

or

[TABLE]

Appendix B Fermion damping rate

We can verify the expression for the fermion damping directly from the optical theorem using the finite temperature cutting rules [65]. This relates the decay rate to the imaginary part of the tree level cross-section which in turn is related to the imaginary part of the self-energy. We can define $\mathcal{M}_{L}$ and $\mathcal{M}_{T}$ as the longitudinal and transverse matrix elements for the tree level scattering processes of an electron emitting a soft photon

[TABLE]

where $\bar{u}(p)$ and $u(k)$ are massless on shell spinors. We can write $\Gamma_{{\bm{k}}}$ in terms of these matrix elements as

[TABLE]

where $k=(k^{0},{\bm{k}})$ and $p_{s}=(sE_{\bm{p}},\bm{p})$. Plugging $\mathcal{M}_{L(T)}$ into B and using the relation $u(k)\bar{u}(k)=\not{k}$ we get

[TABLE]

we find the relation for the self energy and the decay rate given in Eqn. 2.60. We can now write the decay rate $\Gamma_{{\bm{k}}}$ in Coulomb gauge as

[TABLE]

Defining the Lorentz vectors $k=(k^{0},|{\bm{k}}|,0)$, $p_{s}=(sE_{\bm{p}},|\bm{p}|\cos\theta,|\bm{p}|\sin\theta)$, and the transverse polarization vector $\epsilon^{\mu}_{T}({\bm{k}}-\bm{p})=\frac{1}{|{\bm{k}}-\bm{p}|}\left(0,|\bm{p}|\sin\theta,|{\bm{k}}|-\bm{p}\cos\theta\right)$, we get the final expression for $\Gamma_{{\bm{k}}}$ in Coulomb gauge

[TABLE]

Appendix C Full expressions for the scalar kernel

Using the definition in Eq. 4.19, we can write an expression for the scalar kernel in the Bethe-Salpeter equation in terms of the rung matrix elements.

[TABLE]

where $\theta_{{\bm{k}}}$ is the angle between ${\bm{k}}$ and the $x$ axis and $\theta_{{\bm{k}}^{\prime}}$ is the angle between ${\bm{k}}^{\prime}$ and the $x$ axis.

If we make a rotationally symmetric ansatz for the wavefunction, only terms depending on $\theta_{{\bm{k}}^{\prime}}-\theta_{{\bm{k}}}$ remain, which leaves us

[TABLE]

Both rungs can be written in the form

[TABLE]

For the first rung, we can write these coefficients as

[TABLE]

where $\theta_{{\bm{k}},{\bm{k}}^{\prime}}$ is the angle between ${\bm{k}}$ and ${\bm{k}}^{\prime}$. For the second rung, we use the general expression given in 3. Defining $\theta_{\bm{p}}$ as the angle between $\bm{p}$ and $\bar{{\bm{k}}}$ and $\phi$ as the angle between $\bar{{\bm{k}}}$ and $\bm{K}$, the transverse polarization vector with respect to the $\bar{{\bm{k}}}$ axis can be written as

[TABLE]

We use these to define the following Lorentz vectors which are contained in the numerators of the longitudinal and transverse components of the second rung in

[TABLE]

Plugging C.5 into C.6, we can write the spatial components of the Lorentz vectors for the longitudinal polarizations as

[TABLE]

and the transverse polarizations as

[TABLE]

We can now now plug the expressions in C.8 into 3 and evaluate the momentum integrals to get complete expressions for the matrix coefficients of rung II.

Appendix D Feynman rules for ladder diagrams

We derive the Feynman rules for the ladder diagrams in the $\frac{1}{N_{F}}$ expansion. Ignoring thermal corrections, we can start by going into the interaction picture with respect to the the free Hamiltonian. The time evolution operator is given by

[TABLE]

here $A$ is the free photon field and $j_{\mu}$ is the free fermion current $j^{\mu}(t)=\bar{\psi}(t)\gamma^{\mu}\psi(t)$. We can expand this (dropping spatial arguments) as

[TABLE]

for $U^{\dagger}_{I}$ we have

[TABLE]

We start by stating the usual anti-commutation relations for free fermions

[TABLE]

For currents we have

[TABLE]

The current-fermion and current-current commutators can be written in terms of the retarded propagator terms of the retarded propagator is now

[TABLE]

Lastly for free fields, both $j^{\mu}_{a}(t,\bm{x})$ and $\tilde{\psi}(t,\bm{x})$ commute with the photon operator so we have

[TABLE]

D.1 Fermion Ladder

In the interaction picture, $\mathcal{F}(t,\bm{x})$ is given by

[TABLE]

where $U^{\dagger}_{I}\tilde{\psi}(t)U_{I}$ is given by the expansion

[TABLE]

From Eqn. D.1 and using the definition of the retarded propagator, we can immediately see that the zeroth order term in $\mathcal{F}(t,\bm{x})$ is given by

[TABLE]

The first nontrivial term is given by

[TABLE]

Expanding the anti-commutator for each time fold, we get

[TABLE]

Putting everything together gives us

[TABLE]

which we identify as a ladder with one type one rung.

Now we can expand to second order on both time folds. We consider the expansion of the commutator contained inside the anti-commutator on one time fold

[TABLE]

From the first term we get

[TABLE]

which will give us the type II rung. This can be written as

[TABLE]

The second term gives

[TABLE]

This will give two copies of the type I rung or a self energy term for $\psi$.

D.2 Squared Current-Current Commutator

In thr interaction $\mathcal{C}(t,\bm{x})$ is given by

[TABLE]

The zeroth order term can be deduced from the expression in D.6c

[TABLE]

Diagrammatically the first terms each look like a box but have momentum flowing in the opposite directions. The second two terms also have opposite momentum flow but have Wightman functions which are crossed. However, ladders with crossed legs should not contribute to chaos given the pole structure [64].

The first nontrivial in the expansion for the current-current commutator is

[TABLE]

Expanding the commutator on one of the time folds (again supressing spatial arguments) gives us

[TABLE]

The whole first order term will have $16$ terms. However, most contain a photon three point function and will vanish due to Furry’s theorem. The nonvanishing terms will gives us

[TABLE]

where the primes will technically be contracted with the unprimed terms. These correspond to boxes with one rung with momenta flowing in opposite directions. The parts in between the Wightman functions are the same as the first order correction to the fermion ladder. These are pictured in Fig. 7. To second order

[TABLE]

By inspection, the first term will give us terms with two photon rails or a self energy correction. The second term will give the second type of rung appearing in the fermion diagram. This give four terms where the momentum flow in the boxes at either end can go in the opposite direction. The third term can give us a ladder with self energy corrections on the rails.

Appendix E Details of numerical computations

Here we discuss the details of our numerical computation of $\Gamma_{{\bm{k}}}$ and $\lambda_{L}$. Our methods are similar to those used in [64], with some adjustments made to address various issues. The most challenging aspect of the numerics is managing the light-cone singularities that appear throughout the computation, which require additional care. To compute the components of the polarization bubble, we dealt with this issue by a slight off-shift in the grid. From our analysis in App. A, we know that all singularities that appear in the decay rate and the kernel should be integrable, so we can manually remove them during each integration.

The first step in our computation is to compute the analytic continuation of the $xx$ and $yy$ components of the polarization bubble. The imaginary parts are computed by numerical integration of Eqn. A.28 using a standard MATLAB routine. The integral is computed separately for $k^{0}<|{\bm{k}}|$ and $k^{0}>|{\bm{k}}|$ because of the non-analyticity at $k^{0}=|{\bm{k}}|$. The real parts are then obtained by a stable version of Kramers-Kronig given in the Eqn. A.41. These are then used to precompute the values of the spectral function, retarded propagator, and Wightman propagator on $200\times 200$ grids for $k^{0}/T\in[0,12.1]$, $|{\bm{k}}|/T\in[0,12.21]$. These grids are then used to create two dimensional cubic interpolations for all propagators. $\Gamma_{{\bm{k}}}$ is calculated from Eqn. B on a $60$ point grid for $|{\bm{k}}|/T\in[0,12]$ where $60$ grid points are used for the angular integral. To calculate the exponent on varying grid sizes, we use a one dimensional spline interpolation for the decay rate.

The kernel is computed using Eqn. C.2 on a $40\times 40$ grid for $|{\bm{k}}|/T,|{\bm{k}}^{\prime}|/T\in[0,12]$ where $40$ grid points are used for the angular integral. We find that most eigenvalues are complex with negative real parts, but the eigenvalue with the largest real part is completely real. This eigenvalue is identified as the Lyapunov exponent. There are also two eigenvalues which are exactly zero. We also compute the kernel using an interpolation on a grid that is slightly offset and find $\lambda_{L}\approx 4.06$.