A contour for the entanglement entropies in harmonic lattices

TL;DR

This paper develops a contour function for entanglement entropies in harmonic lattices, validating it through numerical analysis and comparison with conformal field theory predictions, especially in one-dimensional massless regimes.

Contribution

It introduces a new contour function for entanglement entropy in harmonic lattices and compares it with theoretical predictions in specific regimes.

Findings

Good agreement with conformal field theory in massless regimes

Effective contour function for single and disjoint intervals

Numerical validation across different boundary conditions

Abstract

We construct a contour function for the entanglement entropies in generic harmonic lattices. In one spatial dimension, numerical analysis are performed by considering harmonic chains with either periodic or Dirichlet boundary conditions. In the massless regime and for some configurations where the subsystem is a single interval, the numerical results for the contour function are compared to the inverse of the local weight function which multiplies the energy-momentum tensor in the corresponding entanglement hamiltonian, found through conformal field theory methods, and a good agreement is observed. A numerical analysis of the contour function for the entanglement entropy is performed also in a massless harmonic chain for a subsystem made by two disjoint intervals.

Equations

$$S_A=-\mathrm{Tr}({\rho }_A\log {\rho }_A),$$

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

$$S_A=-\underset{n\to 1}{\lim }{\partial }_n\mathrm{Tr}{\rho }_A^n=\underset{n\to 1}{\lim }{S}_A^{(n)},$$

$$I_{A_1,A_2}^{(n)}=S_{A_1}^{(n)}+S_{A_2}^{(n)}-S_{A_1\cup A_2}^{(n)},$$

$${\rho }_A=e^{-2\pi K_A},$$

$$S_A^{(n)}=\sum _{i\in A}{s}_A^{(n)}(i),$$

$$s_A^{(n)}(i)⩾0,\forall i\in A.$$

$$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{)}\,,$$

$$\gamma_{A}=\bigg{(}\begin{array}[]{cc}Q&R\\ R^{\textrm{t}}&P\\ \end{array}\bigg{)}\,,$$

$$\gamma_{A}=W^{\textrm{t}}\big{(}D\oplus D\big{)}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{]}\,.$$

$$\sum _{i=1}^{\ell }{p}_k(i)=1,1⩽k⩽\ell .$$

$$s_A^{(n)}(i)=\sum _{k=1}^{\ell }{p}_k(i)s_n({\sigma }_k),$$

$$\Phi_{A}\equiv O^{\textrm{t}}\big{(}D\oplus D\big{)}\,O\,,$$

$$S_{A}^{(n)}=\,{\rm Tr}\big{[}s_{n}(D)\big{]}=\frac{1}{2}\,{\rm Tr}\big{[}s_{n}(D\oplus D)\big{]}=\frac{1}{2}\,{\rm Tr}\big{[}s_{n}(\Phi_{A})\big{]}\,,$$

$$S_A^{(n)}=\sum _{\alpha =1}^{\ell }{s}_A^{(n)}(\alpha ),$$

$$s^{(n)}_{A}(\alpha)=\frac{1}{2}\,{\rm Tr}\big{[}X^{(\alpha)}s_{n}(\Phi_{A})\big{]}\,.$$

$$\displaystyle s^{(n)}_{A}(\alpha)\;=\;\frac{1}{2}\,{\rm Tr}\big{[}X^{(\alpha)}s_{n}(O^{\textrm{t}}(D\oplus D)O)\big{]}=\frac{1}{2}\,{\rm Tr}\big{[}OX^{(\alpha)}O^{\textrm{t}}\,s_{n}(D\oplus D)\big{]}$$

$$\displaystyle\hskip 19.91684pt\;=\;\frac{1}{2}\,{\rm Tr}\big{[}O\big{(}\delta^{(\alpha)}\oplus\delta^{(\alpha)}\big{)}O^{\textrm{t}}\big{(}s_{n}(D)\oplus s_{n}(D)\big{)}\big{]}\,.$$

$$O=\bigg{(}\begin{array}[]{cc}U_{O}&Y_{O}\\ Z_{O}&V_{O}\end{array}\bigg{)}\,.$$

$$s^{(n)}_{A}(\alpha)=\frac{1}{2}\,\Big{(}{\rm Tr}\big{[}\big{(}U_{O}\,\delta^{(\alpha)}U_{O}^{\textrm{t}}+Y_{O}\,\delta^{(\alpha)}Y_{O}^{\textrm{t}}\big{)}s_{n}(D)\big{]}+{\rm Tr}\big{[}\big{(}Z_{O}\,\delta^{(\alpha)}Z_{O}^{\textrm{t}}+V_{O}\,\delta^{(\alpha)}V_{O}^{\textrm{t}}\big{)}s_{n}(D)\big{]}\Big{)}\,.$$

$$s_A^{(n)}(\alpha )=\sum _{k=1}^{\ell }{p}_k(\alpha )s_n({\sigma }_k),$$

$$p_{k}(\alpha)\,=\,\frac{1}{2}\Big{(}\big{[}(U_{O})_{k\alpha}\big{]}^{2}+\big{[}(Y_{O})_{k\alpha}\big{]}^{2}+\big{[}(Z_{O})_{k\alpha}\big{]}^{2}+\big{[}(V_{O})_{k\alpha}\big{]}^{2}\Big{)}\,,$$

$$OO^{\textrm{t}}=\bigg{(}\begin{array}[]{cc}U_{O}U_{O}^{\textrm{t}}+Z_{O}Z_{O}^{\textrm{t}}&U_{O}Y_{O}^{\textrm{t}}+Z_{O}V_{O}^{\textrm{t}}\\ Y_{O}U_{O}^{\textrm{t}}+V_{O}Z_{O}^{\textrm{t}}&Y_{O}Y_{O}^{\textrm{t}}+V_{O}V_{O}^{\textrm{t}}\end{array}\bigg{)}=\bigg{(}\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{0}\\ \boldsymbol{0}&\boldsymbol{1}\end{array}\bigg{)}\,.$$

$$\sum_{\alpha=1}^{\ell}\Big{(}\big{[}(U_{O})_{k\alpha}\big{]}^{2}+\big{[}(Z_{O})_{k\alpha}\big{]}^{2}\Big{)}=1\,,\qquad\sum_{\alpha=1}^{\ell}\Big{(}\big{[}(Y_{O})_{k\alpha}\big{]}^{2}+\big{[}(V_{O})_{k\alpha}\big{]}^{2}\Big{)}=1\,,$$

$$s^{(n)}_{A}(i)=\frac{1}{2}\,{\rm Tr}\big{[}X^{(i)}s_{n}(\Phi_{A})\big{]}\,,$$

$$p_{k}(i)\,=\,\frac{1}{2}\Big{(}\big{[}(U_{O})_{ki}\big{]}^{2}+\big{[}(Y_{O})_{ki}\big{]}^{2}+\big{[}(Z_{O})_{ki}\big{]}^{2}+\big{[}(V_{O})_{ki}\big{]}^{2}\Big{)}\,,$$

$$W=K_{\text{L}}E{K}_{\text{R}},E=e^{\chi }\oplus e^{-\chi },\chi \equiv \text{diag}({\chi }_1,\dots ,{\chi }_{\ell }),$$

Full text

Contents

@starttoc

toc

1 Introduction

Entanglement in many-body quantum systems has attracted a lot of research during the last decade (see [1] for a collection of reviews).

Consider a quantum system whose Hilbert space is bipartite, namely $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. Denoting by $\rho$ the density matrix which characterises the state of the whole system, $A$’s reduced density matrix $\rho_{A}$ is obtained by taking the partial trace of $\rho$ over $\mathcal{H}_{B}$. The entanglement entropy is defined as the Von Neumann entropy of $\rho_{A}$, i.e.

[TABLE]

where the normalization condition ${\rm Tr}\rho_{A}=1$ has been imposed.

Other important quantities providing useful information about the entanglement of the bipartition of the system in the state $\rho$ are the Rényi entropies

[TABLE]

which are parameterised by the integer $n\geqslant 2$. The entanglement entropy (1) can be computed from the Rényi entropies (2) by performing the following analytic continuation

[TABLE]

which is also known as the replica limit for the entanglement entropy. We refer to $S_{A}^{(n)}$ with $n\geqslant 1$ as the entanglement entropies, assuming that $S_{A}^{(1)}\equiv S_{A}$ is the entanglement entropy (1) obtained through the replica limit (3).

The entanglement entropies measure the bipartite entanglement associated to the decomposition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ when the whole system is in a pure state. In this manuscript we consider bipartitions of the Hilbert space associated to spatial bipartitions of the entire space, where $A$ is a spatial region and $B$ its complement.

Entanglement entropies have been largely studied in many quantum systems both in lattice models [2, 3, 4, 5, 6, 7] and in quantum field theories (QFTs). As for the latter class of models, important goals have been achieved in conformal field theories (CFTs) [8, 9, 10], free quantum field theories [11] and holography [12, 13].

A relatively simple spatial bipartition is obtained when $A$ is a simply connected domain. Nonetheless, it is very interesting to study also configurations where $A$ is made by disjoint regions. In the case of two disjoint domains $A_{1}$ and $A_{2}$, we have $A=A_{1}\cup A_{2}$ and an important quantity to introduce is the following combination of entanglement entropies

[TABLE]

where the $n=1$ case corresponds to the mutual information $I_{A_{1},A_{2}}\equiv S_{A_{1}}+S_{A_{2}}-S_{A_{1}\cup A_{2}}=\lim_{n\to 1}I_{A_{1},A_{2}}^{(n)}$, which can be obtained by taking the replica limit (3) of the various terms. The mutual information measures the total amount of correlations between the two disjoint regions [14]. The quantity (4) has been studied both in one [15, 16, 17, 18] and in higher spatial dimensions [11, 19].

In principle, the reduced density matrix $\rho_{A}$ contains more information about the entanglement between $A$ and its complement $B$ than the entanglement entropies $S_{A}^{(n)}$. Since $\rho_{A}$ is a positive semi-definite and hermitian operator, it can be written as

[TABLE]

where $K_{A}$ is known as the entanglement hamiltonian (or modular hamiltonian).

The entanglement hamiltonians for free fermions and free bosons on the lattice have been studied in [3] (see [6, 11] for reviews). In QFTs the entanglement hamiltonians are generically non local operators [20]. Nonetheless, there are interesting special configurations where $K_{A}$ can be expressed in terms of integrals over $A$ of local operators multiplied by proper local weight functions. The most important case belonging to this class is a Lorentz invariant QFT in its ground state with $A$ given by the half space $x_{1}>0$. For this example Bisognano and Wichmann found that the entanglement hamiltonian is the generator of the Lorentz boosts in the $x_{1}$ direction [21].

Considering the special class of QFTs given by CFTs, the crucial result by Bisognano and Wichmann and the conformal symmetry allow to write the entanglement hamiltonians corresponding to other interesting configurations as the integral over the domain $A$ of the component $T_{00}$ of the energy-momentum tensor multiplied by a suitable local weight function [22, 23, 24]. For free fermions in one spatial dimension, also the entanglement hamiltonian of disjoint intervals has been found [18].

In this manuscript we study functions $s_{A}^{(n)}:A\to\mathbb{R}$ which assign a real number to every site in the spatial domain $A$ such that

[TABLE]

where the subsystem $A$ is made by either a connected region or disjoint regions. Since this property naturally leads to interpret $s_{A}^{(n)}(i)$ as a density for the entanglement entropies [25], it is natural to impose also the following positivity condition

[TABLE]

The aim is to construct these functions $s_{A}^{(n)}(i)$ from the methods employed to compute the entanglement entropies, in order to identify the contribution of the $i$-th site to the entanglement entropies. These functions have to fulfil (6), (7) and possibly other proper requirements and they are expected to provide information about the spatial structure within $A$ of the entanglement between $A$ and $B$. The same question can be formulated for QFTs in the continuum, where the discrete sum in (6) becomes an integral over the region $A$ and $s_{A}^{(n)}$ is a function of the position $x\in A$.

A function fulfilling the constraints (6) and (7) for free fermions on the lattice has been studied by Chen and Vidal [26]. As for the harmonic lattices, a proposal which satisfies (6) is based on the results obtained by Botero and Reznik [4] and it has been studied in [27]. Following [26], throughout this manuscript the function $s_{A}^{(n)}(i)$ will be called contour function for the entanglement entropies.

Many functions satisfying (6) and (7) could be constructed and a complete list of properties characterising the contour function for the entanglement entropies in a unique way is not available. In [26], the authors made a step in this direction by proposing three further reqirements beside (6) and (7) involving the contour functions for the entanglement entropies. They are given by: (a) a constraint implementing the consistency with any spatial symmetry of the subsystem; (b) the requirement that the contour integrated over any subregion $G\subseteq A$ must be invariant under local unitary transformations in that subregion; (c) a bound meaning that the contour integrated over any subregion $G\subseteq A$ must be smaller or equal than the entanglement of any factor space of $\mathcal{H}_{A}$ which contains the Hilbert space $\mathcal{H}_{G}$ of that subregion. The authors of [26] constructed a contour function for the entanglement entropy of free fermions on the lattice which fulfils these five requirements. These properties could be helpful in order to characterise the contour function within the set made by the densities of entanglement, which can be naturally defined as the functions satisfying the properties (6) and (7). For instance, a natural function to consider is the flat contour $s_{A}^{(n)}(i)=S_{A}^{(n)}/|A|$, which assigns the same value to all the sites in $A$, being $|A|$ the total number of sites contained in $A$. This function fulfils the constraints (6) and (7), but it does not provide any useful information about the spatial structure of the entanglement within $A$. As for the contour function for the entanglement entropies in the harmonic lattices based on [4], the property (6) holds, but the positivity condition (7) has not been proved [27].

In this manuscript we construct new contour functions for the entanglement entropies in generic harmonic lattices which fulfil (6) and (7). For a specific proposal, we prove that it satisfies also a weaker version of the three further requirements introduced in [26]. A detailed numerical analysis of this contour function is performed for harmonic chains in one spatial dimension.

The continuum limit of the harmonic lattice in the massless regime is the CFT given by the free massless boson and in one spatial dimension it has central charge $c=1$. In two-dimensional CFTs, some time-independent examples have been found where the entanglement hamiltonian $K_{A}$ of a finite interval $A$ can be written as an integral over $A$ of the $T_{00}$ component of the energy-momentum tensor multiplied by a suitable local weight function [22, 23, 24]. Focusing on these configurations, for some examples we find that, in the massless regime and in the scaling limit, the contour function highlighted above is in good agreement with the inverse of the local weight function multiplying $T_{00}$ in the corresponding entanglement hamiltonian, specified to the $c=1$ case.

The manuscript is organised as follows. In §2 the usual method to compute the entanglement entropies in harmonic lattices is reviewed, highlighting the role of the Williamson’s theorem. In §3 a set of contour functions $s_{A}^{(n)}(i)$ satisfying (6) and (7) is constructed. Then, we focus on a specific proposal belonging to this set, showing that it also fulfils a weaker version of the three requirements introduced in [26]; therefore it is our best candidate for the contour function for the entanglement entropies in harmonic lattices. In §4 we consider this contour function in the massless regime and in the scaling limit for various configurations, comparing the lattice data with the inverse of the local weight function which multiplies the local operator $T_{00}$ in the corresponding entanglement hamiltonian (in the $c=1$ case), already derived through a CFT analysis. In §5 alternative constructions are discussed, including the one based on [4]. In §6 we draw our conclusions and mention some open problems. In appendix A some technical issues concerning the Williamson’s theorem are discussed. In appendix B we show that the contour constructed in §3.2 satisfies the properties (a), (b) and (c), in the weaker version discussed in §3.3. In appendix C we collect all the correlators employed in our numerical analysis.

2 Entanglement entropies in the harmonic lattices

In this section we briefly recall the standard method to compute the entanglement entropies in harmonic lattices, which involves the symplectic spectrum of the covariance matrix associated to the reduced density matrix $\rho_{A}$ of the subsystem [2, 4, 7]. Although throughout this manuscript we mainly consider the one-dimensional case, the method described in this section holds for a generic number of spatial dimensions.

The hamiltonian of the one-dimensional harmonic chain with $L$ lattice sites and nearest neighbour spring-like interaction reads

[TABLE]

where we impose either periodic ($\hat{q}_{L}=\hat{q}_{0}$ and $\hat{p}_{L}=\hat{p}_{0}$) or Dirichlet ($\hat{q}_{L}=\hat{q}_{0}=\hat{p}_{L}=\hat{p}_{0}=0$) boundary conditions. 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}$ are satisfied.

By arranging the operators $\hat{q}_{i}$ and $\hat{p}_{i}$ into the vector $\hat{\boldsymbol{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 as $[\hat{r}_{i},\hat{r}_{j}]=\textrm{i}J_{ij}$ (here we set $\hbar=1$), being $J$ the standard symplectic matrix

[TABLE]

where the square matrices $\boldsymbol{1}$ and $\boldsymbol{0}$ are the identity matrix and the matrix made by zeros respectively. Notice that $J^{\textrm{t}}=-J$ and $J^{2}=-\boldsymbol{1}$.

The real symplectic group $\textrm{Sp}(L)$ is made by the $2L\times 2L$ real matrices $M$ such that $MJM^{\textrm{t}}=J$. Given a symplectic matrix $M$, we have that $\textrm{det}(M)=1$ and $M^{\textrm{t}}\in\textrm{Sp}(L)$ [28]. The real symplectic group is very important for the analysis of the harmonic lattices because the linear transformations $\hat{\boldsymbol{r}}\to\hat{\boldsymbol{r}}^{\prime}=M\hat{\boldsymbol{r}}$ constructed through symplectic matrices $M$ preserve the canonical commutation relations.

A canonical rescaling of the variables allows to write the hamiltonian (8) as the hamiltonian of a free boson with mass $\omega$ discretised on a lattice with spacing $a=\sqrt{m/\kappa}$. In the numerical computations on the lattice presented in this manuscript we have set $m=\kappa=1$ without loss of generality. The continuum limit of this model provides the free scalar boson with mass $\omega$ in two spacetime dimensions. It is well known that in the massless case $\omega=0$ this quantum field theory is a CFT with central charge $c=1$.

In this manuscript we deal only with some Gaussian states associated to the harmonic chain (8), which can be either pure (e.g. the ground state) or mixed (e.g. the thermal state) (see [29] for a review). A Gaussian state for the harmonic chain (8) is completely characterised by the correlators $\langle\hat{r}_{i}\rangle$ (first moments) and $\langle\hat{r}_{i}\,\hat{r}_{j}\rangle$ (second moments). A shift in the first moments corresponds to a unitary transformation which preserves the Gaussian nature of the state. We will consider Gaussian states with vanishing first moments, which are fully described by the $2L\times 2L$ covariance matrix $\gamma\equiv\textrm{Re}\langle\boldsymbol{r}\,\boldsymbol{r}^{\textrm{t}}\rangle$ collecting all the second moments [2, 4, 7, 28].

The covariance matrix $\gamma$ is real and symmetric. It can be shown that the uncertainty principle coming from the canonical commutation relations implies that $\gamma+\textrm{i}J/2$ is positive semidefinite and that this result is equivalent to impose that $\gamma$ is positive definite and $\sigma_{k}\geqslant 1/2$ for its symplectic eigenvalues [30, 28]. Moreover, the covariance matrix corresponding to a pure Gaussian state satisfies the relation $(\textrm{i}J\,\gamma)^{2}=\tfrac{1}{4}\,\boldsymbol{1}$ [31]. It is straightforward to observe that the linear map $\hat{\boldsymbol{r}}\to\hat{\boldsymbol{r}}^{\prime}=M\hat{\boldsymbol{r}}$ characterised by the symplectic matrix $M$ induces on the covariance matrix the transformation $\gamma\to\gamma^{\prime}=M\gamma M^{\textrm{t}}$.

In this manuscript we are interested in the entanglement associated to spatial bipartitions where the whole harmonic lattice $A\cup B$ is in the Gaussian state characterised by the covariance matrix $\gamma$. Since we are considering the states with $\langle\hat{r}_{i}\rangle=0$, the group of linear transformations which preserves the Gaussian character of a state is $\textrm{Sp}(L)$.

The reduced density matrix $\rho_{A}$ associated to a spatial subsystem $A$ made by $\ell$ sites corresponds to a mixed state also when the state of the whole system $A\cup B$ is pure. For the harmonic chain (8) and its higher dimensional generalisations, the reduced density matrix $\rho_{A}$ remains Gaussian for any choice of the subsystem $A$. Because of this property, $\rho_{A}$ is characterised by the $2\ell\times 2\ell$ reduced covariance matrix $\gamma_{A}$ obtained by extracting from the covariance matrix $\gamma$ of the entire system the rows and the columns corresponding to the lattice sites which belong to $A$. Let us consider the reduced covariance matrices $\gamma_{A}$ with the following block structure

[TABLE]

where $Q$, $P$ and $R$ are the $\ell\times\ell$ correlation matrices corresponding to the subsystem $A$ (the subindex $A$ has been dropped for these matrices in order to lighten the forthcoming expressions). The reduced covariance matrix $\gamma_{A}$ is real, symmetric and positive definite. We remark that $Q$ and $P$ are symmetric and strictly positive, while $R$ is not constrained in general.

A crucial result for the quantitative analysis of the entanglement in the harmonic lattice is the Williamson’s theorem [32]. This theorem holds for a generic real, symmetric and positive matrices, but in this context we are interested in its consequences for the covariance matrices $\gamma_{A}$.

According to the Williamson’s theorem, a symplectic matrix $W\in\textrm{Sp}(\ell)$ exists such that

[TABLE]

where $D=\textrm{diag}(\sigma_{1},\dots,\sigma_{\ell})$ and $\sigma_{k}>0$. The set $\{\sigma_{1},\dots,\sigma_{\ell}\}$ is given by the positive square roots of the spectrum of the matrix $(\textrm{i}J\gamma_{A})^{2}$ and it is known as the symplectic spectrum of $\gamma_{A}$. The symplectic spectrum is invariant under canonical transformations.

The Williamson’s theorem has been proved in various ways [33], but in our analysis we will employ some steps of the proof found in [34]. In the appendix A we briefly discuss some issues related to the Williamson’s theorem, including also the part of the proof given by [34] that will be useful for our discussions.

The entanglement entropies in the harmonic lattices are obtained from the symplectic spectrum of $\gamma_{A}$ as [2, 4, 7]

[TABLE]

where $s(y)$ and $s_{n}(y)$ are the following analytic functions:

[TABLE]

and

[TABLE]

Notice that $\lim_{n\to 1}s_{n}(y)=s(y)$ for $y\geqslant 1/2$, as expected from the replica limit. This naturally leads us to adopt the notation $s_{1}(y)\equiv s(y)$. Notice that $s_{n}(y)\geqslant 0$ when $y\geqslant 1/2$, being $n\geqslant 1$. Thus, the symplectic eigenvalues $\sigma_{k}=1/2$ do not contribute to the entanglement entropies, i.e. the non vanishing terms in the sums (12) correspond only to the symplectic eigenvalues $\sigma_{k}>1/2$.

3 A contour for the entanglement entropies

In this section we describe the construction of our proposal for the contour function for the entanglement entropies in generic harmonic lattices. In §3.1, by adapting the observations of [26], we give a procedure to write contour functions satisfying (6) and (7) from real orthogonal matrices. In §3.2 we focus on the contour function obtained from an orthogonal matrix related to the Euler decomposition of the symplectic matrix $W$ occurring in the Williamson’s theorem (11). Then, in §3.3 we show that the contour function for the entanglement entropies constructed in §3.2 fulfils also a weaker formulation of the three constraints introduced in [26].

3.1 Contour functions from orthogonal matrices

In §2 we reviewed that the entanglement entropies in the harmonic lattice are given by the sums (12) over the symplectic spectrum of the covariance matrix $\gamma_{A}$ characterising the domain $A$ containing $\ell$ lattice sites.

Following [4, 26], we construct the contour function $s_{A}^{(n)}(i)$ satisfying (6) by assuming that we can associate $\ell$ real numbers $p_{k}(i)$ to every symplectic eigenvalue $\sigma_{k}$ (here $1\leqslant i\leqslant\ell$) such that

[TABLE]

Indeed, by inserting (15) into (12) and inverting the sums, it is immediate to recognize that the entanglement entropies can be written in the form (6) with the contour function given by

[TABLE]

where the function $s_{n}(x)$ is given by (13) for the entanglement entropy and by (14) for the Rényi entropies. We call mode participation function the function $p_{k}(i)$ in (16), as done in [4]. It encodes information about the contribution of the $i$-th site in $A$ to the term associated to the $k$-th eigenvalue $\sigma_{k}$ of the symplectic spectrum in the sums (12) providing the entanglement entropies. If we also require that $p_{k}(i)\geqslant 0$ for all $i$ and $k$, then Eq. (15) naturally leads to interpret the mode participation function as a set of probabilities. In this case the positivity condition (7) is guaranteed because $s_{n}(\sigma_{k})\geqslant 0$ for $\sigma_{k}\geqslant 1/2$. Any set of probabilities $p_{k}(i)$ provides a contour function (16) fulfilling the constraints (6) and (7). However, we are interested in finding mode participation functions which are based on the method underlying the computation of the entanglement entropies.

In order to identify some mode participation functions from (12), let us write the entanglement entropies as traces of suitable matrices. In particular, given a $2\ell\times 2\ell$ real orthogonal matrix $O\in O(2\ell)$, let us introduce

[TABLE]

being $D$ the diagonal matrix containing the symplectic spectrum of the reduced covariance matrix $\gamma_{A}$ (see Eq. (11)). Then, the entanglement entropies (12) can be written as

[TABLE]

where the analytic functions $s_{n}(x)$ are given by (13) and (14). In the last step we have employed (17), the cyclic property of the trace and the fact that the functions $s_{n}(x)$ are analytic.

The matrix $O$ defines a linear mapping sending the Williamson’s mode, labelled by $k$, into another set of modes, that we will label by the index $\alpha$. In order to isolate the contribution to the entanglement entropies (18) due to a specific mode characterised by a fixed value of $\alpha$, let us introduce a family $\{X^{(\alpha)},1\leqslant\alpha\leqslant\ell\}$ of orthogonal projectors. These operators are represented by $2\ell\times 2\ell$ matrices which are symmetric, semi-positive definite and they satisfy the property $X^{(\alpha)}X^{(\beta)}=\delta_{\alpha\beta}X^{(\alpha)}$ and $\sum_{\alpha=1}^{\ell}X^{(\alpha)}=\boldsymbol{1}$. In the base defined by the map $O$, the projector $X^{(\alpha)}$ can be written as $X^{(\alpha)}=\delta^{(\alpha)}\oplus\,\delta^{(\alpha)}$, being $\delta^{(\alpha)}$ the $\ell\times\ell$ matrix whose elements are $\delta^{(\alpha)}_{ab}=\delta_{a\alpha}\delta_{b\alpha}$.

Plugging the identity matrix written in the form $\sum_{\alpha=1}^{\ell}X^{(\alpha)}=\boldsymbol{1}$ into the argument of the trace occurring in the last step of (18) and employing the linearity of the trace, it is straightforward to realise that

[TABLE]

where

[TABLE]

By employing the expression $X^{(\alpha)}=\delta^{(\alpha)}\oplus\,\delta^{(\alpha)}$ into (20), together with the fact that $s_{n}(x)$ are analytic functions and the cyclic property of the trace, it becomes

[TABLE]

In order to write (22) in the form (16) and read the corresponding mode participation function, let us partition the orthogonal matrix $O$ introduced in (17) in four $\ell\times\ell$ blocks

[TABLE]

Plugging this block partitioned matrix into (22), for the contour function we find

[TABLE]

Writing explicitly the four terms occurring in this expression, we find that (24) becomes

[TABLE]

with the mode participation function given by

[TABLE]

which is positive by construction.

In order to check that the mode participation function (26) satisfies (15), let us start from the orthogonality condition $OO^{\textrm{t}}=\boldsymbol{1}$ for the block partitioned matrix (23), i.e.

[TABLE]

By considering the $k$-th element (with $1\leqslant k\leqslant\ell$) and the $(\ell+k)$-th element along the diagonal of (27), we obtain respectively

[TABLE]

which tell us that (26) fulfils the constraint $\sum_{\alpha=1}^{\ell}p_{k}(\alpha)=1$ for any integer $k\in[1,\ell]$.

By employing the orthogonality condition $O^{\textrm{t}}O=\boldsymbol{1}$ instead of (27) and following similar steps, one finds that (26) satisfies also the further property $\sum_{k=1}^{\ell}p_{k}(\alpha)=1$ for any integer $\alpha\in[1,\ell]$.

Since we are interested in constructing contour functions, let us consider those cases where the index $\alpha$ labels the sites of $A$, namely $\alpha=i\in A$, according to the notation adopted throughout this manuscript. Henceforth we will employ the projector $X^{(i)}$ corresponding to the $i$-th site of the region $A$.

Summarising, the entanglement entropies (12) can be written in the form (6) with the following contour function

[TABLE]

which can be expressed as in (16) with the mode participation function given by

[TABLE]

in terms of the elements of the orthogonal matrix (23). Since the mode partition function (30) is positive and fulfils the constraint (15), the corresponding contour function (29) satisfies the properties (6) and (7).

We find it worth remarking that, while the above discussion is based on the fact that the matrix $O$ in (17) is orthogonal, a canonical transformation is implemented by a symplectic matrix. By requiring that $O\in O(2\ell)\cap\textrm{Sp}(\ell)$, further constraints for the blocks in (23) coming from the condition $OJO^{\textrm{t}}=J$ can be employed. An explicit example belonging to this class is considered in the next subsection.

3.2 A proposal based on the Williamson’s theorem and the Euler decomposition

The discussion in §3.1 allows to conclude that, given an orthogonal matrix (23), we can construct the mode participation function (30) and, consequently, the corresponding contour function (16) satisfies (6) and (7).

Some particular orthogonal matrices are more relevant for the physics of our problem, which is encoded in the Gaussian reduced density matrix $\rho_{A}$. Being $\rho_{A}$ fully described by the reduced covariance matrix $\gamma_{A}$ in our cases, let us focus on the orthogonal matrices related to $\gamma_{A}$. In particular, since $p_{k}(i)$ provides the contribution of the $i$-th site in $A$ to the term associated to the $k$-th symplectic eigenvalue in (12), we find it worth looking for a meaningful orthogonal matrix within the linear transformation which relates the canonical variables $(\hat{q}_{i},\hat{p}_{i})$, labelled by the index $i$ of the lattice sites, to the canonical variables labelled by the index $k$ associated to the symplectic spectrum. This particular canonical transformation is implemented by the real symplectic matrix $W$ associated to $\gamma_{A}$ through the Williamson’s theorem (11).

The Euler decomposition (also known as Bloch-Messiah decomposition) [35] of the real symplectic matrix $W$ introduced in (11) reads

[TABLE]

where $\chi_{j}\geqslant 0$ and the real matrices $K_{\textrm{\tiny L}}$ and $K_{\textrm{\tiny R}}$ are symplectic and orthogonal. The non-uniqueness of the decomposition (31) is due only to the freedom to order the elements along the diagonal of $\chi$. The set containing the matrices of the form given by $E$ is a subgroup of $\textrm{Sp}(\ell)$ corresponding to the single-mode squeezing operations. Combining the polar decomposition of the real symplectic matrix $W$ and its Euler decomposition (31), we find that it can be written as follows

[TABLE]

where

[TABLE]

The real matrix $K$ is symplectic and orthogonal, while the real matrices $E_{\textrm{\tiny R}}$ and $E_{\textrm{\tiny L}}$ are symplectic, symmetric and positive definite. The orthogonal matrix $K$ is obtained by removing the squeezing matrix $E$ in the Euler decomposition of $W$, which is the factor making $W$ non-orthogonal. In the following we will employ the polar decomposition $W=KE_{\textrm{\tiny R}}$ and this factorisation is unique [36].

Symplectic matrices which are also orthogonal form a subgroup of $\textrm{Sp}(\ell)$ which is isomorphic to the group of the $\ell\times\ell$ unitary matrices (see e.g. Proposition 2.12 of [28]). These matrices preserve the trace of the covariance matrix and they correspond to passive unitary transformations for its density matrix. Instead, when a symplectic transformation is not orthogonal, the trace of the covariance matrix changes and the corresponding unitary transformations for the density matrix are called active [29].

By employing the properties of the matrices highlighted above, we can easily write $E_{\textrm{\tiny L}}$ and $E_{\textrm{\tiny L}}$ in terms of $W$ as follows

[TABLE]

Also the orthogonal matrix $K$ can be written in terms of $W$ by using $K=E_{\textrm{\tiny L}}^{-1}\,W=WE_{\textrm{\tiny R}}^{-1}$ and the relations in (34). The result reads

[TABLE]

The factorisations (31) and (32) hold for any real symplectic matrix. In our analysis we are interested in the real symplectic matrix entering in the Williamson’s theorem. Among the orthogonal matrices occurring in the decompositions (31) and (32), we think that the matrix $K$ is the most natural one to consider in order to construct a mode participation function. Plugging the polar decomposition $W=K\,E_{\textrm{\tiny R}}$ into (11), we find that

[TABLE]

This is an explicit realisation of (17) with $\Phi_{A}=E_{\textrm{\tiny R}}^{-1}\,\gamma_{A}\,E_{\textrm{\tiny R}}^{-1}$ and $O=K$. Notice that in this case $K$ is also symplectic. The relation (36) tells us that $K$ is the orthogonal matrix which diagonalises the symmetric matrix $E_{\textrm{\tiny R}}^{-1}\,\gamma_{A}\,E_{\textrm{\tiny R}}^{-1}$.

The next step consists in writing $E_{\textrm{\tiny R}}^{-1}$ in (36) in terms of the covariance matrix $\gamma_{A}$. This can be done by employing some steps of the constructive proof of the Williamson’s theorem found in [34], which have been briefly recalled in the appendix A. In particular, given the real, symmetric and positive definite covariance matrix $\gamma_{A}$, one introduces the following antisymmetric matrix

[TABLE]

where $J$ is the standard $2\ell\times 2\ell$ symplectic matrix (9). Being the matrix $\hat{\gamma}_{A}$ antisymmetric, an orthogonal matrix $\widetilde{O}$ exists such that

[TABLE]

where $D=\textrm{diag}(\sigma_{1},\dots,\sigma_{\ell})$ is the diagonal matrix containing the symplectic spectrum of $\gamma_{A}$ introduced in (11). By extracting $\hat{\gamma}_{A}$ from (38), we find that

[TABLE]

which is positive definite, being $MM^{\textrm{t}}$ positive definite for any real invertible matrix $M$. From (39) one obtains $|\hat{\gamma}_{A}|$, which reads

[TABLE]

This relation tells us that the symmetric and positive definite matrix $|\hat{\gamma}_{A}|$ is diagonalised by the orthogonal matrix $\widetilde{O}$ and its spectrum coincides with the symplectic spectrum.

It is worth noticing that (40) provides a realisation of (17) different from (36). Indeed, $\Phi_{A}=|\hat{\gamma}_{A}|$ and $O=\widetilde{O}$ in this case. Let us remark that $\widetilde{O}$ is not necessarily symplectic. In §5 the contour function associated to the orthogonal matrix $\widetilde{O}$ will be discussed.

In order to express $E_{\textrm{\tiny R}}^{-1}$ in (36) in terms of $\gamma_{A}$, we employ a crucial step occurring in the proof of the Williamson’s theorem found in [34] (see also the appendix A), where the symplectic matrix $W$ satisfying (11) is constructed as follows

[TABLE]

By using this expression, the second relation in (34) becomes

[TABLE]

where in the last step (39) has been employed. Thus, $E_{\textrm{\tiny R}}^{-1}$ reads

[TABLE]

This result allows us to write the l.h.s. of (36) from the covariance matrix $\gamma_{A}$.

By specialising (18) to the explicit case given by (36), we can write the entanglement entropies (12) as follows

[TABLE]

where the matrix $E_{\textrm{\tiny R}}^{-1}$ is the function of the covariance matrix $\gamma_{A}$ in (43). The corresponding contour function is obtained by specialising the expression (29) to the case given by (36) and the result reads

[TABLE]

Once the orthogonal matrix $K$ introduced in (33) is written in its block form

[TABLE]

the mode participation function corresponding to the contour function (45) is obtained by specialising (30) to the case $O=K$, namely

[TABLE]

which fulfils the constraint (15) because of the orthogonality of $K$, as shown in general in (27) and (28) for the orthogonal matrix $O$.

The symplectic condition $KJK^{\textrm{t}}=J$ for the block matrix (46) tells us that $U_{K}Y_{K}^{\textrm{t}}$ and $V_{K}Z_{K}^{\textrm{t}}$ are symmetric matrices and that $U_{K}V_{K}^{\textrm{t}}-Y_{K}Z_{K}^{\textrm{t}}=\boldsymbol{1}$. These relations do not simplify (47) in the general case.

The expressions in (45) and (47) are the main result of this section and they provide our proposal for the contour function in generic harmonic lattices. This proposal fulfils the constraints (6) and (7). In the next subsection we prove that it also satisfies three further requirements which correspond to a weaker version of the properties introduced in [26] for the contour functions for the entanglement entropies.

In the remaining part of this section, we discuss the cases where the reduced covariance matrix (10) is block diagonal, namely $R=\boldsymbol{0}$ and therefore $\gamma_{A}=Q\oplus P$, being $Q_{ij}=\langle\hat{q}_{i}\hat{q}_{j}\rangle$ and $P_{ij}=\langle\hat{p}_{i}\hat{p}_{j}\rangle$ the $\ell\times\ell$ correlation matrices restricted to the subsystem $A$. All the examples considered in our numerical analysis belong to this class.

When $\gamma_{A}$ is block diagonal, we have that $(\textrm{i}J\gamma_{A})^{2}=(PQ)\oplus(QP)$. This implies that the symplectic spectrum can be found by first computing the spectrum of $QP$ (which is equal to the spectrum of its transpose $PQ$) and then taking its positive square root. The antisymmetric matrix $\hat{\gamma}_{A}$ defined in (37) simplifies to

[TABLE]

By specifying (39) and (40) to this case, we find that also $|\hat{\gamma}_{A}|$ is block diagonal

[TABLE]

Given a block diagonal covariance matrix $\gamma_{A}$, the symplectic matrix $W$ occurring in the Williamson’s theorem (11) is block diagonal as well. Thus, also the factorisations (31) and (32) for $W$ are made by block diagonal matrices. In particular, considering the symplectic, symmetric and positive definite matrix $E_{\textrm{\tiny R}}^{-1}=\Xi\oplus\Pi$, we have that $\Xi$ are $\Pi$ are symmetric and positive definite matrices. The fact that $E_{\textrm{\tiny R}}^{-1}$ is symplectic becomes the condition $\Xi\,\Pi^{\textrm{t}}=\boldsymbol{1}$, which tells us that $\Xi$ and $\Pi$ are not orthogonal. The matrices $\Xi$ and $\Pi$ can be written in terms of $Q$ and $P$ by plugging (49) into (43) and exploiting the block diagonal structure of $\gamma_{A}$. The result reads

[TABLE]

Also the matrix $E_{\textrm{\tiny R}}^{-1}\gamma_{A}E_{\textrm{\tiny R}}^{-1}$ occurring in the contour function (45) becomes block diagonal in this case. In particular, it reads

[TABLE]

where the blocks on the diagonal are complicated functions of $Q$ and $P$ provided by (50).

Since the matrices involved in the factorisations (31) and (32) are block diagonal, the orthogonal and symplectic $K$ in (46) becomes $K=U_{K}\oplus V_{K}$, namely $Y_{K}=Z_{K}=\boldsymbol{0}$. The orthogonality condition for this $K$ is equivalent to require that both $U_{K}$ and $V_{K}$ are orthogonal, while the symplectic condition leads to the relation $U_{K}V_{K}^{\textrm{t}}=\boldsymbol{1}$. Combining these observations, we can conclude that $U_{K}=V_{K}$. By employing this result and (51), the matrix relation (36) simplifies to

[TABLE]

which tell us that the orthogonal matrix $U_{K}$ diagonalises the symmetric matrix $\Xi\,Q\,\Xi=\Pi\,P\,\Pi$. This observation and (50) allow us to compute $U_{K}$ from the correlation matrices $Q$ and $P$.

The mode participation function for $\gamma_{A}=Q\oplus P$ is obtained by first specialising (47) to this simpler case and then employing $U_{K}=V_{K}$. The result reads

[TABLE]

Henceforth we focus on the contour function for the entanglement entropies in the harmonic chain (8), setting $m=\kappa=1$ without loss of generality. Since all the examples explored in this manuscript have a block diagonal reduced covariance matrix $\gamma_{A}=Q\oplus P$, the mode participation function is given by (53). The boundary conditions imposed at the endpoints of the harmonic chain are crucial in the computation of $p_{k}(i)$ because they determine the correlators entering in the matrices $Q$ and $P$. All the correlators employed in our explicit examples have been collected in the appendix C. In the numerical analysis the symplectic eigenvalues $\{\sigma_{k},\,1\leqslant k\leqslant\ell\}$ have been arranged in decreasing order in terms of the label $k$, namely $\sigma_{1}\geqslant\sigma_{2}\geqslant\dots\geqslant\sigma_{\ell}$.

The correlators providing the matrix elements of $Q$ and $P$ for the harmonic chain with periodic boundary conditions in the vacuum are given by (116) and in the thermodynamic limit $L\to\infty$ they become (120) and (121). Considering this regime, in Fig. 1 we show the mode participation function (53) for a single interval $A$ made by $100$ sites when the mass parameter is very small ($\omega=10^{-8}$). We remark that we cannot set $\omega=0$ for periodic boundary conditions because the correlator $\langle\hat{q}_{i}\hat{q}_{j}\rangle$ diverges when $\omega\to 0$. This divergence is due to the occurrence of the zero mode, which is a consequence of the invariance under translations.

In Fig. 1 the peaks of the mode participation function close to the endpoints of the interval and for small values of $k$, which correspond to the large symplectic eigenvalues, tell us that a large contribution to the entanglement entropies comes from the sites close to the endpoints of the interval. We find it worth noticing that the mode participation function is symmetric with respect to the center of the interval and also the occurrence of a profile $\tilde{k}(i)$ which delimits the modes providing a non vanishing contribution to the mode participation function, namely $p_{k}(i)\simeq 0$ for $k>\tilde{k}(i)$ for the $i$-th site.

In Fig. 2 we have considered the dependence on the mass $\omega$ in $p_{k}(i)$ for the same setup described for Fig. 1. The density plots for the mode participation function correspond to different values of $\omega$: the top left panel contains the same data of Fig. 1, while in the other panels $\omega$ is larger. As the mass increases, the set of modes providing a non vanishing contribution to $p_{k}(i)$ at fixed $i$ is sharper, becoming very localised for large values of $\omega$ (see the bottom right panel). In the massive regime, where an area law occurs, the contribution of the large symplectic eigenvalues, which correspond to small values of $k$ and provide the largest part of the entanglement entropies, is localised close to the endpoints of the interval, as one observes by comparing the top panels of Fig. 2. The symmetry of the profile with respect to the center of the interval is not influenced by the value of the mass, as expected. As for the profile of the function $\tilde{k}(i)$, it clearly depends on $\omega$.

In Fig. 3 we show the contour function for the entanglement entropy for increasing values of the mass $\omega$. As the mass increases, the power law behaviour of the contour function in the massless case becomes exponential in the regions between an endpoint and the center of the interval. In the right panel only the first half of the interval is considered and the collapse of the data tells us that, in the massive regime, there is a domain in $A$ where $s_{A}(i)\propto e^{-b\omega i}$ with the parameter $b>0$ independent of $\ell$.

3.3 Three further requirements

The constraints (6) and (7) are the minimal requirements for the contour function for the entanglement entropies. In [26] other three reasonable properties have been introduced. In the following, we show that (45) satisfies a weaker version of these further constraints.

Before discussing the statements of these properties, let us motivate the fact that we restrict our analysis to the special class of canonical transformations corresponding to symplectic matrices which are also orthogonal.

Given a canonical transformation implemented by a symplectic matrix $M$, the covariance matrix transforms as $\gamma_{A}\to\gamma^{\prime}_{A}=M\gamma_{A}M^{\textrm{t}}$. From this transformation rule and the Williamson’s theorem (11), it is straightforward to realise that \gamma_{A}^{\prime}=(W^{\prime})^{\textrm{t}}\big{(}D\oplus D\big{)}W^{\prime}, where for the symplectic matrix $W^{\prime}$ we have $W^{\prime}=WM^{\textrm{t}}$. This shows that the symplectic spectrum is invariant under canonical transformations. Plugging the polar decompositions $W=KE_{\textrm{\tiny R}}$ and $W^{\prime}=K^{\prime}E^{\prime}_{\textrm{\tiny R}}$ into the relation $W^{\prime}=WM^{\textrm{t}}$, one finds

[TABLE]

The last step is obtained by employing the identity matrix in the form $\boldsymbol{1}=M^{-1}M$, in order to recognise the symmetric matrix $ME_{\textrm{\tiny R}}M^{\textrm{t}}$, which is also positive definite, being $E_{\textrm{\tiny R}}$ positive definite. The last expression in (54) does not provide a polar decomposition for $W^{\prime}$ because the symplectic matrix $KM^{-1}$ is not orthogonal for a generic symplectic matrix $M$, unless we restrict to the class of symplectic matrices $M$ which are also orthogonal. In this case, from (54) we have

[TABLE]

consistently with the properties of the matrices entering in the polar decompositions of $W^{\prime}$ and $W$. By employing (55), for the matrix $E_{\textrm{\tiny R}}^{-1}\gamma_{A}E_{\textrm{\tiny R}}^{-1}$ occurring in the contour function (45) we find

[TABLE]

which tells us that, if we restrict to the class of canonical transformations implemented by $M\in\textrm{Sp}(\ell)\cap O(2\ell)$, then the matrix $E_{\textrm{\tiny R}}^{-1}\gamma_{A}E_{\textrm{\tiny R}}^{-1}$ transforms like $\gamma_{A}$.

In the following (see (a), (b) and (c) below) we introduce a weaker version of the three requirements proposed in [26] beside (6) and (7) by considering only the transformations characterised by the subgroup $\textrm{Sp}(\ell)\cap O(2\ell)$ instead of the entire symplectic group $\textrm{Sp}(\ell)$.

(a) Spatial symmetry. If $\rho_{A}$ is invariant under a transformation relating the sites $i$ and $j$ in the subsystem $A$, then $s_{A}^{(n)}(i)=s_{A}^{(n)}(j)$.

This requirement is due to the possible occurrence of a spatial symmetry, which depends both on the underlying lattice model and on the choice of the spatial region $A$. Typical examples of spatial symmetries could be related to the invariance under translations, rotations or space reflections.

In order to formulate the remaining properties, we need to introduce also the contour $s_{A}^{(n)}(G)$ of a subregion $G\subseteq A$ as follows

[TABLE]

In the special case of $G=A$, from (6) and (57) we find $s_{A}^{(n)}(A)=S_{A}^{(n)}$. The contour $s_{A}^{(n)}(G)$ is clearly additive: for any two non intersecting spatial subsets $G\subsetneq A$ and $\tilde{G}\subsetneq A$ we have $s_{A}^{(n)}(G\cup\tilde{G})=s_{A}^{(n)}(G)+s_{A}^{(n)}(\tilde{G})$. Moreover, the contour $s_{A}^{(n)}(G)$ is monotonous, i.e. for $G\subseteq\tilde{G}\subseteq A$ the inequality $s_{A}^{(n)}(G)\leqslant s_{A}^{(n)}(\tilde{G})$ holds.

(b) Invariance under local unitary transformations. Given a system in the state characterised by the density matrix $\rho$ and a unitary transformation $U_{G}$ acting non trivially only on $G\subseteq A$, denoting by $\rho^{\prime}$ the state of the system after such transformation, the same contour $s_{A}^{(n)}(G)$ should be found for $\rho$ and $\rho^{\prime}$.

The property (b) is motivated by the expectation that the contribution of a subregion $G$ to the entanglement should not be modified by a change of basis restricted to $G$.

(c) A bound. Given a system in the pure state $|\Psi\rangle$ and the bipartition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, let us assume that the further decompositions $\mathcal{H}_{A}=\mathcal{H}_{\Omega_{A}}\otimes\mathcal{H}_{\bar{\Omega}_{A}}$ and $\mathcal{H}_{B}=\mathcal{H}_{\Omega_{B}}\otimes\mathcal{H}_{\bar{\Omega}_{B}}$ lead to the following factorisation of the state

[TABLE]

Considering a subregion $G\subseteq A$ such that $\bigotimes_{i\in G}\mathcal{H}_{i}\subseteq\mathcal{H}_{\Omega_{A}}$, we must have that

[TABLE]

where $S^{(n)}(\Omega_{A})$ are the entanglement entropies corresponding to the reduced density matrix $\rho_{\Omega_{A}}$, obtained by tracing over the degrees of freedom of $\mathcal{H}_{\bar{\Omega}_{A}}\otimes\mathcal{H}_{B}$.

The upper bound discussed in the requirement (c) can be also presented as a lower bound. This observation can be examined by adapting the corresponding discussion made in [26] to our case in a straightforward way. We refer the interested reader to [26] for a more detailed discussion on the motivations leading to this property.

In appendix B we have shown that the contour function (45) fulfils also the constraints (a), (b) and (c) for the restricted class of transformations characterised by symplectic and orthogonal matrices. It would be interesting either to extend this analysis to the whole group of the symplectic matrices or to find physical motivations leading to this restriction.

As emphasized in [26], let us remark that the constraints (6) and (7) together with the three further requirements presented in this subsection (even in the version formulated in [26]) do not characterise the contour function for the entanglement entropies in a unique way. Inequivalent expressions for $s_{A}^{(n)}(i)$ satisfying the above five properties could be found. It would be very interesting to list a set of constraints which identify a unique construction of the contour function for the entanglement entropies.

4 Massless regime and entanglement hamiltonians in 2d CFTs

In this section we study the contour function for the entanglement entropies proposed in §3.2 in the massless regime of the harmonic chains where either periodic or Dirichlet boundary conditions are imposed. For some particular configurations of the subsystem $A$ made by a single interval, the corresponding entanglement hamiltonians in two dimensional CFTs suggest a candidate for the continuum limit of the contour function for the entanglement entropies.

In §4.1 we briefly review the static cases where the entanglement hamiltonians can be written as an integral over the domain $A$ of the $T_{00}$ component of the energy-momentum tensor multiplied by a suitable local weight function. We focus our attention on the results obtained in two-dimensional CFTs, where the weight function entering in the entanglement hamiltonian provides a candidate for the continuum limit of the contour function for the entanglement entropies. In §4.2 we explicitly study the static cases shown in the top and middle panels of Fig. 4, where this analysis can be applied by performing a comparison with the corresponding numerical results from the lattice. In §4.3 we consider the contour function for the entanglement entropy in the massless regime when $A=A_{1}\cup A_{2}$ is made by two disjoint intervals in the infinite line (see the bottom panel in Fig. 4). We emphasize that in this case a reliable candidate for the contour function in the continuum limit coming from a CFT analysis is not known.

4.1 Weight function in some entanglement hamiltonians in 2d CFTs

The entanglement hamiltonian $K_{A}$ defined in (5) provides the reduced density matrix of the subsystem; therefore in principle it contains more information than the entanglement entropies. For free models on the lattice the entanglement hamiltonians have been studied in [3]. In relativistic quantum field theories the entanglement hamiltonian is usually non-local but in certain cases it can be written in terms of integrals of local operators multiplied by suitable local weight functions. In a spacetime with a generic number of dimensions, this happens for very few known examples. Instead, for two-dimensional CFTs the large symmetry allows to enlarge the class of the configurations whose entanglement hamiltonian can be written in a local way as explained above, by including both static and time dependent examples. In the following we will focus on the static cases.

The most important example where the entanglement hamiltonian can be written as an integral of a local operator multiplied by a suitable weight function is the vacuum state of a QFT in flat $d$-dimensional Minkowski space ($\mathbb{R}^{d}$ in euclidean space) with the spatial subsystem $A$ given by the half space $x_{1}>0$ and $x_{i}\in\mathbb{R}$ for $2\leqslant i\leqslant d-1$. For this configuration the following crucial result due to Bisognano and Wichmann holds [21]

[TABLE]

which tells us that $K_{A}$ is the generator of the Lorentz boosts along the $x_{1}$-direction, or of the euclidean rotations around $x_{1}=0$ in the euclidean space.

When the QFT is a $d$-dimensional CFT, a conformal transformation can be employed to map (60) into the entanglement hamiltonian $K_{A}$ for a ball $A$ with radius $R$. The result of this mapping reads [22]

[TABLE]

In the special case of $d=2$, the subsystem $A$ is a single interval $A=(-R,R)$ in the infinite line (see Fig. 4, top left panel) at zero temperature and (61) gives the corresponding entanglement hamiltonian. The main difference between (60) and (61) is the local weight function multiplying $T_{00}$ within the integrand.

For two-dimensional CFTs, the entanglement hamiltonian corresponding to other interesting configurations where the subsystem $A$ is a single interval can be written by conformal mapping the result (60) in a suitable way. In [24] these cases have been studied within a unifying framework that we briefly review in the following.

Consider the two-dimensional euclidean spacetime (which can also have boundaries) describing the physical system. The subsystem $A$ is an interval that can be infinite or reach the boundary of this euclidean space. A suitable regularization procedure consists in removing small discs of radius $\epsilon\ll 1$ around the endpoints of the interval $A$ [8]. When the subsystem $A$ touches the boundary of the spacetime, then only one endpoint has to be regularised in this way. In [24] it has been argued that, after this regularisation, whenever the regularised domain can be conformally mapped into an annulus the entanglement hamiltonian corresponding to the initial configuration can be identified with the generator of the translations around the annulus along the direction orthogonal to the direction which connects the two boundaries of the annulus. Denoting by $w$ the complex variable parameterising the annulus and by $z$ the complex coordinate of the original domain where the regularisation discs have been removed, we have that $w=f(z)$.

All the examples that we consider in this manuscript with $A$ made by a single interval (see §4.2) are static and fall into the class of configurations just described. They are shown in the top and middle panels of Fig. 4. For these configurations, we can map back the generator of the translations around the annulus into the $z$-domain, finding that the entanglement hamiltonian can be written as follows [24]

[TABLE]

The analysis of [24] allows to relate the local weight function $1/f^{\prime}(x)$ multiplying $T_{00}(x)$ to the corresponding entanglement entropies. The result reads

[TABLE]

where $c$ is the central charge of the underlying CFT and $\mathcal{W}$ is the width of the annulus in the $w$-domain. The width $\mathcal{W}$ can be computed from the above mentioned conformal transformation $f(z)$ mapping the $z$-domain into the annulus as follows

[TABLE]

where we have introduced the notation $A_{\epsilon}$ to denote the subsystem $A$ after the removal of the small discs around the endpoints. The width $\mathcal{W}$ in (64) is divergent as $\epsilon\to 0$. The simplest example is a single interval $A=(0,\ell)$ in the infinite line at zero temperature (see Fig. 4, top left panel). In this case the integration domain in (64) is $A_{\epsilon}=(\epsilon,\ell-\epsilon)$ and one finds that $\mathcal{W}=2\log(\ell/\epsilon)$.

The constant $C$ in (63) is subleading as $\epsilon\to 0$ and it is related to the boundary entropy [37] associated to the conformally invariant boundary conditions imposed at the boundaries of the $z$-domain, which include also the boundaries due to the small discs removed during the regularization procedure [38, 24].

The expressions in (63) and (64), which provide the entanglement entropies of those configurations where (62) holds, suggest a natural candidate for the contour function for the entanglement entropies in the scaling limit. Indeed, plugging (64) into (63) and neglecting terms which are infinitesimal as $\epsilon\to 0$, it is straightforward to consider the following function

[TABLE]

where $\ell$ is the length of the interval $A$. We remark that $f(x)$ depends both on the euclidean spacetime defining the model and on the configuration of the interval.

In the following we restrict our analysis to the massless scalar field in two spacetime dimensions, which has $c=1$. Our main goal is to compare the expression (65) with $c=1$ to the scaling limit of the contour function for the entanglement entropies proposed in §3.2 in the special cases of one-dimensional harmonic chains in the massless regime. In particular, we are interested in the space dependent term containing $f^{\prime}(x)$ in (65) and not in the vertical shift characterised by the constant $C$, which is influenced by non universal features (see [39] for an exact computation of this term in a specific spin chain).

4.2 Single interval

The first case we consider is a single interval $A$ of length $\ell$ in the infinite line (top left panel of Fig. 4), when the whole system is in the ground state. Because of the invariance under translations, we can fix the origin in the first endpoint of the interval and therefore $A=(0,\ell)$. The map transforming the configuration obtained by removing small discs of radius $\epsilon$ around the endpoints of $A$ into the annulus is given by

[TABLE]

Taking the derivative of this expression and plugging the result into (65), for $c=1$ we get

[TABLE]

In Fig. 5 we show the contour function for the entanglement entropies constructed in §3.2 when $A$ is a single interval of various lengths in the infinite line. The data have been found by considering the harmonic chain (8) in the thermodynamic regime, whose correlators are (120) and (121). The data are obtained for very small but non vanishing mass $\omega$. Indeed, this parameter cannot be set to zero because the correlator $\braket{\hat{q}_{i}\hat{q}_{j}}$ diverges in this limit. This is due to the occurrence of the zero mode, which is a consequence of the invariance under translations of the model.

In Fig. 5 the panels on the right have the same vertical axis of the corresponding ones on the left, but the parameterisations of the horizontal axis are different. In particular, we have employed the function suggested by the CFT result (67), in order to show more clearly the behaviour of the lattice data nearby the endpoints of the interval. This way to display the data will be adopted also in other figures corresponding to the subsequent examples.

The dashed curves in Fig. 5 are obtained from the CFT formula (67), where the constant $C$ is not universal. In order to fix $C$ through a method applicable to all the examples considered throughout this subsection, we can employ (65). In particular, given the lattice points for the contour function, we impose that \sum_{i=1}^{\ell_{\textrm{\tiny max}}}s_{A}^{(n)}(i)=C+\tfrac{1}{12}(1+1/n)\sum_{i=1}^{\ell_{\textrm{\tiny max}}}f^{\prime}(i-1/2)\big{|}_{\ell_{\textrm{\tiny max}}}, where $\ell_{\textrm{\tiny max}}$ is the size of the largest interval considered in the corresponding numerical analysis and the shifted argument for $f^{\prime}$ is introduced because this function diverges at the endpoints of the interval. A deeper analysis that we leave for future work could lead to an insightful method to fix the constant $C$ in (65) from the numerical data.

A very good collapse of the numerical data corresponding to different values of $\ell$ is observed in the top panels of Fig. 5, already for small intervals. Few points close to the endpoints of the interval, where the contour function diverges, have not been shown because non universal features due to the lattice are expected very close to the endpoints of the interval. The agreement between the CFT expression (67) and the lattice data is very good nearby the endpoints of the interval and gets worse around the center. This agrees with the expectation that the universal part of the entanglement entropies is determined by the regions close to the endpoints of the interval. The disagreement between the CFT formula and the lattice data around the center is due to non universal contributions.

In the next example the whole system is on a semi-infinite line, therefore boundary conditions must be imposed on the physical boundary. The subsystem $A=(0,\ell)$ is a single interval of length $\ell$ placed at the beginning of the semi-infinite line (middle left panel of Fig. 4). The regularizing procedure requires to remove a small disc of radius $\epsilon$ only around the second endpoint of $A$. Mapping the resulting configuration into the annular geometry, it is worth remarking that in this case the boundary conditions on the two boundaries of the annulus may be different. The conformal map reads

[TABLE]

Given this function, for this example the expression (65) specialised to $c=1$ becomes

[TABLE]

As for the numerical analysis on the lattice, for this configuration we consider the harmonic chain defined on a segment with Dirichlet boundary conditions imposed to both its endpoints and then we take the thermodynamic limit. In this case the invariance under translations is broken because of the presence of the physical boundary; therefore the zero mode does not occur in the correlators (124) and (125), which become (128) and (129) respectively in the thermodynamic limit. This fact allows us to set $\omega=0$ in the numerical analysis.

In Fig. 6 we show the contour function for the entanglement entropies described in §3.2 with $\omega=0$. Since we have imposed Dirichlet boundary conditions, the curves for the contour function obtained from the lattice data pass through the origin, as highlighted in the insets of the left panels.

The dashed curves in Fig. 6 (only their positive part are shown) correspond to the CFT formula (69) where the constant $C$ has been fixed as explained above by employing the map (68). The CFT curves nicely reproduce the divergent behaviour of the lattice data nearby the endpoint of the interval, while they cannot capture the lattice data for the contour function around the boundary, which introduces non universal features.

The above analysis can be applied also for a single interval $A$ of length $\ell$ in the infinite line when the whole system is in a thermal state at temperature $T=1/\beta$ (top left panel of Fig. 4). Setting the interval in $A=(0,\ell)$ and removing a small disc around both the endpoints of $A$, the conformal map relating this configuration to the annular domain reads

[TABLE]

Given this expression and fixing $c=1$, one finds that (65) becomes

[TABLE]

which has been written in a form highlighting the dependence on the ratios $x/\ell\in(0,1)$ and $\ell/\beta$ in the universal term.

In Fig. 7 we show the contour function for the entanglement entropy of a single interval of length $\ell$ in the harmonic lattice at finite temperature $T\equiv 1/\beta$ and in the thermodynamic regime. The contour function has been evaluated by employing (16) and (53), with the correlators (122) and (123). Since this model is invariant under translations, the zero mode occurs in the correlators and this prevents us from setting $\omega=0$. The data shown in the figure are obtained for very small non vanishing mass $\omega$.

The dashed curve in Fig. 7 corresponds to the CFT formula (71) specialised to $n=1$, with the vertical shift given by the constant $C$ fixed through the function (70), as explained in the first example of this subsection.

In the remaining two examples we consider finite length systems in their ground state. Different boundary conditions (either periodic or Dirichlet) are imposed.

An interesting configuration to study is given by a single interval $A$ of length $\ell$ in a spatial circle whose length is $L>\ell$ (top right panel of Fig. 4). Setting the system in the finite segment $(0,L)$ and fixing the position of the interval to $A=(0,\ell)$, periodic boundary conditions are imposed at $x=0$ and $x=L$. The system is invariant under translations along the circle. After the removal of the small discs at the endpoints of $A$ in the euclidean spacetime describing this system (an infinite cylinder), the conformal map which sends the resulting domain into the annulus reads

[TABLE]

Applying the formula (65) for this map and for $c=1$, one obtains

[TABLE]

This expression has been written in a form highlighting the fact that the universal term in the r.h.s. is a function of the two ratios $x/\ell\in(0,1)$ and $\ell/L\in(0,1)$.

In Fig. 8 we plot the contour function for the entanglement entropy discussed in §3.2 for a single interval of length $\ell$ in a periodic harmonic chain of finite length $L$ in its ground state. The correlators employed for this numerical analysis have been written in (116). Notice that $\langle\hat{q}_{i}\hat{q}_{j}\rangle$ diverges as $\omega\to 0$ because of the occurrence of the zero mode $k=0$; therefore the mass $\omega$ must be very small but non vanishing. The dashed curve in Fig. 8 has been obtained from the CFT formula (73) in the special case of $n=1$, with the constant $C$ fixed by employing the function in (73) as explained for the previous cases.

Our last example in the class of configurations whose entanglement hamiltonian in the continuum can be written in the form (62), is given by a system in a finite segment $(0,L)$ where the same boundary conditions are imposed at both its endpoints. The subsystem $A=(0,\ell)$ is a single interval of length $\ell<L$ sharing an endpoint with the entire system (middle right panel of Fig. 4). By adapting the procedure described above to this case through the map given in [24], one finds that

[TABLE]

In Fig. 9 we compare this CFT formula, with $n=1$ (dashed curve) and the constant $C$ fixed by adapting to this case the method explained above, to the contour function for the entanglement entropy constructed in §3.2 for a single interval made by $\ell$ sites at the beginning of the massless harmonic chain in a segment of finite length $L$ with Dirichlet boundary conditions imposed at its endpoints. The correlators employed for this numerical analysis are given by (124) and (125) with $\omega=0$. Because of the Dirichlet boundary conditions, the contour obtained from the lattice data pass through the origin, as highlighted in the inset of the panel on the left.

In this case, where the massless regime can be considered without approximation, the agreement between the lattice data and the CFT curve is quite remarkable nearby the second endpoint of the interval and it gets worse close to the boundary. This is expected from the fact that the universal part in the CFT expression (74) does not contain information about the specific boundary conditions imposed at the endpoints of the segment $(0,L)$. This non universal information is encoded in the constant $C$.

4.3 Two disjoint intervals

A very interesting configuration to study involves a subsystem $A=A_{1}\cup A_{2}$ made by the union of two disjoint intervals $A_{1}$ and $A_{2}$. In an infinite system, we denote by $\ell_{1}$ and $\ell_{2}$ the lengths of $A_{1}$ and $A_{2}$ respectively, while $d$ is the distance separating the intervals (see the bottom panel of Fig. 4).

Considering the harmonic chain in the thermodynamic regime and for a very small mass ($\omega=10^{-8}$), in Fig. 10 we show the mode participation function (53) for two disjoint and equal intervals with $\ell_{1}=\ell_{2}=40$ separated by a distance $d$. The panels correspond to increasing values of $d$, starting from $d=0$, when $A_{1}\cup A_{2}$ is a single interval of length $\ell_{1}+\ell_{2}=80$, until a large value of $d$ (in the bottom right panel $d=100$), when the intervals are very far apart. For small $d$ the profile of the mode participation function restricted to one of the two intervals is clearly influenced by the presence of the other interval, while for large distances it becomes qualitatively like the mode participation function of a single interval (see the top left panel of Fig. 2 or the top left panel of Fig. 10).

Given the lattice setup described for Fig. 10, in Fig. 11 we show the contour function for the entanglement entropy of two equal and disjoint intervals ($\ell_{1}=\ell_{2}\equiv\ell$) for various lengths $\ell$ and for two fixed values of the dimensionless ratio $d/\ell$. All the data corresponding to the same value of $d/\ell$ nicely collapse on the same curve.

It would be very interesting to find an analytic function through a CFT analysis such that its integral over $A=A_{1}\cup A_{2}$ provides the entanglement entropy of two disjoint intervals for the massless scalar in two spacetime dimensions. This function, which is not known in the literature, would be a natural candidate to compare against the numerical results for the contour function for the entanglement entropy shown in Fig. 11 (left panel).

Despite the lack of a candidate function derived from CFT methods, we find it worth employing a function which captures the expected divergencies close to the endpoints of the two intervals. This function provides only part of the expected CFT result for the entanglement entropies of two disjoint intervals.

Inspired by the results of [9, 18], let us consider

[TABLE]

where $x\in(u_{1},v_{1})\cup(u_{2},v_{2})$ and we have introduced the following notation

[TABLE]

We stress that the expression obtained by plugging $f^{\prime}(x)$ from (75) into (65) does not provide the expected result for the massless free boson. Indeed, by removing the small discs of radius $\epsilon$ around the four endpoints of the intervals and integrating $f^{\prime}(x)$ from (75) over the remaining domain $x\in A_{\epsilon}=(u_{1}+\epsilon,v_{1}-\epsilon)\cup(u_{2}+\epsilon,v_{2}-\epsilon)$, one obtains

[TABLE]

where

[TABLE]

This is not the results found in [15, 16] for the entanglement entropies of two disjoint interval for this model because an important additive term is missing in (77). In particular, by adopting the notation of [16], the missing term is given by $\tfrac{1}{1-n}\log\mathcal{F}_{n}$, where $\mathcal{F}_{n}$ is a function of the harmonic ratio of the four endpoints which has been computed analytically for $n\geqslant 2$ in [15, 16]. In [16], an explicit expression written as an integral of an analytic function has been found for the $n=1$ case.

In the left panel of Fig. 11, the dashed curves correspond to the r.h.s. of (65) with the function $f(x)$ given by (75). The vertical shift due to the constant $C$ in (65) has been fixed by trying to maximise the overlap with the lattice data. We cannot employ the method adopted in the previous single interval cases because for two disjoint intervals a non trivial function obtained from CFT providing the entanglement entropies once integrated over $A_{1}\cup A_{2}$ is not available for the free massless scalar. Despite the fact that (75) is not expected to provide the correct result, the agreement between the dashed curves and the corresponding lattice data is quite satisfactory as shown in the left panel of Fig. 11.

In order to highlight the missing contribution due to the term $\tfrac{1}{1-n}\log\mathcal{F}_{n}$ in the Rényi entropies, let us introduce the contour function corresponding to the combination of entanglement entropies in (4), which gives the contour function for the mutual information when $n\to 1$. It reads

[TABLE]

where the contour functions are constructed as explained in §3.2. The functions $s^{(n)}_{A_{j}}(i)$ and $\theta_{A_{j}}(i)$ in (79) are respectively the contour function and the Heaviside step function corresponding to the interval $A_{j}$.

In the continuum, let us consider the following CFT expression

[TABLE]

being $c$ is the central charge ($c=1$ for the massless free boson) and

[TABLE]

where the function $f(x)$ is given by (75), the function $f_{A_{j}}(x)$ by (76) and $\theta_{A_{j}}(x)$ is the Heaviside step function with support in $A_{j}$. Because of the invariance under translations on the infinite line, we can set $u_{1}=0$, $v_{1}=\ell_{1}$, $u_{2}=\ell_{1}+d$ and $v_{2}=\ell_{1}+d+\ell_{2}$ in (82). This leads to the following expression

[TABLE]

where $\ell$ is a generic length. For instance, we can set either $\ell\equiv\ell_{1}+d+\ell_{2}$ or $\ell\equiv\ell_{1}+\ell_{2}$.

In the right panel of Fig. 11, by considering the harmonic chain with mass $\omega=10^{-8}$ in the thermodynamic regime, we show the contour function for the mutual information, which is given by (79) with $n=1$. The dashed curves are obtained from (80) and (83), where the constant $C$ has been fixed by trying to maximise the overlap with the corresponding curves found from the lattice data, as done in the left panel for the reason discussed above. The agreement is reasonable for $d/\ell=1/16$ while it gets worse for $d/\ell=1/2$. This is expected from the CFT expression for $\mathcal{F}_{n}$ [16], a function of the harmonic ratio of the four endpoints, which is $(1+d/\ell)^{-2}$ for two equal intervals of length $\ell$. Indeed, $\mathcal{F}_{n}\to 1$ when $d/\ell\to 0$ and $d/\ell\to\infty$, while it reaches its maximum for $(1+d/\ell)^{-2}=1/2$, i.e. when $d/\ell=\sqrt{2}-1$, which is close to $d/\ell=1/2$.

We find it worth remarking that a CFT candidate for the contour function for the entanglement entropies when $A$ is made by disjoint intervals is available in the case of free fermions on the infinite line [18]. In particular, for two disjoint intervals it is given by the function employed to plot the dashed curves in the left panel of Fig. 11. It would be interesting to compare this function with the lattice data coming from the contour function constructed in [26] for free fermions.

5 Alternative proposals

In this section we discuss two alternative constructions for the contour function which are different from the one presented in §3.2. The first one is based again on the symplectic matrix $W$ entering in (11) and it has been inspired by the proof of the Williamson’s theorem found in [34]. The second one is based on the mode participation function proposed by Botero and Reznik [4] and it has been studied more recently in [27]. In §5.3 we describe a deformation of the procedure explained in §3.1 which also provides a positive mode participation function satisfying the sum rule (15).

5.1 A proposal based on a proof of the Williamson’s theorem

In §3 we have shown that, given a real orthogonal matrix (23), the expression (29) provides a natural candidate for the contour function of the entanglement entropies and the corresponding mode participation function is (30). We find it worth focusing our attention on the orthogonal matrices which naturally occur in the analysis of the symplectic spectrum of $\gamma_{A}$. In §3 we discussed the orthogonal and symplectic matrix $K$ in (33) coming from the Euler decomposition of the symplectic matrix $W$ of the Williamson’s theorem (11). As already remarked in §3, another orthogonal matrix naturally related to the symplectic matrix $W$ involved in the Williamson’s theorem (11) is the orthogonal matrix $\widetilde{O}$ defined by (38), which enters in the factorisation (41) of $W$ provided in the constructive proof of the Williamson’s theorem found in [34].

The expression (40) tells us that $\widetilde{O}$ is the orthogonal matrix diagonalising $|\hat{\gamma}_{A}|$. Moreover, by specialising (18) and (29) to this case where $\Phi_{A}=|\hat{\gamma}_{A}|$, one finds respectively that the entanglement entropies can be written as

[TABLE]

and that the corresponding contour function reads

[TABLE]

Decomposing $\widetilde{O}$ in blocks like in (23), i.e.

[TABLE]

we have that the mode participation function in this case is (30) with $O=\widetilde{O}$, namely

[TABLE]

We find it instructive to discuss in some detail the simpler examples of reduced covariance matrices which are block diagonal, namely $\gamma_{A}=Q\oplus P$. In these cases, (41) tells us that $\widetilde{O}=U_{\widetilde{O}}\oplus V_{\widetilde{O}}$ and $W=U\oplus V$ are also block diagonal.

Being $\widetilde{O}$ orthogonal, both the matrices $U_{\widetilde{O}}$ and $V_{\widetilde{O}}$ on its diagonal are orthogonal. By plugging (48) and the block diagonal form of $\widetilde{O}$ into the square of (40), one finds that

[TABLE]

which tell us that $U_{\widetilde{O}}$ and $V_{\widetilde{O}}$ are the orthogonal matrices diagonalising the symmetric matrices $Q^{1/2}\,P\,Q^{1/2}$ and $P^{1/2}\,Q\,P^{1/2}$ respectively. Once $U_{\widetilde{O}}$ and $V_{\widetilde{O}}$ have been computed, the mode participation function in these cases is given by (87) with $Y_{\widetilde{O}}=Z_{\widetilde{O}}=\boldsymbol{0}$, namely

[TABLE]

The contour function $\tilde{s}^{(n)}_{A}(i)$ for the entanglement entropies can be constructed like in (16), with the mode participation function $\tilde{p}_{k}(i)$ instead of $p_{k}(i)$.

We have repeated the numerical analysis performed in §3 and §4 by employing the mode participation function (89). All the examples presented in §4.2 and §4.3 have been considered and basically the same curves have been found, except for the harmonic chain in thermodynamic limit at finite temperature $T=1/\beta$ and for the periodic chain with finite length $L$. The contour function for the entanglement entropy evaluated through the mode participation function (89) in these two cases is shown in Fig. 12 e Fig. 13 respectively. From the right panels of these figures one can clearly observe the disagreement with the corresponding CFT formulas, which provide straight lines whose slope is equal to one.

5.2 A proposal based on the eigenvectors of $(\textrm{i}J\gamma_{A})^{2}$

The first proposal we found in the literature of a mode participation function in harmonic lattices was suggested by Botero and Reznik [4]. The contour function associated to this mode participation function was studied in [27] and in this subsection we further discuss this construction. For simplicity, we focus only on block diagonal covariance matrices $\gamma_{A}$.

When $\gamma_{A}=Q\oplus P$, we have that $(\textrm{i}J\gamma_{A})^{2}=(PQ)\oplus(QP)$ and the symplectic spectrum is given by the positive square root of the spectrum of $QP$, as already mentioned in §3. In [4] the right eigenvectors of $QP$ and $PQ$ have been employed to construct a mode participation function. In the following we show that these eigenvectors can be read from the rows of the symplectic matrix $W$ entering in the Williamson’s theorem (11).

Given a block diagonal covariance matrix $\gamma_{A}=Q\oplus P$, the symplectic matrix $W=U\oplus V$ occurring in (11) is block diagonal as well and the symplectic condition for $W$ is equivalent to the relations $UV^{\textrm{t}}=\boldsymbol{1}=VU^{\textrm{t}}$, which tell us that $U^{-1}=V^{\textrm{t}}$ and $V^{-1}=U^{\textrm{t}}$. In particular, $U$ and $V$ are not orthogonal matrices. Then, specialising (41) to this case, one gets

[TABLE]

where the block diagonal structure of $\widetilde{O}=U_{\widetilde{O}}\oplus V_{\widetilde{O}}$ has been employed. From (48) and the block diagonal structure of $\widetilde{O}$, one finds that the relation (38) specified to this simpler case becomes $U_{\widetilde{O}}\,Q^{1/2}\,P^{1/2}\,V_{\widetilde{O}}^{\textrm{t}}=D$. By using this expression, one can check that $UV^{\textrm{t}}=\boldsymbol{1}$ holds for (90), as expected.

By employing (90) and the fact that $U_{\widetilde{O}}$ and $V_{\widetilde{O}}$ are orthogonal, we can write that

[TABLE]

which do not provide the diagonalization of the real and symmetric matrices $Q$ and $P$ because $U$ and $V$ are invertible but not orthogonal. Expressions relating $Q$ and $P$ to the orthogonal matrices $U_{\widetilde{O}}$ and $V_{\widetilde{O}}$ respectively are obtained by inverting the relations (90). They read $Q^{1/2}=U_{\widetilde{O}}^{\textrm{t}}\,D^{1/2}\,U=U^{\textrm{t}}D^{1/2}\,U_{\widetilde{O}}$ and $P^{1/2}=V_{\widetilde{O}}^{\textrm{t}}\,D^{1/2}\,V=V^{\textrm{t}}D^{1/2}\,V_{\widetilde{O}}$. From these results, (91) and the condition $UV^{\textrm{t}}=\boldsymbol{1}=VU^{\textrm{t}}$, one can check that the expressions in (88) are recovered, as expected.

Finally, by using (91) and $UV^{\textrm{t}}=\boldsymbol{1}=VU^{\textrm{t}}$, one finds

[TABLE]

which tell us that the invertible matrices $U$ and $V$ diagonalise $PQ$ and $QP$ respectively.

The first relation in (92) can be written as $QP\,V^{-1}=V^{-1}D^{2}$, which means that the $k$-th column of $V^{-1}=U^{\textrm{t}}$ is the right eigenvector $\boldsymbol{v}_{k}$ of $QP$ corresponding to the eigenvalue $\sigma_{k}^{2}$, namely $QP\,\boldsymbol{v}_{k}=\sigma_{k}^{2}\,\boldsymbol{v}_{k}$. In the same way, from the second relation in (92) one concludes that the $k$-th column of $U^{-1}=V^{\textrm{t}}$ provides the right eigenvector $\boldsymbol{u}_{k}$ of $PQ$ corresponding to the eigenvalue $\sigma_{k}^{2}$, i.e. $PQ\,\boldsymbol{u}_{k}=\sigma_{k}^{2}\,\boldsymbol{u}_{k}$. Denoting by $v_{k}(i)$ and $u_{k}(i)$ the $i$-th element of $\boldsymbol{v}_{k}$ and $\boldsymbol{u}_{k}$ respectively, in [4] the following expression has been proposed for the mode participation function

[TABLE]

By inverting $UV^{\textrm{t}}=\boldsymbol{1}$, one finds $(V^{-1})^{\textrm{t}}\,U^{-1}=\boldsymbol{1}$, and the diagonal elements of this relation provide the normalization condition $\sum_{k=1}^{\ell}\check{p}_{k}(i)=1$ with $1\leqslant k\leqslant\ell$. Given (93), the corresponding contour function $\check{s}^{(n)}_{A}(i)$ for the entanglement entropies can be constructed like in (16), with the mode participation function $\check{p}_{k}(i)$ instead of $p_{k}(i)$.

For subsystems made by two disjoint intervals, one can also define $\check{\Upsilon}_{A}^{(n)}$ like in (79), by replacing $s^{(n)}_{A}(i)$ with $\check{s}^{(n)}_{A}(i)$.

We remark that the property $\check{p}_{k}(i)\geqslant 0$ is not obvious for the proposal (93), as also emphasised in [27]. Actually, for the configuration where $A$ is made by two disjoint intervals, we have found several cases where $\check{p}_{k}(i)$ becomes negative. For instance, given two disjoint intervals with the same length $\ell_{1}=\ell_{2}=30$ separated by $d=10$ in the harmonic chain in the thermodynamic limit with $\omega=10^{-8}$, we found that the mode participation function (93) for $k=14$ reaches small negative values of order $10^{-4}$ for some sites, while the largest positive ones are of order $10^{-1}$. Nonetheless, the contour function obtained from this $\check{p}_{k}(i)$ is positive. Moreover, we have always found $\check{s}^{(n)}_{A}(i)\geqslant 0$ for all the sites in all the examples that we have considered.

In Fig. 14 we show the contour function constructed through (93) when $A$ is a single interval of length $\ell$ at the beginning of a massless open chain with Dirichlet boundary conditions in the thermodynamic limit. The dashed curve comes from the CFT expression (69) for $n=1$ and it is the same dashed curve showed in the top panels of Fig. 6. Thus, the contour function corresponding to the mode participation function (93) does not seem to provide the CFT expressions discussed in §4.2 in the scaling limit.

In Fig. 15 we have considered the contour function for the mutual information (namely Eq. (79) for $n=1$) for the periodic chain in the thermodynamic regime when $\omega=10^{-8}$, with the subsystem $A$ made by two equal and adjacent intervals of length $\ell$. The results obtained from the mode participation functions (53) (left panel) and (93) (right panel) are shown, together with the CFT expression (80) for $n=1$ (dashed curves). Comparing the insets of the two panels, we notice that the data collapse is better for the contour function constructed in §3.2 and we also observe that $\check{\Upsilon}_{A}^{(n=1)}$ becomes negative close to the two endpoints of the single interval $A_{1}\cup A_{2}$ for large values of $\ell$.

5.3 A deformation of the contour

We find it instructive to discuss briefly a deformation of construction described in §3 through a parameter $\xi\in(0,1)$. The resulting mode participation function provides a contour function which satisfies (6) and (7).

Given two real and positive parameters $\xi$ and $\eta$, let us revisit the construction described in §3.1 by considering the operator \breve{X}^{(i)}\,\equiv\,2\big{[}\xi\,\delta^{(i)}\oplus\eta\,\delta^{(i)}\big{]} associated to the $i$-th site instead of the projector $X^{(i)}$, which is recovered in the special case of $\xi=\eta=1/2$. The contour function in (29) with $X^{(i)}$ replaced by $\breve{X}^{(i)}$ gives

[TABLE]

where $\Phi_{A}$ is (17) with the orthogonal matrix $O$ written in terms of its blocks like in (23). By adapting the steps in (21) to this case, for (94) we find

[TABLE]

where in the last expression we have introduced

[TABLE]

which reduces to the mode participation function (30) for $\xi=\eta=1/2$, as expected.

Being $\xi$ and $\eta$ positive, we have that $\breve{p}_{k}(i)\geqslant 0$ by construction and this implies that $\breve{s}^{(n)}_{A}(i)\geqslant 0$ for the function in (94). Moreover, because of the relations (28) coming from the orthogonality condition $O^{\textrm{t}}O=\boldsymbol{1}$, we find that the property $\sum_{k=1}^{\ell}\breve{p}_{k}(i)=1$ is satisfied only when $\xi+\eta=1$. Thus, the contour function (94) fulfils the constraints (6) and (7) for \breve{X}^{(i)}\,\equiv\,2\big{[}\xi\delta^{(i)}\oplus(1-\xi)\delta^{(i)}\big{]} with $\xi\in(0,1)$.

Whenever $\xi\neq 1/2$, we have that $\sum_{i\in A}\widetilde{X}^{(i)}=2[\xi\,\boldsymbol{1}\oplus(1-\xi)\,\boldsymbol{1}]\neq\boldsymbol{1}$ and that $\breve{X}^{(i)}$ is not a projector because \big{(}\breve{X}^{(i)}\big{)}^{2}\neq\breve{X}^{(i)}. Since in §3.3 the properties of $X^{(i)}$ have been employed in a crucial way, this proposal for the contour function could be ruled out by the requirements introduced in [26]. Notice that, by employing the orthogonality condition $OO^{\textrm{t}}=\boldsymbol{1}$, one finds that the constraint $\sum_{i=1}^{\ell}\breve{p}_{k}(i)=1$ for any integer $k\in[1,\ell]$, which is not required for the mode participation function, can be fulfilled only for $\xi=1/2$.

6 Conclusions

In this paper we studied the contour function for the entanglement entropies in generic harmonic lattices. Our proposal is given by the expression in (45), which can be written also in the form (16) in terms of the particular mode participation function in (47). This proposal fulfils the basic constraints (6) and (7) and also three further requirements corresponding to a weakened version of the properties introduced in [26].

Focusing on one-dimensional harmonic chains with either periodic or Dirichlet boundary conditions in the massless regime, we have studied the configurations whose entanglement hamiltonian found though CFT methods can be written in the form (62), namely as an integral over the interval $A$ of the component $T_{00}$ of the energy-momentum tensor multiplied by a suitable local weight function [22, 23, 24]. Comparing the scaling limit of our contour function for the entanglement entropies with the inverse of the local weight function multiplying $T_{00}$ in the entanglement hamiltonian for various configurations, a good agreement has been observed in the regions of the interval nearby the endpoints, where a universal behaviour is expected. It remains to be understood more quantitatively how to extract the non-universal information contained in the contour function.

We have also considered the contour function for the entanglement entropy for a massless harmonic chain in the thermodynamic regime in some cases where $A$ is made by two disjoint intervals, finding well defined curves from the collapse of the numerical data in the scaling limit. It would be very interesting to find an analytic expression for these curves through a CFT analysis.

Other contour functions for the entanglement entropies satisfying the properties (6) and (7) have been considered (see §5) and we observed that these contour functions do not always agree with the functions coming from the corresponding entanglement hamiltonians obtained through CFT methods. We leave for future work to check whether they fulfil also the three requirements introduced in [26] or at least their weaker version discussed in §3.3.

The contour functions for the entanglement entropies provide a way to characterise the spatial structure of entanglement within the spatial region $A$. Few analysis of the contour functions for bipartite entanglement have been performed in the literature [4, 26, 27] and many issues could be explored in the future.

Finding a complete list of properties which allows to define the contour function for the entanglement entropies in a unique way is a very interesting open problem.

Contour functions are naturally defined also in higher dimensions. Since they give a more refined information on the entanglement structure than entanglement entropy, it could be interesting to study them in higher dimensional systems. As already noticed in [26], the contour function for fermions in two spatial dimensions can discriminate between the ground state of a massive theory and the one of a critical system with a finite number of zero modes. On the contrary, the entanglement entropy scales in both cases with the length of the boundary of the region $A$ and therefore it cannot distinguish the two cases. Moreover, in higher dimensions the shape of the spatial region $A$ plays a crucial role. For example, for critical models in two spatial dimensions, when $A$ has corners logarithmic terms are present which contain universal information [40]. It would be interesting to explore how the contour function encodes this and maybe further information.

We also find it worth studying further the contour function for the free fermions proposed by Chen and Vidal [26]: in the massless case one could compare the lattice results with the CFT expressions found in [24] in the case of $c=1/2$, as done in §4.2 for massless harmonic chain in the case of $c=1$. For the single interval on the infinite line this has been already done in [26] and it would be interesting to explore other configurations. We remark that, for free massless fermions and for subsystems $A$ made by a generic number of disjoint intervals, a natural CFT candidate for scaling limit of the contour function for the entanglement entropies is available in the literature [18].

Finding contour functions which fulfil all the proper requirements also for more interesting lattice models (e.g. the Ising spin chain) is an important open problem for future studies. The construction of contours for the entanglement entropies could be helpful to improve the current understanding about the role played by the entanglement spectrum [41] and by the eigenvectors of the reduced density [42, 6] in the characterisation of the bipartite entanglement. The interesting question to understand is whether and how is possible to find a general way to obtain a contour function directly from the modular hamiltonian. Some results for specific configurations and particular free models have been found in [26] and in this paper, but more complicated configurations and more interesting models (the interacting ones in particular) should be addressed. Recently, the modular hamiltonian for free lattice models has been employed to suggest another way to recover the universal function in (67) for $n=1$ through a lattice calculation [43]. Moreover, also analytical results for the modular hamiltonian of an interval on the infinite line for the free fermion have been found [44].

Tensor networks are very successful in describing the low energy physics of lattice models [45]. The scaling of entanglement in the ground state of lattice systems has been used as input in constructing the tensor networks algorithms which allow to simulate those systems [46]. Studying the spatial distribution of entanglement could potentially improve such tensor networks algorithms.

The problem of providing a list of properties which characterises the contour function in a unique way can be extended also to other measures of the bipartite entanglement. For entanglement measures different from the entanglement entropies explicit constructions of the corresponding contour functions are missing. A measure of the bipartite entanglement for mixed states which has attracted some attention is the logarithmic negativity [47]. A QFT approach to this measure of entanglement has been proposed in [48], with particular focus on CFTs (see also [49] for the massive case). Numerical calculations of the logarithmic negativity have been performed in lattice models, both in one spatial dimension [50] and in two-dimensional lattices [51]. In a forthcoming publication [52] we apply the approach described in this manuscript to study a contour function for the logarithmic negativity in generic harmonic lattices.

Acknowledgements

We are grateful to Marcus Cramer, Mihail Mintchev, Guifre Vidal and in particular to Pasquale Calabrese, Robert Myers and Luca Tagliacozzo for helpful discussions or correspondence. AC thanks SISSA, Trieste, and the University of Strathclyde, Glasgow, for warm hospitality and support during part of this work. ET acknowledges YITP, Kyoto, and NORDITA, Stockholm, for warm hospitality and support during part of this work. AC acknowledges financial support from MINECO (grant MTM2014-54240-P) and Comunidad de Madrid (grant QUITEMAD+-CM, ref. S2013/ICE-2801). AC has been also supported by European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 648913). ET has been supported by the ERC under Starting Grant 279391 EDEQS.

Appendices

Appendix A On the Williamson’s theorem

In this appendix we first review part of the proof of the Williamson’s theorem found in [34], which has been employed to construct the mode participation functions in §3 and in §5.1. Then, some issues related to the uniqueness of the symplectic matrix occurring in the Williamson’s theorem are briefly discussed.

Williamson’s theorem [32]. Given a $2m\times 2m$ real matrix $M$ which is also symmetric and positive definite, a real symplectic matrix $W$ exists such that

[TABLE]

where $D=\textrm{diag}(\sigma_{1},\dots,\sigma_{m})$ with $\sigma_{j}>0$. The sequence $\{\sigma_{1},\dots,\sigma_{m}\}$ is called symplectic spectrum of $M$. The symplectic spectrum can be obtained by taking the modulus of the spectrum of the matrix $JM$.

Various proofs of the Williamson’s theorem have been found [33] and we focus on the one given in [34]. In this proof the matrix $\hat{M}\equiv M^{1/2}JM^{1/2}$ is introduced. It is straightforward to realise that $\hat{M}$ and $JM$ have the same spectrum. Since $\hat{M}$ is real and antisymmetric, an orthogonal matrix $O$ exists such that

[TABLE]

being $D=\textrm{diag}(\sigma_{1},\dots,\sigma_{m})$ with $\sigma_{j}>0$.

By employing the orthogonal matrix $O$ defined in (98), the symplectic matrix $W$ entering in (97) is constructed as follows

[TABLE]

It is straightforward to check that the orthogonality condition for $O$ implies that the matrix $W$ in (99) satisfies (97). Moreover, by employing (98) one easily finds that $W$ in (99) is symplectic. Notice that the non uniqueness of the matrix $O$ is determined by the diagonalization problem defined by (98) and therefore it is related to the degeneracy of the symplectic eigenvalues.

In this manuscript the Williamson’s theorem is applied to the reduced covariance matrix $\gamma_{A}$ characterising the subsystem $A$. The symplectic matrix $W$ corresponding to $\gamma_{A}$ is crucial for the constructions of the contour functions for the entanglement entropies discussed in the main text. Indeed, the Euler decomposition of the symplectic matrix $W$ corresponding to $\gamma_{A}$ is employed to define the mode participation function (47) described in §3.2. Moreover, the orthogonal matrix defined through the relation (98) in the special case of $M=\gamma_{A}$ (see Eq. (38)) is the key object occurring in the construction of the mode participation function (89) discussed in §5.1.

Given a matrix $M$ satisfying the hypothesis of the Williamson’s theorem and assuming that $W$ and $W^{\prime}$ are two symplectic matrices such that (97) holds for both of them with the same $D$, it can proved that, $W^{\prime}\,W^{-1}\equiv\mathcal{U}$ is symplectic and orthogonal (Proposition 8.12 of [28]).

On one side, by considering the decomposition (99) for $W$ and the analogue one for $W^{\prime}$, which involves the same $D$ and $M$ but a different orthogonal matrix $O^{\prime}$, and plugging these decompositions into the relation $W^{\prime}\,W^{-1}\equiv\mathcal{U}$, it is straightforward to observe that O^{\prime}=\big{[}(D^{1/2}\oplus D^{1/2})\,\mathcal{U}\,(D^{-1/2}\oplus D^{-1/2})\big{]}O. On the other side, both $O$ and $O^{\prime}$ fulfil (98), and this leads to observe that $O|\hat{M}|O^{\textrm{t}}=D\oplus D=O^{\prime}|\hat{M}|(O^{\prime})^{\textrm{t}}$. Thus, the orthogonal matrices $O$ and $O^{\prime}$ are related through an orthogonal matrix $\mathcal{V}$ which does not mix eigenspaces corresponding to different eigenvalues in $D\oplus D$. Plugging $O^{\prime}=\mathcal{V}O$ into the relation between $O$ and $O^{\prime}$ written right above, one finds $\mathcal{V}=(D^{1/2}\oplus D^{1/2})\,\mathcal{U}\,(D^{-1/2}\oplus D^{-1/2})$. By first isolating $\mathcal{U}$ in the r.h.s. of the latter equation and then employing that $\mathcal{V}$ commutes with $D^{1/2}\oplus D^{1/2}$, we obtain that $\mathcal{V}=\mathcal{U}$.

The above observations support the claim that, ultimately, the non uniqueness of the orthogonal matrix $O$ depends on the degeneracies within the symplectic spectrum. Nonetheless, this freedom does not influence the contour function. Indeed, in the diagonalization problem defined by $O|\hat{M}|O^{\textrm{t}}=D\oplus D$ there is an obvious degeneracy due to the structure of the diagonal matrix in the right hand side. This allows to mix the operators $\hat{q}_{i}$ and $\hat{p}_{i}$ corresponding to the same site, but such freedom does not change the contour function (see the property (b) in §3.3 in the special case where $G$ is just the $i$-th site). When degeneracies occur in the symplectic spectrum, the contour (16) can be written as $\sum_{\kappa=1}^{\ell^{\prime}<\ell}\mathcal{P}_{\kappa}(i)s_{n}(\sigma_{\kappa})$, where the sum is performed over eigenspaces corresponding to different symplectic eigenvalues and $\mathcal{P}_{\kappa}(i)=\sum_{\kappa^{\prime}}p_{\kappa^{\prime}}(i)$ is the sum of the mode participation function over a base of the eigenspace corresponding to the symplectic eigenvalue $\sigma_{\kappa}$. From (47) one observes that orthogonal transformations mixing eigenvectors within the same eigenspace indexed by $\kappa$ do not change $\mathcal{P}_{\kappa}(i)$.

Let us conclude with an observation about the uniqueness of the matrix $\hat{M}$ introduced in the proof of the Williamson’s theorem found in [34]. In particular, let us try to extend the analysis by considering $\hat{M}_{(a)}\equiv M^{a}JM^{a}$ for a real power $a$, which becomes $\hat{M}$ in the special case of $a=1/2$. The matrices $\hat{M}_{(a)}$ and $JM^{2a}$ have the same spectrum. Since the matrix $\hat{M}_{(a)}$ is real and antisymmetric, a real orthogonal matrix $O_{(a)}$ exists such that

[TABLE]

By employing the orthogonality of $O_{(a)}$, it is straightforward to show that the relation M=W_{(a)}^{\textrm{t}}\big{(}D_{(a)}\oplus D_{(a)}\big{)}W_{(a)} is satisfied for the matrix $W_{(a)}\equiv(D_{(a)}^{-1/2}\oplus D_{(a)}^{-1/2})\,O_{(a)}M^{1/2}$. Then, by requiring that $W_{(a)}$ satisfies the symplectic condition, one finds that $O_{(a)}\hat{M}\,O_{(a)}^{\textrm{t}}$ gives the r.h.s. of (100). This allows to conclude that $O_{(a)}\hat{M}_{(a)}O_{(a)}^{\textrm{t}}=O_{(a)}\hat{M}\,O_{(a)}^{\textrm{t}}$, which leads to $\hat{M}_{(a)}=\hat{M}$ and therefore $a=1/2$.

Appendix B On the properties of the contour function

In this appendix we consider the three properties of the contour function for the entanglement entropies introduced in [26] besides the constraints (6) and (7). Focussing on the harmonic lattices, in the following analysis, which supports and extends the discussion made in §3.3, we show that these three further properties are satisfied by the contour function described in §3.2 if we restrict to the canonical transformations implemented by the matrices $M\in\textrm{Sp}(\ell)\cap O(2\ell)$.

(a) Spatial symmetry. If $\rho_{A}$ is invariant under a transformation relating the sites $i$ and $j$ in the subsystem $A$, then $s_{A}^{(n)}(i)=s_{A}^{(n)}(j)$.

Considering a canonical transformation that maps the $i$-th site into the $j$-th site, for the corresponding symplectic $M$ we have that $X^{(i)}\to X^{(j)}=M^{\textrm{t}}X^{(i)}M$. Since $X^{(j)}$ is a projector (i.e. $(X^{(j)})^{2}=X^{(j)}$), the matrix $M$ is orthogonal. Assuming that this transformation is a symmetry of $\rho_{A}$, from the transformation rule (56) we conclude that E_{\textrm{\tiny R}}^{-1}\gamma_{A}E_{\textrm{\tiny R}}^{-1}=M\big{(}E_{\textrm{\tiny R}}^{-1}\gamma_{A}E_{\textrm{\tiny R}}^{-1}\big{)}M^{\textrm{t}} for this particular $M$. Thus, under this mapping we have

[TABLE]

This shows that the contour (44) fulfils the constraint (a) about the spatial symmetry of the subsystem.

In order to formulate the remaining properties, we need to introduce also the contour $s_{A}^{(n)}(G)$ of a subregion $G\subseteq A$ as follows

[TABLE]

In the special case of $G=A$, from (6) and (103) we find $s_{A}^{(n)}(A)=S_{A}^{(n)}$. The contour $s_{A}^{(n)}(G)$ is clearly additive: for any two non intersecting spatial subsets $G\subsetneq A$ and $\tilde{G}\subsetneq A$ we have $s_{A}^{(n)}(G\cup\tilde{G})=s_{A}^{(n)}(G)+s_{A}^{(n)}(\tilde{G})$. Moreover, the contour $s_{A}^{(n)}(G)$ is monotonous, i.e. for $G\subseteq\tilde{G}\subseteq A$ the inequality $s_{A}^{(n)}(G)\leqslant s_{A}^{(n)}(\tilde{G})$ holds.

(b) Invariance under local unitary transformations. Given a system in the state characterised by the density matrix $\rho$ and a unitary transformation $U_{G}$ acting non trivially only on $G\subseteq A$, denoting by $\rho^{\prime}$ the state of the system after such transformation, the same contour $s_{A}^{(n)}(G)$ should be found for $\rho$ and $\rho^{\prime}$.

In order to discuss the requirement (b), let us introduce the projector for the subregion $G\subseteq A$ as follows

[TABLE]

By employing the realisation of the matrices $X^{(i)}$ described in §3.1, we have that $X^{(G)}$ is the identity on $G$, while it vanishes outside. The unitary transformation $U_{G}$ corresponds to a symplectic matrix $M_{G}\in\textrm{Sp}(\ell)$ which acts non trivially only on $G$. We recall that in our analysis we restrict to symplectic matrices that are also orthogonal, i.e. $M_{G}\in\textrm{Sp}(\ell)\cap O(2\ell)$. Since $M_{G}$ acts non trivially only on $G$, where $X^{(G)}$ is the identity, we have that $[X^{(G)},M_{G}]=0$.

From (45) and (103) we get

[TABLE]

where we have used first that $\boldsymbol{1}=M_{G}^{\textrm{t}}\,M_{G}$, then the cyclic property of the trace and finally the fact that $X^{(G)}$ and $M_{G}$ commute. Thus, also the property (b) is satisfied for the restricted class of unitary transformations $U_{G}$ associated to $M_{G}\in\textrm{Sp}(\ell)\cap O(2\ell)$.

(c) A bound. Given a system in the pure state $|\Psi\rangle$ and the bipartition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, let us assume that the further decompositions $\mathcal{H}_{A}=\mathcal{H}_{\Omega_{A}}\otimes\mathcal{H}_{\bar{\Omega}_{A}}$ and $\mathcal{H}_{B}=\mathcal{H}_{\Omega_{B}}\otimes\mathcal{H}_{\bar{\Omega}_{B}}$ lead to the following factorisation of the state

[TABLE]

Considering a subregion $G\subseteq A$ such that $\bigotimes_{i\in G}\mathcal{H}_{i}\subseteq\mathcal{H}_{\Omega_{A}}$, we must have that

[TABLE]

where $S^{(n)}(\Omega_{A})$ are the entanglement entropies corresponding to the reduced density matrix $\rho_{\Omega_{A}}$, obtained by tracing over the degrees of freedom of $\mathcal{H}_{\bar{\Omega}_{A}}\otimes\mathcal{H}_{B}$.

We refer the interested reader to [26] for a more detailed discussion on the motivations leading to this property.

Assuming that $G$ is made by the first $\ell_{G}<\ell$ sites of $A$, let us order the modes by considering the vector $\hat{\boldsymbol{r}}_{G}$ given by $\{\hat{r}_{i,\alpha}\,|\,i\in G\}$ and $\hat{\boldsymbol{r}}_{\bar{G}}$ which collects the remaining ones $\{\hat{r}_{i,\alpha}\,|\,i\in A\setminus G\}$, being $\hat{r}_{i,1}\equiv\hat{q}_{i,2}$ and $\hat{r}_{i,2}\equiv\hat{p}_{i}$. Similarly, we can consider the modes $\hat{\boldsymbol{w}}_{\Omega_{A}}$ and $\hat{\boldsymbol{w}}_{\bar{\Omega}_{A}}$ given by $\{\hat{w}_{m,\alpha}\,|\,m\in\Omega_{A}\}$ and $\{\hat{w}_{m,\alpha}\,|\,m\in\bar{\Omega}_{A}\}$ respectively. Notice that $\ell_{G}\leqslant\ell_{\Omega_{A}}$, being $2\ell_{\Omega_{A}}$ the number of elements in $\hat{\boldsymbol{w}}_{\Omega_{A}}$. A consequence of the hypothesis $\bigotimes_{i\in G}\mathcal{H}_{i}\subseteq\mathcal{H}_{\Omega_{A}}$ is that the $\hat{\boldsymbol{r}}$ modes are related to the $\hat{\boldsymbol{w}}$ modes through a linear map such that all the modes $\hat{\boldsymbol{r}}_{G}$ are contained into $\hat{\boldsymbol{w}}_{\Omega_{A}}$, namely

[TABLE]

where $V_{\Omega G}$, $V_{\Omega\bar{G}}$ and $V_{\bar{\Omega}\bar{G}}$ are rectangular matrices with proper sizes which partition the $2\ell\times 2\ell$ matrix $V$.

In order to preserve the canonical commutation relations, the matrix $V$ in (110) must be symplectic, but we also impose that $V$ is orthogonal. As for the covariance matrices $\gamma_{\Omega}$ and $\gamma_{A}$ corresponding to the $\hat{\boldsymbol{r}}$ modes and to the $\hat{\boldsymbol{w}}$ modes respectively, from (108) we have that $\gamma_{\Omega}=\gamma_{\Omega_{A}}\oplus\gamma_{\bar{\Omega}_{A}}$. Combining this observation with (110), we obtain that

[TABLE]

From $s^{(n)}_{A}(G)$ in (105) and the polar decomposition (32) for the symplectic matrix occurring in the Williamson’s theorem for $\gamma_{\Omega}$, which provides the matrix E_{\textrm{\tiny R},{\textrm{\tiny\Omega}}}, we find that

[TABLE]

where in (112) we used (56), the step (113) has been obtained from E_{\textrm{\tiny R},{\textrm{\tiny\Omega}}}=E_{\textrm{\tiny R},{\textrm{\tiny\Omega_{A}}}}\oplus E_{\textrm{\tiny R},{\textrm{\tiny\bar{\Omega}{A}}}}, which is a consequence of $\gamma_{\Omega}=\gamma_{\Omega_{A}}\oplus\gamma_{\bar{\Omega}_{A}}$, and (114) comes from $VX^{(G)}V^{\textrm{t}}=V_{\Omega G}\,V_{\Omega G}^{\textrm{t}}\oplus\boldsymbol{0}$. The inequality in (115) is obtained by employing the fact that $V_{\Omega G}\,V_{\Omega G}^{\textrm{t}}$ is a projector. Equivalently, one can employ that from the condition $VV^{\textrm{t}}=\boldsymbol{1}$ one gets $V_{\Omega G}\,V_{\Omega G}^{\textrm{t}}=\boldsymbol{1}-V_{\Omega\bar{G}}\,V_{\Omega\bar{G}}^{\textrm{t}}$. Then, by first plugging the latter expression into (114) and then discarding the term with the positive definite matrix $V_{\Omega\bar{G}}\,V_{\Omega\bar{G}}^{\textrm{t}}$, we find the inequality in (115), which can be written by exploiting that s_{n}(E_{\textrm{\tiny R},{\textrm{\tiny\Omega{A}}}}^{-1}\gamma_{\Omega_{A}}E_{\textrm{\tiny R},{\textrm{\tiny\Omega_{A}}}}^{-1}) is positive definite.

Appendix C Correlators

In this appendix we collect the correlators employed in our numerical analysis. Considering the one-dimensional harmonic chain (8) with $L$ lattice sites, in the following we write the two-point correlators corresponding to either periodic (§C.1) or Dirichlet (§C.2) boundary conditions.

C.1 Periodic chain

The periodic harmonic chain is defined by the hamiltonian (8) with boundary conditions given by $q_{L}=q_{0}$ and $p_{L}=p_{0}$. By introducing the creation and annihilation operators in the standard way, for the vacuum state one finds the following two-point correlators

[TABLE]

where the dispersion relation reads

[TABLE]

It is important to remark that the correlator $\langle\hat{q}_{i}\hat{q}_{j}\rangle$ in (116) is not well defined for $\omega=0$. Indeed, for the zero mode $k=0$ we have $\omega_{0}|_{\omega=0}=0$ and therefore the corresponding term in $\langle\hat{q}_{i}\hat{q}_{j}\rangle$ diverges as $\omega\to 0$. The occurrence of the zero mode is due to the translation invariance of the model with periodic boundary conditions. In order to explore the massless regime, in our numerical analysis we keep $\omega$ non vanishing and such that $1/\omega$ is much larger than the other scales involved in the analysis.

We also consider the periodic harmonic chain in the thermal state with temperature $T=1/\beta$. For this state the correlators entering in the numerical analysis read

[TABLE]

Also in this case the correlator (118) diverges as $\omega\to 0$ because of the occurrence of zero mode at $k=0$.

In the thermodynamic limit $L\to\infty$, the vacuum state correlators (116) become [4]

[TABLE]

where $z\equiv(\omega-\sqrt{\omega^{2}+4\kappa/m}\,)^{2}/(4\kappa/m)$. As for the correlators (118) and (119) at finite temperature, in the thermodynamic limit they become

[TABLE]

with $\omega_{k}=\sqrt{\omega^{2}+(4\kappa/m)\sin^{2}(k/2)}$.

The above two-point functions provide the elements of the matrices $Q$ and $P$ entering in the reduced covariance matrix $\gamma_{A}$, whose symplectic spectrum leads to the entanglement entropies as discussed in §2.

C.2 Open chain with Dirichlet boundary conditions

The open harmonic chain with Dirichlet boundary conditions is defined by the hamiltonian (8) with the conditions $q_{0}=q_{L+1}=0$ and $p_{0}=p_{L+1}=0$ imposed at its endpoints. Because of these boundary conditions, the invariance under translations does not occur.

The vacuum state correlators for this model read [53]

[TABLE]

where the dispersion relation is given by

[TABLE]

We remark that the correlators (124) and (125) are well defined also in the massless regime, i.e. finite quantities are obtained when $\omega=0$. This key feature significantly improves the comparison between the numerical results from the lattice and the corresponding ones obtained through CFT methods.

In the thermodynamic limit $L\to\infty$, the two-point functions (124) and (125) become

[TABLE]

where $\omega_{k}$ has been defined in the text below (123). When $\omega=0$ these integrals can be performed analytically and the results read respectively [48]

[TABLE]

being $\psi(z)$ the digamma function.