Spin wave contribution to entanglement in Heisenberg models

Employing the Holstein–Primakoff transformation, we obtain an approximate description of both the ground state entanglement and the entanglement dynamics in a spin system in terms of spin waves. As an example, we apply the method to a two-dimensional spin system.


Introduction
Extensive attention has been devoted in the last few years to understanding the overlap between condensed matter and quantum information theory [1]. Most of the models studied so far have dealt with one-dimensional spin systems as they are amenable to an exact solution. O'Connors and Wootters showed that in the Heisenberg chain the maximization of the entanglement at zero temperature is related to the energy minimization [2]. At finite temperatures, violation of Bell's inequality can be directly related to the properties of the internal energy [3]. Finite temperatures are not necessarily detrimental to the amount of entanglement present in the system; it has been shown indeed that both temperature and magnetic field can increase it [4,5]. Near a quantum phase transition, entanglement can be classified in the framework of scaling theory [6]- [9] although there are notable differences between non-local quantum and classical correlations. Interest was confined not only to the equilibrium properties of entanglement but also to its dynamical properties (see [10] and references therein).
In two dimensions, an analysis of entanglement is far more difficult as one cannot resort to exact solutions. Roscilde et al [11] studied the entanglement of a class of anisotropic magnetic models by means of Monte Carlo simulation. In the present paper, we would like to address the question of to what extent it is possible to use approximation methods, such as the spinwave analysis, developed in the theory of magnetism [12] to study the entanglement. We discuss separately the pairwise entanglement of the ground state and the connections between entanglement dynamics and spin wave propagation.

The model
The system we consider is a spin-1/2 ferromagnetic lattice with an exchange coupling J in a transverse magnetic field of strength h. The Hamiltonian is given by where S a l are the spin-1/2 matrices (a = x, y, z) relative to site l of a D-dimensional spin lattice, and N is the total number of sites (we assume periodic boundary conditions). This Hamiltonian allows us to study a range of systems from the isotropic XY to the quantum Ising model by varying the anisotropy parameter 0 < γ 1.
We concentrate on the two-site entanglement, which can be evaluated by analysing the twoqubit reduced density matrix ρ (2) ı,  , obtained after tracing out all the spins except those at positions ı and . In this case, the amount of entanglement for an arbitrary mixed state of two qubits can be quantified through the concurrence, C[ρ (2) ], defined as [13] C[ρ (2) R = ρ (2) σ y ⊗ σ y (ρ (2) ) * σ y ⊗ σ y , where λ max is the largest eigenvalue of the matrix √ R. The great advantage of the concurrence over other entanglement measures [14], is that it is directly defined in terms of the density matrix ρ (2) , without any minimization procedure. It follows that C[ρ (2) ] can also be obtained from the two-point correlation functions of the corresponding spins.
In order to evaluate the concurrence it is essential to identify the structure of the two-qubit reduced density matrix ρ (2) . In the basis {|0 = |↑↑ , |1 = |↑↓ , |2 = |↓↑ , |3 = |↓↓ }, we will consider density matrices of the following form (see however [15]) Due to the symmetries of the Hamiltonian, this structure is preserved by the unitary time evolution. For this class of mixed states, the concurrence is and can be expressed in terms of a spin-spin correlation function. Using the notations g αβ =  /2 for correlation functions relative to spins at sites ı and , we obtain for the concurrence For states which are translational invariant δ z = 0, i.e. ρ 11 = ρ 22 in equation (4). Moreover, in the presence of reflection symmetry, g xy = g yx , and therefore ρ 12 is real. These two conditions are satisfied for the ground state (in the absence of symmetry breaking, see below) and for the time evolution starting from a singlet state. We can, then, use the simplified expression The class of models in equation (1) can be diagonalized exactly by means of the Jordan-Wigner transformation [16,17], which maps spins into spinless fermions, in one dimension. Also the correlation functions can be evaluated exactly for both translational invariant [18] and nontranslational invariant states [19]. This is the approach pursued in [10] to analyse the entanglement dynamics.
In the following, we will employ the spin-wave approach to obtain an expression for these correlation functions, and then for the concurrence. This method is approximate, but it can be employed also for higher-dimensional lattices.

Spin-wave approximation
The spin-wave approximation is a standard tool in the theory of magnetism. A detailed presentation can be found, e.g. in [12]. Here, we review its essentials to keep the presentation selfcontained. The procedure consists of three steps. First, we need to identify a suitable reference state (typically, the 'classical' ground state). Then, we obtain the spin waves as single-particle boson excitations defined over this reference state. Finally, we set the occupation of the boson modes as a function of the temperature of the system. In order to simplify the notation, for the ground state entanglement, we consider the case of γ = 0 since the general case with 0 γ 1 leads only to quantitative corrections. In contrast, we will see in the following sections that this is not the case for the dynamics; therefore, we will come back to the full Hamiltonian (1) in section 2.
The classical (i.e. factorized) state, which minimizes the energy, is the zeroth order term in an expansion of the Hamiltonian that will be performed below. For γ = 0, the model shows a rotational invariance in the equatorial plane of the spin space. Due to this symmetry, we can choose the common spin direction in the x-z plane, and set S x = S cos ϕ, S y = 0, S z = S sin ϕ. Here, S is the value of the spin, taken to be 1/2 at the end of the calculations. The value of ϕ depends on the ratio h/J (and on the coordination number z), reflecting the quantum phase transition. To write the results in the simplest way, we define the parameter χ = h/(zJS). Minimization of the average energy with respect to ϕ gives The case ϕ = 0 corresponds to the ordered phase, with spins aligned along the field axis; the case ϕ = 0, called canted phase, is mainly due to spin exchange. We can express the spin operators in terms of Bose operators by performing the Holstein-Primakoff transformation; a subsequent 1/S expansion allows us to obtain a (approximate) description of the quantum fluctuations around the reference state defined above.
To begin with, we need to rotate the spin axes so that the new quantization (i.e. z) axis points in the direction of the reference spin S . Then, we will relate the (rotated) spin operators to boson ones. It is clear that we need to rotate the spin around the y-axis by an angle −ϕ, whose value depends on the parameter χ. Letting the rotated operators to be labelled by a tilde, we have  Using the rotated operators, we perform the transformatioñ whereS ± =S x ± iS y are the spin raising and lowering operators, while the vector index l labels the lattice positions. Equations (10)- (12) can be employed for any value of S to transform the spin Hamiltonian into one for an interacting boson system. We retain only the first orders of a 1/S expansion, corresponding to a quadratic Hamiltonian in the Bose operators: Equation (13) can be easily diagonalized by going to the momentum representation and then performing (only if χ < 1) a Bogoliubov transformation. In this last case we have with where the primed sum runs over oriented nearest-neighbour bonds. By performing the linear transformation with tanh 2θ µ = µ J µ , one gets where the dispersion relation for the spin waves is given by The correlation functions will depend on the mean values of combination of two of the original (direct lattice) boson operators. They are given by where α † µ α µ are given by the usual Bose factors.
The concurrence is evaluated through equation (6). We note that, although this procedure is quite commonly employed to obtain equilibrium properties of a spin system, we will use it here to describe not only the structure of the ground state, but also the dynamical evolution, starting from an out-of-equilibrium initial state. We remark that we cannot expand the concurrence in powers of 1/S and then take the first orders for S = 1/2, since it is well defined only for S = 1/2. Rather, we express the various correlation functions in terms of the bose operators, expand them in powers of 1/S, retaining only the first-order contributions and employ these approximate functions to evaluate the concurrence.

Ground state concurrence
In this section, we relate the entanglement content of the ground state to the properties of the spin waves.

Concurrence for χ > 1
In this case ϕ = 0 ⇒ a † ı a †  = a ı a  = 0, giving a very simple form for the two-site reduced density matrix (we put S = 1/2 in the following) (the mean occupation number n ı ≡ a † ı a ı is independent of the site label ı due to the translation invariance). By means of (7) and (23) the concurrence simplifies considerably and can be cast in the form It is evident how spin-waves are responsible for the entanglement in the system.

Concurrence for χ 1
In this case, the Holstein-Primakoff transformation describes the excitations of the system with respect to a reference state which breaks the rotational symmetry of the Hamiltonian (this is due to the fact that the reference spin direction has a component in the equatorial plane). As a result, the density matrix does not display any of the symmetry properties that led to equation (4) (apart from translation invariance). Nevertheless, the procedure to obtain the concurrence can still be carried out analytically, and the result turns out to be expressible in very simple terms where the average values can be taken from equations (21) and (22). Note that this result reduces to equation (24) for χ > 1.

Entanglement propagation
In this section, we will employ a procedure similar to the one discussed above, to obtain an approximate expression of both the correlation functions and the reduced density matrix in terms of time-dependent Bose operators. To this end, we will first seek the time dependence of the creation and annihilation operators. We will obtain this time evolution in general and then illustrate the results for two analytically solvable cases: (i) the isotropic model (γ = 0) and (ii) the case in which there is no external field, h = 0 in equation (1), with the dynamics of the spin lattice generated only by the neighbour site-spin interaction. This last case will be analysed in the presence of the anisotropy parameter γ = 0. This analysis will give us the opportunity to illustrate the limitations of the overall procedure described so far.
Our aim is to describe the distribution of pairwise entanglement in the lattice as a function of time, starting either from a completely separable state, such as the vacuum state defined below, or from a well localized entangled pair prepared in a singlet state. In the first case, we will see that only the anisotropic XY model (i.e. the one with γ = 0) can give rise to creation of entanglement from the vacuum. In the second case, we will suppose that a singlet state has been prepared (i.e. superimposed over the vacuum) at sites a and b, and look at the entanglement propagation all over the lattice. In both cases, we will take the vacuum state |vac = |↑, ↑, ↑, . . . , ↑ as our reference state. This is the right choice (independently of χ) since we are limiting ourself to the investigation of the time evolution starting from a configuration in which at most one spin is reversed (i.e. in spin wave language, at most one boson is present). Thus, we can suppose that |vac is the best reference state since the state of the system will contain only small deviations (i.e. very few excitations) with respect to it. It is clear that if the number of excitations becomes too high, then a quadratic Hamiltonian in Bose operators is no longer sufficient.
As the reference state is chosen to have all the spins aligned along the z-direction, we can perform the Holstein-Primakoff transformation directly on the original spin operators S a ı . By using equations (10)- (12) in the Hamiltonian, expanding in powers of 1/S and retaining terms up to first order, we get (for S = 1 2 ) Note that the contributions proportional to γ do not preserve the total number of bosons; in spin language, they connect different spin-sectors of the Hilbert space of the lattice. This Hamiltonian can be diagonalized by a Bogoliubov transformation exactly along the line discussed previously (one can notice that the terms proportional to γ here give rise to the same kind of effect as those proportional to (cos 2 ϕ − 1) there). Namely, after switching to the reciprocal lattice, we can introduce the Bogoliubov operator α µ as in equations (17) and (18), but this time with coefficients determined by where, again, the primed sum runs over oriented nearest-neighbour bonds.
In terms of these rotated operators, the Hamiltonian is diagonal (see equation (19)) with spin-wave frequencies given by By using the simple time evolution of the α µ 's and by inverting the Bogoliubov transformations, one can obtain the time evolution for the boson operators relative to the direct lattice sites: with Equation (28) allows us to get the time dependence of the various correlation functions and, therefore, of the reduced density matrix. The two limiting cases below are amenable of a fairly simple solution and will be discussed in some detail.

The isotropic case: γ = 0
For γ = 0, i S z i is a constant of the motion. In terms of spin-wave analysis, if one starts from a one-particle state (like the singlet (|↑ a ↓ b − |↓ a ↑ b )/ √ 2), then the first order in the 1/S expansion gives the exact result, since all other (many-particle) corrections vanish. Since sinh 2θ µ = 0, the time dependent G functions introduced above simplify to The vanishing of G (2) x (t) above implies that the reduced density matrix is given by equation (4), with ρ 03 = ρ 33 = 0 for states belonging to the sector of the total Hilbert space in which just one spin is reversed with respect to the vacuum state. This implies that the concurrence becomes C = 2|ρ 12 |. By using the correlation functions evaluated to order O(1/S), we can express the concurrence between sites ı and  in the very simple form For N → ∞ (and for a D-dimensional hyper-cubic lattice) G (1) (t) is (up to an irrelevant phase factor) where J n are nth-order Bessel functions. The evolution of the concurrence is illustrated in figure 1. The various plots show the concurrence between the spin at site a (initially entangled) and any other spin in the lattice at different times. At the beginning, the concurrence is different from zero only between spins at sites a and b, but, as time goes on, we find that the spin at a becomes entangled with both its neighbourhood and the neighbourhood of site b. We also observe that the entanglement displays many collapses and revivals while diffusing all around the lattice with unit speed (in units of the coupling constant J and the lattice spacing) in both directions.

Zero field case: h = 0
Another case in which it is possible to analytically obtain a closed form for the time evolution of the boson operators is the case where no field is applied. In this case, as will become clear in the following, our approximated spin-wave procedure breaks down on increasing γ.
To first order in 1/S, the Hamiltonian for this case reduces to Equations (28)-(30) are still valid and, in the continuum limit (N → ∞), we have where we defined a scaled time τ = 2Jt 1 − γ 2 . The general form of the reduced density matrix, expressed in terms of boson operators to order O(1/S), is given by where the time dependence of the various averages can be obtained by using equations (28), (36) and (37). We will calculate the concurrence for two different initial states. First, we will describe the homogeneous generation of bi-partite entanglement from the vacuum state and its subsequent diffusion all over the lattice; as a second example, we will examine once again the case of a singlet state imposed on top of the vacuum. In both cases, one has a ı = 0 ∀ ı; therefore, ρ (2) takes again the form of equation (4).

Vacuum as initial state.
Starting from the vacuum, all the results become translation invariant. The coefficient of the reduced density matrix relative to sites ı and , for example, are found to depend only on the distance r =  − ı. For a bi-dimensional lattice and writing r = x, y, we find The contributions coming from higher orders in the 1/S expansion are non-negligible for γ > 1/2. This roughly indicates a validity limit of the 1/S expansion for large γ (further discussions are provided below). The concurrence is shown in figure 2. We observe that, for short times, every site becomes entangled with its neighbourhood; furthermore, the range of entangled partners increases with time (nearly linearly). However, the value of the concurrence with another spin at a given distance eventually decreases for long times as ∼1/τ. This long time behaviour, which is obtained also for a one-dimensional case in the spin-wave treatment, is in contrast with what is found from the exact one-dimensional solution [10], where a small residual bi-partite entanglement is found to persist at long times. This indicates another limit of the approximated picture: the spinwave analysis is not valid for arbitrarily long times. Intuitively, this long-time failure is due to the fact that (since γ = 0) more and more excitations are generated as time goes on; then we cannot expect the long-time behaviour to be accurately depicted by retaining only few spin-waves. Eventually, for large γ, the spin wave approach fails in the description of the dynamics also for short times. This can be explained by noticing that the vacuum average number of spin excitations of a given wave vector µ increases with γ as vac| a † µ a µ |vac = For γ = 1/2, this number is of the order of 0.1 and it seems plausible that for larger values of γ, the number of excitations becomes too large to neglect interaction among them.
For small values of the anisotropy parameter, however, our approach remains valid and, due to the fact that the correlation functions depend only on the scaled time τ = 2Jt 1 − γ 2 , we can observe in-phase oscillations of the concurrence as a function of τ for various values of γ, see figure 3. Again, this feature is lost when higher-order corrections become important.

Singlet initial state.
Another interesting case to analyse is the time evolution of the concurrence for an initial singlet state, superimposed on the vacuum and located at sites a and b, The reduced density matrix is given by equation (38), where the various averages can be evaluated as a function of time using again equations (36) and (37). The expressions for the various matrix elements, however, are quite lengthy, with the only remarkable feature that each of them is given by the sum of two contributions, the first of them being just the corresponding vacuum value. As a result, the concurrence contains three kinds of contributions: first, related to the vacuum concurrence; second, strictly due to the presence of the singlet; and third, due to the interference between the previous two. It turns out that, as a function of time, the concurrence is given in general by a complicated interference pattern between its various components; however, due to the small value of γ, the contribution due to the 'propagation' of the singlet dominates. For example, in figure 4, we show the behaviour of the amount of nearest-neighbour entanglement as a function of time for a square spin lattice.

Conclusions
To summarize, we described the method of spin waves as an approximate tool to evaluate correlation functions, which can be used (through the concurrence) to quantify bipartite entanglement in spin lattices in more than one dimension. We showed that, in many cases, the amount of entanglement can be expressed through few physical properties related to spin waves (such as their average number). Furthermore, we explored the limitations of the method by analysing the anisotropic ferromagnetic XY model in various regimes. We obtained an exact result for one-particle states (such as the singlet) in the absence of anisotropy in the spin-spin interaction (the γ = 0 case), but also showed that the approximation scheme breaks down at long times for other cases and, in particular, for values of the anisotropy parameter γ greater than 1 2 . We argued that this is due to the fact that interaction among spin waves cannot be neglected and, therefore, the 1/S expansion cannot be truncated at a low order in this case.