A Replica method.

Entanglement Entropy of Systems with Spontaneously Broken Continuous Symmetry


We study entanglement properties of systems with spontaneously broken continuous symmetry. We find that in addition to the expected area law behavior, the entanglement entropy contains a subleading contribution which diverges logarithmically with the subsystem size in agreement with the Monte Carlo simulations of A. Kallin et. al. (Phys. Rev. B 84, 165134 (2011)). The coefficient of the logarithm is a universal number given simply by where is the number of Goldstone modes and is the spatial dimension. This term is present even when the subsystem boundary is straight and contains no corners, and its origin lies in the interplay of Goldstone modes and restoration of symmetry in a finite volume. We also compute the “low-energy” part of the entanglement spectrum and show that it has the same characteristic “tower of states” form as the physical low-energy spectrum obtained when a system with spontaneously broken continuous symmetry is placed in a finite volume.

I Introduction

In recent years there has been a lot of theoretical interest in entanglement properties of quantum states of matter. Entanglement has proved to be a useful probe of non-local correlations for both gapless and gapped systems. The most commonly used characterization of entanglement is the entanglement entropy , defined as the von Neumann entropy associated with the reduced density matrix of a subsystem. A close relative of the entanglement entropy is the Renyi entanglement entropy .

To date, the most impressive progress has been made in the study of entanglement entropy of one dimensional critical states. Here, it has been shown(1); (2) that for systems described by conformal field theories (CFT’s) the entanglement entropy of a system of length behaves as,


where is the short distance cut-off. The coefficient is a universal number known as the central charge of the CFT. The subleading constant depends on the system geometry (e.g. the ratio of the subsystem size to the full system size). Although is not fully universal as is clear from the cut-off dependence of (1), it is universal up to an additive constant. The universal behavior of the entanglement entropy (1) has proved useful for extraction of the central charge of the governing CFT in numerical density-matrix renormalization group (DMRG) studies of one-dimensional critical systems.

Our present understanding of entanglement in dimension is far less complete. However, it is generally expected that for both critical and non-critical systems the leading contribution to the entanglement entropy scales as the area of the subsystem boundary , .1 This “area law” contribution is related to short-range entanglement in the vicinity of the boundary, and as a result, the proportionality constant is non-universal. However, for critical scale invariant systems one expects a subleading, fully universal, geometric contribution to the entanglement entropy,


The scaling form (2) is believed to hold for scale invariant systems in with arbitrary smooth subsystem boundary and in with flat subsystem boundary.(5); (6); (7); (8) Additional logarithmic contributions are expected in if the boundary has sharp corners and in if the boundary is curved. We note that all the above stated results/hypotheses on the entanglement entropy apply also to the Renyi entanglement entropy , with the constants and acquiring a dependence on .

The scaling hypothesis (2) relies on the following argument. The entanglement entropy, being a dimensionless quantity, can only depend on ratios of length (or energy) scales. However, in a scale invariant theory, the only two length scales are the total system size and the short distance cut-off , with the corresponding energy scales and , where is the dynamical critical exponent. Thus, any dependence of the entanglement entropy on the system size must come together with the dependence on the short distance cut-off. The variation comes from short-range entanglement in the vicinity of the boundary and is expected to take the form of some local geometric quantity integrated over the boundary area. If the boundary is straight and has no corners, there are no non-trivial local geometric quantities, and the only possible cut-off dependent contribution to the entanglement entropy is proportional to the integral of over the boundary, i.e. the boundary area, in accord with Eq. (2). On the other hand, in two dimensions, if the boundary of the subsystem has corners, each corner can contribute a constant to resulting in an entanglement entropy which depends logarithmically on and hence on . For a generalization of these arguments to curved boundaries, see Refs. (9); (10); (11).

In a recent breakthrough it has become possible to extract the Renyi entanglement entropy of quantum systems using Monte-Carlo simulations.(12) One of the first applications of the method of Ref. (12) has been to study the entanglement entropy of a spin- Heisenberg model on a two-dimensional square lattice.(13) The Monte-Carlo simulations were performed using a torus geometry, with the subsystem being either a cylinder or a square. As expected, the leading contribution to the Renyi entanglement entropy was found to scale linearly with the system size. However, surprisingly, a subleading logarithmic correction was observed for both geometries studied.

The ground state of the Heisenberg model on a two-dimensional square lattice is known to spontaneously break the spin-rotation symmetry to a subgroup. Thus, in the infinite volume limit the ground state is infinitely degenerate and the ground state manifold is a two-dimensional sphere, whose points are labeled by the orientation of the Néel order parameter . The low energy excitations above a particular ground state are two linearly dispersing Goldstone bosons, known as spin-waves. The interactions between the spin-waves are irrelevant in the low-energy limit and one, thus, simply obtains a theory of two free scalar bosons. Naively, this theory is scale invariant and so according to Eq. (2), one does not expect any subleading logarithmic contributions to the entanglement entropy as long as the boundary of the subregion is smooth, e.g. for the cylindrical subregion geometry. For the square subregion geometry, one could attribute the logarithmic correction to the corners of the boundary, however, the Monte-Carlo estimate of the coefficient of the logarithm differs both in sign and by one order of magnitude compared to the previous calculation based on a free bosonic field theory.(8)

The caveat to the above discussion is that a system with a spontaneously broken continuous symmetry is, in fact, not scale invariant at low energy in the same way that a CFT is. The reason for this is that symmetry is always restored in a finite volume and instead of the degenerate vacuum manifold one has a unique ground state and a “tower” of excited states carrying a definite charge under the symmetry group. For instance, in the case of the Heisenberg model with its symmetry breaking, the tower of states can be described by an effective Hamiltonian,


where is the total spin of the system, is the spin-wave velocity, is the spin-stiffness and is the volume of the system. Note that the spacing between the energy levels of the tower scales as . For system dimension , where spontaneous breaking of continuous symmetry is possible, this spacing is much smaller than the spin-wave gap .

Thus, in a system with spontaneous breaking of continuous symmetry, to a given length scale there correspond two infra-red energy scales: the tower of states spacing and the spin-wave gap . The ratio of these energy scales corresponds to the square of the magnitude of order parameter fluctuations associated with a spin-wave mode (normalized by the total size of the order parameter manifold). In general, we expect the entanglement entropy to depend on this ratio, so by dimensional analysis


Note that in the large limit, is much greater than one.2 In this paper, we demonstrate that for a system with symmetry breaking, for a smooth subsystem boundary in (straight boundary in ),


Thus, the entanglement entropy contains an extra subleading term which diverges logarithmically with the system size. The coefficient of this term is directly expressed in terms of the number of Goldstone modes . In particular, for the case of the Heisenberg antiferromagnet , . We note that in addition to the logarithmic correction in Eq. (5) there also appears a finite constant that depends on the system geometry. Unlike in the case of entanglement entropy in one dimension, Eq. (1), where the presence of a logarithm rendered universal only up to an additive contribution, here is fully universal, as all the short-distance physics is absorbed into the order-parameter stiffness and the Goldstone velocity . Given a subsystem geometry, our calculation method allows us to determine numerically; in Fig. 1 we present the result for the geometry studied in the Monte-Carlo simulations of Ref. (13); (19).

Figure 1: The universal constant contribution , Eq. (5), to the Renyi entropy of the non-linear -model in dimension (describing e.g. the square lattice Heisenberg model). In the geometry studied here, the total system is a torus of size , while the subsystem is a cylinder of size . The details of the calculation are described in section III.3. The solid line is a guide to eye.

We also extend the result (5) to the case when and the boundary has corners, or and the boundary is curved. Here we obtain,


The coefficient can be computed in a free theory of Goldstone fields. For instance, for the case with corners,


The sum in Eq. (7) is over the corners of the boundary; denotes the corner log coefficient in a free scalar bosonic theory. This coefficient depends on the corner angle and has been computed in Ref. (8).

In addition to the entanglement entropy, we compute the Renyi entropies , which are also found to satisfy the scaling forms in Eqs. (5), (6). The coefficient of the area law now depends on the replica index , so does the finite size constant and the corner/curvature coefficient . However, the coefficient of the logarithmic correction is found to be independent of the replica index.

We also study the full low energy spectrum of the entanglement Hamiltonian defined in terms of the reduced density matrix of the subsystem as . We find,


Here denotes the spin of the subsystem and , are bosonic annihilation/creation operators. The first term in Eq. (8) has the same form as the tower of states sector (3) of the physical Hamiltonian, while the second term is harmonic. The coefficient to logarithmic accuracy is given by


The gaps in the harmonic sector are found to scale as . Thus, in where spontaneous breaking of continuous symmetry is possible, the entanglement gap in the tower of states sector is parametrically smaller than the gap in the harmonic sector.

We note that recent DMRG studies of the entanglement spectrum in the superfluid phase of the 2d bosonic Hubbard model are in reasonable agreement with our analytical result (8).(15) In particular, Ref. (15) has observed a tower of states structure of the low-lying entanglement spectrum with the entanglement gap scaling inversely with the length of the subsystem boundary, as predicted in Eq. (9). Note that in a superfluid, the system displays spontaneous breaking of symmetry, so the role of is played by the subsystem particle number (more precisely, its deviation from the ground state expectation value ). The higher lying states in the entanglement spectrum were found in Ref. (15) to exhibit a much weaker system size dependence, consistent with our result .

For symmetry, the result (8) also applies to 1d superfluids/XY magnets, which are described by Luttinger liquid theory. For this 1d case, the tower of states gap and the harmonic gap in the entanglement spectrum scale in the same way as . This is not surprising: Luttinger liquid theory is a CFT, and a powerful result on 1d CFTs states that the entanglement spectrum of a segment of length embedded in an infinite line is identical to the physical spectrum of the system on an open segment of length .(1) Thus, for , the first term in Eq. (8) is the tower of states part of the physical spectrum on an open segment and the second term is the spin-wave part of the physical spectrum on an open segment. We stress that such direct full correspondence of physical spectrum and entanglement spectrum is exact only in . We also note that recent DMRG simulations of Ref. (16) have provided a numerical confirmation of this correspondence for a number of systems described by 1d CFTs, including the superfluid phase of the 1d Bose-Hubbard model and the Luttinger-liquid phase of the XXZ chain.

The results of the present paper are based on the analysis of the non-linear -model, which provides the full low-energy description of the system. We find the ground state wave-function for the non-linear -model in a finite volume, compute the reduced density matrix and obtain its spectrum, Eq. (8). From this spectrum, we then compute the entanglement entropy. This procedure is very similar to the correlation matrix technique used to find the entanglement spectrum and entanglement entropy in a free bosonic field theory,(14) except we pay special attention to the compact nature of the order parameter manifold. We also separately confirm the result on entanglement entropy (5) by performing replica method calculations on the non-linear -model in the path-integral formulation (see Appendix). We stress that all our results are exact in the large system limit.

Curiously, to get an intuition about our exact results it is sufficient to consider a simple quantum mechanical model of two coupled quantum rotors and , representing the average order parameter in subsystem and its complement ,


Here is the angular momentum of each rotor. , and denote the volumes of each subsystem and of the total system respectively. We choose the coupling appropriately to reflect the finite order-parameter stiffness of the system. One finds that the entanglement Hamiltonian corresponding to the ground state wave-function of this model, indeed, has a tower of states form and reproduces the logarithmic correction to the entanglement entropy in Eq. (5).

We would like to note that the presence of logarithmic corrections to the entanglement entropy of a Heisenberg antiferromagnet has been theoretically pointed out previously in Ref. (17). The authors of Ref. (17) have used the spin-wave (large ) expansion of the Heisenberg model. The effect of symmetry restoration in a finite volume has been mimicked by an application of a staggered magnetic field , suitably chosen to give a zero net staggered moment. The final step of the calculation of the entanglement entropy has been performed numerically and finite size scaling was used to extract the coefficient of the logarithmic divergence . This value is quite close to our exact result . Furthermore, we have repeated the calculation in Ref. (17) for system sizes up to sites and found that approaches unity (within a percent) as the system size is increased.3 While our calculation of the coefficient agrees with the semi-numerical method of Ref. (17), our approach has the advantage of being exact: in particular, it correctly treats the compactness of the order parameter and does not rely on the expansion.

An attempt to make a connection between the restoration of symmetry in a finite volume and the appearance of subleading logarithmic terms in the entanglement entropy has also been made in Ref. (13). Here the authors started with a mean-field Néel state of a Heisenberg antiferromagnet and averaged it over all orientations of the order parameter. The entanglement entropy of the resulting spin-singlet state was found to be , where is the total number of lattice sites (here and below we drop the constant piece in ). However, this approximation is too crude: it gives the coefficient of the logarithmic divergence to be , instead of our exact result . The physical reason why the approach of Ref. (13) fails is that it does not take into account the existence of spin-waves, while as we have argued above, the presence of both the tower of states and spin-waves is crucial for the logarithmic correction to entanglement entropy in Eq. (5). The estimate of Ref. (13) for an antiferromagnet is physically similar to the exact result one obtains in the case of a free Bose gas where , with - the total number of particles and - the average interparticle distance. On the other hand, for an interacting Bose gas, i.e. a superfluid, our result, Eq. (5), with gives a logarithmic correction . Thus, based on entanglement entropy we can distinguish a free Bose gas, where no Goldstone boson is present, from an interacting superfluid, which has a Goldstone mode.

This paper is organized as follows. As a warm-up, we begin by calculating the entanglement properties of the toy rotor model (10) in section II. In section III, we determine the ground state wave-function of the non-linear -model and use it to compute the entaglement spectrum and the entanglement entropy. Some concluding remarks are presented in section IV. In the appendix, we give an alternative calculation of the entanglement entropy in the non-linear -model using the replica method and show that the result is in complete agreement with of our wave-function calculation in section III.

Ii Warm up: rotor model.

In this section, we calculate the Renyi entropy in the model (10), describing two quantum rotors and of unit length. These rotors are taken to represent the average orientation of the order parameter in subsystems and . We choose the coupling between the rotors in the following way. Recall that when the system is placed in a box of size , the energy cost to twist the order parameter by an angle between the two sides of the box is related to the order-parameter stiffness via,


Thus, we take . Note that this value of is approximate and should be understood in a scaling sense only.

To find the spectrum of the model (10) it is convenient to first work with the Lagrangian formulation,


Let us introduce the average and relative coordinates and via,


Here, , and are unit vectors forming an orthonormal set together with . We expect the fluctuations of the average coordinate to be parametrically slower than those of the relative coordinate . Moreover, we expect the fluctuations of to be small. Therefore, expanding the Lagrangian (12) to leading order in and in derivatives of we obtain,


Here is the total volume and is the reduced volume. We see that the Lagrangians of the average and relative coordinates decouple. The dynamics of the average coordinate are governed by the quantum rotor Lagrangian (15); the associated Hamiltonian


is, indeed, the appropriate “tower of states” Hamiltonian describing the lowest energy excitations of a system with spontaneous symmetry breaking. Here , , are angular momentum operators, i.e. generators of rotations of , and hence of . They can, thus, be identified with the total angular momentum (spin) of the system. We use the short-hand , which is also equal to the Laplacian, , on the sphere.

The dynamics of the relative coordinate are governed by the Lagrangian (16) describing an dimensional harmonic oscillator. The frequency of this harmonic oscillator is given by


In the true physical system, the spectrum of the “relative motion” involves Goldstone modes with a dispersion , where the momentum is quantized in a finite geometry. Our rotor model (10) replaces this multi-mode spectrum with a single dimensional oscillator whose energy (18) is of order of the finite-size gap of the Goldstone modes. It turns out that such a replacement is sufficient for capturing the logarithmic contribution to the entanglement entropy in Eq. (5).

The ground-state wave-function corresponding to the Lagrangian (14) is a product of a wave-function of the rotor Hamiltonian (17) and the ground state wave-function of the harmonic oscillator (16),


Here, is the volume of a unit sphere and


Note that the condition guarantees that the amplitude of the relative fluctuations is, indeed, small.

We proceed to compute the reduced density matrix from the wave-function (19),


The contributions to the integral over above come from . Therefore, is non-negligible only for . We may change variables in Eq. (21) to and, to leading order, take both and to lie in the tangent plane of . The integral over then becomes Gaussian and gives,


Since is just the Laplacian on the sphere, we may use the heat kernel expansion


Thus, as , we may write,


Defining the entanglement Hamiltonian as , we obtain,


Hence, the entanglement Hamiltonian has the same form as the physical tower of states Hamiltonian (17). The entanglement gap in Eq. (25) is found to scale as . In the next section we will perform an analysis of the full non-linear -model, instead of the toy rotor model studied here, revealing that the lowest branch of the exact entanglement Hamiltonian, indeed, has a tower of states form , with given by Eq. (9). Thus, up to a logarithmic factor, the scaling of the entanglement gap in Eq. (25) agrees with the exact result (9). Note that in Eq. (25) is determined via Eq. (20) by the coupling of the rotor model that we simply postulated to scale as . By choosing appropriately, we can force the tower of states spectrum (25) to match the exact result (9).

We can now compute the Renyi entropy of the rotor model. From Eq. (24)


where we’ve used Eq. (23) in the last step. Therefore,


Note that if we wish to interpret the result (28) in terms of the actual physical system rather than the toy rotor model, we must remember that Eq. (28) only captures the contribution of the lowest branch of the entanglement spectrum given by Eq. (25). There will also be a contribution from the higher lying part of the entanglement spectrum. In the next section, we will determine this higher lying part, and show that it contributes a standard area law term to the Renyi entropy, so we identify the Renyi entropy of the rotor model, Eq. (28), with the logarithmic correction in Eq. (5). In addition, since the result (27) depends on (), and hence on the details of the spectrum of “relative fluctuations” (Goldstone modes), it is only logarithmically accurate. In particular, it does not capture the universal geometric constant of Eq. (5).4 However, our full calculation in the next section will allow us to determine , as well.

Iii Wave-function method.

In this section, we calculate the entanglement entropy in the non-linear -model using the wave-function method. The action of the theory is given by,


with , , an -dimensional unit vector. We will set the spin-wave velocity below and restore it in the final results. For simplicity, we consider the theory on a spatial torus, although a system with open boundaries can be treated with minimal modifications.

Following the standard treatment, we write


Here, . and are unit vectors that together form an orthonormal basis: , . The fields are constrained to satisfy,


The vector describes the (slow) fluctuations of the overall direction of the order parameter in the system, while describe the spin-wave fluctuations. Expanding the action (29) to leading order in , we obtain


where is the volume of the system. We remind the reader that higher order terms in the expansion in are irrelevant in the RG sense for . Thus, to leading order the motion of and decouples. The Hamiltonian corresponding to the action (32) is given by,


where desribes the sector and - the sector. The action for is that of a particle on a unit sphere, so


Here are angular momentum operators, i.e. generators of rotations of , and hence of . They can, thus, be identified with the total angular momentum (spin) of the system. We use the short-hand , which is also equal to the Laplacian, , on the sphere.

The action for is quadratic, so the spin-wave Hamiltonian takes the form




and are the momenta allowed by the periodic boundary conditions on the torus. Note that the mode is missing from the sums (35), (36) due to the constraint (31). For future reference, we introduce the static propagator,




Let us now write the ground-state wave-function(al) of the system, . It is given by a product of ground state wave-functions in the and sectors. The ground state in the sector carries zero angular momentum and the corresponding wave-function is a constant. The ground state wave-function in the sector is a Gaussian; its form can be easily deduced by requiring that the wave-function reproduce the correlator (37). We, thus, obtain






In Eq. (39) and below, we ignore the normalization of the wave-function as it will not be essential for our purposes. Thus, the overall ground state wave-function of the system is,


where the last equality holds to leading order in .

Next, we divide our system into region and its complement and compute the reduced density matrix associated with region (at this point we keep the shape of regions and arbitrary).


Here, and denote the values of restricted to regions and , respectively. Recall that we are considering configurations of with small, smooth flucutations about a global direction , i.e. each maps the space into some small patch of the sphere. Since the values of in the two wave-functions in the integrand of Eq. (43) are identified, we can also take , , as well as to lie in the same small patch of the sphere. We then compute patch by patch. By rotational symmetry, it is sufficient to compute in the patch centered at the North pole. Then, write


Here, the index on runs over , and unlike , is unconstrained. Expanding to leading order in ,


Here and below we use the short-hand notation . and denote restricted to regions and , respectively. Similarly, denotes with both arguments restricted to region , denotes with both arguments restricted to region , and denotes with restricted to region and restricted to region . Performing the integral over in Eq. (45),


where is defined with its arguments over region and satisfies,


Note that although defined over the entire system does not have an inverse due to the presence of the zero mode, - the restriction of to region generally does possess an inverse. Now from Eq. (40) we deduce,


where is defined with its arguments over region and satisfies,


Again, even though defined over the entire system does not possess an inverse, generally does. is also defined with its arguments over region and is given by, , with - the volume of region . In other words, is the projector onto the constant mode over region . Now using Eq. (48), takes the form,


Using Eq. (44), we can now rewrite in terms of , ,


Note that invariance is now restored in Eq. (51) so it can be used on the entire sphere (and not just in the vicinity of the North pole).

Next, we would like to deduce the spectrum of the reduced density matrix (entanglement Hamiltonian). To do this, it is convenient to parametrize as,


Where and , , are unit vectors forming an orthonormal basis, , . We take the fields to satisfy the constraint,