Entanglement and fluctuations in the XXZ model with power-law interactions

TL;DR

This paper explores how long-range power-law interactions in the XXZ model affect its ground-state properties, revealing three regimes with distinct behaviors and a continuously varying dynamical exponent that influences correlations and entanglement.

Contribution

It identifies and characterizes three distinct interaction regimes in the long-range XXZ model, including a novel medium-range regime with a tunable dynamical exponent.

Findings

Existence of three regimes depending on decay exponent and dimension

In the medium-range regime, the dynamical exponent varies continuously

Long-range interactions significantly alter correlation and entanglement scaling

Abstract

We investigate the ground-state properties of the XXZ model with $1/r^{\alpha}$ interactions, describing spins interacting with long-range (LR) transverse (XX) ferromagnetic interactions and longitudinal (Z) antiferromagnetic interactions, or hardcore bosons with LR repulsion and hopping. The LR nature of the couplings allows us to quantitatively study the spectral, correlation and entanglement properties of the system by making use of linear spin-wave theory, supplemented with density-matrix renormalization group in one-dimensional systems. Our most important prediction is the existence of three distinct coupling regimes, depending on the decay exponent $\alpha$ and number of dimensions $d$: 1) a short-range regime for $\alpha > d + \sigma_c$ (where $\sigma_c = 1$ in the gapped N\'eel antiferromagnetic phase exhibited by the XXZ model, and $\sigma_c = 2$ in the gapless XY ferromagnetic phase), sharing the same properties as those of finite-range interactions ($\alpha=\infty$); 2) a long-range regime $\alpha < d$, sharing the same properties as those of the infinite-range interactions ($\alpha=0$) in the thermodynamic limit; and 3) a most intriguing medium-range regime for $d < \alpha < d+\sigma_c$, continuously interpolating between the finite-range and the infinite-range behavior. The latter regime is characterized by elementary excitations with a long-wavelength dispersion relation $\omega \approx \Delta_g + ck^z$ in the gapped phase, and $\omega \sim k^z$ in the gapless phase, exhibiting a continuously varying dynamical exponent $z = (\alpha - d) / \sigma_c$. In the gapless phase of the model the $z$ exponent is found to control the scaling of fluctuations, the decay of correlations, and a universal sub-dominant term in the entanglement entropy, leading to a very rich palette of behaviors for ground-state quantum correlations beyond what is known for finite-range interactions.

Tables

Table 1. Table 1: Summary of main results for the long-range ordered XY phase (in d = 2 𝑑 2 d=2 and 3 for all α 𝛼 \alpha , and in d = 1 𝑑 1 d=1 for α < 3 𝛼 3 \alpha<3 ) and the Néel phase, concerning the dispersion relation, the scaling of fluctuations and of entanglement. The scaling laws refer to a hypercubic-lattice geometry with linear size L 𝐿 L , or to a A 𝐴 A subsystem with linear extent L A subscript 𝐿 𝐴 L_{A} . The ‘—’ symbols indicate the absence of predictions in the corresponding regimes: in the XY phase for d = 1 𝑑 1 d=1 the area-law term in the scaling of the entanglement entropy becomes a subdominant constant; and the Néel phase is simply absent for α < d 𝛼 𝑑 \alpha<d (see the phase diagram on Fig. 1 ).

XY phase	SHORT RANGE  ($\alpha >d+2$)	MEDIUM RANGE  ()	LONG RANGE  ()
dynamical exponent $z$     ($\omega \sim k^z$)	$z=1$	$z=\left(\alpha -d\right)/2$	$z=0$
$⟨{\left(\delta S^x\right)}^2⟩$ scaling (total     system)	$L^{\max (d,2)}$	$L^{\max (d,2z)}$	$L^d$
$⟨{\left(\delta S^y\right)}^2⟩$ scaling (total     system)	$L^{d+1}$	$L^{d+z}$	$L^d$
$⟨{\left(\delta S^z\right)}^2⟩$ scaling (on the $A$ subsystem)	$L_A^{d-1}\log L_A$	$L_A^{d-z}$	$L_A^d$
scaling of subsystem entanglement entropy:     dominant term	$L_A^{d-1}$ ($\log L_A$ in $1d$)	$L_A^{d-1}$ ($\log L_A$ in $1d$)	—
scaling of subsystem entanglement entropy:      logarithmic term	$\frac{d-1}{2}\log L_A$	$\frac{d-z}{2}\log L_A$	$\frac{d}{2}\log L_A$

Table 2. Table 2: Decay exponent of the correlation functions for α > d 𝛼 𝑑 \alpha>d . “Long.” and “trans.” stand for longitudinal and transerve (with respect to the direction of the order parameter). In the XY phase, z = min ⁡ [ 1 , ( α − d ) / 2 ] 𝑧 1 𝛼 𝑑 2 z=\min[1,(\alpha-d)/2] is the dynamical exponent of Table 1 . In the Néel phase, the gap induces a short distance exponential decay, while the algebraically decaying interaction as 1 / r α 1 superscript 𝑟 𝛼 1/r^{\alpha} gives a long distance algebraic tail for arbitrary α 𝛼 \alpha . The correlation function ⟨ S i z ​ S j z ⟩ − ⟨ S i z ⟩ ​ ⟨ S j z ⟩ < 0 delimited-⟨⟩ superscript subscript 𝑆 𝑖 𝑧 superscript subscript 𝑆 𝑗 𝑧 delimited-⟨⟩ superscript subscript 𝑆 𝑖 𝑧 delimited-⟨⟩ superscript subscript 𝑆 𝑗 𝑧 0 \langle S_{i}^{z}S_{j}^{z}\rangle-\langle S_{i}^{z}\rangle\langle S_{j}^{z}\rangle<0 for i ≠ j 𝑖 𝑗 i\neq j in the XY phase, and shows a staggered pattern in the Néel phase.

Correlation functions	XY	Néel
$⟨S_i^{\nu }{S}_j^{\nu }⟩-⟨S_i^{\nu }⟩⟨S_j^{\nu }⟩$	$\sim 1/r^{{\eta }_{\nu }}$	$\sim e^{-r/\xi }+1/r^{{\eta }_{\nu }}$
${\eta }_x$	(long.) $2\left(d-z\right)$	(trans.) $\alpha$
${\eta }_y$	(trans.) $d-z$	(trans.) $\alpha$
${\eta }_z$	(trans.) $d+z$	(long.) $2\alpha$

Equations

$${\mathcal{H}}_{\mathrm{XXZ}}=\sum _{i\ne j}\frac{J_0}{r_{ij}^{\alpha }}\left[-(S_i^x{S}_j^x+S_i^y{S}_j^y)+\Delta S_i^z{S}_j^z\right]$$

$$\displaystyle{\cal H}_{\rm XXZ}=\sum_{i\neq j}~{}\frac{J_{0}}{r_{ij}^{\alpha}}~{}\Big{[}$$

$$(\Delta r_{ij}{)}_p=(L/\pi )\sin (r_{ij}^{(p)}\pi /L)r_{ij}^{(p)}=({\mathbf{r}}_i-{\mathbf{r}}_j)\cdot {\mathbf{e}}_p.$$

$${\Psi }_{\mathrm{FM}}$$

$${\Psi }_{\text{Neˊel}}$$

$${\Psi }_{\mathrm{XY}}$$

$$E_{\mathrm{FM}}$$

$$E_{\text{Neˊel}}$$

$$E_{\mathrm{XY}}$$

$${\Delta }_c(\alpha )=-\frac{{\sum }_{i\ne j}{r}_{ij}^{-\alpha }}{{\sum }_{i\ne j}{ϵ}_i{ϵ}_j{r}_{ij}^{-\alpha }}$$

$$S_i^x$$

$$S_i^y$$

$$S_i^z$$

$${\mathcal{H}}^{(2)}=\frac{1}{2}\sum _{\mathbf{k}}\left(\begin{array}{cc}{b}_{\mathbf{k}}^{†}& b_{\mathbf{k}}\end{array}\right)\left(\begin{array}{cc}{A}_{\mathbf{k}}& B_{\mathbf{k}}\\ B_{\mathbf{k}}& A_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{b}_{\mathbf{k}}\\ b_{\mathbf{k}}^{†}\end{array}\right)$$

$$A_{\mathbf{k}}={\gamma }_0+(\Delta -1){\gamma }_{\mathbf{k}}/2B_{\mathbf{k}}=(\Delta +1){\gamma }_{\mathbf{k}}/2$$

$${\gamma }_{\mathbf{k}}=\sum _{\mathbf{r}\ne 0}\frac{e^{i\mathbf{k}\cdot \mathbf{r}}}{r^{\alpha }}$$

$${\mathcal{H}}_{\mathrm{XY}}^{(2)}=\sum _{k\ne 0}{\gamma }_0\sqrt{1-\frac{{\gamma }_{\mathbf{k}}}{{\gamma }_0}}\sqrt{1+\Delta \frac{{\gamma }_{\mathbf{k}}}{{\gamma }_0}}{\beta }_{\mathbf{k}}^{†}{\beta }_{\mathbf{k}}+{\mathcal{H}}_{k=0}.$$

$$E_{\mathbf{k}}=\sqrt{A_{\mathbf{k}}^2-B_{\mathbf{k}}^2}={\gamma }_0\sqrt{1-{\gamma }_{\mathbf{k}}/{\gamma }_0}\sqrt{1+\Delta {\gamma }_{\mathbf{k}}/{\gamma }_0}.$$

$${\gamma }_{\mathbf{k}}={\gamma }_{\mathbf{k}}^{(\mathrm{nn})}+{\int }_{\rho >1}{d}^d\rho \frac{e^{ia\mathbf{k}\cdot \mathbf{\rho }}}{{\rho }^{\alpha }}$$

$${\gamma }_0-{\gamma }_{\mathbf{k}}$$

$$S_i^z$$

$$S_i^y$$

$$S_i^x$$

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

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

$${\mathcal{H}}_{\text{Neˊel}}^{(2)}=\sum _k\sqrt{(\Delta {\gamma }_{\mathbf{K}}+{\gamma }_{\mathbf{k}})(\Delta {\gamma }_{\mathbf{K}}+{\gamma }_{\mathbf{k}-\mathbf{K}})}{\beta }_{\mathbf{k}}^{†}{\beta }_{\mathbf{k}}.$$

$$H_{\mathrm{XXZ}}=\sum _{\mathbf{k}}{\gamma }_{\mathbf{k}}\left(-S_{\mathbf{k}}^x{S}_{-\mathbf{k}}^x-S_{\mathbf{k}}^y{S}_{-\mathbf{k}}^y+\Delta S_{\mathbf{k}}^z{S}_{-\mathbf{k}}^z\right)$$

$$H_{\mathrm{XXZ}}/{\gamma }_0\approx \frac{1}{N}\left[-({\mathbf{S}}_{\mathrm{tot}}{)}^2+(1+\Delta )(S_{\mathrm{tot}}^z{)}^2\right].$$

$$H_{\mathrm{XXZ}}|S,M⟩=\frac{{\gamma }_0}{N}\left[-S(S+1)+(1+\Delta )M^2\right]|S,M⟩.$$

$$|{\Psi }_0⟩=|S=N/2,M=0⟩$$

$$(E_{\mathrm{XY}}-E_0)/N=1/2N\to 0$$

$$⟨S_i^x{S}_j^x⟩-⟨S_i^x⟩⟨S_j^x⟩$$

$$⟨S_i^y{S}_j^y⟩$$

$$⟨S_i^z{S}_j^z⟩$$

$$⟨S_i^y{S}_j^y⟩$$

Full text

1. a short-range regime for $\alpha>d+\sigma_{c}$ (where $\sigma_{c}=1$ in the gapped Néel antiferromagnetic phase exhibited by the XXZ model, and $\sigma_{c}=2$ in the gapless XY ferromagnetic phase), sharing the same properties as those of finite-range interactions ($\alpha=\infty$); 2) a long-range regime $\alpha<d$, sharing the same properties as those of the infinite-range interactions ($\alpha=0$) in the thermodynamic limit; and 3) a most intriguing medium-range regime for $d<\alpha<d+\sigma_{c}$, continuously interpolating between the finite-range and the infinite-range behavior. The latter regime is characterized by elementary excitations with a long-wavelength dispersion relation $\omega\approx\Delta_{g}+ck^{z}$ in the gapped phase, and $\omega\sim k^{z}$ in the gapless phase, exhibiting a continuously varying dynamical exponent $z=(\alpha-d)/\sigma_{c}$. In the gapless phase of the model the $z$ exponent is found to control the scaling of fluctuations, the decay of correlations, and a universal sub-dominant term in the entanglement entropy, leading to a very rich palette of behaviors for ground-state quantum correlations beyond what is known for finite-range interactions.

I Introduction

Long-range (LR) interacting quantum many-body systems have attracted an increasing level of attention in the recent years. The experimental improvements in cooling, controlling and addressing few- to many-body atomic and molecular quantum systems possessing sizable LR interactions Lahayeetal2009 ; Schneideretal2012 ; Monroeetal2014 ; Browaeysetal2016 have triggered intense theoretical efforts, aimed at exploring the possibility that such interactions may stabilize stronger quantum collective phenomena with respect to the case of short-range interactions. In trapped-ion experiments, in particular, it has become possible to engineer Ising and exchange interactions between the spin of the ions decaying as a power-law $1/r^{\alpha}$ of the distance $r$, and with a continuously tunable exponent $\alpha$ ($0<\alpha<3$ brittonetal2012 ). Within this setup, experimentalists were able to observe how the dynamics of correlation spreading after a quantum quench is modified by the long-range interactions with respect to the case of ultracold neutral atoms interacting via a contact potential cheneauetal2012 . Within the context of ultracold neutral gases several groups have attained the quantum degeneracy of atoms possessing a large intrinsic magnetic moment, namely fermionic and bosonic isotopes of Cr griesmaieretal2005 ; dePazetal2013 ; nayloretal2015 , Dy luetal2011 ; luetal2012 and Er aikawaetal2012 ; aikawaetal2014 , and they were able to observe the coherent spin-exchange dynamics in these systems, induced by the large dipole-dipole ($1/r^{3}$) interaction dePazetal2013 ; baieretal2016 . Moreover infinite-range cavity-mediated interactions in a Bose-Einstein condensate have been experimentally demonstrated Baumannetal2010 ; Ritschetal2013 , leading to the spontaneous formation of long-range ordered phases (solid and supersolid). Finally, recent experimental progress in the manipulation of molecular systems with a large electric dipole hazzardetal2013 ; boetal2013 and of ensembles of Rydberg atoms schaussetal2015 ; barredoetal2015 ; Browaeysetal2016 has opened new pespectives for the quantum simulation of spin Hamiltonians with LR couplings.

In parallel to these remarkable experimental achievements, the theoretical efforts have focused on the study of equilibrium and out-of-equilibrium properties of LR interacting quantum lattice models. A large number of works has focused on the peculiar post-quench spreading of correlations eisertetal2013 ; haukeT2013 ; hazzardetal2014 ; schachenmayeretal2015 ; foss-feigetal2015 ; cevolanietal2015 ; maghrebietal2016 ; buyskikhetal2016 ; cevolanietal2016 and entanglement schachenmayeretal2013 ; haukeT2013 ; hazzardetal2014 ; buyskikhetal2016 in LR interacting systems, in connection with the breakdown (or generalizations) of Lieb-Robinson bounds liebR1972 constraining the dynamics eisertetal2013 ; foss-feigetal2015 . Further studies have rather focused on the ground-state properties of these models, where LR interactions can strongly affect the decay of correlation functions duttaB2001 ; laflorencieetal2005 ; peteretal2012 ; vodolaetal2014 , lead to phase transitions duttaB2001 ; capogrosso-sansoneetal2010 ; maghrebietal2015 ; gongetal2016 or modify substantially the entanglement properties vodolaetal2014 . In particular, the breakdown of the Mermin-Wagner theorem for sufficiently LR interactions leads to the possibility of spontaneous breaking of a continuous symmetry even in one-dimensional systems duttaB2001 ; maghrebietal2015 ; gongetal2016 .

In the face of the mounting body of experimental and theoretical results on LR interacting systems, it is of central importance to develop a broad (and possibly exhaustive) picture of the effect of LR interactions on the many-body physics of the system, and particularly so as the strength and decay law of interactions is varied continuously. While this endeavor might be too arduous to pursue for all the models of interest to experiment and theory in recent years, one can adopt a different strategy which focuses on a sufficiently simple model possessing nonetheless a rich phenomenology, and aims at extracting the most salient features and driving principles - in terms of excitations, fluctuations and entanglement - of the various regimes of LR interactions.

In this work we pursue the latter strategy by focusing our attention on the ground-state physics of a paradigmatic model in the theory of magnetism, the XXZ model for $S=1/2$, possessing power-law decaying ($1/r^{\alpha}$) isotropic ferromagnetic interactions for the $x$ and $y$ spin components, and ferro- or antiferromagnetic interactions for the $z$ spin components. The LR XXZ model describes also the physical situation of hardcore bosons with long-range hopping and density-density interactions. This model offers a rich showcase of effects of LR interactions on symmetry breaking phases: 1) in any physical dimension $d$ its phase diagram exhibits an extended ferromagnetic XY phase with breaking of the continuous rotation symmetry in the $xy$ plane, and the apparition of a gapless Goldstone mode. In particular in $d=1$ such a phase exists only thanks to the LR interactions, and it is stabilized by a decay exponent $\alpha<3$ against a Luttinger-liquid (namely gapless disordered) phase; 2) the above phase is competition with a Néel phase exhibiting long-range antiferromagnetic ordering along the $z$ axis and a gapped spectrum. As a consequence the LR XXZ model allows to monitor the effect of LR interactions with a continuously varying exponent $\alpha$ on phases breaking either a continuous or a discrete symmetry. In the XY phase the LR interactions are expected to stabilize against fluctuations the long-range order which already appears for nearest-neighbor interactions (in $d\geq 2$) . This fully justifies to treat such fluctuations as harmonic ones, as done by linear spin-wave (LSW) theory, on which we base most of our analysis. On the other hand the Néel phase is destabilized by the LR interactions due to their frustrated nature. Nonetheless in both phases, as well as at the transition between the two, the assumption of weak quantum fluctuations underlying LSW theory is well verified a posteriori, making our predictions quantitative. The case of antiferromagnetic XY interactions could also be treated in principle by LSW theory – but in this case the power-law decay induces frustration, which is expected to progressively weaken long-range order, making the LSW approach less and less justified as $\alpha$ decreases.

In particular LSW theory puts on the center stage the dispersion relation of free bosonic excitations, which is found in turn to control the scaling of fluctuations (or, equivalently, the decay of correlations), and the scaling of entanglement. In particular a systematic analysis of the long-wavelength properties of the dispersion relation allows to identify three distinct regimes upon varying the $\alpha$ exponent: 1) a short-range regime for $\alpha>d+\sigma_{c}$ ($\sigma_{c}=1$ in the Néel phase and $\sigma_{c}=2$ in the XY phase), whose main properties reproduce those observed in the limit $\alpha=\infty$ of finite-range interactions; 2) a medium-range regime for $d<\alpha<d+\sigma_{c}$, whose scaling properties (in the XY phase) are dominated by a continuously varying dynamical exponent $z=(\alpha-d)/\sigma_{c}$, governing the $k\to 0$ limit of the dispersion relation; 3) and a properly defined long-range regime for $\alpha<d$, possessing dispersionless excitations, and reproducing the properties of the infinite-range limit $\alpha=0$ in the thermodynamic limit. In particular this classification shows that the most interesting regime is the intermediate (medium-range) one, as it is the only one with markedly distinct features with respect to the two extreme limits of finite-range ($\alpha=\infty$) and infinite-range ($\alpha=0$) interactions which have been thoroughly investigated in the past. Importantly, the medium-range regime of the XY phase is also the one in which to frame the dipolar interaction in $d=$1, 2 and 3 (although $\alpha=3$ falls on the boundaries of the medium-range regime for $d=1$ and 3, and in its bulk only for $d=2$).

Despite their harmonic nature within LSW theory, ground-state quantum fluctuations exhibit a very rich structure in terms of scaling properties, and they are associated with an equally complex scaling of the entanglement entropy (EE) of a subsystem. A list of such properties in the three regimes of the LR XXZ model is presented in Table 1, which summarizes the most important results of the present work.

The structure of our paper is as follows: Sec. II introduces the model and is theoretical treatment; Secs. III and  IV discuss the structure of quantum fluctuations and entanglement in the ground state, respectively; conclusions are drawn in Sec. V.

II The Hamiltonian and its low energy properties.

In this Section, we introduce the model Hamiltonian under investigation, and proceed to determine its phase diagram within mean-field theory. We then investigate harmonic quantum fluctuations around the mean-field limit by diagonalizing the quadratic (LSW) Hamiltonian, and obtaining the dispersion relation of elementary excitations. The behavior of the dispersion relation in the limit of a vanishing wavevector $k\to 0$ is controlled by the Fourier transform of the power-lay-decaying interactions, and therefore it fundamentally depends on the decay exponent $\alpha$ as well as on the number of dimensions $d$. We further prove the self-consistency of the LSW approximation, showing that the LSW corrections to the mean-field solution are weak throughout the phase diagram, and even vanishing in the thermodynamic limit for $\alpha<d$.

II.1 Model Hamiltonian

The Hamiltonian of the XXZ model with LR interactions reads

[TABLE]

where the indices $i,j$ run over the nodes of a $N=L^{d}$ hyper-cubic lattice in $d$ spatial dimensions with periodic boundary conditions (unless otherwise specified). $S_{i}^{\beta}$ ($\beta=x,y,z$) are quantum-spin operators attached to each node. In the following we shall specify our discussion to the case of $S=1/2$. We consider ferromagnetic interactions for $xy$ spin components, while the $z$ spin components may have either ferromagnetic ($\Delta<0$) or antiferromagnetic ($\Delta>0$) interactions. The coupling $J_{0}$ simply sets the overall energy scale, and it will be set to one in the following.

The choice of $S=1/2$ and of ferromagnetic $xy$ couplings allows for a meaningful mapping of the above Hamiltonian onto that of hardcore bosons with long-range hopping and interaction

[TABLE]

where $a_{i}$, $a_{i}^{\dagger}$, and $n_{i}=a_{i}^{\dagger}a_{i}$ are hardcore-boson operators (satisfying the relations $\{a_{i},a_{i}^{\dagger}\}=1$ and $[a_{i},a_{j\neq i}^{(\dagger)}]=0$), related to the $S=1/2$ spin operators by the transformation $a_{i}=S_{i}^{-}$, $n_{i}-1/2=S_{i}^{z}$.

Unless explicitely stated, we shall work with periodic boundary conditions, and choose the finite-size regularization $r_{ij}^{2}=\sum_{p=1}^{d}(\Delta r_{ij})_{p}^{2}$, with

[TABLE]

The distance $r_{ij}$ is dimensionless, as it is measured in units of the lattice spacing; this choice leads to $r_{ij}^{-\alpha}\to\delta_{\langle ij\rangle}$ for $\alpha\to\infty$, where $\langle ij\rangle$ denotes a pair of nearest neighbors, namely the $\alpha\to\infty$ limit corresponds to finite-range interactions. On the opposite end, the limit $\alpha\to 0$ produces couplings with an infinite range. Traditionally, the distinction between short- and long-range interactions relies on the convergence properties of the sum $\sum_{j\neq i}1/r_{ij}^{\alpha}$ in the thermodynamic limit. As the latter diverges whenever $\alpha\leq d$, the separation between short- and long-range interactions is set at $\alpha=d$.

II.2 Mean-field phase diagram

In $d=2$ and $d=3$, the ground-state of the Hamiltonian Eq. (1) exhibits three different phases : an Ising ferromagnetic phase (FM) with spins aligned along the $z$ direction – corresponding to an insulating state of hardcore bosons with one particle per site; a Néel ordered phase with spins exhibiting a finite staggered magnetization along the $z$ axis – corresponding to a checkerboard solid of hardcore bosons; and an XY phase where the rotational symmetry about the $z$ axis is spontaneously broken, and spins order ferromagnetically along (for instance) the $x$ axis – corresponding to a superfluid condensate for the hardcore bosons.

The approximate location of the transition lines between the different phases is predicted by a simple mean-field argument, by comparing the energy of the corresponding mean-field states :

[TABLE]

where even (odd) sites are located at positions ${\bm{r}}_{i}$ such that $\sum_{p=1}^{d}\bm{r}_{i}\cdot{\bm{e}}_{p}$ is even (odd). The corresponding energies $E=\langle\Psi|{\cal H}_{\rm XXZ}|\Psi\rangle$ are :

[TABLE]

with $s=1/2$ in the present case. $\epsilon_{i}=1$ ($-1$) if $i$ is an even (odd) site. For $\Delta<-1$, the Ising ferromagnet (FM) has the lowest energy, and in fact it corresponds to the exact ground state: given that this state does not possess quantum correlations nor entanglement, we shall not discuss it any further, and we will restrict our attention to the case $\Delta>-1$. In the latter range, the condition $E_{\textnormal{N\'{e}el}}=E_{\rm XY}$ sets the transition line at

[TABLE]

Notice that, since $r_{ij}^{-\alpha}$ decays with distance between $i$ and $j$, the denominator is always negative, so that the above ratio is positive. For $\Delta>\Delta_{c}(\alpha)$, the system displays Néel order, while for $\Delta<\Delta_{c}(\alpha)$, the in-plane (XY) ferromagnetic order is favored. For $\alpha\to\infty$, the interactions are restricted to nearest neighbor, and one recovers the known result $\Delta_{c}(\infty)=1$. Due to the frustration of the interaction among $z$ spin components inherent to the long-range nature of the couplings, we expect that when $\alpha$ decreases, a larger value of $\Delta$ is necessary to stabilize the Néel order. And indeed, one can predict that at $\alpha<d$, the XY order is always stabler than the Néel order, since $\sum_{i\neq j}\epsilon_{i}\epsilon_{j}r_{ij}^{-\alpha}$ is always finite, while $\sum_{i\neq j}r_{ij}^{-\alpha}$ diverges. One thus obtains the phase diagram shown on Fig. 1. When including harmonic quantum corrections to the mean-field solution (as discussed in the following) the phase diagram preserves its topology in $d=2$ and 3. In $d=1$, on the other hand, the XY phase is destabilized for $\alpha>3$, in compliance with Mermin-Wagner theorem Auerbach-book , as further discussed in Sec. II.4 (the correct phase diagram for $d=1$ including quantum corrections to mean-field theory is to be found in Fig. 4).

II.3 Spin waves and excitation spectrum

In this section we study the low energy properties of the XXZ Hamiltonian through linear spin-wave theory, representing a semi-classical expansion around the mean-field ground-state colettaetal2012 .

II.3.1 XY phase

The mean-field ground-state in the XY phase is a perfect ferromagnet with all spins aligned along e.g. the $x$ axis. We then introduce Holstein-Primakoff (HP) boson operators $b_{i}$ which describe small deviation with respect to this perfect ferromagnet (namely the spin waves) holstein_primakoff ; colettaetal2012 :

[TABLE]

The LSW approximation consists in discarding all terms beyond quadratic in the HP transformation as well as in the resulting Hamiltonian. Its validity relies on the assumption that the quantum fluctuations of the spins are not strong enough to destroy the classical order. Technically, this hypothesis requires the populations of the HP bosons $\langle b_{i}^{\dagger}b_{i}\rangle$ to remain significantly smaller than 1/2 (or $s$ for a generic spin length $s$). A systematic comparison of some predictions of the LSW expansion with exact Monte-Carlo calculations in $d=2$ has been presented in Ref. colettaetal2012 , in the $\alpha=\infty$ (nearest-neighbor) case and at $\Delta=0$, demonstrating a remarkable accuracy for the LSW results. The accuracy of LSW predictions can only be expected to improve upon lowering the $\alpha$ exponent, as long as this does not entail crossing the phase transition separating the XY phase from the Néel phase, as discussed in the previous section. If instead one considers values of $\alpha$ and $\Delta$ close to the transition line in Fig. 1 one may expect the quality of LSW theory to worsen: nonetheless, as discussed in Sec. II.4 the first-order nature of the XY-Néel transition guarantees that our LSW treatment remains justified.

Taking advantage of the translational invariance, the Hamiltonian is conveniently expressed in momentum space. Introducing $b_{\bm{k}}=N^{-1/2}\sum_{j}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}}b_{j}$, and expanding the Hamiltonian up to second order in the $b_{\bm{k}}$ operators, one obtains

[TABLE]

where a constant term has been dropped. We have introduced the coefficients

[TABLE]

and the fundamental geometric coefficient

[TABLE]

corresponding to the Fourier transform of the coupling matrix.

The quadratic Hamiltonian ${\cal H}^{(2)}$ is diagonalized by a Bogoliubov transformation blaizot-ripka $b_{\bm{k}}=u_{\bm{k}}\beta_{\bm{k}}-v_{\bm{k}}\beta_{-\bm{k}}^{\dagger}$, with coefficients $u_{\bm{k}}=(1/\sqrt{2})(A_{\bm{k}}/\sqrt{A_{\bm{k}}^{2}-B_{\bm{k}}^{2}}+1)^{1/2}$ and $v_{\bm{k}}=\sqrt{u_{\bm{k}}^{2}-1}$, to lead to the form:

[TABLE]

In fact, the Bogoliubov transformation only applies to $k\neq 0$, while the $k=0$ sector deserves a special treatment. When $\alpha>d$, $\gamma_{0}=\sum_{r\neq 0}1/r^{\alpha}$ is non-divergent, and therefore the $k=0$ sector provides a microscopic contribution to the Hamiltonian and to the thermodynamic properties 333Note that $\gamma_{\bm{k}}$ converges for any nonzero $k$ and any $\alpha>0$.. For $\alpha>d$, ${\cal H}_{k=0}$ can then be safely ignored in the calculations when taking the thermodynamic limit. The $\beta_{\bm{k}}$ operators describe quasi-particles which, at the LSW level of approximation, correspond to exact eigenmodes of the many-body Hamiltonian. They represent collective fluctuations of the spins above the ground-state, with a dispersion relation

[TABLE]

The LSW approximation is dynamically stable if $E_{\bm{k}}$ is real for any $\bm{k}$. Since $\gamma_{\bm{k}}\leq\gamma_{0}$, the term $\sqrt{1-\gamma_{\bm{k}}/\gamma_{0}}$ is not problematic. As $\gamma_{\bm{k}}$ is maximally negative for $\bm{K}=(\pi,\pi,\dots)$ (see Fig. 2), we must require $\Delta<-\gamma_{0}/\gamma_{\bm{K}}$, coinciding with the condition $\Delta<\Delta_{c}(\alpha)$, Eq. (8), for the mean-field stability of the XY phase; this condition must be supplemented with $\Delta>-1$ to avoid the instability towards the Ising ferromagnetic phase. The resulting spectrum is gapless at $k=0$, corresponding to the Goldstone mode associated with the broken $U(1)$ rotational symmetry around the $z$ axis. An important qualitative difference with the short-range regime is that the usual linear dispersion $\omega\sim k$ for this Goldstone mode, is altered for $d<\alpha<d+2$ into $\omega\sim k^{z}$, with $z=(\alpha-d)/2$, as we shall discuss below. For $\alpha<d$, on the other hand, one observes that $\gamma_{0}\to\infty$, while $\gamma_{\bm{k}}$ remains convergent. As a consequence, the excitation spectrum (rescaled in units of $\gamma_{0}$) becomes dispersionless in the thermodynamic limit. We further discuss the case $\alpha<d$ in Section II.5.

Dynamical exponent.— The geometric coefficient $\gamma_{\bm{k}}$ controls the low energy dispersion relation $E_{\bm{k}}$ (see Eq. (13) and Eq. (19)), and the large distance decay of the spin correlations is also directly controlled by the small $k$ behavior of $\gamma_{\bm{k}}$ (see Appendix B). In order to obtain the correct scaling regimes, independently of the details of the lattice (as expected in the small-$k$ limit), and to recover simultaneously the correct dispersion relation in the limit $\alpha\to\infty$, we found convenient to treat exactly the nearest-neighbor contribution to $\gamma_{\bm{k}}$, and to approximate the rest of the sum by an integral. We thus have :

[TABLE]

where ${\bm{\rho}}={\bm{r}}/a$, $a$ is the lattice constant, and $\gamma_{\bm{k}}^{\rm(nn)}$ ($=2\sum_{i=1}^{d}\cos(k_{i}a)$ on the cubic lattice under the present investigation) is the value of $\gamma_{\bm{k}}$ for $\alpha\to\infty$, i.e. for nearest-neighbor interactions only. For any $\alpha>d$, $\gamma_{0}$ is expressed via a convergent integral. In Appendix A, we show that $\gamma_{\bm{k}}$ has the following behavior at small $k$ :

[TABLE]

As a consequence $\sqrt{\gamma_{0}-\gamma_{\bm{k}}}$ develops a cusp around $k=0$ for $d<\alpha<d+2$, turning into a divergence for $\alpha<d$ (see Fig. 2).

Interestingly, although the calculation of the integral in Eq. (15) requires a different treatment for each value of $d$, the low-$k$ dispersion shows a clear change at the simple $d$-dependent value of the $\alpha$ exponent, $\alpha=d+2$. When $\alpha>d+2$, the qualitative behavior of $\gamma_{0}-\gamma_{\bm{k}}$ is the same as for the short range limit $\alpha\to\infty$, with a prefactor of the $k^{2}$ scaling that depends on $\alpha$ and on the details of the lattice, and which diverges at $\alpha=d+2$.

On Fig. 2, we have plotted $\gamma_{\bm{k}}$ in $d=2$ along the $(k,k)$ diagonal of the Brillouin zone.

Excitation spectrum and sound velocity .— Given the expressions of the excitation spectrum $E_{\bm{k}}$ in the XY phase Eq. (13), and given the discussion of the previous paragraph (see Eq. (16)), the small $k$ behavior of the excitation spectrum is straightforwardly derived.

The result is summarized on Table 1, and illustrated on Fig. 3 at various representative points of the phase diagram. The spectrum is gapless in the whole XY phase, corresponding to the Goldstone mode of the broken U(1) rotational symmetry. Remarkably, the linear dispersion relation of this Goldstone mode, well known in the short range regime [Fig. 3(a)] is recovered only when $\alpha>d+2$, with an $\alpha$- and $\Delta$-dependent sound velocity $\partial E_{\bm{k}}/\partial k|_{k=0}$; while for $d<\alpha<d+2$, the dispersion relation behaves as $k^{z}$, with $z=(\alpha-d)/2<1$, and the group velocity diverges as $k^{z-1}$ at small $k$ [Fig. 3(d)].

II.3.2 Néel phase

In the Néel phase, the Holstein-Primakoff bosons represent small deviations around the reference mean-field state $\Psi_{\textnormal{N\'{e}el}}$, and they are introduced via the following transformation:

[TABLE]

with $\epsilon_{i}=1$ for even sites and $\epsilon_{i}=-1$ for odd sites. The reduction to an effective quadratic Hamiltonian goes along the same line as for the XY ferromagnet. The $A_{\bm{k}},B_{\bm{k}}$ coefficients take the form (with again $\bm{K}=(\pi,\pi,\dots)$)

[TABLE]

resulting in the following quadratic Hamiltonian

[TABLE]

The spectrum is gapped and twofold degenerate (because it is folded onto the smaller magnetic Brillouin zone, which is half the geometric one) with $E_{\bm{k}}=E_{{\bm{k}}-{\bm{K}}}$ (see Fig. 3). The stability requirement that the spectrum be real is equivalent to the condition $\Delta>-\gamma_{0}/\gamma_{\bm{K}}$, which is again the same criterion defining the phase boundary between the Néel phase and the XY ferromagnetic phase. Near $k=0$, the spectrum behaves as $E_{\bm{k}}\approx\Delta_{g}+ck^{z}$, where $\Delta_{g}=\gamma_{\bm{K}}\sqrt{\Delta+1}\sqrt{\Delta+\gamma_{0}/\gamma_{\bm{K}}}$, $c$ is some ($\alpha-$ and $\Delta-$dependent) constant and $z=\min(\alpha-d,2)$. The dispersion relation expected in the short-range regime is therefore recovered when $\alpha>d+2$ [Fig. 3(b)]. Moreover, for $d<\alpha<d+1$, the maximal group velocity (moving from a finite $k$ towards $k=0$) diverges as $v_{G}^{\max}=\max_{\bm{k}}|\bm{\nabla}_{\bm{k}}E_{\bm{k}}|\sim k^{\alpha-d-1}$ [Fig. 3(c)].

II.4 Stability of the spin-wave approximation

In this section, we briefly discuss the self-consistency of the LSW approximation in the different regimes. Self-consistency requires that the modification of the order parameter $m_{\rm SW}$, due to the non-zero population of Holstein-Primakoff bosons, remains small with respect to its classical value $m_{\rm cl}=1/2$. In the XY phase, the order parameter $m$ is the average magnetization along $x$, while it is the staggered magnetization in the Néel phase. According to the Holstein-Primakoff transformations Eq. (9) and (17), one thus has to evaluate $(m_{\rm cl}-m_{\rm SW})/m_{\rm cl}=2\langle b_{i}^{\dagger}b_{i}\rangle=(2/N)\sum_{\bm{k}}n_{\bm{k}}$, with $n_{\bm{k}}=\langle b_{\bm{k}}^{\dagger}b_{\bm{k}}\rangle=v_{\bm{k}}^{2}$. This integral is always finite in the gapped Néel phase, but could possibly diverge in the XY gapless phase due to the small $k$ behavior of $v_{\bm{k}}^{2}\sim 1/E_{\bm{k}}\sim k^{-z}$: this happens if $z\geq d$. Given that $z=1$ for $\alpha\geq d+2$ and $z<1$ otherwise, the only true instance of instability of the LSW approximation is found in $d=1$ for $\alpha>3$ in the XY phase, featuring a logarithmic divergence in the system size. The strongest renormalization of the order parameter due to quantum fluctuations is found for $\alpha\to\infty$ and $\Delta=1$, namely at the Heisenberg point of the XXZ model with nearest-neighbor interactions. Therefore in $d=2$ and $d=3$ the renormalization of the order parameter throughout the $\Delta$-$\alpha$ phase diagram is upped bounded by that of the nearest-neighbor Heisenberg model, namely 40% on the $d=2$ square lattice and 16% on the $d=3$ cubic lattice Mattis-book . This means in particular that the Néel-to-XY quantum phase transition at $\alpha<\infty$ does not possess stronger quantum renormalization than those at the Heisenberg point, and therefore it is still quantitatively described by LSW theory. This is intrinsically due to the first-order nature of the phase transition, which implies that true quantum critical fluctuations and entanglement do not develop.

In $d=1$, the renormalization of the classical order diverges in the thermodynamic limit when $\alpha\to 3^{-}$ – a finite-size system the calculation still delivers a finite renormalization, e.g. for $\alpha=3$ and $L=10^{5}$ the renormalization reaches $86\%$. The breakdown of LSW theory for $\alpha\geq\alpha_{c}=3$ signals a phase transition, absent at the mean-field level, between a long-range-ordered XY ferromagnet (or a superfluid condensate) for $\alpha<\alpha_{c}$ and a quasi-long-range-ordered Luttinger liquid (LL) phase Giamarchi-book for $\alpha>\alpha_{c}$. The phase diagram in $d=1$ is reported in Fig. 4. This quite unusual phase transition, specific to $d=1$, has been studied in more details in maghrebietal2015 ; gongetal2016 beyond the harmonic approximation, while LSW theory can only indicate the existence of such a transition, but it cannot quantitatively describe it. In particular, as shown in Ref. maghrebietal2015 the phase transition occurs at a $\Delta$-dependent value $\alpha_{c}(\Delta)<3$. In view of the limitations of LSW results in this regime, in Sec. IV.3.2 we shall complement them with fully quantitative ones based on the DMRG approach, in order to investigate the evolution of entanglement properties across the XY-LL transition.

II.5 Exact ground state for $\alpha<d$

When $0<\alpha<d$, LSW theory predicts that classical order is not renormalized by quantum fluctuations, irrespective of the dimension $d$ and of the precise value of $\alpha$. Indeed, in this long-range regime – which for $\Delta>-1$ always falls into the XY phase – $\gamma_{0}$ diverges, while $\gamma_{\bm{k}}$ goes to a constant for any fixed nonzero $k$ in the thermodynamic limit. This implies that $A_{\bm{k}}/B_{\bm{k}}\to 0$ at any $\bm{k}\neq 0$, so that $u_{\bm{k}}\to 1$ and $v_{\bm{k}}\to 0$ in the thermodynamic limit. In other words, the spin waves $b_{\bm{k}}^{\dagger}|0\rangle$ become exact eigenstates of the Hamiltonian, whose ground state is the mean-field trial state. One is thus left with an effective dynamics for the ground state which occurs only in the $k=0$ sector (the collective spin ${\bm{S}}_{\rm tot}=\sum_{i}{\bm{S}}_{i}$). One might question the reliability of the LSW approach to draw a definite conclusion on the nature of the ground state. However, we can show that the prediction of LSW theory is essentially exact: the XY ferromagnetic state, while not being the exact ground state on a finite size system, has a vanishingly small energy density above the ground state in the thermodynamic limit. In order to do so, we first rewrite the Hamiltonian as :

[TABLE]

where we have introduced $S_{\bm{k}}^{\beta}=N^{-1/2}\sum_{i}e^{i{\bm{k}}\cdot{\bm{r}}_{i}}S_{i}^{\beta}$. Since $\gamma_{k\neq 0}/\gamma_{0}\to 0$ in the thermodynamic limit, we may keep only the $k=0$ sector, which reproduces the so-called Lipkin-Meshkov-Glick model Latorreetal2005 :

[TABLE]

The latter model is readily diagonalized by the total-spin basis $|S,M\rangle$ to give:

[TABLE]

As long as $\Delta>-1$, the ground-state is

[TABLE]

with an energy $E_{0}=-N/4-1/2$ (since $\gamma_{0}$ diverges for $N\to\infty$, we measure energy in units of $\gamma_{0}$, so that the thermodynamic limit is well defined). The state $|\Psi_{0}\rangle$ – often referred to in the atomic-physics literature as the spin-squeezed state – is obtained from the application of $(\sum_{i}\sigma_{i}^{+})^{N/2}$ to the state with all spins down, where $\sigma_{i}^{+}$ flips the spin at position $i$ if it is a spin down, and gives zero if applied to a spin up. Being the symmetric, equal weight superposition of all possible spin configurations with $N/2$ spins up and $N/2$ spins down, it contains no further correlation between the spins beyond the global constraint that the total spin along $z$ be zero. And, in the thermodynamic limit, the mean-field XY ferromagnet (which does not possess any form of correlation), has a vanishingly small energy density compared to the squeezed state:

[TABLE]

and it has essentially the same properties as the true ground state. Another way to understand this result, following the early works of Anderson Anderson1952 ; Anderson-book , is to realize that the mean-field XY ferromagnet also lies in the $S=N/2$ manifold, but it is a superposition of many states with different values of $M$. Those states differ in energy with respect to the ground state by $E(N/2,M)-E_{0}=(\Delta+1)(M^{2}/N)$ which scales to 0 as $1/N$ — building up the famous “tower of states”. Any superposition of states in the $S=N/2$ sector sufficiently peaked around $M=0$ thus becomes degenerate with the true ground state in the thermodynamic limit. This latter point of view only partly applies to the XY ferromagnet, since the distribution of $M$ has, in that case, a width of order $O(\sqrt{N})$, and hence an excitation energy of order ${\cal O}(1)$, as we have shown above. If the behavior on finite-size systems can strongly deviate from the mean-field one, one can nevertheless conclude that any correlation effect for $\alpha<d$ is a finite size effect. The entanglement properties in the long-range regime will be further discussed in Sec. IV.1

III Structure of quantum correlations and fluctuations in the ground state

III.1 Decay of correlation functions

In this Section, we discuss the LSW predictions for the large-distance decay of the spin correlations in the ground-state. We focus on $\alpha>d$, since all correlations vanish in the thermodynamic limit for $\alpha<d$ – they just stem from the conservation of $S^{z}=\sum_{i}S^{z}_{i}$ (see Section II.5). Our observations are summarized in Table 2, and the calculations leading to the various spin correlations are detailed in Appendix B.

III.1.1 XY phase

In the XY phase, the rotational symmetry about the $z$ axis is spontaneously broken. As a result, the spin correlations in the $xy$ plane are not isotropic : the $S^{y}S^{y}$ spin correlations (transverse to the order parameter), do not possess the same decay as the $S^{x}S^{x}$ spin correlations (longitudinal to the order parameter). The $S^{z}S^{z}$ spin correlations exhibit a third distinct decay behavior. Fig. 5 shows the various spin correlations in the XY phase for $d=2$.

LSW predicts the following behaviors:

[TABLE]

An analytical understanding for these decay exponents can be obtained from the small-$k$ behavior of the structure factor $S_{\bm{k}}^{\beta\beta}=N^{-1}\sum_{i,j}e^{i\bm{k}\cdot(\bm{r}_{j}-\bm{r}_{i})}\langle S_{i}^{\beta}S_{j}^{\beta}\rangle$ ($\beta=x,y,z$). In Appendix B, we show indeed that $S_{\bm{k}}^{zz}\sim k^{z}$ while $S_{\bm{k}}^{yy}\sim k^{-z}$; then the decay exponent of the correlation function can be related to that of the small-$k$ of the corresponding structure factor via the calculation of Appendix A , assuming algebraically decaying correlations. Moreover the fact that the decay exponent of the $S^{x}S^{x}$ correlation function is twice the one of the $S^{y}S^{y}$ correlation function can be understood in that $S^{x}S^{x}$ is quartic in the HP bosons, while $S^{y}S^{y}$ is quadratic (see Appendix B). As the large-distance decay of the spin-spin correlations is controlled by the $\alpha$- and $d$-dependent $z$ exponent, we can conclude that the same decay exhibited for finite-range interactions ($\alpha=\infty$) is recovered over the whole short-range regime $\alpha>d+2$; while in the medium-range regime ($d<\alpha<d+2$) the decay exponents continuously depend on $\alpha$ exhibiting a broad variety of different correlation regimes. Note that the exact ground state state on a finite-size system does not break the U(1) symmetry; and the LSW prediction for the $\eta_{y}$ exponent is the one which reproduces the power-law decay of correlations in the symmetric case.

III.1.2 Néel phase

In the Néel phase, the rotational symmetry about the $z$ axis is preserved, so that the correlations in the $xy$ plane are fully isotropic – a feature reproduced by the LSW approach, as we show in Appendix B. As illustrated in Fig. 5, the spin correlations show an hybrid decay (exponential at short distance, followed by an algebraic decay at long distance):

[TABLE]

with $a$, $a^{\prime}$ some ($\alpha$- and $\Delta$- dependent) constants. The simple exponential decay expected in the finite-range regime is thus recovered only for $\alpha\to\infty$ in a rather nontrivial way, while for any $\alpha<\infty$ the spatial decay of interactions dictates the long-distance spatial decay of correlations. Similar findings were reported earlier on other gapped systems in Refs. peteretal2012 ; vodolaetal2014 . The hybrid decay (exponential followed by a power law) is a common feature of gapped systems with long-range interactions – it has been observed, for instance, in the dynamics following a quench gongetal2014 , and in the ground state of a spin-1 topological phase gongetal2016b .

III.2 Fluctuations of the collective spin

To complement the analysis of the large distance decay of the spin correlations, we offer in this Section a scaling analysis of the associated fluctuations of the collective spin. Considering the collective spin component $S^{\beta}=\sum_{i}S^{\beta}_{i}$ with $\beta=x,y,z$, one has

[TABLE]

(where we have introduced the symbol $\delta O=O-\langle O\rangle$), namely the fluctuations are composed of a local term and of a correlation term. If $\left[S^{\beta},\cal H\right]\neq 0$ - namely for the non-conserved spin components $\beta=x,y$ of the collective spin, the ground state of the Hamiltonian has finite global fluctuations of the collective spin component in question, whose scaling can be readily estimated from the knowledge of the power-law decay of the correlation function $\langle\delta S_{i}^{\beta}\delta S_{j}^{\beta}\rangle\sim r_{ij}^{-\eta_{\beta}}$. Indeed, from Eq. (30) one readily deduces that:

[TABLE]

Clearly, if $d-1-\eta_{\beta}<-1$ ($\eta_{\beta}>d$) , then the scaling of fluctuations is ${\cal O}(L^{d})$, while $d-1-\eta_{\beta}>-1$ ($\eta_{\beta}<d$) will provide a correction to the conventional volume law of fluctuation scaling, namely

[TABLE]

In the Néel phase, $\eta_{x,y}=\alpha>d$, so that corrections to volume scaling of fluctuations do not appear.

In the XY phase, on the other hand, given the results presented in Table 2, the whole medium-range regime is characterized by a violation of the volume law of fluctuations for the $S^{y}$ collective spin, $\langle(\delta S^{y})^{2}\rangle\sim L^{d+z}$, while the longitudinal fluctuations of the order parameter $\langle(\delta S^{x})^{2}\rangle\sim L^{\max(d,2z)}$ could violate a volume law only if $z\geq d/2$, implying $2d<\alpha<d+2$ – a condition which is only satisfied in $d=1$ for $2<\alpha<3$ (see Fig. 6(a)).

When considering instead the $S^{z}$ spin component, which is conserved by the XXZ Hamiltonian, one has that $\langle(\delta S^{z})^{2}\rangle=0$, so that the only meaningful question to ask concerns the scaling of the fluctuations of the $S^{z}_{A}=\sum_{i\in A}S_{i}^{z}$ operator defined on subsystem $A$. In particular the fluctuations of the local $S^{z}_{A}$ spin component can be expressed entirely in terms of correlations between $A$ and $B$, namely

[TABLE]

The scaling of $\langle\delta^{2}S^{z}_{A}\rangle$ is then fundamentally governed by the decay of correlations. Considering that the system is defined on a hyper-torus, and considering an equal $A-B$ bipartition of the system (namely $L_{A},L_{B}\sim L$), one can show that

[TABLE]

Details on the calculation are provided in Appendix C, showing that further logarithmic corrections are expected if $\eta_{z}=d$ or $\eta_{z}=d-1$. Hence in the Néel phase, where $\eta_{z}=2\alpha>2d\geq d+1$, an area-law scaling of $S^{z}$ fluctuations is verified over the whole medium- and short-range regimes in all dimensions $d\geq 1$.

In the XY phase, on the other hand, $\eta_{z}=d+z$, so that the whole of the short-range regime with $z=1$ exhibits a logarithmically violated area law for the $S^{z}$ fluctuations, as observed in the finite-range limit bipartite-fluct ; whereas the medium-range regime with $z<1$ exhibits a continuous violation of the area law, $\langle\delta^{2}S^{z}_{A}\rangle\sim L^{d-z}$, up to a fully developed volume-law scaling in the long-range regime.

Our findings for the scaling of fluctuations are summarized in Table 1. Fig. 6 shows the predicted scaling exponents for the collective-spin fluctuations, in the XY phase compared with the numerical reconstruction of the exponents via direct calculations on systems with sizes $L=10^{4}\dots 10^{5}$ in $d=1$ and $L=10\dots 200$ in $d=2$, showing that finite-size effects may be significant, yet they do not prevent from observing the strong $\alpha$ dependence of the scaling of fluctuations in the medium-range regime.

IV Scaling of the subsystem entanglement entropy

In this section we investigate the entanglement properties of the ground state, focusing on the scaling of the entanglement entropy of a subsystem. The latter is defined starting from the reduced density matrix (RDM) $\rho_{A}$ of subsystem $A$ as a partial trace of the ground-state projector over $B$ (namely the complement of $A$):

[TABLE]

The entanglement entropy (EE) is then defined as the von Neumann entropy of the RDM, $S_{A}=-{\rm Tr}\rho_{A}\log\rho_{A}$. It is useful to recall that any quantum state admits a so-called Schmidt decomposition Peres-book related to the $A$-$B$ bipartition of the system, in the form:

[TABLE]

where the states $|\psi_{A}^{(p)}\rangle$ (resp. $|\phi_{B}^{(p)}\rangle$) form an orthogonal basis of the Hilbert space $H_{A}$ (resp. $H_{B}$) of subsystem $A$ (resp. $B$). The EE is then simply obtained as the Shannon entropy of the Schmidt coefficients $\lambda_{p}$ $S_{A}=-\sum_{p}\lambda_{p}\log\lambda_{p}$.

LSW theory allows for a very efficient calculation of the entanglement properties of the ground state, thanks to the Gaussian form of the RDM of any subsystem. We refer the reader to Sec. IV.2 and to Refs. song ; Luitzetal2015 ; FrerotR2015 , where the calculation of the EE for quadratic bosonic systems is detailed.

IV.1 Entanglement and fluctuations in the long-range regime

When $\alpha<d$ the Schmidt decomposition of the exact ground state in the thermodynamic limit can be determined explicitly, allowing for an exact calculation of the EE Latorreetal2005 . We have argued in Section II.5 that for $\alpha<d$ the ground state is, in the thermodynamic limit, $|\Psi_{0}\rangle=|S=N/2,~{}M=0\rangle$. Introducing $|S=N/2,~{}M=-N/2+p\rangle$ for the symmetric, equal-weight superposition of all states having $p$ spins up and $N-p$ spins down, the Schmidt decomposition of $|\Psi_{0}\rangle\propto(\sum_{i}\sigma_{i}^{+})^{N/2}|S_{A}=\frac{N_{A}}{2},~{}M_{A}=-\frac{N_{A}}{2}\rangle\otimes|S_{B}=\frac{N_{B}}{2},~{}M_{B}=-\frac{N_{B}}{2}\rangle$ for a bipartition into $N_{A}\leq N/2$ and $N_{B}=N-N_{A}$ spins is then simply

[TABLE]

where $\lambda_{p}=\binom{N_{A}}{p}\binom{N_{B}}{N/2-p}/\binom{N}{N/2}$ counts the number of ways (normalized to unity) in which to distribute $p$ up-spins up among $N_{A}$ spins, and $N/2-p$ up-spins among $N_{B}$ spins — the remaining $N/2$ spins being down-spins. The distribution of the subsystem magnetization $M_{A}=S_{A}^{z}$ is also given by $\lambda_{p}$, and it is centered around $S_{A}^{z}=0$ with a width of order $\sqrt{N_{A}}$, hence a variance scaling as $N_{A}\sim L_{A}^{d}$, consistent with the scaling $L_{A}^{d-z}$ and $z=0$ (compare Table 1).

The EE is simply the entropy of the $\lambda_{p}$ distribution, which scales as $(d/2)\log L_{A}$ —indeed, the entropy is of order $\log\Omega$, with $\Omega\sim\sqrt{N_{A}}$ the typical number of nonzero $\lambda_{p}$’s. In this long-range regime, the Schmidt basis $|\psi_{A}^{(i)}\rangle$ is contained in a small subspace of the local Hilbert space (namely that of symmetric superposition states with a fixed magnetization), having dimensions $O(N_{A})$ to be contrasted with the local Hilbert space dimensions ($2^{N_{A}}$). The example at hand highlights the existence of a simple relation between the structure of ground-state entanglement and that of the fluctuations of the subsystem magnetization $S_{A}^{z}$, the latter being the only physical mechanism responsible for entanglement. Even though such a simple relationship cannot be found in the short- and medium-range regimes, one can always say that, if $S_{z}$ is globally conserved in the system, its fluctuations on subsystem $A$ are a sufficient condition for the existence of $AB$ entanglement FrerotR2015 .

IV.2 Entanglement and fluctuations in the medium- and short-range regime

In this section, we discuss the scaling of the EE of a subsystem in the medium- and short-range regime as obtained via LSW theory, contrasting it with the scaling of the fluctuations of the subsystem magnetization $\langle\delta^{2}S_{A}^{z}\rangle$. Our findings concerning the scaling of EE are summarized in Table 1. In the gapped Néel phase, both quantities are found to obey and area law RMP-area-laws , namely to scale as the boundary of $A$, in both regimes. In the XY phase, on the other hand, there is a stark contrast between the scaling of entanglement and that of magnetization fluctuations. Indeed entanglement is always found to obey an area law in both medium- and short-range regimes, while fluctuations exhibit systematic violations of the area-law scaling, as discussed in Sec. III.2. In particular the contrast is very stark in the medium-range regime, in which the prefactor of the area-law scaling of entanglement vanishes progressively as $\alpha$ is reduced towards $d$ (see Fig. 7(b)), while the scaling of $S_{A}^{z}$ fluctuations moves from an area-law one at $\alpha=d+2$ to a volume-law one at $\alpha=d$ [see Fig. 7(a)].

In Fig. 7(b) we fitted the EE as

[TABLE]

with an additional subdominant logarithmic contribution $b\log L$, and a constant term $c$. Fig. 7(d) shows that the area-law coefficient decreases monotonically to zero when $\alpha$ decreases towards $d$. The $b$ coefficient, shown in Fig. 7(c), can be attributed to the Goldstone modes associated with the broken rotational symmetry metlitskiG2011 , and it will be discussed in Section IV.3.1, while the special case of $d=1$ (possessing a logarithmic scaling even in the absence of long-range order) will be analyzed in Sec. IV.3.2.

The strong decrease of the EE upon decreasing $\alpha$ can be understood from the structure of the single-particle entanglement spectrum, and from its comparison to the (bulk) physical spectrum. The single-particle entanglement spectrum is defined starting from the Gaussian structure of the RDM $\rho_{A}$ describing the $A$ subsystem, which in turn is a consequence of the harmonic approximation made within LSW theory FrerotR2015 :

[TABLE]

where ${\cal H}_{A}$ is the so-called entanglement Hamiltonian, ${\bm{b}}=(b_{1},b_{2},...,b_{N_{A}})$ is the vector of Bose operators on the $A$ subsystem, and $h_{A}$ is the single-particle entanglement Hamiltonian. The diagonalization of $h_{A}$ FrerotR2015 brings the quadratic entanglement Hamiltonian to the form ${\cal H}_{A}=\sum_{\bm{k}_{||},n}\tilde{E}_{\bm{k}_{||},n}\beta_{\bm{k}_{||},n}^{\dagger}\beta_{\bm{k}_{||},n}$, where $\bm{k}_{||}$ is the momentum along the $A$-$B$ cut – representing a good quantum number for the half-torus geometry of subsystem $A$ that we adopt here – and $n$ is a further mode index related to the dynamics generated by ${\cal H}_{A}$ transversely to the cut.

The Gaussian RDM $\rho_{A}$ describes a thermal gas of quasiparticles with dispersion relation given by $\tilde{E}_{\bm{k}_{||},n}$ at unit temperature - whose thermal entropy is identical to the EE. Fig. 8 shows that the low-lying modes of the single-particle entanglement spectrum $\tilde{E}_{\bm{k}_{||},n}$ become increasingly stiff as $\alpha$ decreases – in striking analogy with the physical spectrum of spin-wave excitations, whose low-energy density of states behaves as $\rho(\omega)\sim\omega^{d/z-1}$, becoming smaller when $\alpha$ (and $z=(\alpha-d)/2$) decrease. A progressive reduction in the low-energy density of states of the single-particle entanglement spectrum $\tilde{E}_{\bm{k}_{||},n}$ obviously implies a reduction of entropy at fixed temperature for $\rho_{A}$, namely of EE. Conversely, the EE is maximal along the phase-transition line between the XY and Néel phase, as shown on Fig. 9, and this increase can be again associated with the $\tilde{E}_{\bm{k}_{||},n}$ spectrum, which, at the transition, acquires a second zero mode at $\bm{k}_{||}=(\pi,\pi,...)$, besides the one already present at $\bm{k}_{||}=0$, as shown in Fig. 10. The appearance of a second soft mode obviously boosts the low-energy density of states, leading to an increase of the EE FrerotR2016 . Correspondingly, the spin-wave spectrum (also shown in Fig. 10) shows the same appearance of a second soft mode (at $\bm{K}=(\pi,\pi,...)$), which signals the instability of the XY ground state to the appearance of long-range staggered spin order as in the Néel phase. The corresponding increase in the low-energy density of states of the spin-wave excitations implies in turn a stronger quantum correction of the classical ordered moment: hence the simultaneous softening of the spin-wave dispersion relation and of the single-particle entanglement spectrum at the transition is responsible for the striking similarity between the strong enhancement of entanglement (Fig. 9) and of quantum fluctuations of the order parameter (Fig. 1) around the XY-Néel transition.

Finally, Fig. 11 shows the variance of $S_{A}^{z}$ and of $S_{A}^{y}$ across the phase diagram (for a fixed finite size). The evolution of spin fluctuations along the $y$ axis strongly resembles that of the entanglement shown in Fig. 9, although it more clearly reveals the first-order nature of the transition with a jump in the fluctuations properties. The jump is related to a sudden change in the scaling properties of the fluctuations, as detailed in Table 1. A jump is also to be observed in the EE, albeit much weaker as the dominant scaling behavior of entanglement is the same area law on both sides of the transition.

IV.3 Medium- and short-range regime in the XY phase: sub-leading contribution to the entanglement entropy

We now take a closer look at the details of the entanglement scaling in the XY phase, focusing on the observed logarithmic corrections to the area law scaling of EE - as reported in the previous section. We first recall in Section IV.3.1 the origin of this logarithmic contribution in the presence of a broken symmetry, as stemming from the interplay between the tower-of-state entanglement spectrum of a subsystem and the low-lying Goldstone modes coupling two subsystems. In so doing we shall rephrase arguments which have been put forward in earlier works metlitskiG2011 ; kulchytskyyetal2015 , but we will also generalize them to the case of long-range interactions, which add the new ingredient of a continuously varying dynamical exponent $z\leq 1$. Furthermore we shall specialize in Section IV.3.2 our discussion to the case of $d=1$, where the analysis is somewhat more subtle. There the logarithmic correction to the area law becomes the dominant term in the medium-range regime, and interestingly a similar logarithmic scaling persists also in the short-range regime, featuring a Luttinger-liquid (LL) phase with conformal invariance: supplementing our LSW calculations with numerical DMRG data will allow us to study the evolution of entanglement across the XY-LL transition.

IV.3.1 XY phase: logarithmic term from the tower-of-state spectrum

A universal additive logarithmic contribution to the area-law scaling of EE has been predicted to appear in systems breaking a continuous symmetry in their ground state in Ref. metlitskiG2011 , and extensively verified numerically in Refs.Luitzetal2015 ; kulchytskyyetal2015 ; FrerotR2015 . Such a contribution can be traced back to the existence of a low-lying tower of states (ToS) in the entanglement spectrum of a subsystem $A$, akin to the low-lying spectrum of a finite-size system exhibiting spontaneous symmetry breaking in the thermodynamic limit Anderson-book . If isolated from its complement $B$, subsystem $A$ would indeed possess such a ToS in the low-lying excitation spectrum, described by the effective Hamiltonian $H^{\rm(ToS)}={\cal L}^{2}/2I$. Here $\cal L$ stands for the generator of the symmetry spontaneously broken in the thermodynamic limit – in the case of XY symmetry at hand ${\cal L}=S_{A}^{z}$, the generator of rotations in the $xy$ plane. The ToS Hamiltonian describes therefore the angular momentum of a rigid rotor living on a $N$-dimensional sphere, where $N$ is the number of components of the order parameter ($N=2$ in our case), and possessing a moment of inertia $I\sim L_{A}^{d}$, scaling as the volume of the subsystem. The exact same observations apply to subsystem $B$.

The coupling between subsystems $A$ and $B$ creates entanglement between the $A$ and $B$ rigid rotors in the joint ground state. As a consequence of the $A$-$B$ coupling the system possesses $N_{G}=N-1$ Goldstone modes spreading coherently across the two subsystems. In particular the characteristic energy scale for the coupling of the $A$ and $B$ rigid rotors is set by the smallest-wavevector Goldstone mode(s) with $k\sim 1/L$, to which we attribute an energy $\Delta_{G}=\Delta_{G}(L)$. As discussed in Ref. metlitskiG2011 ; kulchytskyyetal2015 , the Hamiltionan of the coupled $A$ and $B$ rotors can then be approximated as that of a harmonic oscillator of frequency $\Delta_{G}/\hbar$ metlitskiG2011 ; kulchytskyyetal2015 . Tracing out subsystem $B$ leads to a density-matrix description of subsystem $A$, in which ToS modes are populated up to an energy of the order of $\Delta_{G}$, namely up to an angular momentum ${\cal L}_{\max}\sim(I\Delta_{G})^{1/2}$; hence, knowing that ToS levels ${\cal L}$ have a degeneracy of order ${\cal L}^{N-2}$, the EE can be estimated by simple state counting as $S_{A}\sim\log\Omega_{A}$, where

[TABLE]

As a consequence the EE receives a contribution from the ToS spectrum of the kind

[TABLE]

As discussed at length in Sec. II.3 in the long-range XXZ model of interest here $\Delta_{G}\sim L^{-z}$, leading then to the result

[TABLE]

The above result generalizes to the case $z\neq 1$ the universal logarithmic term in the EE of a continuous-symmetry-breaking phase first obtained in Ref. metlitskiG2011 .

Within LSW theory there is no guarantee that this additive logarithmic contribution can be accurately captured, especially in view of the fact that the LSW approach cannot describe the ToS sector of the spectrum, as it assumes an explicitly broken symmetry even on finite-size systems. Nevertheless, in the case of linearly dispersing ($z=1$) Goldstone modes it has been shown song ; FrerotR2015 that a careful treatment of the $k=0$ sector of the LSW Hamiltonian enables one to actually capture the universal logarithmic term within LSW theory. This is possible by gapping out the $k=0$ LSW mode with the addition of a small term $-BS_{\rm tot}^{x}$ in the Hamiltonian which stabilizes the ferromagnetic order and introduces a gap of order $\sqrt{B}$ at $k=0$. The choice $B=1/L^{2d}$ allows therefore to mimic the $1/L^{d}$ scaling of the ToS spectrum in the thermodynamic limit. As shown on Fig. 7(b,c), this procedure leads to the appearance of a logarithmic term $b\log L_{A}$ in the EE extracted from the LSW data, whose $b$ prefactor reproduces rather closely the predicted universal form $b=(d-z)/2$. The deviations can be in part attributed to finite-size effects, which, for a given maximum size $L$ used in the fits of the entanglement scaling, become more significant the lower $\alpha$, justifying the modest agreement, observed at low $\alpha$ in Fig. 7(c), between finite-size LSW calculations and the universal prediction.

IV.3.2 From Luttinger liquid behavior to the XY phase in $d=1$

As already commented above, the case $d=1$ is rather special in that the system possesses two ground-state phases in the XY regime upon varying $\alpha$: an XY phase with long-range order for $\alpha<3$, and a Luttinger-liquid (LL) phase with algebraic correlations for $\alpha>3$. In the latter regime the low-energy physics of the system is captured by conformal field theory (CFT) calabreseC2004 , predicting a universal logarithmic violation of the entanglement area law in the form $S_{A}=(nc/6)\log L_{A}$, where $c=1$ is the central charge of the CFT, and $n$ ($=1,2$) is the number of common boundaries between subsystem $A$ and its complement $B$. On the other hand, for $\alpha<3$ conformal invariance breaks down (as $z<1$) and the system develops long-range order. The EE still possesses a universal logarithmic term, as discussed in the previous section, that is independent of the geometry of the $A$-$B$ boundary. It becomes therefore extremely interesting to explore the evolution of the EE across the XY-LL transition to track how the CFT prediction and the ToS contribution evolve into each other.

LSW theory is obviously inadequate to the latter scope, as it breaks down for $\alpha\geq 3$ (and it becomes increasingly inaccurate as $\alpha\to 3^{-}$). Therefore we complement our approach with ground-state DMRG calculations whitedmrg ; white93 , a method that gives easily access to the EE of a partition of the system. We perform the numerical simulations on lattices up to $L=50$ sites implementing both open (OBC) and periodic boundary conditions (PBC) and setting $\Delta=0$. Since the Hamiltonian Eq. (1) commutes with total magnetization we restrict our study to zero magnetization sector. In order to keep a truncation error smaller than $10^{-7}$, throughout the simulations we used up to $M=300$ states. We chose a hopping amplitude of the form $1/|r_{i}-r_{j}|^{\alpha}$ for both OBC and PBC, both in the DMRG and the LSW calculations.

The scaling of half-chain EE is shown in Fig. 12 for OPC and PBC in panels (a) and (b) respectively for different values of $\alpha$ ranging from $0.1$ to $5$. A general form for the dominant term of the EE scaling reads

[TABLE]

where one expects $b_{\rm CFT}=c/6$ and $b_{\rm ToS}=0$ for $\alpha>3$, $b_{\rm CFT}\neq 0$ and $b_{\rm ToS}=(d-z)/2=(3-\alpha)/4$ for $1<\alpha<3$, and $b_{\rm CFT}=0$, $b_{\rm ToS}=1/2$ for $\alpha<1$. The meaning of $b_{\rm CFT}$ in the $1<\alpha<3$ (which lacks conformal invariance) is simply that of the prefactor to a logarithmic term which depends on the number of boundaries $n$. The ability of both LSW theory 444The introduction of OBC in LSW theory requires to diagonalize the quadratic Hamiltonian in real space - a well-established procedure which we implement following Ref. blaizot-ripka . and DMRG to simulate both OBC ($n=1$) and PBC ($n=2$) allows to systematically extract finite-size estimates of the $b$ coefficients as

[TABLE]

The coefficients so extracted are shown in Fig. 12: there we find a relatively good agreement between the theoretical expectations - valid in the infinite-size limit - and the finite-size numerical data, coming both from the quadratic LSW approximation, and from exact DMRG calculations. Finite-size effects appear to be very significant (and especially so at small $\alpha$). Yet the application of a transverse field $B\sim 1/L^{2}$ allows to use LSW theory on finite sizes even in the regime $\alpha>3$, and it is remarkable to observe that even in this regime LSW theory on finite system sizes remains predictive when compared with DMRG.

V Conclusions and perspectives

Inspired by the possibilities of quantum simulation in ultracold atoms, we studied the low energy properties of the $d$-dimensional XXZ Hamiltonian with ferromagnetic couplings in the $xy$ plane and both ferro- and antiferromagnetic along the $z$ axis, both decaying as a power law $1/r^{\alpha}$ with the distance. Linear spin-wave theory proved to be a reliable tool to determine the ground state phase diagram, the low energy excitation spectrum as well as the correlation and entanglement properties in the various phases exhibited by the system: a gapped antiferromagnetic Néel ordered phase, a gapless ferromagnetic XY phase, and a gapped ferromagnetic Ising phase.

In the gapped Néel phase (present only for $\alpha>d$), we identified two fundamentally different regimes of interactions: a short-range regime ($\alpha>d+1$), and a medium-range regime ($d<\alpha<d+1$). Both regimes are similar in most aspects to the $\alpha\to\infty$ (nearest-neighbor) limit, with the main qualitative difference that the spin-spin correlations exhibits a short-distance exponential decay controlled by the gap $\Delta_{g}$ followed by a long-distance power-law decay of the correlations in the ground state, controlled directedly by the exponent $\alpha$. Furthermore, in the medium-range regime ($d<\alpha<d+1$) the short-range regime dispersion relation $E_{k}\sim\Delta_{g}+ck^{2}$ acquires a cusp $E_{k}\sim\Delta_{g}+ck^{\alpha-d}$ at small $k$, while the short-range dispersion relation is recovered at $\alpha>d+2$.

In the XY phase, our calculations identify three fundamentally distinct regimes of interactions: a short-range regime ($\alpha>d+2$) characterized by a dynamical exponent $z=1$ akin to the finite-range limit ($\alpha=\infty$), a medium-range regime ($d<\alpha<d+2$) with a $\alpha$-dependent dynamical exponent $z=(\alpha-d)/2$, and a long-range regime ($\alpha<d$) exhibiting the same properties as the infinite-range regime ($\alpha=0$) in the thermodynamic limit, namely correlations uniquely stemming from finite-size effects. In all regimes the dynamical exponent $z$ is found to control directly the long-distance decay of the spin correlations and the scaling of fluctuations of the collective spin in the ground state. In particular, the medium-range regime of the XY phase exhibits a continuously varying palette of scalings for the collective-spin fluctuations, all of them violating the scalings exhibited in the conventional short-range regime. In particular the local fluctuations of the conserved collective spin ($S^{z}$) on a subsystem exhibit algebraic corrections to the area law, as a result of the coherent exchange of magnetization quanta between two subsystems mediated by the non-local couplings.

The scaling of fluctuations is in sharp contrast with the behavior of entanglement entropy, which in the XY phase is found to always satisfy an area law (plus additive logarithmic corrections), and with an area-law prefactor which generically decreases when the range of interactions increases. This behavior can be traced back to the density of states in the quasi-particle entanglement spectrum (as well as in the physical spectrum), which is lowered upon decreasing $\alpha$ as the dispersion relation of quasi-particle excitations stiffens. Finally, an additive logarithmic correction to the area-law scaling is found in the XY phase, associated with the spontaneous breaking of the continuous U(1) symmetry. The prefactor to the logarithmic term is universal, uniquely depending on the $z$ exponent and the number of components of the order-parameter; our prediction generalizes in a non-trivial manner that of Ref. metlitskiG2011 for short-range interacting systems, and it can be quantitatively tested against accurate density-matrix renormalization group results in $d=1$.

Beside their intrinsic theoretical interest, our results have a direct relevance for ongoing experiments on ultracold atoms in optical lattices. Indeed the Hamiltonian (1) with $\alpha=3$ in $d=2$ is realized by the Mott insulating phase of magnetic atoms dePazetal2013 when imposing the conservation of the magnetization along the quantization axis. Furthermore its $d=1$ implementation can be envisioned in trapped ions PorrasC2004 ; Monroeetal2014 , which also enable to vary continuously the $\alpha$ exponent of the power-law decay of interactions. In particular the existence of a medium-range regime exhibiting an $\alpha$-dependent, “curved” dispersion relation of elementary excitations, as well as an $\alpha$-dependent scaling of fluctuations, represent a prediction which lend itself rather naturally to an experimental test. In the case $\alpha=3$, generalizing the model to antiferromagnetic interactions in the $xy$ plane (while taking $\Delta=0$) would mimic the physics of lattice-trapped Rydberg atoms with resonant interactions barredoetal2015 ; Browaeysetal2016 . The XX model with nearest-neighbor ($\alpha=\infty$) interactions on a bipartite lattice has the same physics irrespective of the sign of the interactions. Therefore, in the case $\alpha>d$ - for which nearest-neighbor interactions dominate the energetics of the system - we may expect that much of the physics observed in the ferromagnetic case carries over to the antiferromagnetic one, as long as one restricts to bipartite lattices. The exploration of antiferromagnetic $xy$ interactions on bipartite lattices, as well as on non-bipartite ones, represents an exciting extension of our present work.

Our work shows that long-range interactions offer a very rich landscape in terms of scaling properties of ground-state quantum fluctuations and subsystem entanglement entropy. Establishing a quantitative link between the two is necessary to connect entanglement to directly measurable properties, yet it appears rather challenging. A strategy we shall pursue in the future is based on the concept of “local entanglement thermodynamics” swingleM2016 , postulating an explicit Ansatz for the entanglement Hamiltonian in the form of the original microscopic Hamiltonian of the system, yet with spatially modulated coupling constants. Finally, the existence of sharply distinct regimes for the dispersion relation of elementary excitations in the model of interest suggests that the non-equilibrium unitary dynamics following a quench will also be extremely rich, revealing an unconventional spreading of correlations which directly reflects the existence of a continuously varying dynamical exponent $z$ Frerotetal2017 .

VI Acknowledgements

We thank D. Vodola for fruitful discussions. PN acknowledges the Institut Universitaire de France and the INFN grant QUANTUM for financial support. This work was supported by ANR (”ArtiQ” project).

Appendix A Fourier transform of the interaction potential

In this Section, we prove the scaling behavior of Eq. (16). To this goal, we analyze the small-$k$ behavior of the integral :

[TABLE]

for $\alpha>d$. More precisely, we will analyse $\gamma_{0}-\gamma_{\bm{k}}$ at small $k$. Since $\gamma_{0}^{(\rm n.n)}-\gamma_{\bm{k}}^{(\rm n.n)}=2\sum_{i=1}^{d}[1-\cos(k_{i}a)]\sim(ka)^{2}$ on the hypercubic lattice we consider, $\gamma_{0}-\gamma_{\bm{k}}$ will always scale at least as fast as $(ka)^{2}$ at small $k$. Note that, although we can obtain the correct scaling for $\gamma_{0}-\gamma_{\bm{k}}$, we are not able to predict the correct prefactors, as they depend on the details of the lattice which are not captured by our continuous approximation. Finally, we note that if $\alpha<d$:

[TABLE]

A.1 $d=1$

In $d=1$, one has to evaluate :

[TABLE]

As $\alpha>1$, the integral diverges at $k\to 0$. One thus integrates by parts to obtain :

[TABLE]

Now, since $1-\alpha<0$, the first term on the r.h.s is $\cos(ak)(ak)^{1-\alpha}/(\alpha-1)$. Secondly, since $x^{1-\alpha}\sin x\sim x^{2-\alpha}$ for $x\to 0$, the integral on the r.h.s converges to a finite value $C$ when $1<\alpha<3$ (note that it diverges for any $\alpha<1$ and it is not defined for $\alpha=1$). Therefore we can establish that, for $1<\alpha<3$ :

[TABLE]

for some constant $C$ (which depends on $\alpha$). Since $2/(\alpha-1)$ is the value of the integral on the l.h.s at $k=0$, and since $\alpha-1<2$, the dominant scaling at small $k$ in $d=1$ is :

[TABLE]

If $\alpha>3$, the integral on the r.h.s of Eq. (48) diverges at $k\to 0$ as $k^{3-\alpha}$, so that one obtains :

[TABLE]

for some constant $C^{\prime}$. Hence we can prove that in $d=1$ :

[TABLE]

A.2 $d=2$

In $d=2$, the integral to be calculated is :

[TABLE]

where ${\cal J}_{n}(x)=(1/2\pi)\int_{0}^{2\pi}e^{i(x\cos\theta-n\theta)}d\theta$ is a Bessel function. As ${\cal J}_{0}(0)=1$, the integral on the r.h.s is divergent for $\alpha>2$, and we perform an integration by parts. Using the fact that ${\cal J}_{0}^{\prime}(x)=-{\cal J}_{1}(x)$, one has :

[TABLE]

The first term on the r.h.s gives $(ak)^{2-\alpha}/(2-\alpha)+O(k^{4-\alpha})$, while, as ${\cal J}_{1}(x)\sim x$ at small $x$, the integral on the r.h.s converges when $2<\alpha<4$ and diverges for $\alpha>4$ (note that as ${\cal J}_{1}(x)\sim\sqrt{2/\pi x}\cos(x-3\pi/4)$ for $x\to\infty$, the integral converges in $+\infty$). Thus, we have proved that, for $2<\alpha<4$ :

[TABLE]

for some constant $C$. As $2\pi/(2-\alpha)$ is the value of the integral on the l.h.s of Eq. (53) at $k=0$, we have thus established that in $d=2$ :

[TABLE]

If $\alpha>4$, the integral on the r.h.s of Eq. (54) diverges as $k^{4-\alpha}$ when $k\to 0$, so that

[TABLE]

for some constant $C^{\prime}$. We have thus established that, in $d=2$ :

[TABLE]

A.3 $d=3$

In $d=3$, we have to evaluate :

[TABLE]

The integral on the r.h.s diverges in $k=0$ for $\alpha>3$, and we integrate it by parts. The remaining integral shows the same divergence in $k=0$, so that a second integration by parts is needed. One obtains :

[TABLE]

The first term on the r.h.s gives the $k=0$ value of the integral to evaluate in Eq. (59). If $3<\alpha<5$, the integral on the r.h.s or Eq. (A.3) converges, so that one obtains the scaling behavior in $d=3$ :

[TABLE]

If $5<\alpha$, the integral on the r.h.s of Eq. (A.3) diverges in $k\to 0$ as $k^{5-\alpha}$, thus the result :

[TABLE]

Appendix B Structure factor of the spin correlations

In this Appendix, we detail all the calculations of the spin correlations in the ground state at the level of the spin-wave approximation. The starting point is the Holstein-Primakoff (HP) representation of the spin operators (we consider a general spin $s$) :

[TABLE]

in the XY ferromagnetic phase (where the reference mean-field state has all spins pointing towards $+X$). In the Néel phase (where the reference state has spins pointing alternatingly towards $\pm Z$), we have :

[TABLE]

where $\epsilon_{i}=+1$ on even sites, and $\epsilon_{i}=-1$ on odd sites. As the LSW approximate ground state, being the ground state of a quadratic Hamiltonian, satisfies Wick’s theorem blaizot-ripka , all the properties of the quantum fluctuations in the ground state are contained in the two-point correlation functions $\langle b_{i}^{\dagger}b_{j}\rangle$ and $\langle b_{i}b_{j}\rangle$. Higher order correlations can be expressed, through Wick’s theorem, as a function of two-point ones blaizot-ripka . If one assumes translational invariance, the correlations are most naturally expressed in momentum space. As the Bogoliubov rotation reads $b_{\bm{k}}=u_{\bm{k}}\beta_{\bm{k}}-v_{\bm{k}}\beta_{-\bm{k}}^{\dagger}$, and given the form of the Hamiltonian Eq. (10) in the main text, one obtains :

[TABLE]

where $v_{\bm{k}}^{2}=(1/2)(A_{\bm{k}}/E_{\bm{k}}-1)$ and $u_{\bm{k}}v_{\bm{k}}=B_{\bm{k}}/2E_{\bm{k}}$. Here, $E_{\bm{k}}=\sqrt{A_{\bm{k}}^{2}-B_{\bm{k}}^{2}}$, and $A_{\bm{k}}$ and $B_{\bm{k}}$ take the expressions in Eqs. (11) and (18).

Notice that none of the correlation functions possesses an imaginary part. An equivalent expression for the correlations in Fourier space is then:

[TABLE]

The other correlations vanish.

XY phase.— In the XY phase, the rotational invariance about the $z$ axis is spontaneously broken. As a consequence, the $S^{y}S^{y}$ and $S^{x}S^{x}$ spin correlations need not have the same structure. At the level of the LSW approach, the $S^{x}S^{x}$ spin correlation is indeed of order 4 in the HP bosons operators, while the $S^{y}S^{y}$ correlation is only of order 2. One may thus expect that the decay exponent of the $S^{x}S^{x}$ spin correlations is twice the decay exponent of the $S^{y}S^{y}$ spin correlations, and this expectation is indeed confirmed numerically. The calculation for the $S^{x}S^{x}$ spin correlations is better formulated in real space :

[TABLE]

where on the last line, we have used Wick’s theorem.

On the other hand, the $S^{y}S^{y}$ spin correlations read (notice that $\langle S_{i}^{y}\rangle=0$) :

[TABLE]

where we have used the fact that the correlation functions are real. One can directly deduce from this formula the expression of the static spin structure factor for the $y$ spin components:

[TABLE]

Finally, the $S^{z}S^{z}$ correlations are:

[TABLE]

and we obtain the structure factor for the $z$ spin components as:

[TABLE]

As $\sqrt{1-\gamma_{\bm{k}}/\gamma_{0}}\sim k^{z}$ at small $k$, we thus obtain the scaling behavior for the structure factors in the XY phase :

[TABLE]

In correspondence with the small-$k$ behavior of the structure factors, we find that the associated correlations decay at large distance as :

[TABLE]

Indeed, if one assumes that correlations decay algebraically with distance, this behavior follows from the calculations of Appendix A, where we proved that

[TABLE]

as long as $-d<\sigma<2$, and FT stands for Fourier transform.

On the other hand, we verified numerically that the $S^{x}S^{x}$ spin correlations decay algebraically with a decay exponent twice as large as the decay exponent of the $S^{y}S^{y}$ correlations, as expected in view of the fact that the $S^{x}S^{x}$ correlation is quartic in the HP bosons, while $S^{y}S^{y}$ is quadratic:

[TABLE]

Néel phase.— In the Néel phase, the rotational symmetry about the $z$ axis is not broken, so that one expects $S_{\bm{k}}^{xx}=S_{\bm{k}}^{yy}$. Given the formula for the HP transformation, the expression of $S_{\bm{k}}^{xx}$ in the Néel phase is exactly the same as the expression of $S_{\bm{k}}^{zz}$ in the XY phase, namely :

[TABLE]

On the other hand, the calculation for the $S^{y}S^{y}$ spin correlations gives (with ${\bm{K}}=(\pi,\pi,\dots)$) :

[TABLE]

from which we can deduce the structure factor :

[TABLE]

Given the expression of $A_{\bm{k}}$ and $B_{\bm{k}}$ in the Néel phase, we thus have $S_{\bm{k}}^{xx}=S_{\bm{k}}^{yy}$, as expected. The structure factor scales as the dispersion relation at small momentum:

[TABLE]

From the expression of the structure factor, one can then reconstruct the spatial decay of the spin correlations. We find an exponential decay at short distance (associated with the gap $\Delta$), and an algebraic decay at large distance, with an exponent $\alpha$, directly controlled by the decay of the spin-spin interaction. Finally, the $S^{z}S^{z}$ spin correlations read:

[TABLE]

which also decays as $1/|\bm{r}_{i}-\bm{r}_{j}|^{\alpha}$ at large separation, and show a staggered pattern.

Appendix C Scaling of the covariance

In this Appendix, we provide details on the calculation leading to the prediction of Eq. (34). Considering two subsystems $A$ and $B$ of linear size $L$ and two observables ${\cal O}_{A/B}=\sum_{i\in A/B}{\cal O}_{i}$, the covariance of ${\cal O}_{A}$ and ${\cal O}_{B}$ is given by

[TABLE]

Assuming that the correlation function decays as a power-law with distance :

[TABLE]

we can relate the scaling behavior of $\langle\delta{\cal O}_{A}\delta{\cal O}_{B}\rangle$ to the exponent $\eta$. For the sake of mathematical simplicity, we take $A$ and $B$ translationally invariant in $d-1$ dimensions (they are thus lines in $d=1$ and cylinders in $d=2$, while in $d=3$, the geometry is more difficult to visualize, since the $A-B$ contact area has the topology of a torus. As our aim is simply to predict the scaling of the covariance with the linear size of $A$ and $B$, their precise shape is irrelevant. The analysis proceeds in two steps 1) we analyse the decay of the correlation contour FrerotR2015

[TABLE]

where $i$ is a site in $A$, whose distance from the $A-B$ boundary is denoted $a$, as illustrated on Fig. 13;

and 2) we deduce the scaling of $\langle\delta{\cal O}_{A}\delta{\cal O}_{B}\rangle=\sum_{i\in A}C_{i}$ from the decay of $C_{i}$ when moving away from the $A-B$ boundary.

Decay of the contour.— We have to evaluate the scaling behavior of

[TABLE]

The last integral converges if $\eta\geq d-1$ and is of order $[(L/x)/(1+u)]^{d-1-\eta}$ otherwise (with logarithmic corrections if $\eta=d-1$). In this latter case ($\eta<d-1$), we find that $C(x)\sim L^{d-\eta}$. If $\eta>d-1$, we have that

[TABLE]

If $\eta>d$, the integral converges and $C(x)\sim x^{d-\eta}$. If $\eta<d$, the integral is of order $(L/x)^{d-\eta}$ and we find again that $C(x)\sim L^{d-\eta}$. We also predict logarithmic corrections if $\eta=d$. To sum up, we have the following behavior for the contour

[TABLE]

with logarithmic corrections if $\eta=d$ or $\eta=d-1$.

Scaling of the covariance.— Summing the contour $C_{i}$ over all sites in $A$ gives a trivial factor of $L^{d-1}$ for the $d-1$ directions transverse to $x$. Integrating the distance $x$ to the $A-B$ boundary gives the final result of Eq. (34) with further logarithmic corrections if $\eta=d$ or $\eta=d-1$.