Entanglement entropy for the valence bond solid phases of two-dimensional dimerized Heisenberg antiferromagnets

TL;DR

This paper calculates entanglement entropies in two-dimensional dimerized Heisenberg antiferromagnets, revealing area law behavior and entanglement spectrum features near quantum phase transitions.

Contribution

It introduces an analytical method combining bond-operator formalism and modified spin-wave theory to compute entanglement entropies in VBS phases of 2D dimerized antiferromagnets.

Findings

Entanglement entropies obey an area law in VBS phases.

Entanglement increases but remains finite approaching the Neel-VBS transition.

Entanglement spectra are characterized for VBS ground states.

Abstract

We calculate the bipartite von Neumann and second R\'enyi entanglement entropies of the ground states of spin-1/2 dimerized Heisenberg antiferromagnets on a square lattice. Two distinct dimerization patterns are considered: columnar and staggered. In both cases, we concentrate on the valence bond solid (VBS) phase and describe such a phase with the bond-operator representation. Within this formalism, the original spin Hamiltonian is mapped into an effective interacting boson model for the triplet excitations. We study the effective Hamiltonian at the harmonic approximation and determine the spectrum of the elementary triplet excitations. We then follow an analytical procedure, which is based on a modified spin-wave theory for finite systems and was originally employed to calculate the entanglement entropies of magnetic ordered phases, and calculate the entanglement entropies of the VBS ground states. In particular, we consider one-dimensional (line) subsystems within the square lattice, a choice that allows us to consider line subsystems with sizes up to $L' = 1000$. We combine such a procedure with the results of the bond-operator formalism at the harmonic level and show that, for both dimerized Heisenberg models, the entanglement entropies of the corresponding VBS ground states obey an area law as expected for gapped phases. For both columnar-dimer and staggered-dimer models, we also show that the entanglement entropies increase but do not diverge as the dimerization decreases and the system approaches the N\'eel--VBS quantum phase transition. Finally, the entanglement spectra associated with the VBS ground states are presented.

Figures

Tables

Table 1. Table 1: Coefficients a 𝑎 a , b 𝑏 b , and c 𝑐 c obtained by fitting the von Neumann entanglement entropies 𝒮 1 subscript 𝒮 1 \mathcal{S}_{1} shown in Figs.  6 (a) and (b) with the curve ( 49 ).

			Columnar				Staggered
$J^{\prime }$		a	b	c		a	b	c
3.40		0.03	7.31e-09	0.04		0.02	1.16e-05	0.04
2.70		0.04	6.18e-09	0.07		0.03	1.07e-05	0.06
2.20		0.06	3.78e-08	0.09		0.04	5.84e-05	0.08

Table 2. Table 2: Coefficients a 𝑎 a , b 𝑏 b , and c 𝑐 c obtained by fitting the second Rényi entanglement entropies 𝒮 2 subscript 𝒮 2 \mathcal{S}_{2} shown in Figs.  6 (a) and (b) with the curve ( 49 ).

			Columnar				Staggered
$J^{\prime }$		a	b	c		a	b	c
3.40		0.01	3.03e-09	0.02		0.01	5.82e-06	0.01
2.70		0.02	8.62e-10	0.03		0.01	6.44e-06	0.02
2.20		0.03	1.39e-08	0.04		0.02	5.08e-08	0.04

Equations

$$\mathcal{S}=\mathcal{S}({\rho }_A)=-\mathrm{Tr}\left({\rho }_A\ln {\rho }_A\right),$$

$${\mathcal{S}}_{\alpha }={\mathcal{S}}_{\alpha }({\rho }_A)=\frac{1}{1-\alpha }\ln \left[\mathrm{Tr}\left({\rho }_A^{\alpha }\right)\right],$$

$${\mathcal{S}}_{\alpha }=a_{\alpha }{L}^{d-1}-\gamma .$$

$${\mathcal{S}}_{\alpha }=a_{\alpha }{L}^{d-1}+\frac{1}{2}{n}_G\ln \left(\frac{{\rho }_s}{v}{L}^{d-1}\right)+{\gamma }_{\alpha }^{ord},$$

$$H=\sum _{⟨ij⟩}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$H=J^{\prime }\sum _{i\in \mathcal{D}}{\mathbf{S}}_i^1\cdot {\mathbf{S}}_i^2+\sum _{\mu \nu }\sum _{i\in \mathcal{D}}\sum _{\tau }{\mathbf{S}}_i^{\mu }\cdot {\mathbf{S}}_{i+\tau }^{\nu }.$$

$$\mbox{\boldmath$\tau$}_{1}=2a\hat{x},\;\;\;\;\;\;\;\;\;\;\;\mbox{\boldmath$\tau$}_{2}=a\hat{y},$$

$$\mbox{\boldmath$\tau$}_{1}=a(\hat{y}-\hat{x}),\;\;\;\;\;\;\;\;\;\;\;\mbox{\boldmath$\tau$}_{2}=a(\hat{x}+\hat{y}),$$

$$H=$$

$$\displaystyle+\sum_{i\in\mathcal{D}}\Bigl{(}{\bf S}_{i}^{1}\cdot{\bf S}_{i+2}^{1}+{\bf S}_{i}^{2}\cdot{\bf S}_{i+2}^{2}+{\bf S}_{i}^{2}\cdot{\bf S}_{i+1}^{1}\Bigr{)},$$

$$H=J^{\prime }\sum _{i\in \mathcal{D}}{\mathbf{S}}_i^1\cdot {\mathbf{S}}_i^2$$

$$\displaystyle+{\bf S}_{i}^{2}\cdot{\bf S}_{i+2}^{1}+{\bf S}_{i}^{2}\cdot{\bf S}_{i+2-1}^{1}\Bigr{)}$$

$$|s⟩$$

$$|t_y⟩$$

$$|s⟩=s^{†}|0⟩,|t_{\alpha }⟩=t_{\alpha }^{†}|0⟩,\alpha =x,y,z.$$

$$s^{†}s+\sum _{\alpha }{t}_{\alpha }^{†}{t}_{\alpha }=1.$$

$$S_{\alpha }^{1,2}=\pm \left(s^{†}{t}_{\alpha }+t_{\alpha }^{†}s\mp i{ϵ}_{\alpha \beta \gamma }{t}_{\beta }^{†}{t}_{\gamma }\right).$$

$$H=H_0+H_2+H_3+H_4.$$

$$H_0=$$

$$H_2=$$

$$H_3=$$

$$H_4=$$

$$g_2(\tau )$$

$$g_3(\tau )$$

$$g_4(\tau )$$

$$g_2(\tau )$$

$$g_3(\tau )$$

$$g_4(\tau )$$

$$-\mu \sum _i\left(s_i^{†}{s}_i+t_{i\alpha }^{†}{t}_{i\alpha }-1\right).$$

$$s_i^{†}=s_i=⟨s_i^{†}⟩=⟨s_i⟩\to \sqrt{N_0}$$

$$t_{i\alpha }^{†}=\frac{1}{\sqrt{N^{\prime }}}\sum _{\mathbf{k}\in \mathrm{BZ}}{e}^{-i\mathbf{k}\cdot {\mathbf{R}}_i}{t}_{\mathbf{k}\alpha }^{†},$$

$$H_0=$$

$$H_2=$$

$$H_3=$$

$$H_4=$$

$$A_{\mathbf{k}}$$

$$B_{\mathbf{k}}$$

Full text

I Introduction

In the last few years, bipartite entanglement entropies have been used to characterise many-body quantum systems rmp08 ; rmp-area-law ; grover13 ; wen-book ; review-nicolas . Such quantities could offer additional information that, in principle, could not be obtained from correlation functions. In particular, bipartite entanglement entropies have been employed to study interacting spin systems song11 ; kallin11 ; lou11 ; grover11 ; luitz14 ; wessel14 ; luitz15 ; alet15 ; melko17 .

The bipartite entanglement entropy for pure states is defined as follows: Consider, for instance, the ground state $|\Psi\rangle$ of a system $S$, a subsystem $A$ (arbitrary size and shape) and its complementary $\bar{A}$ such that $S=A\cup\bar{A}$. The entanglement entropy is defined, for instance, as the von Neumann entropy grover13 ; review-nicolas ,

[TABLE]

where $\rho_{A}={\rm Tr}_{\bar{A}}|\Psi\rangle\langle\Psi|$ is the reduced density matrix of the subsystem $A$. Alternatively, the entanglement entropy is defined as the Rényi entropy

[TABLE]

where the index $\alpha>0$ acts as a weight for the probabilities. In the limit $\alpha\rightarrow 1$, the Rényi entropy (2) reduces to the von Neumann entropy (1).

An important issue is the scaling of the $T=0$ entanglement entropy when the system $S$ is in the thermodynamic limit. For instance, for gapped systems in spatial dimensions $d>1$ described by local Hamiltonians, it is found that the entanglement entropy (2) assumes the general form wen-book

[TABLE]

Here the leading term, that depends on the size $L^{d-1}$ of the boundary between the subsystem $A$ and its complementary $\bar{A}$, is the so-called area law rmp-area-law . The coefficient $a_{\alpha}$ is a non-universal constant. The second term $\gamma>0$ is a universal constant known as the topological entanglement entropy kitaev06 . Such a quantity indicates whether the ground state has non-trivial topological order wen-book . An example of a topologically ordered phase (a phase that cannot be characterised by a local order parameter) is a gapped Z2 spin liquid with $\gamma=\ln 2$ melko17 .

The entanglement entropy has also been used to characterize gapless systems. In particular, it was found that the entanglement entropies of the ground state of the spin-$1/2$ antiferromagnet (AFM) Heisenberg model with nearest-neighbor interactions on a square lattice obey an area law with additive logarithmic corrections song11 ; kallin11 . Later, Metliski and Grover grover11 analytically calculated the Rényi entanglement entropy (2) of a phase that spontaneously breaks a continuous symmetry and, for a corner-free subsystem, it was showed that

[TABLE]

where $\rho_{s}$ is the spin stiffness, $v$ is the velocity of the $n_{G}$ Goldstone modes, and $\gamma_{\alpha}^{ord}$ is a non-universal constant. Interestingly, the coefficient of the additional logarithmic correction to the area law is equal to one-half the number of Goldstone modes $n_{G}$. Recently, the entanglement entropy of the Néel phase of spin-$1/2$ square lattice Heisenberg AFMs has been calculated within a modified spin-wave theory for finite systems luitz15 ; alet15 and the obtained results are in good agreement with Eq. (4).

In this paper, we calculate the von Neumann (1) and second ($\alpha=2$) Rényi (2) entanglement entropies of the ground states of spin-$1/2$ dimerized Heisenberg antiferromagnets on a square lattice focusing on (quantum paramagnet) valence bond solid (VBS) phases. We describe the VBS phases within the bond-operator representation bond-op , a formalism that allow us to map the original spin Hamiltonian into an effective interacting boson model for the triplet (triplon) excitations. We then consider such an effective boson Hamiltonian at the harmonic approximation and determine the bipartite entanglement entropies via a procedure similar to the modified spin-wave theory for finite systems song11 ; luitz15 ; alet15 employed for symmetry broken phases. In particular, we consider one-dimensional (line) subsystems of size $L^{\prime}$ within the square lattice and analytically calculate the bipartite entanglement entropies. Such a procedure also allow us to derive the corresponding entanglement spectra.

I.1 Overview of the results

In the first part of the paper (Secs. II–IV), we study two square lattice dimerized Heisenberg AFMs with columnar [Fig. 1(a)] and staggered [Fig. 1(b)] dimerization patterns. We calculate the dispersion relation of the elementary (triplon) excitations of the VBS phases of the two dimer-models (Fig. 4) within the bond-operator formalism at the (mean-field) harmonic level. The triplon energy gaps (Fig. 5) and the quantum critical points where the Néel–VBS quantum phase transition (QPT) takes place (Sec. IV.1) are determined.

In the second part of the paper (Sec. V), we calculate the bipartite von Neumann (1) and second Rényi (2) entanglement entropies of the VBS ground states of both dimer-models. It is shown here that the combination of the bond-operator results at the harmonic approximation with an approach similar to the one used in Refs. song11 ; luitz15 ; alet15 for magnetic ordered phases provides the area law behaviour for the entanglement entropies, a behaviour expected for gapped phases. This is indeed our main result. Importantly, our results are derived for line (chain) subsystems $A$, a choice that allows us to determine the entanglement entropies for very large subsystem sizes. Furthermore, we show that the entanglement entropies seem to not diverge as the system approaches the Néel–VBS quantum phase transition, but only reaches a maximum value (Fig. 7). Finally, we show that the corresponding entanglement spectra for the VBS phases are gapped even when close to the Néel–VBS quantum critical point (Fig. 8).

I.2 Outline

Our paper is organized as follows: In Sec. II, we introduce the square lattice columnar and staggered dimerized Heisenberg AFMs considered in our study. In Sec. III, the bond-operator representation bond-op for spin operators is briefly summarized and the effective interacting boson models corresponding to the two dimerized Heisenberg antiferromagnets are derived. The analysis of the effective boson models within the harmonic approximation, in particular, the determination of the triplet (triplon) excitation spectra, is presented in Sec. IV. In Sec. V, we briefly review the procedure employed in Refs. song11 ; luitz15 ; alet15 for the calculation of the entanglement entropies of magnetic ordered phases, determine the bipartite entanglement entropies for the VBS phases of the columnar-dimer and staggered-dimer models, and discuss the corresponding entanglement spectra. Finally, in Sec. VI, we provide a brief summary of our main findings. A short discussion about the classical dimerized Heisenberg AFMs and some technical details of the scheme adopted for the calculation of the entanglement entropies are included in the three Appendices.

II Square lattice dimerized antiferromagnets

Let us consider the dimerized AFM Heisenberg model on a square lattice:

[TABLE]

where ${\bf S}_{i}$ is a spin-$1/2$ operator at site $i$ and the nearest-neighbor exchange couplings $J_{ij}=J>0$ and $J^{\prime}>0$ are arranged according to the columnar and staggered patterns respectively shown in Figs. 1(a) and (b). Hereafter, we set $J=1$.

The Hamiltonian (5) is an interesting model system to study quantum phase transitions vojta-qpt since its ground state depends on the (intra-dimer) exchange coupling $J^{\prime}$: for $J^{\prime}\sim 1$, the ground state has semiclassical Néel magnetic long-range order while, for $J^{\prime}\gg 1$, a quantum paramagnetic (disordered) phase sets in and the ground state is given by a VBS of short singlets as illustrated in Figs. 1(c) and (d). The Néel–VBS quantum phase transition takes place at the critical couplings $J^{\prime}_{c}=1.9096(2)$ (columnar) wenzel09 and $J^{\prime}_{c}=2.5196(2)$ (staggered) wenzel08 . According to the quantum-to-classical mapping, this QPT should be in the same universality class of the classical Heisenberg model in $(2+1)$-dimensions, the so-called $O(3)$ universality class doretto11 . However, quantum Monte Carlo (QMC) results wenzel08 indicated that such a scenario applies only to the columnar-dimer model [Fig. 1(a)]: For the staggered-dimer model [Fig. 1(b)], it was found that the critical exponents showed small deviations from the ones of the $O(3)$ universality class, an interesting feature that motivated further investigations doretto11 ; jiang09 ; jiang12 ; yasuda13 ; sandvik18 . It was then proposed doretto11 that the critical exponents of the staggered-dimer model are indeed the ones of the $O(3)$ universality class, but with anomalously large corrections to scaling related to cubic triplet interactions, see Eq. (23) below. Such a scenario was later confirmed by QMC calculations jiang12 ; yasuda13 ; sandvik18 . This interesting feature of the critical behaviour of the columnar-dimer and staggered-dimer models found in finite-size QMC simulations is also a motivation for our study. In the following, we will focus on the region $J^{\prime}>1$ of the phase diagram of the model (5).

In order to describe the VBS phases of the columnar-dimer and staggered-dimer Heisenberg AFMs, it is useful to rewrite the Hamiltonian (5) in terms of the underline lattices defined by the strong couplings $J^{\prime}$:

[TABLE]

Here $i$ indicates a site of the dimerized lattice $\mathcal{D}$, which has two spins per unit cell labeled by the indices $\mu$ and $\nu=1,2$, see Figs. 1(a) and (b). The index $\tau=1,2$ corresponds to the dimer nearest-neighbor vectors \mbox{\boldmath\tau}_{i}: for the columnar-dimer model [Fig. 1(a)], we have

[TABLE]

whereas, for the staggered-dimer model [Fig. 1(b)],

[TABLE]

with $a$ being the lattice spacing of the original square lattice (in the following we set $a=1$). In terms of the nearest-neighbor vectors (7) and (8), the Hamiltonian (6) can be explicitly written as

[TABLE]

for the columnar-dimer model, and

[TABLE]

for the staggered-dimer model.

III Bond-operator representation

The VBS phases of the dimerized Heisenberg AFMs (5) can be described within the bond-operator representation for spin operators bond-op . In the following, we briefly summarize this formalism.

We start by considering two spins-$1/2$: ${\bf S}^{1}$ and ${\bf S}^{2}$. The Hilbert space of the system is made out of a singlet state $|s\rangle$ and three triplet states $|t_{\alpha}\rangle$:

[TABLE]

It is possible to define a set of boson operators $s^{\dagger}$ and $t^{\dagger}_{\alpha}$ with $\alpha=x$, $y$, $z$ that creates the states (11) out of a fictitious vacuum $|0\rangle$, namely,

[TABLE]

The unphysical states of the enlarged Hilbert space are removed via the introduction of the constraint

[TABLE]

We then calculate the matrix elements of each component of the two spins operators within the basis $|s\rangle$ and $|t_{\alpha}\rangle$, i.e., we determine $\langle s|S^{1}_{\alpha}|s\rangle$, $\langle s|S^{1}_{\alpha}|t_{\beta}\rangle$, $\ldots$, and therefore, based on the obtained results, conclude that the components of the spin operators ${\bf S}^{1}$ and ${\bf S}^{2}$ can be written in terms of the boson operators $s^{\dagger}$ and $t^{\dagger}_{\alpha}$ as

[TABLE]

Here the indices $\alpha,\beta,\gamma=x,y,z$, $\epsilon_{\alpha\beta\gamma}$ is the completely antisymmetric tensor with $\epsilon_{xyz}=1$, and the summation convention over repeated indices is implied.

The bond-operator representation (14) can be generalized to the lattice case, allowing us to express the dimerized Heisenberg models (9) and (10) in terms of the boson operators $s_{i}^{\dagger}$ and $t^{\dagger}_{i\alpha}$.

III.1 Effective boson models

Substituting Eq. (14) generalized to the lattice case into the Hamiltonian (9) of the columnar-dimer model, we find that the Hamiltonian can be written as

[TABLE]

Here the $H_{n}$ terms have $n$ triplet operators and are given by

[TABLE]

with summation convention over repeated indices implied. The $g_{i}(\tau)$ functions are defined as

[TABLE]

with the dimer nearest-neighbor vectors \mbox{\boldmath\tau}_{i} given by Eq. (7). A similar expression is found for the Hamiltonian (10) of the staggered-dimer model, but now the $g_{i}(\tau)$ functions read

[TABLE]

with the \mbox{\boldmath\tau}_{i} vectors defined as in Eq. (8). One should note that only the last term of the Hamiltonian (9) contributes to the cubic term $H_{3}$ whereas, for the staggered-dimer model, all three nonlocal terms of the Hamiltonian (10) provide a nonvanishing contribution for $H_{3}$, a feature that can be understood on symmetry grounds, see Sec. II.C from Ref doretto11 . Finally, the constraint (13) is taken into account on average via a Lagrange multiplier $\mu$, i.e., we add the following term to the Hamiltonian (15)

[TABLE]

Within the bond-operator formalism, the VBS ground states illustrated in Figs. 1(c) and (d) can be viewed as a condensate of the singlets $s_{i}$. We then set

[TABLE]

in the Hamiltonian (15) and arrive at an effective boson Hamiltonian only in terms of the triplet $t_{i\alpha}$ boson operators. Both constants $N_{0}$ and $\mu$ will be self-consistently calculated for a fixed value of the exchange coupling $J^{\prime}$.

Performing a Fourier transform,

[TABLE]

where ${\bf R}_{i}$ is a vector of the dimerized lattice, $N^{\prime}=N/2$ is the number of dimers ($N$ is the number of sites of the original square lattice), and the momentum sum runs over the corresponding dimerized first Brillouin zones, see Fig. 2, one shows that the four terms (16) of the Hamiltonian (15) assume the form

[TABLE]

For the columnar-dimer model, the coefficients $A_{\bf k}$, $B_{\bf k}$, $\xi_{\bf k}$, and $\gamma_{\bf k}$ are given by

[TABLE]

while, for the staggered-dimer model, we have

[TABLE]

It is important to mention that the bond-operator approach to VBS phases is quite similar to the Holstein–Primakoff one to magnetic ordered phases, but while the latter considers fluctuations (spin–waves) above a semiclassical magnetic ordered state, the former describes excitations above a quantum paramagnetic state. Such similarity will be useful in the calculation of the entanglement entropies below. Further comparisons between the two approaches can be found in Sec. II.B from Ref. doretto12 .

The procedure discussed above for the description of a VBS phase within the bond-operator formalism follows the lines of Refs. bond-op ; doretto12 . Such a scheme is slightly distinct from the previous bond-operator study doretto11 of the dimerized Heisenberg models (5), which is based on the procedure discussed in Ref. kotov98 : in this case, it is assumed that the boson operators $t^{\dagger}_{i\alpha}$ create triplet excitations out of a singlet background $|\Psi_{0}\rangle=\prod_{i}s^{\dagger}_{i}|0\rangle$; the equivalent of Eq. (19) reads $s_{i}^{\dagger}=s_{i}=1$ and the constraint (13) becomes an inequality, $\sum_{\alpha}t^{\dagger}_{i\alpha}t_{i\alpha}\leq 1$, which is implemented via an on-site triplet–triplet repulsion term added to the Hamiltonian. For both dimer-models at the harmonic approximation, it is found that the Néel–VBS QPT takes place at the critical coupling $J^{\prime}_{c}=3$. As shown below, the procedure implement in our work provides better results for $J^{\prime}_{c}$ at the (lowest-order) harmonic approximation.

IV Harmonic approximation

In this section, we consider the effective boson model (15) in the lowest-order approximation, the so-called harmonic approximation. In this case, we keep the terms of the Hamiltonian (15) up to quadratic order in the triplet boson operators $t_{{\bf k}\alpha}$, namely,

[TABLE]

Since the Hamiltonian (27) is quadratic in the triplet operators $t_{{\bf k}\alpha}$, it can be diagonalized via a Bogoliubov transformation

[TABLE]

It is then easy to show that the Hamiltonian (27) assumes the form

[TABLE]

where

[TABLE]

is the ground state energy,

[TABLE]

is the energy of the triplet (triplon) excitations above the VBS ground state, and the coefficients $u_{\bf k}$ and $v_{\bf k}$ of the Bogoliubov transformation (28) are given by

[TABLE]

Finally, we would like to quote some triplet-triplet ground state expectation values that will be useful in the determination of the entanglement entropy, see Sec. V.1 below. With the aid of Eq. (28) and using the fact that the ground state of the Hamiltonian (29) is the vacuum for the boson operators $b$, one easily shows that

[TABLE]

IV.1 Self-consistent equations

The ground state energy (30) and the triplon excitation spectrum (31) are fully determined once we calculate the constants $N_{0}$ and $\mu$ for a fixed value of the exchange coupling $J^{\prime}$. By minimizing the ground state energy (30) with respect to $\mu$ and $N_{0}$, we find a system of self-consistent equations,

[TABLE]

which are numerically solved.

The numerical solutions of Eq. (34), i.e., the behaviour of $N_{0}$ and $\mu$ in terms of $J^{\prime}$, are respectively shown in Figs. 3(a) and (b). As expected, see discussion in Sec. II, we find solutions for the system of self-consistent equations (34) only for $J^{\prime}\geq 1.70$ (columnar-dimer) and $J^{\prime}\geq 2.00$ (staggered-dimer). One sees that (i) $N_{0}\rightarrow 1$ as the coupling $J^{\prime}$ increases (system deep in the VBS phase) and (ii) $N_{0}$ decreases as $J^{\prime}$ approaches the Néel–VBS QPT.

Figure 4 shows the triplon excitation spectra (31) of the columnar-dimer [Fig. 4(a)] and the staggered-dimer [Fig. 4(b)] models for three different values of $J^{\prime}$. One sees that the triplon excitation spectrum is gapped for both models and that the triplon energy gap $\Delta$ decreases as $J^{\prime}$ approaches the Néel–VBS QPT (see details below). For the columnar-dimer model, the triplon gap $\Delta$ is located at the $Y$ point [see Fig. 2(a)] while, for the staggered-dimer model, the triplon gap is located at the centre of the first Brillouin zone, the $\Gamma$ point [see Fig. 2(b)]. For both dimer-models, it is possible to show that the momentum associated with the triplon gap $\Delta$ is equal to the ordering wave vector ${\bf Q}$ of the corresponding Néel magnetic long-range ordered phase that sets in for $J^{\prime}<J^{\prime}_{c}$, see Appendix A for details.

The behaviour of the triplon gaps $\Delta$ as a function of $J^{\prime}$ are displayed in Fig. 5. Again, one notices that the triplon gaps $\Delta$ close as the systems reach the Néel–VBS quantum critical points. In order to estimate the critical coupling $J^{\prime}_{c}$, we follow the lines of Ref. doretto12 , i.e., we assume a continuous Néel–VBS QPT, fit the data with the curve

[TABLE]

and then consider the condition $\Delta=0$. Following such a procedure, within the harmonic approximation, the critical couplings are $J^{\prime}_{c}=1.61$ (columnar) and $J^{\prime}_{c}=1.93$ (staggered), which are in quite reasonable agreement with the ones determined via QMC calculations, namely, $J^{\prime}_{c}=1.9096(2)$ (columnar) wenzel09 and $J^{\prime}_{c}=2.5196(2)$ (staggered) wenzel08 . Such an agreement is expected due to the small number of triplets $t$ in the VBS ground state, $1-N_{0}$, a quantity that could be taken as a control parameter within the bond-operator formalism. More accurate results for $J^{\prime}_{c}$ can be obtained within the bond-operator formalism by perturbatively including the cubic (23) and quartic (24) terms as done, e.g., in Ref. doretto12 . Finally, one should also mention that the critical couplings found here are in better agreement with the QMC simulations than the ones ($J_{c}^{\prime}=3$) obtained in the previous bond-operator study doretto11 at the same approximation level.

V Entanglement entropies

In this section, we calculate the von Neumann (1) and the second ($\alpha=2$) Rényi (2) entanglement entropies for the VBS ground states of both columnar-dimer and staggered-dimer models. Both quantities are good measures of entanglement, but the second Rényi entanglement entropy is easier to numerically determined grover13 . For instance, Helmes and Wessel calculated the second Rényi entanglement entropy of a two-dimensional bilayer Heisenberg AFM wessel14 based on a QMC procedure introduced in Ref. humeniuk12 . Since we would like to compare our analytical results with future numerical ones, it is interesting to determined both the the von Neumann (1) and the second ($\alpha=2$) Rényi (2) entanglement entropies within our scheme.

Due to the similarities between the description of VBS phases within the bond-operator formalism and the description of magnetic ordered phases within spin–wave theory (see Sec. III.1), we follow the lines of Refs. song11 ; luitz15 ; alet15 , where bipartite entanglement entropies for the Néel phase of two-dimensional Heisenberg AFMs are determined via a modified spin–wave theory for finite systems. For completeness, in the following we briefly outline such a scheme which is indeed based on Refs. peschel01 ; peschel03 ; peschel09 ; barthel06 ; frerot15 .

Let us consider a $d$-dimensional system $S$ described by a generic quadratic Hamiltonian frerot15

[TABLE]

where $a_{n}$ is a boson operator associated with the site $n$ of a $d$-dimensional hypercubic lattice with $N$ sites and $A_{nm}$ and $B_{nm}$ are $N\times N$ matrices. We divide the system $S$ into a subsystem $A$ with $N_{A}<N$ sites and its complementary $\bar{A}$ such that $S=A\cup\bar{A}$, see Sec. I. It is possible to show, e.g., with the aid of coherent states barthel06 ; peschel01 , that the reduced density matrix $\rho_{A}$ of the subsystem $A$ assumes the form peschel01 ; peschel03 ; peschel09 ; barthel06 ; frerot15

[TABLE]

where $\mathcal{K}$ is a normalization constant and $\mathcal{H}_{E}$ is the so-called entanglement Hamiltonian frerot15

[TABLE]

with $\mathcal{A}_{ij}$ and $\mathcal{B}_{ij}$ being $N_{A}\times N_{A}$ matrices. Notice that both Hamiltonians $H$ and $\mathcal{H}_{E}$ have the same quadratic form and the latter is restricted to the sites $i$ and $j$ of the subsystem $A$. Due to this similarity, the Hamiltonians $H$ and $\mathcal{H}_{E}$ can then be diagonalized by the same Bogoliubov transformation. In particular, we have

[TABLE]

where the energies $\epsilon_{k}$ define the entanglement spectrum review-nicolas ; haldane-es and the boson operators $a_{i}$ and $b_{k}$ are related by a Bogoliubov transformation, see Eq. (55). Since the reduced density matrix $\rho_{A}$ has the form (37), the von Neumann entanglement entropy (1) is simply given by the expression of the thermal entropy, i.e.,

[TABLE]

where $n_{k}=1/[\exp(\epsilon_{k})-1]$ is the occupation of the $k$ mode. Therefore, once the entanglement spectrum $\epsilon_{k}$ is known, the bipartite von Neumann entanglement entropy (1) is determined. Similar considerations hold for fermionic systems described by Hamiltonians of the form (36) barthel06 ; peschel01 ; cheong04 . For a more general expression for the entanglement Hamiltonian, we refer the reader, e.g., to Ref. calabrese18 , where an approximate entanglement Hamiltonian for a general lattice model is derived based on a lattice version of the so-called Bisognano-Wichmann theorem.

Instead of performing the partial trace described above to calculate the matrix elements $\mathcal{A}_{ij}$ and $\mathcal{B}_{ij}$ and then find the entanglement entropy (40), we can alternatively determine $\mathcal{S}$ from single-particle Green’s functions associated with the Hamiltonian (36) barthel06 ; peschel01 ; peschel03 ; peschel09 . Indeed, the entanglement spectrum $\epsilon_{k}$ is related to the eigenvalues of the so-called correlation matrix $C$, which is defined as barthel06 ; peschel01 ; peschel03 ; peschel09

[TABLE]

Here $i$, $j$, and $s$ refer to sites of the subsystem $A$ and $f_{ij}$ and $g_{ij}$ are single-particle Green’s functions,

[TABLE]

As discussed in details in Appendix B, one shows that the $N_{A}$ eigenvalues $\mu^{2}_{k}$ of the correlation matrix $C$ can be written in terms of the entanglement spectrum $\epsilon_{k}$ as

[TABLE]

Substituting Eq. (43) into the expression (40), one shows that the von Neumann entanglement entropy (1) reads

[TABLE]

Similarly, one finds that the Rényi entanglement entropies (2) assume the form song11 ; luitz15 ; alet15

[TABLE]

Therefore, the bipartite entanglement entropies are completely determined, once the eigenvalues $\mu^{2}_{k}$ of the correlation matrix $C$ are known.

Within linear spin–wave theory, the effective boson model that describes the Néel phase of an Heisenberg AFM has the same form as the Hamiltonian (36). Due to such similarity, the procedure discussed above was employed to calculated bipartite entanglement entropies of magnetic ordered phases of Heisenberg AFMs song11 ; luitz15 ; alet15 . In particular, the single–particle Green’s functions (42) can be easily calculated within linear spin–wave theory. Notice that the same considerations hold for the description of the VBS phases of the columnar-dimer and staggered-dimer models within the bond-operator formalism at the harmonic approximation, see Eq. (27). Therefore, in the next section, we apply the scheme described above for the VBS phases of the dimerized Heisenberg models (9) and (10).

V.1 Entanglement entropies of the VBS phases

To determine the bipartite entanglement entropies for the VBS phases of the columnar-dimer and the staggered-dimer models, we considerer a line subsystem $A$, i.e., an one-dimensional spin chain with size $L^{\prime}=2(N_{A}-1)$, as illustrated in Fig. 1(c). Such a partition is quite interesting, since it allows us to reach very large system sizes luitz15 in addition to analytically determine the entanglement entropies alet15 . Indeed, a line subsystem has been employed to study interacting spin systems luitz14 ; luitz15 ; alet15 and, in particular, it provides luitz15 a prefactor for the logarithmic term in Eq. (4) in good agreement with the analytical results of Metliski and Grover grover11 .

The matrix elements (41) of the correlation matrix $C$ are easily calculated. From Eqs. (33) and (42), one shows that

[TABLE]

where ${\bf R}_{i}=2i\hat{x}$, with $i=1,2,\ldots,N_{A}$, is a vector of the (underline) dimerized lattice of the line subsystem $A$, the coefficients $A_{\bf k}$ and $B_{\bf k}$ are given by Eqs. (25) and (26) respectively for the columnar-dimer and staggered-dimer models, $\Omega_{\bf k}$ is the triplon excitation energy (31), and the momentum sums run over the corresponding first Brillouin zones. Notice that the correlation matrix $C$ is completely determined by the coefficients $A_{\bf k}$ and $B_{\bf k}$ of the effective boson model (27).

In principle, the eigenvalues $\mu^{2}_{m}$ of the correlation matrix $C$ are numerically calculated, see, e.g., Refs. song11 ; alet15 . However, for a one-dimensional (line) subsystem $A$, the eigenvalues of the correlation matrix can be analytically determined luitz15 ; alet15 since the correlation matrix $C$ is indeed a circulant matrix gray : In this case, the eigenvalues $\mu^{2}_{m}$ are given by the Fourier transform of the first line of the correlation matrix $C$. For both dimer models, one finds (see appendix C for details)

[TABLE]

where the index $m=1,2,\cdots,N_{A}$ is related to the momentum $k_{x}$ parallel to the system-subsystem boundary,

[TABLE]

with $N^{\prime}=N_{A}N_{y}$ and $N_{A}=(L^{\prime}+2)/2$. Therefore, for a line subsystem $A$, the $N_{A}$ eigenvalues $\mu^{2}_{m}$ of the correlation matrix $C$ can be easily expressed in terms of the coefficients $A_{{\bf k}}$ and $B_{{\bf k}}$ of the effective boson model (27). Once the sum over the momentum component $k_{y}$ is evaluated by changing it to an integral, the entanglement entropies follow from Eqs. (44) and (45).

From Eqs. (47) and (48), one clearly sees how the finite-size nature of the subsystem $A$ is included in the calculation of the entanglement entropy within our approach. Recall that the coefficients $A_{{\bf k}}$ and $B_{{\bf k}}$ are obtained within the bond-operator method, a formalism suitable to describe the VBS phases in the thermodynamic limit. Similar considerations hold for Refs. song11 ; luitz15 ; alet15 , where the entanglement entropies for magnetic ordered phases are calculated via a modified spin–wave theory. One should note that in Refs. song11 ; luitz15 ; alet15 , however, an additional information about the subsystem size is encoded in a staggered magnetic field $h$ that is added to restore the spin rotational symmetry of the finite-size subsystem $A$. The value of $h$ is determined by imposing that, at each site $i$ of the lattice, the $z$-component of the spin operator $\langle S^{z}_{i}\rangle=0$ [see also Sec. VII from Ref. frerot15 for an alternative procedure]. As discussed in Ref. alet15 , this staggered magnetic field $h$ is an important ingredient to find the prefactor of the logarithmic correction to the area law proportional to the number of Goldstone modes [second term of Eq. (4)]. Here, for the VBS phases, such an additional magnetic field is not necessary, since these phases preserve the spin rotational symmetry.

Figure 6 shows the von Neumann $\mathcal{S}_{1}$ and second Rényi $\mathcal{S}_{2}$ bipartite entanglement entropies in terms of the subsystem size $L^{\prime}$ for the VBS ground states of the columnar-dimer [Fig. 6(a)] and the staggered-dimer [Fig. 6(b)] models. We consider one-dimensional subsystems $A$ with sizes up to $L^{\prime}=1000$ and show the results for three different values of the exchange coupling $J^{\prime}$. Notice that $\mathcal{S}_{1}$ is larger than $\mathcal{S}_{2}$ for the same value of $J^{\prime}$, a feature that has been found for the magnetic ordered phase of Heisenberg AFMs song11 ; luitz15 ; alet15 . Moreover, both entropies $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are dominated by an area law as expected for two-dimensional gapped phases rmp-area-law ; wen-book . Indeed, we fit the data shown in Figs. 6(a) and (b) with the curve

[TABLE]

and, for the three values of the exchange coupling $J^{\prime}$, we find that $b<10^{-5}$, see tables 1 and 2 for details. Finally, one sees that the prefactor $b$ of the logarithmic term is smaller for the columnar-dimer model than for the staggered-dimer one for both entanglement entropies $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$.

Although both bipartite entanglement entropies $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ increase as $J^{\prime}$ decreases, it seems that they do not diverge as $J^{\prime}$ approaches the Néel–VBS quantum critical point. We illustrated such a behaviour in Fig. 7, where it is shown the von Neumann entanglement entropy $\mathcal{S}_{1}$ as a function of the exchange coupling $J^{\prime}$ for a subsystem $A$ with size $L^{\prime}=400$. For both columnar-dimer and staggered-dimer models, one sees that $\mathcal{S}_{1}$ has the same qualitatively behaviour, although its is larger for the columnar-dimer model than for the staggered-dimer one. For both dimer models, $\mathcal{S}_{1}$ reaches a maximum value at the smallest exchange coupling $J^{\prime}$ determined via the numerical solutions of the self-consistent equations (34) [$J^{\prime}=1.70$ (columnar-dimer) and $J^{\prime}=2.00$ (staggered-dimer), see Sec. IV.1], a feature that indicates a possible absence of divergence at the quantum critical point. One should mention that such an absence of divergence of the entanglement entropy at criticality was previously observed in the Néel–VBS QPT of a two-dimensional bilayer Heisenberg AFM wessel14 and in the superfluid-Mott insulator QPT of a two-dimensional Bose-Hubbard model frerot16 .

Finally, one also sees in Fig. 7 that, for larger values of $J^{\prime}$, the entanglement entropy $\mathcal{S}_{1}$ slowly decreases. Indeed, one finds, e.g., for the columnar-dimer model, that $\mathcal{S}_{1}=0.57$, $0.27$, and $0.16$ respectively for $J^{\prime}=20.0$, $30.0$, and $40.0$. Such a behaviour is expected since, as the (intra-dimer) exchange coupling $J^{\prime}$ increases, the dimers get more and more isolated and, therefore, $\mathcal{S}_{1}$ should vanish in this limit.

V.2 Entanglement spectra for the VBS phases

In addition to the bipartite entanglement entropies $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, the procedure employed in our work allow us to easily calculate the entanglement spectrum $\epsilon_{k_{x}}$ as defined in Eq. (39). Notice that once the eigenvalues $\mu^{2}_{m}=\mu^{2}_{k_{x}}$ of the correlation matrix $C$ are known, the entanglement spectrum $\epsilon_{k_{x}}$ follows from Eq. (43).

In Fig. 8, we show the entanglement spectra of the columnar-dimer [Fig. 8(a)] and staggered-dimer [Fig. 8(b)] models for four different values of the exchange coupling $J^{\prime}$ (a line subsystem $A$ with size $L^{\prime}=1000$ is considered). One sees that the entanglement spectra of the two dimer models are qualitatively similar, although the bandwidth is larger for the staggered-dimer model than for the columnar-dimer one at the same value of the exchange coupling $J^{\prime}$. For both dimer models, the energy gap of the entanglement spectra is at $k_{x}=0$ and the gap decreases as $J^{\prime}$ approaches the critical coupling $J^{\prime}_{c}$. Differently from the triplon spectrum $\Omega_{\bf k}$, the gap of the entanglement spectrum does not close as the system approaches the Néel–VBS quantum critical point: we find $\epsilon_{k_{x}=0}=1.77$ ($J^{\prime}=1.70$, columnar-dimer) and $\epsilon_{k_{x}=0}=2.47$ ($J^{\prime}=2.00$, staggered-dimer), compare with Fig. 5. Such a feature is in contrast with the behaviour of the two-dimensional Bose-Hubbard model which also displays a QPT between a gapped (Mott-insulator) and a gapless (superfluid) phases frerot16 : it was found that the gap of the entanglement spectrum closes at the superfluid–Mott insulator QPT driven by the ratio $t/U$ between the tunneling amplitude $t$ and the on-site repulsion $U$ at integer (fixed) filling. Finally, for larger values of the exchange coupling $J^{\prime}$, the entanglement spectrum is almost flat. We exemplify this feature in Fig. 8, where the entanglement spectra for both dimer models with exchange coupling $J^{\prime}=10.0$ are shown. This behaviour is indeed in agreement with the fact that the dimers are almost isolated for $J^{\prime}\gg 1$, see discussion at the end of Sec. V.1.

After the proposal of Li and Haldane haldane-es that the low-lying entanglement spectrum can be used to identify topological order, a series of papers has been devoted to study the entanglement properties of two-dimensional topological phases, see, e.g., Ref. [2] from Ref. alba13 . On the other hand, the more conventional phases realized in two-dimensional systems have been received less attention review-nicolas . In the latter case, an interesting result is due to Alba et al. alba13 who showed that the entanglement spectrum of the Mott-insulator phase of the square lattice two-dimensional Bose-Hubbard model is dominated by degrees of freedom located at the system-subsystem boundary. In particular, the entanglement spectrum can be interpreted as the spectrum of a (boundary) tight-binding model whose sites are at the system-subsystem boundary. Assuming that this is indeed a quite general feature of a gapped phase, we expect that the entanglement spectra shown in Figs. 8(a) and (b), which are derived for line subsystems $A$, might be characteristic of the columnar-dimer and staggered-dimer models, i.e., such a features might be found in entanglement studies regardless the subsystem shape.

VI Summary

In this paper, we have studied the (quantum paramagnet) VBS phases of the columnar-dimer and staggered-dimer Heisenberg AFMs on a square lattice within the bond-operator formalism at the harmonic approximation. In particular, these results, combined with a procedure employed in Refs. song11 ; luitz15 ; alet15 for magnetic ordered phases, allowed us to calculate the bipartite von Neumann and second Rényi entanglement entropies for the VBS ground states of the two dimer models. Choosing an one-dimensional (line) subsystem $A$, this formalism provides the area law behaviour for the entanglement entropies as expected for gapped phases.

It would be interesting to apply the bond-operator based approach discussed here, e.g., for rectangular strip and square subsystems $A$ [see, e.g., Fig. 1 from Ref. alet15 ]. Such studies would allow us to check whether our results for the entanglement entropies depend on the shape of the subsystem $A$. In this case, however, the eigenvalues $\mu^{2}_{k}$ of the correlation matrix $C$ should be numerically determined. Moreover, it would be important to determine the effects of the cubic $H_{3}$ [Eq. (23)] and quartic $H_{4}$ [Eq. (24)] terms of the effective boson model (15) on the entanglement entropies. These two terms could be perturbatively considered as done, e.g., in Ref. doretto12 . Recall the role of the cubic term in the distinction between the two dimer models as discussed in Sec. II. As expected, the mean-field results obtained here are qualitatively similar for both dimer models. We intend to performed these two studies in a future publication.

Finally, it would also be interesting to compare the entanglement spectra derived here with the ones determined via density matrix renormalization group (DMRG) calculations as done, e.g., for the square lattice Bose–Hubbard model alba13 . However, as far as we know, such DMRG data for the columnar-dimer and staggered-dimer models are not available at the moment.

Acknowledgements.

We thank E. Miranda for helpful discussions. L.S.G.L. kindly acknowledges the financial support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

Appendix A Classical dimerized Heisenberg models

The classical phases of the columnar-dimer [Eq. (9)] and staggered-dimer [Eq. (10)] models can be determined by parametrizing the spins ${\bf S}^{1}_{i}$ and ${\bf S}^{2}_{i}$ [see Figs. 1(a) and (b)] as

[TABLE]

where ${\bf Q}$ is the ordering wave vector, ${\bf R}_{i}$ is a vector of the dimerized lattice, and the set of unit vectors $\hat{e}_{i}$ obeys the following relations

[TABLE]

Substituting Eq. (LABEL:s1s2) into the Hamiltonian (9) of the columnar-dimer model, we obtain the energy $E$ as a function of the components $Q_{x}$ and $Q_{y}$ of the ordering wave vector and the angle $\theta$, namely

[TABLE]

The ground state energy follows from the minimization of Eq. (51) with respect to the parameters $Q_{x}$, $Q_{y}$, and $\theta$. For $J^{\prime}>0$, we find ${\bf Q}=(0,\pi)$ and $\theta=\pi$, i.e.,

[TABLE]

With the aid of Fig. 1(a), one easily sees that the configuration (52) corresponds to a collinear Néel phase.

Similarly, substituting Eq. (LABEL:s1s2) into the Hamiltonian (10) of the staggered-dimer model, we arrive at

[TABLE]

where the dimer nearest-neighbor vectors $\tau$ are given by Eq. (8). In this case, we find that ${\bf Q}=(0,0)$ and $\theta=\pi$, i.e.,

[TABLE]

Appendix B Entanglement spectrum and the correlation matrix

In this section, we derive the relation (43) between the eigenvalues $\mu^{2}_{k}$ of the correlation matrix $C$ and the entanglement spectrum $\epsilon_{k}$.

As mentioned in Sec. V, the generic quadratic Hamiltonian (36) and the corresponding entanglement Hamiltonian (38) can be diagonalized by the same Bogoliubov transformation. Let us consider

[TABLE]

where the coefficients $\phi_{ki}$ and $\psi_{ki}$ are assumed to be real for simplicity. Due to the bosonic algebra of the $b_{k}$ operators, the coefficients $\phi_{ki}$ and $\psi_{ki}$ obey the relations

[TABLE]

which imply that

[TABLE]

For sites $i$ and $j$ associated with the subsystem $A$, we have

[TABLE]

where the states $|\xi_{A}\rangle$ and $|\xi_{\bar{A}}\rangle$ are respectively associated with the subsystem $A$ and its complementary $\bar{A}$ as defined in Sec. I. Moreover, since the Bogoliubov transformation (55) diagonalizes the entanglement Hamiltonian $\mathcal{H}_{E}$ [Eq. (39)] and the reduced density matriz $\rho_{A}$ has the form (37), we have

[TABLE]

Therefore, from Eqs. (57) and (58) and the inverse of the transformation (55), one shows that

[TABLE]

and, similarly,

[TABLE]

The single-particle Green’s functions (42) assume the form

[TABLE]

where $i$ and $j$ $\in\,A$. With the aid of the orthogonality condition (56), one shows that

[TABLE]

The above equation can be written in a matrix form

[TABLE]

where the elements of the $N_{A}\times N_{A}$ matrices $G^{++}$ and $G^{--}$ are given by

[TABLE]

and the vectors $\hat{\psi}_{k}$ and $\hat{\phi}_{k}$ are defined as

[TABLE]

From Eq. (63), we find the eigenvalue equation for the correlation matrix C,

[TABLE]

that provides the relation (43) between the eigenvalues $\mu^{2}_{k}$ of the correlation matrix $C$ and the entanglement spectrum $\epsilon_{k}$.

Finally, one notices that

[TABLE]

which is the bosonic version of Eq. (16) from Ref. peschel03 written in a slightly different notation.

Appendix C Eigenvalues of the correlation matrix for one-dimensional (line) subsystem

For an arbitrary subsystem $A$, the correlation matrix (41) satisfies the property $C_{ij}=C_{|i-j|}$, as one can easily see from Eq. (46), indicating that the correlation matrix $C$ is a Toeplitz matrix gray . In particular, for an one-dimensional subsystem $A$ with periodic boundary conditions, one finds that the correlation matrix $C$ is indeed a circulant matrix: in this case, each row of the matrix can be obtained from the first row by a shift of the matrix elements gray . Due to this translational property, the eigenvalues of a circulant matrix can be obtained by a discrete Fourier transform of its first row elements.

Let us consider the one-dimensional subsystem $A$ shown in Fig. 1(c), i.e., a chain of size $L^{\prime}$ and $N_{A}=(L^{\prime}+2)/2$ sites. Since the vectors of the underline dimerized lattice are ${\bf R}_{i}=2i\hat{x}$, with $i=1,2,\ldots,N_{A}$, the expression (46) for the single-particle Green’s function $f_{ij}$ and $g_{ij}$ can be written as

[TABLE]

where the functions $\alpha_{k_{x}}$ and $\beta_{k_{x}}$ are defined as

[TABLE]

The eigenvalues $\mu^{2}_{k}$ of the correlation matrix C are given by the discrete Fourier transform gray

[TABLE]

where $C_{0j}$ are the elements of the first row of the correlation matrix (41). Using the convolution property $(f*g)(x)=(g*f)(x)$, one shows that the matrix elements (41) can be written as

[TABLE]

where $s-i=x$, $j-i=l$, and $j-s=l-x$. Therefore, Eq. (68) assumes the form

[TABLE]

Finally, from Eqs. (66), (69), and (70) and with the aid of the Fourier transform property $\mathcal{F}\{(f*f)(x)\}_{l}=\mathcal{F}\{(f)(x)\}_{l}\mathcal{F}\{(f)(x)\}_{l}$, one shows that

[TABLE]

with $k_{x}$ given by Eq. (48).