Symmetry resolved entanglement in two-dimensional systems via dimensional reduction

Given a bipartition of a system in a pure state |ψ⟩ into A ∪ B, the reduced density matrix (RDM) ρ_A of the subsystem A is defined by tracing over the degrees of freedom of the subsystem B, i.e.

$$\begin{equation}{\rho }_{\mathrm{A}}={\mathrm{Tr}}_{\mathrm{B}}\enspace \rho ,\end{equation} \tag{ 1 }$$

where ρ = |ψ⟩⟨ψ| is the density matrix of the entire system. A measure of the bipartite entanglement is given by the Rényi entanglement entropies S_n, defined as

$$\begin{equation}{S}_{n}\equiv \frac{1}{1-n}\enspace \mathrm{log}\enspace \mathrm{Tr}\enspace {\rho }_{\mathrm{A}}^{n},\end{equation} \tag{ 2 }$$

whose limit n → 1 is the von Neumann entanglement entropy, i.e., S₁ = −Tr ρ_A log ρ_A. These entanglement entropies have been investigated for a large number of extended quantum systems in several different physical situations and with many different techniques (see e.g. references [4–6] as reviews).

Below we only report some results that we will use in the following sections. We start by considering the one-dimensional tight binding model, i.e., the free spinless fermions described by the Hamiltonian

$$\begin{equation}{H}_{\mathrm{F}\mathrm{F}}=-\frac{1}{2}\sum _{i}\left({c}_{i+1}^{{\dagger}}{c}_{i}+{c}_{i}^{{\dagger}}{c}_{i+1}\right)+\sum _{i}\mu {c}_{i}^{{\dagger}}{c}_{i},\end{equation} \tag{ 3 }$$

where μ is the chemical potential, c_i and ${c}_{i}^{{\dagger}}$ the ladder operators of the fermions obeying standard anticommutation relations $\left\{{c}_{i},{c}_{j}^{{\dagger}}\right\}={\delta }_{ij}$. We only focus on the ground state here. When |μ| < 1 the theory is gapless. The Jordan Wigner transformation maps the model to the spin-1/2 XX chain in a magnetic field. We consider the subsystem A to be an interval made of ℓ consecutive sites. For large ℓ, the asymptotic scaling of the entanglement entropies is given by [48, 49]

$$\begin{equation}{S}_{n}=\frac{1+{n}^{-1}}{6}\enspace \mathrm{log}\left(2\ell \enspace \mathrm{sin}\enspace {k}_{\mathrm{F}}\right)+{{\Upsilon}}_{n},\end{equation} \tag{ 4 }$$

where k_F = arccos(μ/2) is the Fermi momentum and

$$\begin{equation}{{\Upsilon}}_{n}=\frac{n+1}{n}{\int }_{0}^{\infty }\frac{\mathrm{d}t}{t}\left[\frac{1}{1-{n}^{-2}}\left(\frac{1}{n\enspace \mathrm{sinh}\enspace t/n}-\frac{1}{\mathrm{sinh}\enspace t}\right)\frac{1}{\mathrm{sinh}\enspace t}-\frac{{\mathrm{e}}^{-2t}}{6}\right].\end{equation} \tag{ 5 }$$

The leading logarithmic term in equation (4) is universal and follows from conformal field theory [7–10]; in contrast, the non-universal constant ϒ_n has been derived using the Fisher–Hartwig conjecture [48]. For a finite system of length L with PBC's, the same form also holds replacing ℓ with $\frac{L}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi \ell }{L}$ [9].

Exact results are available also for free bosonic systems on the lattice, i.e., for the harmonic chain. In this case, one exploits the Baxter corner transfer matrix (CTM) approach [50]. A chain of oscillators of mass M = 1 with frequency ω₀, coupled together by springs with elastic constant k (which, without loss of generality, we assume to be k = 1 − ω₀), is described by the Hamiltonian

$$\begin{equation}{H}_{\mathrm{B}}=\frac{1}{2}\sum _{i}{p}_{i}^{2}+{\omega }_{0}^{2}{q}_{i}^{2}+k{\left({q}_{i+1}-{q}_{i}\right)}^{2},\end{equation} \tag{ 6 }$$

where the p_i and q_i satisfy the canonical commutation relations [q_i, q_j] = [p_i, p_j] = 0 and [q_i, p_j] = iδ_ij. A canonical transformation of variables (p_i, q_i) allows us to have only one relevant system parameter in equation (6), i.e., ${\omega }_{0}^{2}/k$ or, equivalently, (1 − k)²/k [51].

In references [52–54] it has been shown that a harmonic chain is related to a two-dimensional classical Gaussian model. Such correspondence is at the basis of the CTM approach which allows us to write the RDM as the partition function of the two-dimensional classic model. For the case of an infinite subsystem size (namely, a semi-infinite line), this correspondence allows us to express ρ_A as (up to a prefactor) [54–56]

$$\begin{equation}{\rho }_{\mathrm{A}}\sim {\mathrm{e}}^{-{\mathcal{H}}_{\mathrm{C}\mathrm{T}\mathrm{M}}},\end{equation} \tag{ 7 }$$

where ${\mathcal{H}}_{\mathrm{C}\mathrm{T}\mathrm{M}}$ is an effective Hamiltonian which, due to the Gaussian nature of the model, can be diagonalised. This means that all the eigenvalues of the RDM can be determined exactly. In particular, in reference [54] it was shown that

$$\begin{equation}{\mathcal{H}}_{\mathrm{C}\mathrm{T}\mathrm{M}}=\sum _{j=0}^{\infty }{\epsilon}\left(2j+1\right){\beta }_{j}^{{\dagger}}{\beta }_{j},\quad {\epsilon}=\frac{\pi I\left(\sqrt{1-{\kappa }^{2}}\right)}{I\left(\kappa \right)},\end{equation} \tag{ 8 }$$

where I(κ) is the complete elliptic integral of the first kind, κ can be defined in terms of the system parameters as $\kappa /{\left(1-\kappa \right)}^{2}=k/{\omega }_{0}^{2}$ and ${\beta }_{j},{\beta }_{j}^{{\dagger}}$ are bosonic ladder operators. From equation (2), we easily read off the Rényi entropies of a harmonic chain as [9]

$$\begin{equation}{S}_{n}=\sum _{j=0}^{\infty }n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}}\right]-\sum _{j=0}^{\infty }\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}}\right].\end{equation} \tag{ 9 }$$

In the critical regime in which ω₀ → 0 (i.e., → 0), we recover the field theory result [9]

$$\begin{equation}{S}_{n}\simeq \frac{1+{n}^{-1}}{12}\enspace \mathrm{log}\enspace \xi .\end{equation} \tag{ 10 }$$

where $\xi \sim {\omega }_{0}^{-1}$ is the correlation length (the inverse gap) of the system. Although all the previous formulas are valid for a semi-infinite line, we have that for a finite subsystem of length ℓ, as long as ℓ ≫ ξ, the clustering of the RDM implies that it becomes the product of two independent ones at the two boundaries [57] and hence that the Rényi entropies are just the double of the one above for a single boundary, with exponentially suppressed corrections in ℓ/ξ [58, 59].

In the following, we will mainly be interested in systems with a U(1) internal symmetry. We then consider a complex bosonic theory which on the lattice is a chain of complex oscillators (we dub complex or double harmonic chain). The latter is the sum of two real harmonic chains in the variables (p⁽¹⁾, q⁽¹⁾) and (p⁽²⁾, q⁽²⁾), i.e.

$$\begin{equation}{H}_{\mathrm{C}\mathrm{B}}\left({p}^{\left(1\right)}+\mathrm{i}{p}^{\left(2\right)},{q}^{\left(1\right)}+\mathrm{i}{q}^{\left(2\right)}\right)={H}_{\mathrm{B}}\left({p}^{\left(1\right)},{q}^{\left(1\right)}\right)+{H}_{\mathrm{B}}\left({p}^{\left(2\right)},{q}^{\left(2\right)}\right).\end{equation} \tag{ 11 }$$

For the double chain the entanglement Hamiltonian is the sum of two ${\mathcal{H}}_{\mathrm{C}\mathrm{T}\mathrm{M}}$ of the form (8). For each of them we introduce the associated ladder operators, β_1,i and β_2,i. Since these operators commute, the RDM factorises as

$$\begin{equation}{\rho }_{\mathrm{A}}={\rho }_{\mathrm{A}}^{{\beta }_{1}}\bigotimes{\rho }_{\mathrm{A}}^{{\beta }_{2}},\end{equation} \tag{ 12 }$$

where we denoted the RDM for β_1,i and β_2,i with ${\rho }_{\mathrm{A}}^{{\beta }_{1}}$ and ${\rho }_{\mathrm{A}}^{{\beta }_{2}}$ respectively. Therefore the Rényi entropies for a complex chain are just the double of those for a real chain.

Let us focus now on systems with a U(1) symmetry and with an associated conserved charge denoted by Q. Moreover, we assume it is possible to write Q as the sum of local operators, i.e., Q = ∑_iQ_i, where Q_i is the local contribution. Taking the trace over B of [ρ, Q] = 0, and defining Q_A = ∑_i∈AQ_i, we find that [ρ_A, Q_A] = 0. This implies that ρ_A is block-diagonal and each block corresponds to a different charge sector labelled by the eigenvalue q of Q_A

$$\begin{equation}{\rho }_{\mathrm{A}}={\oplus }_{q}{{\Pi}}_{q}{\rho }_{\mathrm{A}},\end{equation} \tag{ 13 }$$

where Π_q is the projector onto the subspace of states of region A with charge q (with a slight abuse of notation we just wrote Π_qρ_A instead of Π_qρ_AΠ_q, being sure that will not generate any confusion). Such decomposition defines the (normalised) RDM in the q-sector, ρ_A(q) ≡ Π_qρ_A/Tr(Π_qρ_A). The generalised moments

$$\begin{equation}{\mathcal{Z}}_{n}\left(q\right)\equiv \mathrm{Tr}\left({{\Pi}}_{q}{\rho }_{\mathrm{A}}^{n}\right),\end{equation} \tag{ 14 }$$

are useful to write down the symmetry resolved Rényi and von Neumann entropies, meaning the entropies associated to each ρ_A(q) that are given by

$$\begin{equation}{S}_{n}\left(q\right)\equiv \frac{1}{1-n}\enspace \mathrm{log}\enspace \mathrm{Tr}{\left({\rho }_{\mathrm{A}}\left(q\right)\right)}^{n}=\frac{1}{1-n}\enspace \mathrm{log}\left[\frac{{\mathcal{Z}}_{n}\left(q\right)}{{\mathcal{Z}}_{1}^{n}\left(q\right)}\right],\quad {S}_{1}\left(q\right)=\underset{n\to 1}{\mathrm{lim}}\enspace {S}_{n}\left(q\right).\end{equation} \tag{ 15 }$$

Notice that ${\mathcal{Z}}_{1}\left(q\right)$ is the probability p(q) of finding q in a measurement of Q_A, i.e., p(q) = Tr(Π_qρ_A). The total von Neumann entanglement entropy S₁ and the symmetry resolved ones satisfy the relation [47, 60]

$$\begin{equation}{S}_{1}=\sum _{q}p\left(q\right){S}_{1}\left(q\right)-\sum _{q}p\left(q\right)\mathrm{log}\enspace p\left(q\right).\end{equation} \tag{ 16 }$$

equation (16) has a clear physical interpretation. The first contribution is known as configurational entanglement entropy (S^c) and depends on the entropy of each charge sector, weighted with its probability. The second contribution is the fluctuation entanglement entropy (S^f) which is due, as the name says, to the fluctuations of the charge within the subsystem. The fluctuation entropy S^f, also known as number entropy, has been studied in many different context [47, 61–65]. More complicated formulas can also be written for the Rényi entropies [66].

Finally, we can define the (normalised) charged moments of ρ_A as

$$\begin{equation}{Z}_{n}\left(\alpha \right)\equiv \mathrm{Tr}\enspace {\rho }_{\mathrm{A}}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{Q}_{\mathrm{A}}\alpha },\end{equation} \tag{ 17 }$$

which, importantly, are related to the generalised moments (14) by Fourier transform, i.e., [40]

$$\begin{equation}{\mathcal{Z}}_{n}\left(q\right)={\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }{Z}_{n}\left(\alpha \right).\end{equation} \tag{ 18 }$$

Recently, such relation allowed to derive interesting results about the different symmetry-resolved contributions for CFTs, free gapped and gapless systems of bosons and fermions, integrable spin chains and disordered systems. As for the total entropy, we will review only the results relevant for our purposes, while a plethora of others can be found in the literature [39–46, 66–78].

In the one-dimensional tight binding model, the generalised Fisher–Hartwig conjecture has been used to obtain the asymptotic behaviour of the symmetry resolved entropies at leading and subleading orders [45, 46]. The charged moments are

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}\left(\alpha \right)=\frac{\mathrm{i}{k}_{\mathrm{F}}\ell }{\pi }\alpha -\left[\frac{1}{6}\left(n-\frac{1}{n}\right)+\frac{2}{n}{\left(\frac{\alpha }{2\pi }\right)}^{2}\right]\mathrm{log}\enspace 2\ell \enspace \mathrm{sin}\enspace {k}_{\mathrm{F}}+{\Upsilon}\left(n,\alpha \right)+o\left(1\right),\end{equation} \tag{ 19 }$$

where ϒ(n, α) is a real and even function of α defined as

$$\begin{equation}{\Upsilon}\left(n,\alpha \right)=ni{\int }_{-\infty }^{\infty }\mathrm{d}w\left[\mathrm{tanh}\left(\pi w\right)-\mathrm{tanh}\left(\pi nw+\mathrm{i}\alpha /2\right)\right]\mathrm{log}\enspace \frac{{\Gamma}\left(\frac{1}{2}+\mathrm{i}w\right)}{{\Gamma}\left(\frac{1}{2}-\mathrm{i}w\right)}.\end{equation} \tag{ 20 }$$

Taking the Fourier transforms and expanding for large ℓ, one gets for the symmetry resolved entropies at the leading orders [45]

$$\begin{equation}{S}_{n}\left(q\right)={S}_{n}-\frac{1}{2}\enspace \mathrm{ln}\left(\frac{2}{\pi }\enspace \mathrm{ln}\left(2\ell \enspace \mathrm{sin}\enspace {k}_{\mathrm{F}}\right)\right)+\frac{\mathrm{ln}\enspace n}{2\left(1-n\right)}+o\left(1\right).\end{equation} \tag{ 21 }$$

The fact that up to order O(1) the symmetry resolved entropies do not depend on q has been dubbed equipartition of entanglement [43]. The first term breaking equipartition appears at order O(1/(log ℓ)²) [45].

The same quantities have been also investigated for off-critical quantum bosonic chains through the Baxter's CTM for the bipartition in two semi-infinite systems [44] (and generalised to finite subsystems in [67]). The approach of the previous subsection can, in fact, be adapted to the computation of the symmetry resolved entropies in the non-critical complex harmonic chain which (in contrast with its real analogue) possesses a U(1) symmetry. The starting point is to write the charge operator Q_A in terms of the β₁'s and β₂'s ladder operators associated to the bosons of the two (real) chains (see equation (12)), as [44]

$$\begin{equation}{Q}_{\mathrm{A}}=\sum _{j=0}^{\infty }{\beta }_{1,j}^{{\dagger}}{\beta }_{1,j}-{\beta }_{2,j}^{{\dagger}}{\beta }_{2,j}.\end{equation} \tag{ 22 }$$

For the charged moments, we need to compute $\mathrm{Tr}\left({\rho }_{\mathrm{A}}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{Q}_{\mathrm{A}}\alpha }\right)$, but using the form in equation (22) for Q_A, the trace factorises as

$$\begin{equation}{Z}_{n}\left(\alpha \right)=\mathrm{Tr}\enspace {\rho }_{\mathrm{A}}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{Q}_{\mathrm{A}}\alpha }=\mathrm{Tr}\left[{\left({\rho }_{\mathrm{A}}^{{\beta }_{1}}\right)}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{N}_{\mathrm{A}}^{{\beta }_{1}}\alpha }\right]{\times}\left[\mathrm{Tr}{\left({\rho }_{\mathrm{A}}^{{\beta }_{2}}\right)}^{n}\enspace {\mathrm{e}}^{-\mathrm{i}{N}_{\mathrm{A}}^{{\beta }_{2}}\alpha }\right],\end{equation} \tag{ 23 }$$

where ${N}_{\mathrm{A}}^{{\beta }_{1}}={\sum }_{j\in \mathrm{A}}{\beta }_{1,j}^{{\dagger}}{\beta }_{1,j}$ and ${N}_{\mathrm{A}}^{{\beta }_{2}}={\sum }_{j\in \mathrm{A}}{\beta }_{2,j}^{{\dagger}}{\beta }_{2,j}$. The two factors are equal, except for the sign of α. Therefore, the charged moments for the complex harmonic chain are [44]

$$\begin{align}\mathrm{log}\enspace {Z}_{n}\left(\alpha \right)=\sum _{j=0}^{\infty }2n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}}\right]-\sum _{j=0}^{\infty }\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}+\mathrm{i}\alpha }\right]\\ \quad \quad \quad \quad \quad \quad -\sum _{j=0}^{\infty }\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}-\mathrm{i}\alpha }\right].\hfill \end{align} \tag{ 24 }$$

In the critical region → 0, equation (24) reduces to

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}\left(\alpha \right)\simeq \left[\frac{1}{n}{\left(\frac{\alpha }{2\pi }\right)}^{2}-\frac{\vert \alpha \vert }{2\pi n}+\frac{1}{6n}-\frac{n}{6}\right]\mathrm{log}\enspace \xi +O\left(1\right).\end{equation} \tag{ 25 }$$

A Fourier transform allows to get the generalised moments and, eventually, the symmetry resolved entropies. We only report the final result, which reads

$$\begin{equation}{S}_{n}\left(q\right)=\frac{2}{1-n}\sum _{k=1}^{\infty }\left[n\enspace \mathrm{log}\left(1-{\mathrm{e}}^{-2{\epsilon}k}\right)-\mathrm{log}\left(1-{\mathrm{e}}^{-2n{\epsilon}k}\right)\right]+\frac{1}{1-n}\enspace \mathrm{log}\enspace \frac{{{\Phi}}_{q}\left({\mathrm{e}}^{-n{\epsilon}}\right)}{{\left({{\Phi}}_{q}\left({\mathrm{e}}^{-{\epsilon}}\right)\right)}^{n}},\end{equation} \tag{ 26 }$$

with

$$\begin{equation}{{\Phi}}_{q}\left(u\right)=\sum _{k=0}^{\infty }{\left(-1\right)}^{k}{u}^{{k}^{2}+k+\vert q\vert \left(2k+1\right)}.\end{equation} \tag{ 27 }$$

While in general there is no entanglement equipartition for these bosonic chains, in the critical limits → 0 one has

$$\begin{equation}{S}_{n}\left(q\right)=\frac{1}{1-n}\enspace \mathrm{log}\enspace \frac{{\mathcal{Z}}_{n}\left(q\right)}{{\mathcal{Z}}_{1}^{n}\left(q\right)}={S}_{n}\left(q=0\right)+\frac{n{{\epsilon}}^{2}{q}^{2}}{2}+O\left({{\epsilon}}^{3}\right),\end{equation} \tag{ 28 }$$

and equipartition is recovered at the leading order in . Again, all the previous formulas are valid for A being a semi-infinite line. The results for a finite interval are obtained exploiting the clustering of the RDM, following, e.g., references [57, 67].

Let us consider a quadratic fermionic system on a two-dimensional square lattice with isotropic hopping between nearest-neighbour sites. It is described by the following Hamiltonian

$$\begin{equation}{H}_{\mathrm{F}\mathrm{F}}=-\frac{1}{2}\sum _{\left\langle \mathbf{i},\mathbf{j}\right\rangle }\left({c}_{\mathbf{i}}^{{\dagger}}{c}_{\mathbf{j}}+{c}_{\mathbf{j}}^{{\dagger}}{c}_{\mathbf{i}}\right)+\mu \sum _{\mathbf{i}}{c}_{\mathbf{i}}^{{\dagger}}{c}_{\mathbf{i}},\end{equation} \tag{ 29 }$$

where μ is the chemical potential for the spinless fermions c_i, with i = (i₁, i₂) a vector identifying a given lattice site, and ⟨i, j⟩ stands for nearest neighbours. Specifically, we consider a set of N coupled identical parallel chains, hence N is the finite length along one direction (say the y-axis). In the other direction, say the x-axis, the system is either infinite or finite with length L. PBC's are imposed along the y-axis. The subsystem A is a (periodic) strip of length ℓ along the x-axis, (see figure 1).

Given the special geometry we consider, we can take the Fourier transform along the transverse y direction. The partial Fourier transforms ${\tilde {c}}_{{j}_{1},r}$ and its inverse are

$$\begin{equation}{\tilde {c}}_{{j}_{1},r}=\frac{1}{\sqrt{N}}\sum _{j=0}^{N-1}{c}_{{j}_{1},j}\enspace {\mathrm{e}}^{-2\pi ijr/N},\quad {c}_{{j}_{1},{j}_{2}}=\frac{1}{\sqrt{N}}\sum _{r=0}^{N-1}{\tilde {c}}_{{j}_{1},r}\enspace {\mathrm{e}}^{2\pi i{j}_{2}r/N},\end{equation} \tag{ 30 }$$

leading to the Hamiltonian in mixed space-momentum representation

$$\begin{equation}{H}_{\mathrm{F}\mathrm{F}}=\sum _{r=0}^{N-1}{H}_{{k}_{y}^{\left(r\right)}}.\end{equation} \tag{ 31 }$$

The operator ${H}_{{k}_{y}^{\left(r\right)}}$ is the Hamiltonian in the ${k}_{y}^{\left(r\right)}=\frac{2\pi r}{N}$ transverse momentum sector:

$$\begin{equation}{H}_{{k}_{y}^{\left(r\right)}}=-\frac{1}{2}\sum _{i=1}^{L}\left({\tilde {c}}_{i,r}^{{\dagger}}{\tilde {c}}_{i+1,r}+\mathrm{h}.\mathrm{c}.\right)+\sum _{i}{\mu }_{r}{\tilde {c}}_{i,r}^{{\dagger}}{\tilde {c}}_{i,r},\end{equation} \tag{ 32 }$$

where

$$\begin{equation}{\mu }_{r}=\mu -\mathrm{cos}\enspace {k}_{y}^{\left(r\right)},\end{equation} \tag{ 33 }$$

and L is the length of the chain along the x-axis. In this way, the Hamiltonian is mapped to a sum of N independent one-dimensional chains with chemical potential μ_r depending on the transverse momentum ${k}_{y}^{\left(r\right)}$.

We focus on the critical regime of the whole 2D system, which is in attained for 0 < μ < 2. In terms of the one dimensional systems, this constraint on μ means that all transverse modes with |μ_r| < 1 are critical, while the others are not. This inequality is satisfied for

$$\begin{equation}r\in {{\Omega}}_{\mu }=\left[0,\frac{\mathrm{arccos}\left(\mu -1\right)N}{2\pi },\left[\cup \right],N\left(1-\frac{\mathrm{arccos}\left(\mu -1\right)}{2\pi }\right),N-1\right].\end{equation} \tag{ 34 }$$

The inner extremes of the intervals are not part of Ω_μ. The case μ = 0 deserves particular attention: when dealing with a finite number of chains, also the mode r = 0 has to be removed from Ω_μ. This difference is irrelevant in the limit N → ∞ when the fraction of critical chains is simply given by $\frac{\mathrm{arccos}\left(\mu -1\right)}{\pi }$.

Since the Hamiltonian is a sum of different sectors, the ground state density matrix factorises and so does the RDM

$$\begin{equation}{\rho }_{\mathrm{A}}=\underset{r=1}{\overset{N}{\bigotimes}}{\rho }_{{k}_{y}^{\left(r\right)}}^{\mathrm{A}}=\underset{r\in {{\Omega}}_{\mu }}{\bigotimes}{\rho }_{{k}_{y}^{\left(r\right)}}^{\mathrm{A}}.\end{equation} \tag{ 35 }$$

In the last equality, we stress that the only relevant modes are the ones corresponding to critical 1D chains. The blocks corresponding to non-critical chains are projectors on the 1D vacuum state, i.e. without fermions. As a consequence, hereafter, we only take into account the gapless modes, which belong to Ω_μ. The RDM ${\rho }_{{k}_{y}^{\left(r\right)}}^{\mathrm{A}}$ of the 1D subsystem associated to the rth mode can be written as [51, 79, 80]

$$\begin{equation}{\rho }_{{k}_{y}^{\left(r\right)}}^{\mathrm{A}}=\mathrm{det}\enspace {C}_{{k}_{y}^{\left(r\right)}}\enspace \mathrm{exp}\left(\sum _{i,j}{\left[\mathrm{log}\left({C}_{{k}_{y}^{\left(r\right)}}^{-1}-1\right)\right]}_{i,j}{\tilde {c}}_{i,r}^{{\dagger}}{\tilde {c}}_{j,r}\right),\end{equation} \tag{ 36 }$$

where the matrix ${C}_{{k}_{y}^{\left(r\right)}}\equiv \langle {\tilde {c}}_{i,r}^{{\dagger}}{\tilde {c}}_{j,r}\rangle$ is the correlation matrix restricted to the rth subsystem A. The entanglement entropy is easily expressed in terms of the eigenvalues of such correlation matrix (see appendix A for further details).

We start by considering the model in the thermodynamic limit in the longitudinal (x) direction, i.e., L → ∞ (figure 1, left panel). For the ground-state of N infinite chains the correlation matrix C of the whole (two-dimensional) subsystem can be written as

$$\begin{equation}\mathbf{C}={\oplus }_{r}{C}_{{k}_{y}^{\left(r\right)}},\end{equation} \tag{ 37 }$$

where ${C}_{{k}_{y}^{\left(r\right)}}$ reads

$$\begin{equation}{C}_{{k}_{y}^{\left(r\right)}}\left(i,j\right)=\frac{\mathrm{sin}\enspace {k}_{r}^{\mathrm{F}}\left(i-j\right)}{\pi \left(i-j\right)},\quad {k}_{r}^{\mathrm{F}}=\mathrm{arccos}\enspace {\mu }_{r},\end{equation} \tag{ 38 }$$

as a function of the Fermi momentum ${k}_{r}^{\mathrm{F}}$ of each rth chain (μ_r is given in equation (33)). This is due to the factorisation of the Hilbert space into the different modes, which corresponds to a block diagonal structure of the correlation matrix: each block is associated to a transverse mode and, as a consequence, to a given 1D ground state.

From the structure of the Hamiltonian in equation (31), the Rényi entropies can be computed by invoking the one-dimensional results discussed in section 2.1: the entanglement entropy is additive on tensor products and therefore decomposes as

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(\underset{r\in {{\Omega}}_{\mu }}{\bigotimes}{\rho }_{{k}_{y}^{\left(r\right)}}^{\mathrm{A}}\right)=\sum _{r\in {{\Omega}}_{\mu }}{S}_{n,r}^{1\mathrm{D}},\quad {S}_{n,r}^{1\mathrm{D}}=\frac{1}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\left(2\ell \enspace \mathrm{sin}\enspace {k}_{r}^{\mathrm{F}}\right)+{{\Upsilon}}_{n}+o\left(1\right),\end{equation} \tag{ 39 }$$

where ϒ_n is in equation (5).

Note that in our setting, the Fermi momentum of each transverse mode-chain can be explicitly written down as

$$\begin{equation}\mathrm{sin}\enspace {k}_{r}^{\mathrm{F}}=\sqrt{1-{\left(\mu -\mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}}.\end{equation} \tag{ 40 }$$

Plugging this relation into equation (39) we get

$$\begin{align}{S}_{n}^{2\mathrm{D}}=\frac{{f}_{N}\left(\mu \right)N}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\left(2\ell \right)+{f}_{N}\left(\mu \right)N{{\Upsilon}}_{n}+\frac{1}{12}\left(1+\frac{1}{n}\right)\\ \quad \quad \quad {\times}\sum _{r\in {{\Omega}}_{\mu }}\mathrm{log}\left[1-{\left(\mu -\mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}\right],\hfill \end{align} \tag{ 41 }$$

where f_N(μ) denotes the fraction of critical modes (i.e., the number of modes belonging to Ω_μ divided by N). It is useful to define also the quantity

$$\begin{equation}{A}_{N}\left(\mu \right)=\frac{1}{2N}\sum _{r\in {{\Omega}}_{\mu }}\mathrm{log}\left[1-{\left(\mu -\mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}\right],\end{equation} \tag{ 42 }$$

so that we have

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{1}{6}\left(1+\frac{1}{n}\right)\left({f}_{N}\left(\mu \right)\mathrm{log}\enspace 2\ell +{A}_{N}\left(\mu \right)\right)N+N{f}_{N}\left(\mu \right){{\Upsilon}}_{n},\end{equation} \tag{ 43 }$$

As aforementioned, when N → ∞ the prefactor of the logarithmic term simply becomes

$$\begin{equation}N{f}_{\infty }\left(\mu \right)=N\frac{\mathrm{arccos}\left(\mu -1\right)}{\pi }.\end{equation} \tag{ 44 }$$

In the left panel of figure 2 we report f_N(μ) as function of N for a few values of μ, showing the approach to N → ∞.

Zoom In Zoom Out Reset image size

In the right panel of figure 2 we report a similar plot for A_N(μ), as function of N for four different values of μ. As N increases, it approaches an asymptotic value that can be explicitly calculated. In fact, in the limit of large N, the sum in equation (41) turns into

$$\begin{align}\frac{1}{2}\sum _{r\in {{\Omega}}_{\mu }}\mathrm{log}\left\vert 1-{\left(\mu -\mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}\right\vert & \to \frac{N}{2\pi }{\int }_{0}^{\mathrm{arccos}\left(\mu -1\right)}\enspace \mathrm{d}x\enspace \mathrm{log}\left(1-{\left(\mu -\mathrm{cos}\left(x\right)\right)}^{2}\right)\\ & \quad \enspace -\frac{1-\mathrm{log}\left(2\pi \sqrt{\mu \left(2-\mu \right)}\right)+\mathrm{log}\enspace N}{2},\hfill \end{align} \tag{ 45 }$$

where we have subtracted the (divergent) contribution from the upper extreme of integration, corresponding to the modes r = {f_∞(μ)N, (1 − f_∞(μ))N}, which are excluded from the sum in the left-hand side (see equation (34)). The explicit computation of the integral gives

$$\begin{align}\hfill & \frac{N}{2\pi }{\int }_{0}^{\mathrm{arccos}\left(\mu -1\right)}\enspace \mathrm{d}x\enspace \mathrm{log}\left(1-{\left(\mu -\mathrm{cos}\left(x\right)\right)}^{2}\right)\hfill \\ \hfill & \quad =-\frac{N}{2\pi }\left[\pi \enspace \mathrm{log}\left(1+4\mu +2{\mu }^{2}-2\left(1+\mu \right)\sqrt{{\mu }^{2}+2\mu }\right)\right.\hfill \\ \hfill & \quad \quad +\mathrm{arccos}\left(\mu -1\right)\mathrm{log}\left(4\left(1+\mu +\sqrt{{\mu }^{2}+2\mu }\right)\right)\hfill \\ \hfill & \quad \quad +\left.\mathrm{Im}\left({\mathrm{L}\mathrm{i}}_{2}\left({\mathrm{e}}^{2\mathrm{i}\mathrm{arccos}\left(\mu -1\right)}\right)+2{\mathrm{L}\mathrm{i}}_{2}\left({\mathrm{e}}^{\mathrm{i}\mathrm{arccos}\left(\mu -1\right)}\left(1+\mu +\sqrt{{\mu }^{2}+2\mu }\right)\right)\right)\right],\hfill \end{align} \tag{ 46 }$$

where Li₂ is the dilogarithmic function ${\mathrm{L}\mathrm{i}}_{2}\equiv {\sum }_{k=1}^{\infty }\frac{{z}^{k}}{{k}^{2}}$. Once again, the case μ = 0 deserves particular attention because also the divergence coming from the lower extreme of integration in (45) has to be subtracted (i.e., the limits μ → 0 and N → ∞ do not commute). Thus, one has to carefully perform a Taylor expansion of the integrand around both extremes of integration. The final result is

$$\begin{equation}{A}_{N}\left(0\right)\to {A}_{\infty }\left(0\right)=-\mathrm{log}\enspace 2-\frac{2-2\enspace \mathrm{log}\enspace \pi +2\enspace \mathrm{log}\enspace N}{N}.\end{equation} \tag{ 47 }$$

The logarithmic correction for small values of N is evident in figure 2 for all values of μ, but it is more pronounced for μ = 0, as clear from the analytic expressions. Hence the total entropy for large N is

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{{f}_{\infty }\left(\mu \right)N}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\enspace 2\ell +\frac{N}{6}\left(1+\frac{1}{n}\right){A}_{\infty }\left(\mu \right)+{f}_{\infty }\left(\mu \right)N{{\Upsilon}}_{n},\end{equation} \tag{ 48 }$$

which we recall is valid at order o(ℓ⁰) and O(N⁰). We will see that to have a good agreement with numerical data at finite but large ℓ, it is needed to keep the log N contribution in A_N(μ).

When both ℓ and N are large, it is useful to look at the special case of the subsystem A being a square strip with N = ℓ, when equation (48) is rewritten as

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{{f}_{\infty }\left(\mu \right)}{6}\left(1+\frac{1}{n}\right)\ell \enspace \mathrm{log}\enspace 2\ell +\ell \frac{1}{6}\left(1+\frac{1}{n}\right){A}_{\infty }\left(\mu \right)+{f}_{\infty }\left(\mu \right)\ell {{\Upsilon}}_{n}+O\left({\ell }^{0}\right).\end{equation} \tag{ 49 }$$

(Actually, any choice of N and ℓ proportional to each other, N = aℓ, would be equivalent, with just an overall factor a, but for simplicity let us just think to a = 1.) Let us briefly comment equation (49). It shows the expected logarithmic correction to the area law and our derivation gives a clearer understanding of such behaviour: it is a simple consequence of the fact that we are dealing with an extensive a number of critical chains, i.e. proportional to N = ℓ, whose entropy obeys a logarithmic scaling so that each of them contributes proportionally to log ℓ to the total entropy. Moreover, it also agrees with the result obtained by the application of Widom conjecture (see, e.g., [11–13]) that provides an explicit formula for the prefactor of the leading term of the entanglement entropy of free fermions in any dimension, i.e., S₁ = Cℓ^d−1 log ℓ + O(ℓ^d−1), with C given by

$$\begin{equation}C=\frac{1}{12\left(2\pi \right)}{\int }_{\partial {\Lambda}}{\int }_{\partial {\Gamma}\left(\mu \right)}\vert {n}_{x}\cdot {n}_{p}\vert \mathrm{d}{S}_{x}\mathrm{d}{S}_{p},\end{equation} \tag{ 50 }$$

where Λ is the considered subsystem with volume normalised to one, Γ(μ) is the volume in momentum space enclosed by the Fermi surface, n_p, n_x are the unit normals to the boundaries of these volumes and the integration is carried over the surface of both domains. In the case of interest for this paper, given the compactification along the y direction, the Fermi surface is defined by the solutions of S_p = 0 (with S_p = μ − cos k_x − cos k_y). By performing the line integrals, equation (50) becomes

$$\begin{equation}C=\frac{\mathrm{arccos}\left(\mu -1\right)}{3\pi },\end{equation} \tag{ 51 }$$

in agreement with equation (49). We stress that while the leading terms in the two approaches are identical, the dimensional reduction provides an explicit prediction also for the subleading term proportional to ℓ, as in equation (49), which cannot be derived by Widom conjecture.

3.2.1. Some generalisations

The same approach is straightforwardly adapted to the computation of the entanglement entropies in the case of Dirichlet (open) boundary conditions (DBC's) along the transverse direction (y-axis), i.e., imposing ${c}_{{j}_{1},0}={c}_{{j}_{1},N}=0$. Although these boundary conditions break the translational invariance in the transverse direction, one can use the Fourier sine transform (rather than the standard one). The only final difference is that the set of modes Ω_μ in (34) corresponding to critical chains will now start from r = 1 (instead of r = 0). The same strategy applies when the total system is a finite block of L sites along the x-direction, with PBC's (see the right panel in figure 1). In this case, the only difference is that the scaling of the one-dimensional Rényi entropies for a system with PBC's reads

$$\begin{equation}{S}_{n,j}^{1\mathrm{D}}=\frac{1+n}{6n}\enspace \mathrm{log}\left(\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\mathrm{sin}\enspace {k}_{j}^{\mathrm{F}}\right)+{{\Upsilon}}_{n}.\end{equation} \tag{ 52 }$$

Therefore, we have for any finite N

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{n+1}{6n}\left[{f}_{N}\left(\mu \right)\mathrm{log}\left(\frac{L}{2\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)+{A}_{N}\left(\mu \right)\right]N+{f}_{N}\left(\mu \right)N{{\Upsilon}}_{n},\end{equation} \tag{ 53 }$$

and similarly for large N with f_N(μ) → f_∞(μ) and A_N(μ) → A_∞(μ).

3.2.2. Numerical checks

We now benchmark the results for the total entropies against exact numerical calculations obtained by the free-fermion techniques reported in the appendix A. In figure 3 we report the numerical data of the Rényi entropies for different values of the index n and chemical potential μ, both for infinite (panel (a)) and finite (panel (b)) system size. We fix the transverse direction N to be equal to the longitudinal subsystem length ℓ, so that the subsystem A is a square with PBC in the transverse direction. We also properly choose the values of μ and ℓ such that ℓf_ℓ(μ) is an integer number to eliminate effects due to partial fillings of modes. The theoretical predictions for the leading scaling in equations (49) and (53) are also reported for comparison. These include both the leading term and the subleading one proportional to the area (2ℓ) between the subsystem A and the rest of the system.

Zoom In Zoom Out Reset image size

It is evident that the analytical results correctly describe the data. We also report (as dashed lines) the sole leading universal behaviour ∝ℓ log ℓ: this universal term alone does not match the data for these values of ℓ, highlighting the importance of the subleading terms ∝ℓ that we calculated analytically here for the first time. In the panel (c) of the same figure, we plot the data for the von Neumann entropy where we subtracted the leading term f_ℓ(μ)ℓ log ℓ to show the non-universal subleading terms found in equation (49) alone.

In figure 3, subleading oscillating corrections for n ≠ 1 are visible. These are easily understood as a consequence of the well studied unusual corrections to the scaling in 1D [49, 81–83], which are present for generic bipartitions (see e.g. [14]). Anyhow, in our special case we can exploit the dimensional reduction also to derive exact predictions for these corrections in 2D. In 1D, for the tight-binding model, they behave like [49, 82]

$$\begin{equation}{d}_{n}^{1\mathrm{D}}\left(\ell \right)\equiv {S}_{n}^{1\mathrm{D}}\left(\ell \right)-{S}_{n}^{1\mathrm{D},\left(0\right)}\left(\ell \right)={f}_{n}\enspace \mathrm{cos}\left(2\ell {k}_{\mathrm{F}}\right)\vert 2\ell \enspace \mathrm{sin}\enspace {k}_{\mathrm{F}}{\vert }^{-2/n},\end{equation} \tag{ 54 }$$

where ${S}_{n}^{1\mathrm{D},\left(0\right)}\left(\ell \right)$ is the (leading) prediction from the Fisher–Hartwig conjecture in equation (4), (which contains both the CFT prediction and the non universal constant term ϒ_n) and the amplitude is

$$\begin{equation}{f}_{n}=\frac{2}{1-n}\frac{{{\Gamma}}^{2}\left(\left(1+{n}^{-1}\right)/2\right)}{{{\Gamma}}^{2}\left(\left(1-{n}^{-1}\right)/2\right)}.\end{equation} \tag{ 55 }$$

By dimensional reduction we have that each chain with Hamiltonian ${H}_{{k}_{y}^{\left(r\right)}}$ in equation (31) has corrections given by equation (54) with the appropriate Fermi momentum. Summing over the contributions given by each mode in equation (54), we obtain (for the case N = ℓ of the figure)

$$\begin{equation}{d}_{n}^{2\mathrm{D}}\left(\ell \right)\equiv {S}_{n}^{2\mathrm{D}}\left(\ell \right)-{S}_{n}^{2\mathrm{D},\left(0\right)}\left(\ell \right)={f}_{n}\sum _{j\in {{\Omega}}_{\mu }}\mathrm{cos}\left(2\ell {k}_{j}^{\mathrm{F}}\right)\vert 2\ell \enspace \mathrm{sin}\enspace {k}_{j}^{\mathrm{F}}{\vert }^{-2/n},\end{equation} \tag{ 56 }$$

where ${S}_{n}^{2\mathrm{D},\left(0\right)}$ stands for the leading terms in equation (43). The accuracy of these subleading corrections is tested against numerical data in figure 4 for three different values of the chemical potential, finding perfect agreement. Several frequencies corresponding to the various ${k}_{j}^{\mathrm{F}}$ are clearly visible in the figure.

Zoom In Zoom Out Reset image size

The same dimensional reduction technique can further be used to compute the symmetry resolved entanglement entropies. Indeed, from equation (29) the particle number $Q={\sum }_{\mathbf{i}}{c}_{\mathbf{i}}^{{\dagger}}{c}_{\mathbf{i}}$ is a conserved U(1) charge of the model in arbitrary dimension. The strategy is exactly as before: we consider a finite system in the transverse direction with PBC and so reduce to a one-dimensional problem for the charged moments and then, via Fourier transform, we get the symmetry resolved entropies.

3.3.1. Charged moments

Because of the factorisation of the RDM (35) and of the additivity of the conserved charge, we can rewrite

$$\begin{equation}{\rho }_{\mathrm{A}}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{Q}_{\mathrm{A}}\alpha }=\underset{r\in {{\Omega}}_{\mu }}{\bigotimes}{\rho }_{A,{k}_{y}^{\left(r\right)}}^{n}\enspace {\mathrm{e}}^{\mathrm{i}{Q}_{\mathrm{A}}^{\left(r\right)}\alpha },\end{equation} \tag{ 57 }$$

where ${Q}_{\mathrm{A}}^{\left(r\right)}$ is the charge operator restricted to the rth transverse mode. This factorisation allows us to rewrite in terms of the one-dimensional results for the charged moment

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)=\sum _{r\in {{\Omega}}_{\mu }}\mathrm{log}\enspace {Z}_{n,r}^{1\mathrm{D}}\left(\alpha \right),\end{equation} \tag{ 58 }$$

and, using the explicit 1D result equation (19), the sum is performed as

$$\begin{equation}\mathrm{log}{Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq \mathrm{i}\overline{q}\alpha -\left[\frac{1}{6}\left(n-\frac{1}{n}\right)+\frac{2}{n}{\left(\frac{\alpha }{2\pi }\right)}^{2}\right]\left({f}_{N}\left(\mu \right)\mathrm{log}2\ell +{A}_{N}\left(\mu \right)\right)N+N{f}_{N}\left(\mu \right){\Upsilon}\left(n,\alpha \right).\end{equation} \tag{ 59 }$$

The first term in equation (59) is purely imaginary and it is the average number of particle within A, for large N explicitly given by

$$\begin{equation}\overline{q}=\frac{N\ell }{{\pi }^{2}}{\int }_{0}^{\mathrm{arccos}\left(\mu -1\right)}\mathrm{d}x\enspace \mathrm{arccos}\left(\mu -\mathrm{cos}\enspace x\right).\end{equation} \tag{ 60 }$$

It is extensive in the subsystem volume (Nℓ), as it should, and at half-filling, μ = 0, it reproduces the simple result $\overline{q}=N\ell /2$.

In equation (59), it is useful to write ϒ(n, α) as

$$\begin{equation}{\Upsilon}\left(n,\alpha \right)={\Upsilon}\left(n\right)+\gamma \left(n\right){\alpha }^{2}+{\epsilon}\left(n,\alpha \right),\quad {\epsilon}\left(n,\alpha \right)=O\left({\alpha }^{4}\right),\end{equation} \tag{ 61 }$$

where

$$\begin{equation}\gamma \left(n\right)=\frac{ni}{4}{\int }_{-\infty }^{\infty }\mathrm{d}w\left[{\mathrm{tanh}}^{3}\left(\pi nw\right)-\mathrm{tanh}\left(\pi nw\right)\right]\mathrm{log}\enspace \frac{{\Gamma}\left(\frac{1}{2}+\mathrm{i}w\right)}{{\Gamma}\left(\frac{1}{2}-\mathrm{i}w\right)}.\end{equation} \tag{ 62 }$$

In reference [45] it has been shown that the quadratic approximation of equation (61) is appropriate for many of applications since (n, α) ≪ γ(n)α². In particular, this approximation allows us for an explicit analytic computation of the symmetry resolved moments ${\mathcal{Z}}_{n}\left(q\right)$. Therefore, hereafter we will keep only the terms up to O(α²) and we rewrite (59) in the compact form as:

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq \mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(0\right)+\mathrm{i}\overline{q}\alpha -{\alpha }^{2}\left({\mathcal{B}}_{n}{f}_{N}\left(\mu \right)\mathrm{log}\enspace 2\ell +{\mathcal{C}}_{n}\right)N,\end{equation} \tag{ 63 }$$

with

$$\begin{equation}\begin{aligned}\hfill & {\mathcal{B}}_{n}=\frac{1}{2{\pi }^{2}n},\hfill \\ \hfill & {\mathcal{C}}_{n}=\frac{{A}_{N}\left(\mu \right)}{2{\pi }^{2}n}-{f}_{N}\left(\mu \right)\gamma \left(n\right).\hfill \end{aligned}\end{equation} \tag{ 64 }$$

In figure 5 we report the numerical data both for the real and the imaginary part of $\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ for different values of n and α. Here the system is an infinite cylinder of circumference ℓ and the subsystem A is again a periodic square strip with longitudinal length equal to ℓ. Here and throughout this section, the values of μ and ℓ are chosen such that ℓf_ℓ(μ) is an integer number. Moreover, when μ = 0, we focus on the case ℓ even. The theoretical prediction in equation (59) is also reported for comparison, showing that the analytical result correctly describes the data as long as |α| < π (as well known already in 1D, see e.g. [45, 46]). The oscillating corrections to the scaling become relevant when α moves close to ±π. The reason is that some terms in the generalised Fisher–Hartwig approach become larger and equation (19) is not a good approximation at the considered intermediate values of ℓ [45, 46]. In the figure it is evident that these oscillations arise also for n = 1, contrarily to what happens for α = 0.

Zoom In Zoom Out Reset image size

We can also study the subleading oscillatory behaviour exploiting the one-dimensional results [45] (valid for −π < α < π), i.e.,

$$\begin{align}\hfill {d}_{n}^{1\mathrm{D},j}\left(\alpha ,\ell \right)& \equiv \mathrm{log}\enspace {Z}_{n}^{1\mathrm{D}}\left(\alpha \right)-\mathrm{log}\enspace {Z}_{n}^{1\mathrm{D},\left(0\right)}\left(\alpha \right)=\hfill \\ \hfill & ={\mathrm{e}}^{-2\mathrm{i}{k}_{j}^{\mathrm{F}}\ell }{\left(2\ell \enspace \mathrm{sin}\enspace {k}_{j}^{\mathrm{F}}\right)}^{-\frac{2}{n}\left(1-\frac{\alpha }{\pi }\right)}{f}_{n}^{\left(1\right)}\left(\alpha \right)+{\mathrm{e}}^{2\mathrm{i}{k}_{j}^{\mathrm{F}}\ell }{\left(2\ell \enspace \mathrm{sin}\enspace {k}_{j}^{\mathrm{F}}\right)}^{-\frac{2}{n}\left(1+\frac{\alpha }{\pi }\right)}{f}_{n}^{\left(2\right)}\left(\alpha \right)+\dots ,\hfill \end{align} \tag{ 65 }$$

where ${Z}_{n}^{1\mathrm{D},\left(0\right)}\left(\alpha \right)$ is the (leading) prediction of the generalised Fisher–Hartwig conjecture, equation (19), and

$$\begin{equation}{f}_{n}^{\left(1\right)}\left(\alpha \right)=\frac{{{\Gamma}}^{2}\left(\frac{1}{2}+\frac{1}{2n}-\frac{\alpha }{2\pi n}\right)}{{{\Gamma}}^{2}\left(\frac{1}{2}-\frac{1}{2n}+\frac{\alpha }{2\pi n}\right)},\quad {f}_{n}^{\left(2\right)}\left(\alpha \right)=\frac{{{\Gamma}}^{2}\left(\frac{1}{2}+\frac{1}{2n}+\frac{\alpha }{2\pi n}\right)}{{{\Gamma}}^{2}\left(\frac{1}{2}-\frac{1}{2n}-\frac{\alpha }{2\pi n}\right)}.\end{equation} \tag{ 66 }$$

In 2D, the subleading oscillatory behaviour is easily obtained summing the contributions for each mode given by equation (65), resulting for N = ℓ in

$$\begin{equation}{d}_{n}^{2\mathrm{D}}\left(\alpha ,\ell \right)\equiv \mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)-\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D},\left(0\right)}\left(\alpha \right)=\sum _{j\in {{\Omega}}_{\mu }}{d}_{n}^{1\mathrm{D},j}\left(\ell ,\alpha \right),\end{equation} \tag{ 67 }$$

which is compared with numerical data in figure 6, finding perfect agreement. For μ = 0, the oscillatory behaviour of the imaginary part of the charged moments vanishes.

Zoom In Zoom Out Reset image size

A simple but interesting generalisation of the calculation we just presented concerns the geometry of a torus as depicted in the right of figure 1. The longitudinal size of the system is L. The charged moments are again obtained by summing up the contribution of the different transverse modes as critical 1D chains, using the finite size form with the chord length. Summing up the contributions of the the transverse modes we get

$$\begin{align}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)& \simeq \mathrm{i}\alpha \overline{q}-\left[\frac{1}{6}\left(n-\frac{1}{n}\right)+\frac{2}{n}{\left(\frac{\alpha }{2\pi }\right)}^{2}\right]N\left[{f}_{N}\left(\mu \right)\mathrm{log}\left[\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right]\right.\\ & \quad \left.+{A}_{N}\left(\mu \right)\right]+N{f}_{N}\left(\mu \right){\Upsilon}\left(n,\alpha \right).\hfill \end{align} \tag{ 68 }$$

The accuracy of this prediction is tested in figure 7 against exact numerical calculations.

Zoom In Zoom Out Reset image size

3.3.2. Symmetry resolution

We now can compute the Fourier transform ${\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)$ of the charged moments using the leading order terms of ${Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ taking into account the effect of the non-universal pieces. This Fourier transform is

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)={\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }{Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq {Z}_{n}^{2\mathrm{D}}\left(0\right){\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}\left(q-\overline{q}\right)\alpha -{\alpha }^{2}{b}_{n}},\end{equation} \tag{ 69 }$$

where the coefficient of the quadratic term is

$$\begin{equation}{b}_{n}={\mathcal{B}}_{n}{f}_{N}\left(\mu \right)\mathrm{log}\enspace 2\ell +{\mathcal{C}}_{n}N.\end{equation} \tag{ 70 }$$

For large subsystem size ℓ and/or N, we can treat the integral by means of the saddle point approximation and use as domain of integration [−∞, +∞], getting

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)\simeq {Z}_{n}^{2\mathrm{D}}\left(0\right){\mathrm{e}}^{-\frac{{\left(q-\overline{q}\right)}^{2}}{4N\left({\mathcal{B}}_{n}{f}_{N}\left(\mu \right)\mathrm{log}2\ell +{\mathcal{C}}_{n}\right)}}\sqrt{\frac{1}{4\pi N\left({\mathcal{B}}_{n}{f}_{N}\left(\mu \right)\mathrm{log}\enspace 2\ell +{\mathcal{C}}_{n}\right)}},\end{equation} \tag{ 71 }$$

where ${Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ is given in equation (59) and we report it again for completeness in a coincise form for α = 0

$$\begin{equation}{Z}_{n}^{2\mathrm{D}}\left(0\right)={\left({\left(2\ell \right)}^{{f}_{N}\left(\mu \right)}{\mathrm{e}}^{{A}_{N}\left(\mu \right)}\right)}^{-\frac{1}{6}\left(n-\frac{1}{n}\right)N}\enspace {\mathrm{e}}^{{\Upsilon}\left(n\right)N{f}_{N}\left(\mu \right)}.\end{equation} \tag{ 72 }$$

In full analogy to the 1D case, the probability distribution functions given by these moments are still Gaussian with mean $\overline{q}$ and variance that for large ℓ and N grows as $\sqrt{N\enspace \mathrm{log}\enspace \ell }$. An equivalent result was already obtained in references [46, 62] for Fermi gases in arbitrary dimension using the Widom's conjecture. The novelty of this formula is an exact prediction for the coefficient ${\mathcal{C}}_{n}$ that renormalises the variance at order O(ℓ) and, as we will see, will play a crucial role for an accurate computation of the symmetry resolved entropies.

Let us briefly discuss the terms that have been neglected in the derivation of equation (71) which are the same as in 1D [45]. The main approximation is to ignore (n, α) in equation (62) which induces a correction going like 1/(N log ℓ). The subleading corrections to ${Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ in equation (56) only induce power-law corrections and are subdominant compared to one above. Finally the corrections coming from having replaced the extremes of integration ±π with ±∞ are really small: they decay as ${\mathrm{e}}^{-{\pi }^{2}{b}_{n}}/{b}_{n}$, i.e., exponentially in N.

The accuracy of equation (71) is checked for different values of μ in figure 8 where we report the numerically calculated Fourier transforms and the analytical prediction. It is evident from the data in the main frames and in the insets that both the ℓ and the q dependence of ${\mathcal{Z}}_{n}\left(q\right)$ is perfectly captured by our approximation.

Zoom In Zoom Out Reset image size

With these ingredients at our disposal, we are ready to compute the asymptotic behaviour of the symmetry resolved entanglement, given by

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(q\right)=\frac{1}{1-n}\enspace \mathrm{log}\left[\frac{{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)}{{\mathcal{Z}}_{1}^{2\mathrm{D}}{\left(q\right)}^{n}}\right]\simeq \frac{1}{1-n}\enspace \mathrm{log}\enspace \frac{{Z}_{n}^{2\mathrm{D}}\left(0\right)}{{\left({Z}_{1}^{2\mathrm{D}}\left(0\right)\right)}^{n}}\frac{{\mathrm{e}}^{-\frac{{\left(q-\overline{q}\right)}^{2}}{4{b}_{n}}}}{{\mathrm{e}}^{-\frac{n{\left(q-\overline{q}\right)}^{2}}{4{b}_{1}}}}\frac{{\left(4\pi {b}_{n}\right)}^{-1/2}}{{\left(4\pi {b}_{1}\right)}^{-n/2}}.\end{equation} \tag{ 73 }$$

After some simple algebra, we can write

$$\begin{align}\hfill {S}_{n}^{2\mathrm{D}}\left(q\right)& ={S}_{n}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left(\frac{2N}{\pi }\left({f}_{\infty }\left(\mu \right)\mathrm{log}\left(2\ell \right)+{f}_{\infty }\left(\mu \right){\delta }_{n}+{A}_{\infty }\left(\mu \right)\right)\right)+\frac{\mathrm{log}\enspace n}{2\left(1-n\right)}\hfill \\ \hfill & \quad +\enspace {\left(q-\overline{q}\right)}^{2}{\pi }^{4}\frac{n}{1-n}\frac{\left(\gamma \left(1\right)-n\gamma \left(n\right)\right)}{N{\left[{f}_{\infty \left(\mu \right)\mathrm{log}\left(2\ell \right)+{f}_{\infty }\left(\mu \right){\kappa }_{n}}+{A}_{\infty }\left(\mu \right)\right]}^{2}}+\cdots \enspace ,\hfill \end{align} \tag{ 74 }$$

where ${S}_{n}^{2\mathrm{D}}$ is the total Rényi entropy,

$$\begin{equation}{\delta }_{n}=-\frac{2{\pi }^{2}n\left(\gamma \left(n\right)-\gamma \left(1\right)\right)}{1-n},\end{equation} \tag{ 75 }$$

and

$$\begin{equation}{\kappa }_{n}=-{\pi }^{2}\left(\gamma \left(1\right)+n\gamma \left(n\right)\right).\end{equation} \tag{ 76 }$$

equation (74) provides the leading behaviour for large ℓ and N as well as the non-universal additive constants, and a q-dependent subleading correction which scales as N⁻¹(log ℓ)⁻². Such correction provides the first term in the expansion for large ℓ and N which depends on the symmetry sector. As in the corresponding 1D calculation [45], it can be calculated from the subleading terms of the variance of ${\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)$, in particular the additive non-universal constant ${\mathcal{C}}_{n}$ in equation (64). So while few leading terms satisfy the equipartition of entanglement, we can precisely identify the first term that breaks it. Taking now the limit for n → 1 of (74), we get the von Neumann entropy

$$\begin{align}\hfill {S}_{1}^{2\mathrm{D}}\left(q\right)& ={S}_{1}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left(\frac{2N}{\pi }\left({f}_{\infty }\left(\mu \right)\mathrm{log}\left(2\ell \right)+{f}_{\infty }\left(\mu \right){\delta }_{1}+{A}_{\infty }\left(\mu \right)\right)\right)-\frac{1}{2}+\hfill \\ \hfill & \quad +\enspace {\left(q-\overline{q}\right)}^{2}{\pi }^{4}\frac{\left(\gamma \left(1\right)+{\gamma }^{\prime }\left(1\right)\right)}{N{\left[{f}_{\infty }\left(\mu \right)\mathrm{log}\left(2\ell \right)+{f}_{\infty }\left(\mu \right){\kappa }_{1}+{A}_{\infty }\left(\mu \right)\right]}^{2}}+\cdots \hfill \end{align} \tag{ 77 }$$

These predictions for the symmetry resolved entanglement are compared with the numerical data in figure 9. In the left panel we consider the scaling with ℓ of ${S}_{n}\left(\overline{q}\right)$ and it is evident that the numerical data perfectly match with the theoretical prediction in equations (74) and (77). The corrections in $\left(q-\overline{q}\right)$ are suppressed as 1/(N(log ℓ)²) and the curves in the right panel seem to be on top of each other on the scale of the plot. In order to appreciate their distance, in the inset we report the differences with ${S}_{n}\left(\overline{q}\right)$ (focussing on n = 1) and we show that they are well described by our prediction. The agreement is excellent even for relatively small values of ℓ and N of the order of 20.

Zoom In Zoom Out Reset image size

In the one-dimensional case, the corrections of order log(log ℓ) coming from the symmetry resolved entanglement exactly cancel in the total entanglement entropy, when summing to the fluctuation entanglement as in equation (16). The same occurs also in 2D. In fact, the fluctuation entanglement in our case is given by

$$\begin{align}{S}^{2\mathrm{D},f}=-{\int }_{q}{\mathcal{Z}}_{1}\left(q\right)\mathrm{log}\enspace {\mathcal{Z}}_{1}\left(q\right)\simeq \frac{1}{2}\left(1+\mathrm{log}\enspace 4\pi {b}_{1}\right)\\ \quad \quad \enspace =\frac{1}{2}+\frac{1}{2}\enspace \mathrm{log}\left(\frac{2}{\pi }N{f}_{\infty }\left(\mu \right)\mathrm{log}\enspace \ell \right)+O\left({\left(\mathrm{log}\enspace \ell \right)}^{-1}\right).\hfill \end{align} \tag{ 78 }$$

From this equation, it is clear that the first two leading terms in S^2D,f cancel exactly with the corresponding ones in the symmetry resolved entanglement in equation (77).

We quickly discuss now what happens in the case of the torus geometry: we only report the final results, since the calculations are only a slight modification of the previous ones. The Fourier transform of equation (68) is

$$\begin{equation}{\mathcal{Z}}_{n}\left(q\right)\simeq {Z}_{n}^{2\mathrm{D}}\left(0\right)\frac{{\mathrm{e}}^{-\frac{{\left(q-\overline{q}\right)}^{2}}{4N\left[\frac{1}{2{\pi }^{2}n}\left({f}_{N}\left(\mu \right)\mathrm{log}\left(\frac{2L}{\pi }\mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)+{A}_{N}\left(\mu \right)\right)-{f}_{N}\left(\mu \right)\gamma \left(n\right)\right]}}}{\sqrt{4N\pi \left[\frac{1}{2{\pi }^{2}n}\left({f}_{N}\left(\mu \right)\mathrm{log}\left(\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)+{A}_{N}\left(\mu \right)\right)-{f}_{N}\left(\mu \right)\gamma \left(n\right)\right]}}\end{equation} \tag{ 79 }$$

where

$$\begin{equation}{Z}_{n}^{2\mathrm{D}}\left(0\right)={\mathrm{e}}^{{\Upsilon}\left(n\right)N{f}_{N}\left(\mu \right)}{\left({\mathrm{e}}^{{A}_{N}\left(\mu \right)}{\left(\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)}^{{f}_{N}\left(\mu \right)}\right)}^{-\frac{1}{6}\left(n-\frac{1}{n}\right)N}.\end{equation} \tag{ 80 }$$

The symmetry resolved entropies are then easily worked out as

$$\begin{align}\hfill {S}_{n}^{2\mathrm{D}}\left(q\right)& ={S}_{n}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left[\frac{2N}{\pi }\left({f}_{\infty }\left(\mu \right)\mathrm{log}\left(\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)+{f}_{\infty }\left(\mu \right){\delta }_{n}+{A}_{\infty }\left(\mu \right)\right)\right]+\frac{\mathrm{log}\enspace n}{2\left(1-n\right)}\hfill \\ \hfill & \quad +\enspace {\left(q-\overline{q}\right)}^{2}{\pi }^{4}\frac{n}{1-n}\frac{\left(\gamma \left(1\right)-n\gamma \left(n\right)\right)}{N{\left[{f}_{\infty }\left(\mu \right)\mathrm{log}\left(\frac{2L}{\pi }\enspace \mathrm{sin}\left(\frac{\pi \ell }{L}\right)\right)+{f}_{\infty }\left(\mu \right){\kappa }_{n}+{A}_{\infty }\left(\mu \right)\right]}^{2}}+\cdots \enspace .\hfill \end{align} \tag{ 81 }$$

These results are tested with numerics in figure 10.

Zoom In Zoom Out Reset image size

Let us examine a two-dimensional system of L × N real coupled oscillators, where L and N are the lengths along the x- (longitudinal) and y-direction (transverse), respectively (see figure 1). The 2D Hamiltonian describing a real 2D square lattice of harmonic oscillators is

$$\begin{equation}{H}_{\mathrm{B}}=\frac{1}{2}\sum _{x=1}^{L}\sum _{y=1}^{N}\left[{p}_{x,y}^{2}+{\omega }_{0}^{2}{q}_{x,y}^{2}+{\kappa }_{x}{\left({q}_{x+1,y}-{q}_{x,y}\right)}^{2}+{\kappa }_{y}{\left({q}_{x,y+1}-{q}_{x,y}\right)}^{2}\right],\end{equation} \tag{ 82 }$$

where q_x,y, p_x,y and ω₀ are coordinate, momentum and self-frequency of the oscillator at site (x, y) while κ_x and κ_y are the nearest-neighbour couplings. As in the one-dimensional case, the 2D lattice of complex oscillators is

$$\begin{equation}{H}_{\mathrm{C}\mathrm{B}}\left({p}^{\left(1\right)}+\mathrm{i}{p}^{\left(2\right)},{q}^{\left(1\right)}+\mathrm{i}{q}^{\left(2\right)}\right)={H}_{\mathrm{B}}\left({p}^{\left(1\right)},{q}^{\left(1\right)}\right)+{H}_{\mathrm{B}}\left({p}^{\left(2\right)},{q}^{\left(2\right)}\right).\end{equation} \tag{ 83 }$$

If we define $\mathbf{p}=\frac{{p}^{\left(1\right)}+\mathrm{i}{p}^{\left(2\right)}}{\sqrt{2}}$ and $\mathbf{q}=\frac{{q}^{\left(1\right)}+\mathrm{i}{q}^{\left(2\right)}}{\sqrt{2}}$, equation (83) becomes

$$\begin{align}\hfill {H}_{\mathrm{C}\mathrm{B}}& =\sum _{x=1}^{L}\sum _{y=1}^{N}\left[{\mathbf{p}}_{x,y}^{{\dagger}}{\mathbf{p}}_{x,y}+{\omega }_{0}^{2}{\mathbf{q}}_{x,y}^{{\dagger}}{\mathbf{q}}_{x,y}++{\kappa }_{x}{\left({\mathbf{q}}_{x+1,y}-{\mathbf{q}}_{x,y}\right)}^{{\dagger}}\left({\mathbf{q}}_{x+1,y}-{\mathbf{q}}_{x,y}\right)\right.\hfill \\ \hfill & \quad \left.+\enspace {\kappa }_{y}{\left({\mathbf{q}}_{x,y+1}-{\mathbf{q}}_{x,y}\right)}^{{\dagger}}\left({\mathbf{q}}_{x,y+1}-{\mathbf{q}}_{x,y}\right)\right].\hfill \end{align} \tag{ 84 }$$

Imposing PBC's along the y-direction, we can exploit the translational invariance and use Fourier transform in the transverse direction, to get the mixed space-momentum representation

$$\begin{equation}{\mathbf{q}}_{x,y}=\frac{1}{\sqrt{N}}\sum _{r=0}^{N-1}{\tilde {q}}_{x,r}\enspace {\mathrm{e}}^{2\pi \mathrm{i}ry/N},\end{equation} \tag{ 85 }$$

and similarly for p_x,y. We set κ_y = κ_x = 1 to shorten the notation (indeed they can be absorbed by a canonical transformation). The Hamiltonian (83) then becomes

$$\begin{equation}{H}_{\mathrm{C}\mathrm{B}}=\sum _{x=1}^{L}\sum _{r=0}^{N-1}{\tilde {p}}_{x,r}^{{\dagger}}{\tilde {p}}_{x,r}+{\omega }_{r}^{2}{\tilde {q}}_{x,r}^{{\dagger}}{\tilde {q}}_{x,r}+{\left({\tilde {q}}_{x+1,r}-{\tilde {q}}_{x,r}\right)}^{{\dagger}}\left({\tilde {q}}_{x+1,r}-{\tilde {q}}_{x,r}\right).\end{equation} \tag{ 86 }$$

where

$$\begin{equation}{\omega }_{r}^{2}={\omega }_{0}^{2}+4\enspace {\mathrm{sin}}^{2}\enspace \frac{\pi r}{N}.\end{equation} \tag{ 87 }$$

Because of the additivity of the independent transverse chain modes in equation (86), the entanglement entropy can be computed by using the 1D results, as we did for free fermions.

As already discussed in section 2.1, the ground-state reduced density matrices of each 1D chain is obtained by means of corner transfer matrices, for the bipartition of the infinite chain in two halves. Therefore, the entanglement spectrum associated to the r-mode/chain along the y direction is given by equation (8), specialised to the frequency ω_r, i.e., the eigenvalues of the entanglement Hamiltonian are now given by

$$\begin{equation}{{\epsilon}}_{j}^{\left(r\right)}={{\epsilon}}_{r}\left(2j+1\right),\end{equation} \tag{ 88 }$$

where the energy levels are

$$\begin{equation}{{\epsilon}}_{r}=\frac{\pi I\left(\sqrt{1-{\kappa }_{r}^{2}}\right)}{I\left({\kappa }_{r}\right)},\quad {\kappa }_{r}=\frac{1}{2}\left(2+{\omega }_{r}^{2}-{\omega }_{r}\sqrt{4+{\omega }_{r}^{2}}\right).\end{equation} \tag{ 89 }$$

The parameter κ_r is obtained by solving the equation ${\omega }_{r}^{2}={\left(1-{\kappa }_{r}\right)}^{2}/{\kappa }_{r}$.

In our semi infinite strip, each mode gives a contribution to the entanglement entropy which can be computed through the CTM approach. Since they are independent, such contributions simply add up leading to

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{2}{1-n}\sum _{r=0}^{N-1}\sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){{\epsilon}}_{r}}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}}\right]\right).\end{equation} \tag{ 90 }$$

This result is valid for arbitrary integer N. We can now take the limit of large transverse direction. The sum becomes an integral in ζ = r/N, we can write (_r → (ζ))

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{2N}{1-n}{\int }_{0}^{1}\mathrm{d}\zeta \sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}\left(\zeta \right)}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}\left(\zeta \right)}\right]\right),\end{equation} \tag{ 91 }$$

and in the limit n → 1

$$\begin{equation}{S}_{1}^{2\mathrm{D}}=2N{\int }_{0}^{1}\mathrm{d}\zeta \sum _{j=0}^{\infty }\left(\frac{\left(2j+1\right){\epsilon}\left(\zeta \right)}{{\mathrm{e}}^{\left(2j+1\right){\epsilon}\left(\zeta \right)}-1}-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}\left(\zeta \right)}\right]\right).\end{equation} \tag{ 92 }$$

These results for the entanglement entropies are numerically tested in panel (a) of figure 11 using the free-boson techniques reported in the appendix A. The considered subsystem is a strip, periodic in the transverse direction (hence of length N) and of longitudinal size equal to ℓ, but such that ℓ is much larger than the correlation length of the system (of order ${\omega }_{0}^{-1}$) and so the entanglement is just the double of the one for a semi-infinite subsystem. We have fixed ℓ = 50 which is large enough for the considered values of ω₀. Equations (91) and (92) perfectly predict the prefactor of the area-law term, in all cases when the thermodynamic limit along the transverse direction is a good approximation.

Zoom In Zoom Out Reset image size

We now discuss the critical regime ω₀ → 0. Here we do not set ω₀ = 0 from the beginning, but we take a very small ω₀ and then take the large ℓ limit. The two limits are known to not commute in 1D [84]. Although the technique to obtain the 2D results from 1D one is the same for bosons and fermions, the physics is very different. Indeed, while for fermions the N chains in the Hamiltonian (31) are all critical, just with renormalised chemical potentials (33), for free bosons only the zero-mode chain is critical and all the other have a gap given by equation (87) that does not close as ω₀ → 0. This different behaviour is the origin of the logarithmic multiplicative correction to the area law for massless fermions, while massless bosons follow a strict area law, with additive logarithmic corrections. While these physical results are well known in the literature (see, e.g., [20]) we find their explanation with dimensional reduction particularly clear.

We can now sum the contributions of the various transverse modes to get the total entanglement entropy. For the zero-mode with r = 0 we take the result from the massive Klein–Gordon theory [84]. Summing up the various contributions, we have for the two-dimensional lattice complex oscillators

$$\begin{align}{S}_{n}^{2\mathrm{D}}=\frac{n+1}{3n}\enspace \mathrm{log}\enspace \ell +n\enspace \mathrm{log}\left(-\mathrm{log}\left({\omega }_{0}\ell \right)\right)+{c}_{n}^{1\mathrm{D}}\\ \quad \quad \enspace +\frac{2}{1-n}\sum _{r=1}^{N-1}\sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){{\epsilon}}_{r}}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}}\right]\right).\hfill \end{align} \tag{ 93 }$$

Here the first line is the zero gapless transverse mode and the second is the sum over all massive ones. The additive constant ${c}_{n}^{1\mathrm{D}}$ is non-universal and is not predicted by field theory; we will fix it numerically with a standard fit of the 1D system. All the chains with r > 0 give a O(1) contribution in ℓ since, for large enough ℓ, it holds $\ell \gg {\omega }_{r}^{-1}$; hence they give rise to an area-law scaling (i.e., ∝N). The panel (b) of figure 11 confirms the accuracy of the prediction (93) for the critical regime as a function of N.

4.1.1. Some generalisations

As in the fermionic case, let us mention that this technique can also be applied when imposing DBC's along the transverse direction, i.e., q_i,0 = q_i,N = p_i,0 = p_i,N = 0. Because of the breaking of translational invariance, we can simply use the Fourier sine transform

$$\begin{equation}{\mathbf{q}}_{x,y}=\sqrt{\frac{2}{N}}\sum _{r=1}^{N-1}{\tilde {q}}_{x,r}\enspace \mathrm{sin}\left(\frac{\pi ry}{N}\right),\quad {\tilde {q}}_{x,r}=\sqrt{\frac{2}{N}}\sum _{y=1}^{N-1}{\mathbf{q}}_{x,y}\enspace \mathrm{sin}\left(\frac{\pi ry}{N}\right).\end{equation} \tag{ 94 }$$

(and similarly for ${\tilde {p}}_{x,y}$). The key difference with respect to the periodic case is that the frequencies of the transverse modes are

$$\begin{equation}{\omega }_{r}^{2}={\omega }_{0}^{2}+4\enspace {\mathrm{sin}}^{2}\enspace \frac{\pi r}{2N},\quad r=1,\enspace \dots ,\enspace N-1.\end{equation} \tag{ 95 }$$

Thus, within these BC, the frequencies are all different from zero, even for ω₀ = 0. The Rényi entanglement is

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{2}{1-n}\sum _{r=1}^{N-1}\sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){{\epsilon}}_{r}}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}}\right]\right).\end{equation} \tag{ 96 }$$

We can now take the limit of large N, similarly to what done in equation (91) for ω₀ > 0, to get

$$\begin{align}\hfill {S}_{n}^{2\mathrm{D}}& =\frac{2N}{1-n}{\int }_{0}^{1}\mathrm{d}\zeta \sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}\left(\zeta \right)}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}\left(\zeta \right)}\right]\right)\hfill \\ \hfill & \quad -\enspace \frac{2}{1-n}\sum _{j=0}^{\infty }\left(n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){{\epsilon}}_{0}}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{0}}\right]\right),\hfill \end{align} \tag{ 97 }$$

where we need to subtract the contribution from the zero mode, since in equation (96) the sum starts from r = 1 rather than 0. The accuracy of equation (96) is checked by numerics in the panel (c) of figure 11 in which the agreement is perfect.

Here we compute the contributions to the entanglement entropy coming from the different U(1) symmetry sectors for the 2D lattice of oscillators. The conserved charge reduced to the subsystem is just the 2D generalisation of Q_A in equation (22). To get the 2D results for the strip geometry, we use dimensional reduction and the 1D findings of reference [44] through the CTM approach.

4.2.1. Charged moments

For the 2D lattice and for a subsystem being a periodic strip, the charged moments are obtained as a sum of the N independent chains as

$$\begin{align}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)& =\sum _{r=0}^{N-1}\sum _{j=0}^{\infty }\left[2n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){{\epsilon}}_{r}}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}+\mathrm{i}\alpha }\right]\right.\\ & \quad \left.-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}-\mathrm{i}\alpha }\right]\right],\hfill \end{align} \tag{ 98 }$$

and, taking the limit of large N

$$\begin{align}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)& =N{\int }_{0}^{1}\mathrm{d}\zeta \sum _{j=0}^{\infty }\left[2n\enspace \mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right){\epsilon}\left(\zeta \right)}\right]-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}\left(\zeta \right)+\mathrm{i}\alpha }\right]+\right.\\ & \quad \left.-\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{\epsilon}\left(\zeta \right)-\mathrm{i}\alpha }\right]\right].\hfill \end{align} \tag{ 99 }$$

Notice that $\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ is real and even in α. We plot $\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ as a function of α and N in figure 12 (left panels) together with the corresponding numerical data (for the numerical details see the appendix A). The agreement is perfect for all considered values of n, α, N, and ω₀. We find that for all α, $\mathrm{log}\enspace {Z}_{1}^{2\mathrm{D}}\left(\alpha \right)$ is a monotonously increasing function of the self frequency of the oscillators, ω₀. As a function of α, they have a single maximum at α = 0. The plots as a function of N show that the integral in equation (99) well predicts the prefactor of the area-law term of Z_n(α) in the massive case, $N\gg {\omega }_{0}^{-1}$. To obtain the data in the figure we fix ℓ = 50 which is much larger than the correlation length at all considered ω₀; hence, by cluster decomposition, they approach the double of the prediction (99) (see discussion in section 2.1 below equation (10)). We also analyse the scaling of the charged moments for DBC's along the transverse direction. The corresponding results are also displayed in the right panels of figure 12. The computation through the dimensional reduction perfectly works for every ω₀.

Zoom In Zoom Out Reset image size

In the critical regime, ω₀ → 0, the subsystem is a finite strip of longitudinal length ℓ and the resulting pattern is similar to the case encountered when α = 0. Using equation (25) and assuming that ℓ ≫ ξ, we obtain

$$\begin{align}\mathrm{log}\enspace \frac{{Z}_{n}^{2\mathrm{D}}\left(\alpha \right)}{{Z}_{n}^{2\mathrm{D}}\left(0\right)}& \simeq \frac{2}{n}\left[{\left(\frac{\alpha }{2\pi }\right)}^{2}-\frac{\vert \alpha \vert }{2\pi }\right]\enspace \mathrm{log}\ell +-\sum _{r=1}^{N-1}\sum _{j=0}^{\infty }\left[\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}+\mathrm{i}\alpha }\right]\right.\\ & \quad \left.+\mathrm{log}\left[1-{\mathrm{e}}^{-\left(2j+1\right)n{{\epsilon}}_{r}-\mathrm{i}\alpha }\right]\right],\hfill \end{align} \tag{ 100 }$$

where the second line is an additive (ℓ-independent) term that represents the contribution of the chains with r > 0. Figure 13 shows that the agreement of equation (100) with numerical data is better as ω₀ is smaller, i.e., when the approximation of a critical regime is valid. Also it works better for n closer to 1. Note that we had to consider values of ω₀ much smaller as compared to the same calculation for α = 0 to fit the numerics with the analytical prediction for the critical regime. This is likely due to the lack of a more detailed knowledge of the subleading corrections to ${Z}_{n}^{\left(1\mathrm{D}\right)}\left(\alpha \right)$ in the critical regime (which instead we have for free fermions).

Zoom In Zoom Out Reset image size

4.2.2. Symmetry resolution

In order to get the symmetry resolution it is convenient to first rewrite $\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ in the limit N → ∞ as

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)=N{\int }_{0}^{1}\mathrm{d}\zeta \mathrm{log}\left[\frac{{\theta }_{4}\left(0\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}{{\theta }_{4}\left(\frac{\alpha }{2}\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}\prod _{j=1}^{\infty }\frac{{\left(1-{\mathrm{e}}^{-\left(2j-1\right){\epsilon}\left(\zeta \right)}\right)}^{2n}}{{\left(1-{\mathrm{e}}^{-\left(2j-1\right)n{\epsilon}\left(\zeta \right)}\right)}^{2}}\right]\equiv N{f}_{n}\left(\alpha \right),\end{equation} \tag{ 101 }$$

where θ₄ is one of the Jacobi theta function

$$\begin{equation}{\theta }_{4}\left(z\vert u\right)=\sum _{k=-\infty }^{\infty }{\left(-1\right)}^{k}\enspace {u}^{{k}^{2}}\enspace {\mathrm{e}}^{2\mathrm{i}kz}.\end{equation} \tag{ 102 }$$

Then, the Fourier transform ${\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)$ is

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)={\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }\prod _{r=0}^{N-1}{Z}_{n,r}^{1\mathrm{D}}\left(\alpha \right),\end{equation} \tag{ 103 }$$

i.e., it is the convolution of the Fourier transforms ${\mathcal{Z}}_{n,r}^{1\mathrm{D}}\left(q\right)$ of ${Z}_{n,r}^{1\mathrm{D}}\left(\alpha \right)$. This formula can be easily evaluated for any finite N, even very large. In order to test its accuracy we take the Fourier transform of the numerical data for ${Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ in the previous section and compare it with equation (103). The results are shown in figure 14 where the symmetry resolved moments are plotted both against N and q for different values of ω₀. The agreement is excellent. Note that ${\mathcal{Z}}_{n}\left(q\right)$ is peaked at q = 0, which is the average charge in the subsystem.

Zoom In Zoom Out Reset image size

For large N, we can use equation (101) so that equation (103) can be rewritten as

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)\simeq {\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }\enspace {\mathrm{e}}^{N{f}_{n}\left(\alpha \right)}.\end{equation} \tag{ 104 }$$

Given that we are interested in the large N limit, the integral may be performed by saddle point method, with the only maximum of f_n(α) in α = 0, as we can see in figure 12. Therefore, the integral becomes

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)\simeq {\mathrm{e}}^{N{f}_{n}\left(0\right)}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }\enspace {\mathrm{e}}^{-\frac{N{\alpha }^{2}}{2}{\int }_{0}^{1}\mathrm{d}\zeta \frac{{\theta }_{4}^{{\prime\prime}}\left(0\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}{4{\theta }_{4}\left(0\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}}={Z}_{n}^{2\mathrm{D}}\left(0\right)\frac{{\mathrm{e}}^{-\frac{{q}^{2}}{2g\left(n\right)N}}}{\sqrt{2\pi Ng\left(n\right)}},\end{equation} \tag{ 105 }$$

where we defined

$$\begin{equation}g\left(n\right)\equiv {\int }_{0}^{1}\mathrm{d}\zeta \frac{{\theta }_{4}^{{\prime\prime}}\left(0\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}{4{\theta }_{4}\left(0\vert {\mathrm{e}}^{-n{\epsilon}\left(\zeta \right)}\right)}.\end{equation} \tag{ 106 }$$

The probability distributions given by these moments are Gaussian with mean $\overline{q}=0$ and variance that grows as $\sqrt{N}$. Unfortunately it is difficult to test equation (105) against numerical calculations because we would need rather large values of N. We instead checked that indeed equation (103) converges for large N to (105). Anyhow, we will show the corresponding plot only for the symmetry resolved entropies below.

The last step now is to use equation (103) to calculate the symmetry resolved entropies

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(q\right)=\frac{1}{1-n}\enspace \mathrm{log}\left[\frac{{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)}{{\mathcal{Z}}_{1}^{2\mathrm{D}}{\left(q\right)}^{n}}\right],\end{equation} \tag{ 107 }$$

whose limit n → 1 is the symmetry resolved von Neumann entropy

$$\begin{equation}{S}_{1}^{2\mathrm{D}}\left(q\right)\simeq -\frac{{\partial }_{n}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right){\vert }_{n=1}}{{\mathcal{Z}}_{1}^{2\mathrm{D}}\left(q\right)}+\mathrm{log}\enspace {\mathcal{Z}}_{1}^{2\mathrm{D}}\left(q\right).\end{equation} \tag{ 108 }$$

Using the previously obtained ${\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)$, we have analytic predictions valid for any N. In figure 15 we test the accuracy of this prediction for the entropy in each symmetry sector for N as large as 200. The figure clearly shows that for these relatively small values of N, the equipartition of the entanglement does not hold, even though for n = 1 data start becoming parallel to each other, suggesting a possible onset of equipartition.

Zoom In Zoom Out Reset image size

In order to understand if and how equipartition is attained at larger values of N, we work out the large N limit. In the limit N → ∞ plugging equation (105) into equation (107), we obtain

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(q\right)=\frac{1}{1-n}\enspace \mathrm{log}\frac{{Z}_{n}^{2\mathrm{D}}\left(0\right)}{{\left({Z}_{1}^{2\mathrm{D}}\left(0\right)\right)}^{n}}\enspace {\mathrm{e}}^{-\frac{{q}^{2}}{2N}\left(\frac{1}{g\left(n\right)}-\frac{n}{g\left(1\right)}\right)}\frac{{\left(2\pi Ng\left(1\right)\right)}^{n/2}}{{\left(2\pi Ng\left(n\right)\right)}^{1/2}}.\end{equation} \tag{ 109 }$$

The first ratio in equation (109) just gives the total Rényi entropy of order n, while the non-trivial dependence on n of g(n) is responsible for the breaking of the equipartition of the entanglement. After some algebra, we obtain

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(q\right)={S}_{n}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left(2\pi N\right)-\frac{1}{2\left(1-n\right)}\enspace \mathrm{log}\enspace \frac{g\left(n\right)}{g{\left(1\right)}^{n}}-\frac{{q}^{2}}{2\left(1-n\right)N}\left(\frac{1}{g\left(n\right)}-\frac{n}{g\left(1\right)}\right),\end{equation} \tag{ 110 }$$

whose limit n → 1 is

$$\begin{equation}{S}_{1}^{2\mathrm{D}}\left(q\right)={S}_{1}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left(2\pi N\right)-\frac{{q}^{2}}{2N}\frac{{g}^{\prime }\left(1\right)+g\left(1\right)}{g{\left(1\right)}^{2}}+\frac{1}{2}\left(\frac{{g}^{\prime }\left(1\right)}{g\left(1\right)}-\mathrm{log}\enspace g\left(1\right)\right).\end{equation} \tag{ 111 }$$

Hence, we have shown that the leading terms in the expansion for large N satisfy the equipartition of entanglement. The first term breaking it is at order 1/N and has an amplitude proportional to q².

Unfortunately, as already mentioned, it is difficult to test numerically the validity of equation (110) because it requires too large value of N. A posteriori, the reason of this peculiar behaviour is easily understood from equations (110) and (111): the prefactor of the equipartition breaking term multiplying q²/N is −103.485... for n = 1 and −1793.66... for n = 2, very large in both cases. Hence, we should get to values of N of order of thousands in order to see equipartition and this is not simply done numerically. What instead we can easily do is to test that for large N the analytic prediction (107) tends indeed to the predicted asymptotic behaviour (110). This is shown in the left of figure 15 where we see that very large values of N are required to recover the asymptotic behaviour, especially for large values of q and n. Hence equipartition is attained for larger and larger values of N as q and n grow, as very clear from the figure.

From the structure of the Hamiltonian in equation (31), the Rényi entropies decompose as

$$\begin{equation}{S}_{n}^{2\mathrm{D}}\left(\underset{r=0}{\overset{N-1}{\bigotimes}}{\rho }_{{k}_{y}^{\left(r\right)},\mathrm{A}}^{n}\right)=\sum _{r=0}^{N-1}{S}_{n,r}^{1\mathrm{d}}\end{equation} \tag{ 120 }$$

where ${S}_{n,r}^{1\mathrm{D}}$ is given in equation (39).

In this setting, the Fermi momentum of each transverse mode-chain is

$$\begin{equation}\mathrm{sin}\enspace {k}_{r}^{\mathrm{F}}=\sqrt{1-{\left(\frac{\mu }{2}-\frac{1}{2}\enspace \mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}}.\end{equation} \tag{ 121 }$$

Plugging this relation into equation (120) we get

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{N}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\left(2\ell \right)+N{{\Upsilon}}_{n}+\frac{1}{12}\left(1+\frac{1}{n}\right)\sum _{r=1}^{N}\mathrm{log}\left[1-{\left(\frac{\mu }{2}-\frac{1}{2}\enspace \mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}\right].\end{equation} \tag{ 122 }$$

For any given N, the sum over r can been performed using elementary trigonometric identities which lead to

$$\begin{align}{S}_{n}^{2\mathrm{D}}& =\frac{N}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\enspace \frac{\ell }{2}+\frac{1}{12}\left(1+\frac{1}{n}\right)\mathrm{log}\left[4\left({T}_{N}\left(2-\mu \right)-{\left(-1\right)}^{N}\right)\right.\\ & \quad {\times}\left.\left({T}_{N}\left(2+\mu \right)-1\right)\right]+N{{\Upsilon}}_{n},\hfill \end{align} \tag{ 123 }$$

where T_N is the Nth Chebyshev polynomial. This formula is valid for any finite N. It is useful to define a quantity analogous to equation (42), i.e.

$$\begin{equation}{A}_{N}\left(\mu \right)=\frac{1}{2N}\enspace \mathrm{log}\left[4\left({T}_{N}\left(2-\mu \right)-{\left(-1\right)}^{N}\right)\left({T}_{N}\left(2+\mu \right)-1\right)\right],\end{equation} \tag{ 124 }$$

so that we have

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{N}{6}\left(1+\frac{1}{n}\right)\left(\mathrm{log}\enspace \frac{\ell }{2}+{A}_{N}\left(\mu \right)\right)+N{{\Upsilon}}_{n},\end{equation} \tag{ 125 }$$

The function A_N(μ) is plotted as function of N in figure 16: as N increases, it approaches an asymptotic value more quickly than in the isotropic case (see figure 2). In the limit of large N, the sum in equation (122) turns into an integral

$$\begin{equation}\frac{1}{6}\sum _{r=1}^{N}\mathrm{log}\left\vert 1-{\left(\frac{\mu }{2}-\frac{1}{2}\enspace \mathrm{cos}\left(\frac{2\pi r}{N}\right)\right)}^{2}\right\vert \to \frac{N}{12\pi }{\int }_{0}^{2\pi }\mathrm{d}x\enspace \mathrm{log}\left(1-\frac{1}{4}{\left(\mu -\mathrm{cos}\left(x\right)\right)}^{2}\right).\end{equation} \tag{ 126 }$$

which, can be performed (e.g., using the Taylor expansion of the logarithm and resumming the series). The final result is (easier than the analogous in equation (46))

$$\begin{equation}N{A}_{\infty }\left(\mu \right)=\frac{N}{2}\enspace \mathrm{log}\left[\left(2+\mu +\sqrt{{\mu }^{2}+4\mu +3}\right)\left(2-\mu +\sqrt{{\mu }^{2}-4\mu +3}\right)\right].\end{equation} \tag{ 127 }$$

Hence the total entropy for large N is

$$\begin{equation}{S}_{n}^{2\mathrm{D}}=\frac{N}{6}\left(1+\frac{1}{n}\right)\mathrm{log}\enspace \frac{\ell }{2}+\frac{N}{6}\left(1+\frac{1}{n}\right){A}_{\infty }\left(\mu \right)+N{{\Upsilon}}_{n},\end{equation} \tag{ 128 }$$

which shows the expected logarithmic correction as a consequence of the fact that we are dealing with N critical chains.

Zoom In Zoom Out Reset image size

As for the isotropic setting, the same approach is straightforwardly adapted to the computation of the entanglement entropies in the case of DBC's along the transverse direction (y-axis), where the only final difference is that, in equation (120), we have to sum over N − 1 modes, rather than N.

Moreover, also in this case the same strategy applies when the total system is a finite block of L sites along the x-direction with PBC's, by replacing ℓ with $\frac{L}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi \ell }{L}$.

The results for the total entropies are checked against exact numerical calculations in figure 17, where we report the numerical data of the Rényi entropies for different values of the index n and chemical potential μ, both for infinite (panel (a)) and finite (panel (b)) system size. It is evident that the analytical results correctly describe the data, not only through the sole leading universal behaviour ∝ℓ log ℓ, but, above all, through to the subleading terms ∝ℓ, as showed in the panel (c) of the same figure. As in the isotropic setting, the subleading oscillating corrections for n ≠ 1 are described by the sum over the oscillating contributions given by each mode in equation (54).

Zoom In Zoom Out Reset image size

As already discussed in appendix B.2, the same dimensional reduction technique can further be used to compute the symmetry resolved entanglement entropies of a system in which (in the current setting) all transverse modes correspond to 1D critical chains. Let us start with the computation of the charged moments. The factorisation of equation (57) allows us to write

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)=\sum _{r=0}^{N-1}\mathrm{log}\enspace {Z}_{n,r}^{1\mathrm{D}}\left(\alpha \right),\end{equation} \tag{ 129 }$$

and, using the explicit 1D result equation (19), the sum is performed as

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq \mathrm{i}\overline{q}\alpha -\left[\frac{1}{6}\left(n-\frac{1}{n}\right)+\frac{2}{n}{\left(\frac{\alpha }{2\pi }\right)}^{2}\right]\left(\mathrm{log}\enspace \frac{\ell }{2}+{A}_{N}\left(\mu \right)\right)N+N{\Upsilon}\left(n,\alpha \right).\end{equation} \tag{ 130 }$$

The first purely imaginary term in equation (130) is the average number of particle within A, for large N explicitly given by

$$\begin{equation}\overline{q}=\frac{N\ell }{2{\pi }^{2}}{\int }_{0}^{2\pi }\mathrm{d}x\enspace \mathrm{arccos}\left(\frac{\mu -\mathrm{cos}\enspace x}{2}\right).\end{equation} \tag{ 131 }$$

It is extensive in the subsystem volume (Nℓ), as it should, and at half-filling, μ = 0, it reproduces the simple result $\overline{q}=N\ell /2$.

Keeping only the terms up to O(α²), we rewrite (130) in the compact form as:

$$\begin{equation}\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq \mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(0\right)+\mathrm{i}\overline{q}\alpha -{\alpha }^{2}\left({\mathcal{B}}_{n}\enspace \mathrm{log}\enspace \ell +{\mathcal{C}}_{n}\right)N,\end{equation} \tag{ 132 }$$

with ${\mathcal{B}}_{n}$ given in equation (64) and

$$\begin{equation}{\mathcal{C}}_{n}=\frac{1}{2{\pi }^{2}n}\left({A}_{N}\left(\mu \right)-\mathrm{log}\enspace 2\right)-\gamma \left(n\right).\end{equation} \tag{ 133 }$$

In figure 18 we report the numerical data both for the real and the imaginary part of $\mathrm{log}\enspace {Z}_{n}^{2\mathrm{D}}\left(\alpha \right)$ for different values of n and α with the theoretical prediction in equation (130) showing that the analytical result correctly describes the data as long as |α| < π. In the same figure, we also test the usual generalisation of the geometry of a torus using the finite size form with the chord length. Also in this case, the subleading oscillatory behaviour is easily obtained by the sum over contributions for each mode given by equation (65).

Zoom In Zoom Out Reset image size

We now can compute the Fourier transform ${\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)$ of the charged moments, i.e.

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)={\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}q\alpha }{Z}_{n}^{2\mathrm{D}}\left(\alpha \right)\simeq {Z}_{n}^{2\mathrm{D}}\left(0\right){\int }_{-\pi }^{\pi }\frac{\mathrm{d}\alpha }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}\left(q-\overline{q}\right)\alpha -{\alpha }^{2}{b}_{n}},\end{equation} \tag{ 134 }$$

where the coefficient of the quadratic term is

$$\begin{equation}{b}_{n}={\mathcal{B}}_{n}N\enspace \mathrm{log}\enspace \ell +{\mathcal{C}}_{n}N.\end{equation} \tag{ 135 }$$

By means of the saddle point approximation we get a Gaussian distribution function with mean $\overline{q}$ and variance growing as $\sqrt{N\enspace \mathrm{log}\enspace \ell }$, i.e.

$$\begin{equation}{\mathcal{Z}}_{n}^{2\mathrm{D}}\left(q\right)\simeq {Z}_{n}^{2\mathrm{D}}\left(0\right){\mathrm{e}}^{-\frac{{\left(q-\overline{q}\right)}^{2}}{4\left({\mathcal{B}}_{n}N\mathrm{log}\ell +{\mathcal{C}}_{n}N\right)}}\sqrt{\frac{1}{4\pi \left({\mathcal{B}}_{n}N\enspace \mathrm{log}\enspace \ell +{\mathcal{C}}_{n}N\right)}},\end{equation} \tag{ 136 }$$

whose accuracy is checked in figure 19.

Zoom In Zoom Out Reset image size

Finally, the asymptotic behaviour of the symmetry resolved entanglement, given by equation (73), is

$$\begin{align}{S}_{n}^{2\mathrm{D}}\left(q\right)& ={S}_{n}^{2\mathrm{D}}-\frac{1}{2}\enspace \mathrm{log}\left(\frac{2N}{\pi }\left(\mathrm{log}\left(\ell /2\right)+{\delta }_{n}+{A}_{\infty }\left(\mu \right)\right)\right)+\frac{\mathrm{log}\enspace n}{2\left(1-n\right)}\\ & \quad +{\left(q-\overline{q}\right)}^{2}{\pi }^{4}\frac{n}{1-n}\frac{\left(\gamma \left(1\right)-n\gamma \left(n\right)\right)}{N{\left[\mathrm{log}\left(\ell /2\right)+{\kappa }_{n}+{A}_{\infty }\left(\mu \right)\right]}^{2}}+\cdots \enspace ,\hfill \end{align} \tag{ 137 }$$

where ${S}_{n}^{2\mathrm{D}}$ is the total Rényi entropy, δ_n and κ_n are respectively defined in equations (75) and (76) (see figure 20).

Zoom In Zoom Out Reset image size