Long-range Ising and Kitaev Models: Phases, Correlations and Edge Modes

We analyze the quantum phases, correlation functions and edge modes for a class of spin-1/2 and fermionic models related to the 1D Ising chain in the presence of a transverse field. These models are the Ising chain with anti-ferromagnetic long-range interactions that decay with distance $r$ as $1/r^\alpha$, as well as a related class of fermionic Hamiltonians that generalise the Kitaev chain, where both the hopping and pairing terms are long-range and their relative strength can be varied. For these models, we provide the phase diagram for all exponents $\alpha$, based on an analysis of the entanglement entropy, the decay of correlation functions, and the edge modes in the case of open chains. We demonstrate that violations of the area law can occur for $\alpha \lesssim1$, while connected correlation functions can decay with a hybrid exponential and power-law behaviour, with a power that is $\alpha$-dependent. Interestingly, for the fermionic models we provide an exact analytical derivation for the decay of the correlation functions at every $\alpha$. Along the critical lines, for all models breaking of conformal symmetry is argued at low enough $\alpha$. For the fermionic models we show that the edge modes, massless for $\alpha \gtrsim 1$, can acquire a mass for $\alpha<1$. The mass of these modes can be tuned by varying the relative strength of the kinetic and pairing terms in the Hamiltonian. Interestingly, for the Ising chain a similar edge localization appears for the first and second excited states on the paramagnetic side of the phase diagram, where edge modes are not expected. We argue that, at least for the fermionic chains, these massive states correspond to the appearance of new phases, notably approached via quantum phase transitions without mass gap closure. Finally, we discuss the possibility to detect some of these effects in experiments with cold trapped ions.


Introduction
Topological superconductors and insulators have generated enormous interest in recent years as they correspond to examples of novel quantum phases that are not captured by the familiar Ginzburg-Landau theory of phase transitions. Breakthrough experiments have already led to the observation of symmetry protected topological phases both in condensed-matter systems [1] and atomic, molecular, and optical physics [2,3]. While topological phases are finding applications in fields as diverse as photonics and spintronics, the recent probable observation of Majorana modes [4][5][6][7][8][9][10] in solid-state materials represents the first major step towards the realization of topological quantum computing.
Majorana modes are non-dispersive states with zero energy. In [11], Kitaev has shown that these modes can exist localized at the edges of a one-dimensional superconductor made of spinless fermions with short-range (SR) p-wave pairing. This model is solvable and the underlying lattice Hamiltonian can be mapped exactly onto the well-known Ising chain in a transverse field in one-dimension. For SR interactions, the latter is a text-book example of Hamiltonian displaying a quantum phase transition, here from an ordered (anti-)ferromagnetic phase to a disordered paramagnetic one. Following earlier theoretical works [12,13], recent experiments with cold trapped ions have generated enormous interest by demonstrating that long-range (LR) Ising-type Hamiltonians arise as the effective description for the dynamics of the internal states of cold trapped ions, acting as (pseudo-)spins with two or, recently three, internal states. In these experiments, effective spin interactions are generated by a laser-induced manipulation of the vibrational modes of the ion chain [13][14][15][16][17], which are naturally long ranged. The resulting Ising-type interactions are antiferromagnetic and decay with distance r as a power-law a r 1 , with an adjustable exponent α usually in the range   a 0 3.5. In experiments with cold ions, the quantum state of individual particles can be prepared and measured. As a result, both the static and dynamical properties of the many-body system are accessible. Recent experiments have led to the observation of instances of interaction-induced frustration [18], non-local propagation of correlations [19][20][21][22][23] and entanglement in a quantum many-body system [16,24]. Very recently, spectroscopy experiments have focused on the precise determination of the excited states of LR models [25].
The experimental works described above are based on the understanding of the phase diagram of the Isingchain in a transverse field, which is known exactly for SR interactions only. In a seminal work [26], Koffel, Lewenstein and Tagliacozzo have explored the phase diagram of this system with LR interactions in the parameter range  a 0.5. The results were intriguing: (i) the connected correlation functions decay with a power-law tail even within the gapped paramagnetic phase, at odds with conventional wisdom inherited from SR models and consistent with earlier predictions for other quantum models with LR interactions [12,[21][22][23]27]. Crucially, (ii) the entanglement entropy, usually a constant within gapped phases, seems to scale logarithmically with the system size within the paramagnetic phase for sufficiently small  a 1 as well as sublogarithmically for a > 1. This is remarkable as would signal a violation of the so-called 'area law', dictating the behavior of the entanglement entropy in SR quantum mechanical systems. These studies also confirm (iii) the persistence of antiferromagnetic and paramagnetic orders with decreasing α.
Research in the area of topological phases with LR interactions is very active, and several possible experimental realizations have been recently proposed. In particular, Kitaev chains with non-local hopping and pairing may be realized in solid state architectures with so-called helical Shiba chains, made of magnetic impurities on an s-wave superconductor [28,29]. For atomic and molecular systems, key implementations of topological phases have been proposed with polar molecules, dipolar ground state atomic quantum gases and Rydberg excited atoms [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. In addition, the famous Haldane phase may be soon realized in cold ion experiments with three internal states per ion [48,49], simulating spin-1 particles. For this latter model, very recent theoretical work [50] has demonstrated that major features of symmetry-protected topological order can persist for LR interactions. It remains an open question to determine the validity of these results for generic symmetry-protected topological phases with LR interactions.
In this work, we analyze the quantum phases, correlation functions and edge-mode localization of a class of spin-1/2 and fermionic models related to the one-dimensional Ising chain in the presence of a transverse field. These models are the Ising chain with anti-ferromagnetic LR interactions, as discussed in [26], as well as a class of Hamiltonians corresponding to a generalization of the Kitaev chain, where both the hopping and pairing terms are LR with an algebraic decay a r 1 , and their relative strength can be varied. For these models, we provide the phase diagram for all exponents α, based on an analysis of: (i) the entanglement entropy; (ii) the decay of correlation functions in all phases; (iii) the mass and the localization properties of the edge modes when the chains are open.
In the case of the long-range Ising (LRI) chain we utilize numerical calculations based on the density-matrixrenormalization group (DMRG) method [51,52], while the long-range Kitaev-type (LRK) models remain exactly solvable for all α, allowing for analytical calculation.
The following results are obtained for all models.
(i) A violation of the area law for the entanglement entropy occurs in gapped regions with  a 1. For  a 1, no violation is found.
(ii) For any finite α, connected correlation functions within the gapped phases display a hybrid decay that is exponential at short distances and algebraic at long ones. The power of the algebraic decay, however, as well as the extension of the two decay regimes, depends on α. However, when  a 1, the connected correlation functions show a purely algebraic decay.
(iii) For the LRK models, we provide an exact analytical expression for the decay of correlation functions within the gapped phases that describes the hybrid behavior with distance mentioned above and explains its origin.
(iv) Along the critical lines, we demonstrate that conformal symmetry is broken for sufficiently small α, by analyzing the finite size scalings of the von Neumann entropy and of the energy density for the ground states, as well as the behavior of the dynamical exponent with α around the minima of the spectrum.
(v) We find the existence of two kinds of edge modes: massless and massive. For the LRK models, massless (Majorana) modes, as previously found in [53], appear in the antiferromagnetic region of the phase diagram for large α. The antiferromagnetic phase for the LRK models is defined in analogy with that of the shortrange Ising chain. The massive modes, instead, are entirely new and are found to appear in a broad area of the antiferromagnetic phase for  a 1 and when we choose the unbalance ò between the strengths of the hopping and pairing terms to be different from1. This results suggests, for  a 1, a restoration of the  2 symmetry associated with the (absence of) ground state degeneracy, as well as a possible transition to a novel  2 symmetric phase and without mass gap closure. Interestingly, if we choose  = 1, the massless modes survive for all a > 0 and they are exponentially localized at the edge of the system.
(vi) For the LRI chain in the antiferromagnetic phase, edge modes are massless for all α, and, up to numerical precision are exponentially localized at the edge of the chain, in contrast, e.g. to [50]. However, surprisingly we find in the paramagnetic phase alocalization of excited, gapped, energy eigenstates for  a 1, which for  a 1 are instead delocalized in the bulk.
(vii) We finally discuss the persistence of some of the LR effects discussed above (e.g., hybrid decay of correlation functions and edge mode localization) in small chains of up to 30 sites, as relevant to current experiments.
The paper is organized as follows. In section 2 we introduce the model Hamiltonians that we consider in this work (section 2.1), the observables that are used to characterize the various phases (section 2.2), and present the corresponding phase diagrams (section 2.3). In particular, in section 2.4 we discuss the critical lines of the Ising and Kitaev models, and argue that conformal symmetry is broken for sufficiently small α. In section 3 we provide an analytic calculation of the correlation functions for the fermionic Hamiltonians that explains the hybrid exponential and algebraic decay observed in these LR models (sections 3.1 and 3.2). In section 3.4 we provide a numerical comparison with results for the LRI chain, displaying similar behavior. In section 4 we analyze the edge modes in the LRI and LRK chains. In particular, in section 4.1 we analyze the properties of gapless Majorana modes that are found in the anti-ferromagnetic phases of the LRK models for  a 1. In section 4.2, instead, we demonstrate that the edge modes can become massive for  a 1, signalling a transition to a new phase. We discuss similar results obtained for some excited states that get localized on the edges in the paramagnetic phase of the LRI model for  a 1. In section 5 we discuss the observability of some of the results above in small chains of a few tens of particles, as relevant for cold ions experiments. Finally, section 6 discusses conclusions and outlook.

Model Hamiltonians and quantum phases
In this section we introduce the model Hamiltonians that we consider in this work and present the corresponding phase diagrams that we compute based on results from the entanglement entropy, decay of correlation functions, spectrum of excitations, and edge mode localization, as discussed in detail in the following sections. where s n j (n = x y z , , ) are Pauli matrices for a spin-1/2 at site j on a chain of length L. The first term on the right-hand side of equation (1) describes spin-spin interactions that we choose antiferromagnetic (AM) with q > sin 0 (or equivalently q p < < 0 ). The second term describes the coupling of individual spins to an external field pointing in the z-direction. Thus, while the first term favours an antiferromagnetic phase with spins pointing along the x direction, the second terms favours a paramagnetic (PM) phase where all spins align along z. In the case of SR interactions (i.e., for a  ¥) the Hamiltonian equation (1) is exactly solvable and a quantum phase transition between these two phases is known to occur at q p = 4 c . Reference [26] has shown numerically that a quantum phase transition separating the AM and the PM phases survives for all finite  a 0.5. Below, we are interested in exploring the phase diagram of equation (1) for all α and θ. where † a j describes the creation operator for a fermionic particle at site j, and = † n a a j j j . The Hamiltonians equation (2) represent generalizations of the Kitaev chain for spinless fermions with superconducting p-wave pairing, where both the hopping and pairing terms decay with distance algebraically with exponentα. Here, all energies are expressed in dimensionless units and the parameter ò governs the unbalance between the hopping and pairing terms.
In the SR limit a  ¥, equation (

von Neumann entropy
Entanglement measures are routinely used to characterize the critical properties of strongly correlated quantum many-body systems [55]. A key example is the von Neumann entropy  ℓ that we employ in this work. For a system of L sites that is partitioned into two disjoint subsystems A and B containing ℓ and L−ℓ sites, respectively,  ℓ is defined as Tr log , 5 2 where r ℓ is the reduced density matrix of the subsystem A. Two general behaviors of  ℓ are known for the ground states of one-dimensional SR interacting systems. Within gapped phases,  ℓ saturates to a constant value independently of ℓ and thus obeys to the so-called area law [56]. On the contrary,  ℓ diverges with ℓ for critical gapless phases and, for conformally invariant systems, satisfies the universal law [57]: 2 s i n for the case of open boundaries. Here, a is a non-universal constant and c is the central charge of the theory. The latter characterizes the universality class of the gapless phase [58,59].  ℓ , and thus the central charge, can in principle be directly computed by numerical techniques, such as density matrix renormalization group [51,52] as well as by means of analytical methods for quadratic Hamiltonians [60][61][62][63].
In the case of LR models it has been shown [26,53,56] that the divergence of  ℓ in equation (6) can also occur for gapped phases, corresponding to a so-called violation of the area law. Since this violation is found to be logarithmic, an effective central charge c eff may be defined also within the gapped phases and used to characterize the main features of the phase diagram for LR interactions [26].

Correlation functions
The various quantum phases can be characterized by the decay of two-point correlation functions with distance. For the LRI model, we are interested in the connected correlations x in the LRI model, respectively, via the Jordan-Wigner transformation given above [54]. Since the LRK models are quadratic, the functions equations (8) and(9) can be directly obtained from the one-point correlations † a a i j and † † a a i j via Wick's theorem. In particular, one finds For SR interactions, the connected correlations above are known to decay exponentially (algebraically) with distance within the gapped (gapless) phases. Surprisingly, for several models with LR interactions it was reported that algebraic decay of correlations can coexist with an initial exponential decay within gapped phases [12,22,53]. An analytic understanding of this effect has so far proven elusive.
In the following we use the decay of correlation functions to characterize the various phases. In particular, for the LRI chain we provide extensive numerical results using DMRG techniques, while for the LRK chains we exploit the integrability of the models to derive an analytic expression of the behavior of correlation functions in all parameter regimes.

Edge states and edge gaps
Localized edge states within topological phases have attracted much interest over the last decade, largely because of possible applications in schemes for topological quantum computing [64]. In particular, [11] has shown that localized states arise at the edge of a 1D superconductor with p-wave SR pairing interactions (described e.g. by the limit a  ¥ of equation (2)). The existence of these localized states is related with the spontaneous breaking of the discrete  2 symmetry associated with the parity of the fermion number: when this  2 symmetry is broken, two degenerate ground states with different parity appear. Here, they will be labeled by + 0 and -0 in the even and odd parity sectors, respectively.
As Hamiltonians equations (1) and (2) are equivalent in the limit a  ¥, the same breaking of  2 symmetry (now related to spin flips along x direction) described above occurs also for the LRI chain, resulting in the presence of two degenerate edge states in this model. From the discussion above it turns our that the analysis of the edge modes can be utilized to characterize the quantum phases of the system.
In this work, we identify the localized states by directly computing their wavefunctions and masses, which can be achieved either numerically for equation (1) using DMRG simulations or exactly for equations (2). In order to accomplish this task for the LRI model, it is useful to exploit the Jordan-Wigner transformation equation (4) to define new fermionic operators (similar to equation (4)) from spin operators s  j . We then compute the wave-functions j ( ) j 1,2 of the massless edge modes as [65] Here, the states + 0 and -0 are the ground states in the even and odd parity sector also for the LRI model, respectively. The mass of the two modes j ( ) j 1,2 , also known as edge gap (in order to distinguish it from the usual mass gap, i.e. the energy difference between the ground state and the first excited bulk state, see [50,66]), is defined as the difference --+ E E 0 0 , with  E 0 being the energy of the state  0 . In the following we will also be interested in characterizing the localization of massive edge states that are found in the phase of LRI where the  2 symmetry is preserved, and thus where a unique ground state + 0 occurs in the even-parity sector. This localization arises for the first two excited states -1 and -2 for  a 1, which are instead delocalized in the bulk for a  1 (see section 4.2). Their wavefunctions read Here the mass of the mode ( ) w j s is defined as the difference D = - The difference between equations (12)- (14) is that in the first two expressions the second term on the right-hand side of the equations is nonzero because of the zero-energy condition [65].
For the LRK models, the wavefunctions for both the massless and massive modes can be extracted following [54]. The latter describes a technique for the exact diagonalization of a generic fermionic quadratic Hamiltonian of the form Here, L n are single-particle energies (ordered as with energy L n , while g nj and h nj are the wavefunctions of h n and h † n .

Phase diagrams of LRI and LRK models
In this section, we present the phase diagrams for the LR models equations (1) and (2), obtained from an analysis of the observables described above. The results for the LRI model were obtained numerically via DMRG techniques for chains of a length L up to L = 200. For all calculations we utilized up to 128 local basis states and 10 finite-size sweeps [51,52]. The discarded error on the sum of the eigenvalues of the reduced density matrix was always less than 10 −8 . For the LRK models all results were obtained (semi-)analytically.  (1). Hamiltonian equation (1) is invariant under the transformation q q cos cos and thus the phase diagram is symmetric around q p = 2.
We find that for  a 1, c eff is zero everywhere except along two critical lines. By comparing with results for the energy gap (not shown), we find that the critical lines separate two gapped regions (denoted as PM1 and AM in figure 1(a)) that for a  ¥ correspond to the known paramagnetic and antiferromagnetic phases of the SR model. Similar to [26], we find that the behavior of the full correlation functions s s i x j x and s s i z j z is consistent with the persistence of paramagnetic and antiferromagnetic orders for all α. However, different from the SR model, we find that the connected correlation functions decay with distance with a hybrid behavior that is exponential at short distances and algebraic at long ones. An example is shown in figure 2(a) for C R xx 1, in the PM1 phase. Surprisingly, we find numerically that the exponent g x of the long-distance decay for C R xx 1, displays three difference behaviors: (i) for a > 2 it fulfills g a = x , consistent with the results of [12,26]. However, (ii) for a < < 1 2 we obtain a hybrid exponential and algebraic decay with a different γ that depends linearly on α with a slope consistent, with~( ) 0.55 5 and (iii) for  a 1 we observe numerically a curve compatible with a pure algebraic decay, with an α-dependence of g x that is linear with slope~( ) 0.25 2 . The fitted exponent g x is shown in figure 2(c).
For  a 1 in the paramagnetic regions of the phase diagram denoted as PM2 we find that the effective central charge grows continuously with decreasing α from zero to a finite value that appears to be θ-dependent and has a maximum of order 1 for a = 0 and q p » 2. An example for q p = 0.2 is plotted in figure 1(d) (blue squares). As mentioned above, in this PM2 region, the correlation function C R xx 1, is found numerically to decay as an almost pure power-law.
The energy spectrum changes in this region PM2 compared to the case PM1: the energy gaps D ( ) E 1 and D ( ) E 2 between the ground state and the first excited states in the odd parity sector increase with decreasing α, as shown in figure 8(b), and the wave-functions of the two lowest-energy excited states -1 and -2 , defined in section 2.2.3, become localized at the edges of the chain (figure 7(b)).
In the antiferromagnetic phase, denoted as AM, the effective central charge is instead zero for all α [26]. For a > 1 the connected correlation functions C R zz 1, display a clear algebraic decay at long-distances (see the example in figure 6(a) below), while our numerical results do not allow for establishing whether an initial exponential decay is also present, as expected. The ground-state is found to be doubly degenerate for all α. This degeneracy is due to the spontaneous breaking of the  2 spin-flip symmetry [58,67], and is related to the two modes that are localized at the edges of the chain as in the short-range Ising model 5 . While a LR power-law tail may be present [50], the localization of these modes is here found to be consistent with exponential up to numerical precision (see figure 7(b) below). We come back to this pointbelow.

LR Kitaev models.
The phase diagram of the LRK model equation (2) for  = 0 is reported in figure 1(b), where we plot the effective central charge c eff defined in section 2.2.1 as a function of the angle θ and the power α of the decay of the pairing term. In this case the invariance of equation (2) under q q cos cos is lost for any finite α and the phase diagram is not symmetric around q p = 2. Figure 1(b) shows that for  a 1 and q p < < 0 , two phases exist that are denoted as PM and AM1 separated by two critical lines. In the limit a  ¥ these phases correspond to the paramagnetic and antiferromagnetic phases of the LRK model and are gapped. Consistently, figure 1(b) shows that the effective   along the critical lines, as expected from general results for SR systems [57].
The PM and AM1 phases are distinguished by different asymptotic values of the correlation functions S( ) i j , defined in section 2.2.2. In the region denoted as AM1, S( ) i j , has a finite value for - ¥ |i j , while S( ) i j , decays for - ¥ |i j within the PM phase. Similar to the situation in the LRI model (see above), the decay to zero of S( ) i j , in the PM phase shows a hybrid exponential and power-law behavior with distance. This is shown for a particular value of θ in figure 2(b), where we find numerically that the exponent g S for the power-law tail of S( ) i j , equals g a = z when a > 1. For a < 1, however, the exponential part becomes numerically unobservable, and S( ) i j , decays essentially algebraically within the PM phase with an exponent that grows to 2 for a  0 (see figure 2(b)).
Remarkably, we show below in section 3 that the hybrid exponential and algebraic behavior described above can be obtained analytically in all phases for several correlation functions, such as the one-body and the densitydensity correlation functions. In particular, the leading contribution to the one-body correlation function † where the pre-factors  a q , and  a q , are derived below in section 3. The algebraic part of the decay of ( ) 2 for a > 1 and a < 1, respectively. In section 4.1, we show that a similar hybrid exponential and power-law decay is found for the localization of the edge modes within the antiferromagnetic phase AM1: for  a 1, the edge modes are massless, as expected from the SR Kitaev model. However, for  a 1 the edge modes acquire a finite mass, i.e., become gapped. This is shown in figure 1(c), where we plot the edge gap L 0 defined in section 2.  , respectively. A Bogoliubov transformation brings equation (20) in diagonal form as Here, the new fermionic operators h h -( ) † , k k are given in terms of the operators The ground state of equation (21) is the vacuum of the h k fermions and has an energy density  figure 1(a)). In [53] it was demonstrated for a related model with LR pairing only that this behavior corresponds to an exotic change for the decay of density-density correlation functions: For  a 1 their oscillations mimic those of a Luttinger liquid. Here, we find a similar behavior (not shown). The increasing of c eff , below a = 1, is also found in the very recent work [70] where the scaling of the von Neumann entropy in the thermodynamic limit is analytically analized. This anomalous behavior of c eff points towards a breaking of conformal symmetry along the critical line, which we analyze further below.
The breaking of conformal symmetry can be inferred also by analyzing the scaling of the energy density a ( ) e L , 0 with the system size L [53]: for a conformally invariant theory the following relation must hold [ does not satisfy the scaling law (24), which implies a breaking of conformal symmetry. This behavior also implies that the quantum phase transition between the PM and the AM1 phase for q p < 2 and  a 2 is in a different universality class from that of the SR Kitaev (Ising) model. We notice that, even if conformal symmetry is broken, the von Neumann entropy predicts a value for c eff which tends to 1 as α goes to zero, compatible with the observed decay of density-density correlation functions.
We further confirm the breaking of conformal symmetry for the fermionic models by looking at the behavior of their low-energy spectra. The dispersion relation of a conformal field theory is linear in the momentum k, implying a dynamical exponent z = 1 [59]. Consistently, by expanding the dispersion relation l a ( ) k for the long range Kitaev models on the critical line with q p > 2, for a > 2 we find lã ( ) k k, for  k 0. However, for a < < 1 2 we obtain the scaling lã a-( ) k k 1 . This latter scaling implies a dynamical exponent a =z 1 that varies continuously with α and is different from that of a conformal field theory. This would imply that the quantum phase transition between the PM and the AM1 phase for q p > 2 and a < 2 is in a different universality class from that of the SR Kitaev model. The appearance of a new universality class due to LR interactions is also found in [71,72]. Incidentally, we notice that linearity of the spectrum around the minimum is only a necessary condition for the persistence of conformal invariance: indeed along the critical line at q p < 2 conformal symmetry breaking arises even if the low-energy spectrum is linear for every α.

LR Ising model.
We locate the critical line (for q p < 2, the other being symmetric) of the LRI model numerically by using two complementary ways that agree up to finite-size effects. Firstly, we determine the points in the phase diagram of figure 1(a) where the energy gap between the ground state and the first excited state reaches its minimum. Then, we compute the effective central charge c eff for different system sizes L from equation (6) and determine the precise values of α and θ for which its value does not depend on L [73,74]. Examples of this latter technique, which is found to be particular precise, are presented in figure 3 for different α. We notice however that this method does not allow us to extract a precise value for c eff when  a 0.2, since within our numerical results lines with different L do not cross at a single point in this region.
On the critical line, we find that for  a 1 (black solid lines in figure 1) c eff is equal to 1/2 as expected for the central charge of the critical SR Ising model. However, for  a 1 (red dashed lines in figure 1), c increases continuously up to a value of order 1 as shown in figure 3(b). We argue that on this line the conformal invariance of the model is broken as the found values of c do not coincide with the discrete set allowed for the known conformal field theories [58,67]. Based on the mismatch between the predictions for c from the von Neumann entropy and the ground state energy density found in the previous subsection for the LRK model, we cannot exclude here a conformal symmetry breaking also in a certain range for α above a = 1. However our results do not allow us to provide a final answer, since our DMRG calculations for the LRI model cannot reach sufficiently large sizes to perform a satisfying finite-size scaling for the energy density.

Correlation functions for the LRK
In this section, we present an analytical calculation of the one-body correlation functions for the LRK models. The latter display a hybrid exponential and algebraic decay with distance that is explained by exploiting the integrability of the models. Higher-order correlation functions, such as the density-density correlations, are readily obtained from these correlations via Wick's theorem (see section 2.3.2 and below). respectively. In the following, we focus first on the one-body correlation function † a a R 0 and come back to the anomalous correlation † † a a R 0 in section 3.3 below. In order to evaluate equation (25), we use the Cauchy theorem applied to the contour in the complex plane drawn in figure 4 and = + z k y i , and where we have chosen  = 0. In equation (27) we have neglected the contributions from ^and  ¢^as they vanish when  ¥ M . As we explain below, the integrations over the lines   and  p 0,2 (and thus momenta p  k and  k 0 and p 2 , respectively) are responsible for the exponential and algebraic behavior observed in these models, respectively.
in the limit  y 0 [75], and reads The decay constant ξ is related to the zeroes of the denominator of  a ( ) z and is obtained by solving the equation Two cases must be distinguished: Li e 0 admits a solution, since the function -- Li e 0 admits a solution, for the same reason as above. In the following we focus, without loss of generality, on the first case, where ξ is solution of Li e . Notably for a = 0 solutions exist only for q p equation (31) does not admit any solutions, which implies the absence of exponential decay. This is in contrast with, e.g., the expected behavior of correlation functions within gapped phases for SR models. Figure 4(b) shows the decay constant x 1 as a function of α for two different values of θ. In particular, for q p = 0.2 and a  ¥, x 1 tends to the SR value (x q  | | log tan ), as expected. However, for a  0 we find that x 1 tends to zero, essentially linearly with α. As explained below, this can result in the non-observability of the exponential dependence of correlation functions for  a 1. Notice however that even if ξ is finite when  a 1 for q p = 0.4 , the exponential part of the correlation functions is still unobservable.

Algebraic decay
The sum of the integrals on the lines  0 (where  While, e.g., the phase diagram of figure 1(b) demonstrates the persistence of individual paramagnetic and antiferromagnetic phases with varying α, the analytic expressions equations (34) and(35) for the one-body correlation function clearly show that different regions are in fact present within each phase.

Hybrid decay and other correlations
The two contributions from equations (29) and (34) sum up to give the hybrid exponential-algebraic behavior of equation (19) valid for all α. Figure 5(a) shows that the analytical results are in perfect agreement with a numerical solution of equation (25). The anomalous correlation † † a a R 0 can be computed along the same lines as before and is given by ⎧ and, from equations (34) and (36), the leading part of its LR power law tail is found to be ⎧ We notice that a hybrid exponential and algebraic behavior similar to that described above has been already observed numerically in certain spin and fermionic models, e.g., in [12,26,50,53]. This behavior is characteristic of the non-local interactions and from our analysis appears to be largely unrelated to the presence or absence of a gap in the spectrum. In fact, we have shown here for Hamiltonians equations (2) that this hybrid decay does not require exotic properties of the spectrum (see [21] and [22]), rather is due to the different contributions of momenta p = k (as for the SR limit) and p = k 0, 2 respectively. In particular, the latter momenta are responsible for the LR algebraic decay.
Finally, as mentioned before, we notice that (i) the contribution to the imaginary part in equation (32) is due to a ( ) Li e y and disappears in the limit a  ¥. This implies a simple exponential decay of correlations as i j , in figure 2, since the initial exponential decay is too small to be observed, as expected from the discussion above. Moreover, no universal behavior for the decay exponents is identifiable for C i j xx , and C i j zz , in this region. The exponents for C i j xx , and S( ) i j , differ here in general, probably because the contribution of the string operators in H LRI becomes more relevant.
In the AM phase, instead, we find that the decay of C i j zz , (figure 6(a)) for  a 1 displays an algebraic tail with an exponent compatible with g a = z , mimicking C i j xx , in the PM phase (a very precise estimate is forbidden by the decay oscillations between even and odd sites). These results are consistent with those of [12], obtained for the case a = 3. Notably the decay exponent is here always different from the value a 2 analytically computed for ( ) g i j , 2 in equation (39). This discrepancy is again probably due to the role of the string operators, here quantitatively relevant even for  a 1. For  a 1 again no universal behavior for g z is found.

Edge modes properties 4.1. Massless edge modes
It is known that the SR Kitaev chain hosts modes localized at the edges [11] in the ordered phase for p q p < < 4 3 4 . These modes are fermionic and massless (in the limit  ¥ L ), thus they have Majorana nature [76] and are a consequence of the topological non-triviality of the ordered phase (see e.g. [77] and references therein).
For the LRK models of equation (2), Majorana massless modes are found for  a 1 in the AM1 region of the phase diagram of figure 1(b). Plotting the square of the wavefunction y j 0 corresponding to the zero edge gap L 0 , defined in section 2.2.3, we find for  ¹ 1 in equation (2) a hybrid exponential and algebraic decay with the distance from one edge of the chain. The exponent of the algebraic decay of y j 0 2 is found to be equal to a 2 . This behavior is similar to that observed in [53] in the presence of LR pairing only. Interestingly, here the algebraic tail of the Majorana modes can be tuned to completely disappear by changing the parameter ò that fixes the unbalance between the hopping and the pairing terms in the Hamiltonian of equation (2). Figure 7(a) shows an example for q p = 0.4 and a = 3. When  ¹ 1 the hybrid behavior is fully visible and the state y j 0 decays in the bulk with an algebraic tail. However, by approaching the value  = 1 the algebraic tail of y j 0 decreases and eventually disappears. As a result, for  = 1 the wave function y j 0 becomes exponentially localized at one edge. We find that this exponential localization for  = 1 is present also in the parameter region a < 1. Moreover, the edge gap L 0 scales to zero exponentially with the system size L for all α as figure 8(c) shows.
For the LR Ising Hamiltonian of equation (1), edge modes appear in the AM phase for every a > 0. They have zero mass, since the edge gap - Examples of these edge modes for different α are given in figure 7(b), where we plot the square of the wavefunction j ( ) j 1 defined in equation (12). We find numerically that j ( )

Massive edge modes
If we extend the analysis of the LR Kitaev Hamiltonian of equation (2) to different ò and sufficiently small α, a totally new situation arises for the edge gap L 0 and the edge modes: in the region denoted as AM2 in the phase diagram of figure 1(b), we find that L 0 (which is zero for  a 1) becomes nonzero for  a 1 and  ¹ 1, also in the thermodynamical limit. This case is shown in figure 1(e), where we plot the edge gap L 0 as function of α and θ for  = 10. Between  = 1 and  = 10 we checked a continuous increase of the extension of the region AM2 and no transitions in between. Similarly in figure 8(a) we plot L 0 together with several other single-particle energies as a function of α and for a fixed θ. These energies have been computed as described in section 2.2.3. In figure 8(a) the mass of the edge mode L 0 is easily recognizable for all ò, since it is separated from all bulk modes by a finite gap.
Consistently with the discussion above, for  a 1, L 0 is zero as expected from the SR model, so that two degenerate ground states exist as the  2 symmetry of the model is spontaneously broken. However, surprisingly, for  a 1 and  ¹ 1, we find that L 0 becomes finite and thus the ground state is unique. This indicates that the  2 symmetry of the model (which is broken for  a 1 in the AM1 phase) is restored for  a 1. As a consequence, the region AM2, where the  2 symmetry is restored, must be separated from AM1 by a quantum phase transition, even if no closure for the mass gap arises in the bulk. The wavefunction of the lowest massive L 0 state is now given by the matrix element g 0i defined in section 2.2.3. By plotting the probability density g i 0 2 , we now obtain a localization on the edges that is symmetric with respect to the middle of the chain. This probability density decays algebraically when approaching half of the chain, as is clearly seen in figure 7(d).
A similar wavefunction localization at the edges of the system is found also for the LRI model in the PM2 region in the phase diagram of figure 1(a) for  a 1, where the  2 symmetry is preserved. However, while for the LRK massive edge modes originate from Majorana edge modes present at  a 1, for the LRI the edge localization arises for excited states of the bulk spectrum in the region PM1. For  a 1, these states are degenerate, as shown in figure 8(b), and separated from the third excited state by a gap that is finite in the thermodynamic limit. Because of this degeneracy, we consider the probability density = + with ( ) w j 1,2 wavefunctions of |( ) 1, 2 defined in section 2.2.3. A typical situation is depicted in figure 7(d), where we plot p j as function of the lattice site j for different values of α. For  a 1, p j is oscillating and delocalized in the bulk, while it is localized exponentially at the edges for  a 1 and is symmetric with respect to half of the chain.
We leave as an open question whether the edge localization here signals the appearance of a new phase (and without mass-gap closure) with preserved  2 symmetry, similar to the LRK models above.

Observability in current experiments
Recent experiments with cold ions have made possible the realization of LR Ising-type Hamiltonians as equation (1) with   a 0 3.5 [15,16,18,20]. In these experiments, both static and dynamical spin-spin correlations, as well as the spectrum of quasi-particle excitations [25], can be measured with extreme precision, which in principle could allow for an analysis of some of the observables discussed above. For example, the observation of long-distance algebraic of correlations, as well as spectroscopic signatures of the formation of localized excited edge modes for  a 1 could allow for the precise determination of the properties of these LR models. Figure 8. (a) First 45 low energy states of the LRK models for L = 500 and q p = 0.7 as function of α. It is possible to distinguish the L 0 mode (solid points for three different ò) separated from the L > n 0 bulk modes (light green triangles for  = 10). When  a 1 the vanishing of L 0 for  ¹ 1 signals the spontaneous breaking of the  2 symmetry of LRK and defines the AM1 phase (see also the phase diagram in figure 1(b)). When L > 0 0 the  2 symmetry is restored and the phase AM2 appears. (b) Main panel: D ( ) E 1 (empty circles) and D ( ) E 2 (solid triangles) for the LRI as functions of the inverse of the length of the chain L for q p = 0.159 and for three α reported in the plot. These two states are degenerate in the  ¥ L limit, as finite-size scaling shows. Inset: D ( ) E 1 as function of α for q p = 0.159 . When α decreases, D ( ) E 1 increases, showing a change for the energy spectrum in the region PM2 of the phase diagram of figure 1(a). (c) Edge gap L 0 for (squares)  = 0.5 and two different α as function of the system size L. When  = 1, the edge gap scales exponentially to zero for all α. For comparison, the case (crosses)  = 2 and a = 3 is also reported. There the edge mass displays an hybrid exponential and power-law decay.
One key aspect of experiments, however, is that experimentally attainable lengths for ion chains are currently limited to at most few tens of ions. It is thus natural to ask whether the characteristic long-distance decay of correlation functions described above can be observed in systems of such length. To explore this issue, figure 6 (right panel) shows the correlation C R xx 1, for Hamiltonian equation (1) in the PM1 and PM2 phases for a chain of L = 30 sites with open boundary conditions and for different α. For a > 2, the initial exponential decay dominates the correlations for  R 10, while a comparatively small algebraic tail is found for  R 10. For  a 1, however, the exponential part has essentially disappeared and the decay is purely algebraic at all distances, as expected from the discussion of sections 3.4. This fundamental change of behavior around a  1 may be observable. We note, however, that the exponent g x of the algebraic decay is here different from that presented in figure 2(a), due to strong finite size effects in these systems. We find similar results for the correlation C R zz 1, . On the other hand, the emergence of massive edge modes in the LRI chain could be a convenient diagnostic of the change of nature of the paramagnetic phase for  a 1.

Summary and outlook
In this work we have analyzed the phase diagram of the long range anti-ferromagnetic Ising chain and of a class of fermionic Hamiltonian of the Kitaev type, with LR pairing and hopping. We have clarified in what regions of the phase diagram violation of the area law occurs, and have provided numerical evidence and exact analytical results for the observed hybrid decay of correlation functions, which are found to decay exponentially at short range and algebraically at long range, for all α. We have further demonstrated the breaking of conformal symmetry along the critical lines in both models at low enough α. Most interestingly, for the fermionic models we have demonstrated for the first time that the topological edge modes can become massive for sufficiently small values of  a 1. This implies the existence of a transition to a novel phase without closure of the mass gap, to the known phase with massless Majorana modes for  a 1. We conjecture that the possibility of a phase transition with nonzero mass gap is due to the peculiar behavior of LR correlations, showing power-law tails also when the gap does not vanishes. Similarly, we have found that excited bulk states in the paramagnetic phase of the Ising model can become localized at the edges of the chain for  a 1.
This works may open several exciting research directions. The first question concerns the nature and topological properties of the proposed new phase of the Kitaev model with  a 1, and of its localized edge modes. We conjecture that these massive edge modes are due to the hybridization of the Majorana modes at small α, due to the bulk overlap between their wave functions, whose decay is slower and slower for decreasing α. This aspect will be the subject of future studies. Another important open question is whether the appearance of massive edge modes may be connected also to the violation of the area-law for the entanglement entropy in these models.
In general, these results represent counter-examples for the topological properties of existing topological models with LR interactions, as recently analyzed in [50]. The question of whether a possible universal behavior exists for topological models with LR interactions is thus still wide open.