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

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

$$\begin{eqnarray}&&{H}_{\mathrm{LRI}}=\mathrm{sin}\theta \displaystyle \sum _{i=1;j\gt i}^{L}\displaystyle \frac{{\sigma }_{i}^{x}{\sigma }_{j}^{x}}{{| i-j| }^{\alpha }}+\mathrm{cos}\theta \displaystyle \sum _{i=1}^{L}{\sigma }_{i}^{z},\end{eqnarray} \tag{ 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 $\mathrm{sin}\theta \gt 0$ (or equivalently $0\lt \theta \lt \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 }_{{\rm{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 α and θ.

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

$$\begin{eqnarray}&&{H}_{\mathrm{LRK}}=\mathrm{sin}\theta \displaystyle \sum _{i,j}\displaystyle \frac{{a}_{i}^{\dagger }{a}_{j}+(1+\epsilon )\;{a}_{i}{a}_{j}+{\rm{h}}.{\rm{c}}.}{{| i-j| }^{\alpha }}+2\mathrm{cos}\theta \displaystyle \sum _{i}{n}_{i},\end{eqnarray} \tag{ 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 α. Here, all energies are expressed in dimensionless units and the parameter 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]

$$\begin{eqnarray}&&{\sigma }_{j}^{+}={a}_{j}^{\dagger }\;\mathrm{exp}({\rm{i}}\;\pi \displaystyle \sum _{{\ell }\lt j}{n}_{{\ell }});\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\sigma }_{j}^{z}=2{n}_{j}-1,\end{eqnarray} \tag{ 4 }$$

with ${\sigma }_{j}^{+}=({\sigma }_{j}^{x}+{\rm{i}}{\sigma }_{j}^{y})/2$. However, at finite α this identification does not hold anymore due to the contributions of the string operators $\mathrm{exp}({\rm{i}}\pi {\displaystyle \sum }_{{\ell }\lt j}{n}_{{\ell }})$ in equation (4). In particular, unlike the LRI, equation (2) remains exactly solvable, allowing for analytic solutions at any finite α.

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 ${{ \mathcal 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 ℓ and L − ℓ sites, respectively, ${{ \mathcal S }}_{{\ell }}$ is defined as

$$\begin{eqnarray}&&{{ \mathcal S }}_{{\ell }}=-\mathrm{Tr}\;{\rho }_{{\ell }}{\mathrm{log}}_{2}{\rho }_{{\ell }},\end{eqnarray} \tag{ 5 }$$

where ${\rho }_{{\ell }}$ is the reduced density matrix of the subsystem A.

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

$$\begin{eqnarray}&&{{ \mathcal S }}_{{\ell }}=a+\displaystyle \frac{c}{6}{\mathrm{log}}_{2}d({\ell },L),\end{eqnarray} \tag{ 6 }$$

with $d({\ell },L)=(2L/\pi )\mathrm{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]. ${{ \mathcal 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 ${{ \mathcal 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}_{{\rm{eff}}}$ may be defined also within the gapped phases and used to characterize the main features of the phase diagram for LR interactions [26].

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

$$\begin{eqnarray}&&{C}_{i,j}^{\nu \nu }=\langle {\sigma }_{i}^{\nu }{\sigma }_{j}^{\nu }\rangle -\langle {\sigma }_{i}^{\nu }\rangle \langle {\sigma }_{j}^{\nu }\rangle ,\end{eqnarray} \tag{ 7 }$$

with $\nu =x,z$.

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

$$\begin{eqnarray}&&{g}_{2}(i,j)=\langle {n}_{i}{n}_{j}\rangle -\langle {n}_{i}\rangle \langle {n}_{j}\rangle \;\end{eqnarray} \tag{ 8 }$$

and the function

$$\begin{eqnarray}&&{\rm{\Sigma }}(i,j)=\langle ({a}_{i}^{\dagger }+{a}_{i})\mathrm{exp}(\pi {\rm{i}}\displaystyle \sum _{{\ell }=i+1}^{j-1}{a}_{{\ell }}^{\dagger }{a}_{{\ell }})({a}_{j}^{\dagger }+{a}_{j})\rangle \end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{g}_{2}(i,j)={| \langle {a}_{i}^{\dagger }{a}_{j}^{\dagger }\rangle | }^{2}-{| \langle {a}_{i}^{\dagger }{a}_{j}\rangle | }^{2}\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}{\rm{\Sigma }}(i,j)=\mathrm{det}(\begin{array}{ccc}{G}_{i,i+1} & ... & {G}_{i,j}\\ \vdots & \ddots & \\ {G}_{j-1,i+1} & ... & {G}_{j-1,j}\end{array})\end{eqnarray} \tag{ 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.

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 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 ${c}_{j}^{\dagger }={\sigma }_{j}^{+}\mathrm{exp}({\rm{i}}\pi {\displaystyle \sum }_{{\ell }\lt j}{\sigma }_{{\ell }}^{+}{\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]

$$\begin{eqnarray}&&{\varphi }_{j}^{(1)}=\langle {0}^{-}| {c}_{j}^{\dagger }| {0}^{+}\rangle +\langle {0}^{+}| {c}_{j}^{\dagger }| {0}^{-}\rangle ,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\varphi }_{j}^{(2)}={\rm{i}}(\langle {0}^{-}| {c}_{j}^{\dagger }| {0}^{+}\rangle -\langle {0}^{+}| {c}_{j}^{\dagger }| {0}^{-}\rangle ).\end{eqnarray} \tag{ 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}^{\pm }\rangle }$ being the energy of the state $| {0}^{\pm }\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

$$\begin{eqnarray}&&{w}_{j}^{(1)}=\langle {1}^{-}| {c}_{j}^{\dagger }| {0}^{+}\rangle \qquad {w}_{j}^{(2)}=\langle {2}^{-}| {c}_{j}^{\dagger }| {0}^{+}\rangle .\end{eqnarray} \tag{ 14 }$$

Here the mass of the mode ${w}_{j}^{(s)}$ is defined as the difference ${\rm{\Delta }}E(s)={E}_{| {s}^{-}\rangle }-{E}_{| {0}^{+}\rangle }$ 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

$$\begin{eqnarray}&&H=\displaystyle \sum _{i,j=1}^{L}[{A}_{{ij}}{a}_{i}^{\dagger }{a}_{j}+\displaystyle \frac{1}{2}({a}_{i}^{\dagger }{B}_{{ij}}{a}_{j}^{\dagger }+{\rm{h}}.{\rm{c}}.)],\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}&&H=\displaystyle \sum _{n=0}^{L-1}{{\rm{\Lambda }}}_{n}{\eta }_{n}^{\dagger }{\eta }_{n}.\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{n}=\displaystyle \sum _{{ij}}{\psi }_{{nj}}({A}_{{ji}}+{B}_{{ji}}){\phi }_{{in}}.\end{eqnarray} \tag{ 17 }$$

Here, ${{\rm{\Lambda }}}_{n}$ are single-particle energies (ordered as ${{\rm{\Lambda }}}_{0}\leqslant {{\rm{\Lambda }}}_{1}\leqslant \cdots \leqslant {{\rm{\Lambda }}}_{L-1}$) and ${\eta }_{n}$ are fermionic operators defined by the following Bogoliubov transformation

$$\begin{eqnarray}&&{\eta }_{n}=\displaystyle \sum _{j}{g}_{{nj}}{a}_{j}+{h}_{{nj}}{a}_{j}^{\dagger }.\end{eqnarray} \tag{ 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 ${\rm{i}}({\eta }_{n}-{\eta }_{n}^{\dagger })/2$ with energy ${{\rm{\Lambda }}}_{n}$, while g_nj and h_nj are the wavefunctions of ${\eta }_{n}$ and ${\eta }_{n}^{\dagger }$.

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

Zoom In Zoom Out Reset image size

We find that for $\alpha \gtrsim 1$, ${c}_{{\rm{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 α. 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 \gt 2$ it fulfills ${\gamma }_{x}=\alpha$, consistent with the results of [12, 26]. However, (ii) for $1\lt \alpha \lt 2$ we obtain a hybrid exponential and algebraic decay with a different γ that depends linearly on α 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 α-dependence of ${\gamma }_{x}$ that is linear with slope $\sim 0.25(2)$. The fitted exponent ${\gamma }_{x}$ is shown in figure 2(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 α from zero to a finite value that appears to be θ-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.

Zoom In Zoom Out Reset image size

The energy spectrum changes in this region PM2 compared to the case PM1: the energy gaps ${\rm{\Delta }}E(1)$ and ${\rm{\Delta }}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}^{-}\rangle$ and $| {2}^{-}\rangle$, defined in section 2.2.3, become localized at the edges of the chain (figure 7(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the antiferromagnetic phase, denoted as AM, the effective central charge is instead zero for all α [26]. For $\alpha \gt 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 α. 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 ⁵ . 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.

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}_{{\rm{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 $\mathrm{cos}\theta \to -\mathrm{cos}\theta$ is lost for any finite α and the phase diagram is not symmetric around $\theta =\pi /2$.

Figure 1(b) shows that for $\alpha \gtrsim 1$ and $0\lt \theta \lt \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 central charge ${c}_{{\rm{eff}}}$ is zero within the phases for all $\alpha \gtrsim 1$ (see figure 1(d) for an example), while ${c}_{{\rm{eff}}}\ne 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 ${\rm{\Sigma }}(i,j)$ defined in section 2.2.2. In the region denoted as AM1, ${\rm{\Sigma }}(i,j)$ has a finite value for $| i-j| \to \infty$, while ${\rm{\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 ${\rm{\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 θ in figure 2(b), where we find numerically that the exponent ${\gamma }_{{\rm{\Sigma }}}$ for the power-law tail of ${\rm{\Sigma }}(i,j)$ equals ${\gamma }_{z}=\alpha$ when $\alpha \gt 1$. For $\alpha \lt 1$, however, the exponential part becomes numerically unobservable, and ${\rm{\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

$$\begin{eqnarray}\langle {a}_{R}^{\dagger }\;{a}_{0}\rangle ={{ \mathcal A }}_{\alpha ,\theta }\cdot \displaystyle \frac{{(-1)}^{R}{{\rm{e}}}^{-\xi R}}{\sqrt{R}}+{{ \mathcal B }}_{\alpha ,\theta }\cdot \{\begin{array}{cc}\displaystyle \frac{1}{{R}^{\alpha +1}} & \alpha \gt 2,\\ \displaystyle \frac{1}{{R}^{2\alpha -1}} & 1\lt \alpha \lt 2,\\ \displaystyle \frac{1}{{R}^{2-\alpha }} & 0\lt \alpha \lt 1,\end{array}\end{eqnarray} \tag{ 19 }$$

where the pre-factors ${{ \mathcal A }}_{\alpha ,\theta }$ and ${{ \mathcal 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 \gt 1$ and $\alpha \lt 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 ${{\rm{\Lambda }}}_{0}$ defined in section 2.2.3 as a function of α for a few values of the parameter 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.

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

$$\begin{eqnarray}{H}_{\mathrm{LRK}}=\mathrm{sin}\theta \displaystyle \sum _{k}(\begin{array}{cc}{a}_{k}^{\dagger } & {a}_{-k}\end{array})(\begin{array}{cc}{g}_{\alpha }(k)+\mathrm{cot}\theta & -{\rm{i}}(1+\epsilon ){f}_{\alpha }(k)\\ {\rm{i}}(1+\epsilon ){f}_{\alpha }(k) & -({g}_{\alpha }(k)+\mathrm{cot}\theta )\end{array})(\begin{array}{c}{a}_{k}\\ {a}_{-k}^{\dagger }\end{array}),\end{eqnarray} \tag{ 20 }$$

where a_k reads ${a}_{k}=\displaystyle \frac{1}{\sqrt{L}}{\displaystyle \sum }_{j}{{\rm{e}}}^{{\rm{i}}{kj}}{{\rm{a}}}_{j},$ $k=2\pi n/L$ is the lattice momentum with $n=0,...,L-1$, and the functions ${f}_{\alpha }(k)$ and ${g}_{\alpha }(k)$ read ${f}_{\alpha }(k)={\displaystyle \sum }_{{\ell }}\mathrm{sin}(k{\ell })/{{\ell }}^{\alpha }$ and ${g}_{\alpha }(k)={\displaystyle \sum }_{{\ell }}\mathrm{cos}(k{\ell })/{{\ell }}^{\alpha }$, respectively.

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

$$\begin{eqnarray}&&{H}_{\mathrm{LRK}}=\displaystyle \sum _{k}{\lambda }_{\alpha }(k)({\eta }_{k}^{\dagger }{\eta }_{k}-1/2),\end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray}&&{\lambda }_{\alpha }(k)=2\mathrm{sin}\theta \sqrt{{({g}_{\alpha }(k)+\mathrm{cot}\theta )}^{2}+{(1+\epsilon )}^{2}{f}_{\alpha }{(k)}^{2}}.\end{eqnarray} \tag{ 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{eqnarray}&&(\begin{array}{c}{\eta }_{k}^{\dagger }\\ {\eta }_{-k}\end{array})=(\begin{array}{cc}\mathrm{sin}{\beta }_{k} & \mathrm{icos}{\beta }_{k}\\ \mathrm{icos}{\beta }_{k} & \mathrm{sin}{\beta }_{k}\end{array})(\begin{array}{c}{a}_{k}^{\dagger }\\ {a}_{-k}\end{array}),\end{eqnarray} \tag{ 23 }$$

where $\mathrm{tan}(2{\beta }_{k})=(1+\epsilon )f(k)/[g(k)+\mathrm{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)=-{\displaystyle \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 $\mathrm{cot}\theta +{g}_{\alpha }(0)=0$ and on the line $\mathrm{cot}\theta +{g}_{\alpha }(\pi )=0;$ (ii) for a system in the thermodynamic limit $L\to \infty$ instead ${f}_{\alpha }(k)={\mathfrak{I}}{\mathrm{Li}}_{\alpha }({{\rm{e}}}^{{\rm{i}}k})$ and ${g}_{\alpha }(k)={\mathfrak{R}}{\mathrm{Li}}_{\alpha }({{\rm{e}}}^{{\rm{i}}k})$, ${\mathrm{Li}}_{\alpha }(x)$ being the polylogarithm of order α [69], and the critical line $\mathrm{cot}\theta +{g}_{\alpha }(0)=0$ ends at the point $\theta =\pi$ and $\alpha =1$.

For the critical line with $\theta \lt \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}_{\mathrm{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}_{\mathrm{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}_{\mathrm{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]

$$\begin{eqnarray}&&{e}_{0}(\alpha ,L)={e}_{\infty }(\alpha )-\displaystyle \frac{\pi \;{v}_{F}\;c}{6{L}^{2}},\end{eqnarray} \tag{ 24 }$$

where ${e}_{\infty }(\alpha )$ is the energy density in the thermodynamic limit, ${v}_{{\rm{F}}}=| 2\mathrm{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 \lt \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}_{\mathrm{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 ${\lambda }_{\alpha }(k)$ for the long range Kitaev models on the critical line with $\theta \gt \pi /2$, for $\alpha \gt 2$ we find ${\lambda }_{\alpha }(k)\sim k$, for $k\to 0$. However, for $1\lt \alpha \lt 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 α 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 \gt \pi /2$ and $\alpha \lt 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 \lt \pi /2$ conformal symmetry breaking arises even if the low-energy spectrum is linear for every α.

We locate the critical line (for $\theta \lt \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}_{{\rm{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}_{\mathrm{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}_{{\rm{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 α 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.