Entanglement Dynamics after a Quench in Ising Field Theory: A Branch Point Twist Field Approach

TL;DR

This paper develops a perturbative approach using branch point twist fields to analyze the time evolution of entanglement entropy after a mass quench in the 1+1D Ising field theory, with results confirmed by lattice simulations.

Contribution

It introduces a novel perturbative method for calculating entanglement dynamics in massive quantum field theories post-quench, extending the twist field approach to time-dependent problems.

Findings

Linear growth of Rènyi entropies at large times

Presence of oscillatory contributions with frequency 2m

Good agreement with lattice numerical results

Abstract

We extend the branch point twist field approach for the calculation of entanglement entropies to time-dependent problems in 1+1-dimensional massive quantum field theories. We focus on the simplest example: a mass quench in the Ising field theory from initial mass $m_0$ to final mass $m$. The main analytical results are obtained from a perturbative expansion of the twist field one-point function in the post-quench quasi-particle basis. The expected linear growth of the R\'enyi entropies at large times $mt\gg 1$ emerges from a perturbative calculation at second order. We also show that the R\'enyi and von Neumann entropies, in infinite volume, contain subleading oscillatory contributions of frequency $2m$ and amplitude proportional to $(mt)^{-3/2}$. The oscillatory terms are correctly predicted by an alternative perturbation series, in the pre-quench quasi-particle basis, which we also discuss. A comparison to lattice numerical calculations carried out on an Ising chain in the scaling limit shows very good agreement with the quantum field theory predictions. We also find evidence of clustering of twist field correlators which implies that the entanglement entropies are proportional to the number of subsystem boundary points.

Tables

Table 1. Table 1: UV central charge and operator dimensions from the scaling limit extrapolations in the ground state with m = 0.04 𝑚 0.04 m=0.04 and ℓ = 128 ℓ 128 \ell=128 .

	$c$	${\Delta }_{\mu }$	${\Delta }_2$	${\Delta }_3$	${\Delta }_4$
Theory	0.5	0.125	0.0625	0.11111	0.15625
Fit	0.50195(3)	0.124969(4)	0.06258(1)	0.11002(9)	0.1525(3)

Equations

$${\rho }_A:={\mathrm{Tr}}_B(e^{-itH(\lambda )}|0⟩⟨0|e^{itH(\lambda )}),$$

$$S_n(t):=\frac{\log \mathrm{Tr}{\rho }_A^n}{1-n},$$

$$H_{\text{Ising}}(h)=-J\sum _{i=1}^N\left({\sigma }_i^x{\sigma }_{i+1}^x+h{\sigma }_i^z\right),$$

$$m:=2J|h-1|,v:=2Ja,$$

$$S_n(t)=\frac{\log \left({\epsilon }^{{\Delta }_n}{}_n⟨0|\mathcal{T}(0,t)|0{⟩}_n\right)}{1-n},$$

$${\tau }_n:={}_n⟨0|\mathcal{T}(0,0)|0{⟩}_n,$$

$$\underset{\ell \to \infty }{\lim }{}_n⟨0|\mathcal{T}(0,t){\mathcal{T}}^{†}(\ell ,t)|0{⟩}_n\sim {}_n⟨0|\mathcal{T}(0,t)|0{⟩}_n^2.$$

$$K(\theta )=i\tan \left[\frac{1}{2}{\tan }^{-1}(\sinh \theta )-\frac{1}{2}{\tan }^{-1}\left(\frac{m}{m_0}\sinh \theta \right)\right]:=i\hat{K}(\theta ),$$

$${}_n⟨0|\mathcal{T}(0,t)|0{⟩}_n={\tilde{\tau }}_n^{\prime }\exp \left[-\frac{n{\Gamma }^{\prime }mt}{2}-\frac{n{\mu }^2}{64\pi mt}-\frac{\mu }{8\sqrt{\pi }n}\frac{\cos \frac{\pi }{2n}}{{\sin }^2\frac{\pi }{2n}}\frac{\cos (2mt-\frac{\pi }{4})}{(mt{)}^{\frac{3}{2}}}+\cdots \right],$$

$$\mu :=1-\frac{m}{m_0}=-\frac{\delta m}{m_0},\mathrm{with}\delta m:=m-m_0.$$

$${\tilde{\tau }}_n^{\prime }={\tilde{\tau }}_n{e}^{A+O(K^3)}\mathrm{and}{\Gamma }^{\prime }=\Gamma +O(K^4),$$

$$\Gamma :=2{\int }_0^{\infty }\frac{d\theta }{\pi }{\hat{K}}^2(\theta )\sinh \theta ,$$

$$A:=\frac{1}{2\sin \frac{\pi }{n}}{\int }_{-\infty }^{\infty }\frac{d\theta }{2\pi }{\hat{K}}^2(\theta ).$$

$$S_n(t)=\frac{\log ({\epsilon }^{{\Delta }_n}{\tilde{\tau }}_n^{\prime })}{1-n}+\frac{\Gamma nmt}{2(n-1)}+\frac{n{\mu }^2}{64\pi mt(n-1)}+\frac{\mu }{8\sqrt{\pi }n}\frac{\cos \frac{\pi }{2n}}{{\sin }^2\frac{\pi }{2n}}\frac{\cos (2mt-\frac{\pi }{4})}{(n-1)(mt{)}^{\frac{3}{2}}}+O(t^{-3}).$$

$$\Gamma =\frac{{\mu }^2}{3\pi }+O({\mu }^3),A=\frac{{\mu }^2}{24\pi \sin \frac{\pi }{n}}+O({\mu }^3).$$

$$\underset{n\to 1}{\lim }\frac{\mu }{8\sqrt{\pi }n(n-1)}\frac{\cos \frac{\pi }{2n}}{{\sin }^2\frac{\pi }{2n}}\frac{\cos (2mt-\frac{\pi }{4})}{(mt{)}^{\frac{3}{2}}}=\frac{\mu }{4}\frac{\sqrt{\pi }\cos (2mt-\frac{\pi }{4})}{4(mt{)}^{3/2}}.$$

$$S_n(0)=\frac{\log \left({\epsilon }^{{\Delta }_n}{\tau }_n\right)}{1-n},$$

$${\Delta }_n=\frac{c}{12}\left(n-\frac{1}{n}\right),$$

$$S_n(0)=-\frac{c}{12}\left(n+\frac{1}{n}\right)\log (m_{0}a)+O(1),$$

$$|\Omega ⟩=\exp \left[{\int }_0^{\infty }\frac{d\theta }{2\pi }K(\theta )a^{†}(-\theta )a^{†}(\theta )\right]|\tilde{0}⟩.$$

$$|\Omega {⟩}_n=\exp \left[\sum _{j=1}^n{\int }_0^{\infty }\frac{d\theta }{2\pi }K(\theta )a_j^{†}(-\theta )a_j^{†}(\theta )\right]|\tilde{0}{⟩}_n,$$

$$S_n(t)=\frac{1}{1-n}\log \left(\frac{{\epsilon }^{{\Delta }_n}{}_n⟨\Omega |\mathcal{T}(0,t)|\Omega {⟩}_n}{{}_n⟨\Omega |\Omega {⟩}_n}\right).$$

$${}_n⟨\Omega |\mathcal{T}(0,t)|\Omega {⟩}_n:={\tilde{\tau }}_n\sum _{k_1,k_2=0}^{\infty }{C}_{2k_1,2k_2}(t),$$

$${\tilde{\tau }}_n{C}_{2k_1,2k_2}(t)=\frac{1}{k_1!k_2!}\sum _{j_1,\dots ,j_{k_1}=1}^n\sum _{p_1,\dots ,p_{k_2}=1}^n$$

$$\times \left[\prod _{s=1}^{k_1}{\int }_0^{\infty }\frac{d{\theta }_s^{\prime }}{2\pi }K({\theta }_s^{\prime }{)}^{\ast }{e}^{2itE({\theta }_s^{\prime })}\right]\left[\prod _{r=1}^{k_2}{\int }_0^{\infty }\frac{d{\theta }_r}{2\pi }K({\theta }_r)e^{-2itE({\theta }_r)}\right]$$

$$\times {}_{n;j_1{j}_1\dots j_{k_1}{j}_{k_1}}⟨{\theta }_1^{\prime },-{\theta }_1^{\prime },\dots ,{\theta }_{k_1}^{\prime },-{\theta }_{k_1}^{\prime }|\mathcal{T}(0,0)|-{\theta }_{k_2},{\theta }_{k_2},\dots ,-{\theta }_1,{\theta }_1{⟩}_{p_{k_2}{p}_{k_2}\dots p_1{p}_1;n},$$

$${}_n⟨\Omega |\Omega {⟩}_n:=\sum _{k=0}^{\infty }{Z}_{2k},$$

$$Z_{2k}=\frac{1}{(k!{)}^2}\sum _{j_1,\dots ,j_k=1}^n\sum _{p_1,\dots ,p_k=1}^n\left[\prod _{s=1}^k{\int }_0^{\infty }\frac{d{\theta }_s^{\prime }d{\theta }_s}{(2\pi {)}^2}K({\theta }_s^{\prime }{)}^{\ast }K({\theta }_s)\right]$$

$$\times {}_{j_1{j}_1\dots j_k{j}_k}⟨{\theta }_1^{\prime },-{\theta }_1^{\prime },\dots ,{\theta }_k^{\prime },-{\theta }_k^{\prime }|-{\theta }_k,{\theta }_k,\dots ,-{\theta }_1,{\theta }_1{⟩}_{p_k{p}_k\dots p_1{p}_1}\mathrm{for}k>0,$$

$$\frac{{}_n⟨\Omega |\mathcal{T}(0,t)|\Omega {⟩}_n}{{}_n⟨\Omega |\Omega {⟩}_n}:={\tilde{\tau }}_n\sum _{k_1,k_2=0}^{\infty }{D}_{2k_{1}2k_2}(t),$$

$$D_{2k_1,2k_2}(t)=\sum _{p=0}^{\min (k_1,k_2)}{\tilde{Z}}_{2p}{C}_{2(k_1-p),2(k_2-p)}(t),$$

$${}_{n;j_1}⟨{\theta }_1|\mathcal{T}(0,0)|{\theta }_2{⟩}_{j_2;n}$$

Full text

**Entanglement Dynamics after a Quench in Ising Field Theory:

A Branch Point Twist Field Approach**

Olalla A. Castro-Alvaredo*♡, Máté Lencsés♢, István M. Szécsényi♠, and Jacopo Viti♣*

${}^{\heartsuit\,\spadesuit}$ Department of Mathematics, City, University of London, 10 Northampton Square EC1V 0HB, UK

${}^{\diamondsuit\,\clubsuit}$International Institute of Physics, UFRN, Campos Universitário, Lagoa Nova 59078-970 Natal, Brazil

♣ECT, UFRN, Campos Universitário, Lagoa Nova 59078-970 Natal, Brazil

We extend the branch point twist field approach for the calculation of entanglement entropies to time-dependent problems in 1+1-dimensional massive quantum field theories. We focus on the simplest example: a mass quench in the Ising field theory from initial mass $m_{0}$ to final mass $m$. The main analytical results are obtained from a perturbative expansion of the twist field one-point function in the post-quench quasi-particle basis. The expected linear growth of the Rényi entropies at large times $mt\gg 1$ emerges from a perturbative calculation at second order. We also show that the Rényi and von Neumann entropies, in infinite volume, contain subleading oscillatory contributions of frequency $2m$ and amplitude proportional to $(mt)^{-3/2}$. The oscillatory terms are correctly predicted by an alternative perturbation series, in the pre-quench quasi-particle basis, which we also discuss. A comparison to lattice numerical calculations carried out on an Ising chain in the scaling limit shows very good agreement with the quantum field theory predictions. We also find evidence of clustering of twist field correlators which implies that the entanglement entropies are proportional to the number of subsystem boundary points.

Keywords: Entanglement Entropy, Quench Dynamics, Ising Model, Integrability, Branch Point Twist Fields, Perturbation Theory

♡ o.castro-alvaredo@city.ac.uk

♢ mate.lencses@gmail.com

♠ istvan.szecsenyi@city.ac.uk

♣ viti.jacopo@gmail.com

1 Introduction

Out-of-equilibrium many-body quantum dynamics is one of the most active and challenging research areas in low-dimensional physics; see the special issue [1] and in particular the review [2]. A typical setup triggering out-of-equilibrium evolution from an initial equilibrium state is that of a quantum quench [3]. In a quench protocol, a quantum system is prepared at $t<0$ in the ground state, denoted by $|0\rangle$, of a Hamiltonian $H(\lambda_{0})$ which depends on a parameter $\lambda_{0}$. At $t=0$, the parameter $\lambda_{0}$ is suddently changed to a new value $\lambda\not=\lambda_{0}$ and the unitary time evolution for positive times is governed by the new Hamiltonian $H(\lambda)$. The state of the system at time $t$ may be then formally written as $e^{-itH(\lambda)}|0\rangle$.

In this context, the evolution of the bipartite, or von Neumann, entanglement entropy following a quantum quench has been intensively studied; see [4] for a review and references therein. Consider a space bipartition of a 1+1-dimensional quantum system as sketched in fig. 1 and assume that regions $A$ and $B$ are semi-infinite. Then the entanglement entropy associated to region $A$ after a quench may be expressed as $S(t)=-{\rm Tr}_{A}(\rho_{A}\log\rho_{A})$ where formally

[TABLE]

is the reduced density matrix associated to subsystem $A$. Since the regions are semi-infinite, the entropies will not explicitly depend on the subsystem’s length. Another set of entanglement measures is provided by the Rényi entropies which are defined as

[TABLE]

and have the property $\lim_{n\rightarrow 1}S_{n}(t)=S(t)$. It is in fact these Rényi entropies which we will mostly be studying in this manuscript.

The universal features of the evolution of entanglement after a quench at a critical point described by conformal field theory have been studied in [5, 3, 6]. In these works an intuitive picture was put forward, namely one based on the production of highly entangled quasi-particle pairs of opposite momenta right after the quench. These then propagate in space-time until a critical time $t_{\text{sat}}=\frac{\ell}{2v}$, where $\ell$ is the size of the subsystem and $v$ is the propagation velocity. In this region, the entanglement entropy grows linearly in time. For $t_{\text{sat}}>\frac{\ell}{2v}$ the entanglement saturates to a value proportional to the subsystem’s size $\ell$. These features were later demonstrated analytically for the XY chain in a transverse magnetic field in [7], where the exact coefficients of the terms linear in $t$ and in $\ell$ were computed. Note however that for our configuration of two semi-infinite regions $t_{\text{sat}}\rightarrow\infty$. Thus we expect the entanglement to continue to grow linearly in time for all times; this must be beared in mind when comparing analytic results with lattice numerical calculations.

Although the entanglement evolution after a quench has been studied for many physical lattice models, both free [7, 8, 9, 10, 5, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] and interacting [21, 22, 23, 24], complete analytic derivations of such linear behaviour remained elusive so-far. In $1+1$ dimensions, they are still based either on conformal field theory techniques [5, 25] or conjectural large space-time asymptotics for block Töplitz matrices [7, 8, 9]. An exception is represented by models with random unitary evolution such as random circuits [26, 27].

More recently, it has been shown that for gapped systems, the entanglement dynamics after a quench features other non-trivial effects. In particular the time-dependence of the entanglement entropy can show subleading corrections that might qualitatively alter the leading linear increase at large times. For instance, the studies [28, 29, 30] have shown that quasi-particle confinement in a linear potential can lead to oscillatory behaviour in time, as well as suppression of linear growth for sufficiently large times.

The purpose of this paper is then twofold: first we will provide a general quantum field theory framework to analyse entanglement dynamics in massive systems. Secondly, we will provide evidence that subleading oscillatory terms are actually a common feature of entanglement dynamics in infinite volume. To this end, we will focus on one of the simplest and best known theories: the Ising field theory. We may regard this as the scaling limit of the Ising chain described by the Hamiltonian

[TABLE]

where $J>0$, and $h$ is known as the transverse field. The Ising spin chain has a quantum critical point, with a gapless spectrum, at $h=1$, which separates a paramagnetic phase $(h>1)$ from a ferromagnetic phase ($h<1$). The two phases are related by a Kramers-Wannier duality transformation, which interchanges the spin with the disorder field.

Near the critical point, for $|h-1|\ll 1$ it is possible to define the scaling limit by taking $J\rightarrow\infty$ and $a\rightarrow 0$, where $a$ is the lattice spacing, while keeping

[TABLE]

fixed and finite [31]. In the scaling limit, the low energy excitations of (3) are then relativistic real non-interacting fermions with positive mass $m$, while the speed of light is fixed to $v$.

Finally, a word is due on the techniques that we will be using in this paper. We will exploit the well-known relationship between Rényi entropies and correlation functions of branch point twist fields [32, 33, 34]. For the simple configuration of fig. 1 this means that we will be computing a one-point function of a branch point twist field $\mathcal{T}(x,t)$ and studying its time dependence after the quench. Branch point twist fields are defined on a replicated quantum field theory containing $n$ identical copies of the original theory, and have been interpreted as symmetry fields associated to the cyclic permutation of the copies in [34]. Explicitly, the Rényi entropies at time $t$ are given by

[TABLE]

where $\varepsilon$ is some UV cut-off; ${\cal T}(0,t):=e^{iH(\lambda)t}{\cal T}(0,0)e^{-iH(\lambda)t}$ is the time-evolved twist field in the Heisenberg picture; $n$ is the replica number; $|0\rangle_{n}$ is the ground state of the replica theory before the quench that is, for coupling constant $\lambda_{0}$ (corresponding to mass gap $m_{0}$); $\Delta_{n}$ is the scaling dimension of the twist field at criticality. For $t=0$, (5) gives the Rényi entropies at equilibrium in terms of

[TABLE]

the Vacuum Expectation Value (VEV) of the branch point twist field, which by dimensional analysis must be proportional to $m_{0}^{\Delta_{n}}$. In particular, the UV cut-off $\varepsilon$ is chosen in such a way that no finite $O(1)$ term appears at $t=0$ on the right hand side of (5).

When comparing our field theoretical predictions for (5) with lattice calculations in the Ising chain in the scaling limit, we will be actually comparing (5) with a similar quantity involving a two-point function of branch point twist fields ${}_{n}\langle 0|{\cal T}(0,t){\cal T}^{\dagger}(\ell,t)|0\rangle_{n}$: that is the entanglement entropy of an interval of length $\ell$. To allow for comparison, we will take the length of such interval to be very large, in which case we expect clustering of the two-point function to occur, namely, the factorization

[TABLE]

Thus our results for (5) will generally give half the values obtained from computations involving a large but finite interval $\ell$. We indeed confirm this in section 6 of this paper.

This paper is organized as follows: in section 2 we present a summary of our analytical results. In section 3 we review the field theoretical tools that we have used: a time-dependent formulation of the branch-point twist field approach [34] for the calculation of entanglement measures. In particular, we present an expansion of the twist field one-point function in the post-quench quasi-particle basis following the route traced in [31] for the order parameter. In section 4, we derive the main analytical results. In section 5 we further generalize the perturbative approach to the quench dynamics put forward in [35] to the calculation of the twist field one-point function. We show that for sufficiently small quenches, these results are in agreement with the main outcome of section 4. In section 6 we present a detailed test of our field theoretical predictions against lattice results obtained in the scaling limit. Finally, we conclude in section 7. An appendix with additional numerical lattice results completes the paper.

2 Summary of the Main Results

Consider the Ising field theory with mass scale $m_{0}$ and a quench that changes it to a new value $m$. Let us also introduce the function ($\theta\in\mathbb{R}$)

[TABLE]

whose meaning we discuss in section 3. Long time after the quench, namely for $mt\gg 1$ the expectation value of the branch point twist field is conjectured to be

[TABLE]

where

[TABLE]

The ellipsis in (9) denote terms that are subleading with respect to $t^{-3/2}$ for large times. The parameters $\tilde{\tau}_{n}^{\prime}$ and $\Gamma^{\prime}$ in (9) are calculated perturbatively in the function $K(\theta)$ whose absolute value is then assumed small for $\theta$ real. In particular

[TABLE]

where $\tilde{\tau}_{n}$ is the expectation value of the branch point twist field in the post-quench ground state $|\tilde{0}\rangle_{n}$, similar to the definition (6) but with mass gap $m$. The decay rate $\Gamma$ and the constant $A$ in (11) are

[TABLE]

and

[TABLE]

We will provide a full derivation of the twist field one-point function up to $O(K^{2})$ and conjecture its general form in (9) following an analogous calculation as for the one-point function of the spin operator after a mass quench [31] in Ising field theory. In particular, we will show that the decay rate $\Gamma$ in (12) is the same as for the spin operator [31] up to the second-order corrections in the function $K$. The oscillatory contribution in (9) has also the same frequency and power law in $mt$ as for the spin operator, albeit with a different $n$-dependent overall coefficient.

From the explicit calculation of the one-point function of the branch-point twist field, we can derive an exact expression at $O(K^{2})$ for the Rényi entropies (5) at large times after the quench which is given by

[TABLE]

Since for $|\mu|\ll 1$, $\hat{K}(\theta)$ is $O(\mu)$ and

[TABLE]

Eq. (14) can also be viewed as a large-time expansion of a perturbative series in the quench parameter $\mu$. In particular, due to (15), eq. (14) is exact up to second-order terms in $\mu$. This claim will be checked explicitly against analytical and numerical lattice results for the Rényi entropies in the scaling limit in section 6.

Furthermore, notice that the oscillating term is $O(\mu)$. Indeed, under the assumption $|\mu|\ll 1$ and provided the replacement $m\rightarrow m_{0}$ in the frequency, it can be also derived within a first-order perturbative approach to the quench dynamics [35, 36]. We will postpone details of this alternative derivation to section 5. Remarkably, this first-order perturbative approach is not enough to capture the leading large-time asymptotics of the Rényi entropies. The latter is governed by the decay rate $\Gamma$ and is therefore an $O(\mu^{2})$ effect.

Another interesting feature of (14) is that the limit $n\rightarrow 1$ is only well-defined for the oscillatory term; for the von Neumann entropy we obtain in particular

[TABLE]

The same limit is obviously ill-defined for all the other contributions in (14). The reason for this has to do with the way the $O(\mu^{2})$ terms are computed, namely from branch point twist fields which are a priori only defined for $n\ \in\mathbb{N}\setminus\{0,1\}$. Taking the limit generally requires an understanding of the analytic continuation to $n\in\mathbb{R}$ of the one-point function, which is non-trivial, particularly for higher particle form factor contributions (e.g. precisely the ones that give rise to the problematic terms in (14)). We have not carried out this limit here, but we have found good agreement between (16) and lattice calculation in the scaling limit for the Neumann entropy. We report these comparisons in section 6

3 Review of the Main Techniques

In this section we review the main techniques that we have employed in order to derive the results of the previous section: branch point twist fields in relation to entanglement measures in the Ising field theory and the expansion of the one-point function of a local operator in the post-quench quasi-particle basis, as developed in [31]. Results obtained from a perturbative approach [35, 36] in the quench parameter will be given in section 5.

3.1 Branch Point Twist Fields and Entanglement

The main properties of branch point twist fields were described at length in [34] and the subsequent review article [37]. In quantum field theory, it has been known for some time [32, 33, 34] that, for integer $n$, the Rényi entropies in (2) may be expressed in terms of correlation functions of branch point twist fields, with the number of twist field insertions equalling the number of boundary points of the subsystems under consideration. This means that the entanglement entropy of a semi-infinite region is simply given by the one-point function of the branch point twist field. Eq. (5) can be thought as the obvious time-dependent generalization of the setting in [34] at equilibrium.

Before embarking into the study of the time-dependent one-point function of the twist field, it is useful to recall some of its properties at equilibrium. At equilibrium, the $n^{\mathrm{th}}$ Rényi entropy of a semi-infinite system in the ground state of $H(\lambda_{0})$ is given by (5) at $t=0$, namely

[TABLE]

Obviously in the ground state $|\tilde{0}\rangle_{n}$ of the post-quench Hamiltonian $H({\lambda})$ (with a mass gap $m$) (17) applies by replacing $\tau_{n}\rightarrow\tilde{\tau}_{n}$. The power $\Delta_{n}$ is the scaling dimension of the branch point twist field at criticality, which is given by

[TABLE]

where $c$ is the central charge of the underlying CFT [38, 39, 32]; for instance, $c=\frac{1}{2}$ for Ising field theory. The parameter $n$, which was already introduced in section 1, is the number of copies of the replicated Hilbert space of the quantum field theory, upon which the branch point twist field acts. Therefore, in such a replicated theory, the pre- and post-quench ground states $|0\rangle_{n}$ and $|\tilde{0}\rangle_{n}$ are tensor products of $n$ copies of the physical ground states defined in the introduction. The same construction carries over for the time-dependent case.

Finally, it is useful to remember that when comparing with lattice results for the Ising spin chain, the natural choice for the cut-off is the lattice spacing $a$. The UV cut-off $\varepsilon$ and $a$ are related by a model-dependent (i.e. non-universal) proportionality constant. Therefore on the lattice (17) reads [32]

[TABLE]

where $O(1)$ denotes non-universal terms that are finite or vanish in the scaling limit $a\rightarrow 0$. The leading logarithmic lattice spacing dependence in (19) can be used to extract the twist field scaling dimension, alias the central charge of the UV fixed point, in lattice numerical calculations. Actually the quality of such an extrapolation provides a useful measure of how close the numerical calculation is to the scaling regime of the lattice model; see section 6 and appendix A.

3.2 Expansion of the Time-Dependent One-Point Function in the Post-Quench Basis

In this section we review the approach first employed in [31] to study relaxation dynamics of a local operator in the Ising field theory after a mass quench. The technique needs two inputs: an expansion of the initial state into eigenstates of the post-quench Hamiltonian and the matrix elements of the local operator one is interested in, between states of the post-quench quasi-particle basis. In principle the method is applicable also to interacting post-quench theories, provided that such analytical data are known; in particular the post-quench theory, considered in infinite volume and for all times, should be integrable. See for instance [40] for recent activity devoted to overlap calculations.

In the specific case of the Ising mass quench, the non-normalized initial state $|\Omega\rangle:=\sqrt{\langle\Omega|\Omega\rangle}|0\rangle$, $|0\rangle$ being the ground state of the pre-quench Hamiltonian, can be exactly expressed in terms of eigenstates of the post-quench Hamiltonian as

[TABLE]

Notice therefore that $|\tilde{0}\rangle$ is the vacuum of the post quench Ising field theory (i.e. with mass gap $m$); $a^{\dagger}(\theta)$ is the fermionic creation operator, and $K(\theta)$ is the function given earlier in (8). The integral in (20) is over the so-called rapidity which parametrizes the energy ($E$) and momentum ($P$) of the one-particle state $|\theta\rangle:=a^{\dagger}(\theta)|\tilde{0}\rangle$ as follows: $E=m\cosh\theta$ and $P=m\sinh\theta$. The normalization of the one-particle states is $\langle\theta|\theta^{\prime}\rangle=2\pi\delta(\theta-\theta^{\prime})$.

Quenches leading to states with the structure (20) where studied in detail in [41, 42]. A derivation of the function (8) is given in Appendix A of [42]. These states have the same structure of the boundary states first described by Ghoshal and Zamolodchikov [43]. Their structure neatly fits with the quasi-particle picture put forward in [5, 3, 6] as the initial state (20) can be regarded as a coherent superposition of particle pairs, also known as a squeezed coherent state. Exact solvability of the quench dynamics, which is generally not possible, has been also related [44] to initial states analogous to (20), see for instance [35].

In the $n$-copy theory, this simply generalises to

[TABLE]

where $a_{j}^{\dagger}(\theta)$ is the fermionic creation operator in copy $j$. We denote by $|\theta_{1},\ldots,\theta_{k}\rangle_{j_{1},\ldots,j_{k};n}$ an element of an orthonormal basis in the replicated ($in$ or $out$) Hilbert space consisting of $k$ particles with rapidities $\theta_{i}$ and copy labels $j_{i}$, $i=1,\ldots,k$. The energy and momentum of multi-particle states are the sum of the energies and momenta of their one-particle constituents.

In such a framework, the Rényi entropies after the quench can be written as

[TABLE]

Substituting the representation (21) of the replicated initial state into (22), both numerator and denominator admit a formal expansion as sums of integrals of matrix elements in the post quench basis. Borrowing notations from [31], we will write these series as

[TABLE]

with

[TABLE]

and analogously

[TABLE]

where now

[TABLE]

and $Z_{0}=1$. The ratio in (22) can be then expanded formally in powers of the function $K$

[TABLE]

with

[TABLE]

where $\tilde{Z}_{2p}$ are the expansion coefficients of the inverse of the norm, i.e. $\sum_{k,p=0}^{\infty}Z_{2k}\tilde{Z}_{2p}=1$. In section 4 we present the calculation up to $O(K^{2})$.

The matrix elements of the twist field in (24) can be related to the so-called elementary form factors [45, 46, 34], see (35) in the next section. The transformation that relates the two functions is called crossing. Consider for instance the matrix element ${}_{n;j_{1}}\langle\theta_{1}|\mathcal{T}(0,0)|\theta_{2}\rangle_{j_{2};n}$. This can be written as

[TABLE]

where $\theta_{12}:=\theta_{1}-\theta_{2}$, $\eta$ is a small positive parameter and $F_{2}^{j_{1}j_{2}}(\theta)$ will be given in (33). This relation can be generalized to matrix elements involving states with larger number of particles [46]. The shift by $i\eta$ makes the function $F_{2}^{j_{1}j_{2}}(\theta)$ on the right hand side of (29) regular for $\theta\rightarrow i\pi$. There are however additional sources of divergences related to the normalization of the asymptotic states in infinite volume, see the $\delta$ function in (29).

These infinite volume singularities are expected to be cancelled by similar singularities in the denominator in (26) in the combination as (28). The precise way in which this cancellation occurs has been the object of much investigation over the past decade and a rigorous understanding now exists. That is, to consider the theory in finite volume $V$ and use the volume as regulator [47, 48]. However, this rigorous approach is rather involved and for this reason some simpler methods, have also been developed. In [31] a regularization scheme known as $\kappa$-regularization [49, 50] was used. The technique requires to shift the coinciding rapidities by a real value $\kappa$ (or several values $\kappa_{i}$ for multi-particle states) so that the singularities are avoided. Then introduce a smooth function $P(\kappa)$ which is strongly peaked about $\kappa=0$ with the properties

[TABLE]

Of course, there are many functions that would meet the criteria above but one expects that in the infinite volume limit $V\rightarrow\infty$ they will all lead to the same finite result. A natural choice also employed in [31] is a gaussian $P(\kappa)=Ve^{-\pi\kappa^{2}V^{2}}$. For instance for a two-particle form factor the regularization would be implemented as

[TABLE]

After using the crossing relation (29), it is possible to isolate the infinite volume divergences, coming from the normalization of the states, and the leading contribution for $V\rightarrow\infty$, by expanding the integrand as a series about $\kappa=0$. Applications can be found in [31] and in section 4.

3.3 Form Factors of the Branch Point Twist Field in the Ising Field Theory

In this section, we will finally recall the necessary results for the form factors of the twist field in the Ising field theory. The explicit form of the form factors is needed to evaluate the numerator of (22). In the replicated Ising field theory fermionic particles have an extra copy index $j=1,\dots,n$. There is also an internal $\mathbb{Z}_{2}$ symmetry in each copy which implies, for $\mathbb{Z}_{2}$ even fields, such as the twist field, that only even-particle form factors are non-vanishing. As already mentioned in the previous section, let $|\theta_{1},\ldots,\theta_{k}\rangle_{j_{1},\ldots,j_{k};n}$ be an asymptotic in state of the replicated theory consisting of $k$ particles with rapidities $\theta_{i}$ and copy labels $j_{i}$, $i=1,\ldots,k$. We further assume $\theta_{1}>\theta_{2}>\cdots>\theta_{k}$. The two-particle twist field form factor is defined as [34]

[TABLE]

and is given by

[TABLE]

with

[TABLE]

and $\tau_{n}$ was defined in (6). Here we have used the fact that the twist field is a Lorentz scalar and therefore the form factor depends only on the rapidity difference, rather than two separate rapidities. Due to the free nature of the theory the $k$-particle form factors are given in terms of Pfaffians [51],

[TABLE]

where $\mathrm{Pf}$ is the Pfaffian of the matrix $W$ (i.e. $\mathrm{Pf}^{2}(W)=\det(W)$), which in turn is defined as

[TABLE]

In practice, (35) implies that, if we call the two-particle form factor in (36) a contraction, $k$-particle form factors ($k$ even) are obtained as sums of products of contractions as prescribed by the Wick theorem for fermionic fields. Due to the particular monodromy properties of the branch point twist field discussed in [34], all form factors can be ultimately expressed in terms of form factors involving only one copy of the theory. In particular

[TABLE]

for $j_{1}\geq j_{2}\geq\cdots\geq j_{k}$. For the Ising field theory this means that the two-particle form factor of particles in the first copy is effectively the building block for any other form factor. For this reason it is useful to adopt a simpler notation for this form factor. We then define the normalized two-particle form factor

[TABLE]

All formulas presented in this section are valid when considering the post-quench ground state $|\tilde{0}\rangle_{n}$, if one replaces $\tau_{n}$ with $\tilde{\tau}_{n}$.

4 Rényi Entropies after a Mass Quench: Field Theory Results

As outlined at the beginning of section 3.2, the calculation is organized as a perturbation series in powers of the function $K$ introduced in (8). In principle, the final result is not limited to $\delta m\ll 1$, see (10), provided $K(\theta)$ is sufficiently small for $\theta\in\mathbb{R}$. Physically this is equivalent to truncating the series in (20) to a few multi-particle states. In this section we will fill in the details of the derivation of (14).

4.1 Contributions at $O(K)$

Apart from a trivial $K$-independent term, corresponding to $C_{00}=D_{00}=1$ in (23), the leading term in the $K$ expansion of (22) is $O(K)$ and given by

[TABLE]

where we have used (38) and (33). Notice that the expansion of the denominator in (25) starts as $1+O(K^{2})$, therefore, see (28), $C_{2,0}+C_{0,2}=D_{2,0}+D_{0,2}$. At large times, according to stationary phase analysis, we can expand the integrand in (39) close to $\theta=0$ and observe that

[TABLE]

with $\mu$ defined by (10). By retaining only contributions up to $O(K)$, the one-point function of the twist field is then for $mt\gg 1$

[TABLE]

As anticipated in section 1 for $|\mu|\ll 1$ the same result can be derived from a perturbation theory approach [35]; see section 5 for details. We will show in a subsequent section that the terms above are in fact just the first two contributions to the expansion of an exponential, hence the expression (9).

4.2 Contributions at $O(K^{2})$

The $O(K^{2})$ contributions are considerably more involved and provide a first indication that Rényi entropies after the quench grow linearly in time. Taking into account numerator and denominator in (22), the $O(K^{2})$ contributions in the expansion of the one-point function are given by, see again (28)

[TABLE]

and

[TABLE]

4.2.1 The Contribution $D_{2,2}$

We start analysig $D_{2,2}$ in (42); from (24) one has

[TABLE]

The second equality follows from permutation symmetry of the replicas, and we also defined

[TABLE]

Finally, $C_{0,0}=1$ and from (25) it follows

[TABLE]

To further manipulate (44), we exploit the crossing relation [46]

[TABLE]

which generalises (29) to four-particle states. In (47) and hereafter, we used the notation: $\theta_{\pm}:=\theta\pm i\pi\,$. After substituting (47) into (44) we regroup the result into three terms: $C_{2,2}:=C_{2,2}^{(0)}+C_{2,2}^{(2)}+C_{2,2}^{(4)}$. The integrand of $C_{2,2}^{(0)}$ contains the first line in (47), the integrand of $C_{2,2}^{(2)}$ contains the second and third line in (47) while the integrand of $C_{2,2}^{(4)}$ includes the four-particle form factor in the last line of (47). Now it is easy to see that $C_{2,2}^{(0)}=Z_{2}C_{00}$ and therefore the only non-vanishing contribution at $O(K^{2})$ is, see (42),

[TABLE]

The double integral $C_{2,2}^{(2)}$, after integrating the delta function over $\theta^{\prime}$ and exploiting the symmetries $M(-x,y;t)=-M(x,y;t)$ and $M(-x,-y;t)=M(x,y;t)$, can be rewritten as

[TABLE]

Notice that the sum over $j$ in (44) reduces in this case to only one term, due to the Kronecker delta in (47). In the $\kappa$ regularization scheme, eq. (49) should be first integrated with the measure $P(\kappa)$, discussed in section 3, and then the outcome of the integration evaluated in the limit $V\rightarrow\infty$ and $\eta\rightarrow 0$. In practice, one expands (49) in a power series in $\kappa$ close to $\kappa=-i\eta/2$ and observes that $\int d\kappa P(\kappa)\kappa^{n}=O(V^{-n})$, therefore in the infinite volume limit only terms that are singular or finite for $\kappa,\eta\rightarrow 0$ contribute to the final result. Actually, when summed up at a given order in $K$, divergent terms in $\kappa$ should cancel consistently. Expanding the function $C_{2,2}^{(2)}$ around $\kappa=-i\eta/2$ we obtain

[TABLE]

The third and fourth terms vanish by symmetry, while the first one which is divergent in the limit $\kappa,\eta\rightarrow 0$, will be cancelled by an opposite contribution coming from $C_{2,2}^{(4)}$. In conclusion only the second time-independent term in (50), contributes to the final result for the twist field one-point function. Such a constant was called $A$ in (13).

Let us then finally analyze $C_{2,2}^{(4)}$. This is a double integral weighted by the function $M(\theta^{\prime},\theta;t)$ of the four-particle form factor in (47). For the Ising model such a form factor is obtained, see the definition (35), applying the Wick theorem as

[TABLE]

Using (37) and in particular $F_{2}^{1j}(\theta)=\tilde{\tau}_{n}f(2\pi i(j-1)-\theta)=-\tilde{\tau}_{n}f(\theta-2\pi i(j-1))$ for $j\neq 1$, we can rewrite the two-particle form factors in (51) in terms of the elementary function $f$, see section 3. The sums over $j$, needed to construct $C_{2,2}^{4}$, see (44), can be then performed by using the identity

[TABLE]

that can be found for instance in the Appendix of [51]. This gives

[TABLE]

The two lines in (53) have to be finally integrated over the rapidities $\theta$ and $\theta^{\prime}$ to obtain the function $C_{2,2}^{4}(t)$. We define $I(t)$ to be the result of integrating the second line in (53) (i.e. the function inside the square bracket) and $I^{\prime}(t)$ to be the result of integrating the first line (i.e. the product of functions $f$). In this way $C_{2,2}^{4}(t)=I(t)+I^{\prime}(t)$; we start by analyzing $I^{\prime}(t)$ which is simply

[TABLE]

$C_{2,0}$ given in (39); the result follows from $f(\theta)=-f(\theta)^{*}$ for $\theta\in\mathbb{R}$. Notice that, according to the discussion in (16), for large times $I^{\prime}(t)=O(t^{-3})$.

The remaining integral to complete our calculation of $D_{2,2}(t)$ is $I(t)$. After substituting the explicit form for the function $G$, given in (52) into (53) and using $M(x,y;t)=-M(x,-y;t)$ it can be eventually rewritten as

[TABLE]

where we have introduced the function

[TABLE]

which is regular along the integration contour in the variable $\theta$ in the limit $\eta_{1,2}\to 0$; also $\eta_{12}:=\eta_{1}-\eta_{2}$.

The denominator in (55) has poles at $\theta=\theta^{\prime}-\kappa-i\eta_{1}$ and $\theta=\theta^{\prime}+\kappa+i\eta_{2}$. To calculate the $\kappa$-regularized part of the integral and evaluate the $\eta_{1,2}\to 0$ limit, we modify the integration contour for $\theta$ to be the sum of the contours

[TABLE]

where $s$ and $\phi$ are parameters chosen carefully. We have that $s<\theta^{\prime}-\kappa$, $\eta_{1}<\phi$, and $\phi$ has to be smaller than the position of the branch point in the function $\hat{K}$. When shifting the contour from the real axis to $\mathcal{C}_{3}$ we encounter a pole at $\theta=\theta^{\prime}-\kappa-i\eta_{1}$ and pick up the residue contribution, in a clockwise direction, with the value

[TABLE]

where $\eta:=\eta_{1}+\eta_{2}$. Expanding the integrand in (58) around $\kappa=-i\eta/2$, sending $\eta\rightarrow 0$, and calling $\theta$ the integration variable, we have

[TABLE]

The first term in (59) exactly cancels the two-particle form factor singularity, i.e. the first term on the right hand side of (50). The second term is remarkably linear in time with coefficient $-\frac{nm\Gamma}{2}$ and $\Gamma$ given in (12). The third term vanishes due to $\hat{K}^{2}(0)=\hat{K}^{2}(\infty)=0$. We can then finally write

[TABLE]

where the last integral is now well defined as there are no singularities along the contour of integration of the rapidity $\theta$.

It is also possible to extract the large time limit of the integral

[TABLE]

which appears in (60). The integrand of (61) has a double pole on the real axis of the variable $\theta$, however this can be cured, without spoiling convergence at infinity, by taking a double derivative with respect to time. After taking the double derivative the integration contour for $\theta$ can be lifted back to the real axis. We are then led to consider the large $t$ asymptotics of the following double integral

[TABLE]

This can be done by standard application of the stationary phase approximation for two-dimensional integrals. There is only one stationary point at $\theta=\theta^{\prime}=0$; Taylor-expanding the integrand about $\theta=\theta^{\prime}=0$ gives at leading order

[TABLE]

Integrating back twice we obtain the desired asymptotic for the integral $R(t)$ in (60) which is

[TABLE]

Notice that when integrating back, we are setting to zero possible terms $O(t)$ and $O(1)$, due to the asymptotic of the original integral. In summary, we have shown that

[TABLE]

We will revisit this result in subsection 4.3 where we argue that these contributions are nothing but the first non-trivial term in the expansion of the exponential featuring in (9).

4.2.2 The Contributions $D_{0,4}$ and $D_{4,0}$

Finally we analyze the contributions $D_{0,4}$ and $D_{4,0}$. Since $D_{4,0}=D_{0,4}^{*}$, we focus only on $D_{0,4}$, which is given by $D_{0,4}=C_{0,4}$ with

[TABLE]

and we defined, analogously to (45)

[TABLE]

The four-particle form factor in (66) can be decomposed as in (51) by applying Wick’s theorem and the sum over the index $j$ performed by recalling (53). By repeating steps similar to those employed in the section 4.2, the integral $C_{0,4}$ can be written as a sum of two terms, namely $C_{0,4}(t):=C_{0,4}^{(1)}(t)+C_{0,4}^{(2)}(t)$. In particular, one obtains

[TABLE]

and

[TABLE]

where $H(x,y)$ is the same function given in (56). By applying the stationary phase approximation we can estimate the large time limit of  (69). This gives another $O(t^{-3})$ contribution (since (68) is also of $O(t^{-3})$), namely

[TABLE]

with

[TABLE]

This closes our calculation of the branch point twist field one-point function at $O(K^{2})$.

4.2.3 The Complete Formula at $O(K^{2})$ for the Twist Field One-point Function

We can finally summarize the result for the twist field one-point function up to $O(K^{2})$

[TABLE]

where $A$ is given in (13), $C_{2,0}(t)+C_{0,2}(t)$ in (39), $\Gamma$ in (12), $R(t)$ in (61) and $C_{0,4}^{(2)}(t)=[C_{4,0}^{(2)}(t)]^{*}$ in (69). Note also that

[TABLE]

with $I^{\prime}(t)$ given in (54) and the other terms in (68). The presence of the contributions $1+C_{2,0}(t)+C_{0,2}(t)+\frac{(C_{2,0}(t)+C_{0,2}(t))^{2}}{2}$ suggests that these and higher-order terms may arise from the exponentiation of $C_{2,0}(t)+C_{0,2}(t)$. Indeed, it is possible to argue that all terms in the expansion (72) exponentiate. We will give a simple argument towards this conclusion in subsection 4.3.

Reinserting the cut-off dependence $\varepsilon^{\Delta_{n}}$ in (72), taking the logarithm, expanding its argument up to $O(K^{2})$ and dividing by $1-n$ we obtain the following expression up to $O(K^{2})$ for the Rényi entropies

[TABLE]

where $\tilde{S}_{0,n}$ is the Rényi entropy in the ground state of the post-quench Hamiltonian i.e.

[TABLE]

Expanding (74) for $mt\gg 1$, we obtain the result quoted in (14).

4.3 An Argument Towards Exponentiation at Higher Orders

A simple combinatorial argument can be provided to show that all the terms in the expansion (72) exponentiate. In other words, they result from the expansion of an exponential at order $O(K^{2})$. The exponent will receive $O(K^{3})$ and higher corrections which we will not investigate in this paper. Note that in [31] an argument was given for the exponentiation of the term $-\Gamma mt$ (the equivalent of our $-\frac{\Gamma nmt}{2}$ term but for the order parameter). An entirely similar argument can be given for the branch point twist field to show the exponentiation of this term. However, we find that exponentiation is a much more general feature of the one-point function, extending to other terms at $O(K^{2})$ as well. As we will see, our calculation does not use any special properties of the branch point twist field form factors, apart from their Pfaffian structure. Therefore we expect the same exponentiation to occur for the order parameter $\frac{\langle\Omega|\sigma(0,t)|\Omega\rangle}{\langle\Omega|\Omega\rangle}$.

Examining the generalization of the crossing relation (47) and the Wick contraction nature of the form factor expressions (35) and (36), it is natural to expand the $C_{2k,2l}(t)$ functions as sums of products of connected contributions, that is

[TABLE]

where $C_{2i,2j}^{c}(t)$ are related to integrals of “connected” matrix elements, which are defined recursively from the condition of not being factorizable into other connected expressions. The $n_{i,j}$ are non-negative integers that satisfy the constraints $\sum_{i,j=0}^{\infty}i\,n_{i,j}=k$ and $\sum_{i,j=0}^{\infty}j\,n_{i,j}=l$. By inverting the expansion (76), for the first few connected terms we get for instance

[TABLE]

These new combinations of terms are immediately recognizable from our earlier computation. For instance, $C_{0,4}^{c}(t)$ is nothing but $C_{0,4}^{(2)}(t)$ defined in (69), and $C_{2,2}^{c}(t)$ is obtained from $C_{2,2}(t)$ after subtracting the term $I^{\prime}(t)$ defined in (54). The norm of the initial state $|\Omega\rangle_{n}$ admits an analogous expansion

[TABLE]

where $\sum_{i=0}^{\infty}i\,\tilde{n}_{i}=k$. To calculate the regular terms of the one-point function $D_{2k,2l}(t)$, see eq. (28), we need to calculate the inverse of the norm defined by the condition $\sum_{km,p=0}^{\infty}Z_{2k}\tilde{Z}_{2p}=1$. From observation of the first few terms of the inverse of the norm, we expect its connected expansion to have the form

[TABLE]

In the following, we are only focusing on terms of the one-point function, that contain connected matrix elements of at most $O(K^{2})$, i.e. we consider only terms where powers $n_{i,j}=0$ for $i+j>2$ and $\tilde{n}_{i}=0$ for $i>1$. With this assumption (76) takes the form

[TABLE]

where

[TABLE]

with $\lfloor.\rfloor$ denoting the integer part. Plugging these formulas and (81) into $D_{2(k+l),2k}$ using the form (28), and exchanging the order of the summations leads to

[TABLE]

where the combination $D_{2,2}^{c}(t)=C_{2,2}^{c}(t)-Z_{2}^{c}$ is both connected and regular. Similar results hold for $D_{2k,2(k+l)}(t)$ and $D_{2k,2k}(t)$, hence the one-point function (27) takes the form

[TABLE]

Further manipulation of the order and range of the summations, allows us to write the one-point function as

[TABLE]

where $D_{2i,0}^{c}(t)=C_{2i,0}^{c}(t)$, $D_{0,2i}^{c}(t)=C_{0,2i}^{c}(t)$, since these terms are regular without any subtraction. Note that the terms in the exponent are precisely those inside the bracket in (74). With this we showed, using the assumption (81), that the one-point function exponentiates up to $O(K^{2})$ terms. The expression is regular, since all the singularities are cancelled as explained in detail in the previous sections.

Given the simplicity of the Ising form factors, we expect exponentiation to occur also at higher orders in $K$, and we are planning to investigate this further in the future.

5 Perturbation Theory in the Quench Parameter

Integrable model perturbation theory was developed in [52] for the study of integrable models subject to a small integrability-breaking perturbation. The case considered there was translation invariant in time. In a non-equilibrium protocol, such as a quench, the field theory action is no longer time-translation invariant. In [35] it was then observed that requiring factorization of the scattering at all times for such an action is consistent only if the latter is free. An approach to tackle the quench problem was also proposed in which the state in the Heisenberg picture after the quench could be expanded perturbatively in the quench parameter over the pre-quench quasi-particle basis. The approach requires the pre-quench theory to be integrable but allows for considering integrability breaking protocols.

Let us first review the main results of [35]. Consider an integrable quantum field theory with ground state $|0\rangle$ and action $\mathcal{A}_{0}$. At time $t=0$ the system is quenched and from $t=0$ onwards it is described by the new action

[TABLE]

where $\Psi(x,t)$ is some local field. In the interaction picture, with respect to the Hamiltonian of the pre-quench theory, the state of the system at infinite time after the quench is the time ordered exponential

[TABLE]

The state $|\psi_{0}\rangle$ in (88) can then be expanded perturbatively in $\lambda$ over the basis of the out-states of the pre-quench theory. $k$-particle states of this type are denoted by $|\theta_{1},\ldots,\theta_{k}\rangle_{\text{out}}$, with $\theta_{1}<\theta_{2}<\cdots<\theta_{k}$, being the rapidities. It can then be shown that integrability of the pre-quench theory allows for relaxing the constraint of ordering on the rapidities in the expansion over the out-states. In fact, the expansion [35]

[TABLE]

represents the state in the pre-quench basis in the Heisenberg picture at all times after the quench, up to first order in $\lambda$. $E_{0}(\theta)=m_{0}\cosh\theta$ and $P_{0}(\theta)=m_{0}\sinh\theta$ are the pre-quench energy and momenta of the particles, and

[TABLE]

is a $k$-particle form factor of the local field $\Psi$, calculated in the pre-quench quasi-particle basis. This state can then be employed to compute perturbative corrections to the one-point function of any local field $\Phi$ after the quench. These are found to be

[TABLE]

where [36]

[TABLE]

is a constant which is introduced to ensure that $\delta\langle\Phi(0)\rangle=0$ at first order in perturbation theory.

5.1 Perturbation Theory for the Entanglement Entropy

In order to calculate the quantity $\delta\langle{\cal T}(t)\rangle_{n}$, defined similarly to (91), we shall work in a replica version of (87); however this introduces a few changes. As discussed in section 3 particles are labelled by a replica index $j=1,\ldots,n$ and we will denote the replicated normalized pre-quench ground state by $|0\rangle_{n}=\otimes_{n}|0\rangle$. By repeating the steps leading to (89), it follows that the first order expansion of the state of the system in the replica theory after the quench is

[TABLE]

The expression (93) is essentially identical to (89) except for the prefactor $n$, which takes into account the sum over the replicas. Such a sum is however trivial since the local operator $\Psi$ when insterted in the $j$-th replica has only non-vanishing form factors among particles with copy index $j$. Similarly, the generalization of (91) and (92) for the twist field is also straightforward and given by

[TABLE]

where

[TABLE]

and $F_{k}$ are the form factors defined in (35), see section 3. The entanglement entropy may then be computed at first order in perturbation theory as

[TABLE]

Finally, we can define the first order correction to the Rényi entropies as

[TABLE]

Notice that in the perturbative approach, the pre-quench VEV (i.e. $\tau_{n}$) appears at the denominator of (97).

5.2 Entanglement Entropy Oscillations after a Small Mass Quench

Let us now evaluate (97) for a mass quench in the Ising field theory. In this case the field $\Psi(x,t)$ is the energy field, denoted by $\varepsilon(x,t)$, which has only a non-vanishing two-particle form factor with the pre-quench basis. Indeed, the pre-quench action $\mathcal{A}_{0}$ in (87) is obtained by perturbing the conformal invariant UV fixed point by the energy operator itself. The two-particle form factor, suitably normalized reads

[TABLE]

With the normalization choice for the energy form factor given in (98), we can directly identify [36] $\lambda$ in (94) with $\delta m\ll 1$, given in (10). From (94)–(96), and also recalling (33), the first order correction in $\delta m$ to the Rényi entropies after the quench can be easily calculated:

[TABLE]

The constant $C_{\mathcal{T}}^{n}$ can be determined exactly in this case [36] and it turns out to be

[TABLE]

where $\Delta_{n}$ is the scaling dimension of the twist field in (18). Observing that at first order in the quench parameter

[TABLE]

then from (96) and (99), the large time limit at first order in perturbation theory for the Rényi entropies finally follows

[TABLE]

Eq. (102) reproduces the main result in (14), up to $O(\mu)$ as anticipated in section 2. By expanding $\hat{K}$ around $\mu=0$, it is actually easy to verify that eq. (99) coincides with the first order of (74) at all times.

6 Lattice Results and Numerical Study in the Scaling Limit

In this section we present a detailed comparison of the field theory results obtained in section 4 against lattice numerical calculations in the Ising spin chain with Hamiltonian

[TABLE]

In the following, as already anticipated in the introduction, the lattice spacing will be denoted by $a$. In a lattice model, it is not possible to access directly the Rényi entropies $S_{n}(t)$ for a semi-infinite interval after a quench. Numerical techniques, based on the correlation matrix, are however known for calculating the Rényi entropies $S_{n}^{L}(t)$ of a subsystem of $L$ neighbouring sites with physical length $\ell:=La$, embedded into an infinite system (i.e. in the limit $N\rightarrow\infty$ in (103)). To extract the semi-infinite Rényi entropies, which we determined analytically in section 4, we then assume the validity of the same clustering property that holds for the spin-operator two-point function. The clustering property translates for the Rényi entropies into

[TABLE]

As also discussed in section 1, from a field theoretical perspective (see for instance (7)) eq. (104) is a consequence of locality of the branch point twist field, nevertheless it constitutes a non-trivial and new prediction when applied to the lattice model.

6.1 Correlation Matrix

For the Ising spin chain, the time evolution of correlation functions and entropies can be calculated using the restricted correlation matrix of a subsystem $A$ of $L$ sites [53, 5], embedded into an infinite system

[TABLE]

where111Note that the matrix $\Gamma_{L}^{A}$ has nothing to do with the decay rate $\Gamma$ introduced earlier in (12). Both these notations have been previously used in the literature so we maintain them here.

[TABLE]

and, for the Ising chain, we have

[TABLE]

Here, we already rewritten the transverse fields $h_{0}$ and $h$ in terms of the pre- and post-quench masses defined in the scaling field theory by: $h_{0}=1+am_{0}$, $h=1+am$. Since $m_{0},m$ are positive, the lattice calculations will be performed in the paramagnetic phase, the results should hold also in the ferromagnetic phase by duality. We also set the speed of light to $v=1$, therefore, according to (4), $J=\frac{1}{2a}$.

The matrix $\Gamma_{L}^{A}$ has $2L$ purely imaginary eigenvalues $\pm\mathrm{i}\nu_{k},\,k=1,\dots,L$, and the $2^{L}$ eigenvalues of the reduced density matrix $\rho_{A}$ matrix have the form

[TABLE]

where $a^{(j)}_{k}\in\{0,1\}$. A straightforward calculation gives the Rényi entropies for an interval of length $L$

[TABLE]

where the polynomials are

[TABLE]

Therefore the Rényi entropies can be easily calculated numerically by diagonalizing the correlation matrix.

6.2 Linear Growth in the Scaling Limit

Exact lattice results are available [7] for the leading, linear in time, contribution to the Rényi entropies $S_{n}^{L}(t)$. The linear growth is obtained in the regime $1\ll t\ll L$, while for $t\gg L$, according to a semi-classical quasi-particle picture the Rényi entropies saturate to a value proportional to the size $L$ of the interval.

We will then compare the field theoretical results of section 4 valid up to the second order in the quench parameter $\mu$, with the scaling limit of the lattice predictions in [7]. Computationally, see again (4) and the remarks below (107), the scaling limit is defined as follows: replace lattice quantities according to (107), then introduce the continuum momentum variable $p$ substituting $\varphi:=pa$, and $0\leq p\leq\frac{2\pi}{a}$, and eventually take the limit $a\rightarrow 0$. The mathematical operation will be denoted by the shorthand notation $\lim_{\text{scal}}$. For instance, it is easy to verify that

[TABLE]

The main result in [7], evaluated for $L\rightarrow\infty$ then reads

[TABLE]

with $\epsilon^{\prime}_{\varphi}:=\frac{d\epsilon_{\varphi}}{d\varphi}$ and $P_{n}$ as in (110). The $lin$ subscript indicates that the formula only captures the linear growth part of the entanglement. With the definition

[TABLE]

the scaling limit of (112) is thus

[TABLE]

Expanding (114) for small quenches (i.e. $m=m_{0}+\delta m$) and substituting $p=m_{0}\sinh\theta$, which is consistent at second order in $\delta m$, we find

[TABLE]

where we used the definition (10). It then follows, as expected according to (104), that the result in (115) is precisely twice the leading large-time asympotics obtained expanding (14) for a small quench, see in particular (15). By considering the limit $n\rightarrow 1$ in (112), the scaling limit of the von Neumann entropy turns out to be

[TABLE]

and expanding for small quenches

[TABLE]

In the scaling limit, and for small quenches, the von Neumann entropy determined in [7] is dominated by a term $O(\delta m^{2}\log\delta m)$ and therefore is not analytic in the quench parameter. This unexpected result provides another indication that the limit $n\rightarrow 1$ in the Rényi entropies does not commute with a pertubative expansion in $\delta m$. Incidentally, the emergence of logarithmic corrections to the expectation values of certain fields in the massive Ising field theory is compatible with previous studies, such as  [54]. It is generally due to ambiguities in the definition of some local operators: For instance in the Ising field theory, the energy field can be regarded as a linear combination of the usual fermion bilinear and a term proportional to the identity field with proportionality constant equal to the mass.

6.3 Numerical Evaluation of the Correlation Matrix

In this section we test the field theory results against the numerical results on the lattice, where one can directly diagonalize the correlation matrix (105) and calculate the entropies as shown in section 6.1.

We expect that the results match in the scaling limit defined in the previous section and for small quenches. Therefore we chose the transverse field to be close to the critical value i.e. $m,m_{0}\ll 1$ and the quench $\delta m\ll 1$. Then the scaling limit can be carried out by decreasing $a$ while keeping the physical subsystem size $\ell=La$ fixed. Then one can extrapolate to $a=0$ taking into account corrections to the scaling of the entropies as discussed in Appendix A.

This method gives the entanglement entropies of a finite subsystem, which is proportional to the logarithm of the two point function of branch point twist fields. Therefore one needs to consider large enough subsystem sizes in order to observe clustering, namely factorization into a one-point function squared. All our numerical results show excellent agreement with analytical predictions up to a factor two due to clustering. We note that for larger subsystems and times the evaluation of (105) gets harder due to the highly oscillatory integrands in (106).

As recalled in (19), in a massive field theory the logarithmic divergence of the von Neumann entropy is encoded in the term

[TABLE]

where $c$ is the central charge of the ultraviolet CFT and the $O(1)$ corrections are discussed in Appendix A. The scaling limit technique can be then used to extract the central charge $c$. For instance at values of the mass $m=0.04$ we obtained $c=0.50195(3)$, which is very close to the theoretical value $c=\frac{1}{2}$. The central charge extrapolation provides a means to numerically probe the scaling regime of the Ising spin chain. We found that differences of entropies calculated at different times do not have any divergences in the scaling limit as expected. More details on the numerical results corresponding to the scaling limit can be found in Appendix A.

6.3.1 Saturation and Oscillations

As was already pointed out in [5] for a finite subsystem size the entanglement entropy saturates to a constant after a finite time. The saturation constant is linear in the subsystem size. The authors studied large quenches, where the leading behaviour is the linear growth, and there was no trace of oscillatory behaviour. The left panel of fig. 2 shows the time evolution of the von Neumann entropy of different subsytems for a small quench with fixed lattice spacing. It is clear that also after saturation the entropy continues to oscillate. The baseline of the oscillations saturates as well, but apart from this offset, the functional form is predicted well by the field theory formula (39), which in principle is not supposed to be valid for $t\gg L$.

For a comparison we plot the entanglement entropy together with (39) shifted by an arbitrary constant. In the right panel of fig. 2 we plot the differences of the entanglement entropy, calculated at different subsystem sizes. After the saturation the curves are equally spaced, which shows that the oscillations do not depend on the subsystem size.

From fig. 2 we can draw several conclusions:

• •

Before saturation, the values of the entanglement entropy are independent of the subsystem’s size. This demonstrates the clustering of the two point function of twist fields.

• •

The saturation times and saturation values are equally spaced for different subsystem sizes (with fixed difference in the size). This shows the $\propto L$ behaviour of the saturation values, which was already discussed in [5, 7].

• •

The oscillations are present, independently of the linear growth and the saturation. After the saturation sets in, the shape of the oscillations is the same for different subsystem sizes. Moreover, they persist for large times and are well reproduced by the formula (39) up to a constant offset, and a factor of two. As already discussed, the factor of two is the result of the clustering.

6.3.2 Linear Growth and Oscillations

To observe linear growth in time and test the field theory result (74), one needs larger subsystem sizes in order to prevent saturation within the time window. Fig. 3 shows the time evolution of the Rényi entropies after a mass quench. The theoretical prediction (74) is in remarkably good agreement with the numerical data. One can also see how the successive contributions of the different terms in (74) improve accuracy.

Within our field theoretical approach, as mentioned in section 2, for the von Neumann entropy we can obtain only the oscillatory behaviour. Additional contributions are also expected to be present in this case. One can then take the time derivative in order to eliminate the time independent offset and subtract the value predicted by (117) to eliminate the linear growth, which, in turn, produces an offset in the time derivative. The results can be seen in fig. 4. The agreement with the field theory prediction is again very good except for the small $t$ region. Notice that the corrections coming from $C^{(2)}_{40}(t)+C^{(2)}_{04}(t)$ and $R(t)$ can not be calculated at $n=1$ and, although subleading for large time, they might affect the small time behaviour.

7 Conclusion

In this paper we have presented an analytic derivation of the leading large-time post-quench dynamics of entanglement in the massive Ising field theory. We considered in particular a global quench resulting from a sudden change in the mass of the fermionic particle, from an initial value $m_{0}$ at time $t=0$ to a subsequent value $m$ for $t>0$. For the first time in a dynamical context for massive quantum field theories, we have employed the branch point twist field approach [34] in our computations. We have computed the Rényi entropies of a semi-infinite interval, which are proportional to the logarithm of the one-point function of a branch point twist field in a replica quantum field theory. In particular, the twist field one-point function has been computed exactly up to $O(K^{2})$, in the post-quench quasi-particle expansion of the initial state by employing a regularization scheme for the infinite volume divergences discussed in [31].

Such an expansion can be also recast as a perturbative series in the quench parameter $\delta m:=m-m_{0}$, and is then exact up to $O(\delta m^{2})$. At first order in the quench parameter $\delta m$ the result for the twist field one-point function can also be recovered from a perturbative expansion in the pre-quench quasi-particle basis by generalizing the approach introduced in [35]. We demonstrated, moreover, that crucial effects of the relaxation dynamics, such as linear growth of entanglement must manifest as second order corrections in such a perturbative expansion. The main conclusions from the analytic results can be then summarized as follows:

• •

We showed the presence of a contribution to the Rényi entropies which grows linearly in time with slope $\frac{n\Gamma m}{2(n-1)}$, where $\Gamma$ is up to $O(K^{2})$ exactly the decay rate of the spin operator after a mass quench found in [31].

• •

The Rényi and von Neumann entropies contain contributions which are oscillatory, with frequency of oscillation $2m$ and amplitude proportional to $(mt)^{-3/2}$ for large-time. Those contributions are of first order in $K$ and can be also obtained from a perturbative expansion in the quench parameter $\delta m$. Our result implies that oscillations in the entanglement entropies are not produced by finite size-effects, as for instance stated in [7], but are rather inherent properties of those quantities. The same reasoning will apply to the constant shift of the Rényi entropies, which we also determined in section 4.

• •

We have provided a simple argument to show that, up to a constant normalization by the VEV of the twist field, the one-point function can be expressed as the exponential of a Laurent polynomial in the variables $(mt)^{\frac{p}{2}}$ for $p\leq 2$ and $p\neq 1$. We have shown this at $O(K^{2})$ and expect to generalize this conclusion to higher orders in the future. Interestingly, the arguments leading to exponentiation of the one-point function also apply to the order parameter discussed in [31].

As expected, our field-theoretical results for the linear large-time behaviour of the Rényi entropies reproduce the scaling limit of the formulae found in [7] up to $O(\delta m^{2})$. The field-theoretical expansion, however, extends the lattice results to intermediate times and is confirmed with remarkable accuracy by numerical calculations directly in the scaling limit. Comparison between analytic and numerical results also shows that the two-point function of branch point twist fields after a global quench satisfies clustering, as previously observed for the spin field [8, 9, 31]. Thus the entanglement entropies are proportional to the number of subsystem boundary points, just as in equilibrium situations.

It would be interesting to use twist fields and the approaches discussed in section 3.2 and section 5 to consider other (small) global quenches. In particular, quenches that drive the theory away from an integrable point or to a different interacting integrable model. A particular case is the quench of the longitudinal magnetic field in the Ising spin chain while fixing the transverse field at its critical value $h=1$. In the scaling limit, this corresponds to a mass quench in the so-called minimal $E_{8}$ Toda theory. The prediction is then that the entanglement entropies will oscillate with frequencies that are directly the quasi-particle masses of the $E_{8}$ field theory. An analogous phenomenon has been observed numerically in the Ising spin chain [28, 29, 30], for a quench of the longitudinal magnetic field, but in the ferromagnetic phase $h<1$. The presence of oscillations in the entanglement entropies and their slower (linear) growth in time were ascribed to the confinement of the kinks.

Finally, it would be useful to develop a quasi-particle interpretation/derivation of the oscillatory contributions to entanglement. Even though our results are restricted to the Ising field theory, the emergence of such oscillations in the context of form factor expansions seems very natural. This suggests that it is a universal feature of quenches in gapped theories.

A powerful unifying picture emerges from our work: the dynamics of entanglement and that of correlators of local fields after a global quench are not fundamentally distinct. Rather, the dynamics of entanglement is just the dynamics of correlators of a particular field, the branch point twist field. As a consequence, the large time linear growth of entanglement emerging from the quasi-particle picture of [5] is nothing but the exponential decay of correlators (in our case, the one-point function) at large time after the quench. Suggestively, out-of-equilibrium dynamics where no indication of linear growth of entanglement is observed, such as in the presence of confinement [28, 29, 30], could signal that certain local observables fail to relax exponentially fast at large times. The present work extends the seminal results of [3] out-of-criticality and, for a very simple model, provides further evidence of the rich and interesting dynamics of correlators in out-of-equilibrium massive QFT.

Acknowledgements:

We are indebted to Aldo Delfino, Fabian Essler and Dirk Schuricht for discussions and helpful e-mail exchanges regarding their works [35, 31], and to Gábor Takács for sharing his insights into the wider interpretation and applicability of our results.

OCA and IMSZ are also grateful to Benjamin Doyon and Cecilia De Fazio for many useful discussions and feedback. They gratefully acknowledge support from EPSRC through the standard proposal “Entanglement Measures, Twist Fields, and Partition Functions in Quantum Field Theory” under reference number EP/P006108/1 and from the International Institute of Physics in Natal (Brazil). This project was initiated during the Workshop on Transport in Strongly Correlated Systems held there in the summer of 2018. OCA’s research was partly supported by an Emmy Noether Visiting Fellowship of the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science.

ML is grateful to Fábio Novaes for general discussions. ML and JV acknowlegde financial support from the Brazilian ministries MEC and MCTIC.

Appendix A Numerics in the Scaling Limit

A.1 Extrapolations and Extraction of the Central Charge and Scaling Dimensions

In this appendix we present further details on the numerical scaling limit discussed in section 6.3. The central charge and the operator scaling dimensions can be extracted from the evaluation of the restricted correlation function at $t=0$. We fix $m$ to a certain value in the scaling field theory, and change $a$ and $L$ in such a way that the physical subsystem size $\ell=aL$ is kept fixed. Then one can fit the lattice spacing dependence of the logarithm of the one-point function of the disorder operator and the Rényi entropies with the following functional forms

[TABLE]

For the Rényi entropies unusual corrections are present [55]. For the disorder operator one assumes standard corrections, more on this operator can be found in section A.2.

The coefficient of the $\log a$ term corresponds to the scaling dimension of the operators or in the case of the von Neumann entropy the central charge of the UV CFT. We carried out the fit with the following parameters: $m=0.04,\,\ell=128,\,a=1/4,1/7,1/8\dots 1/20$. The results are summarized in Table 1.

The agreement for the central charge, the dimension of the disorder operator, and of the twist field with $n=2,3,4$ is excellent (see equation (18)). Note that the central charge can be extracted with better precision by calculating the von Neumann entropy at the critical point, with fixed lattice spacing and changing the number of sites in the subsystem, based on the logarithmic violation of the area law. In our case we extract the UV central charge away from the critical point, therefore we have less precision.

Using the above fits one can extrapolate to $a=0$. At different times, we used the same set of lattice spacings to carry out the extrapolations. Note that in our comparison we subtract the post-quench ground state entropy from the numerical results. The logarithmic singularity cancels from these differences, therefore we did not include the logarithm when extrapolating these quantities.

Note that going closer to the critical point would require more computational power since one has to increase the subsystem size correspondingly, therefore the size of the correlation matrix increases. The calculation of the intergrals (106) gets also more difficult. One also has to make sure, that the chosen subsystem size is large enough for the clustering of the two-point functions. For the quench studied in this paper we checked this using the saturation of the post-quench entropies. We found that for $a=1$ and $L\approx 120$ the entropies are saturated up to $O(10^{-6})$, therefore we claim that our numerical results for the entropies have errors of this order.

A.2 Scaling Limit and the Disorder Operator

In [8, 9] the decaying exponential characterizing the post-quench behaviour of the spin operator was found to be

[TABLE]

therefore, in the scaling limit, see (113)

[TABLE]

For small quenches, eq. (122) becomes

[TABLE]

that agrees with the field theory result presented in [31], when expanded up to the second order in the quench parameter.

In [9] the authors determine the time dependence of the order parameter by calculating the determinant of the correlation matrix. However, their definition of the correlation matrix is slightly different. In particular, in our notations, their Töplitz matrix starts with $\Pi_{-1}$ in the upper left corner. From [56] one can see that this can be absorbed into the redefinition of $h\rightarrow 1/h$, realizing the Kramers–Wannier duality. Therefore calculating the determinant of (105) gives the square of the two-point function of the disorder operator, up to a constant:

[TABLE]

Using the fitting procedure outlined in section A.1 we obtained $\Delta_{\mu}\approx 0.12497$, which is very close to the theoretical value $\frac{1}{8}=0.125$.

If the separation is large, the two-point function of $\mu$ clusters, just as for the order parameter in the ferromagnetic phase [9]

[TABLE]

Therefore for large enough separations one can get access to the one-point function. It can be also seen that in the scaling limit $\tilde{\mu}(t)=\log\langle\tilde{0}|\mu(0,t)|\tilde{0}\rangle-\log\langle\tilde{0}|\mu(0,0)|\tilde{0}\rangle$ has a finite limit. Based on the Kramers–Wannier duality the formulas of [31] for the order parameter in the ferromagnetic phase can be directly used to test the disorder operator in the paramagnetic phase. Such a comparison can be seen in fig. 5. The agreement is excellent. Note that in the case of the order/disorder operator there is no offset at $O(K^{2})$ [31].