Entanglement hamiltonians in 1D free lattice models after a global quantum quench

TL;DR

This paper investigates how the entanglement Hamiltonian evolves over time after a global quantum quench in 1D free lattice models, revealing insights into entanglement dynamics and spectrum gaps.

Contribution

It provides a detailed analysis of the temporal evolution of entanglement Hamiltonians and spectra in free lattice models following a global quench, combining numerical and conformal field theory approaches.

Findings

Entanglement spectrum gaps evolve distinctly over time.

Contours for entanglement entropy are characterized and analyzed.

Conformal field theory and quasi-particle models explain the entanglement dynamics.

Abstract

We study the temporal evolution of the entanglement hamiltonian of an interval after a global quantum quenchin free lattice models in one spatial dimension. In a harmonic chain we explore a quench of the frequency parameter. In a chain of free fermions at half filling we consider the evolution of the ground state of a fully dimerised chain through the homogeneous hamiltonian. We focus on critical evolution hamiltonians. The temporal evolutions of the gaps in the entanglement spectrum are analysed. The entanglement hamiltonians in these models are characterised by matrices that provide also contours for the entanglement entropies. The temporal evolution of these contours for the entanglement entropy is studied, also by employing existing conformal field theory results for the semi-infinite line and the quasi-particle picture for the global quench.

Figures

Equations

$${\rho }_A=\frac{e^{-{\hat{K}}_A}}{{\mathcal{Z}}_A}$$

$$\hat{H}=\sum _{i=0}^{L-1}\left(\frac{1}{2m}{\hat{p}}_i^2+\frac{m{\omega }^2}{2}{\hat{q}}_i^2+\frac{\kappa }{2}({\hat{q}}_{i+1}-{\hat{q}}_i{)}^2\right)$$

$$J\equiv\bigg{(}\begin{array}[]{cc}\boldsymbol{0}&\boldsymbol{1}\\ -\,\boldsymbol{1}&\boldsymbol{0}\\ \end{array}\bigg{)}$$

$$\widehat{K}_{A}=\frac{1}{2}\,\boldsymbol{\hat{r}}^{\textrm{t}}H_{A}\,\boldsymbol{\hat{r}}\qquad\boldsymbol{\hat{r}}=\bigg{(}\begin{array}[]{c}\boldsymbol{\hat{q}}\\ \boldsymbol{\hat{p}}\end{array}\bigg{)}$$

$${\gamma }_A=W^{\text{t}}{\mathcal{D}}_{\text{d}}W$$

$$S_A=\sum _{k=1}^{\ell }s({\sigma }_k)S_A^{(n)}=\sum _{k=1}^{\ell }{s}_n({\sigma }_k)$$

$$s(y)\equiv (y+1/2)\log (y+1/2)-(y-1/2)\log (y-1/2)$$

$$s_{n}(y)\equiv\frac{1}{n-1}\,\log\big{[}(y+1/2)^{n}-(y-1/2)^{n}\big{]}.$$

$$(\text{i}J{\gamma }_A{)}^2=W^{-1}{\mathcal{D}}_{\text{d}}^{2}W(\text{i}{\gamma }_{A}J{)}^2={\tilde{W}}^{-1}{\mathcal{D}}_{\text{d}}^2\tilde{W}$$

$$\tilde{W}\equiv J^{\text{t}}WJ=W^{-\text{t}}.$$

$$H_A=W_H^{\text{t}}{\mathcal{E}}_{\text{d}}{W}_H$$

$$\mathcal{E}=2\,\textrm{arccoth}(2\,\mathcal{D})=\log\!\bigg{(}\frac{\mathcal{D}+1/2}{\mathcal{D}-1/2}\bigg{)}\,.$$

$$H_A=2\text{i}J\text{arccoth}(2\text{i}{\gamma }_{A}J)=2\text{i}\text{arccoth}(2\text{i}J{\gamma }_A)J.$$

$$W_H=\tilde{W}$$

$$H_{A}=J^{\textrm{t}}\gamma_{A}\,J\;h\big{(}\sqrt{(\textrm{i}\gamma_{A}J)^{2}}\,\big{)}=h\big{(}\sqrt{(\textrm{i}J\gamma_{A})^{2}}\,\big{)}\,J^{\textrm{t}}\gamma_{A}\,J$$

$$\widehat{K}_{A}=\frac{1}{2}\;\boldsymbol{\hat{s}}^{\textrm{t}}\,\mathcal{E}_{\textrm{\tiny d}}\,\boldsymbol{\hat{s}}\;\;\qquad\;\;\widetilde{W}\boldsymbol{\hat{r}}\equiv\boldsymbol{\hat{s}}\equiv\bigg{(}\!\begin{array}[]{c}\hat{\boldsymbol{\mathfrak{q}}}\\ \hat{\boldsymbol{\mathfrak{p}}}\end{array}\!\bigg{)}$$

$$\boldsymbol{\hat{b}}\equiv\Omega^{-1}\boldsymbol{\hat{s}}\;\;\;\qquad\;\;\;\hat{\mathfrak{b}}_{k}\equiv\frac{\hat{\mathfrak{q}}_{k}+\textrm{i}\,\hat{\mathfrak{p}}_{k}}{\sqrt{2}}\qquad\Omega\equiv\frac{1}{\sqrt{2}}\bigg{(}\!\!\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{1}\\ -\textrm{i}\boldsymbol{1}&\textrm{i}\boldsymbol{1}\\ \end{array}\!\bigg{)}$$

$$\widehat{K}_{A}=\frac{1}{2}\;\boldsymbol{\hat{b}}^{\textrm{t}}\bigg{(}\!\begin{array}[]{cc}\boldsymbol{0}&\mathcal{E}\\ \mathcal{E}&\boldsymbol{0}\\ \end{array}\!\bigg{)}\,\boldsymbol{\hat{b}}=\sum_{k=1}^{\ell}\varepsilon_{k}\!\left(\hat{\mathfrak{b}}_{k}^{\dagger}\,\hat{\mathfrak{b}}_{k}+\frac{1}{2}\,\right).$$

$${\rho }_A=\sum _{\mathbf{n}}\frac{e^{-{\sum }_{k=1}^{\ell }{\epsilon }_k(n_k+1/2)}}{{\mathcal{Z}}_A}|\mathbf{n}⟩⟨\mathbf{n}|\equiv \sum _{\mathbf{n}}{\lambda }_{\mathbf{n}}|\mathbf{n}⟩⟨\mathbf{n}|$$

$$\mathcal{Z}_{A}=\sum_{\boldsymbol{m}}\langle\boldsymbol{m}|\bigg{(}\!\sum_{\boldsymbol{n}}e^{-\!\sum_{k=1}^{\ell}\varepsilon_{k}(n_{k}+1/2)}|\boldsymbol{n}\rangle\langle\boldsymbol{n}|\bigg{)}|\boldsymbol{m}\rangle=\prod_{k=1}^{\ell}\frac{e^{-\varepsilon_{k}/2}}{1-e^{-\varepsilon_{k}}}=\prod_{k=1}^{\ell}\sqrt{\sigma_{k}^{2}-1/4}$$

$${\lambda }_{\mathbf{n}}=\prod _{k=1}^{\ell }\left[\left(1-e^{-{\epsilon }_k}\right)e^{-n_k{\epsilon }_k}\right]=\prod _{k=1}^{\ell }\frac{1}{{\sigma }_k+1/2}{\left(\frac{{\sigma }_k-1/2}{{\sigma }_k+1/2}\right)}^{n_k}.$$

$${\lambda }_{\mathbf{n}}={\lambda }_{\text{max}}{e}^{-{\sum }_{k=1}^{\ell }{n}_k{\epsilon }_k}={\lambda }_{\text{max}}\prod _{k=1}^{\ell }{\left(\frac{{\sigma }_k-1/2}{{\sigma }_k+1/2}\right)}^{n_k}.$$

$$g_{\mathbf{m},\mathbf{n}}\equiv \log {\lambda }_{\mathbf{m}}-\log {\lambda }_{\mathbf{n}}=\sum _{k=1}^{\ell }(n_k-m_k){\epsilon }_k=\sum _{k=1}^{\ell }(m_k-n_k)\log \left(\frac{{\sigma }_k-1/2}{{\sigma }_k+1/2}\right).$$

$$\gamma_{A}\equiv\bigg{(}\begin{array}[]{cc}Q&R\\ R^{\textrm{t}}&P\\ \end{array}\bigg{)}\;\;\;\qquad\;\;\;H_{A}\equiv\bigg{(}\begin{array}[]{cc}M&E\\ E^{\textrm{t}}&N\\ \end{array}\bigg{)}\,.$$

$$W\equiv\bigg{(}\begin{array}[]{cc}U&Y\\ Z&V\\ \end{array}\bigg{)}\;\;\;\qquad\;\;\;\widetilde{W}\equiv\bigg{(}\begin{array}[]{cc}\widetilde{U}&\widetilde{Y}\\ \widetilde{Z}&\widetilde{V}\\ \end{array}\bigg{)}$$

$$\bigg{(}\begin{array}[]{cc}\widetilde{U}&\widetilde{Y}\\ \widetilde{Z}&\widetilde{V}\\ \end{array}\bigg{)}=\bigg{(}\begin{array}[]{cc}V&\!\!-Z\\ -Y&\!\!U\\ \end{array}\bigg{)}\,.$$

$$\gamma_{A}=\bigg{(}\begin{array}[]{cc}U^{\textrm{t}}\,\mathcal{D}\,U+Z^{\textrm{t}}\,\mathcal{D}\,Z&U^{\textrm{t}}\,\mathcal{D}\,Y+Z^{\textrm{t}}\,\mathcal{D}\,V\\ Y^{\textrm{t}}\,\mathcal{D}\,U+V^{\textrm{t}}\,\mathcal{D}\,Z&Y^{\textrm{t}}\,\mathcal{D}\,Y+V^{\textrm{t}}\,\mathcal{D}\,V\end{array}\bigg{)}$$

$$H_{A}=\bigg{(}\begin{array}[]{cc}V^{\textrm{t}}\,\mathcal{E}\,V+Y^{\textrm{t}}\,\mathcal{E}\,Y&-\,V^{\textrm{t}}\,\mathcal{E}\,Z-Y^{\textrm{t}}\,\mathcal{E}\,U\\ -\,Z^{\textrm{t}}\,\mathcal{E}\,V-U^{\textrm{t}}\,\mathcal{E}\,Y&Z^{\textrm{t}}\,\mathcal{E}\,Z+U^{\textrm{t}}\,\mathcal{E}\,U\end{array}\bigg{)}\,.$$

$$S_A^{(n)}=\sum _{i\in A}{s}_A^{(n)}(i)s_A^{(n)}(i)⩾0.$$

$$\sum _{i=1}^{\ell }{p}_k(i)=1p_k(i)⩾0.$$

Full text

Contents

@starttoc

toc

1 Introduction

The entanglement associated to the state of a quantum system and to a bipartition of its Hilbert space has attracted a lot of attention during the last decade. Many techniques have been developed to quantify entanglement in quantum many-body systems, quantum field theories and quantum gravity [1, 2, 3, 4, 5, 6, 7].

Given a bipartition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ of the Hilbert space and a state characterised by the density matrix $\rho$, the reduced density matrix associated to $\mathcal{H}_{A}$ is $\rho_{A}\equiv\textrm{Tr}_{\mathcal{H}_{B}}\rho$. In this manuscript we focus on bipartitions of Hilbert spaces that come from spatial bipartitions; hence $A$ denotes a spatial region and $B$ its complement. [8]. The reduced density matrix is a positive semi-definite hermitian operator normalised to $\textrm{Tr}_{\mathcal{H}_{A}}\rho_{A}=1$ (hereafter the trace will be always over $\mathcal{H}_{A}$). These properties allow to write $\rho_{A}$ as

[TABLE]

where the hermitian operator $\widehat{K}_{A}$ is the entanglement hamiltonian (also called modular hamiltonian) and $\mathcal{Z}_{A}={\rm Tr}(e^{-\widehat{K}_{A}})$.

The moments of the reduced density matrix ${\rm Tr}\rho_{A}^{n}$, that are parameterised by the integer power $n\geqslant 2$, provide the Rényi entropies $S_{A}^{(n)}\equiv\tfrac{1}{1-n}\log({\rm Tr}\rho_{A}^{n})$. The entanglement entropy is defined as the Von Neumann entropy of the reduced density matrix $S_{A}=-{\rm Tr}(\rho_{A}\log\rho_{A})$. A very important property of the Rényi entropies is that $S_{A}^{(n)}\to S_{A}$ as $n\to 1$. This limit can be very complicated because it requires to perform an analytic continuation in $n$. The entanglement entropies are $S_{A}^{(n)}$ for $n\geqslant 1$, where $S_{A}^{(1)}\equiv S_{A}$ is assumed. These scalar quantities depend only on the eigenvalues of $\rho_{A}$ (the entanglement spectrum) [9]; hence the entanglement hamiltonian contains more information about the entanglement of the bipartition with respect to the corresponding entanglement entropies.

In quantum field theories, Bisognano and Wichmann found the explicit expression of the entanglement hamiltonian of half space through the stress tensor of the model [10]. In conformal field theories (CFT), this formula has been employed as the starting point to obtain entanglement hamiltonians for other interesting configurations [11, 12, 13, 14, 15, 16, 17, 18]. Crucial tools in this analysis are the methods developed for CFT with boundary [19]. In quantum gravity models in the context of the gauge/gravity correspondence, elaborating on the holographic prescription for the entanglement entropy [20, 21], also the entanglement hamiltonian has been explored [22].

In many body quantum systems, entanglement hamiltonians have been studied for bosonic Gaussian states in harmonic lattices and for fermionic Gaussian states in lattice models of free fermions [1, 2, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. The entanglement hamiltonian is a quadratic operator in these free models; hence it is fully characterised by a matrix. The results available in the literature are mostly based on numerical analysis of this matrix. Analytic expressions for the entanglement hamiltonian in the thermodynamic limit have been obtained in [31] for the static configuration where $A$ is a single block in the infinite chain of free fermions in the ground state. The matrix characterising this entanglement hamiltonian at half filling has non vanishing odd diagonals at any allowed distance from the main diagonal. The continuum limit of this entanglement hamiltonian has been studied in [34] and the expected CFT result of [11, 12] has been derived through this method. In spin chain systems, entanglement hamiltonians have been studied in [35].

The entanglement spectrum is another very insightful quantity to consider in order to get insights about the bipartite entanglement. It has been explored by employing both CFT methods [15, 36] and numerical techniques [37, 38]. The entanglement spectrum is not enough to construct the entanglement hamiltonian.

The contours for the entanglement entropies are simpler to obtain with respect to the entanglement hamiltonian and they are not fully determined by the entanglement spectrum as well. The contours for the entanglement entropies are positive functions on the subsystem $A$ such that their integral over $A$ provides the entanglement entropies. The aim of these scalar quantities is to describe the contribution of the various sites in $A$ to the entanglement entropies of the bipartition and the entanglement spectrum is not enough to determine them. In free lattice models, contours for the entanglement entropies have been constructed and numerical analysis have been performed [39, 40, 41, 42, 16]. For some static configurations described by CFT, they are related to the corresponding entanglement hamiltonians [42]. Density functions for the entanglement entropy have been explored also in [43]. The contours for the entanglement entropies are special density functions because further constraints are required for them [40, 42]. Nonetheless, a complete list of properties that allows to select a unique contour for the entanglement entropies is not known.

The bipartite entanglement has been studied during the last decade also to explore the out of equilibrium dynamics of quantum systems. Global quantum quenches are insightful processes where the system is prepared in the ground state $|\psi_{0}\rangle$ of a translationally invariant hamiltonian and at $t=0$ a sudden change modifies the hamiltonian keeping the translation invariance. The non trivial dynamics is due to the fact that the initial state $|\psi_{0}\rangle$ is not an eigenvector of the hamiltonian determining its unitary evolution. The temporal evolutions of various quantities after global quantum quenches and also other kinds of quenches have been studied in quantum field theories, in many-body systems on the lattice [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61] and also in quantum gravity through the gauge/gravity correspondence [21, 62] (see the reviews [63, 64, 65] for an extensive list of references).

In this manuscript we study the temporal evolution of the entanglement hamiltonian and of a contour for the entanglement entropy after a global quantum quench in free lattice models. For simplicity, we consider, in one spatial dimension, a harmonic chain and a chain of free fermions. In the harmonic chain, we explore a quench of the frequency parameter such that the unitary evolution is governed by the massless hamiltonian. In the chain of free fermions, we consider the global quench introduced in [49, 1], where the system is prepared in the ground state of a fully dimerised chain while the evolution hamiltonian is fully homogeneous.

This manuscript is organised as follows. In §2 and §3 we discuss the construction of the entanglement hamiltonians and of the contours for the entanglement entropies in the models of interest. In §4 we provide a naive formula for the contour function obtained through CFT results. The numerical results for the harmonic chain and for the chain of free fermions are described in §5 and §6 respectively. In §7 we draw some conclusions. The Appendices A, B, C and D contain technical details about some issues discussed in the main text.

2 Harmonic lattices

In this section we review the construction of the matrix characterising the entanglement hamiltonian for Gaussian states in harmonic lattices. We focus on one dimensional chains, but the discussion can be easily extended to harmonic lattices in higher dimensions.

The hamiltonian of the harmonic chain with nearest neighbour spring-like interaction reads

[TABLE]

where $L$ is the number of lattice sites and the hermitian operators $\hat{q}_{i}$ and $\hat{p}_{i}$ satisfy the canonical commutation relations $[\hat{q}_{i},\hat{q}_{j}]=[\hat{p}_{i},\hat{p}_{j}]=0$ and $[\hat{q}_{i},\hat{p}_{j}]=\textrm{i}\delta_{ij}$ (throughout this manuscript $\hbar=1$). The boundary conditions, that are crucial to determine the expression of the correlators, do not influence explicitly the following discussion. By arranging the position and momentum operators into the vector $\boldsymbol{\hat{r}}\equiv(\hat{q}_{1},\dots,\hat{q}_{L},\hat{p}_{1},\dots,\hat{p}_{L})^{\textrm{t}}$, the canonical commutation relations can be written in the form $[\hat{r}_{i},\hat{r}_{j}]=\textrm{i}J_{ij}$, being $J$ the standard symplectic matrix

[TABLE]

where $\boldsymbol{1}$ is the $L\times L$ identity matrix and $\boldsymbol{0}$ is the $L\times L$ matrix made by zeros. Notice that $J^{\textrm{t}}=-J$ and $J^{2}=-\boldsymbol{1}$. We also need that $J^{\textrm{t}}(a\oplus b)J=b\oplus a$.

The linear transformations $\boldsymbol{\hat{r}}\to\boldsymbol{\hat{r}}^{\prime}=S\boldsymbol{\hat{r}}$ preserving the canonical commutation relations define the real symplectic group $\textrm{Sp}(L)$, made by the real $2L\times 2L$ matrices $S$ that satisfy the relation $SJS^{\textrm{t}}=J$ [66]. Given a symplectic matrix $S$, we have that $\textrm{det}(S)=1$, $S^{\textrm{t}}\in\textrm{Sp}(L)$ and $S^{-1}=JS^{\textrm{t}}J^{-1}$. The above observations imply that $S^{-\textrm{t}}=J^{\textrm{t}}SJ$, where we have adopted the notation $M^{-\textrm{t}}\equiv(M^{\textrm{t}})^{-1}$.

By employing a canonical transformation, the hamiltonian (2) can be written as the hamiltonian of a free boson with mass $\omega$ discretised on a lattice with spacing $a=\sqrt{m/\kappa}$. This implies that we can set $m=\kappa=1$ without loss of generality. The continuum limit of this model gives the free scalar boson with mass $\omega$ in two spacetime dimensions. In the massless regime, this quantum field theory is a conformal field theory (CFT) with central charge $c=1$.

In this manuscript we focus on Gaussian states of (2), which are completely characterised by the correlators $\langle\hat{r}_{i}\rangle$ (first moments) and $\langle\hat{r}_{i}\,\hat{r}_{j}\rangle$ (second moments). Since a shift of the first moments corresponds to a unitary transformation that preserves the Gaussian nature of the state, let us consider Gaussian states having $\langle\hat{r}_{i}\rangle=0$. Thus, the second moments fully describe these states and they can be collected into the $2L\times 2L$ covariance matrix $\gamma\equiv\textrm{Re}\langle\boldsymbol{\hat{r}}\,\boldsymbol{\hat{r}}^{\textrm{t}}\rangle$, which is a real, symmetric and positive definite matrix [5, 6, 67, 39, 66]. A canonical transformation $\boldsymbol{\hat{r}}\to\boldsymbol{\hat{r}}^{\prime}=S\boldsymbol{\hat{r}}$ characterised by the symplectic matrix $S$ induces the transformation $\gamma\to\gamma^{\prime}=S\gamma S^{\textrm{t}}$ on the covariance matrix. The covariance matrix of a pure Gaussian state satisfies the relation $(\textrm{i}J\,\gamma)^{2}=\tfrac{1}{4}\,\boldsymbol{1}$.

Given the harmonic chain (2) in a Gaussian state $\rho$ characterised by the covariance matrix $\gamma$, let us introduce a bipartition of the Hilbert space $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ corresponding to a spatial bipartition $A\cup B$. The reduced density matrix $\rho_{A}\equiv\textrm{Tr}_{\mathcal{H}_{B}}\rho$ associated to the spatial subsystem $A$ characterises a mixed state also when the whole system $A\cup B$ is in a pure state. For the harmonic chain (2) and its higher dimensional generalisations, $\rho_{A}$ remains Gaussian for any choice of $A$. This implies that $\rho_{A}$ is fully described by the reduced covariance matrix $\gamma_{A}\equiv\textrm{Re}\langle\boldsymbol{\hat{r}}\,\boldsymbol{\hat{r}}^{\textrm{t}}\rangle|_{A}$, obtained by extracting from $\gamma$ of the entire system the rows and the columns corresponding to the lattice sites belonging to the subsystem $A$. The reduced covariance matrix $\gamma_{A}$ is real, symmetric and positive definite. In the numerical analysis of this manuscript we consider only the case where $A$ is an interval made by $\ell$ sites, hence $\gamma_{A}$ is a $2\ell\times 2\ell$ matrix. Nonetheless, the considerations reported in this section hold for a generic number of spatial dimensions, once $\ell$ is understood as the number of sites in the subsystem $A$.

2.1 Entanglement hamiltonian

In the harmonic chain, the entanglement hamiltonian $\widehat{K}_{A}$ corresponding to a region $A$ introduced in (1) is a quadratic hermitian operator, hence it can be written as [23, 2]

[TABLE]

where the $2\ell$ dimensional vector $\boldsymbol{\hat{r}}$ collects the position and the momentum operators $\hat{q}_{i}$ and $\hat{p}_{i}$ with $i\in A$. The matrix $H_{A}$, that fully determines the entanglement hamiltonian in (4), is real, symmetric and positive definite. This guarantees that $\widehat{K}_{A}$ is hermitian, like the operators $\hat{q}_{i}$ and $\hat{p}_{i}$. In the following we often call $H_{A}$ the entanglement hamiltonian matrix.

A very important tool employed to quantify the bipartite entanglement in harmonic lattices is the Williamson’s theorem [68, 69, 70]. It holds for any real, symmetric and positive matrix having even order, but in our analysis we employ it for the $2\ell\times 2\ell$ matrices $\gamma_{A}$ and $H_{A}$. Considering the reduced covariance matrix $\gamma_{A}$ first, the Williamson’s theorem guarantees that we can construct a symplectic matrix $W\in\textrm{Sp}(\ell)$ such that

[TABLE]

where $\mathcal{D}_{\textrm{d}}\equiv\mathcal{D}\oplus\mathcal{D}$ and the diagonal matrix $\mathcal{D}=\textrm{diag}(\sigma_{1},\dots,\sigma_{\ell})$ collects the symplectic eigenvalues $\sigma_{k}>0$ of $\gamma_{A}$. The symplectic eigenvalues are uniquely determined up to permutation and constitute the symplectic spectrum. We will refer to the r.h.s. of (5) as the Williamson’s decomposition of $\gamma_{A}$. We will choose a decreasing ordering for the symplectic eigenvalues.

The symplectic spectrum of $\gamma_{A}$ provides the entanglement entropy $S_{A}$ and the Rényi entropies $S_{A}^{(n)}$ as follows [5, 39, 67, 6, 25, 1]

[TABLE]

where $s(y)$ and $s_{n}(y)$ are the analytic functions given by

[TABLE]

and

[TABLE]

The parameter $n$ is an integer $n\geqslant 2$. Performing an analytic continuation of this integer parameter, the entanglement entropy $S_{A}$ can be obtained as $S_{A}^{(1)}\equiv\lim_{n\to 1}S_{A}^{(n)}$, that is the replica limit for the entanglement entropy. The scalar quantities $S_{A}^{(n)}$ with $n\geqslant 1$ are usually called entanglement entropies, assuming that $S_{A}^{(1)}\equiv S_{A}$.

The symplectic spectrum of $\gamma_{A}$ is obtained by diagonalising $(\textrm{i}J\gamma_{A})^{2}$ or $(\textrm{i}\gamma_{A}J)^{2}$. Indeed, by employing that $W$ is symplectic and that $J^{\textrm{t}}\mathcal{D}^{r}_{\textrm{\tiny d}}J=\mathcal{D}^{r}_{\textrm{\tiny d}}$ for any non negative integer $r$, it is not difficult to realise that

[TABLE]

where we have introduced

[TABLE]

Being $W\,\widetilde{W}^{\textrm{t}}=\widetilde{W}^{\textrm{t}}\,W=\boldsymbol{1}$, the matrix $W$ is not orthogonal.

The Williamson’s theorem can be employed to decompose also the entanglement hamiltonian matrix $H_{A}$ defined in (4), namely

[TABLE]

where $W_{H}\in\textrm{Sp}(\ell)$ and $\mathcal{E}_{\textrm{d}}\equiv\mathcal{E}\oplus\mathcal{E}$, being $\mathcal{E}=\textrm{diag}(\varepsilon_{1},\dots,\varepsilon_{\ell})$ the symplectic spectrum of $H_{A}$. The symplectic eigenvalues $\varepsilon_{k}$ are often called single particle entanglement energies.

The symplectic spectrum of $H_{A}$ is related to the symplectic spectrum of $\gamma_{A}$ as follows

[TABLE]

Being $\mathcal{E}$ and $\mathcal{D}$ diagonal matrices, we have $\varepsilon_{k}=2\,\textrm{arccoth}(2\sigma_{k})=\log[(\sigma_{k}+1/2)/(\sigma_{k}-1/2)]$, whose inverse reads $\sigma_{k}=\coth(\varepsilon_{k}/2)/2=[(e^{\varepsilon_{k}}+1)/(e^{\varepsilon_{k}}-1)]/2$ for $1\leqslant k\leqslant\ell$. Let us remark that the uncertainty principle leads to $\sigma_{k}\geqslant 1/2$, and this implies that $\varepsilon_{k}>0$.

The determinant of $H_{A}$ is determined only by the symplectic spectrum. Indeed, by employing the Williamson’s decomposition (11) and the fact that a symplectic matrix has determinant equal to one, one finds $\textrm{det}\,H_{A}=\textrm{det}\,\mathcal{E}_{\textrm{\tiny d}}=\prod_{k=1}^{\ell}\varepsilon_{k}^{2}$.

In [28] the entanglement hamiltonian matrix $H_{A}$ has been expressed in terms of the reduced covariance matrix $\gamma_{A}$ as follows

[TABLE]

This relation can be obtained from (11) by employing (12) and

[TABLE]

where $\widetilde{W}$ is given by (10). The details of this derivation have been reported in the Appendix A. Notice that, if $B=A^{-1}(\Lambda\oplus\Lambda)A$ for some diagonal matrix $\Lambda$, then also $B=(J^{\textrm{t}}A)^{-1}(\Lambda\oplus\Lambda)(J^{\textrm{t}}A)$. This implies that we can replace $W$ with $J^{\textrm{t}}W$, $-W$ or $-J^{\textrm{t}}W$ in (9), and similarly for $\widetilde{W}$. From (5) and (11), one observes that $\gamma_{A}$ and $H_{A}$ are not affected by this ambiguity.

The expression (13) can be also written as

[TABLE]

where we have introduced the function $h(y)\equiv y^{-1}\log[(y+1/2)/(y-1/2)]$. The equivalence between (15) and (13) is shown in Appendix A. Notice that the expression (13) or its equivalent form (15) cannot be applied for pure states, where $(2\textrm{i}J\gamma)^{2}=\boldsymbol{1}$.

By employing (5), (10), (11) and (14), one finds \gamma_{A}\,H_{A}=W^{\textrm{t}}\big{(}\mathcal{D}_{\textrm{\tiny d}}\,\mathcal{E}_{\textrm{\tiny d}}\big{)}\widetilde{W} and its transpose H_{A}\,\gamma_{A}=\widetilde{W}^{\textrm{t}}\big{(}\mathcal{D}_{\textrm{\tiny d}}\,\mathcal{E}_{\textrm{\tiny d}}\big{)}W. This implies $\textrm{Tr}((H_{A}\gamma_{A})^{n})=2\sum_{k}(\sigma_{k}\,\varepsilon_{k})^{n}$ for $n\geqslant 1$. Thus, for a generic function $f(x)$ we have $\textrm{Tr}(f(H_{A}\,\gamma_{A}))=\textrm{Tr}(f(\gamma_{A}\,H_{A}))=2\sum_{k}f(\sigma_{k}\,\varepsilon_{k})$.

The above discussion tells that the entanglement hamiltonian matrix $H_{A}$ can be obtained from $\gamma_{A}$ either through the explicit expressions in (11) and (13) or by finding the symplectic spectrum of the reduced covariance matrix $\gamma_{A}$ and the symplectic matrix $W$ first and then use (11), (12) and (14). In the latter procedure, the needed data can be obtained through a standard diagonalisation, as indicates in (9). In our numerical analysis, this step is not straightforward because the matrix $(\textrm{i}J\gamma_{A})^{2}$ to diagonalise is not symmetric. In the Appendix B we describe the procedure that we have employed to construct the symplectic matrix $W$, that is based on the Cholesky decomposition of $\gamma_{A}$.

The Williamson’s theorem (11) and the relation (14) allow to write entanglement hamiltonian (4) in the following diagonal form

[TABLE]

where the hermitian operators $\hat{\mathfrak{q}}_{k}$ and $\hat{\mathfrak{p}}_{k}$ collected in the vector $\boldsymbol{s}$ satisfy the canonical commutation relations, being $\widetilde{W}$ symplectic. It is convenient to introduce the annihilation and creation operators $\boldsymbol{\hat{b}}^{\textrm{t}}\equiv(\hat{\mathfrak{b}}_{1}\,\dots\,\hat{\mathfrak{b}}_{\ell}\;\,\hat{\mathfrak{b}}_{1}^{\dagger}\,\dots\,\hat{\mathfrak{b}}_{\ell}^{\dagger}\,)$ as follows

[TABLE]

where the $2\ell\times 2\ell$ matrix $\Omega$ is unitary. Since $\Omega^{-1}(\textrm{i}J)\Omega^{-\textrm{t}}=J$, it is straightforward to check that $[\hat{b}_{i},\hat{b}_{j}]=J_{ij}$. In terms of these operators, the entanglement hamiltonian (16) becomes

[TABLE]

This operator has a well known form; hence it can be treated in the standard way. By introducing the eigenstates $|n_{k}\rangle$ of the occupation number operator $\mathfrak{b}_{k}^{\dagger}\mathfrak{b}_{k}$, whose eigenvalues are given by non negative integers $n_{k}$, the reduced density matrix (1) can be written through the projectors $|\boldsymbol{n}\rangle\equiv\bigotimes_{k=1}^{\ell}|n_{k}\rangle$ as follows [71]

[TABLE]

where $n_{k}$ is the $k$-th element of the $\ell$-dimensional vector $\boldsymbol{n}$. The coefficients $\lambda_{\boldsymbol{n}}$ provide the entanglement spectrum. The normalisation condition ${\rm Tr}\rho_{A}=1$ leads to

[TABLE]

where the orthonormality of the states $|\boldsymbol{n}\rangle$ has been employed to write $\mathcal{Z}_{A}$ as the product of $\ell$ geometric series (one for each $n_{k}$). Combining (19) and (20), the generic element of the entanglement spectrum can be expressed in terms of the symplectic spectrum $\mathcal{D}$ as follows

[TABLE]

The largest eigenvalue is given by $\lambda_{\text{\tiny max}}=\lambda_{\boldsymbol{0}}$, namely it has $n_{k}=0$ for all $k$. Its explicit expression reads $\lambda_{\text{\tiny max}}=\prod_{k=1}^{\ell}(1-e^{-\varepsilon_{k}})=\prod_{k=1}^{\ell}(\sigma_{k}+1/2)^{-1}$. Thus, (21) can be written as

[TABLE]

The gap between two eigenvalues in the entanglement spectrum characterised by the vectors $\boldsymbol{m}$ and $\boldsymbol{n}$ is given by

[TABLE]

The gaps with respect to the largest eigenvalue are $\log\lambda_{\text{\tiny max}}-\log\lambda_{\boldsymbol{n}}=\sum_{k=1}^{\ell}\varepsilon_{k}\,n_{k}$ and correspond to the special case $\boldsymbol{m}=\boldsymbol{0}$ in (23). By arranging the symplectic spectrum in decreasing order $\sigma_{1}\geqslant\sigma_{2}\geqslant\dots\geqslant\sigma_{\ell}$, we have $\varepsilon_{1}\leqslant\varepsilon_{2}\leqslant\dots\leqslant\varepsilon_{\ell}$ for the single particle entanglement energies. In order to find the gaps $0<g_{1}<g_{2}<\dots$ with respect to the largest eigenvalue, one can first introduce \mathcal{B}_{p}^{(r)}\equiv\big{\{}\sum_{j=1}^{p}\varepsilon_{k_{j}}\,,1\leqslant k_{j}\leqslant r\big{\}} for $r\geqslant 1$ and then compute the $r$-th gap as $g_{r}=\textrm{min}[(\cup_{p=1}^{r}\mathcal{B}_{p}^{(r)})\setminus\{g_{j},1\leqslant j\leqslant r-1\}]$. This procedure is not optimal, but it allows to get $g_{r}$ easily also for high values of $r$. The smallest gaps are given by $g_{1}=\varepsilon_{1}$, $g_{2}=\textrm{min}\{2\varepsilon_{1},\varepsilon_{2}\}$ and g_{3}=\textrm{min}\big{\{}3\varepsilon_{1},\textrm{max}\{2\varepsilon_{1},\varepsilon_{2}\},\varepsilon_{3}\big{\}}. The curves in Fig. 7 have been found by employing the above expressions.

It is useful to partition the symmetric matrices $\gamma_{A}$ and $H_{A}$ into $\ell\times\ell$ blocks as follows

[TABLE]

The blocks of $\gamma_{A}$ are the correlators $Q_{ij}=\langle\hat{q}_{i}\,\hat{q}_{j}\rangle$, $P_{ij}=\langle\hat{p}_{i}\,\hat{p}_{j}\rangle$ and $R_{ij}=\textrm{Re}\langle\hat{q}_{i}\,\hat{p}_{j}\rangle$, with $1\leqslant i\leqslant\ell$ and $1\leqslant j\leqslant\ell$. Correspondingly, the partitions in $\ell\times\ell$ blocks of the symplectic matrices $W$ and $\widetilde{W}$ occurring in the Williamson’s decompositions of $\gamma_{A}$ and $H_{A}$ (see (5), (11) and (10)) read

[TABLE]

and therefore the relation (10) becomes

[TABLE]

These partitions in $\ell\times\ell$ blocks and the relation (14) lead to write the Williamson’s decompositions (5) and (11) of $\gamma_{A}$ and $H_{A}$ respectively as

[TABLE]

and

[TABLE]

Since the symplectic spectrum $\mathcal{E}$ of $H_{A}$ is related to the symplectic spectrum $\mathcal{D}$ of $\gamma_{A}$ through (12), the expression (28) provides the entanglement hamiltonian matrix $H_{A}$ in terms of the blocks of the matrices occurring in the Williamson’s decomposition of $\gamma_{A}$.

In the Appendix C we discuss the special case where $\gamma_{A}$ and $H_{A}$ are block diagonal. This class is very important because it includes the entanglement hamiltonians matrices corresponding to static configurations [2, 29, 30, 33].

2.2 A contour for the entanglement entropies

The contour for the entanglement entropies is a function of the position inside the region $A$ (and of the time, when the system is out of equilibrium) that has been introduced to understand the spatial structure of the bipartite entanglement in pure states. For free lattice models, it has been studied in [39, 40, 41, 42, 16]. Like the entanglement hamiltonian, also the contour for the entanglement entropies cannot be determined from the entanglement spectrum only.

The contour function for the entanglement entropies $s_{A}^{(n)}:A\to\mathbb{R}$ must satisfy the following two conditions

[TABLE]

In [40] other three requirements have been introduced and a contour function for the entanglement entropy of free fermions has been constructed, verifying that it fulfils all these five properties. This set of constraints does not fix uniquely the contour function.

In the harmonic lattices, a contour function has been proposed in [42], that satisfies the two conditions in (29) and a weakened version of the other three conditions proposed in [40]. In this manuscript we will employ this prescription. Another expression fulfilling (29) has been suggested in [42], but the remaining three properties have not been proved for this proposal. A third contour function has been considered in [39, 41]: it satisfies the first condition in (29) but numerical violations of the positivity condition $s^{(n)}_{A}(i)\geqslant 0$ have been observed for some configurations [42].

Let us consider the case where $A$ is an interval made by $\ell$ sites, for simplicity. A contour function $s^{(n)}_{A}(i)$ can be constructed by associating $\ell$ real numbers $p_{k}(i)$ to every symplectic eigenvalue $\sigma_{k}$, where $1\leqslant i,k\leqslant\ell$. The function $p_{k}(i)$ is often called mode participation function [39] and it fulfils the following conditions

[TABLE]

A mode participation function $p_{k}(i)$ allows to write the entanglement entropies (6) as in (29) with

[TABLE]

being $s_{n}(y)$ the functions defined in (7) and (8) for $n\geqslant 1$.

We find it natural to consider also the contribution to the entanglement entropies given by an interval $(i_{1},i_{2})\subset A$, namely

[TABLE]

The entanglement entropies $S^{(n)}_{A}$ correspond to the special case where $i_{1}$ and $i_{2}$ are the sites at the boundary of $A$.

We are interested in the temporal evolution of the above quantities after a global quench and it is often useful to study their increments with respect to their initial value. In particular, for $\mathcal{S}^{(n)}_{A}(i_{1},i_{2};t)$ given by (32), let us introduce

[TABLE]

In harmonic lattices, the contour function (31) is not unique because any mode participation function fulfilling (30) provides a contour function (31) that satisfies (29). Furthermore, any $2\ell\times 2\ell$ orthogonal matrix $O$ can be employed to construct a mode participation function $p_{k}(i)$. In order to make contact with the reduced covariance matrix $\gamma_{A}$, in [42] the Euler decomposition of the symplectic matrix $W$ in (5) has been employed to construct an orthogonal matrix $K$ leading to a reasonable candidate for the contour function $s_{A}^{(n)}(i)$. This orthogonal matrix can be written in terms of $W$ as follows

[TABLE]

Partitioning this $2\ell\times 2\ell$ orthogonal matrix into $\ell\times\ell$ blocks

[TABLE]

the contour function for the entanglement entropies in harmonic lattices proposed in [42] is (31) with the mode participation function given by

[TABLE]

Thus, the symplectic spectrum of $\gamma_{A}$ and the symplectic matrix $W$ occurring in its Williamson’s decomposition (5) provide both the entanglement hamiltonian matrix $H_{A}$ of [28] (see e.g. (28)) and the contour function for the entanglement entropies of [42].

Let us remark that, although the discussion about the contour function $s_{A}^{(n)}(i)$ has been done for an interval, it is straightforward to extend the above formulas to higher spatial dimensions or to subregions $A$ made by many disconnected components.

3 Free fermionic lattices

In this section we consider a fermionic system described by a quadratic hamiltonian. In §3.1 we study the entanglement hamiltonian of a subsystem, that is quadratic as well. In §3.2 we review the construction of the contour for the entanglement entropy in these systems proposed in [40]. Relevant special cases are discussed in §3.3.

3.1 Entanglement hamiltonian

In a chain of free fermions described by a quadratic hamiltonian and in a Gaussian state, the entanglement hamiltonian $\widehat{K}_{A}$ takes the following quadratic form [25]

[TABLE]

where for the operators $\hat{c}_{i}$ and $\hat{c}_{i}^{\dagger}$ the canonical anticommutation relations $\{\hat{c}_{i},\hat{c}_{j}^{\dagger}\}=\delta_{i,j}$ and $\{\hat{c}_{i},\hat{c}_{j}\}=\{\hat{c}_{i}^{\dagger},\hat{c}_{j}^{\dagger}\}=0$ hold. The complex matrices in (37) satify $D=-B^{\ast}$ and $G=-F^{\ast}$ to ensure that $\widehat{K}_{A}$ is hermitian. Collecting the operators into the $2\ell$ dimensional vector $\boldsymbol{\hat{a}}$ defined by $\boldsymbol{\hat{a}}^{\textrm{t}}\equiv(\hat{c}_{1}\,\dots\,\hat{c}_{\ell}\;\,\hat{c}_{1}^{\dagger}\,\dots\,\hat{c}_{\ell}^{\dagger}\,)$, the canonical anticommutation relations become $\{\hat{a}_{i},\hat{a}_{j}\}=I_{i,j}$, being $I$ the $2\ell\times 2\ell$ symmetric matrix with the $\ell\times\ell$ matrix $\boldsymbol{0}$ along the diagonal and the $\ell\times\ell$ identity matrix off diagonal. Standard manipulations allow to write the entanglement hamiltonian (37) in the following form [25]

[TABLE]

up to constant terms, where the $\ell\times\ell$ matrices $T$ and $X$ are

[TABLE]

and satisfy the conditions $T^{\dagger}=T$ and $X^{\textrm{t}}=-X$.

The matrix $\Omega$ defined in (17) allows to introduce the vector $\boldsymbol{\hat{r}}\equiv\Omega\,\boldsymbol{\hat{a}}$, whose elements are hermitian operators that satisfy $\{\hat{r}_{i},\hat{r}_{j}\}=\delta_{i,j}$, being $\Omega\,I\,\Omega^{\textrm{t}}=\boldsymbol{1}$. In terms of these Majorana operators, the entanglement hamiltonian (38) becomes

[TABLE]

where $H_{A}$ is the following matrix

[TABLE]

Since $H_{A}$ is $2\ell\times 2\ell$ real and antisymmetric, a $2\ell\times 2\ell$ real orthogonal matrix $\widetilde{O}$ exists such that

[TABLE]

where $\mathcal{E}=\textrm{diag}\{\varepsilon_{1},\dots,\varepsilon_{\ell}\}$ with $\varepsilon_{k}>0$ [72]. The $\varepsilon_{k}$ are sometimes called single particle entanglement energies. From (42) one observes that $\widetilde{O}\,(\textrm{i}H_{A})^{2}\,\widetilde{O}^{\textrm{t}}=\mathcal{E}^{2}\oplus\mathcal{E}^{2}$. Thus, the real symmetric matrix $-H_{A}^{2}$ has a degenerate spectrum given by $\{\,\varepsilon_{k}^{2}\,,1\leqslant k\leqslant\ell\,\}$ and is diagonalised by the orthogonal matrix $\widetilde{O}$.

The second crucial $2\ell\times 2\ell$ real antisymmetric matrix in the analysis of the fermionic Gaussian state is the following correlation matrix [73, 75, 74, 76]

[TABLE]

which gives $\langle\hat{r}_{i}\hat{r}_{j}\rangle=[\,\textrm{i}\,(\Gamma_{A})_{i,j}+\delta_{i,j}\,]/2$. The same theorem employed above for $H_{A}$ allows to claim that a $2\ell\times 2\ell$ real orthogonal matrix $O$ exists such that

[TABLE]

where \mathcal{N}\equiv\textrm{diag}\big{\{}\nu_{1},\dots,\nu_{\ell}\big{\}} with $\nu_{k}>0$.

The single particle entanglement spectrum $\mathcal{E}$ and the orthogonal matrix $\widetilde{O}$ in (42) are related to the diagonal matrix $\mathcal{E}$ and the orthogonal matrix $O$ in (44) as follows [25]

[TABLE]

where the first expression is equivalent to $\varepsilon_{k}=2\,\textrm{arctanh}(\nu_{k})=\log[(1+\nu_{k})/(1-\nu_{k})]$ with $1\leqslant k\leqslant\ell$. By employing (44), the relations in (45) lead to111Notice that, for an odd function $f_{\textrm{\tiny odd}}(x)\equiv\sum_{p\geqslant 0}t_{p}\,x^{2p+1}$, by using (44), one finds

f_{\textrm{\tiny odd}}\big{(}O\,\Gamma_{A}\,O^{\textrm{t}}\big{)}=\bigg{(}\!\!\begin{array}[]{cc}\boldsymbol{0}\!&\tilde{f}_{\textrm{\tiny odd}}(\mathcal{N})\\ -\,\tilde{f}_{\textrm{\tiny odd}}(\mathcal{N})\!&\boldsymbol{0}\\ \end{array}\!\bigg{)}

(46) where $\tilde{f}_{\textrm{\tiny odd}}(x)\equiv\sum_{p\geqslant 0}(-1)^{p}\,t_{p}\,x^{2p+1}$.

[TABLE]

This formula plays the same role of (13) in the bosonic case (see §2.1).

The relation (42) leads to introduce $\boldsymbol{\hat{d}}\equiv O\,\boldsymbol{\hat{r}}$, whose elements satisfy $\{\hat{d}_{i},\hat{d}_{j}\}=\delta_{i,j}$, being $O$ orthogonal. Furthermore, it is convenient to define fermionic creation and annihilation operators $\hat{\mathfrak{f}}_{k}^{\dagger}$ and $\hat{\mathfrak{f}}_{k}$ as $\boldsymbol{\hat{f}}\equiv\Omega^{-1}\boldsymbol{\hat{d}}$, where $\boldsymbol{\hat{f}}^{\textrm{t}}=(\,\hat{\mathfrak{f}}_{1}\,\dots\,\hat{\mathfrak{f}}_{\ell}\;\,\hat{\mathfrak{f}}_{1}^{\dagger}\,\dots\,\hat{\mathfrak{f}}_{\ell}^{\dagger}\,)$, whose elements satisfy $\{\hat{f}_{i},\hat{f}_{j}\}=I_{i,j}$. Combining the linear maps introduced above, one obtains

[TABLE]

$\Omega^{-1}\,O\,\Omega$ is unitary, being $O$ real orthogonal and $\Omega$ unitary. The relation (42) and the anticommutation relations of $\boldsymbol{\hat{f}}$ allow to write the entanglement hamiltonian (40) in terms of the operators in $\boldsymbol{\hat{d}}$ and $\boldsymbol{\hat{\mathfrak{f}}}$ respectively as

[TABLE]

that makes explicit the role of the single particle entanglement energies.

Let us introduce the partition in $\ell\times\ell$ blocks of the orthogonal matrix defined in (44), namely

[TABLE]

where the blocks are constrained by the orthogonality of $O$ as follows

[TABLE]

In terms of the blocks of $O$ in (50), the unitary matrix in (48) becomes

[TABLE]

By employing that $\boldsymbol{\hat{d}}\equiv O\,\Omega\,\boldsymbol{\hat{a}}$ in (49) and comparing the resulting expression with (38), one writes $T$ and $X$ in terms of the blocks of $O$ in (50) as follows

[TABLE]

The expression (49) for the entanglement hamiltonian naturally leads to introduce the fermionic occupation number operators $\hat{\mathfrak{f}}_{k}^{\dagger}\,\hat{\mathfrak{f}}_{k}$, whose eigenvalues and eigenvectors allow to write the reduced density matrix as $\rho_{A}=\sum_{\boldsymbol{n}}\lambda_{\boldsymbol{n}}|\boldsymbol{n}\rangle\langle\boldsymbol{n}|$. Thus, one can perform an analysis similar to the one reported for the harmonic lattices in §2.1, with the crucial difference that the elements of the vector $\boldsymbol{n}$ are $n_{k}\in\{0,1\}$ in these fermionic models. In particular, the reduced density matrix can be written as

[TABLE]

whose normalisation condition ${\rm Tr}\rho_{A}=1$ provides the following normalisation constant

[TABLE]

Plugging this expression in (56), one obtains the entanglement spectrum

[TABLE]

where the largest eigenvalue corresponds to $n_{k}=0$ for all $k$, being $\varepsilon_{k}>0$. This gives $\lambda_{\text{\tiny max}}=\prod_{k=1}^{\ell}(1+e^{-\varepsilon_{k}})^{-1}$, that has been used to get the last expression in (57). The gap between two eigenvalues can be written as follows

[TABLE]

We find it convenient to consider the gaps with respect to the largest eigenvalue. These gaps and their temporal evolution after a global quench are discussed in §6.

3.2 A contour for the entanglement entropies

An exhaustive analysis of a contour function for the entanglement entropy in chains of free fermions has been carried out in [40]. We employ this prescription in our analysis and refer the reader interested in further details to this reference.

From (44) and the matrix $\Omega$ introduced in (17), it is straightforward to find

[TABLE]

Similarly, considering (42) and (45), one also obtains

[TABLE]

Thus, by employing (59), for a generic function $f$, we have [5, 40]

[TABLE]

The entanglement entropy can be obtained by specifying this formula to $f(x)=-\,x\log(x)$, while the $n$-th moment ${\rm Tr}\rho_{A}^{n}$ of the reduced density matrix corresponds to $f(x)=x^{n}$. We remark that $\nu_{k}<1$ [75].

A natural contour for the entanglement entropies is obtained as explained in §2.2 (see (31)) with the mode participation function given by [40]

[TABLE]

In [40] it has been shown that the contour function constructed from (62) satisfies other three properties beside (29) that we do not discuss in this manuscript. All these constraints can be part of a more detailed definition (still unknown) that could lead to identify the contour function for the entanglement entropies in a unique way.

3.3 Special cases

It is worth focussing the expressions discussed in §3.1 and §3.2 in some special cases.

First, let us consider the class of entanglement hamiltonians having

[TABLE]

that includes the entanglement hamiltonian after the global quench discussed in §6.

The unitary matrix (52) that provides the map (48) becomes the diagonal matrix $\Omega^{-1}\,O\,\Omega=(U_{O}+\textrm{i}\,Z_{O})\oplus(U_{O}-\textrm{i}\,Z_{O})$, with $U_{O}+\textrm{i}\,Z_{O}$ unitary. The latter condition can be obtained also by specialising the orthogonality condition (51) to this case. Furthermore, (3.1) and (3.1) simplify respectively to

[TABLE]

where we remind that $T$ is hermitian. Hence, (38) becomes a fermionic hamiltonian with hopping terms given by [25]

[TABLE]

and (41) simplifies to

[TABLE]

The hermitian matrix $T$ in (65) can be diagonalised by a unitary matrix $\widetilde{U}$ and it has real eigenvalues $\eta_{k}$, i.e. $\widetilde{U}\,T\,\widetilde{U}^{\dagger}=\textrm{diag}\{\eta_{i},\dots\eta_{\ell}\}$. This decomposition allows to write

[TABLE]

Plugging this result into (65), the entanglement hamiltonian becomes $\widehat{K}_{A}=\sum_{k=1}^{\ell}\eta_{k}\,\hat{\mathfrak{f}}^{\dagger}_{k}\,\hat{\mathfrak{f}}_{k}$, where we have defined $\boldsymbol{\hat{\mathfrak{f}}}=\widetilde{U}\,\boldsymbol{\hat{c}}$. Comparing this result with the map in (48), that is diagonal in this special case and it is written in the text below (63), we conclude that $\widetilde{U}=U_{O}+\textrm{i}\,Z_{O}$. These observations have been done in [25], where it has also found that

[TABLE]

being $C_{A}$ the correlation matrix whose generic element is $(C_{A})_{i,j}=\langle\hat{c}_{i}^{\dagger}\hat{c}_{j}\rangle$, with $i,j\in A$. From (68), we have that the eigenvalues $\eta_{k}$ of $T$ are related to the eigenvalues $\zeta_{k}$ of $C_{A}$ as follows

[TABLE]

Comparing with §3.1, one finds $\varepsilon_{k}=|\eta_{k}|$ and $\nu_{k}=|2\zeta_{k}-1|$.

The Rényi entropies are obtained from the eigenvalues of $C_{A}$ as [1, 2, 4, 5, 6]

[TABLE]

The limit $n\to 1^{+}$ of these expressions provides the entanglement entropy as

[TABLE]

As for the contour for the entanglement entropies, in the special case (63) the mode participation function (62) simplifies to

[TABLE]

As consistency check, let us observe that the $i$-th element of the diagonal of the matrix relation $\widetilde{U}^{\dagger}\widetilde{U}=\boldsymbol{1}$ gives the condition $\sum_{k=1}^{\ell}p_{k}(i)=1$ for $1\leqslant i\leqslant\ell$.

Another relevant special case corresponds to $Y_{O}=Z_{O}=\boldsymbol{0}$. For these entanglement hamiltonian matrices the unitary matrix (52) reduces to

[TABLE]

where $U_{O}$ and $V_{O}$ are orthogonal. In this case (3.1) and (3.1) become respectively

[TABLE]

Let us remarks that in the formulas described above we have not specified the explicit form of the two-point correlators, hence they can be used to study more complicated configurations like e.g. the ones where $A$ is the union of many disjoint blocks of sites.

Finally, an interesting class of entanglement hamiltonian matrices included in the one defined by (63) corresponds to $V_{O}=U_{O}$ and $Y_{O}=Z_{O}=\boldsymbol{0}$. In this case the unitary matrix (52) becomes the diagonal matrix $\Omega^{-1}\,O\,\Omega=U_{O}\oplus U_{O}$, where $U_{O}$ is orthogonal. Furthermore, the expression of $T$ in (64) simplifies to $T=U_{O}^{\textrm{t}}\,\mathcal{E}\,U_{O}$. For instance, at equilibrium, the entanglement hamiltonians of an interval in the ground state considered in [31, 32] belongs to this class.

4 Insights from CFT

In this section we employ the analytic results obtained for the entanglement hamiltonian of a semi-infinite line after a global quench in CFT [15] to get some insights about the qualitative behaviour of the entanglement spectrum and of the contour for the entanglement entropies.

The CFT analysis performed in [44] (see also the recent review [63]) leads to the following linear growth of the entanglement entropies before the saturation

[TABLE]

(comparing the notation with [63], we have $\tau_{0}=4\tau_{0,\textrm{\tiny there}}$), where $c$ is the central charge of the model. This gives $S_{A}\simeq 2\pi c\,t/(3\tau_{0})$ for the entanglement entropy and, by fitting this linear growth, we can get $\tau_{0}$ numerically. A slight dependence on $n$ has been observed in $\tau_{0}$ (see e.g. [55]). We find worth remarking that a factor of $2$ in (75) is due to the fact that the interval has two endpoints. Taking the limit $n\to\infty$ in (75) one obtains $-\log\lambda_{\textrm{\tiny max}}\simeq\pi c\,t/(3\tau_{0})$. The parameter $\tau_{0}$ encodes some features of the initial state (we refer to [47] for a complete discussion).

The temporal evolution of the entanglement hamiltonian of a semi-infinite line after a global quantum quench in CFT has been studied in [15] by employing methods and results of CFT with boundaries [19]. This analysis provides also the entanglement spectrum and, in particular, the above linear growth for $-\log\lambda_{\textrm{\tiny max}}$ is recovered up to a factor of $2$, which is due to the fact that the semi-infinite line has only one endpoint. This result suggests to explore the qualitative behaviour of the entanglement spectrum before the saturation at $t/\ell\simeq 1/2$ by taking the analytic CFT expressions for the semi-infinite line and introducing properly a factor of $2$ to take into account the occurrence of two endpoints. As for the gaps in the entanglement spectrum, from [15] we obtain

[TABLE]

which is $1/2$ of the corresponding quantity for the semi-infinite line. In (76), the coefficients $\Delta_{a}$ are the non vanishing conformal dimensions (including also the ones of the descendants) of a boundary CFT with the proper boundary condition, as discussed in [15]. We find it worth considering also the ratios $g_{b,0}/g_{a,0}=\Delta_{b}/\Delta_{a}$ because they are independent of $\tau_{0}$. From (75) and (76) it is straightforward to construct another expression where $\tau_{0}$ does not occur, namely

[TABLE]

The temporal evolutions of (76), of (77) and of the ratios of the entanglement gaps with respect to the first one after the global quenches that we consider in this manuscript are shown in Fig. 7 and Fig. 20. From these data we notice that the best agreement with CFT is observed for the ratios of the entanglement gaps.

The above analysis can be easily adapted to the temporal evolution after a local quantum quenches [48] by employing the corresponding results of [15]. The crucial difference is that logarithmic growths occur in these cases. For instance, for local quenches (76) has $\log t$ instead of $t$ in the denominator. Instead, (77) holds also for local quenches.

In the remaining part of this section we exploit the CFT results for the entanglement hamiltonian of the semi-infinite line after a global quench [15] to get insights about the temporal evolution of the contour function for the entanglement entropies of an interval.

For some particular one-dimensional spatial bipartitions at equilibrium (e.g. an interval in the infinite line at zero and finite temperature, or in a finite periodic system for the ground state) a natural candidate for the contour function for the entanglement entropies in CFT can be obtained from the weight function occurring in the entanglement hamiltonian [42]. By applying this idea to the half line after a global quench, one finds the following CFT ansatz for the contour function [15]

[TABLE]

where $n\geqslant 1$ and

[TABLE]

It is straightforward to observe that $\tau_{0}\,s^{(n)}_{\textrm{h.l.}}(x,t)$ depends on $x/\tau_{0}$, $t/\tau_{0}$ and the central charge $c$. Furthermore, $\mathcal{F}_{\textrm{h.l.}}(x,t)=1/x+O(1)$ as $x\to 0^{+}$, independently of time. This implies that the expected logarithmic divergence of $S_{A}^{(n)}$ as the UV cutoff $\epsilon\to 0$ is independent of time as well. Indeed, this divergence can be obtained by removing an infinitesimal disk around the endpoint of the half line and integrating the contour function (78) from $\epsilon$ to some point inside the half line. Instead, by integrating (78) in the entire half line $(\epsilon,+\infty)$ with $\epsilon/\tau_{0}\ll 1$, one obtains [15]

[TABLE]

that gives the linear growth found in [44], namely $\tfrac{\pi c}{3\tau_{0}}t$ when $t\gg\tau_{0}$. The r.h.s. of (80) provides only the leading term of the entanglement entropy as $\epsilon\to 0$. A subleading $O(1)$ term comes from the contribution of the conformal boundary states introduced through the regularisation procedure [15, 77] (see [42] for a discussion on this term in relation to the contour function at equilibrium).

Considering $t/\tau_{0}$ fixed in (78), we find that $s^{(n)}_{\textrm{h.l.}}(x,t)\sim e^{-2\pi x/\tau_{0}}$ as $x\to+\infty$ and that the height of the horizontal plateau is given by $\tfrac{\pi\,c}{12\tau_{0}}(1+\tfrac{1}{n})$ (this can be done by taking $t\to\infty$ first and then $x\to+\infty$).

A typical temporal evolution of the contour function for the entanglement entropy $s_{\textrm{h.l.}}(x,t)\equiv s^{(1)}_{\textrm{h.l.}}(x,t)$ is shown in the left panel of Fig. 1. The linear divergence as $x\to 0^{+}$ and the exponential decay as $x\to+\infty$ are independent of time. These two regimes are connected for the intermediate values of $x$ by a smooth step at $x\simeq 2t$ whose height is $\tfrac{\pi\,c}{6\,\tau_{0}}$. The occurrence of this step moving in time has been first observed in the numerical analysis performed in [40].

In this manuscript we are interested in the contour function for the entanglement entropies of an interval $A$ of length $\ell$ in the infinite line after a global quantum quench. A CFT analysis for this case is not available in the literature. Hence, let us consider the naive contour function obtained by superposing the two contour functions (78) associated to the half lines $(0,+\infty)$ and $(-\infty,\ell)$, namely

[TABLE]

with $\mathcal{F}_{\textrm{h.l.}}(x,t)$ given by (79). The observations made above about $\mathcal{F}_{\textrm{h.l.}}(x,t)$ lead to notice that $\tau_{0}\,s^{(n)}_{A,0}(x,t)$ obtained from (81) is a function of $t/\ell$, $x/\ell$ and $\tau_{0}/\ell$; and that $s^{(n)}_{A,0}(x,t)$ displays a linear divergence for $x\to 0^{+}$ and for $x\to\ell^{-}$, with the same time independent coefficient. For $t=0$, the naive analytic result (81) gives

[TABLE]

and in the asymptotic regime of long time we find

[TABLE]

The naive contour function for the entanglement entropy is given by (81) with $n=1$ and its typical temporal evolution is shown in the right panel of Fig. 1. This evolution can be explained qualitatively in terms of two fronts with height $\tfrac{\pi c}{6\tau_{0}}$ that start at the two different endpoints of the interval and then propagate in opposite directions with velocity equal to $2$. The two fronts cross and superpose at $t/\ell\simeq 1/4$ and each front takes $t\simeq\ell/2$ to travel across the entire interval. This qualitative features have been first observed in [40] by analysing the numerical data obtained in a chain of free fermions. Further numerical data supporting this picture are presented in §5 and §6.

Given two points $0<x_{1}<x_{2}<\ell$ inside the interval $A$, let us integrate the naive contour function (81) in the interval $(x_{1},x_{2})$. The result reads

[TABLE]

which is a function of $t/\ell$, $x_{1}/\ell$, $x_{2}/\ell$ and $\tau_{0}/\ell$. Subtracting to (84) its value at $t=0$, one finds

[TABLE]

where we have introduced the following function

[TABLE]

In Fig. 11 and Fig. 22 this formula (where the value of $\tau_{0}$ is obtained by fitting the linear growth of the entanglement entropy) gives the grey dashed-dotted curves, that can be compared against numerical data corresponding to particular quenches in a harmonic chain (see §5) and in a chain of free fermions (see §6) respectively.

The expression in (85) provides the entanglement entropies when $(x_{1},x_{2})=(\epsilon,\ell-\epsilon)$. In this case, we have

[TABLE]

Taking $t\gg\tau_{0}$ and $\ell\gg\tau_{0}$ first, and then considering the regimes $t\ll\ell/2$ and $t\gg\ell/2$, this expression simplifies respectively to

[TABLE]

in agreement with [44].

Let us remark that the results reported above have been obtained by employing in a naive way the CFT analytic expressions of [15] corresponding to the half line to study the case of a finite interval, assuming that the dynamics before the thermalisation at $t/\ell\simeq 1/2$ is governed by a neighbourhood of the endpoints, where the interval is indistinguishable from the half line. It would be interesting to extend the analysis of [15] for the quantum quenches to a finite interval.

5 Interval in a harmonic chain

In this section we study the temporal evolution of the entanglement hamiltonian matrix $H_{A}$ of an interval $A$ in an infinite harmonic chain after a global quench of the frequency parameter.

The hamiltonian (2) with periodic boundary conditions $\hat{q}_{L}=\hat{q}_{0}$ and $\hat{p}_{L}=\hat{p}_{0}$ can be diagonalised in the Fourier space by introducing the annihilation and creation operators $a_{k}$ and $a_{k}^{\dagger}$ in the standard way. This procedure can be understood also in terms of canonical transformations through the Williamson’s theorem applied to the matrix defining the quadratic hamiltonian (2). The diagonalised form of (2) is given by

[TABLE]

where the dispersion relation reads

[TABLE]

The quench protocol is the following. The system is prepared in the ground state $|\psi_{0}\rangle$ of the hamiltonian (2) with periodic boundary conditions and non vanishing mass $\omega=\omega_{0}\neq 0$, whose diagonalised form is

[TABLE]

being the dispersion relation $\omega_{0,k}$ given by (90) with $\omega=\omega_{0}$. At $t=0$ the frequency parameter is suddenly quenched to a different value $\omega\neq\omega_{0}$; hence the unitary time evolution of $|\psi_{0}\rangle$ for $t>0$ is determined by the hamiltonian (2) with $\omega\neq\omega_{0}$, namely

[TABLE]

The temporal evolution of several quantities after this quench have been studied in various papers (see the reviews [63, 64] for an exhaustive list of references).

We are interested in the temporal evolution of the entanglement hamiltonian matrix $H_{A}$ discussed in §2.1 after the global quench of the frequency parameter described above. The entanglement hamiltonian matrix $H_{A}(t)$ can be evaluated once the reduced covariance matrix $\gamma_{A}(t)$ is known (see (13), (15) and (28)). Thus, the first step is to construct the covariance matrix $\gamma(t)$ of the entire system, whose blocks are given by the following two-point correlators

[TABLE]

where the time evolved operators in the Heisenberg picture $\hat{q}_{j}(t)$ and $\hat{q}_{j}(t)$ read

[TABLE]

For periodic boundary conditions, the explicit expressions of the correlators in (93) in the thermodynamic limit read [46]

[TABLE]

where the functions $Q_{k}(t)$, $P_{k}(t)$ and $R_{k}(t)$ containing the dependence on $\omega$, $\omega_{0}$ and $t$ are defined as follows

[TABLE]

where $\omega_{k}\equiv\sqrt{\omega^{2}+\tfrac{4\kappa}{m}[\sin(k/2)]^{2}}$ and $\omega_{0,k}$ is given by the same expression with $\omega$ replaced by $\omega_{0}$. At $t=0$ these functions become respectively

[TABLE]

Plugging these expressions into (95), one recovers the correlators associated to the ground state of $\widehat{H}_{0}$, as expected.

It is important to remark that for $t>0$ all the modes labelled by $k$ in (95) and (96) provide a finite contribution when $\omega=0$. The zero mode, i.e. the one corresponding to $k=0$, deserves a particular attention because, for static configurations, the occurrence of this mode leads to a divergence in the correlator $\langle\hat{q}_{i}\hat{q}_{j}\rangle$ as $\omega\to 0^{+}$. This is related to the translation invariance of the periodic chain. Instead, the correlators (95) are finite in the massless regime; hence, we are allowed to consider a global quench where the evolution hamiltonian has $\omega=0$.

The time evolved state (92) is also Gaussian, hence it can be described through the formalism discussed in §2.1. In particular, the correlators (95) provide the blocks of the covariance matrix $\gamma(t)$ after the global quench. The correlators (95) or their modifications have been employed to study different aspects of entanglement in harmonic lattices (see e.g. [52, 55, 59]).

5.1 Numerical results for the entanglement hamiltonian matrix

We are interested in the temporal evolution of the entanglement hamiltonian matrix $H_{A}(t)$ corresponding to an interval $A$ made by $\ell$ sites after the global quench of the frequency parameter. The reduced covariance matrix $\gamma_{A}(t)$ and the entanglement hamiltonian matrix $H_{A}(t)$ can be partitioned into $\ell\times\ell$ blocks (see (24)) that depend on $t>0$, namely

[TABLE]

The relation between $H_{A}(t)$ and the reduced covariance matrix $\gamma_{A}(t)$ can be expressed as in (13), (15) or (28).

In our numerical analysis, the entanglement hamiltonian matrix $H_{A}(t)$ has been constructed by employing the Williamson’s decomposition (11) for $H_{A}$ and (14), where the symplectic matrix $W$ entering in the Williamson’s decomposition (5) for $\gamma_{A}$ has been obtained as explained in the Appendix B. The symplectic eigenvalues of $H_{A}(t)$ have been found from the symplectic spectrum of $\gamma_{A}$ through the relation (12). The correlators employed to define $\gamma_{A}(t)$ have been introduced in (95). In this manuscript we mainly consider the global quench of the frequency parameter described above with $\omega_{0}=1$ and $\omega=0$. Some data corresponding to an initial state having $\omega_{0}=5$ are also provided.

The relation (12), that provides the symplectic spectrum of $H_{A}$ in terms of the one of $\gamma_{A}$, requires $\sigma_{k}>1/2$. Since many elements of the symplectic spectrum of $\gamma_{A}$ are very close to $1/2$, the software approximates them to $1/2$, spoiling the applicability of (12) to get the symplectic spectrum of $H_{A}$. This forces to obtain $\sigma_{k}$’s with very high numerical precision. The number of digits depends on various parameters of the configuration like $\ell$ and $t$. In our numerical analysis for the harmonic chain we worked with precisions between $200$ and $1000$ digits.

In Fig. 2 we show the temporal evolution of the entanglement hamiltonian matrix $H_{A}(t)$ after the global quench written in the form (98), i.e. in terms of its $\ell\times\ell$ blocks. The main feature to highlight with respect to the entanglement hamiltonian matrices corresponding to static configurations [2] is that the off diagonal block $E$ is not vanishing and therefore it contributes in a non trivial way to determine the out of equilibrium dynamics. The magnitudes of the elements of $H_{A}(t)$ become very small for long time. Indeed, a smaller scale has been chosen in the three bottom panels of Fig. 2 to show the structure of the non vanishing elements. At early times only the diagonal blocks are non vanishing and mostly around their main diagonals. As time evolves, also the off diagonal elements of the blocks become relevant. The most evident feature that can be observed from the time evolution in Fig. 2 is the occurrence of bands in the different blocks of $H_{A}(t)$ whose widths increase with time. We consider these widths in more detail below, during the discussion of Fig. 4, Fig. 5 and Fig. 6.

The reflection symmetry of the configuration with respect to the center of the interval leads to a corresponding reflection symmetry with respect to the center in any given diagonal of the blocks $M(t)$ and $N(t)$ at any fixed time. More explicitly, we checked numerically that $M_{i,i}=M_{\ell-i+1,\ell-i+1}$ and $N_{i,i}=N_{\ell-i+1,\ell-i+1}$ with $1\leqslant i\leqslant\ell$ for the main diagonals and that $M_{i,i+p}=M_{\ell-p-i+1,\ell-i+1}$ and $N_{i,i+p}=N_{\ell-p-i+1,\ell-i+1}$ along the $p$-th diagonal, with $p>0$ and $1\leqslant i\leqslant\ell-p$. As for the off diagonal block $E(t)$, we observe a symmetry with respect to the center of the block, namely $E_{i,j}=E_{\ell-i+1,\ell-j+1}$ for $1\leqslant i,j\leqslant\ell$.

The temporal evolution of the blocks composing $H_{A}(t)$ is qualitatively different. As for the symmetric $\ell\times\ell$ matrix $M(t)$ in (98), the largest contributions come from the main diagonal, that decrease as time evolves and remain positive during all the evolution, reaching a stationary curve for long times (see the top left panel of Fig. 3). The first diagonals are mostly negative and in some point become positive at some time (see the top right panel of Fig. 3 and the corresponding inset). In the top left panel of Fig. 4 we show the antidiagonal of $M(t)$, that displays oscillations around zero.

Interesting features can be observed also for the temporal evolution of the symmetric block $N(t)$ in (98). The diagonal of $N(t)$ is shown in the middle left panel of Fig. 3: first it grows until a maximal curve and then it relaxes to a positive curve for long time. The elements of $N(t)$ are mostly positive but some of them become negative for large times (see e.g. the inset of the top right panel in Fig. 4 and the one of the middle right panel in Fig. 3). The block $N(t)$ is the only one where the evolution of the main diagonal and of the first diagonal are qualitatively similar: a smooth positive wedge grows from the center of the interval until $t/\ell\simeq 0.2$, then its tip decreases and gets smooth forming a plateau whose height decreases in time.

The evolution of the block $E(t)$ is qualitatively different from the one of $M(t)$ and $N(t)$ because this matrix is not symmetric. Furthermore, it vanishes for $t=0$ and it seems that its contribution is zero also for long times. The elements of this block are mostly negative but without a definite sign for all the times (see e.g. its first diagonal in the bottom right panel of Fig. 3). From Fig. 2, we observe the formation of two mostly negative bands close to the diagonal for small times having the shape of two parallelograms (see also the top right panel in Fig. 18) whose height along the antidiagonal increases in time until $t/\ell\simeq 0.5$, when they vanish. As for the main diagonal of $E(t)$, a smooth negative wedge develops until it reached a minimum value; then it goes back to zero in a peculiar way; indeed a small plateau is formed whose width increases and then decreases, leading to two fronts that move in the opposite directions towards the endpoints of the interval (see the bottom left panel of Fig. 3).

In Fig. 4 we show the temporal evolution of the antidiagonals of the blocks composing the entanglement hamiltonian matrix $H_{A}(t)$ in (98). Although the temporal evolutions of the antidiagonals in the three blocks look quite different, the main common feature is the fact that at a fixed time they become very small (and not vanishing) around the same points. The distance $2d_{0}$ between these points provides the total width of the bands occurring at any given time in each block (see Fig. 2). The width $d_{0}$ increases linearly in time with velocity equal to one, as shown in the top left panel of Fig. 18. This feature has been highlighted in Fig. 5, where each panel shows all the three antidiagonals together at a given time.

As for the block $M(t)$, its main diagonal (top left panel of Fig. 3) is not captured in the top left panel of Fig. 4 because $\ell$ is even. In $N(t)$ both the temporal evolutions of the diagonal (middle left panel in Fig. 3) and of the antidiagonal (top right panel in Fig. 4) display a plateau but we remark that the behaviour before the formation of this plateau is very different.

It is instructive to compare the temporal evolution of the blocks in entanglement hamiltonian matrix $H_{A}(t)$ with the one of the blocks in the reduced covariance matrix $\gamma_{A}(t)$ (see (98)). The main difference is the fact that the blocks of $H_{A}(t)$ are not homogeneous along their diagonals (see Fig. 2 and Fig. 3), while the blocks of $\gamma_{A}(t)$ are constant along any given diagonal (see the bottom right panel of Fig. 6, that corresponds to a given time), as it can be easily inferred from (95). In Fig. 6 we consider the temporal evolution of the antidiagonals of the blocks in $\gamma_{A}(t)$. Comparing this figure with Fig. 5, one concludes that also the width of the bands occurring in the three blocks of $\gamma_{A}(t)$ increase linearly in time with velocity equal to one. Nonetheless, let us stress that the shape of the antidiagonals occurring in $H_{A}(t)$ and in $\gamma_{A}(t)$ is very different.

As first quantitative check, we have compared the linear growth of the entanglement entropy obtained from our numerical data for the harmonic chain with the CFT formula $S_{A}\simeq 2\pi c\,t/(3\tau_{0})$ [44, 63] (see also §4) with $c=1$, finding $\tau_{0}\simeq 3.79$ (which agrees e.g. with the numerical value obtained in [42])

5.2 Entanglement spectrum

A numerical analysis for some quantities related to the entanglement spectrum has been reported in Fig. 7. The top left panel shows the symplectic spectrum of $H_{A}$ for various times and two different values of $\ell$. For the low-lying part of this spectrum, i.e. small $k$’s, we have $\varepsilon_{k}\simeq\varepsilon_{k+1}$ and this degeneracy disappears for high values of $k$. At small times, $\varepsilon_{k}$ is linear in terms of $k$ while for long time the curve of $\varepsilon_{k}$ bends towards zero, although it never vanishes. It would be interesting to consider $\varepsilon_{k}/\ell$ in the limit of large values of $\ell$. In the top right panel of Fig. 7 we consider the temporal evolution of the largest eigenvalue of the entanglement spectrum (see the text below (21)). Since $S_{A}^{(n)}\to-\log\lambda_{\textrm{\tiny max}}$ as $n\to\infty$, this curve can be compared with the CFT results for the Rényi entropies [44, 63] and a good agreement is found. Fitting the linear growth of $-\log\lambda_{\textrm{\tiny max}}$ we get a slope of $\tfrac{1}{2.1}$ times the slope of the entanglement entropy, while the one predicted by (75) is $\tfrac{1}{2}$.

In the remaining panels of Fig. 7 we explore the temporal evolutions of the gaps in the entanglement spectrum. The data about these evolutions display two distinct temporal regimes separated by $t/\ell\simeq 1/2$. For $t/\ell<1/2$, the CFT result (76), obtained from the analysis of [15] for the semi-infinite line, predicts linear growths in time for the inverse of the gaps with slopes proportional to the conformal spectrum (including also the dimensions of the descendants) allowed by the proper conformal boundary conditions. The numerical data in the middle left panel of Fig. 7 display these linear growths for $t/\ell<1/2$. By fitting the slopes of these linear growths through (76) with $\tau_{0}\simeq 3.79$ obtained above from the linear growth of the entanglement entropy, we find $\Delta_{1}=1.022$, $\Delta_{2}=2.045$ and $\Delta_{3}=3.058$. In order to reduce the influence of the initial state encoded in $\tau_{0}$, in the middle right panel and in the bottom panel of Fig. 7 we consider the temporal evolutions of $g_{r}S_{A}$ and the ratios $g_{r}/g_{1}$ respectively, which should be independent of $\tau_{0}$ according to (75) and (76). As for $g_{r}S_{A}$, it is evident that curves corresponding to different values of $\ell$ do not collapse; hence more values of $\ell$ (possibly also larger than the ones considered here) are needed in order to make comparisons with CFT results.

The temporal evolutions of the ratios $g_{r}/g_{1}$ of the entanglement gaps in the bottom panel of Fig. 7 display interesting features. For $t/\ell<1/2$ the curves having different $\ell$’s collapse forming plateaux whose heights are given by strictly positive integers. This result agrees with the fact that the underlying CFT contains the primary $\partial_{z}\phi$ and its descendants.

In Fig. 8 we show numerical data for the temporal evolution of the gaps in the entanglement spectrum when the initial state has $\omega_{0}=5$, in order to explore the robustness of the observations made above under changes of the initial state. For this $\omega_{0}$ we consider small values of $\ell$ such that the product $\omega_{0}\ell$ gets the same values corresponding to the data shown in Fig. 7. The qualitative behaviour of $1/g_{r}$ is the same one observed in Fig. 7 (the slopes of the linear growths are different, as expected, being (76) dependent on the initial state through $\tau_{0}$). In the temporal evolution of $g_{r}S_{A}$ we observe oscillations that are due to the small values of $\ell$ (indeed, they do not occur in the middle right panel of Fig. 7). Interestingly, the temporal evolutions of the ratios $g_{r}/g_{1}$ coincide with the ones reported in Fig. 7, meaning that this quantity displays some independence on the initial state. It would be instructive to consider other values of $\omega_{0}$ and higher gaps in order to understand better how much the temporal evolutions of $g_{r}/g_{1}$ are robust under modifications of the initial state. In the bottom right panel of Fig. 8 we have considered also long times and from these data we can identify two regimes: $t/\ell\leqslant 1/2$ and the long time regime. In both these temporal regimes we observe plateaux having the same heights. Thus, the ratios in the CFT spectrum can be read also from the long time regime.

5.3 A contour function from the quasi-particle picture

In [44] a quasi-particle picture has been introduced to explain the temporal evolution of the entanglement entropy after a global quantum quench. The underlying idea is based on the fact that the initial state has very high energy with respect to the ground state of the hamiltonian governing the temporal evolution; hence it can be seen as a source of quasi-particle excitations. In particular, in one spatial dimension, it is assumed that at $t=0$ each point of the space emits two quasi-particles with opposite momenta $p$ and $-p$ according to certain probability distribution that depends on both the initial state and the evolution hamiltonian. Only the particles emitted at the same point are entangled and all the points of the space emit the quasi-particles in the same way. For $t>0$, the positions of the quasi-particles emitted at the same point $x$ are $x+v_{p}t$ and $x-v_{p}t$, being $v_{p}>0$ and $v_{-p}=-v_{p}$ for $p>0$.

Considering a spatial bipartition $A\cup B$ and two points $x_{1}\in A$ and $x_{2}\in B$; at time $t$ they are entangled only if they are reached simultaneously by two quasi-particle emitted from the same point $x$ at $t=0$. The bipartite entanglement between $A$ and $B$ is obtained by summing the contributions of all the points $x$ fulfilling this condition. In particular, for the entanglement entropy we have [44]

[TABLE]

where $\tilde{s}(p)$ is obtained by multiplying the momentum distribution function and the contribution of the pair of quasi-particles with momenta $p$ and $-p$ to the entanglement entropy. The integration domain of the allowed momenta is model dependent. Performing the spatial integrals in (99), one finds

[TABLE]

When a maximum velocity $v_{\textrm{\tiny max}}$ exists, $S_{A}$ grows linearly in time for $v_{\mathrm{max}}t\leqslant\ell/2$, while $S_{A}\propto\ell$ for $v_{\mathrm{max}}t\gg\ell/2$. In a CFT we have $|v_{p}|=v_{\mathrm{max}}=1$ for all $p$ [44, 63]. Notice that (100) does not take into account the initial value of the entanglement entropy.

The expression (100) of the entanglement entropy obtained from this quasi-particle picture provides some functions $s_{A}(x,t)$ such that $S_{A}(t)=\int_{A}s_{A}(x,t)dx$. Let us introduce three positive functions $f_{0}(x)$, $f_{1}(x,p,t)$ and $f_{2}(x,p,t)$ fulfilling the following conditions

[TABLE]

Taking into account the initial value $S_{A}|_{t=0}$ in (100) and employing the constraints (101), it is straightforward to write the entanglement entropy as $S_{A}(t)=\int_{A}s_{A}(x,t)\,dx$ with

[TABLE]

Since our configuration is symmetric with respect to the center of the interval, it is natural to require that $f_{0}(x)=f_{0}(\ell-x)$ and $f_{i}(x,p,t)=f_{i}(\ell-x,p,t)$ for $i=1,2$. As mentioned above, (99) is obtained by assuming that the emission of the quasi-particles is spatially homogeneous, and this assumption provides a significant restriction to the form of the functions in (102).

Let us consider the infinitesimal contribution to the entanglement entropy $dS_{A}(t)=s_{A}(x,t)dx$ provided by an infinitesimal interval $(x-dx/2,x+dx/2)$ centered in a point $x\in A$. At $t>0$, this quantity is proportional to the number of quasi-particles (i) that are in this infinitesimal interval (ii) whose entangled quasi-particle is in $B$. The quasi-particles fulfilling these conditions have been emitted at $t=0$ from the position $x-|v_{p}|t$ or from $x+|v_{p}|t$. For $v_{p}\,t\leq\ell/2$, the quasi-particles emitted at $x-|v_{p}|t$ contribute when $2|v_{p}|\,t>x$, while the ones emitted at $x+|v_{p}|t$ matter when $2|v_{p}|\,t>\ell-x$. Instead, for $v_{p}\,t>\ell/2$ the quasi-particles coming from both $x-|v_{p}|t$ and $x+|v_{p}|t$ contribute. These considerations lead to the following expressions

[TABLE]

where $\varphi_{j}$ are non negative, $\varphi_{2}(\ell-x)=\varphi_{2}(x)$ and $\Theta$ is the Heaviside step function. The spatial homogeneity in the quasi-particle production leads to drastic simplifications given by $\varphi_{1}(x)=\varphi_{2}(x)=1$ identically.

The functions in (103) fulfil (101) and the symmetry conditions introduced in the text below (102). Thus, the contour function (102) becomes

[TABLE]

Integrating this expression between two generic points $0<x_{1}<x_{2}<\ell$, we obtain

[TABLE]

The function $\tilde{s}(p)$ can be computed by employing the fact that the density of thermodynamic entropy in the stationary state coincides with the one of the entanglement entropy in (100) [60].

For the global quench in the harmonic chain that we are considering, the stationary values of local observables can be described by a Generalised Gibbs Ensemble [50, 57] (we refer the interested reader to the review [65] for an extensive list of references). The velocity $v_{p}$ of the quasi-particles can be computed as $v_{p}=\partial_{p}\omega_{p}$, being $\omega_{p}$ the dispersion relation of the model. For the harmonic chain, in [60, 46] it has been found that

[TABLE]

where $\omega_{p}\equiv\sqrt{\omega^{2}+\tfrac{4\kappa}{m}[\sin(p/2)]^{2}}$ (see (90)) and $\omega_{0,p}$ is obtained by replacing $\omega$ with $\omega_{0}$ in $\omega_{p}$. The dispersion relation (90) provides the velocity of each momentum mode as follows

[TABLE]

Let us remark that the above results based on the quasi-particle picture hold for any value of the frequencies $\omega_{0}$ and $\omega$, while in our numerical analysis $\omega_{0}=1$ and $\omega=0$. It would be interesting to consider also other values for these parameters.

In Fig. 9 we show the temporal evolution of the contour for the entanglement entropy evaluated through the prescription of [42], outlined in §2.2. This qualitative behaviour, which has been observed also in [40] for fermionic chains after a global quench, has been obtained through the naive analysis performed in §4 by employing the CFT formulas of [15] (see the right panel of Fig. 1). The curves in Fig. 9 can be interpreted in terms of two fronts starting from the endpoints and travelling in the opposite directions towards the center of the interval. Each front has a plateau, that is not exactly horizontal in the numerical data. Since the two fronts have velocity equal to $2$ (see also §4), they cross each other around $t/\ell\simeq 1/4$ and superpose until they reach the opposite endpoint around $t/\ell\simeq 1/2$. Notice that the data shown in the two panels of Fig. 9 do not overlap, meaning that larger values of $\ell$ are needed to obtain a prediction for this curve in the limit of large $\ell$. This prediction can be done from the data obtained for the global quench in the chain of free fermions considered in §6 (see Fig. 21). We remark that the data in Fig. 9 display a divergence close to the endpoints of the interval that is independent of time.

The contour function for the entanglement entropy obtained through the quasi-particle picture (namely by employing (104), (106) and (107)) is shown in Fig. 10, where in the left panel $f_{0}(x)=0$ identically, while $f_{0}(x)=s^{(1)}_{A,0}(x,0)$ in (82) has been chosen in the right panel. The latter choice leads to reproduce also the linear divergencies of the contour function for the entanglement entropy near the endpoints of the interval. The quasi-particle picture formula (104) captures in a better way some features of the numerical data with respect to the corresponding CFT expression (81). For instance, the curve in Fig. 10 having $t/\ell=1/4$ is not flat in the middle of the interval and that the local maximum of the curve having $1/4<t/\ell<1$ increases with time, while in the right panel of Fig. 1 it is constant and equal to half of its asymptotic value for $t\to\infty$.

In Fig. 11 we consider the quantity introduced in (33) with $n=1$ for three different choices of $(i_{1},i_{2})$: the case with $i_{1}=0$ (top panel), the case with fixed $0<i_{1}<\ell/2$ and $i_{2}$ variable (middle panel) and the case with $(i_{1},i_{2})$ in the middle of the interval $A$ (bottom panel). The data points (obtained through (31), (32), (33) and (36)) are compared against the corresponding formula coming from the quasi-particle picture expression (105) (black dashed curves), where the terms containing $f_{0}$ cancel, and also against the naive CFT expression (85) (grey dashed-dotted curves). Notice that the CFT curves and the ones obtained from the quasi-particle picture change slope at the same values of $t/\ell$. At the beginning of the evolution, $\Delta\mathcal{S}_{A}$ grows linearly with a slope that depends on whether $(i_{1},i_{2})$ contains one or two endpoints of $A$ (notice that the red curves in the top and bottom panels correspond to the entanglement entropy), while it vanishes whenever the endpoints of $A$ do not belong to $(i_{1},i_{2})$. For $t/\ell>1/2$ all the pairs of quasi-particles have contributed and the curves reach a constant value proportional to the size of the interval $(i_{1},i_{2})$.

For $t/\ell<1/2$, all the curves in Fig. 11 exhibit a piecewise linear behaviour and, for a given configuration, the slope changes at most four times in correspondence to values of $t/\ell$ that depend on $(i_{1},i_{2})$ and that are captured both by the naive CFT expression (85) and by the quasi-particle picture expression (105). In the top panel of Fig. 11 the changes of slope occur at $t/\ell=x_{0}/(2\ell)$, $t/\ell=(\ell-x_{0})/(2\ell)$ and $t/\ell=1/2$; in the middle panel at $t/\ell=x_{1}/2$, $t/\ell=(a_{1}+x_{1})/(2\ell)$, $t/\ell=(\ell-a_{1}-x_{1})/(2\ell)$ and $t/\ell=(\ell-x_{1})/(2\ell)$ (in the data shown correspond to $x_{1}=\ell/5$); and in the bottom panel at $t/\ell=(\ell/2-a_{0})/(2\ell)$ and $t/\ell=(\ell/2+a_{0})/(2\ell)$. These values come from the quasi-particle picture. First we observe that the segments composing any piecewise linear curve are horizontal or have a positive slope that can take two values such that one is twice the other one. The quasi-particles can entangle $(i_{1},i_{2})$ with $B$ by crossing one or both the endpoints of $A$. Denoting by $\mathfrak{n}_{R}$ ($\mathfrak{n}_{L}$) the number of particles entangling $(i_{1},i_{2})$ with $B$ across the right (left) endpoint of $A$, we have that the segment is horizontal whenever both $\mathfrak{n}_{R}$ and $\mathfrak{n}_{L}$ are constant in time. Instead, when only either $\mathfrak{n}_{R}$ or $\mathfrak{n}_{L}$ is increasing in time, the corresponding segment has a positive slope, that becomes twice this value whenever both $\mathfrak{n}_{R}$ and $\mathfrak{n}_{L}$ are increasing in time. By adapting these considerations to the configurations for $(i_{1},i_{2})$ considered in the three panels of Fig. 11, one obtains the above values of $t/\ell$ corresponding to the changes in the slope of the piecewise linear curves.

6 Interval in a chain of free fermions

In this section we study the temporal evolution of the entanglement hamiltonian matrix after the global quench in a chain of free fermions introduced in [49, 1].

Let us consider the following inhomogeneous hamiltonian written in terms of the fermionic creation and annihilation operators $c_{n}$ and ${c}_{n}^{\dagger}$, (that satisfy the standard anticommutation relations $\{c_{n}^{\dagger},c_{m}^{\dagger}\}=\{c_{n},c_{m}\}=0$ and $\{c_{n},c^{\dagger}_{m}\}=\delta_{m,n}$) [49]

[TABLE]

where $t_{2n}=1$ and $t_{2n+1}=0$, namely only pairs of sites are coupled (dimerized chain). The system is half filled and prepared in the ground state $|\psi_{0}\rangle$ of $H_{0}$.

At $t=0$ the inhomogeneity is removed and the unitary time evolution of $|\psi_{0}\rangle$ is governed by the translation invariant hopping hamiltonian given by

[TABLE]

which is also known as the tight binding model at half filling.

In [49] the analytic expression of the correlation matrix after the global quench has been computed. Its generic element $C_{i,j}(t)\equiv\langle c_{i}^{\dagger}(t)\,c_{j}(t)\rangle$ reads

[TABLE]

where $J_{\nu}(x)$ is the Bessel function of the first kind and

[TABLE]

By employing that $J_{-\nu}(x)=(-1)^{\nu}J_{\nu}(x)$, we have that $\textrm{Re}[C_{i,j}(t)]$ is symmetric and that $\textrm{Im}[C_{i,j}(t)]$ is antisymmetric. Furthermore, we remark that $C_{i,j}(t)\to C^{(\infty)}_{i,j}$ as $t\to+\infty$; hence the correlation matrix is real in the asymptotic regime of long time.

6.1 Entanglement hamiltonian matrix

Considering an interval $A$ containing $\ell$ sites labelled by $1\leqslant i\leqslant\ell$, the generic element of the $\ell\times\ell$ correlation matrix $C_{A}(t)$ is obtained by just restricting to the rows and columns corresponding to the sites in $A$, namely $C_{A}(t)_{i,j}=C(t)_{i,j}$ in (110) with $1\leqslant i,j\leqslant\ell$. The entanglement hamiltonian of $A$ after this global quench is the operator (65) with the matrix $T$ given by (68), i.e. [25]

[TABLE]

This implies that the relation between the eigenvalues $\eta_{k}$ of $T$ and the eigenvalues $\zeta_{k}$ of $C_{A}$ reads $\eta_{k}=\log(1/\zeta_{k}-1)$, as already discussed in §3.3. For a generic value of $t$, the entanglement hamiltonian matrix $T$ in (112) is a complex $\ell\times\ell$ matrix. Since $T$ is hermitian at any $t$, its real part is symmetric and its imaginary part is antisymmetric. Thus, we can focus on the $p$-th diagonal of these two matrices with $p\geqslant 0$.

The eigenvalues $\zeta_{k}\in(0,1)$ and many of them lie exponentially close to 0 and 1. This forces us to work with very high precision in order to get finite results for (112). In particular, the largest interval we considered has $\ell=600$ sites and we employed between 900 and 2500 digits, depending on the value of $t$. The needed working precision decreases as time increases.

Also for this global quench, the configuration given by an interval in the infinite line is symmetric with respect to the center of the interval and this leads to a symmetry in the elements of the $p$-th diagonal. In particular, the odd diagonals of $T$ are symmetric, while the even ones are antisymmetric, namely $T_{i,i+p}=T_{\ell-i-p+1,\ell-i-p+1}$ for odd $p$ and $T_{i,i+p}=-T_{\ell-i-p+1,\ell-i-p+1}$ for even $p$. The main diagonal vanishes identically.

In Fig. 12 and Fig. 13 we show respectively the real part of $T$ and the imaginary part of $T$ of an interval with $\ell=400$ for nine values of $t$. For $t\to\infty$ the imaginary part of $T$ vanishes (a detailed analysis of this asymptotic regime is performed in §6.4). We also checked numerically that the imaginary part of $T$ vanishes for $t\to 0$, as expected. In the complex matrix $T$ the odd diagonals are real, while the even diagonals are purely imaginary, like in the correlation matrix (110).

As time evolves, the amplitude of the elements of $T$ decreases; indeed, in Fig. 12 and Fig. 13 a zoom is needed for large values of $t$ in order to appreciate the fact that $T$ is non vanishing. It is straightforward to observe that at the beginning the main contribution is localised on the diagonals close to the main one and that, during the temporal evolution, also the diagonals corresponding to longer range interactions become more important.

In Fig. 14 we show the temporal evolution of the first four diagonals in the entanglement hamiltonian matrix. Notice that $T_{i,i+1}$ is always positive, while in the other panels also negative values occur. We find it worth remarking that all the curves displayed in this figure vanish at the endpoints of the interval. The spatial inhomogeneity of these diagonals is a characteristic feature of the entanglement hamiltonian matrix, if compared with the correlation matrix restricted to the interval (see Fig. 19).

It is worth studying the limit $\ell\to\infty$ of the entanglement hamiltonian matrix, as done [31, 32, 29] for a static configurations (see [34] for the continuum limit). In Fig. 15 we consider increasing values of $\ell$ for some fixed values of $t/\ell$. While the data have reached the asymptotic curve when $t/\ell=0.1$ and $t/\ell=0.5$, it seems that larger values of $\ell$ are needed to observe this collapse for higher values of $t/\ell$.

In Fig. 12 and Fig. 13 it is straightforward to observe that the diagonals providing a long range interaction become more relevant as time evolves by forming a band around the main diagonal whose width $2d_{0}$ increases with time. A similar feature has been highlighted for the global quench in the harmonic chain discussed in §5.1. This observation naturally leads to consider the antidiagonals in $\textrm{Re}(T)$ and $\textrm{Im}(T)$, that are shown in Fig. 16 and Fig. 17, making evident the occurrence of a band and the increasing of its width during the temporal evolution. The width $2d_{0}$ of this band can be defined by taking the distance between the two minima for $\textrm{Re}(T)$ and the distance between the minimum and the maximum for $\textrm{Im}(T)$. In the bottom panels of Fig. 18 we show that $d_{0}(t)$ defined in this way increases linearly with velocity equal to one as time evolves.

A similar band, whose width increases linearly with velocity equal to one as well, occurs also in the temporal evolution of $\textrm{Re}(C_{A})$ and $\textrm{Im}(C_{A})$ whose antidiagonals display a behaviour like the one shown in Fig. 16 and Fig. 17 for the entanglement hamiltonian matrix (see Fig. 19 for $t/\ell=0.2$). Thus, the main difference between the temporal evolutions of $T$ and of $C_{A}$ is given by the inhomogeneity of the former one along the diagonals, that can be observed by comparing the top panels in Fig. 19 with the corresponding panels in Fig. 12 and Fig. 13.

6.2 Entanglement spectrum

Also for this global quench we find it worth considering the temporal evolution of some quantities related to the entanglement spectrum. In the top left panel of Fig. 20, the eigenvalues $\eta_{k}=\log(1/\zeta_{k}-1)$ of the hermitian matrix $T$ obtained from the eigenvalues $\zeta_{k}$ of the reduced correlation matrix $C_{A}(t)$ are shown for different values of time. The solid black curve corresponds to the asymptotic curve for $t\to\infty$ obtained in [49]. The eigenvalues $\eta_{k}$ occur in the expression (67) of the entanglement hamiltonian matrix $T$. Since $\eta_{k}=-\eta_{\ell-k}$ and only even values of $\ell$ are considered, we find it convenient to label these eigenvalues as $\eta_{k}$ with $k=\pm\frac{1}{2},\,\pm\frac{3}{2},\dots\pm\frac{\ell-1}{2}$. This gives $\eta_{k}=-\eta_{-k}$. Since $\widehat{K}_{A}=\sum_{k}\eta_{k}\,\hat{\mathfrak{f}}^{\dagger}_{k}\,\hat{\mathfrak{f}}_{k}$, by adapting to this model the steps that provide (57) for the harmonic chain, one finds that the largest eigenvalue in the entanglement spectrum is given by the configuration where all the modes with negative $\eta_{k}$ are occupied. This leads to write the largest eigenvalue as \lambda_{\text{\tiny max}}=\big{[}\prod_{k}(1+e^{-\eta_{k}})\big{]}^{-1}\prod_{k<0}\,e^{-\eta_{k}}, whose temporal evolution for $\ell=200$ and $\ell=400$ is shown in the top right panel of Fig. 20.

The gaps in the entanglement spectrum can be computed from the single particle entanglement energies $\eta_{k}$ by considering the particle-hole excitations with respect to the Fermi level. For the first gaps $0<g_{1}<g_{2}<\dots$, this analysis gives $g_{1}=2\eta_{\frac{1}{2}}$, $g_{2}=\eta_{\frac{1}{2}}+\eta_{\frac{3}{2}}$ and g_{3}=\textrm{min}\big{\{}\eta_{\frac{1}{2}}+\eta_{\frac{5}{2}}\,,2\eta_{\frac{3}{2}}\big{\}}. The expressions for the higher gaps become difficult to write as $r$ increases, but this can be done numerically in a systematic way. For these entanglement gaps we carry out an analysis similar to the one performed for the corresponding quantities in the harmonic chain (see Fig. 7), whose results are collected in the middle and bottom panels of Fig. 20. Also for this global quench the data about the temporal evolution of the entanglement gaps display two distinct regimes separated by $t/\ell\simeq 1/2$. In the middle left panel of Fig. 20 we show the temporal evolution of the inverse of the entanglement gaps, in order to highlight the linear growths predicted by the CFT formula (76) before $t/\ell\simeq 1/2$. The numerical value for $\tau_{0}$ is obtained by fitting the linear growth of the entanglement entropy through the CFT formula $S_{A}\simeq 2\pi c\,t/(3\tau_{0})$ of [44, 63] (see §4) with $c=1$ for the free massless Dirac fermion, finding $\tau_{0}\simeq 3.26$. This allows to fit the slope of the linear growths in the middle left panel of Fig. 20 through (76), finding $\Delta_{1}=1.224$, $\Delta_{2}=2.448$ and $\Delta_{3}=3.662$. We find it curious that the linear growths in the middle left panels of Fig. 7 and Fig. 20 basically coincide, despite the diversity of the underlying models (notice that the values of $\tau_{0}$ are different in the two quenches).

As done for the harmonic chain, in order to reduce the influence of the parameter $\tau_{0}$, we consider the temporal evolutions of $g_{r}S_{A}$ and of the ratios $g_{r}/g_{1}$. As for $g_{r}S_{A}$ (middle right panel of Fig. 20), the curves corresponding to different values of $\ell$ collapse much better than the ones obtained for the quench in the harmonic chain (middle right panel of Fig. 7) because the values of $\ell$ considered for this fermionic chain are large enough. Nonetheless, the zooms in the insets show that these curves are distinguishable.

Also for this quench the temporal evolutions of the ratios $g_{r}/g_{1}$ in the entanglement spectrum (bottom panel of Fig. 20) display the most interesting features. When $t/\ell<1/2$, curves having different $\ell$’s nicely collapse identifying neatly plateaux that correspond to strictly positive integers. The same positive integers are obtained also in the asymptotic regime of long time and large $\ell$. We checked numerically this result by plugging into the code employed to study $g_{r}/g_{1}$ the asymptotic values for $\eta_{k}$ found in [49], that will be reported later in §6.4 (see (117)). The plateaux before $t/\ell\simeq 1/2$ and for long time should correspond to the ratios between the conformal dimensions in the spectrum of the underlying CFT, and the data agree with the expected results for a massless Dirac fermion with free boundary conditions. We remark that in our analysis we considered only the particle-hole excitations. It would be interesting to improve this numerical analysis by including higher gaps or by considering other configurations.

In some spin chains at the critical point, the entanglement spectrum for an interval at equilibrium has been studied numerically in [37] finding the conformal spectrum of a boundary CFT with free boundary conditions. In the Ising model, the entanglement spectrum after a global quench has been considered in [56] and also in this model the gaps close as time evolves when the evolution is determined by a critical hamiltonian.

6.3 A contour function from the quasi-particle picture

In order to construct a contour function for the entanglement entropy by employing the quasi-particle picture, we need to know the dispersion relation of the model after the quench and the entropy density in momentum space of the quasi-particles (see §5.3), which are quantities corresponding to the entire system. The dispersion relation of the hamiltonian (109), which determines the evolution of the system after the quench, is given by $\omega(\theta)=\cos\theta$, where $-\pi\leqslant\theta\leqslant\pi$ [73]. In order to compute the entropy density, let us consider the asymptotic state of the system after the quench. For $t\to\infty$, the correlators are given by (111) with $i$ and $j$ labelling two generic sites of the infinite chain. Considering a periodic chain made by an even number $L$ of sites, a circulant matrix $\widetilde{C}^{(\infty)}_{i,j}$ with $i,j\in[-L/2,L/2]$ is obtained, whose eigenvalues are $\zeta_{k}=\left[1+\cos(2\pi k/L)\right]/2$, with $k\in[-L/2,L/2]$. Plugging these eigenvalues into (71), one finds the asymptotic entropy

[TABLE]

After introducing $\theta=2\pi k/L$, we take the limit $L\to\infty$ of $S^{(\infty)}/L$ with $S^{(\infty)}$ given by (113), and the result provides the entropy density $s^{(\infty)}(\theta)$ as follows

[TABLE]

These observations lead to construct a contour function from the quasi-particle picture by employing $s^{(\infty)}(\theta)$ and $v(\theta)\equiv\omega^{\prime}(\theta)=\sin\theta$ into (104), being $\theta\in[-\pi,\pi]$.

In Fig. 21 we show the temporal evolution of the contour for the entanglement entropy evaluated through the prescription of [40] discussed in §3.3, namely by employing (29), (30), (70) and (72). Notice that the data corresponding to different values of $\ell$ collapse on the same curves. The black dashed curves are obtained from (104) with $f_{0}(x)=0$ and (114) as explained above: they display perfect agreement with the data, except in the neighbourhoods of the endpoints. Similarly to the case of the harmonic chain described in §5.3, the divergencies of the contour function close to the endpoints can be reproduced by choosing a suitable expression for $f_{0}(x)$ (e.g. (82) with $n=1$ and the proper $\tau_{0}$, as done for the harmonic chain in the right panel of Fig. 10). The qualitative behaviour of the curves in Fig. 21 is the same observed in [40] for a different global quench and it is roughly reprodiced also by the naive CFT formula (81) discussed in §4 (see the right panel in Fig. 1). The blue dashed line in Fig. 21 corresponds to the asymptotic regime $t\to\infty$ discussed below in §6.4. Also for this global quench the numerical data display linear divergencies close to the endpoints of the interval that are independent of time.

In Fig. 22 we show the temporal evolution of (33) with $n=1$ for the same choices of $(i_{1},i_{2})$ considered in Fig. 11. The qualitative behaviour of the curves is similar to the corresponding curves obtained for the global quench in the harmonic chain, hence the qualitative arguments reported §5.3 about the changes of the slopes holds also in this case. On the other hand, for $t/\ell>1/2$ the asymptotic values for long times are different from the ones observed in the harmonic chain, as expected from the fact that the asymptotic value of the entanglement entropy depends both on the model and on the initial state. We find worth remarking that for the chain of free fermions we have access to larger values of $\ell$ with respect to the harmonic chain, hence the curves corresponding to different values of $\ell$ display a better collapse than the ones shown in Fig. 11 for the harmonic chain, capturing the asymptotic curves for $\ell\to\infty$. The black dashed lines, obtained from the quasi-particle picture analysis by using (32), (33) and (104), display a very good agreement with the numerical data. The data are also compared against the corresponding naive CFT expression (85) (grey dotted-dashed curves), where $\tau_{0}=3.26$ has been used.

6.4 Long time regime

In the asymptotic regime $t\to\infty$, the imaginary part of the correlation matrix (110) vanishes, as already remarked above. In particular, $C_{i,j}(t)\to C^{(\infty)}_{i,j}$ being $C^{(\infty)}_{i,j}$ the real symmetric and tridiagonal matrix defined in (111), which has the same value along a given diagonal. Because of the simple structure of $C^{(\infty)}_{i,j}$, analytic expressions for its eigenvalues and for the corresponding eigenvectors can be written. The eigenvalues of $C^{(\infty)}$ read [49]

[TABLE]

In order to study the matrix occurring in the entanglement hamiltonian, we also need the eigenvectors corresponding to the eigenvalues (115), that are given by

[TABLE]

where the normalisation condition $\sum_{r}u_{k}(r)^{2}=1$ has been imposed. Then, the orthogonal matrix $\widetilde{U}$ having the eigenvector $u_{k}(r)$ in (116) as $k$-th column can be constructed and, by employing this matrix in (67), the generic element of the $\ell\times\ell$ entanglement hamiltonian matrix reads

[TABLE]

where $1\leqslant i,j\leqslant\ell$ and the single particle entanglement energies have been obtained from (68) and (115). Notice that $\sin(r\theta_{k})=0$ when $r=0$ or $r=L+1$. The matrix $T$ in (117) is symmetric, as expected; hence we can focus on its $p$-th diagonal with $0\leqslant p\leqslant\ell-1$, whose generic element is

[TABLE]

Notice that $\theta_{\ell-k+1}=\pi-\theta_{k}$, and this implies $\eta_{\,\ell-k+1}=-\eta_{\,k}$. By using this property and splitting the sum in (118) into a sum from $1$ to $\ell/2$ and a sum from $\ell/2+1$ to $\ell$, one can show that the even diagonals vanish, namely $T_{i,i+2p}=0$ for non negative integers $p$. Furthermore, simple trigonometric identities allow to observe that the $p$-th diagonal is symmetric with respect to its middle point, namely $T_{\ell-i-p+1,\ell-i+1}=T_{i,i+p}$.

In the top panels of Fig. 23, the expression (118) has been checked against numerical data corresponding to large values of $t$ for $p=1$ (left) and $p=3$ (right). The agreement is satisfactory, but the insets show that larger values of $t$ should be taken in order to improve it.

The contour for the entanglement entropies can be written by employing (31), (72) and (116), together with the observations made above. This gives

[TABLE]

where $s_{n}(\zeta_{k})$ is obtained by evaluating (70) or (71) for the $\zeta_{k}$’s given in (115). The expression (119) corresponds to the black dashed curve in Fig. 21 and to the solid lines in the bottom panel of Fig. 23. Agreement is observed with the data points obtained for large values of time.

It is useful to consider the limit $\ell\to\infty$ of the above expressions holding in the asymptotic regime of long time. As the entanglement hamiltonian matrix element (118), by neglecting the highly oscillating terms coming from $\cos(2i\theta_{k})$ and $\sin(2i\theta_{k})$, for $i$ far enough from $1$ and $\ell-p$, one finds (see also section 7.1 of [1])

[TABLE]

where $\eta(\theta)\equiv 2\log[\tan(\theta/2)]$ (see (117)).

The limit (120) is checked in the top panels of Fig. 24 for $p=1$ (left) and $p=3$ (right).

As for the entanglement entropy, from (71) and the eigenvalues in (115), one obtains

[TABLE]

that is the flat value reached for $t\to\infty$ both in the inset and in the main plot of Fig. 21. Also the horizontal dashed line in the bottom panel of Fig. 24 corresponds to (121). The consistency of this result can be checked by employing (119), that for $\ell\to\infty$ gives

[TABLE]

where we used that $1/2$ is the mean value constant coming from $[\sin(i\theta_{k})]^{2}$. For the Rényi entropies, from (70) one finds

[TABLE]

that can be evaluated for explicit values of $n$, finding e.g. $\log(24-16\sqrt{2})$ for $n=2$ and $\log(4/3)$ for $n=3$.

The previous analysis has provided the values of the entanglement hamiltonian and of the contour for the entanglement entropies in the central part of the interval. In the following we improve this analysis by capturing also the behaviour of these quantities close to one of the endpoints. Here we only mention the results, reporting the details of their derivation in the Appendix D222These results and their derivation reported in the Appendix D have been obtained by Ingo Peschel. We are grateful to him for allowing us to include his analysis in this manuscript. .

When $\ell\to\infty$, the expression in (117) becomes

[TABLE]

where $1\leqslant i\leqslant\ell$ and $1\leqslant j\leqslant\ell$. This implies that the generic element of the $p$-th diagonal is $T_{i,i+p}=0$ when $p$ is even and $T_{i,i+p}=-4i/[p(2i+p)]$ when $p$ is odd. Taking $i\to\infty$ in these expressions, one finds $T_{i,i+p}\approx-2/p+1/i+O(1/i^{2})$ for odd $p$; hence (120) is recovered. A numerical check of (124) is reported in the top panels of Fig. 25.

The contour for the entanglement entropy can be analysed in a similar way. Taking $\ell\to\infty$ in (119) with $n=1$, one obtains

[TABLE]

being $i=1,2,3,\dots$ up to infinity. The limit $i\to\infty$ of this expression gives $s_{A}(i)\to\log 4-1+O(1/i^{3})$, which allows to recover (122) and also to find the power law decay $O(1/i^{3})$ for the contour function $s_{A}(i)$ close to the endpoint. The expression (125) has been checked numerically in the bottom panel of Fig. 25.

In this analysis the interval becomes a semi-infinite line; hence we can study the behaviour of the entanglement hamiltonian matrix and of the contour function only around one of the endpoints of the interval.

7 Conclusions

We studied the temporal evolution of the entanglement hamiltonian of an interval after a global quantum quench in two simple free lattice models: the harmonic chain and a chain of free fermions. In the harmonic chain (in the thermodynamic limit) we considered a quench of the frequency parameter such that the evolution hamiltonian is massless. In the chain of free fermions at half filling, we explored the global quench introduced in [49], where the initial state is the ground state of a dimerised chain and the evolution hamiltonian is determined by the homogenous hopping hamiltonian. In these free models, the time dependent entanglement hamiltonian is a quadratic operator; hence it is fully determined by a matrix that can be written explicitly as a function of the matrix whose elements are the two point correlators [23, 24, 2, 28, 25]. The entanglement hamiltonian matrix can be decomposed in terms of a diagonal matrix containing the single particle entanglement spectrum and another matrix that is symplectic for the harmonic chain or orthogonal for the chain of free fermions. The single particle entanglement spectrum and the matrix occurring in this decomposition can be employed to construct also a contour for the entanglement entropies [40, 42], that encodes information about the spatial structure of the bipartite entanglement. For static configurations, some entanglement hamiltonian matrices in the harmonic chain and in a chain of free fermions have been explored in [23, 24, 2, 28, 25, 29, 30, 33, 31, 32, 34].

In this manuscript we explored the temporal evolution of the entanglement hamiltonian matrices and of the contours for the entanglement entropy of an interval made by $\ell$ sites after the global quenches mentioned above. During the temporal evolution, the entanglement hamiltonian matrix is not block diagonal in the harmonic chain and it is complex in the chain of free fermions. All the diagonals of the $\ell\times\ell$ blocks composing the entanglement hamiltonian matrix in the harmonic chain and of the $\ell\times\ell$ real and imaginary parts of the entanglement hamiltonian matrix in the chain of free fermions (with the obvious exception of the ones that vanish at half filling) display a non trivial temporal evolution. Bands of diagonals around the main diagonal in these $\ell\times\ell$ matrices can be identified whose width growths linearly with velocity equal to one. This is observed also in the covariance matrix for the harmonic chain and in the correlation matrix for the chain of free fermions, where all the elements of each diagonal are equal.

The analytic results obtained in CFT for the entanglement hamiltonian of a semi-infinite line [15] have been exploited to write expressions for the finite interval in the continuum limit that reproduce qualitatively the behaviour of the numerical data of some quantities in some regimes. In particular, the linear growth of the gaps in the entanglement spectrum before $t/\ell\simeq 1/2$ and the qualitative temporal evolution of the contour for the entanglement entropy have been obtained. In our numerical analysis we have also observed that the temporal evolution of the ratios of the gaps in the entanglement spectrum can be employed to read the ratios in the conformal spectrum of the underlying CFT when the evolution hamiltonian is critical. This observation holds both for $t/\ell<1/2$ (as predicted in [15]) and in the asymptotic regime of long time. Furthermore, it seems true also for different initial states. It would be interesting to further explore this idea by considering global quenches governed by different critical hamiltonians.

The quasi-particle picture of [44] has been employed to obtain semi-empirical analytic expressions of the contour for the entanglement entropy, finding reasonable agreements with the numerical data, in particular for the chain of free fermions, where larger subsystems have been studied. For the harmonic chain, the formula for the contour is valid also for the global quenches where the mass in the evolution hamiltonian is non vanishing; hence it would be interesting to test it also for these protocols. For the global quench in the chain of free fermions, the entanglement hamiltonian matrix and the contour for the entanglement entropy have been studied in the asymptotic regime of $t\to\infty$.

The analysis of entanglement hamiltonians and the contours for the entanglement entropies can be extended in many interesting directions. For instance, it would be useful to find some analytic expressions describing the asymptotic regime of long time for the global quench considered here in the harmonic chain, like the ones discussed in §6.4 in the chain of free fermions. Other kinds of global quenches can be considered or also quenches of different nature, like the local quenches [48]. An interesting development is given by more complicated bipartitions, like the one where $A$ is made by disjoint intervals [13, 15, 33] whose entanglement entropies have been studied in CFT, lattice models [78] and holography [79]. Connections could be established also with some recent results in harmonic lattices [80].

It is worth studying the temporal evolutions of the entanglement hamiltonians and of the contours for the entanglement entropies after global and local quenches also in interacting models and in higher dimensions. Indeed, in this manuscript we have highlighted the role of the entanglement gaps to obtain the conformal dimensions of the underlying CFT when the evolution hamiltonian is critical; hence it would be interesting to test this idea also in interacting models. In the context of the gauge/gravity correspondence, the bit threads approach to the holographic entanglement entropy [81] suggests a possible gravitational dual of the contour for the entanglement entropy [82].

The entanglement hamiltonians and their dynamics are very important also because they are the building blocks to explore other interesting quantities that capture different aspects of entanglement. For instance, the bipartite entanglement of mixed states can be quantified through the partial transpose and the logarithmic negativity [83], which have been studied also in quantum field theories [84] and in some lattice models [85]. Other interesting directions to explore concern excited states [86] and quantities involving them, like e.g. the relative entropy [87].

Acknowledgements

It is our pleasure to thank Andrea Coser, Viktor Eisler, Ingo Peschel, Giulia Piccitto and Luca Tagliacozzo for important insights. We are also grateful to Horacio Casini, Fabian Essler, Paul Fendley, Juan Maldacena, Sara Murciano, Giuseppe Mussardo, Mohammad Ali Rajabpour, Germán Sierra, Hubert Saleur and Jacopo Surace for useful discussions. We thank an anonymous referee for observations that have led us to the contour function from the quasi-particle picture discussed in §6.3. This work started at the Galileo Galilei Institute (GGI) in Florence during the programme Entanglement in Quantum Systems in June and July 2018: we acknowledge GGI for financial support and the stimulating environment. RA thanks Associazione di Fondazioni e di Casse di Risparmio (ACRI) for financial support.

Appendices

Appendix A On the Williamson’s decomposition of $H_{A}$

In the harmonic lattices described in §2, the entanglement hamiltonian matrix (13) and its equivalent form (15) can be obtained from its Williamson’s decomposition (11). In this Appendix we discuss the details of this derivation.

By employing the Williamson’s decomposition (5) of the reduced covariance matrix $\gamma_{A}$, one finds

[TABLE]

Then, since $(\textrm{i}\mathcal{D}_{\textrm{\tiny d}}J)^{2}=(\textrm{i}J\mathcal{D}_{\textrm{\tiny d}})^{2}=\mathcal{D}_{\textrm{\tiny d}}^{2}$, we have that

[TABLE]

Considering a generic even function $f_{\textrm{\tiny e}}(x)$, the first expression in (127) leads to $f_{\textrm{\tiny e}}(\mathcal{D}_{\textrm{\tiny d}})=f_{\textrm{\tiny e}}(\textrm{i}\mathcal{D}_{\textrm{\tiny d}}J)=f_{\textrm{\tiny e}}(\textrm{i}J\mathcal{D}_{\textrm{\tiny d}})$. By employing (126) into this result, one obtains

[TABLE]

Similarly, for a generic odd function $f_{\textrm{\tiny o}}(x)$, the second expression in (127) implies that $f_{\textrm{\tiny o}}(\mathcal{D}_{\textrm{\tiny d}})=f_{\textrm{\tiny o}}(\textrm{i}\mathcal{D}_{\textrm{\tiny d}}J)\,\textrm{i}J=\textrm{i}J\,f_{\textrm{\tiny o}}(\textrm{i}J\mathcal{D}_{\textrm{\tiny d}})$. This result combined with (126) leads to

[TABLE]

Another useful relation can be derived by noticing that $(\textrm{i}J\gamma_{A})^{n}=J\,(\textrm{i}\gamma_{A}J)^{n}J^{-1}$ for any integer $n\geqslant 0$. This observation implies that, for a generic function $f(x)$, we have

[TABLE]

The relation (12) gives $\mathcal{E}_{\textrm{\tiny d}}=2\,\textrm{arccoth}(2\mathcal{D}_{\textrm{\tiny d}})$, that is an odd function of $\mathcal{D}_{\textrm{\tiny d}}$; hence, by specifying (129) to this case, we find that

[TABLE]

By isolating the matrix within the square brackets, one obtains the first expression in (13) with $H_{A}$ written through its Williamson’s decomposition (11) and $W_{H}$ given by (14). The second expression in (13) can be easily found from the first one by using (130).

Let us observe that, from the relations (11) and (10), we have $J^{\textrm{t}}H_{A}J=W^{\textrm{t}}\mathcal{E}_{\textrm{\tiny d}}W$, which tells that the Williamson’s decompositions of $\gamma_{A}$ and $J^{\textrm{t}}H_{A}\,J$ involve the same symplectic matrix $W$.

An equivalent way to express the entanglement hamiltonian matrix $H_{A}$ in terms of $\gamma_{A}$ is given by (15). Indeed, from the Williamson’s decomposition of $H_{A}$ in (11) with $W_{H}$ given by (14) and (10), by using that $\mathcal{E}_{\textrm{\tiny d}}=h(\mathcal{D}_{\textrm{\tiny d}})\mathcal{D}_{\textrm{\tiny d}}$, with $h(y)$ defined in the text below (15), one finds

[TABLE]

In the second step the identity matrix $\boldsymbol{1}=\widetilde{W}\,W^{\textrm{t}}$ has been properly inserted in order to recognise the Williamson’s decomposition (5) for $\gamma_{A}$ and the relation (128) for the even function $h(y)$, that allow to write the last expression of (132) in the form reported in the last step of (15). Similarly, by inserting the identity matrix $\boldsymbol{1}=(W^{\textrm{t}}J)^{-1}(W^{\textrm{t}}J)$ into the first expression in (132) in the proper way, we find that it can written as follows

[TABLE]

where, by using (5) and (128), one recognises the first expression in (15).

Another equivalent expression for the entanglement hamiltonian matrix $H_{A}$ can be found by observing that the relation (12) between $\mathcal{D}$ and $\mathcal{E}$ gives

[TABLE]

By employing these relations, (5), (10) and (11), we arrive to

[TABLE]

Plugging these expressions into (134), one easily finds

[TABLE]

and

[TABLE]

which provides another form to express the relation between the reduced covariance matrix $\gamma_{A}$ and the entanglement hamiltonian matrix $H_{A}$.

Appendix B On the numerical determination of the symplectic matrix $W$

The symplectic matrix $W$ occurring in the Williamson’s decomposition (5) of the reduced covariance matrix $\gamma_{A}$ plays a crucial role throughout this manuscript. In this Appendix we discuss the method employed to construct $W$ numerically.

The reduced covariance matrix $\gamma_{A}$ is symmetric and positive definite; hence, according to the Cholesky decomposition, there exists a unique way to write $\gamma_{A}$ as follows

[TABLE]

where $L_{\gamma}$ is a real lower triangular matrix. Plugging (138) into (9), one obtains

[TABLE]

The first expression in this sequence of equalities is a real symmetric matrix that can be diagonalised by an orthogonal matrix $W_{\gamma}$, namely $(JL_{\gamma})^{\textrm{t}}\,\gamma_{A}(JL_{\gamma})=W_{\gamma}^{-1}\Lambda\,W_{\gamma}$, being $\Lambda$ a diagonal matrix. Comparing this formula with (139), we find that $\Lambda=\mathcal{D}_{\textrm{\tiny d}}^{2}$ and $W_{\gamma}=\mathcal{N}\,WL_{\gamma}^{-\textrm{t}}$, where $\mathcal{N}$ is an invertible diagonal matrix. By imposing that $W_{\gamma}$ is orthogonal, one obtains $\mathcal{N}=\mathcal{D}_{\textrm{\tiny d}}^{1/2}$. Thus, the matrix $W$ can be written as

[TABLE]

By employing the fact that $W_{\gamma}$ is orthogonal and the Cholesky decomposition (138), it is straightforward to verify that (140) satisfies (5), as expected. The fact that (140) is symplectic can be checked by following the procedure adopted to verify this property in the proof of the Williamson’s theorem given in [70]. As consistency check, one can repeat this analysis starting from the last expression in (139) and this leads (140) again.

In our numerical analysis we have constructed $W$ by employing (140) and Wolfram Mathematica. In particular, given the reduced covariance matrix $\gamma_{A}$, the software provides the unique matrix $L_{\gamma}$ and the orthogonal matrix $W_{\gamma}$ that diagonalises the symmetric matrix $(JL_{\gamma})^{\textrm{t}}\,\gamma_{A}\,(JL_{\gamma})$. Given these two matrices and the symplectic spectrum obtained from (9), the symplectic matrix $W$ is constructed as in (140).

Appendix C Block diagonal entanglement hamiltonians in harmonic lattices

In the harmonic lattices and for static configurations, both the reduced covariance matrix $\gamma_{A}$ and the entanglement hamiltonian matrix $H_{A}$ are block diagonal [2, 29, 30, 33]. In this Appendix we specify the results discussed in §2.1 to the simpler case.

When $\gamma_{A}=Q\oplus P$ is block diagonal, namely $R=\boldsymbol{0}$ in (24), the symplectic matrix $W$ in the Williamson’s decomposition (5) of $\gamma_{A}$ is block diagonal as well, namely $W=U\oplus V$; hence the Williamson’s decomposition of $\gamma_{A}$ can be written as

[TABLE]

The condition that $W$ and $W^{\textrm{t}}$ are symplectic leads respectively to

[TABLE]

telling that $U^{-1}=V^{\textrm{t}}$ and $V^{-1}=U^{\textrm{t}}$. In particular, $U$ and $V$ are not orthogonal matrices. The relations (141) do not provide the diagonalization of the real and symmetric matrices $Q$ and $P$ because $U$ and $V$ are invertible but not orthogonal. Moreover, being $W=U\oplus V$, we have that (10) becomes $\widetilde{W}=V\oplus U$.

The matrices (9) to diagonalise in order to find the symplectic spectrum simplify to $(\textrm{i}J\gamma_{A})^{2}=(PQ)\oplus(QP)$ and $(\textrm{i}\gamma_{A}J)^{2}=(QP)\oplus(PQ)$ when $\gamma_{A}$ is block diagonal. From this observation and (142), the first expression in (9) gives

[TABLE]

When $\gamma_{A}$ is block diagonal, the entanglement hamiltonian matrix $H_{A}$ is block diagonal as well. In particular, from the above remarks and $J^{\textrm{t}}\gamma_{A}J=P\oplus Q$, one observes that the expressions in (15) reduce to [2]

[TABLE]

The equivalence of these two expressions can be verified by transposing one of them and employing that $M$ and $N$ are symmetric combined with $(\sqrt{QP}\,)^{\textrm{t}}=\sqrt{PQ}$, that can be easily obtained from the fact that also $Q$ and $P$ are symmetric. Another useful expression for the block diagonal entanglement hamiltonian matrix $H_{A}$ can be obtained by specifying (28) to this simpler case where $W=U\oplus V$ (i.e. $Y=Z=\boldsymbol{0}$ in (25)). This gives

[TABLE]

Notice that, for a block diagonal $\gamma_{A}$, in (138) we have $L_{\gamma}=L_{Q}\oplus L_{P}$; therefore also the expressions reported in the Appendix B simplify.

Appendix D On the regime of long time and large $\ell$ for the fermionic chain

In this appendix we discuss the derivation of the expressions (124) and (125) about the limit $\ell\to\infty$ of the regime of long time for the quench of the fermionic chain (see §6.4).

For a given function $f$, the large $\ell$ limit of the sum $\frac{1}{\ell+1}\sum_{k=1}^{\ell}f(\theta_{k})$ gives $\int_{0}^{\pi}\frac{d\theta}{\pi}f(\theta)$. Thus, for (117) in this limit we get

[TABLE]

where $\eta(\theta)$ has been defined in the text below (120). Splitting (147) as a sum of two integrals and focussing on one of them, we have

[TABLE]

where an integration by parts has been performed, observing that the boundary terms vanish. The expression (124) is obtained by employing (148) into (147).

As for the contour for the entanglement entropy, in the limit $\ell\to\infty$ we find that (119) becomes

[TABLE]

where the expression for $s(\theta)$ introduced in (122) can be written as

[TABLE]

Plugging this expression into (149), an integral to compute is given by

[TABLE]

where an integration by parts has been performed (the corresponding boundary terms vanish). The remaining integral becomes a sum of two simpler integrals once the product $\cos(\theta)\,\cos(2\theta i)$ is written as a sum of two single cosine functions. These two integrals in (149) can be evaluated as follows

[TABLE]

where $\mu=2\,i\pm 1$ is an odd positive integer. By employing (151) and (152) into (149), the expression (125) is obtained.

As consistency check for (125), let us recover the value of the entanglement entropy found in (121) by summing (125) over the sites of the interval $A$. Since $A$ has $\ell$ sites, the first term of this sum gives $\ell\,(\,\log 4-1\,)$. Then, by using that $\sum_{i=1}^{\infty}\frac{1}{2\,i\,(4\,i^{2}-1)}\,=\,\frac{1}{2}(\log 4-1)$ and taking into account that the interval has two endpoints, we obtain $\sum_{i\in A}s_{A}(i)=(\ell+1)(\log 4-1)$, which is the result reported in (121).