# Long-range Ising and Kitaev models: phases, correlations and edge modes

Davide Vodola, Luca Lepori, Elisa Ercolessi, Guido Pupillo

► To cite this version:

Davide Vodola, Luca Lepori, Elisa Ercolessi, Guido Pupillo. Long-range Ising and Kitaev models: phases, correlations and edge modes. New Journal of Physics, 2015, 18 (1), pp.015001. 10.1088/1367-2630/18/1/015001. [hal-04834110](https://hal.science/hal-04834110v1)

HAL Id: hal-04834110

https://hal.science/hal-04834110v1

Submitted on 12 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Long-range Ising and Kitaev models: phases,
correlations and edge modes

To cite this article: Davide Vodola et al 2016 *New J. Phys.* **18** 015001

View the [article online](#) for updates and enhancements.

Related content

- Singular dynamics and emergence of
  nonlocality in long-range quantum models
  L Lepori, A Trombettoni and D Vodola

- Long-range topological insulators and
weakened bulk-boundary correspondence
L Lepori and L Dell'Anna

- Disordered Kitaev chains with long-range
pairing
Xiaoming Cai

Recent citations

- Phase diagram of the quantum Ising
model with long-range interactions on an
infinite-cylinder triangular lattice
S. N. Saadatmand et al

- Local and nonlocal order parameters in
the Kitaev chain
Gennady Y. Chitov

- Edge insulating topological phases in a
  two-dimensional superconductor with long-
  range pairing
  L. Lepori et al
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

27 August 2015

REVISED

18 November 2015

ACCEPTED FOR PUBLICATION

23 November 2015

PUBLISHED

18 December 2015

Content from this work
may be used under the
terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Long-range Ising and Kitaev models: phases, correlations and edge modes

Davide Vodola¹,⁴, Luca Lepori¹,²,⁴, Elisa Ercolessi³ and Guido Pupillo¹

¹ icFRC, IPCMS (UMR 7504) and ISIS (UMR 7006), Université de Strasbourg and CNRS, F-67000 Strasbourg, France

² Dipartimento di Fisica e Astronomia, Università di Padova, Via Marzolo 8, I-35131 Padova, Italy

³ Dipartimento di Fisica e Astronomia, Università di Bologna and INFN, Via Irnerio 46, I-40127 Bologna, Italy

⁴ These two authors contributed equally to the present work.

E-mail: pupillo@unistra.fr

**Keywords:** long-range interactions, spin chains, Kitaev chains, area law violation, edge modes, Majorana fermions

## Abstract

We analyze the quantum phases, correlation functions and edge modes for 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 long-range interactions that decay with distance $r$ as $1/r^\alpha$, as well as a related class of fermionic Hamiltonians that generalize 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 \lesssim 1$, while connected correlation functions can decay with a hybrid exponential and power-law behavior, 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.

## 1. 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–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–17], which are naturally long ranged. The resulting Ising-type interactions are antiferromagnetic and decay with distance $r$ as a power-law $1/r^\alpha$, with an adjustable exponent $\alpha$ usually in the range $0 \lesssim \alpha \lesssim 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–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 Ising-chain 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 $\alpha \gtrsim 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–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 $\alpha \lesssim 1$ as well as sub-logarithmically for $\alpha > 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 $\alpha$.

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–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 $1/r^\alpha$, and their relative strength can be varied.

For these models, we provide the phase diagram for all exponents $\alpha$, 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-matrix-renormalization group (DMRG) method [51, 52], while the long-range Kitaev-type (LRK) models remain exactly solvable for all $\alpha$, 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 $\alpha \lesssim 1$. For $\alpha \gtrsim 1$, no violation is found.

(ii) For *any finite* $\alpha$, 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 $\alpha$. However, when $\alpha \lesssim 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 $\alpha$, 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 $\alpha$ 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 $\alpha$. The antiferromagnetic phase for the LRK models is defined in analogy with that of the short-range Ising chain. The massive modes, instead, are entirely new and are found to appear in a broad area of the antiferromagnetic phase for $\alpha \lesssim 1$ and when we choose the unbalance $\epsilon$ between the strengths of the hopping and pairing terms to be different from 1. This results suggests, for $\alpha \lesssim 1$, a restoration of the $\mathbb{Z}_2$ symmetry associated with the (absence of) ground state degeneracy, as well as a possible transition to a novel $\mathbb{Z}_2$ symmetric phase and without mass gap closure. Interestingly, if we choose $\epsilon = 1$, the massless modes survive for all $\alpha > 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 $\alpha$, 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 a localization of excited, gapped, energy eigenstates for $\alpha \lesssim 1$, which for $\alpha \gtrsim 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 $\alpha$. 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 $\alpha \gtrsim 1$. In section 4.2, instead, we demonstrate that the edge modes can become massive for $\alpha \lesssim 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 $\alpha \lesssim 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.

## 2. 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.

### 2.1. Model Hamiltonians

#### 2.1.1. LR Ising model

In this work, we are interested in Ising-type Hamiltonians with LR interacting terms. The LRI Hamiltonian [26] reads

$$H_{\text{LRI}} = \sin \theta \sum_{i=1; j>i}^{L} \frac{\sigma_i^x \sigma_j^x}{|i-j|^\alpha} + \cos \theta \sum_{i=1}^{L} \sigma_i^z, \quad (1)$$

where $\sigma_j^\nu$ ($\nu = 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 $\sin \theta > 0$ (or equivalently $0 < \theta < \pi$). 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 $\alpha \to \infty$) the Hamiltonian equation (1) is exactly solvable and a quantum phase transition between these two phases is known to occur at $\theta_c = \pi/4$. Reference [26] has shown numerically that a quantum phase transition separating the AM and the PM phases survives for all finite $\alpha \gtrsim 0.5$. Below, we are interested in exploring the phase diagram of equation (1) for all $\alpha$ and $\theta$.
### 2.1.2. LRK chains

Related to the LRI chain, in the following we introduce and analyze a class of fermionic Hamiltonians of the form

$$H_{\text{LRK}} = \sin \theta \sum_{i,j} \frac{a_i^\dagger a_j + (1 + \epsilon) a_i a_j + \text{h. c.}}{|i-j|^\alpha} + 2 \cos \theta \sum_i n_i, \quad (2)$$

where $a_j^\dagger$ describes the creation operator for a fermionic particle at site j, and $n_j = a_j^\dagger a_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 $\alpha$. Here, all energies are expressed in dimensionless units and the parameter $\epsilon$ governs the unbalance between the hopping and pairing terms.

In the SR limit $\alpha \to \infty$, equation (2) maps into the Ising chain equation (1) via the Jordan-Wigner transformation [54]

$$\sigma_j^\dagger = a_j^\dagger \exp \left( i \pi \sum_{\ell<j} n_\ell \right); \quad (3)$$

$$\sigma_j^z = 2n_j - 1, \quad (4)$$

with $\sigma_j^\dagger = (\sigma_j^x + i\sigma_j^y)/2$. However, at finite $\alpha$ this identification does not hold anymore due to the contributions of the string operators $\exp(i\pi \sum_{\ell<j} n_\ell)$ in equation (4). In particular, unlike the LRI, equation (2) remains exactly solvable, allowing for analytic solutions at any finite $\alpha$.

## 2.2. Observables

### 2.2.1. 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 $S_\ell$ that we employ in this work. For a system of L sites that is partitioned into two disjoint subsystems A and B containing $\ell$ and $L-\ell$ sites, respectively, $S_\ell$ is defined as

$$S_{\ell} = -\text{Tr} \rho_{\ell} \log_2 \rho_{\ell}, \quad (5)$$

where $\rho_\ell$ is the reduced density matrix of the subsystem A.

Two general behaviors of $S_\ell$ are known for the ground states of one-dimensional SR interacting systems. Within gapped phases, $S_\ell$ saturates to a constant value independently of $\ell$ and thus obeys to the so-called area law [56]. On the contrary, $S_\ell$ diverges with $\ell$ for critical gapless phases and, for conformally invariant systems, satisfies the universal law [57]:

$$S_{\ell} = a + \frac{c}{6} \log_2 d(\ell, L), \quad (6)$$

with $d(\ell, L) = (2L/\pi)\sin(\pi\ell/L)$ 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]. $S_\ell$, 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–63].

In the case of LR models it has been shown [26, 53, 56] that the divergence of $S_\ell$ 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_{\text{eff}}$ may be defined also within the gapped phases and used to characterize the main features of the phase diagram for LR interactions [26].

### 2.2.2. 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

$$C_{i,j}^{\nu\nu} = \langle \sigma_i^{\nu} \sigma_j^{\nu} \rangle - \langle \sigma_i^{\nu} \rangle \langle \sigma_j^{\nu} \rangle, \quad (7)$$

with $\nu = x, z$.

For the LRK models we are interested both in the density–density correlation function

$$g_2(i, j) = \langle n_i n_j \rangle - \langle n_i \rangle \langle n_j \rangle \quad (8)$$
and the function

$$ \Sigma(i, j) = \left\langle \left( a_i^\dagger + a_i \right) \exp \left( \pi i \sum_{\ell=1}^{j-1} a_\ell^\dagger a_\ell \right) \left( a_j^\dagger + a_j \right) \right\rangle \quad (9) $$

that correspond to the functions $C_{i,j}^{zz}$ and $\langle \sigma_i^x \sigma_j^x \rangle$ 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 $\langle a_i^\dagger a_j \rangle$ and $\langle a_i^\dagger a_j^\dagger \rangle$ via Wick's theorem. In particular, one finds

$$ g_2(i, j) = \left| \langle a_i^\dagger a_j^\dagger \rangle \right|^2 - \left| \langle a_i^\dagger a_j \rangle \right|^2 \quad (10) $$

and

$$ \Sigma(i, j) = \det \begin{pmatrix} G_{i,i+1} & \cdots & G_{i,j} \\ \vdots & \ddots & \vdots \\ G_{j-1,i+1} & \cdots & G_{j-1,j} \end{pmatrix} \quad (11) $$

with $G_{m,n} = \delta_{mn} + 2 \langle a_m^\dagger a_n^\dagger \rangle + 2 \langle a_m^\dagger a_n \rangle$.

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.

### 2.2.3. 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 $\alpha \to \infty$ of equation (2)). The existence of these localized states is related with the spontaneous breaking of the discrete $\mathbb{Z}_2$ symmetry associated with the parity of the fermion number: when this $\mathbb{Z}_2$ symmetry is broken, two degenerate ground states with different parity appear. Here, they will be labeled by $|0^+\rangle$ and $|0^-\rangle$ in the even and odd parity sectors, respectively.

As Hamiltonians equations (1) and (2) are equivalent in the limit $\alpha \to \infty$, the same breaking of $\mathbb{Z}_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 out 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 $c_j^\dagger = \sigma_j^\dagger \exp(i\pi \sum_{\ell<j} \sigma_\ell^\dagger \sigma_\ell^-)$ (similar to equation (4)) from spin operators $\sigma_j^\pm$. We then compute the wave-functions $\varphi_j^{(1,2)}$ of the massless edge modes as [65]

$$ \varphi_j^{(1)} = \langle 0^- | c_j^\dagger | 0^+ \rangle + \langle 0^+ | c_j^\dagger | 0^- \rangle, \quad (12) $$

$$ \varphi_j^{(2)} = i \left( \langle 0^- | c_j^\dagger | 0^+ \rangle - \langle 0^+ | c_j^\dagger | 0^- \rangle \right). \quad (13) $$

Here, the states $|0^+\rangle$ and $|0^-\rangle$ are the ground states in the even and odd parity sector also for the LRI model, respectively. The mass of the two modes $\varphi_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_{|0^-\rangle} - E_{|0^+\rangle}$, with $E_{|0^+\rangle}$ being the energy of the state $|0^+\rangle$.

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 $\mathbb{Z}_2$ symmetry is preserved, and thus where a unique ground state $|0^+\rangle$ occurs in the even-parity sector. This localization arises for the first two excited states $|1^-\rangle$ and $|2^-\rangle$ for $\alpha \lesssim 1$, which are instead delocalized in the bulk for $\alpha \gg 1$ (see section 4.2). Their wavefunctions read
$$w_j^{(1)} = \langle 1^- | c_j^\dagger | 0^+ \rangle \qquad w_j^{(2)} = \langle 2^- | c_j^\dagger | 0^+ \rangle. \qquad (14)$$

Here the mass of the mode $w_j^{(s)}$ is defined as the difference $\Delta E(s) = E_{|s^-} - E_{|0^+}$ with $s = 1, 2$. 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

$$H = \sum_{i,j=1}^{L} \left[ A_{ij} a_i^\dagger a_j + \frac{1}{2} (a_i^\dagger B_{ij} a_j^\dagger + \text{h. c.}) \right], \qquad (15)$$

with $A = [A_{ij}]$ and $B = [B_{ij}]$ real matrices. Equation (15) can be cast in diagonal form as

$$H = \sum_{n=0}^{L-1} \Lambda_n \eta_n^\dagger \eta_n. \qquad (16)$$

by a singular value decomposition of the matrix $A + B$:

$$\Lambda_n = \sum_{ij} \psi_{nj} (A_{ji} + B_{ji}) \phi_{in}. \qquad (17)$$

Here, $\Lambda_n$ are single-particle energies (ordered as $\Lambda_0 \le \Lambda_1 \le \dots \le \Lambda_{L-1}$) and $\eta_n$ are fermionic operators defined by the following Bogoliubov transformation

$$\eta_n = \sum_j g_{nj} a_j + h_{nj} a_j^\dagger. \qquad (18)$$

The matrix elements $\phi_{nj} = g_{nj} + h_{nj}$ and $\psi_{nj} = g_{nj} - h_{nj}$ can be directly identified as the wavefunctions of the two Majorana modes $(\eta_n + \eta_n^\dagger)/2$ and $i(\eta_n - \eta_n^\dagger)/2$ with energy $\Lambda_n$, while $g_{nj}$ and $h_{nj}$ are the wavefunctions of $\eta_n$ and $\eta_n^\dagger$.

## 2.3. 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.

### 2.3.1. LR Ising model.

Our results for the phase diagram of the LRI model are summarized in figure 1(a), where we plot the effective central charge $c_{eff}$ defined in section 2.2.1 as a function of the angle $\theta$ and the power $\alpha$ of the antiferromagnetic term in equation (1). Hamiltonian equation (1) is invariant under the transformation $\cos \theta \to -\cos \theta$ and thus the phase diagram is symmetric around $\theta = \pi/2$.

We find that for $\alpha \gtrsim 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 $\alpha \to \infty$ 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 $\langle \sigma_i^x \sigma_j^x \rangle$ and $\langle \sigma_i^z \sigma_j^z \rangle$ is consistent with the persistence of paramagnetic and antiferromagnetic orders for all $\alpha$. 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_{1,R}^{xx}$ in the PM1 phase. Surprisingly, we find numerically that the exponent $\gamma_x$ of the long-distance decay for $C_{1,R}^{xx}$ displays three difference behaviors: (i) for $\alpha > 2$ it fulfills $\gamma_x = \alpha$, consistent with the results of [12, 26]. However, (ii) for $1 < \alpha < 2$ we obtain a hybrid exponential and algebraic decay with a different $\gamma$ that depends linearly on $\alpha$ with a slope consistent, with $\sim 0.55(5)$ and (iii) for $\alpha \lesssim 1$ we observe numerically a curve compatible with a pure algebraic decay, with an $\alpha$-dependence of $\gamma_x$ that is linear with slope $\sim 0.25(2)$. The fitted exponent $\gamma_x$ is shown in figure 2(c).

For $\alpha \lesssim 1$ in the paramagnetic regions of the phase diagram denoted as PM2 we find that the effective central charge grows continuously with decreasing $\alpha$ from zero to a finite value that appears to be $\theta$-dependent and has a maximum of order 1 for $\alpha = 0$ and $\theta \approx \pi/2$. An example for $\theta = 0.2\pi$ is plotted in figure 1(d) (blue squares). As mentioned above, in this PM2 region, the correlation function $C_{1,R}^{xx}$ 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 $\Delta E(1)$ and $\Delta E(2)$ between the ground state and the first excited states in the odd parity sector increase with decreasing $\alpha$, as
<figure><img src="image_4.png" /></figure>

**Figure 1.** (a) Phase diagram of the LRI model equation (1) from the effective central charge $c_{\text{eff}}$ extracted from the entanglement entropy $S_{\ell}$ (see equation (6)) for a system of $L = 100$ sites, as a function of $\alpha$ and $\theta$. The phase diagram is symmetric with respect to $\theta = \pi/2$. PM1 and AM denote gapped paramagnetic and antiferromagnetic phases, respectively, with $c_{\text{eff}} = 0$. In the paramagnetic gapped region PM2, $c_{\text{eff}} > 0$ (see panel (c) blue squares) and massive edge modes (see section 4.2 and figure 7(d)) are found. The continuous black lines correspond to critical transition lines with central charge $c = 1/2$. The dashed red lines signal an increase of the central charge from $1/2$ up to $c \approx 1$ with decreasing $\alpha$ (see figure 3(b) blue triangles). (b) Phase diagram of the LRK model equation (2) from $c_{\text{eff}}$ extracted from $S_{\ell}$ in the thermodynamic limit as a function of $\alpha$ and $\theta$, and for $\epsilon = 0$. AM1 and PM are gapped phases with $c_{\text{eff}} = 0$ for $\alpha \gtrsim 1$ and $c_{\text{eff}} \neq 0$ for $\alpha \lesssim 1$ (see panel (c) green circles). AM1 shows massless edge modes with a hybrid exponential and power-law decay with distance. AM2 shows edge states with a finite mass $\Lambda_0$, which increases with decreasing $\alpha$ (see panel (c) and figure 8(a)). The continuous black lines correspond to critical transition lines with central charge $c = 1/2$. The dashed red lines signal an increase of the central charge from $1/2$ up to $c \approx 1$ with decreasing $\alpha$. (c) Effective central charge $c_{\text{eff}}$ for the LR Ising (blue squares) and the LRK models (green circles) for $\theta = 0.2\pi$ as function of $\alpha$. (d) Scaling with the system size $L$ of the mass $\Lambda_0$ of the edge modes within the AM2 phase, for $\theta = 0.75\pi$ and for several $\alpha < 1$. (e) Edge gap $\Lambda_0$ of the LRK model for a system of $L = 400$ sites and $\epsilon = 10$. In the region denoted by 'edge gap' the $\mathbb{Z}_2$ symmetry of the model is restored and the AM2 phase appears.

shown in figure 8(b), and the wave-functions of the two lowest-energy excited states $|1^-\rangle$ and $|2^-\rangle$, 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 $\alpha$ [26]. For $\alpha > 1$ the connected correlation functions $C_{1,R}^{zz}$ 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 $\alpha$. This degeneracy is due to the spontaneous breaking of the $\mathbb{Z}_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.

### 2.3.2. LR Kitaev models.

The phase diagram of the LRK model equation (2) for $\epsilon = 0$ is reported in figure 1(b), where we plot the effective central charge $c_{\text{eff}}$ defined in section 2.2.1 as a function of the angle $\theta$ and the power $\alpha$ of the decay of the pairing term. In this case the invariance of equation (2) under $\cos \theta \to -\cos \theta$ is lost for any finite $\alpha$ and the phase diagram is not symmetric around $\theta = \pi/2$.

Figure 1(b) shows that for $\alpha \gtrsim 1$ and $0 < \theta < \pi$, two phases exist that are denoted as PM and AM1 separated by two critical lines. In the limit $\alpha \to \infty$ these phases correspond to the paramagnetic and antiferromagnetic phases of the LRK model and are gapped. Consistently, figure 1(b) shows that the effective

<footnote>$^5$The presence of massless edge modes in the AM phase of the LRI model is not a sign of symmetry-protected topological order, as it is discussed for the SR Ising model in, e.g., [68].</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) $C_{1,R}^{xx}$ correlation equation (7) for $H_{Ising}$ ($\theta = 0.2\pi$ and $L = 60$), showing the hybrid exponential and power-law behavior for $\alpha \gtrsim 1$ and the purely power-law for $\alpha \lesssim 1$. The oscillations and the bending of $C_{1,R}^{xx}$ for $R \gtrsim 50$ are due to finite size effects of the DMRG simulations. (b) $\Sigma(1, R)$ correlation equation (9) for $H_{LRK}$ at $\theta = 0.2\pi$ and $L = 100$ (PM phase), showing the same hybrid behavior and the same decaying exponent as $C_{1,R}^{xx}$ for $\alpha > 1$. For $\alpha < 1$ instead the power-law tails have exponents differing from $C_{1,R}^{xx}$. (c) Decaying exponent $\gamma_x$ of the algebraic tail of the $C_{1,R}^{xx}$ correlation fitted as $1/R^{\gamma_x}$. Three different behaviors for $\gamma_x$, corresponding to the three dashed black lines are numerically observed (see section 2.3.1). (d) Decaying exponent $\gamma_\Sigma$ of the algebraic tail of the $\Sigma(1, R)$ correlation fitted as $1/R^{\gamma_\Sigma}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Effective central charge extracted from the scaling with $L$ of the von Neumann entropy (6) as function of $\theta$ for two values of $\alpha$, plotted for different system sizes $L = 60, 80, 100$. We located the critical point (indicated by red triangles) where the effective central charge does not depend on $L$. (Main panel) scaling at $\alpha = 2.0$: the critical point is $\theta_c \approx 0.302\pi$ and the central charge $c$ is compatible with $1/2$. (Inset) scaling at $\alpha = 0.7$: the critical point is $\theta_c \approx 0.376\pi$ and the central charge $c$ is compatible with $c \approx 0.66$. (b) Central charge on the critical lines of the LRI model (blue triangles) obtained with the same method as in panel (a) and of the LRK model (green diamonds) obtained by finite size scaling from the expression for entanglement entropy in equation (6).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) Integration contour adopted to evaluate the integral (25). The red dashed line is the branch cut of the square root in the denominator of the integrand in (25). (b) Decay rate 1/ξ of the exponential part of the correlation of equation (29) as function of α for two values of θ. The red lines represent the value of 1/ξ for α → ∞.</figcaption></figure>

central charge $c_{\text{eff}}$ is zero within the phases for all $\alpha \gtrsim 1$ (see figure 1(d) for an example), while $c_{\text{eff}} \neq 0$ 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 $\Sigma(i, j)$ defined in section 2.2.2. In the region denoted as AM1, $\Sigma(i, j)$ has a finite value for $|i - j| \to \infty$, while $\Sigma(i, j)$ decays for $|i - j| \to \infty$ within the PM phase. Similar to the situation in the LRI model (see above), the decay to zero of $\Sigma(i, j)$ in the PM phase shows a hybrid exponential and power-law behavior with distance. This is shown for a particular value of $\theta$ in figure 2(b), where we find numerically that the exponent $\gamma_{\Sigma}$ for the power-law tail of $\Sigma(i, j)$ equals $\gamma_z = \alpha$ when $\alpha > 1$. For $\alpha < 1$, however, the exponential part becomes numerically unobservable, and $\Sigma(i, j)$ decays essentially algebraically within the PM phase with an exponent that grows to 2 for $\alpha \to 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 density–density correlation functions. In particular, the leading contribution to the one-body correlation function $\langle a_i^\dagger a_j \rangle$ reads

$$ \langle a_R^\dagger a_0 \rangle = A_{\alpha, \theta} \cdot \frac{(-1)^R e^{-\xi R}}{\sqrt{R}} + B_{\alpha, \theta} \cdot \begin{cases} \frac{1}{R^{\alpha+1}} & \alpha > 2, \\ \frac{1}{R^{2\alpha-1}} & 1 < \alpha < 2, \\ \frac{1}{R^{2-\alpha}} & 0 < \alpha < 1, \end{cases} \quad (19) $$

where the pre-factors $A_{\alpha, \theta}$ and $B_{\alpha, \theta}$ are derived below in section 3. The algebraic part of the decay of $g_2(R)$ is instead found to be $g_2(R) \sim R^{-2\alpha}$ and $g_2(R) \sim R^{-2}$ for $\alpha > 1$ and $\alpha < 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 $\alpha \gtrsim 1$, the edge modes are massless, as expected from the SR Kitaev model. However, for $\alpha \lesssim 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 $\Lambda_0$ defined in section 2.2.3 as a function of $\alpha$ for a few values of the parameter $\epsilon$ of Hamiltonian equation (2): while for $\alpha \gtrsim 1$ the gap scales to zero with the system size $L$ as $1/L^\alpha$, for $\alpha \lesssim 1$ it remains finite and for $\alpha \sim 0$ can be of order unity. The presence of the gap removes the degeneracy of the ground-state, signaling a new phase for this class of topological LR models. This latter phase is denoted as AM2 in figure 1.

## 2.4. Critical lines

### 2.4.1. LR Kitaev models.

For the LRK models, the critical lines can be computed exactly as the Hamiltonians in equation (2) are integrable. A Fourier transform of the fermionic operators $a_j$ takes equation (2) to the form
$$H_{LRK} = \sin \theta \sum_k (a_k^\dagger a_{-k}) \left( \begin{array}{cc} g_\alpha(k) + \cot \theta & -i(1+\epsilon)f_\alpha(k) \\ i(1+\epsilon)f_\alpha(k) & -(g_\alpha(k) + \cot \theta) \end{array} \right) \left( \begin{array}{c} a_k \\ a_{-k}^\dagger \end{array} \right), \qquad (20)$$

where $a_k$ reads $a_k = \frac{1}{\sqrt{L}} \sum_j e^{ik_j a_j}$, $k = 2\pi n/L$ is the lattice momentum with $n = 0, \dots, L-1$, and the functions $f_\alpha(k)$ and $g_\alpha(k)$ read $f_\alpha(k) = \sum_l \sin(kl)/l^\alpha$ and $g_\alpha(k) = \sum_l \cos(kl)/l^\alpha$, respectively.

A Bogoliubov transformation brings equation (20) in diagonal form as

$$H_{LRK} = \sum_k \lambda_\alpha(k) (\eta_k^\dagger \eta_k - 1/2), \qquad (21)$$

with

$$\lambda_{\alpha}(k) = 2 \sin \theta \sqrt{(g_{\alpha}(k) + \cot \theta)^2 + (1 + \epsilon)^2 f_{\alpha}(k)^2}. \qquad (22)$$

Here, the new fermionic operators $(\eta_k^\dagger, \eta_{-k})$ are given in terms of the operators $(a_k^\dagger, a_{-k})$ by

$$\begin{pmatrix} \eta_k^\dagger \\ \eta_{-k} \end{pmatrix} = \begin{pmatrix} \sin \beta_k & \operatorname{icos}\beta_k \\ \operatorname{icos}\beta_k & \sin \beta_k \end{pmatrix} \begin{pmatrix} a_k^\dagger \\ a_{-k} \end{pmatrix}, \qquad (23)$$

where $\tan(2\beta_k) = (1 + \epsilon)f(k)/[g(k) + \cot \theta]$. The ground state of equation (21) is the vacuum of the $\eta_k$ fermions and has an energy density $e_0(\alpha, L) = -\sum_k \lambda_\alpha(k)/(2L)$.

The critical lines can be computed from the dispersion relation equation (22) as follows: (i) for a finite system with $L$ sites, equation (22) is zero on the line $\cot \theta + g_\alpha(0) = 0$ and on the line $\cot \theta + g_\alpha(\pi) = 0$; (ii) for a system in the thermodynamic limit $L \to \infty$ instead $f_\alpha(k) = \Re \mathrm{Li}_\alpha(e^{ik})$ and $g_\alpha(k) = \Re \mathrm{Li}_\alpha(e^{ik})$, $\mathrm{Li}_\alpha(x)$ being the polylogarithm of order $\alpha$ [69], and the critical line $\cot \theta + g_\alpha(0) = 0$ ends at the point $\theta = \pi$ and $\alpha = 1$.

For the critical line with $\theta < \pi/2$ we compute the value of the central charge $c$ by two methods: (i) by fitting the von Neumann entropy equation (6) and (ii) by studying the finite size corrections to the ground state energy density $e_0(\alpha, L)$ (see [53] for a similar model).

The results for $c$ obtained from the scaling of the von Neumann entropy equation (6) are reported in figure 3(b). For $\alpha \gtrsim 1$, we find $c = 1/2$ as expected from the SR model. For $\alpha \lesssim 1$, however, $c_{\text{eff}}$ increases up to values of order one (see red dashed line in 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 $\alpha \lesssim 1$ their oscillations mimic those of a Luttinger liquid. Here, we find a similar behavior (not shown). The increasing of $c_{\text{eff}}$, below $\alpha = 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_{\text{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 $e_0(\alpha, L)$ with the system size $L$ [53]: for a conformally invariant theory the following relation must hold [59]

$$e_0(\alpha, L) = e_\infty(\alpha) - \frac{\pi v_F c}{6L^2}, \qquad (24)$$

where $e_\infty(\alpha)$ is the energy density in the thermodynamic limit, $v_F = |2 \sin \theta \mathrm{Li}_{\alpha-1}(-1)|$ is the Fermi velocity and $c$ is the central charge of the conformal theory. We analyzed numerically $e_0(\alpha, L)$, finding that relation (24) works properly only for sufficiently large $\alpha \gtrsim 2$. This results in a value $c \simeq 1/2$ for $\alpha \gtrsim 2$, as expected from results for the SR Kitaev chain. Conversely, for $\alpha \lesssim 2$, $e_0(\alpha, L)$ 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 $\theta < \pi/2$ and $\alpha \lesssim 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_{\text{eff}}$ which tends to 1 as $\alpha$ 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 $\lambda_\alpha(k)$ for the long range Kitaev models on the critical line with $\theta > \pi/2$, for $\alpha > 2$ we find $\lambda_\alpha(k) \sim k$, for $k \to 0$. However, for $1 < \alpha < 2$ we obtain the scaling $\lambda_\alpha(k) \sim k^{\alpha-1}$. This latter scaling implies a dynamical exponent $z = \alpha - 1$ that varies continuously with $\alpha$ 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 $\theta > \pi/2$ and $\alpha < 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 $\theta < \pi/2$ conformal symmetry breaking arises even if the low-energy spectrum is linear for every $\alpha$.
### 2.4.2. LR Ising model.

We locate the critical line (for $\theta < \pi/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_{\text{eff}}$ for different system sizes $L$ from equation (6) and determine the precise values of $\alpha$ and $\theta$ 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 $\alpha$. We notice however that this method does not allow us to extract a precise value for $c_{\text{eff}}$ when $\alpha \lesssim 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 $\alpha \gtrsim 1$ (black solid lines in figure 1) $c_{\text{eff}}$ is equal to $1/2$ as expected for the central charge of the critical SR Ising model. However, for $\alpha \lesssim 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 $\alpha$ above $\alpha = 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.

## 3. 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).

The one-body correlation functions $\langle a_R^\dagger a_0 \rangle$ and $\langle a_R^\dagger a_0^\dagger \rangle$ read

$$ \langle a_R^\dagger a_0 \rangle = - \frac{\sin \theta}{2\pi} \int_0^{2\pi} dk e^{ikR} \left( \frac{g_\alpha(k) + \cot \theta}{2\lambda_\alpha(k)} \right) \quad (25) $$

and

$$ \langle a_R^\dagger a_0^\dagger \rangle = \frac{\mathrm{isin}\theta}{2\pi} \int_0^{2\pi} dk e^{ikR} \left( \frac{f_\alpha(k)}{2\lambda_\alpha(k)} \right), \quad (26) $$

respectively. In the following, we focus first on the one-body correlation function $\langle a_R^\dagger a_0 \rangle$ and come back to the anomalous correlation $\langle a_R^\dagger a_0^\dagger \rangle$ 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

$$ \langle a_R^\dagger a_0 \rangle = - \frac{1}{2\pi} \left( \int_{C_0} + \int_{C_-} + \int_{C_+} + \int_{C_{2\pi}} \right) e^{izR} G_\alpha(z) \, dz \quad (27) $$

with $G_\alpha(z) = [g_\alpha(z) + \cot \theta]/[2\sqrt{(\cot \theta + g_\alpha(z))^2 + f_\alpha^2(z)}]$ and $z = k + iy$, and where we have chosen $\epsilon = 0$. In equation (27) we have neglected the contributions from $C_-$ and $C'_-$ as they vanish when $M \to \infty$. As we explain below, the integrations over the lines $L_\pm$ and $C_{0,2\pi}$ (and thus momenta $k \simeq \pi$ and $k \simeq 0$ and $2\pi$, respectively) are responsible for the exponential and algebraic behavior observed in these models, respectively.

### 3.1. Exponential decay

The sum of the integrals on the lines $C_-$ (where $z = \pi^- + iy$) and $C_+$ (where $z = \pi^+ + iy$) of figure 4 gives

$$ I_{\alpha}^{\exp}(R) = - \frac{e^{i\pi R} e^{-\xi R}}{\pi} \int_{0}^{\infty} e^{-yR} G_{\alpha}(\pi + i(y + \xi)) \, dy. \quad (28) $$

Equation (28) displays an exponential behavior with a decay constant $1/\xi$. The appearance of this quantity is due to the square root in the denominator of $G_\alpha(z)$, yielding a branch cut from $z_1 = \pi + i\xi$ to $\infty$. The leading term in equation (28) is obtained by integrating $G_\alpha(\pi + i(y + \xi))$ in the limit $y \to 0$ [75], and reads
$$I_{\alpha}^{\exp}(R) = (-1)^{R} A_{\alpha, \theta} \frac{e^{-\xi R}}{\sqrt{R}}, \qquad (29)$$

with

$$A_{\alpha, \theta} = - \frac{\left( 2 \cot \theta + \mathrm{Li}_{\alpha}(-e^{\xi}) + \mathrm{Li}_{\alpha}(-e^{-\xi}) \right)}{4\sqrt{\pi} \left| \mathrm{Li}_{\alpha-1}(-e^{\xi}) \right|^{1/2} \left[ \cot \theta + \mathrm{Li}_{\alpha}(-e^{-\xi}) \right]^{1/2}}. \qquad (30)$$

The decay constant $\xi$ is related to the zeroes of the denominator of $G_{\alpha}(z)$ and is obtained by solving the equation

$$[\cot \theta + \mathrm{Li}_{\alpha}(-e^{-\xi})][\cot \theta + \mathrm{Li}_{\alpha}(-e^{\xi})] = 0. \qquad (31)$$

Two cases must be distinguished: if $\cot \theta < |\mathrm{Li}_{\alpha}(-1)|$, the equation $[\cot \theta + \mathrm{Li}_{\alpha}(-e^{-\xi})] = 0$ admits a solution, since the function $-\mathrm{Li}_{\alpha}(-e^{-\xi})$ is always decreasing for $\xi > 0$. If instead $\cot \theta > |\mathrm{Li}_{\alpha}(-1)|$, then $[\cot \theta + \mathrm{Li}_{\alpha}(-e^{\xi})] = 0$ admits a solution, for the same reason as above. In the following we focus, without loss of generality, on the first case, where $\xi$ is solution of $\cot \theta = -\mathrm{Li}_{\alpha}(-e^{-\xi})$. Notably for $\alpha = 0$ solutions exist only for $\theta > \pi/4$ [in which case, $\xi = \pm \log(\tan \theta - 1)$ for $\cot \theta \le 1/2$]. For $\theta < \pi/4$, instead, 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 $1/\xi$ as a function of $\alpha$ for two different values of $\theta$. In particular, for $\theta = 0.2\pi$ and $\alpha \to \infty$, $1/\xi$ tends to the SR value ($\xi \to |\log|\tan \theta||$), as expected. However, for $\alpha \to 0$ we find that $1/\xi$ tends to zero, essentially linearly with $\alpha$. As explained below, this can result in the non-observability of the exponential dependence of correlation functions for $\alpha \lesssim 1$. Notice however that even if $\xi$ is finite when $\alpha \lesssim 1$ for $\theta = 0.4\pi$, the exponential part of the correlation functions is still unobservable.

## 3.2. Algebraic decay

The sum of the integrals on the lines $C_0$ (where $z = \eta + iy$) and $C_{2\pi}$ (where $z = 2\pi - \eta + iy$) of figure 4 gives

$$I_{\alpha}^{\mathrm{pow}}(R) = \frac{1}{\pi} \int_{0}^{\infty} e^{-yR} \mathfrak{I}G_{\alpha}(iy) dy, \qquad (32)$$

after sending $\eta \to 0$. The leading contribution for $R \to \infty$ to the integral in equation (32) can be computed again by integrating the imaginary part $\mathfrak{I}G_{\alpha}(iy)$ for $y \to 0$. By exploiting the following series expansion of the polylogarithm [69]

$$\mathrm{Li}_{\alpha}(e^{\pm y}) = \Gamma(1-\alpha)(\mp y)^{\alpha-1} + \sum_{j=0}^{\infty} \frac{\zeta(\alpha-j)}{j!} (\pm y)^j, \qquad (33)$$

one obtains

$$I_{\alpha}^{\mathrm{pow}}(R) = B_{\alpha, \theta} \cdot \begin{cases} \frac{1}{R^{\alpha+1}} & \alpha > 2, \\ \frac{1}{R^{2\alpha-1}} & 1 < \alpha < 2, \\ \frac{1}{R^{2-\alpha}} & \alpha < 1 \end{cases} \qquad (34)$$

with

$$B_{\alpha, \theta} = \begin{cases} \frac{\alpha \zeta(\alpha - 1)}{(\zeta(\alpha) + \cot \theta)^2} & \alpha > 2, \\ \frac{(\cos(\pi\alpha) - 1)\Gamma(1 - \alpha)\Gamma(2\alpha - 1)}{4\Gamma(\alpha)(\cot(\theta) + \zeta(\alpha))^2} & 1 < \alpha < 2, \\ \frac{\cos^3\left(\frac{\pi\alpha}{2}\right)\Gamma(2 - \alpha)(\cot(\theta) + \zeta(\alpha))}{\pi \Gamma(1 - \alpha)} & \alpha < 1. \end{cases} \qquad (35)$$

While, e.g., the phase diagram of figure 1(b) demonstrates the persistence of individual paramagnetic and antiferromagnetic phases with varying $\alpha$, 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.
<figure><img src="image_8.png" /><figcaption>Figure 5. (a) Comparison between numerical solutions of equation (25) (symbols) and the analytic expressions equations (29) and (34) (solid lines) for $\theta = 0.239\pi$ and different $\alpha$. The inset refers to $\alpha = 0.45$, where the exponential part is negligible. (b) $C_{R_0+1,R_0+R}^{zz}$ correlation for $H_{LR1}$ at $\theta = 0.207\pi$, $L = 100$, $R_0 = L/4$ and different $\alpha$. The power-law part of $C_{R_0+1,R_0+R}^{zz}$ shows the same decay exponent $\gamma_z = 2\alpha$ for $\alpha > 1$, as $g_2(0, R)$ for $H_{LRK}$ (equation (39), solid black lines).</figcaption></figure>

### 3.3. 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 $\alpha$. Figure 5(a) shows that the analytical results are in perfect agreement with a numerical solution of equation (25).

The anomalous correlation $\langle a_R^\dagger a_0^\dagger \rangle$ can be computed along the same lines as before and is given by

$$ \langle a_R^\dagger a_0^\dagger \rangle = A_{\alpha, \theta}^{(a)} \frac{(-1)^R e^{-\xi R}}{\sqrt{R}} + B_{\alpha, \theta}^{(a)} \begin{cases} \frac{1}{R^\alpha} & \alpha > 1, \\ \frac{1}{R} & 0 < \alpha < 1 \end{cases} \quad (36) $$

with

$$ A_{\alpha, \theta}^{(a)} = \frac{\mathrm{Li}_{\alpha}(-e^{-\xi}) - \mathrm{Li}_{\alpha}(-e^{\xi})}{4\sqrt{\pi} |\mathrm{Li}_{\alpha-1}(-e^{\xi})|^{1/2} (\cot \theta + \mathrm{Li}_{\alpha}(-e^{-\xi}))^{1/2}} \quad (37) $$

and

$$ B_{\alpha, \theta}^{(a)} = \begin{cases} -\frac{1}{4(\cot \theta + \zeta(\alpha))} & \alpha > 1, \\ -\frac{1}{2\pi} \cos \frac{\pi\alpha}{2} & 0 < \alpha < 1. \end{cases} \quad (38) $$

From $\langle a_R^\dagger a_0 \rangle$ and $\langle a_R^\dagger a_0^\dagger \rangle$ we can compute other local correlations by Wick's theorem. For instance, the density-density correlation function $g_2(i, j)$ reads $g_2(i, j) = |\langle a_i^\dagger a_j^\dagger \rangle|^2 - |\langle a_i^\dagger a_j \rangle|^2$ and, from equations (34) and (36), the leading part of its LR power law tail is found to be

$$ g_2(0, R) \sim \begin{cases} \frac{1}{R^{2\alpha}} & \alpha > 1, \\ \frac{1}{R^2} & 0 < \alpha < 1. \end{cases} \quad (39) $$

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 $k = \pi$ (as for the SR limit) and $k = 0, 2\pi$ 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 $\mathrm{Li}_\alpha(e^y)$ and disappears in the limit $\alpha \to \infty$. This implies a simple exponential decay of correlations as
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a): $C_{R_0+1,R_0+R}^{zz}$ correlation for the LRI for $\theta = 0.35\pi$ (AM phase) and $L = 200$, computed from $R_0 = 100$. The data are fitted by the solid lines $\sim 1/R^\gamma$, $\gamma_z$ reported in the plot. (b): $C_{1,R}^{xx}$ correlation at $\theta = 0.207\pi$ (PM phases) for a system of $L = 30$ sites.</figcaption></figure>

expected for a SR model. Conversely, (ii) the exponential decay is negligible when $\alpha \lesssim 1$, as shown in the inset of figure 5(a), implying an essentially pure algebraic decay.

### 3.4. Comparison with correlations of the LRI model

In this section we compare the decay with distance of the correlation functions $g_2(i, j)$ and $\Sigma(i, j)$ of the LRK models with those of the correlations $C_{i,j}^{zz}$ and $C_{i,j}^{xx}$ of the LRI model, respectively, since they are related by a Jordan-Wigner transformation (see, e.g., section 2.2.2). The correlation functions $C_{i,j}^{zz}$ and $C_{i,j}^{xx}$ have been studied within the PM phase in [26] for $\alpha \gtrsim 0.5$. There, for $\alpha \gtrsim 1$ it was found that the long-distance behavior is characterized by an algebraic decay with exponents $\gamma_z = 2\alpha$ and $\gamma_x = \alpha$ for the two correlations, respectively.

For $\alpha \gtrsim 1$, our own calculations for $C_{i,j}^{xx}$ and $C_{i,j}^{zz}$ are reported in figures 2(b) and 5(b), respectively. There, we compare the decay of correlations with that of the corresponding correlation functions in the LRK models, $g_2(i, j)$ (see also equation (39)) and $\Sigma(i, j)$, respectively, showing very good quantitative agreement with the (semi-)analytic results. In agreement with [12, 53], an hybrid exponential and power-law behavior is found for $C_{i,j}^{xx}$ (as well as for $\Sigma(i, j)$ and $g_2(i, j)$).

Conversely, when $\alpha \lesssim 1$ only an algebraic tail is visible for the decay of $C_{i,j}^{zz}$ in figure 5 and $C_{i,j}^{xx}$ [as well as $\Sigma(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 $\Sigma(i, j)$ differ here in general, probably because the contribution of the string operators in $H_{\text{LRI}}$ becomes more relevant.

In the AM phase, instead, we find that the decay of $C_{i,j}^{zz}$ (figure 6(a)) for $\alpha \gtrsim 1$ displays an algebraic tail with an exponent compatible with $\gamma_z = \alpha$, 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 $\alpha = 3$. Notably the decay exponent is here always different from the value $2\alpha$ analytically computed for $g_2(i, j)$ in equation (39). This discrepancy is again probably due to the role of the string operators, here quantitatively relevant even for $\alpha \gtrsim 1$. For $\alpha \lesssim 1$ again no universal behavior for $\gamma_z$ is found.

## 4. 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 $\pi/4 < \theta < 3\pi/4$. These modes are fermionic and massless (in the limit $L \to \infty$), 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 $\alpha \gtrsim 1$ in the AM1 region of the phase diagram of figure 1(b). Plotting the square of the wavefunction $\psi_{0j}$ corresponding to the zero edge gap $\Lambda_0$, defined in section 2.2.3, we find for $\epsilon \neq 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 $|\psi_{0j}|^2$ is found to be equal to $2\alpha$. This behavior is similar to that observed in [53] in the presence of LR pairing only.
<figure><img src="image_10.png" /></figure>

**Figure 7.** (a) LRK model: spatial probability distribution for the Majorana edge mode $\psi_{0j}$ in the AM1 phase for $\theta = 0.35\pi$, $\alpha = 3$ and different $\epsilon$. When $\epsilon = 1$, the mode is exponentially localized at one edge, while a power-law tail gradually appears when $\epsilon \neq 1$. (b) LRI model: spatial probability distribution for the mode $\varphi_j^{(1)}$ as function of the lattice site $j$ in the AM phase for different $\alpha$ and $\theta = 0.45\pi$. The decay with $j$ is numerically found to be exponential. (c) LRK model: spatial probability distribution for the massive edge mode in the AM2 phase for $\theta = 0.7\pi$ for different $\alpha$ and $\epsilon$. The decay with distance is purely power law and the probability distribution is symmetric with respect to half of the chain. (d) LRI model: probability density $p_j = |w_j^{(1)}|^2 + |w_j^{(2)}|^2$ of the two degenerate excited states $|1^-\rangle$ and $|2^-\rangle$ in the PM2 phase. When $\alpha \gtrsim 1$, $p_j$ spreads into the bulk of the chain. When $\alpha \lesssim 1$, $p_j$ is exponentially localized near a single edge. In panel (b) and (d), we plot only the points compatible with the DRMG errors.

Interestingly, here the algebraic tail of the Majorana modes can be tuned to completely disappear by changing the parameter $\epsilon$ that fixes the unbalance between the hopping and the pairing terms in the Hamiltonian of equation (2). Figure 7(a) shows an example for $\theta = 0.4\pi$ and $\alpha = 3$. When $\epsilon \neq 1$ the hybrid behavior is fully visible and the state $\psi_{0j}$ decays in the bulk with an algebraic tail. However, by approaching the value $\epsilon = 1$ the algebraic tail of $\psi_{0j}$ decreases and eventually disappears. As a result, for $\epsilon = 1$ the wave function $\psi_{0j}$ becomes exponentially localized at one edge. We find that this exponential localization for $\epsilon = 1$ is present also in the parameter region $\alpha < 1$. Moreover, the edge gap $\Lambda_0$ scales to zero exponentially with the system size $L$ for all $\alpha$ as figure 8(c) shows.

For the LR Ising Hamiltonian of equation (1), edge modes appear in the AM phase for every $\alpha > 0$. They have zero mass, since the edge gap $E_{|0^-\rangle} - E_{|0^+\rangle}$ (defined in section 2.2.3) vanishes. Examples of these edge modes for different $\alpha$ are given in figure 7(b), where we plot the square of the wavefunction $\varphi_j^{(1)}$ defined in equation (12). We find numerically that $\varphi_j^{(1)}$ decays exponentially with the distance from the edge of the chain for all values of $\alpha$, as well as $\varphi_j^{(2)}$ (not shown) from the opposite edge.

## 4.2. Massive edge modes

If we extend the analysis of the LR Kitaev Hamiltonian of equation (2) to different $\epsilon$ and sufficiently small $\alpha$, a totally new situation arises for the edge gap $\Lambda_0$ and the edge modes: in the region denoted as AM2 in the phase diagram of figure 1(b), we find that $\Lambda_0$ (which is zero for $\alpha \gtrsim 1$) becomes nonzero for $\alpha \lesssim 1$ and $\epsilon \neq 1$, also in the thermodynamical limit.

This case is shown in figure 1(e), where we plot the edge gap $\Lambda_0$ as function of $\alpha$ and $\theta$ for $\epsilon = 10$. Between $\epsilon = 1$ and $\epsilon = 10$ we checked a continuous increase of the extension of the region AM2 and no transitions in
<figure><img src="image_11.png" /></figure>

**Figure 8.** (a) First 45 low energy states of the LRK models for $L = 500$ and $\theta = 0.7\pi$ as function of $\alpha$. It is possible to distinguish the $\Lambda_0$ mode (solid points for three different $\epsilon$) separated from the $\Lambda_{n>0}$ bulk modes (light green triangles for $\epsilon = 10$). When $\alpha \gtrsim 1$ the vanishing of $\Lambda_0$ for $\epsilon \neq 1$ signals the spontaneous breaking of the $\mathbb{Z}_2$ symmetry of LRK and defines the AM1 phase (see also the phase diagram in figure 1(b)). When $\Lambda_0 > 0$ the $\mathbb{Z}_2$ symmetry is restored and the phase AM2 appears. (b) Main panel: $\Delta E(1)$ (empty circles) and $\Delta E(2)$ (solid triangles) for the LRI as functions of the inverse of the length of the chain $L$ for $\theta = 0.159\pi$ and for three $\alpha$ reported in the plot. These two states are degenerate in the $L \to \infty$ limit, as finite-size scaling shows. Inset: $\Delta E(1)$ as function of $\alpha$ for $\theta = 0.159\pi$. When $\alpha$ decreases, $\Delta 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 $\Lambda_0$ for (squares) $\epsilon = 0.5$ and two different $\alpha$ as function of the system size $L$. When $\epsilon = 1$, the edge gap scales exponentially to zero for all $\alpha$. For comparison, the case (crosses) $\epsilon = 2$ and $\alpha = 3$ is also reported. There the edge mass displays an hybrid exponential and power-law decay.

between. Similarly in figure 8(a) we plot $\Lambda_0$ together with several other single-particle energies as a function of $\alpha$ and for a fixed $\theta$. These energies have been computed as described in section 2.2.3. In figure 8(a) the mass of the edge mode $\Lambda_0$ is easily recognizable for all $\epsilon$, since it is separated from all bulk modes by a finite gap.

Consistently with the discussion above, for $\alpha \gtrsim 1$, $\Lambda_0$ is zero as expected from the SR model, so that two degenerate ground states exist as the $\mathbb{Z}_2$ symmetry of the model is spontaneously broken. However, surprisingly, for $\alpha \lesssim 1$ and $\epsilon \neq 1$, we find that $\Lambda_0$ becomes finite and thus the ground state is unique. This indicates that the $\mathbb{Z}_2$ symmetry of the model (which is broken for $\alpha \gtrsim 1$ in the AM1 phase) is restored for $\alpha \lesssim 1$. As a consequence, the region AM2, where the $\mathbb{Z}_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 $\Lambda_0$ state is now given by the matrix element $g_{0j}$ defined in section 2.2.3. By plotting the probability density $|g_{0j}|^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 $\alpha \lesssim 1$, where the $\mathbb{Z}_2$ symmetry is preserved. However, while for the LRK massive edge modes originate from Majorana edge modes present at $\alpha \gtrsim 1$, for the LRI the edge localization arises for excited states of the bulk spectrum in the region PM1. For $\alpha \lesssim 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 $p_j = |w_j^{(1)}|^2 + |w_j^{(2)}|^2$, with $w_j^{(1,2)}$ wavefunctions of $|(1, 2)\rangle$ 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 $\alpha$. For $\alpha \gtrsim 1$, $p_j$ is oscillating and delocalized in the bulk, while it is localized exponentially at the edges for $\alpha \lesssim 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 $\mathbb{Z}_2$ symmetry, similar to the LRK models above.

## 5. Observability in current experiments

Recent experiments with cold ions have made possible the realization of LR Ising-type Hamiltonians as equation (1) with $0 \lesssim \alpha \lesssim 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 $\alpha \lesssim 1$ could allow for the precise determination of the properties of these LR models.
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_{1,R}^{xx}$ for Hamiltonian equation (1) in the PM1 and PM2 phases for a chain of $L = 30$ sites with open boundary conditions and for different $\alpha$. For $\alpha > 2$, the initial exponential decay dominates the correlations for $R \lesssim 10$, while a comparatively small algebraic tail is found for $R \gtrsim 10$. For $\alpha \lesssim 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 $\alpha \approx 1$ may be observable. We note, however, that the exponent $\gamma_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_{1,R}^{zz}$.

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 $\alpha \lesssim 1$.

## 6. 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 $\alpha$. We have further demonstrated the breaking of conformal symmetry along the critical lines in both models at low enough $\alpha$. 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 $\alpha \lesssim 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 $\alpha \gtrsim 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 $\alpha \lesssim 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 $\alpha \lesssim 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 $\alpha$, due to the bulk overlap between their wave functions, whose decay is slower and slower for decreasing $\alpha$. 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.

## Acknowledgments

We are pleased to thank Alessio Celi, Domenico Giuliano, Alexey Gorshkov, Miguel Angel Martin-Delgado, Andrea Trombettoni and Oscar Viyuela Garcia for fruitful discussions and Fabio Ortolani for help with the DMRG code. We acknowledge support by the ERC-St Grant ColdSIM (No 307688), EOARD, UdS via Labex NIE and IdEX, RYSQ, computing time at HPC-UdS.

## References

[1] Hsieh D, Qian D, Wray L, Xia Y, Hor Y S, Cava R J and Hasan M Z 2008 *Nature* **452** 970–4

[2] Hafezi M, Mittal S, Fan J, Migdall A and Taylor J M 2013 *Nat. Photon.* **7** 1001–5

[3] Jotzu G, Messer M, Desbuquois R, Lebrat M, Uehlinger T, Greif D and Esslinger T 2014 *Nature* **515** 237–40

[4] Mourik V, Zuo K, Frolov S M, Plissard S R, Bakkers E P A M and Kouwenhoven L P 2012 *Science* **336** 1003

[5] Franz M 2013 *Nat. Nanotechnol.* **8** 149–52

[6] Nadj-Perge S, Drozdov I K, Li J, Chen H, Jeon S, Seo J, MacDonald A H, Bernevig B A and Yazdani A 2014 *Science* **346** 602–7

[7] Deng M T, Yu C L, Huang G Y, Larsson M, Caroff P and Xu H Q 2012 *Nano Lett.* **12** 6414–9

[8] Das A, Ronen Y, Most Y, Oreg Y, Heiblum M and Shtrikman H 2012 *Nat. Phys.* **8** 887–95

[9] Rokhinson L P, Liu X and Furdyna J K 2012 *Nat. Phys.* **8** 795–9

[10] Finck A D K, Van Harlingen D J, Mohseni P K, Jung K and Li X 2013 *Phys. Rev. Lett.* **110** 126406

[11] Kitaev A Y 2001 *Phys.—Usp.* **44** 131

[12] Deng X L, Porras D and Cirac J I 2005 *Phys. Rev. A* **72** 063407
[13] Schneider C, Porras D and Schaetz T 2012 *Rep. Prog. Phys.* **75** 024401

[14] Friedenauer A, Schmitz H, Glueckert J T, Porras D and Schaetz T 2008 *Nat. Phys.* **4** 757–61

[15] Britton J W, Sawyer B C, Keith A C, Wang C C J, Freericks J K, Uys H, Biercuk M J and Bollinger J J 2012 *Nature* **484** 489–92

[16] Jurcevic P, Lanyon B P, Hauke P, Hempel C, Zoller P, Blatt R and Roos C F 2014 *Nature* **511** 202–5

[17] Bermudez A, Schaetz T and Plenio M B 2013 *Phys. Rev. Lett.* **110** 110502

[18] Islam R, Senko C, Campbell W C, Korenblit S, Smith J, Lee A, Edwards E E, Wang C C J, Freericks J K and Monroe C 2013 *Science* **340** 583–7

[19] Hauke P and Tagliacozzo L 2013 *Phys. Rev. Lett.* **111** 207202

[20] Richerme P, Gong Z X, Lee A, Senko C, Smith J, Foss-Feig M, Michalakis S, Gorshkov A V and Monroe C 2014 *Nature* **511** 198–201

[21] Gong Z X, Foss-Feig M, Michalakis S and Gorshkov A V 2014 *Phys. Rev. Lett.* **113** 030602

[22] Foss-Feig M, Gong Z X, Clark C W and Gorshkov A V 2015 *Phys. Rev. Lett.* **114** 157201

[23] Cevolani L, Carleo G and Sanchez-Palencia L 2015 *Phys. Rev. A* **92** 041603

[24] Schachenmayer J, Lanyon B P, Roos C F and Daley A J 2013 *Phys. Rev. X* **3** 031015

[25] Jurcevic P, Hauke P, Maier C, Hempel C, Lanyon B P, Blatt R and Roos C F 2015 *Phys. Rev. Lett.* **115** 100501

[26] Koffel T, Lewenstein M and Tagliacozzo L 2012 *Phys. Rev. Lett.* **109** 267203

[27] Hazzard K R A, van den Worm M, Foss-Feig M, Manmana S R, Dalla Torre E G, Pfau T, Kastner M and Rey A M 2014 *Phys. Rev. A* **90** 063622

[28] Pientka F, Glazman L I and von Oppen F 2013 *Phys. Rev. B* **88** 155420

[29] Pientka F, Glazman L I and von Oppen F 2014 *Phys. Rev. B* **89** 180505

[30] Duan L M, Demler E and Lukin M D 2003 *Phys. Rev. Lett.* **91** 090402

[31] García-Ripoll J J, Martin-Delgado M A and Cirac J I 2004 *Phys. Rev. Lett.* **93** 250405

[32] Baranov M A, Osterloh K and Lewenstein M 2005 *Phys. Rev. Lett.* **94** 070404

[33] Cooper N R, Rezayi E H and Simon S H 2005 *Phys. Rev. Lett.* **95** 200402

[34] Micheli A, Brennen G K and Zoller P 2006 *Nat. Phys.* **2** 341–7

[35] Brennen G K, Micheli A and Zoller P 2007 *New J. Phys.* **9** 138

[36] Lahaye T, Menotti C, Santos L, Lewenstein M and Pfau T 2009 *Rep. Prog. Phys.* **72** 126401

[37] Cooper N R and Shlyapnikov G V 2009 *Phys. Rev. Lett.* **103** 155302

[38] Weimer H, Muller M, Lesanovsky I, Zoller P and Buchler H P 2010 *Nat. Phys.* **6** 382–8

[39] Levinsen J, Cooper N R and Shlyapnikov G V 2011 *Phys. Rev. A* **84** 013603

[40] Baranov M A, Dalmonte M, Pupillo G and Zoller P 2012 *Chem. Rev.* **112** 5012–61

[41] Yao N Y, Laumann C R, Gorshkov A V, Bennett S D, Demler E, Zoller P and Lukin M D 2012 *Phys. Rev. Lett.* **109** 266804

[42] Peter D, Griesmaier A, Pfau T and Büchler H P 2013 *Phys. Rev. Lett.* **110** 145303

[43] Yao N Y, Gorshkov A V, Laumann C R, Läuchli A M, Ye J and Lukin M D 2013 *Phys. Rev. Lett.* **110** 185302

[44] Manmana S R, Stoudenmire E M, Hazzard K R A, Rey A M and Gorshkov A V 2013 *Phys. Rev. B* **87** 081106

[45] Gorshkov A V, Hazzard K R and Rey A M 2013 *Mol. Phys.* **111** 1908–16

[46] Maghrebi M F, Yao N Y, Hafezi M, Pohl T, Firstenberg O and Gorshkov A V 2015 *Phys. Rev. A* **91** 033838

[47] Yao N Y, Bennett S D, Laumann C R, Lev B L and Gorshkov A V 2015 *Phys. Rev. A* **92** 033609

[48] Cohen I and Retzker A 2014 *Phys. Rev. Lett.* **112** 040503

[49] Senko C, Richerme P, Smith J, Lee A, Cohen I, Retzker A and Monroe C 2015 *Phys. Rev. X* **5** 021026

[50] Gong Z X, Maghrebi M F, Hu A, Wall M L, Foss-Feig M and Gorshkov A V 2015 arXiv:1505.03146

[51] White S R 1992 *Phys. Rev. Lett.* **69** 2863–6

[52] Schollwöck U 2005 *Phys. Rev. Lett.* **77** 259–315

[53] Vodola D, Lepori L, Ercolessi E, Gorshkov A V and Pupillo G 2014 *Phys. Rev. Lett.* **113** 156402

[54] Lieb E, Schultz T and Mattis D 1961 *Ann. Phys.*, NY **16** 407

[55] Amico L, Fazio R, Osterloh A and Vedral V 2008 *Rev. Mod. Phys.* **80** 517

[56] Eisert J, Cramer M and Plenio M 2010 *Rev. Mod. Phys.* **82** 277

[57] Calabrese P and Cardy J 2004 *J. Stat. Mech.* P06002

[58] di Francesco P, Mathieu P and Senechal D 1997 *Conformal Field Theory* (Berlin: Springer)

[59] Henkel M 1999 *Conformal Invariance and Critical Phenomena* (New York: Springer)

[60] Truong T T and Peschel I 1989 Z. *Phys.* B **75** 119

[61] Peschel I, Kaulke M and Legeza O 1999 *Ann. Phys.* (Leipzig) **8** 153

[62] Peschel I 2003 *J. Phys. A: Math. Gen.* **36** L205

[63] Peschel I 2012 *Braz. J. Phys.* **42** 267

[64] Nayak C, Simon S H, Stern A, Freedman M and Das Sarma S 2008 *Rev. Mod. Phys.* **80** 1083–159

[65] Stoudenmire E M, Alicea J, Starykh O A and Fisher M P 2011 *Phys. Rev. B* **84** 014503

[66] Chan Y H, Chiu C K and Sun K 2015 *Phys. Rev. B* **92** 104514

[67] Mussardo G 2010 *An Introduction to Exactly Solved Models in Statistical Physics Statistical Field Theory* (New York: Oxford University Press)

[68] Greiter M, Schnells V and Thomale R 2014 *Ann. Phys.*, NY **351** 1026–33

[69] Olver F W J, Lozier D W, Boisvert R F and Clark C W 2010 *NIST Handbook of Mathematical Functions* (Cambridge: Cambridge University Press)

[70] Ares F, Esteve J G, Falceto F and de Queiroz A R 2015 *Phys. Rev. A* **92** 042334

[71] Maghrebi M F, Gong Z X and Gorshkov A V 2015 arXiv:1510.01325

[72] Gong Z X, Maghrebi M F, Hu A, Foss-Feig M, Richerme P, Monroe C and Gorshkov A V 2015 arXiv:1510.02108

[73] Campos Venuti L, Degli Esposti Boschi C, Roncaglia M and Scaramucci A 2006 *Phys. Rev. A* **73** 010303

[74] Roncaglia M, Campos Venuti L and Degli Esposti Boschi C 2008 *Phys. Rev. B* **77** 155413

[75] Ablowitz M J and Fokas A T S 2003 *Complex Variables Introduction and Applications* (Cambridge: Cambridge University Press)

[76] Wilczek F 2009 *Nat. Phys.* **5** 614–8

[77] Fidkowski L and Kitaev A 2011 *Phys. Rev. B* **83** 075103