# A Regularized $(XP)^{2}$ Model

Zhao-Feng<EMAIL_ADDRESS>Yu-Qi<EMAIL_ADDRESS>CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China

###### Abstract

We investigate a dynamic model described by the classical Hamiltonian
$H(x,p)=(x^{2}+a^{2})(p^{2}+a^{2})$ , where $a^{2}>0$, in classical, semi-
classical, and quantum mechanics. In the high-energy $E$ limit, the phase path
resembles that of the $(XP)^{2}$ model. However, the non-zero value of $a$
acts as a regulator, removing the singularities that appear in the region
where $x,p\sim 0$ , resulting in a discrete spectrum characterized by a
logarithmic increase in state density. Classical solutions are described by
elliptic functions, with the period being determined by elliptic integrals. We
present three different forms of quantized Hamiltonians, and reformulate them
into the standard Schrödinger equation with $\cosh 2x$-like potentials.
Numerical evaluations of the spectra for these forms are carried out and
reveal minor differences in energy levels. Among them, one interesting form
possesses Hamiltonian in the Schrödinger equation that is identical to its
classical version. In such scenarios, the eigenvalue equations can be
expressed as the vanishing of the Mathieu functions’ value at $i\infty$
points, and furthermore, the Mathieu functions can be represented as the wave
functions.

##  I Introduction

The Hilbert-Polya conjecture posits that the non-trivial zeros of the Riemann
zeta function align with the eigenvalues of a self-adjoint operator, thereby
linking a mathematical conundrum to a physics challenge. While this conjecture
presents a promising avenue to validate the Riemann Hypothesis (RH) riemann:1
; riemann:2 ; riemann:3 ; riemann:4 ; riemann:5 ; riemann:6 ; riemann:7 ;
riemann:8 ; riemann:9 ; riemann:10 ; riemann:11 , identifying such a self-
adjoint operator whose spectrum mirrors the imaginary part of the zeros of the
Riemann zeta function, proves challenging.

Significant progress was made by Berry, Keating, and Connesberry:1999 ;
connes:1999  in 1999 when they proposed the $XP$ model. The primary model is
singular as $x,p\sim 0$ and predicts a continuous spectrum. Imposing
constraints on the $XP$ model in the phase spaces, Berry, Keating and Connes
derived discrete spectra ; notably, the state density scales with $\log E$ at
high energy limits in semi-classical estimations. In quantum physics, this
logarithmic association holds particular intrigue, potentially hinting at
scaling violations. The allure of the regularized $XP$ models lies in their
spectral behavior, which resonates with that of the Riemann zeros.

While a full quantum description of the regularized Berry, Keating and
Connesberry:1999 ; connes:1999  $XP$ model is absent, some several modified
$XP$ models with full quantum description are proposed. In 2011, Berry,
Keating and Sierra showed that the quantization of classical Hamiltonian
$H=x(p+1/p)$, and $H=(x+1/x)(p+1/p)$ contain the smooth approximation of the
Riemann zerossierra:1 ; berry:2011 , and later Sierra generalized the
Hamiltonians in terms of the family of Hamiltonians
$H=U(x)p+V(x)/p$sierra:2012 , capturing the dynamics of a massive particle in
a relativistic spacetime, with metrics informed by functions U and V.
Subsequent works inspired by the $XP$ model have delved deeper into models
associated with the Riemann zeros bender:2017 . In all above model, the
dynamic variable $p$ appears as nonlocal functions of $p$, which are
accommodated in conventional quantum mechanics..

In this paper, we propose a $(XP)^{2}$-like model, with classical Hamiltonian
defined as $H=(x^{2}+a^{2})(p^{2}+a^{2})$. This model offers an alternative
regularization to the traditional $XP$ model, where parameter $a$ acts as a
regulator removing the singularities for $x,p$ near zero, thus, a discrete
spectrum is expected after imposing quantization. Semiclassical estimation
shows that the state density is proportional to $\log E$ in the high energy
limit. The Hamiltonian in our model is local functions of $x$ and $p$, thus
quantization of this model is straightforward, although it suffers from the
operator ordering ambiguities. We explore multiple versions of the quantized
Hamiltonian. Our numerical analyses depict minor spectral differences in high
energy limit. Through variable transformations, we elucidate that both the
classical Hamiltonian and its quantum counterparts can be framed within the
standard Schrödinger equation with $\cosh^{2}$-like potentials. We furnish
numerical comparisons across various versions and interpret their physical
implications.

This paper unfolds as follows: Section II elaborates on the classical
solutions. Section III delves into the semi-classical implications on spectra.
Section IV explore multiple versions of the quantized Hamiltonian and
furnishes numerical calculations on spectra and wave functions. Section V
bridges our findings with conventional Schrödinger equations. Lastly, Section
V contributes conclusions.

##  II  Classical solutions

We start by introducing a classical Hamiltonian defined as:

$H=(x^{2}+a^{2})(p^{2}+a^{2})\,$ (1)

where $a$ is real number. For $a=0$, this Hamiltonian reduces to the square of
the $XP$ model, and consequently, its dynamics closely resemble those of the
conventional $XP$ model. The non-zero $a$ acts as a regulator, constraining
the system’s path within a finite region in the phase space.

For a more general form

$H=(x^{2}+a^{2})(p^{2}+b^{2})\,$ (2)

with $a\neq b$, one may impose a canonical transformation to rescale $x$, $p$
as $x^{\prime}=rx$, $p^{\prime}=p/r$ respectively, using a real number
$r=\sqrt{|a|/|b|}$ . This transformation allows us to rewrite the Hamiltonian
in the form of equation Eq.(1).

$H=(x^{\prime 2}+|ab|)(p^{\prime 2}+|ab|)\,,$ (3)

which aligns with the form of Eq.(1).

The classical canonical Hamilton equations of motion are given by

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial p}=2p(x^{2}+a^{2})\;\;\;,$ (4) $\displaystyle\dot{p}$
$\displaystyle=$ $\displaystyle-\frac{\partial H}{\partial
x}=-2x(p^{2}+a^{2})\;\;\;.$ (5)

This system conserves energy, with paths in phase space adhering to a constant
energy surface given by $H_{E}=(x^{2}+a^{2})(p^{2}+a^{2})$. For a specified
energy $H_{E}$, $x$ and $p$ move periodically with respect to time $t$, as
depicted Fig.1. Setting $H_{E}=E^{2}+a^{4}$ , these equations simplify to:

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\sqrt{(x^{2}+a^{2})(E^{2}-a^{2}x^{2})}\;\;\;,$ (6)
$\displaystyle\dot{p}$ $\displaystyle=$
$\displaystyle-\sqrt{(p^{2}+a^{2})(E^{2}-a^{2}p^{2})}\;\;\;,$ (7)

in which $x(t)$, $p(p)$ are decoupled and then can be expressed by the
elliptical functions ${\rm sn(u|m)}$ as

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle\frac{E}{a}\,{\rm
sn}\left(2\,a^{2}\,t\Big{|}-\frac{E^{2}}{a^{4}}\right)\,,$ (8) $\displaystyle
p(t)$ $\displaystyle=$ $\displaystyle\frac{E}{a}\,{\rm
sn}\left(2\,a^{2}\,\left(t+\frac{T_{E}}{4}\right)\Big{|}-\frac{E^{2}}{a^{4}}\right)\,,$
(9)

where $T_{E}$ denotes the period time, characterized by the complete elliptic
integral of the first kind, $\mathcal{K}(x)$:

$T_{E}=\frac{2}{a^{2}}{\displaystyle\int_{0}^{\frac{\pi}{2}}d\theta\,\frac{1}{\sqrt{1+\frac{E^{2}}{a^{4}}\sin^{2}\theta}}}\equiv\frac{2}{a^{2}}\,\mathcal{K}\left(-\frac{E^{2}}{a^{4}}\right)\;.$
(10)

Here we have set initial condition $x(t=0)=0$. The solution of $x(t)$ and
$p(t)$ retain the same form as expected but $\pi/4$ phase difference.In Fig.2,
the allowed phase space is shown, bounded by the classical path.

Figure 1: Classical trajectory of $H=(x^{2}+a^{2})(p^{2}+a^{2})$, when
$x_{0}=0,E=3,x_{0}=0,p_{0}=\sqrt{E^{2}-1}$.

In the high energy limit, $E\to\infty$, the leading term of $T_{E}$ goes to

$\lim_{E\to\infty}T_{E}=\lim_{E\to\infty}\frac{2}{a^{2}}\,{\cal
K}\left(-\frac{E^{2}}{a^{4}}\right)=\frac{2}{E}\ln\left(\frac{4E}{a^{2}}\right)\;\,.$
(11)

Thus in this limit, $T_{E}$ exhibits a proportionality to the inverse power of
$E$ and the logarithmic factor $\ln\left(\frac{4E}{a^{2}}\right)$, which can
be interpreted as scaling violation, and the parameter $a$ acts as a
regulator, ensuring $T_{E}$ remains finite. This result can be interpreted
through Eq. (4) . Actually, in the high energy limit,

$\dot{x}=2P(x^{2}+a^{2})\simeq 2Px^{2}\simeq 2E\,x\;,$ (12)

from which we infer that :

$T_{E}\simeq\frac{2}{E}\int_{a}^{E/a}\frac{dx}{x}\simeq\frac{2}{E}\ln\frac{E}{a^{2}}\;,$
(13)

which reproduce the logarithmic $E$ term in (11). In the next section, we will
see that this logarithmic factor relates to the state density logarithmic
growth with the energy increasing it high energy limit.

Figure 2: The region between the two curves describes the allowed phase space
with area A bounded by the classical path.

Finally, we explore a classical canonical transformation from $x,p\to u,v$,
which preserves the classical Poisson bracket, represented as:

$\displaystyle x$ $\displaystyle=$ $\displaystyle a\,\sinh u\,,$ (14)
$\displaystyle p$ $\displaystyle=$ $\displaystyle\frac{v}{a\,\cosh u}\,\,.$
(15)

Following the transformation, the Hamiltonian can be reframed as:

$H(u,v)=v^{2}+a^{4}\cosh^{2}u=v^{2}+a^{4}\sinh^{2}u+a^{4}\;.$ (16)

This Hamiltonian corresponds to the standard form for a particle with half-
unit mass, moving one-dimensionally in the potential
$a^{4}\cosh^{2}u=a^{4}\sinh^{2}u+a^{4}$.

The findings indicate that, at the classical mechanics level, the Hamiltonians
presented in Eq.(1) and Eq.(16) depict identical dynamics.

##  III  Semi-classical estimation for the spectrum

The semi-classical estimation serves as a bridge between classical and quantum
descriptions, offering insights into fundamental features of quantum
representation. We employ the semi-classical approximation to estimate the
spectrum of the Hamiltonians presented in (1) and in (16) and to delve into
their high-energy asymptotic behavior.

Considering the classical action variable:

$I=\frac{1}{2\pi}\oint pdq,$ (17)

which is the counter area of the phase space trajectory with given energy $E$,
over $2\pi$. This action variable remains invariant under classical canonical
transformations and is a function of $E$.

The one-dimensional periodic motion suggests a discrete energy spectrum. Thus,
the semi-classical approximation can be employed to estimate $N(E)$, the
number of the quantum states below the energy $E$:

$N(E)\equiv{I\over\hbar}\;.$ (18)

Consequently, for the Hamiltonians defined in (1) and in (16), we deduce:

$N(E)={4\over\hbar}\int_{0}^{\frac{E}{a}}\sqrt{\frac{E^{2}-a^{2}x^{2}}{x^{2}+a^{2}}}dx\;.$
(19)

This integral can be represented using the incomplete elliptic integral of the
second kind $\mathcal{E}\left(\phi|m\right)$:

$N(E)=-{4i\over\hbar}\enspace
E\mathcal{E}\left(i\sinh^{-1}\frac{E}{a^{2}}\Big{|}\frac{a^{4}}{E^{2}}\right)\;.$
(20)

It can also be defined using the complete elliptic integrals of the first kind
$\mathcal{K}(u)$ and the second kind $\mathcal{E}(u)$:

$N(E)={4\over\hbar}\left[\frac{E^{2}+a^{4}}{a^{2}}\mathcal{K}\left(-E^{2}/a^{4}\right)-a^{2}\mathcal{E}\left(-E^{2}/a^{4}\right)\right]\;.$
(21)

In the high-energy limit, the energy gap between two adjacent bound states can
be approximated as:

$\Delta E=\hbar\frac{dE}{dI}\Delta
n=\frac{\hbar}{2E}\frac{dH}{dI}=\frac{\omega_{E}\hbar}{2E}=\frac{\pi\hbar}{T_{E}E}\;\,,$
(22)

where ${\omega_{E}}$ denotes the circular frequency for periodic motion at
energy $E$. Referring to Eq. (11), in high energy limit,

$\Delta E=\frac{\pi\hbar}{\displaystyle 2\ln\frac{4E}{a^{2}}}\;\,$ (23)

Therefore, from this semi-classic estimation we see that in the high-energy
limit, the energy gap between two adjacent bound states is proportional to the
inverse power of $E$.

##  IV Quantum descriptions

The quantum mechanics description for the models associated with the classical
Hamiltonian presented in Eq.(1) is introduced next. Through the process of
canonical quantization, the classical variables $x$ and $p$ in the Hamiltonian
are substituted with their operator counterparts $\hat{x}$ and $\hat{p}$,
which obey the Poisson bracket $[\hat{x},\hat{p}]=i\hbar$. To simplify the
expressions, we often set $\hbar=1$. This resetting can be realized through
the rescaling: $\hat{x}\to\sqrt{\hbar}\hat{x}$,
$\hat{p}\to\sqrt{\hbar}\hat{p}$, $a\to\hbar a$, $E\to\hbar E$,
$H\to\hbar^{2}H$. The ensuing Poisson bracket is $[\hat{x},\hat{p}]=i$. Notice
that classic Hamiltonian given in (1) contains multiplication of $x^{2}$ and
$p^{2}$. Its quantization suffers from operator ordering ambiguities since
$\hat{x}$ and $\hat{p}$ do not commutate within quantum mechanics. This
implies the potential for various quantum Hamiltonian forms corresponding to
the classical Hamiltonian in Eq.(1). We enumerate several representative self-
conjugate Hamiltonian forms below:

Form I preserves the symmetry of $\hat{x},\hat{p}$ exchange:

$H_{1}(\hat{x},\hat{p})=\frac{1}{2}\left[\left(\hat{x}^{2}+a^{2}\right)\left(\hat{p}^{2}+a^{2}\right)+\left(\hat{p}^{2}+a^{2}\right)\left(\hat{x}^{2}+a^{2}\right)\right]\;.$
(24)

This form can also be expressed as:

$H^{\prime}_{1}(\hat{x},\hat{p})=\frac{1}{2}\left(A{A}^{\dagger}+A^{\dagger}A\right)\,,$
(25)

with

$A\equiv\frac{1}{2}\left(\hat{x}\hat{p}+\hat{p}\hat{x}\right)+ia\left(\hat{x}+\hat{p}\right)-a^{2}\;,$
(26)

and

$H^{\prime\prime}_{5}(\hat{x},\hat{p})=\frac{1}{2}\left(B{B}^{\dagger}+B^{\dagger}B\right)\,,$
(27)

with

$B\equiv\frac{1}{2}\left(\hat{x}\hat{p}+\hat{p}\hat{x}\right)+ia\left(\hat{x}-\hat{p}\right)-a^{2}\;.$
(28)

The Hamiltonians $H^{\prime}_{1}(\hat{x},\hat{p})$ and
$H^{\prime\prime}_{1}(\hat{x},\hat{p})$ only differ from
$H_{1}(\hat{x},\hat{p})$ by a constant value:

$H^{\prime}_{1}(\hat{x},\hat{p})=H^{\prime\prime}_{1}(\hat{x},\hat{p})=H_{1}(\hat{x},\hat{p})+\frac{3}{4}\;.$
(29)

This implies that $H_{1}(\hat{x},\hat{p})$ $H_{4}(\hat{x},\hat{p})$
$H_{5}(\hat{x},\hat{p})$ describe identical quantum systems. Therefore, for
Form I, the focus is primarily on $H_{1}(\hat{x},\hat{p})$.

Form II, an alternative quantum version, is presented as:

$H_{2}(\hat{x},\hat{p})=\sqrt[4]{\hat{x}^{2}+a^{2}}\,\hat{p}\,\sqrt{\hat{x}^{2}+a^{2}}\,\hat{p}\,\sqrt[4]{\hat{x}^{2}+a^{2}}+a^{2}\left(\hat{x}^{2}+a^{2}\right)\;.$
(30)

Contrastingly, this form breaks $X,P$ exchange symmetry, yet it offers a more
streamlined Schrödinger equation upon a variable transformation.

Yet another quantum Hamiltonian, denoted as Form III, is described as:

$H_{3}(\hat{x},\hat{p})=\sqrt{\hat{x}^{2}+a^{2}}\left(\hat{p}^{2}+a^{2}\right)\sqrt{\hat{x}^{2}+a^{2}}\;.$
(31)

This form also breaks the $X,P$ exchange symmetry.

Given three kinds of quantum Hamiltonians above, one can derive differential
equations for the wave-functions in the configuration space using the
substitution $\hat{x}\to x$, $\hat{p}\to{\displaystyle-i\frac{d}{dx}}$.

$H_{I}(x,-i\frac{d}{dx})\;\psi_{I}(x)=(E^{2}+a^{4})\;\psi_{I}(x)\;\;\;\;(I=1,2,3)\;.$
(32)

We discuss potential solutions to these equations. Initially, for the
$x\to\infty$ limit, the equations are predominantly influenced by terms like
$x^{2}(-{\displaystyle\frac{d^{2}}{dx^{2}}}+a^{2})$, leading to the asymptotic
behavior:

$\psi_{I}(x)\;\stackrel{{\scriptstyle x\to\infty}}{{\longrightarrow}}\;e^{\pm
ax}\;\;\;\;\;(I=1,2,3)\;,$ (33)

For arbitrary real values of $E$. However, only when $H_{E}=E^{2}+a^{4}$ acts
as an eigenvalue and the corresponding solution $\psi_{I}(x)$ serves as the
eigenstate in Hilbert space, must it adhere to:

$\displaystyle\psi_{I}(x)\;\stackrel{{\scriptstyle
x\to+\infty}}{{\longrightarrow}}\;e^{-ax}\;,$ (34)
$\displaystyle\psi_{I}(x)\;\stackrel{{\scriptstyle
x\to-\infty}}{{\longrightarrow}}\;e^{+ax}\;.$ (35)

It’s evident that the Hamiltonians $H_{I}(x)\;(I=1,2,3)$ remain invariant
under parity transformations, i.e.
$H_{I}(\hat{x},\hat{p})=H_{I}(-\hat{x},-\hat{p})$. This suggests that physical
wave functions possess definitive parities:

It implies that the physical wave functions have definite parities:

$\psi_{I}(x)=\pm\;\psi_{I}(-x)\;\;\;\;\;(I=1,2,3)\;.$ (36)

With respect to all the aforementioned quantum Hamiltonians, Eq. (32) can be
restated as a secondary ordinary differential equation. Specific
representations are:

For form I quantum Hamiltonian (24), the differential equation can be written
as

$(x^{2}+a^{2})\;\psi_{1}^{\prime\prime}(x)+2x\;\psi_{1}^{\prime}(x)+(E^{2}+1-a^{2}\,x^{2})\;\psi_{1}(x)=0\;.$
(37)

For form II quantum Hamiltonian (30), the differential equation can be
expressed as

$\sqrt[4]{x^{2}+a^{2}}\;\frac{d}{dx}\left[\sqrt{x^{2}+a^{2}}\frac{d}{dx}\left(\sqrt[4]{x^{2}+a^{2}}\,\psi_{2}(x)\right)\right]+(E^{2}-a^{2}\,x^{2})\;\psi_{2}(x)=0\;.$
(38)

Setting $\phi_{2}(x)=\left(\sqrt[4]{x^{2}+a^{2}}\,\psi_{2}(x)\right)$,
$\phi_{2}(x)$ then satisfies the following secondary ordinary differential
equations

$(x^{2}+a^{2})\;\phi_{2}^{\prime\prime}(x)+x\;\phi_{2}^{\prime}(x)+(E^{2}-a^{2}\,x^{2})\;\phi_{2}(x)=0\;.$
(39)

For form III quantum Hamiltonian (31), the differential equation can be reads

$\sqrt{x^{2}+a^{2}}\;\frac{d^{2}}{dx^{2}}\left(\sqrt{x^{2}+a^{2}}\,\psi_{3}(x)\right)+(E^{2}+a^{4})\;\psi_{3}(x)=0\;.$
(40)

Setting $\phi_{3}(x)=\left(\sqrt{x^{2}+a^{2}}\,\psi_{3}(x)\right)$, it then
satisfies the following secondary ordinary differential equations

$(x^{2}+a^{2})\;\phi_{3}^{\prime\prime}(x)+(E^{2}-a^{2}\,x^{2})\;\phi_{3}(x)=0\;.$
(41)

Upon comparing Eqs.(37, 39, 41) corresponding to the distinct Hamiltonians
(24, 30, 31), it’s evident they can be consolidated into a general form:

$(x^{2}+a^{2})\;\phi_{I}^{\prime\prime}(x)+(3-I)x\;\phi_{I}^{\prime}(x)+(E^{2}+c_{I}-a^{2}\,x^{2})\;\phi_{I}(x)=0\;\;\;\;\;\;(I=1,2,3)\;,$
(42)

assuming $\phi_{1}(x)=\psi_{1}(x)$ and $c_{1}=1,c_{2}=c_{3}=0$.

Employing the boundary conditions from Eq.(34), one can numerically solve
these equations to obtain eigenvalue $H_{E}=E^{2}+a^{4}$ and the associated
wave functions.

The parity property in Eq. (36) can streamline calculations. For instance,
parity-even solutions yield

$\phi^{\prime}(0)=0\;,$

while parity-odd solutions produce

$\phi(0)=0\;.$

Numerical analyses reveal that exact solutions closely align with semi-
classical energy level estimations for $E$, especially as $E\to\infty$. We
showcase the first 10 energy and the 40th-50th eigenvalues for $a=1$ in
table.1, and in table.2. And discrepancies between the first 20 energy
eigenvalues of H1, H2, and H3 concerning the semi-classical energy eigenvalues
are shown in Fig.3.

Figure 3: Discrepancies between the energy eigenvalues of H1, H2, and H3
concerning the semi-classical energy eigenvalues. The black points correspond
to the 1st to 20th energy eigenvalues of H1, the red points represent the same
for H2, and the blue points respresent that for H3. Table 1: The first ten
Energy eigenvalues of H1, H2, H3 and the semi-classical estimation

Hamiltonian N | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
H1 | 0.80879 | 1.82879 | 2.58598 | 3.24636 | 3.85276 | 4.42303 | 4.96681 | 5.49002 | 5.99659 | 6.48931
H2 | 1.10064 | 1.99182 | 2.70651 | 3.34474 | 3.93687 | 4.49703 | 5.0332 | 5.55043 | 6.05216 | 6.54086
H3 | 0.93685 | 1.95745 | 2.69087 | 3.33685 | 3.93281 | 4.49517 | 5.03272 | 5.55085 | 6.05319 | 6.54232
Semi-classical | 1.05493 | 1.97222 | 2.69374 | 3.33520 | 3.92923 | 4.49064 | 5.02770 | 5.54559 | 6.04784 | 6.53696

Table 2: The 40th-50th ten Energy eigenvalues of H1, H2, H3 and the semi-
classical estimation

Hamiltonian N | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50
---|---|---|---|---|---|---|---|---|---|---
H1 | 18.68292 | 19.04667 | 19.40881 | 19.76939 | 20.12845 | 20.48605 | 20.84223 | 21.19702 | 21.55047 | 21.90261
H2 | 18.70143 | 19.06483 | 19.42664 | 19.78689 | 20.14565 | 20.50296 | 20.85885 | 21.21337 | 21.56656 | 21.91845
H3 | 18.73068 | 19.03024 | 19.41945 | 19.82975 | 20.13614 | 20.47012 | 20.87308 | 21.25504 | 21.53384 | 21.90041
Semi-classical | 18.70039 | 19.06382 | 19.42565 | 19.78593 | 20.14471 | 20.50204 | 20.85795 | 21.21249 | 21.56569 | 21.91759

##  V Schrödinger forms

In this section, we demonstrate that all quantum Hamiltonians specified in
Eqs. (24), (30), (31) can be recast into a standard Schrödinger form. In this
form, the Hamiltonian is expressed as a sum of kinetic term and a potential
term. While the kinetic terms consistently remain the same, the potential
components exhibit variations.

We first impose a variable substitution, Eq.(14), $x=a\sinh u$ on the the
differential equations (37) (39) (41). The derivative term then expressed as:

$\frac{d\phi}{dx}=\frac{du}{dx}\;\frac{d\phi}{du}=\frac{1}{a\cosh
u}\;\frac{d\phi}{du}\;.$ (43)

With (14) and (43), Eq.(37) can be rewritten as:

$\psi_{1}^{\prime\prime}(u)+\tanh
u\,\psi_{1}^{\prime}(u)+(E^{2}+1-a^{4}\sinh^{2}u)\psi_{1}(u)=0\;.$ (44)

Giving $\varphi_{1}(u)=\sqrt{\cosh(u)}\;\psi_{1}(u)$, it reduces to the
following differential equation:

$\varphi_{1}^{\prime\prime}(u)+\left(\frac{1}{2}+E^{2}+\frac{1}{4}\tanh^{2}u-a^{4}\sinh^{2}u\right)\varphi_{1}(u)=0\;,$
(45)

in which while the second derivative term is persists, the first derivative
term is eliminated. Consequently, $\varphi_{1}(u)$ adheres to the standard
Schrödinger equation with a half unit mass and a characteristic effective
potentia

$V_{1}(u)=-\frac{1}{2}-\frac{1}{4}\tanh^{2}u+a^{4}\sinh^{2}u\;.$ (46)

We then derive the Schrödinger form of the Hamiltonia form (30). By applying
Eqs.(14,43), Eq.(39) can be rewritten as:

$\psi_{2}^{\prime\prime}(u)+(E^{2}-a^{4}\sinh^{2}u)\psi_{2}(u)=0\;.$ (47)

Therefore, $\psi_{2}(u)$ complies with the standard Schrödinger equation, once
again sporting a half unit mass but associated with an alternative effective
potential:

$V_{2}(u)=a^{4}\sinh^{2}(u)\;.$ (48)

We emphasis here that this Schrödinger Hamiltonian remains exactly the same
form as the classic one Eq.(16).

Lastly, addressing the Hamiltonian form from Eq.(31) and applying Eqs.(14,43),
Eq.(41) can be cast into:

$\phi_{3}^{\prime\prime}(u)-\tanh
u\,\phi_{3}^{\prime}(u)+(E^{2}-a^{4}\sinh^{2}u)\phi_{3}(u)=0\;.$ (49)

Introducing $\varphi_{3}(u)=\phi_{3}(u)/\sqrt{\cosh u}$, it conforms to the
differential equation:

$\varphi_{3}^{\prime\prime}(u)+\left(\frac{1}{2}+E^{2}-\frac{3}{4}\tanh^{2}u-a^{4}\sinh^{2}u\right)\varphi_{3}(u)=0\;.$
(50)

Again, the second derivative term is retained while the first one vanishes.
This ensures that $\varphi_{3}(u)$ aligns with the standard Schrödinger
equation, furnished with half unit mass, and steered by its own unique
effective potential:

$V_{3}(u)=-\frac{1}{2}+\frac{3}{4}\tanh^{2}u+a^{4}\sinh^{2}u\;.$ (51)

Figure 4: The effective potentials $V_{1}(u)$,$V_{2}(u)$ and $V_{3}(u)$, when
$a=1$.

We show the three kind of effective potentials $V_{1}(u),V_{2}(u)$ and
$V_{3}(u)$in Fig.4 ,when $a=1$. Obversely, those Schrödinger forms differ only
by the potentials. For all value of, the differences among them are minor as
shown in Fig.4. Especially, when $u\to\pm\infty$, all potentials $V_{1}(u)$,
$V_{2}(u)$, $V_{3}(u)$ given in Eqs. (46)(48)(51) exponentially increase with
$u$: $\exp^{2|u|}$, with constant differences.

Such Schrödinger equations can be tackled as eigenvalue problems. For an
arbitrary real number $E$, these solutions are divergent when $u\to\pm\infty$.
Actually, as discussed in the last section in Eq.(33), when $u\to\pm\infty$,
following either solutions are possible:

$\varphi_{I}(u)\;\stackrel{{\scriptstyle
u\to\infty}}{{\longrightarrow}}\;e^{\pm\cosh u}\;\;\;\;(I=1,2,3)\;.$ (52)

When and only when $H_{E}=E^{2}+a^{4}$ is the eigenvalue, the following
convengent solution survives

$\displaystyle\varphi_{I}(u)\;\stackrel{{\scriptstyle
u\to\pm\infty}}{{\longrightarrow}}\;e^{-\cosh
u}\;{\longrightarrow}\;0\;\;\;(I=1,2,3)\;,$ (53)

Then the solution as the wave-function is normalizable in Hilbert space.

Given the parity invariance of these Hamiltonians, eigenstates are distinctly
marked by definitive parities, either parity-odd or parity-even. Implementing
this condition alongside boundary condition:

$\displaystyle\varphi_{I}(\infty)=0\;\;\;(I=1,2,3)\;,$ (54)

one can numerically resolve the Schrödinger equations, procuring both
eigenvalues and wave-functions.

For Eq.(47), solutions can be gained by Mathieu functions. Original Mathieu
functions serve as solutions for the differential equation:

$\frac{d^{2}\,y}{d\,w^{2}}+(c-2q\cos(2w))y=0\;,$ (55)

with parameters $c,q$, and variable $w$ all being real numbers. It is nothing
but the Schrödinger equation with unit mass and periodic potential $2\cos 2w$.
Two independent solutions exist: $y_{1}(w)={\cal M}c(c,q,w)$ and
$y_{2}(w)={\cal M}s(c,q,w)$, corresponding to parity-even and parity-odd
solutions respectively.

Notice that Eq. (50) are related to Eq. (55), by an analytic continuation via
setting variable $w=iu$ and parameters $a^{4}=q$, $E^{2}=c$, resulting in
their solutions are connected by the same analytic continuation.

Specifically, the parity odd and parity even solutions for Eq. (47) are given
by

$\displaystyle\varphi_{2}^{-}(u)={\mathcal{M}}s(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},iu)\;,$
(56)
$\displaystyle\varphi_{2}^{+}(u)={\mathcal{M}c}(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},iu)\;,$
(57)

respectively.

We then can apply boundary conditions Eq. (54) expressed in the following
forms

$\displaystyle\varphi_{2}^{-}(\infty)={\mathcal{M}}s(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},i\infty)=\varphi_{2}^{+}(\infty)={\mathcal{M}c}(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},i\infty)=0\;,$
(58)

to determine the eigenvalues $H_{E}=E^{2}+a^{4}$. Therefore they represent the
eigenvalue equations. The solutions
$\varphi_{2}^{-}(u)={\mathcal{M}}s(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},iu)$
and
$\varphi_{2}^{+}=(u)={\mathcal{M}}s(E^{2}-\frac{a^{4}}{2},-\frac{a^{4}}{4},iu)$
are the corresponding wavefunctions for parity odd and parity even,
respectively.

##  VI conclusions

In this paper, we have studied a dynamic model with the classical Hamiltonian
$H(x,p)=(x^{2}+a^{2})(p^{2}+a^{2})$ $where(a^{2}>0)$ within the frameworks of
classic, semi-classic, and quantum mechanism, respectively. At high energy $E$
limit, the path in phase space resembles that of the $(XP)^{2}$ model.

Historically, the $XP$ model has been heralded as an intriguing paradigm that
may potentially shed light on certain facets of the Riemann Hypothesis,
especially when contextualized within the purview of the Hilbert-Polya
conjecture. The primary $XP$ model yields a continuous spectrum. However, with
appropriate regulation, the model exhibits a discrete spectrum. This spectrum
is marked by a state density that increases logarithmically, analogous to the
distribution of the Riemann zeros in the semi-classical approximation. A
comprehensive quantum depiction of such a regulated $XP$ model is
conspicuously absent in the existing literature. In this paper, we introduce a
regulator for the $(XP)^{2}$ model. Our work sheds some lights on the
quantization of the regulated $XP$ model, and further work on this subject is
currently in preparation.

We have explored three versions of the quantization model and conducted
numerical calculations for each. In presenting their corresponding Schrödinger
forms, we found that quantum Hamiltonian form II (30) is particularly
intriguing for its Hamiltonian appearing in Schrödinger Eq. (44) is exactly
the same with the classic version Eq. (16). In such cases,the eigenvalue
equations can be expressed as conditions for the vanishing of the Mathieu
functions’ values at $i\infty$ points. Furthermore, these Mathieu functions
can be presented as the wavefunctions of the system.

## References

  * (1) B. Riemann, ”On the number of primes less than a given magnitude”, translated by H. M. Edwards from ”Ueber dieAnzahl der Primzahlen unter einer gegebenen Grosse” (1859).
  * (2) H.M. Edwards, Riemann’s Zeta Function, Academic Press, New York, 1974.
  * (3) E.C. Titchmarsh, The Theory of the Riemann Zeta Function, Oxford Univ. Press, Oxford, 1986.
  * (4) H. Davenport, Multiplicative Number Theory, Grad. Texts in Math., vol. 74, Springer-Verlag, New York, 2000.
  * (5) E. Bombieri, Problems of the Millennium: The Riemann Hypothesis, Clay Mathematics Institute, (2000).
  * (6) P. Sarnak, Problems of the Millennium: The Riemann Hypothesis, Clay Mathematics Institute, (2004).
  * (7) B. Conrey, ”The Riemann Hypothesis, Notices of AMS, March, 341 (2003).
  * (8) A. Ivisc, ”On some reasons for doubting the Riemann hypothesis”, arXiv:math/0311162.
  * (9) G. J. Chaitin, ”Thoughts on the Riemann hypothesis”, arXiv:math/0306042.
  * (10) P. Borwein, S. Choi, B. Rooney and A. Weirathmueller (Eds.). The Riemann Hypothesis. A Resource for the Ancionado and Virtuoso Alike, CMS Books in Mathematics, Springer, Canada, 2008.
  * (11) A. Connes, ”An essay on the Riemann Hypothesis”, arXiv:1509.05576.
  * (12) M. V. Berry, J. P. Keating, SIAM Review 41, 236, 1999.
  * (13) A. Connes, Selecta Mathematica New Series 5 29, (1999); math.NT/9811068.
  * (14) G. Sierra and J. Rodrspguez-Laguna, ”The H = xp model revisited and the Riemann zeros”. Phys. Rev. Lett. 106, 200201 (2011); arXiv:1102.5356.
  * (15) M. V. Berry and J. P. Keating, ”A compact hamiltonian with the same asymptotic mean spectral density as the Riemann zeros”, J. Phys. A: Math. Theor. 44, 285203 (2011).
  * (16) G. Sierra, ”General covariant xp models and the Riemann zeros”, J. Phys. A: Math. Theor. 45 055209 (2012); arXiv:1110.3203.
  * (17) Carl M. Bender, Dorje C. Brody, and Markus P. Mueller,PRL 118, 130201 (2017)


# Universal signatures of Majorana zero modes in critical Kitaev chains

Nicolas Laflorencie  Laboratoire de Physique Théorique, Université de
Toulouse, CNRS, UPS, France

###### Abstract

Many topological or critical aspects of the Kitaev chain are well known, with
several classic results. In contrast, the study of the critical behavior of
the strong Majorana zero modes (MZM) has been overlooked. Here we introduce
two topological markers which, surprisingly, exhibit non-trivial signatures
over the entire (1+1) Ising critical line. We first analytically compute the
MZM fidelity ${\cal{F}}_{\rm MZM}$–a measure of the MZM mapping between parity
sectors. It takes a universal value along the (1+1) Ising critical line,
${\cal{F}}_{\rm MZM}=\sqrt{8}/\pi$, independent of the energy. We also obtain
an exact analytical result for the critical MZM occupation number ${{\cal
N}}_{\rm MZM}$ which depends on the Catalan’s constant ${\cal G}\approx
0.91596559$, for both the ground-state (${{\cal N}}_{\rm
MZM}=1/2-4{\cal{G}}/\pi^{2}\approx 0.12877$) and the first excited state
(${{\cal N}}_{\rm MZM}=1/2+(8-4{\cal{G}})/\pi^{2}\approx 0.93934$). We further
compute finite-size corrections which identically vanish for the special ratio
$\Delta/t=\sqrt{2}-1$ between pairing and hopping in the critical Kitaev
chain.

##  I Introduction

Generalities— Much attention has recently turned towards emerging Majorana
bound states in certain condensed matter systems Lutchyn  _et al._ (2010);
Oreg  _et al._ (2010); Sarma  _et al._ (2015); Marra (2022); Mi (2022); Dvir
(2023); Yazdani  _et al._ (2023). One of the simplest toy-model hosting such a
fascinating physics is an exactly solvable quantum chain model, solved in the
magnetic language quite some time ago by Lieb, Schulz, Mattis Lieb  _et al._
(1961) and Pfeuty Pfeuty (1970). Nevertheless, a decisive step was later taken
thanks to the seminal work of Kitaev Kitaev (2001) who realized that non-
trivial topological properties could emerge in such a simple quantum model
when rephrased in fermionic language, nowadays referred to as the Kitaev chain

$\displaystyle{\cal{H}}_{\rm K}$ $\displaystyle=$
$\displaystyle-\sum_{j}\left([t_{j}^{\vphantom{\dagger}}c_{j}^{\dagger}c_{j+1}^{\vphantom{\dagger}}+\Delta_{j}^{\vphantom{\dagger}}c_{j}^{\dagger}c_{j+1}^{{\dagger}}+{\rm{h.c.}}]-\mu_{j}^{\vphantom{\dagger}}c_{j}^{\dagger}c_{j}^{\vphantom{\dagger}}\right)~{}~{}~{}~{}~{}$
(1)

which describes a non-interacting p-wave superconducting wire with hopping
$t_{j}$, pairing $\Delta_{j}$ and potential $\mu_{j}$. Equivalently it
represents the original Lieb  _et al._ (1961) spin chain XY Hamiltonian

${\cal{H}}_{\rm
K}=-\sum_{j}\left(X_{j}^{\vphantom{x}}\sigma_{j}^{x}\sigma_{j+1}^{x}+Y_{j}^{\vphantom{y}}\sigma_{j}^{y}\sigma_{j+1}^{y}+h_{j}^{\vphantom{z}}\sigma_{j}^{z}\right),$
(2)

with couplings $X_{j}=\frac{t_{j}+\Delta_{j}}{2}$,
$Y_{j}=\frac{t_{j}-\Delta_{j}}{2}$, and fields $h_{j}=\frac{\mu_{j}}{2}$.

The total number of fermions ${{N}}_{\rm
f}=\sum_{j}c_{j}^{\dagger}c_{j}^{\vphantom{\dagger}}$
($\sum_{j}(\sigma_{j}^{z}+1)/2$ in the spin language) is not fixed (unless
$\Delta_{j}=X_{j}-Y_{j}=0$), but its parity is conserved:
${{\mathbb{P}}}=\left(-1\right)^{{{N}}_{\rm f}}=\prod_{j}\sigma_{j}^{z}$ has
eigenvalues $p=\pm 1$ and commutes with ${\cal{H}}_{\rm K}$. This yields a
global $\mathbb{Z}_{2}$ symmetry to the problem and one can therefore group
the eigenstates in two distinct parity sectors. The spontaneous breaking of
$\mathbb{Z}_{2}$ symmetry may occur in the thermodynamic limit, associated
with magnetic order $\langle\sigma^{x}\rangle\neq 0$. This standard long-range
order takes a non-local form in fermionic language, with "topological"
unpaired zero-energy Majorana edge states localized at the boundaries Kitaev
(2001).

Strong Majorana zero modes— The idea of strong Majorana zero-mode (MZM),
popularized by Fendley Fendley (2012, 2016) and collaborators Kemp  _et al._
(2017); Else  _et al._ (2017); Olund  _et al._ (2023), goes beyond low energy
Alexandradinata  _et al._ (2016); Wada  _et al._ (2021) as the whole many-body
spectrum is involved, see inset (i) of Fig. 1. A strong MZM operator $\psi$
has the following key properties (assume a $L$ site open chain): (1) it
commutes with the Hamiltonian (at least for $L\gg 1$) $[{{\cal{H}}},\psi]\to
0$; (2) it anti-commutes with the discrete symmetry (here the parity
$\left\\{{{\mathbb{P}}},\psi\right\\}=0$); and (3) it is normalizable
$\psi^{\dagger}\psi^{\vphantom{\dagger}}=\psi^{2}=1$.

Introducing two Majorana fermions at each site
$a_{j}=c^{\dagger}_{j}+c^{\vphantom{\dagger}}_{j}$ and
$b_{j}={\rm{i}}(c^{\dagger}_{j}-c^{\vphantom{\dagger}}_{j})$, the above XY-
Kitaev chain model rewrites

${\cal{H}}_{\rm
K}={\rm{i}}\sum_{j}\left(X_{j}b_{j}a_{j+1}-Y_{j}a_{j}b_{j+1}+h_{j}a_{j}b_{j}\right).$
(3)

Assuming site-independent couplings ($h/X\geq 0$, $X>Y\geq 0$), one can
construct two such strong MZM operators

$\psi_{a}=\frac{1}{N_{a}}\sum_{j=1}^{L}\Theta_{j}^{a}a_{j}\quad{\rm{and}}\quad\psi_{b}=\frac{1}{N_{b}}\sum_{j=1}^{L}\Theta_{j}^{b}b_{L+1-j},$
(4)

which both commute with ${\cal{H}}_{\rm K}$ in the $L\to\infty$ limit under
the simple condition $h<X+Y$ Chepiga and Laflorencie (2023a), see Fig. 1 (a).
They decay exponentially away from the left and right boundaries

$\left|\Theta_{j}^{a,b}\right|\propto\exp\left(-\frac{j-1}{\xi_{\rm
zm}}\right).$ (5)

The MZM localization length diverges if $X+Y-h\to 0^{+}$, following $\xi_{\rm
zm}\approx\frac{X-Y}{X+Y-h}$ Chepiga and Laflorencie (2023a), thus ensuring
their normalization $\psi_{a,b}^{2}=1$, with a finite norm

$N_{a,b}=\sqrt{{\sum_{j=1}\left|\Theta_{j}^{a,b}\right|^{2}}}<\infty\quad{\rm{if}}\quad
h<X+Y.$ (6)

The strong character of the MZMs in Eq. (4) is rooted in the anti-commutation
property $\left\\{{{\mathbb{P}}},\psi_{a,b}\right\\}=0$ not (a). Indeed, when
combined with $[{\cal H}_{\rm K},\psi_{a,b}]=0$, we can build the following
non-local (bilocalized at both edges) Dirac fermion operator

$\Psi^{\dagger}=\frac{1}{2}\left(\psi_{a}-{\rm{i}}\psi_{b}\right)$ (7)

which creates a zero-energy fermion, and provides a mapping between the two
parity sectors for all states. In the thermodynamic limit, the topological
regime is characterized by

Figure 1: (a) Phase diagram of the clean Kitaev-XY chain model Eqs. (1)-(2).
The Ising transition line $h=X+Y$ (white line) separates topological and
trivial regimes. The color map describes the inverse norm squared of the MZM
Eq. (6) with $\mathsf{N}_{0}=N_{a,b}$. The commensurate-incommensurate line
inside the topological regime is shown by the dotted line. Panels (b-c) are
scans along the $Y=0$ line of the phase diagram (corresponding to the
transverse-field Ising model) showing exact diagonalization (ED) results for
open chains of various sizes $L$, as indicated on the plot. (b) The MZM
occupation number ${\cal{N}}_{\rm MZM}=\langle\Psi^{\dagger}\Psi\rangle$,
computed in the ground-state (full lines) and the first excited state (dotted
lines), takes asymptotic values 0 or 1 in the ordered phase ($h<X$), $0.5$ on
the disordered side ($h>X$), and non-trivial universal values at the quantum
critical point (QCP) ${{\cal N}}_{\rm MZM}=1/2-4{\cal{G}}/\pi^{2}\approx
0.1287731$ for the GS and ${{\cal N}}_{\rm
MZM}=1/2+(8-4{\cal{G}})/\pi^{2}\approx 0.939342$ for the first excited state,
where ${\cal G}\approx 0.915966$ is the Catalan’s constant (see text). (c) The
MZM fidelity ${\cal{F}}_{\rm MZM}$ Eq. (9) is shown across the transition for
the same system sizes. It takes asymptotic values 1 or 0 in the two phases,
and a non-trivial universal number ${\cal{F}}_{\rm MZM}^{\rm
critical}=\sqrt{8}/\pi\approx 0.900316$ at the Ising QCP. The three insets
(i-iii) show the bottom of the many-body spectrum for $L=16$ in the three
regimes, resolved in term of the parity quantum number $p=\pm 1$.

${\langle\,n_{p}\,|}\Psi^{\dagger}\Psi^{\vphantom{\dagger}}{|\,n_{p}\,\rangle}=\begin{cases}1&\Rightarrow\quad\Psi{|\,n_{p}\,\rangle}={|\,n_{-p}\,\rangle}\\\
0&\Rightarrow\quad\Psi^{\dagger}{|\,n_{p}\,\rangle}={|\,n_{-p}\,\rangle},\end{cases}$
(8)

for any many-body eigenstate ${|\,n_{p}\,\rangle}$ of parity $p=\pm 1$, where
all the energies are pairwise degenerate $E_{n}^{+}=E_{n}^{-}$.

Fig. 1 summarizes these classic results for the clean Kitaev-XY chain, where
two topological markers are shown: the zero-mode occupation number ${{\cal
N}}_{\rm MZM}=\langle\Psi^{\dagger}\Psi^{\vphantom{\dagger}}\rangle$ for both
the ground-state (GS) and the first excited state (FES), and the MZM fidelity,
defined by

${\cal{F}}_{\rm
MZM}^{(n)}=\frac{1}{2}{\langle\,n_{-p}\,|}\psi_{a}\pm{\rm{i}}\psi_{b}{|\,n_{p}\,\rangle},$
(9)

which quantifies the connection between the two parity sectors via the MZM
mapping Eq. (8).

Main results and paper outline— In this work we present analytical
calculations of these two topological quantities. Along the Ising quantum
critical line $h=X+Y$, despite its non-topological nature, we surprisingly
find a finite universal value for the fidelity ${\cal{F}}_{\rm
MZM}^{(n)}=\sqrt{8}/\pi,\,\forall n$, as well as for the two lowest MZM
occupation numbers, ${{\cal N}}_{\rm MZM}=1/2-4{\cal{G}}/\pi^{2}\approx
0.1287731$ for the ground-state, and ${{\cal N}}_{\rm
MZM}=1/2+(8-4{\cal{G}})/\pi^{2}\approx 0.939342$ for the first excited state,
where ${\cal G}\approx 0.915966$ is the Catalan’s constant Catalan (1865). In
the rest of the paper, we present the analytical derivations, that we then
carefully check numerically using large scale exact diagonalization up to
$L\sim 10^{4}$ lattice sites. We also prove analytically the universality of
our results along the Ising quantum critical line of the clean Kitaev-XY chain
model. Interestingly, we identify a special critical point at $h=X+Y=\sqrt{2}$
where finite corrections vanish completely.

##  II Majorana zero mode fidelity in the clean Kitaev chain

###  II.1 Relation between the fidelity and the parity gap

We first consider a uniform chain of $L$ sites with free ends. The parity gap
is the energy difference within the even-odd doublet, defined by

$\Delta_{\rm
parity}^{(n)}={\cal{S}}_{n}\Bigl{(}\langle\,n_{-}\,|{\cal{H}}|\,n_{-}\,\rangle-\langle\,n_{+}\,|{\cal{H}}|\,n_{+}\,\rangle\Bigr{)},$
(10)

where the sign ${\cal S}_{n}$ ensures that $\Delta_{\rm parity}^{(n)}$ is
positive. In fact, this term has a physical meaning: it reflects the sign of
the end-to-end correlations between the boundary spins not (b).

The even-odd mapping is rigorously given by

$\frac{\psi_{a}+{\rm{i}}{\cal
S}_{n}\mathbb{P}\psi_{b}}{2}|\,n_{p}\,\rangle={\cal F}_{\rm
MZM}^{(n)}|\,n_{-p}\,\rangle-\sum_{n^{\prime}\neq
n}\alpha_{n^{\prime}}|\,n^{\prime}_{-p}\,\rangle,$ (11)

yielding the following expression for the parity gap

$\Delta_{\rm parity}^{(n)}=\frac{{\cal{S}}_{n}\langle\,n_{-}\,|\left[{\cal
H},\psi_{a}\right]|\,n_{+}\,\rangle-{\rm{i}}\langle\,n_{-}\,|\left[{\cal
H},\psi_{b}\right]|\,n_{+}\,\rangle}{2{\cal F}_{\rm
MZM}^{(n)}}~{}~{}~{}~{}~{}\,\,$ (12)

where one recognizes the commutators between the Hamiltonian and the MZMs
$\left[{\cal{H}},\psi_{a,b}\right]$. This expression will be used below to get
exact forms of the MZM fidelity.

###  II.2 The transverse field Ising chain: analytical results

The special point where hopping equals pairing in the fermionic Kitaev chain
Eq. (1) corresponds to the paradigmatic transverse field Ising (TFI) chain
model which is a cornerstone of quantum statistical physics Pfeuty (1970);
Suzuki  _et al._ (2013):

$\displaystyle{\cal{H}}_{\rm TFI}$ $\displaystyle=$
$\displaystyle-\sum_{j}\left(X\sigma_{j}^{x}\sigma_{j+1}^{x}+h\sigma_{j}^{z}\right)$
(13) $\displaystyle=$
$\displaystyle{\rm{i}}\sum_{j}\left(Xb_{j}a_{j+1}+h_{j}a_{j}b_{j}\right).$

Its relative simple form will allow us to provide some details on the
analytical calculation of ${\cal{F}}_{\rm MZM}^{(n)}$. The central quantities
are the MZM commutators with $\cal H_{\rm TFI}$, here given by

$\left[{\cal{H}}_{\rm
TFI},\psi_{a}\right]=-\frac{2{\rm{i}}X}{{{N}}_{a}}\Gamma^{L}b_{L},\quad\left[{\cal{H}}_{\rm
TFI},\psi_{b}\right]=\frac{2{\rm{i}}X}{{{N}}_{b}}\Gamma^{L}a_{1},$ (14)

where $\Gamma=h/X$ is the control parameter for the transition. If one
expresses the boundary operators in the spin language $a_{1}=\sigma_{1}^{x}$
and $b_{L}=-{\rm{i}}\mathbb{P}\sigma_{L}^{x}$, the parity gap Eq. (12) can be
rewritten as follows

$\displaystyle\Delta_{\rm parity}^{(n)}$ $\displaystyle=$
$\displaystyle\frac{X}{{\cal F}_{\rm
MZM}^{(n)}}\Bigl{(}\frac{m_{1}^{s}}{N_{b}}+\frac{m_{L}^{s}}{N_{a}}\Bigr{)}\Gamma^{L}.$
(15)

Here we have introduced the surface magnetization Peschel (1984); Karevski
(2000)
$m_{1,L}^{s}=\left|\langle\,n_{-}\,|\sigma_{1,L}^{x}|\,n_{+}\,\rangle\right|$
whose magnitude is independent of the state $n$ in the free-fermion case, so
is the parity gap (see A+). Moreover, the MZM norms Eq. (6), both equal for
clean chains, are given by

$N_{a,b}=\sqrt{\frac{1-\Gamma^{2L}}{1-\Gamma^{2}}}\equiv{\mathsf{N}_{0}},$
(16)

which hence leads to the following expression for the (energy-independent) MZM
fidelity of finite clean TFI chains

${\cal{F}}_{\rm MZM}(L)=\frac{2Xm^{s}}{\Delta_{0}\mathsf{N}_{0}}\Gamma^{L},$
(17)

where $\Delta_{0}$ is the lowest energy gap, $\mathsf{N}_{0}=N_{a,b}$ is the
MZM norm in Eq. (16), and $m^{s}=m_{1,L}^{s}$ is the surface magnetization,
similar at each boundary of a clean chain. Eq. (17) turns out to be valid
across the entire phase diagram for finite chains, as we discuss now for the
three physical regimes.

####  II.2.1 Disordered phase ${\Gamma>1}$

The trivial (disordered) phase for $h>J$ has a finite energy gap above the GS,
$\Delta_{0}=2X(\Gamma-1)$, a diverging MZM norm
${\mathsf{N}}_{0}\sim\Gamma^{L-1}$, and a power-law vanishing surface
magnetization $m^{s}\sim L^{-3/2}$ Pfeuty (1970). This, when put together,
leads to the following power-law decay for the MZM fidelity in the
topologically trivial phase

${\cal{F}}_{\rm MZM}^{\rm trivial}(L)\sim L^{-3/2},$ (18)

which asymptotically matches the expected vanishing in the disordered phase.
However, the relatively slow algebraic decay is not a trivial result, in a
regime where one would have naively expected a faster exponential decay, see
below Sec. II.3 for a numerical check of this result.

####  II.2.2 Ordered regime ${\Gamma<1}$

In contrast, in the topological (ordered) regime one can show that
$m^{s}=1/\mathsf{N}_{0}=\sqrt{1-\Gamma^{2}}$. Moreover, the gap is
exponentially small Pfeuty (1970); Cabrera and Jullien (1987); Campostrini
_et al._ (2015)
$\Delta_{0}=2X(1-\Gamma^{2})\Gamma^{L}+{\cal{O}}(\Gamma^{2L})$, which yields

${\cal{F}}_{\rm MZM}^{\rm topo}(L)=1-{\cal{O}}(\Gamma^{2L}),$ (19)

as expected from the strong MZM definition. The finite-size corrections
$\sim\exp(-2L/\xi_{\rm zm})$ are exponentially small, controlled by the MZM
localization length

$\xi_{\rm zm}=\frac{1}{\ln\left(1/\Gamma\right)}.$ (20)

Figure 2: Finite-size scaling of the MZM fidelity for the TFI chain model
across the three different regimes. Comparison between analytical predictions
(full lines) and exact diagonalization (ED) results (symbols) obtained from
free-fermion calculations, see Sec. II.3. (a) In the topological regime, one
verifies the very fast convergence ${\cal{F}}_{\rm MZM}^{\rm topo}\to 1$ with
$L$, following Eq. (19), displayed for various values of $\Gamma=h/X<1$: lines
are $\Gamma^{2L}$ and symbols show ED data. (b) At criticality $\Gamma=1$, the
$1/L$ convergence towards ${\sqrt{8}}/{\pi}$ shows perfect agreement between
ED data (symbols) and the analytical prediction Eq. (23) (red line). (c) In
the disordered regime $\Gamma>1$ the expected analytical scaling Eq. (18)
provides a perfect description of ED data.

####  II.2.3 Quantum critical point ${\Gamma=1}$

Perhaps the most remarkable result concerns the quantum critical point (QCP)
itself. Indeed, for $\Gamma=1$ both the finite-size gap and the surface
magnetization vanish, $\Delta_{0}\sim 1/L$ and $m^{s}\sim 1/\sqrt{L}$, while
the MZM norm $\mathsf{N}_{0}=\sqrt{L}$. This implies for Eq. (17) a finite
critical fidelity, a results already visible in Fig. 1 (c) where a finite-size
crossing occurs. One can be more precise, using the critical scaling of the
gap Cabrera and Jullien (1987)

$\Delta_{0}(L)=Xk_{\rm min}(L),$ (21)

where $k_{\rm min}$ is the lowest mode (determined from the boundary
conditions) that one can write

$k_{\rm min}(L)=\frac{\pi}{L+\ell_{\rm eff}}.$ (22)

OBC imply Cabrera and Jullien (1987) $\sin k(L+1)=-\sin kL$, which yields
$\ell_{\rm eff}=\frac{1}{2}$.

The critical scaling of the surface magnetization can be obtained when the QCP
is approached from the ordered side $\Gamma\to 1^{-}$, where
$m^{s}\approx\sqrt{2/\xi_{\rm zm}}$ for large enough MZM localization length.
At the QCP where $\xi_{\rm zm}$ is formally infinite, it is standard to
replace this length scale by the lattice size $L$. Interestingly we
numerically find that the correct length scale to use is precisely the one
which enters in the wave vector $k_{\rm min}$, i.e. $L+\ell_{\rm eff}=L+1/2$
(see Appendix B). When plugged into Eq. (17), we then arrive at a rather
simple expression for the critical fidelity

${\cal{F}}_{\rm MZM}^{\rm
critical}(L)=\frac{\sqrt{8}}{\pi}+\left(\pi\sqrt{2}L\right)^{-1}+{\cal{O}}\left(L^{-2}\right),$
(23)

which converges to $\sqrt{8}/\pi\approx 0.90032$, with finite-size algebraic
corrections $\sim L^{-1}$.

###  II.3 Numerical results

These analytical results can be checked against exact diagonalization
simulations of the free fermionic Hamiltonian that we perform for open chains
of large length $L$. Fig. 2 shows such exact finite-size computations of the
MZM fidelity for the three different regimes where we nicely observe that the
analytical predictions for the trivial regime Eq. (18), the topological phase
Eq. (19), and at criticality Eq. (23), all perfectly agrees with the exact
numerics. In Appendix A we provide some details about the free-fermion exact
diagonalization method, and how the MZM operators and fidelity are obtained.

##  III Universal fidelity at Ising criticality

###  III.1 Analytical results

We now discuss how universal this result is: we repeat a similar calculation
for the more generic Kitaev-XY model Eq. (3), moving away from the $t=\Delta$
point along the Ising quantum critical line at $h=X+Y$, see Fig. 1 (a).

####  III.1.1 Derivation of the critical fidelity

We restart from the definition of the parity gap Eq. (12). The first important
quantity to compute are the commutators
$\left[{\cal{H}}_{K},\psi_{a,b}\right]$, which are straightforward to get at
criticality

$\left[{\cal{H}},\psi_{a}\right]=\frac{-2{\rm{i}}X^{2}}{{{N}}_{a}(X-Y)}b_{L}\,;\,\left[{\cal{H}},\psi_{b}\right]=\frac{2{\rm{i}}X^{2}}{{{N}}_{b}(X-Y)}a_{1},$
(24)

where the norm is

$N_{a,b}={\mathsf{N}}_{0}=\frac{X}{X-Y}\sqrt{L-\ell^{\prime}_{\rm eff}}$ (25)

with (see Appendix C.1.2)

$\ell^{\prime}_{\rm eff}=\frac{y(2+y)}{1-y^{2}},\quad(y=Y/X).$ (26)

Injecting this expression in Eq. (12) we arrive at

${\cal{F}}_{\rm MZM}^{\rm
critical}(L)=\frac{2Xm^{s}}{\Delta_{0}\sqrt{L-\ell^{\prime}_{\rm eff}}}\,,$
(27)

where as before $\Delta_{0}$ is the lowest gap, and $m^{s}$ the surface
magnetization, both evaluated along the Ising quantum critical line. There,
one can use Campostrini  _et al._ (2014)

$\Delta_{0}(L)\approx(X-Y)k_{\rm min}=\frac{(X-Y)\pi}{L+\ell_{\rm eff}}$ (28)

and similar to the TFI case (see Appendix B.3)

$m^{s}(L)\approx\frac{X-Y}{X}\sqrt{\frac{2}{L+\ell_{\rm eff}}},$ (29)

which gives after replacing in Eq. (27)

$\displaystyle{\cal{F}}_{\rm MZM}^{\rm critical}(L)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{8}}{\pi}\sqrt{\frac{L+\ell_{\rm
eff}}{L-\ell^{\prime}_{\rm eff}}}$ (30) $\displaystyle=$
$\displaystyle\frac{\sqrt{8}}{\pi}+\frac{\sqrt{2}\left(\ell_{\rm
eff}+\ell^{\prime}_{\rm
eff}\right)}{\pi}L^{-1}+{\cal{O}}\left(L^{-2}\right).~{}~{}~{}$ (31)

Figure 3: Finite-size scaling of the MZM fidelity for the Kitaev-XY model Eq.
(3). Free-fermion ED results are shown as a function of $Y/X$ and $h/X$, see
the phase diagram in Fig. 1 (a). (a-b) ${\cal{F}}_{\rm MZM}$ across the
transition along two horizontal cuts (a) or vertical cuts (b) for 3
representative system sizes: one sees the curves crossing at $\sqrt{8}/\pi$ in
all cases. (c) ${\cal{F}}_{\rm MZM}$ is shown along the critical line $h=X+Y$
as a function of $Y$ for various $L$. One sees the finite-size convergence
towards $\sqrt{8}/\pi$ with a vanishing of finite-size corrections at
$Y/X=\sqrt{2}-1$ (dashed vertical line). Panel (d) magnifies this regime where
ED data for the largest chains (symbols $L=512,\ldots,16384$) are compared to
the analytical expression Eq. (30) (full line) with $\ell^{\prime}_{\rm eff}$
given by Eq. (26) and $\ell_{\rm eff}$ by Eq. (32). Panel (e) tests the
validity of the ansatz Campostrini  _et al._ (2014) Eq. (32) for $\ell_{\rm
eff}$ using two estimates: The lowest energy gap $\Delta_{0}(L)$ Eq. (28) and
the surface magnetization $m^{2}(L)$ Eq. (29). The agreement bewteeen ED data
(symbols) and the analytical expression Eq. (32) (line) gets clearly better
when $L$ grows.

####  III.1.2 Finite-size corrections towards the universal value

Here a few comments are in order. First we find that the asymptotic critical
value of the MZM fidelity $\sqrt{8}/\pi$ turns out to be universal along the
Ising quantum critical line. Then, the leading finite size corrections are
proportional to $L^{-1}$, with a prefactor which involves $\ell^{\prime}_{\rm
eff}$ and $\ell_{\rm eff}$. We can check that the TFI result Eq. (23) is
perfectly recovered using $\ell_{\rm eff}=1/2$ and $\ell^{\prime}_{\rm eff}=0$
at $Y=0$. For finite $Y$ however, it is more cumbersome to explicitly compute
$\ell^{\prime}_{\rm eff}$. Interestingly, Ref. Campostrini  _et al._ (2014)
gave the following ansatz for this effective length shift

$\ell_{\rm eff}=\frac{1}{2}-\frac{y(3+y)}{1-y^{2}}.$ (32)

Plugging this ansatz in Eq. (30) we see that the finite size corrections
change sign and cancel out exactly for $y+1=\sqrt{2}$, i.e. when
$h_{c}/X=\sqrt{2}$.

###  III.2 Numerical results

####  III.2.1 Ising transitions

The MZM fidelity is obtained from ED calculations, for various chains with
OBC, typically ranging from $L=16$ to $L=16384$ sites. Fig. 3 (a-b) show
${\cal{F}}_{\rm MZM}$ for various cuts in the phase diagram of the Kitaev-XY
Hamiltonian, see Fig. 1 (a). Exactly as was observed for the TFI chain in Fig.
1 (c), here also the transition between trivial (${\cal{F}}_{\rm MZM}\to 0$)
and topological (${\cal{F}}_{\rm MZM}\to 1$) regimes is signalled by a finite-
size crossing at ${\cal{F}}_{\rm MZM}=\sqrt{8}/\pi$, thus confirming this
universal number along the entire Ising critical line of the XY-Kitaev model.

####  III.2.2 Finite-size convergence

Fig. 3 (c-d) show how the critical fidelity approaches this universal value.
In panel (c) one sees a clear convergence with $L$ to ${\sqrt{8}}/{\pi}$,
which features distinct finite-size effects as a function of $y=Y/X$. Indeed,
as expected from the previous part, the finite-size corrections change sign
for $y=\sqrt{2}-1$. This is better seen in panel (d) in a log scale where the
numerical data perfectly compare to the analytical expression Eq. (30). In
terms of the Kitaev chain parameters, this occurs for a particular ratio
between pairing and hopping $\Delta/t=\sqrt{2}-1$.

####  III.2.3 Effective length scales

The vanishing of finite-size corrections is due to the fact that the effective
length scales emerging in the norm $\ell^{\prime}_{\rm eff}$ and from the gap
$\ell_{\rm eff}$ cancel each other in Eq. (30), provided that $\ell_{\rm
eff}+\ell^{\prime}_{\rm eff}=0$. While the analytical expression for
$\ell^{\prime}_{\rm eff}$ in Eq. (26) can be derived exactly from the
normalization of the zero mode operator (see Appendix C.1.2), we have tested
the ansatz Eq. (32) proposed by Campostrini and co-workers in Ref. Campostrini
_et al._ (2014), which indeed perfectly matches our ED data. This is shown in
Fig. 3 (e) where $1/2-\ell_{\rm eff}=y(3+y)/(1-y^{2})$ is plotted together
with ED estimates extracted from the gap Eq. (28) and the surface
magnetization Eq. (29). The agreement is excellent and gets clearly better at
large sizes when $y\to 0$.

##  IV The zero-mode occupation

###  IV.1 Simple expectation in the topological regime

The MZM occupation number for an eigenstate ${|\,n_{p}\,\rangle}$ is

$\displaystyle{\cal{N}}_{\rm
MZM}^{(n_{p})}=\langle\Psi^{\dagger}\Psi\rangle_{n_{p}}=\frac{1}{2}\Bigl{(}1+{\rm{i}}\langle\psi_{a}\psi_{b}\rangle_{n_{p}}\Bigr{)},$
(33)

In the topological regime, the parity fidelity is ${\cal F}_{\rm MZM}=1$ and
therefore Eq. (11) becomes

$|\,n_{-p}\,\rangle=\frac{1}{2}\left(\psi_{a}+{\rm{i}}{\cal
S}_{n}\mathbb{P}\psi_{b}\right)|\,n_{p}\,\rangle,$ (34)

where ${\cal{S}}_{n}=\pm 1$ encodes the sign of the correlation between the
boundary spins in the state ${|\,n_{p}\,\rangle}$. It is then straightforward
to show that for any eigenstate ${|\,n_{p}\,\rangle}$ the MZM occupation

${\cal{N}}_{\rm MZM}^{(n_{p})}=\frac{1}{2}\left(1-p{\cal{S}}_{n}\right).$ (35)

The 4 possible cases are summarized in Table 1.

Parity $p$ | End-to-end correlation $\cal S$ | MZM occupation ${\cal N}_{\rm MZM}$
---|---|---
+1 |  $-1$ (AF) | 1
$-1$ | +1 (FM) | 1
+1 | +1 (FM) | 0
$-1$ |  $-1$ (AF) | 0
Table 1: The occupation ${\cal N}_{\rm MZM}^{(n_{p})}=0~{}{\rm or}~{}1$ in the
topological regime, depends on the parity $p$ and the sign of edge spins
${\cal{S}}_{n}={\rm{sgn}}\left(C^{\rm end}_{xx}\right)$.

A simple example is the ferromagnetic TFI chain which, in the limit of large
coupling and small field limit $\Gamma=h/J\ll 1$, displays eigenstates with
cat-states forms (in the $\\{\sigma^{x}\\}$ basis)

$\displaystyle{|\,n_{p}^{\rm(FM)}\,\rangle}$ $\displaystyle\approx$
$\displaystyle\frac{{|\,\hskip-2.84544pt\uparrow\uparrow\downarrow\uparrow\uparrow\ldots\downarrow\uparrow\,\rangle}+p{|\,\hskip-2.84544pt\downarrow\downarrow\uparrow\downarrow\downarrow\ldots\uparrow\downarrow\,\rangle}}{\sqrt{2}}$
(36) $\displaystyle{|\,n_{p}^{\rm(AF)}\,\rangle}$ $\displaystyle\approx$
$\displaystyle\frac{{|\,\hskip-2.84544pt\uparrow\uparrow\downarrow\uparrow\uparrow\ldots\downarrow\downarrow\,\rangle}+p{|\,\hskip-2.84544pt\downarrow\downarrow\uparrow\downarrow\downarrow\ldots\uparrow\uparrow\,\rangle}}{\sqrt{2}},$
(37)

where one sees that the sole difference concerns the edge spin correlations:
ferromagnetic (FM) or antiferromagnetic (AF). Since the MZMs are essentially
localized on the boundary operators $\psi_{a}\approx a_{1}=\sigma_{1}^{x}$ and
$\psi_{b}\approx b_{L}=-{\rm{i}}{\mathbb{P}}\sigma_{L}^{x}$, one can rewrite
the MZM occupation Eq. (33) as follows

${\cal{N}}_{\rm
MZM}^{(n_{p})}\approx\frac{1}{2}\left({1}-\langle{\mathbb{P}}\sigma_{1}^{x}\sigma_{L}^{x}\rangle_{n_{p}}\right),$
(38)

which nicely matches Eq. (35) for the above cat-state situation.

However, in the general case the MZM correlator
${\rm{i}}\langle\psi_{a}\psi_{b}\rangle_{n_{p}}$ is not a simple object to
directly evaluate, except very deep in the topological regime where it reduces
to the end-to-end correlation, which for cat-states is
$\langle\sigma_{1}^{x}\sigma_{L}^{x}\rangle_{n_{p}}=\pm 1$. Away from the
$\Gamma\to 0$ limit, it is quite interesting to notice that despite its much
more complicated non-local form

$\displaystyle{\rm{i}}\langle\psi_{a}\psi_{b}\rangle_{n_{p}}$ $\displaystyle=$
$\displaystyle\frac{{\rm{i}}}{N_{a}N_{b}}\sum_{j,j^{\prime}}\Theta_{j}^{a}\Theta_{j^{\prime}}^{b}\langle
a_{j}b_{L+1-j^{\prime}}\rangle_{n_{p}},$ (39)

it is expected to be ${\rm{i}}\langle\psi_{a}\psi_{b}\rangle_{n_{p}}=\pm 1$.
This will no longer be true in the non-topological regime where the norms
$N_{a,b}$ diverge, and therefore $\langle\psi_{a}\psi_{b}\rangle_{n_{p}}\to
0$, yielding ${\cal{N}}_{\rm MZM}\to 1/2$, as observed in Fig. 1 (b).

###  IV.2 Free fermion formulation

The non-interacting Kitaev chain Hamiltonian is easily diagonalized by a
Bogoliubov transformation (see Appendix), and takes the following quadratic
form

${\cal{H}}_{\rm
K}=2\sum_{m=1}^{L}\epsilon_{m}\left(\phi^{\dagger}_{m}\phi_{m}^{\vphantom{\dagger}}-\frac{1}{2}\right),$
(40)

where $\phi_{m}$ are new fermionic modes, and the single particle energies are
such that $0\leq\epsilon_{1}\leq\epsilon_{2}\leq\ldots\epsilon_{L}$. In the
topological regime, for large enough chain length $L$ one expects
${\cal{F}}_{\rm MZM}=1$, which implies that
$\Psi^{\dagger}=\phi^{\dagger}_{1}$. This will essentially be true for system
sizes much larger that the correlation length $\xi_{\rm zm}$ (i.e. the
localization length of the MZM). However, when $\xi_{\rm zm}$ is not small as
compared to the chain length $L$, the situation becomes quite interesting.
Indeed, in this case the fidelity being different from unity, it is convenient
to rewrite $\Psi^{\dagger}$ as

$\displaystyle\Psi^{\dagger}$ $\displaystyle=$ $\displaystyle{\cal{F}}_{\rm
MZM}\phi_{1}^{\dagger}$ (41) $\displaystyle+$ $\displaystyle\sum_{m\geq
2}\left(\frac{{\cal{A}}_{m}+{\cal{B}}_{m}}{2}\right)\phi^{\dagger}_{m}+\left(\frac{{\cal{A}}_{m}-{\cal{B}}_{m}}{2}\right)\phi^{\vphantom{\dagger}}_{m},$

where ${\cal{A}}_{m}$ and ${\cal{B}}_{m}$ depends on both the MZM coefficients
$\Theta_{j}^{a,b}$ and the Bogoliubov transformation (see Appendix). This
allows us to express the MZM occupation number as follows

$\displaystyle{{\cal N}}_{\rm MZM}^{(n_{p})}$ $\displaystyle=$
$\displaystyle{\cal{F}}_{\rm
MZM}^{2}\langle\phi^{\dagger}_{1}\phi_{1}^{\vphantom{\dagger}}\rangle_{n_{p}}$
$\displaystyle+$ $\displaystyle\sum_{m\geq
2}\left(\frac{1}{4}\left[{\cal{A}}_{m}-{\cal{B}}_{m}\right]^{2}+{\cal{A}}_{m}{\cal{B}}_{m}\langle\phi^{\dagger}_{m}\phi_{m}^{\vphantom{\dagger}}\rangle_{n_{p}}\right).~{}~{}~{}~{}~{}~{}$

Hence, since many-body eigenpairs
$\left\\{{|\,n_{p}\,\rangle}\,;\,{|\,n_{-p}\,\rangle}\right\\}$ only differ by
their occupation of the lowest single-particle mode $m=1$, the difference in
their MZM occupation is simply given by

$\displaystyle\delta{\cal{N}}_{\rm MZM}^{(n_{p})}$ $\displaystyle=$
$\displaystyle{\langle\,n_{p}\,|}\Psi^{\dagger}\Psi{|\,n_{p}\,\rangle}-{\langle\,n_{-p}\,|}\Psi^{\dagger}\Psi{|\,n_{-p}\,\rangle}$
(43) $\displaystyle=$ $\displaystyle-p{\cal{S}}_{n}{\cal{F}}_{\rm MZM}^{2},$
(44)

where the sign prefactor depends on both the parity $p$ and the the sign
${\cal{S}}_{n}$ of the end-to-end spin correlation. This result is generally
true across the entire phase diagram (even in the trivial regime, but only for
finite $L$ because otherwise the strong MZM operator is not defined anymore).

This can be numerically checked, see Fig. 4 where the MZM occupations have
been computed for random high-energy eigenstates, with or without the $m=1$
mode occupied. While ${\cal{N}}_{\rm MZM}\to 0~{}{\rm{or}}~{}1$ in the
topological regime for all energies, results are much more spread in the
trivial regime where ${\cal{N}}_{\rm MZM}\to 1/2$ slowly with increasing $L$.
We also nicely check in Fig. 4 (c) that the difference Eq. (43) does not
depend on the energy and is given by the square of the fidelity Eq. (44).

Figure 4: ED results for the MZM occupation number ${\cal{N}}_{\rm MZM}$
computed for the TFI model ($Y=0$) against $h/X$, shown for 4 different system
sizes. All symbols correspond to several hundreds of random (high-energy)
eigenstates of the form
${|\,n_{p}\,\rangle}=\prod_{m=1}^{L}\Upsilon_{m}\phi^{\dagger}_{m}{|\,\rm
GS\,\rangle}$ where $\Upsilon_{m}=0~{}{\rm{or}}~{}1$ with probablity $1/2$,
with the lowest fermionic mode $m=1$ being either (a) occupied or empty (b).
The full blue line shows results for the two lowest energy states with $L=128$
sites: (a) for the first excited state ${|\,\rm
FES\,\rangle}=\phi^{\dagger}_{1}{|\,\rm GS\,\rangle}$, and (b) for the ground-
state ${|\,\rm GS\,\rangle}$. (c) The absolute value of the difference Eq.
(43) is independent of the states, and the analytical prediction Eq. (44) is
perfectly verified in the inset.

###  IV.3 Analytical derivation at criticality

Apart from the simple exact relation Eq. (44) which links the difference
between MZM occupations and the fidelity, there is no easy way to analytically
evaluate ${\cal{N}}_{\rm MZM}$ for a generic eigenstate having the following
form in the fermionic basis

${|\,n_{p}\,\rangle}=\prod_{m~{}{\rm{occupied}}}\phi^{\dagger}_{m}{|\,\rm
GS\,\rangle}.$ (45)

Nevertheless, at criticality, using Eq. (33) and Eq. (39), it is
straightforward to arrive for the TFI chain model at

$\displaystyle{\cal{N}}_{\rm critical}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigl{(}1+{\rm{i}}\langle\psi_{a}\psi_{b}\rangle\Bigr{)},$
$\displaystyle=$
$\displaystyle\frac{1}{2}\Bigl{(}1+\frac{\rm{i}}{L}\sum_{i=1}^{L}\sum_{j=1}^{L}\langle
a_{i}b_{j}\rangle\Bigr{)}=\frac{1}{2}+\frac{1}{L}\sum_{ij}\langle
c^{\dagger}_{i}c^{\vphantom{\dagger}}_{j}\rangle,$

where we have used the simple form of the critical MZM at the Ising QCP of the
TFI chain

$\psi_{a}=\frac{1}{\sqrt{L}}\sum_{j=1}^{L}a_{j}\quad{\rm{and}}\quad\psi_{b}=\frac{1}{\sqrt{L}}\sum_{j=1}^{L}b_{j}.$
(47)

It turns out that one can obtain a rather simple analytical expression for the
above sum, building on the fact that an $L$-sites XY chain is equivalent to
two decoupled Ising chains with $L/2$ sites Fisher (1994); Iglói and Juhász
(2008). After a bit of manipulation we arrive at

$\displaystyle\sum_{ij}\langle
c^{\dagger}_{i}c^{\vphantom{\dagger}}_{j}\rangle$ $\displaystyle\approx$
$\displaystyle-\frac{2}{L}\sum_{j=1}^{L/2}\frac{j}{\sin\left(\frac{\pi
j}{L}\right)}$ (48) $\displaystyle\xrightarrow[L\to\infty]{}$
$\displaystyle-\frac{2}{\pi^{2}}\int_{0}^{\pi/2}\frac{x}{\sin(x)}{\rm{d}}x$
(49)

where one recognizes the Catalan’s constant $\cal G$ Catalan (1865), given by

${\cal{G}}=\frac{1}{2}\int_{0}^{\pi/2}\frac{x}{\sin(x)}{\rm{d}}x\approx
0.915965594177\ldots$ (50)

Therefore, the GS occupation of the MZM at criticality is expected to be

${\cal{N}}_{\rm
critical}^{\rm(GS)}=\frac{1}{2}-\frac{4}{\pi^{2}}{\cal{G}}\approx
0.128773127289\ldots,$ (51)

for the GS, and for the first excited state ${|\,\rm
FES\,\rangle}=\phi^{\dagger}_{1}{|\,\rm GS\,\rangle}$

${\cal{N}}_{\rm critical}^{\rm(FES)}={\cal{N}}_{\rm
critical}^{\rm(GS)}+{\cal{F}}_{\rm MZM}^{2}\approx 0.939342596428\ldots$ (52)

###  IV.4 Finite size numerics

These analytical predictions for the critical MZM occupation numbers are
numerically checked, both in the GS and the FES. Fig. 5 shows ED results for
critical TFI chains, up to $L=16384$ sites, where we nicely observe a finite-
size convergence to the predicted values.

Figure 5: Critical MZM occupation number ${\cal{N}}_{\rm MZM}$ computed for
the GS (red, left axis) and the FES (blue, right) of the TFI model ($Y=0$) at
criticality ($h=X$). Finite-size ED results (+) are perfectly fitted with a
simple order 2 polynomial form as indicated on the plot (line); the inifinite-
size extrapolations are also indicated on the graph ($\times$).

Interestingly, as before with the fidelity, here we also find along the
critical line an exact vanishing of the finite-size corrections for the MZM
occupations which occurs at the same special point: $\Delta/t=\sqrt{2}-1$ in
the Kitaev chain language (equivalently $Y/X=\sqrt{2}-1$ and $h/X=\sqrt{2}$
for the XY spin chain).

In addition, we numerically observe a weak finite-size drift of this special
point $Y^{*}(L)$ where the finite-size corrections vanish. This is shown in
Fig. 6 where one sees the convergence of $Y^{*}(L)$ towards its asymptotic
value $\sqrt{2}-1$ with $L$ as a power-law $\sim 1/L$ for both quantities.

Figure 6: Main panel: finite-size convergence of ${\cal{F}}_{\rm MZM}^{\rm
critical}$ towards $\sqrt{8}/\pi$ along the critical Ising line. The vanishing
of the finite-size corrections occurs at a length-dependent coupling
$Y^{*}(L)/X$, which quickly converges to $\sqrt{2}-1$ with increasing $L$.
This convergence of $Y^{*}$ vs. $L$ ($X$ is set to 1) is reported in the Inset
for both quantities ${\cal{F}}_{\rm MZM}^{\rm critical}$ and ${\cal{N}}_{\rm
critical}^{(GS)}$. Black lines are fits of the form $A/L$ with
$A_{\cal{F}}\approx 1.608$ and $A_{\cal N}\approx 0.6317$. The convergence is
even faster for $\cal N$.

##  V Conclusions and discussions

###  V.1 Summary of the main results

In this work, we have revisited the topological transition occurring in the
paradigmatic Kitaev chain model. We introduced two topological markers to
probe the Majorana zero mode (MZM) physics: the MZM fidelity ${\cal F}_{\rm
MZM}$, see Eq. (9), and the MZM occupation number ${\cal{N}}_{\rm MZM}$, Eq.
(8). Whereas these two quantities have long been known to take simple values
in both trivial and topological regimes, we have analytically shown that both
markers take a non-trivial universal value, constant along the quantum (1+1)
Ising critical line of the XY-Kitaev chain model, as summarized in Tab. 2.

| Trivial | Critical | Topological
---|---|---|---
${\cal F}_{\rm MZM}$ | 0 | $\frac{\sqrt{8}}{\pi}\approx 0.9003163$ | 1
${\cal{N}}_{\rm MZM}^{\rm(GS)}$ | 0.5 | $\frac{1}{2}-\frac{4{\cal{G}}}{\pi^{2}}\approx 0.1287731$ | 0
${\cal{N}}_{\rm MZM}^{\rm(FES)}$ | 0.5 | $\frac{1}{2}+\frac{8-4{\cal{G}}}{\pi^{2}}\approx 0.9393426$ | 1
Table 2: Summary of the values taken by the two topological markers studied in
this paper: the MZM fidelity ${\cal F}_{\rm MZM}$, see Eq. (9), and the MZM
occupation number ${\cal{N}}_{\rm MZM}$, Eq. (8), for the ground-state (GS)
and the first excited state (FES). They all take non-trivial values at the
Ising quantum critical point between trivial and topological regimes, values
that are universal for (1+1) Ising criticality.

The MZM fidelity, independent of the energy, has been analytically computed by
establishing a connection with the surface magnetization $m^{s}$, the lowest-
energy gap $\Delta_{0}$, and the norm of the commutator
${\cal{C}}_{{}_{\Psi}}=\Bigl{|}\left[{\cal{H}}_{K},\Psi\right]\Bigr{|}$
between the Hamiltonian and the zero-energy Dirac fermion
($\Psi=\left(\psi_{a}+{\rm{i}}\psi_{b}\right)/2$ built from left and right
MZMs), such that ${\cal{F}}_{\rm MZM}={m^{s}{\cal
C}_{{}_{\Psi}}}/{\Delta_{0}}$. These three quantities exhibit algebraic
behavior at criticality, which compensates and gives rise to a universal
value, that we strikingly find to remain constant along the critical Ising
line: ${\cal{F}}_{\rm MZM}^{\rm critical}=\sqrt{8}/\pi$. We have also computed
the finite-size corrections to this asymptotic result: we analytically
identified a singular point along the Ising critical line where such
corrections vanish completely for the special value $\Delta/t=\sqrt{2}-1$ of
the ratio between pairing and hopping in the Kitaev chain. These predictions
were successfully compared with exact diagonalization results (up to several
thousand lattice sites).

Concerning the MZM occupation number ${\cal{N}}_{\rm MZM}$ we have obtained
exact results, universal along the Ising critical line, for two first energy
levels, see Tab. 2, and also for the zero-mode occupation difference
$\delta{\cal{N}}_{\rm MZM}^{(n_{p})}$ between all many-body eigenpairs
$\left\\{{|\,n_{p}\,\rangle}\,;\,{|\,n_{-p}\,\rangle}\right\\}$ that is given
by ${\cal{F}}_{\rm MZM}^{2}$, see Eq. (44).

###  V.2 Open questions and possible future directions

The rather simple analytical expression of the MZM operators Eq. (4) has
clearly helped us to derive exact results at criticality. In addition, the
free-fermion nature allowed to numerically verify our analytical predictions
with great accuracy using very large Kitaev chains, up to $L=16384$ lattice
sites. Nevertheless, there are still many open questions and directions for
which we can imagine several extensions to go beyond the case of non-
interacting clean Kitaev chains.

A first natural development concerns the effects of the environment, such as
in open systems described by non-Hermitian models Kawabata  _et al._ (2019);
Bergholtz  _et al._ (2021), or when quenched disorder is added directly to the
Kitaev Hamiltonian Motrunich  _et al._ (2001); Bravyi and König (2012);
DeGottardi  _et al._ (2013); Hegde and Vishveshwara (2016); Gergs  _et al._
(2016); Monthus (2018a); Lieu  _et al._ (2018); Monthus (2018a); Levy and
Goldstein (2019); Francica  _et al._ (2023). This last case has a particularly
long and dense history, especially for the magnetic version of the problem,
going back to the seminal work of D. S. Fisher Fisher (1992, 1995). Subsequent
progress has been made Young and Rieger (1996); Fisher and Young (1998); Iglói
and Rieger (1998); Refael and Moore (2004); Iglói and Monthus (2005);
Bonesteel and Yang (2007); Fidkowski  _et al._ (2008); Vojta  _et al._ (2009);
Laflorencie (2022); Zabalo  _et al._ (2023), giving a fairly complete
description of the non-interacting random problem, and especially of the very
unusual properties of the infinite randomness criticality fixed point (IRFP)
Fisher (1995, 1994). It would therefore be quite interesting and relevant to
revisit the Kitaev chain model in the presence of quenched disorder using the
topological markers introduced in the present paper, extending the results
obtained at the clean (1+1) Ising critical point to the physics of the IRFP.

Another very important ingredient concerns the effect of interactions on MZMs
Fidkowski and Kitaev (2010, 2011); Hassler and Schuricht (2012); Thomale  _et
al._ (2013); Katsura  _et al._ (2015). In particular it is known that the
construction of MZM operators is a very difficult task in the general
interacting case Goldstein and Chamon (2012); Kells (2015); Fendley (2016);
Kemp  _et al._ (2017); Else  _et al._ (2017); Więckowski  _et al._ (2018);
Yates  _et al._ (2020a). However, a distinction should be made between the
integrable case, e.g. for the XYZ model where an exact construction has been
shown possible in the clean case Fendley (2016), and the non-integrable case,
e.g. provided by the interacting Ising-Majorana chain model Kemp  _et al._
(2017); Else  _et al._ (2017); Kemp (2019); Więckowski and Ptok (2019); Kemp
_et al._ (2020); Yates  _et al._ (2020a); Pandey  _et al._ (2023) where the
MZMs are only almost strong Else  _et al._ (2017); Kemp  _et al._ (2017);
Yates  _et al._ (2020a, b). Nevertheless, it would be very interesting to
consider the possibility of checking in one way or another the universality of
our results against finite interaction along the self-dual Ising critical line
of the interacting Majorana chain model Rahmani  _et al._ (2015); Chepiga and
Laflorencie (2023a).

A clearly challenging direction touches the combined effects of disorder and
interactions which brings a full set of very exciting questions, some of them
being related to the celebrated many-body localization (MBL) problem Alet and
Laflorencie (2018); Abanin  _et al._ (2019). For instance it was recently
found using large-scale DMRG simulations that interactions are not relevant to
the IRFP of the random TFI chain model at zero temperature Chepiga and
Laflorencie (2023b), while Monthus had previously shown Monthus (2018b) that a
strong-disorder RG treatment generates higher-order terms that prevent a
conclusion, in contrast to the XXZ case Fisher (1994). The disordered XYZ
chain was also shown to display very rich physics both at zero temperature
Roberts and Motrunich (2021), and at high energy Slagle  _et al._ (2016); You
_et al._ (2016). In such a context, the possible existence and stability of
MZMs and their fate at criticality in the presence of both disorder and
interactions remains a fascinating subject, as it has been little discussed
and has shown contrasting conclusions Kells  _et al._ (2018); Sahay  _et al._
(2021); Wahl  _et al._ (2022); Laflorencie  _et al._ (2022); Hannukainen  _et
al._ (2023); Moudgalya  _et al._ (2020).

Finally, it is also worth mentioning the very interesting case of driven
systems Khemani  _et al._ (2016); Yates  _et al._ (2018); Bauer  _et al._
(2019) which has strong experimental relevance Yates and Mitra (2021); Mi
(2022); Harle  _et al._ (2023); Ling  _et al._ (2023).

###### Acknowledgements.

I would like to express my sincere thanks to Natalia Chepiga for collaborating
on related matters. I also acknowledge some discussions with Paul Fendley and
Jack Kemp. This work was granted access to the HPC resources of CALMIP center
under the allocation 2022-P0677 as well as GENCI (grant x2021050225).

##  Appendix A Additional free-fermion results

###  A.1 Exact diagonalization of the Kitaev model

In the general case where translational invariance is absent, we use the Nambu
formalism Young and Rieger (1996) in order to solve the (not necessarily
homogeneous) free-fermion Kiteav chain Hamiltonian Eq. (1). Introducing
$\varphi^{\dagger}=\left(c_{1}^{\dagger}\cdots c_{L}^{\dagger}\,c_{1}\cdots
c_{L}\right)$, we arrive at the simple matrix representation of the
Hamiltonian

${\cal{H}}_{\rm K}=\varphi^{\dagger}\begin{pmatrix}M&{\widetilde{M}}\\\
-{\widetilde{M}}&-M\\\ \end{pmatrix}\varphi.$ (53)

$M$ and ${\widetilde{M}}$ are $L\times L$ matrices, such that

$\displaystyle M_{j,j}$ $\displaystyle=$
$\displaystyle\frac{\mu_{j}}{2}=h_{j},$ $\displaystyle M_{j,j+1}$
$\displaystyle=$ $\displaystyle
M_{j+1,j}=-\frac{t_{j}}{2}=-\frac{X_{j}+Y_{j}}{2},$ (54)

while $M_{i,j}=0$ elsewhere (unless periodic boundary conditions are used:
$M_{1,L}=M_{L,1}=p{t_{L}}/{2}$, where $p=\pm 1$ is the fermionic parity). The
matrix ${\widetilde{M}}$ instead is anti-symmetric:

${\widetilde{M}}_{j,j+1}=-{\widetilde{M}}_{j+1,j}=-\frac{\Delta_{j}}{2}=-\frac{X_{j}-Y_{j}}{2},$
(55)

and ${\widetilde{M}}_{i,j}=0\,\,{\rm{otherwise}}$. The free-fermion problem
can be solved by diagonalizing the $2L\times 2L$ Hamiltonian matrix Eq. (53),
and we then rewrite the quadratic Hamiltonian

${\cal{H}}_{\rm
K}=2\sum_{m=1}^{L}\epsilon_{m}\left(\phi^{\dagger}_{m}\phi_{m}^{\vphantom{\dagger}}-\frac{1}{2}\right),$
(56)

with single particle energies
$0\leq\epsilon_{1}\leq\epsilon_{2}\leq\ldots\epsilon_{L}$.

###  A.2 Majorana zero mode operators and fidelity

The new fermionic modes are given by the Bogoliubov transformation
$\phi_{m}^{\dagger}=\sum_{i=j}^{L}\left(u^{m}_{j}\,c_{j}^{\dagger}+v^{m}_{j}\,c_{j}^{\vphantom{\dagger}}\right)$,
with real $u^{m}_{j}$ and $v^{m}_{j}$. The inverse transformation is
$c_{j}^{{\dagger}}=\sum_{m=1}^{L}\left(u_{j}^{m}\phi_{m}^{{\dagger}}+v_{j}^{m}\phi_{m}^{\vphantom{\dagger}}\right)$,
from which one can express the MZM operators Eq. (4)

$\displaystyle\psi_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{a}}\sum_{j=1}^{L}\Theta_{j}^{a}a_{j}=\sum_{m=1}^{L}{\cal{A}}_{m}\left(\phi^{\dagger}_{m}+\phi^{\vphantom{\dagger}}_{m}\right)$
(57) $\displaystyle\psi_{b}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{b}}\sum_{j=1}^{L}\Theta_{j}^{b}b_{L+1-j}={\rm{i}}\sum_{m=1}^{L}{\cal{B}}_{m}\left(\phi^{\dagger}_{m}-\phi^{\vphantom{\dagger}}_{m}\right),$
(58)

where

$\displaystyle{\cal{A}}_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{a}}\sum_{j=1}^{L}\Theta_{j}^{a}\left(u^{m}_{j}+v^{m}_{j}\right)$
(59) $\displaystyle{\rm{and}}\quad{\cal{B}}_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{b}}\sum_{j=1}^{L}\Theta_{j}^{b}\left(u^{m}_{L+1-j}-v^{m}_{L+1-j}\right).$
(60)

The zero-energy Dirac fermion creation operator Eq. (7) is therefore given by

$\Psi^{\dagger}=\sum_{m=1}^{L}\left(\frac{{\cal{A}}_{m}+{\cal{B}}_{m}}{2}\right)\phi^{\dagger}_{m}+\left(\frac{{\cal{A}}_{m}-{\cal{B}}_{m}}{2}\right)\phi^{\vphantom{\dagger}}_{m},$
(61)

from what we can simply express the MZM fidelity. Indeed, since many-body
eigenpairs only differ by their occupation of the lowest single-particle mode
$m=1$ one has

${\cal{F}}_{\rm MZM}=\frac{{\cal{A}}_{1}+{\cal{B}}_{1}}{2},$ (62)

which further simplifies in the clean case where left and right boundaries are
equivalent, yielding

${\cal{F}}_{\rm
MZM}=\frac{1}{{\mathsf{N}}_{0}}\sum_{j=1}^{L}\Theta_{j}^{a}\left(u^{1}_{j}+v^{1}_{j}\right).$
(63)

##  Appendix B Surface magnetization

###  B.1 Definition

The surface magnetization Peschel (1984); Karevski (2000) is defined by

${{m}}_{1,L}^{s}=\langle\,n_{-}\,|\sigma_{1,L}^{x}|\,n_{+}\,\rangle.$ (64)

The boundary magnetization can be expressed as follows

$\displaystyle\sigma_{1}^{x}$ $\displaystyle=$ $\displaystyle
a_{1}=c_{1}^{\dagger}+c_{1}^{\vphantom{\dagger}}$ (65) $\displaystyle=$
$\displaystyle\sum_{m=1}^{L}\left(u^{m}_{1}+v^{m}_{1}\right)\left(\phi^{\dagger}_{m}+\phi^{\vphantom{\dagger}}_{m}\right),$

and

$\displaystyle\sigma_{L}^{x}$ $\displaystyle=$
$\displaystyle{\rm{i}}{\mathbb{P}}b_{L}=-{\mathbb{P}}\left(c_{L}^{\dagger}-c_{L}^{\vphantom{\dagger}}\right)$
(66) $\displaystyle=$
$\displaystyle-{\mathbb{P}}\sum_{m=1}^{L}\left(u^{m}_{L}-v^{m}_{L}\right)\left(\phi^{\dagger}_{m}-\phi^{\vphantom{\dagger}}_{m}\right).$

###  B.2 Invariance across the many-body spectrum

Since two partner states ${|\,n_{p}\,\rangle}$ and ${|\,n_{-p}\,\rangle}$ only
differ by their occupation of the lowest single-particle mode $m=1$, one gets

$m_{1}^{s}=u^{1}_{1}+v^{1}_{1},$ (67)

which we take positive by convention. On the other hand, at the right boundary
we have

$m_{L}^{s}=\left(u^{1}_{L}-v^{1}_{L}\right)\langle\,n_{-}\,|\left(\phi^{\dagger}_{m}-\phi^{\vphantom{\dagger}}_{m}\right)|\,n_{+}\,\rangle,$
(68)

whose sign depends on the state. Indeed, if the even ($p=+1$) many-body
eigenstate $|\,n_{+}\,\rangle$ has the lowest single-particle mode $m=1$ which
is empty (like the GS for instance), we have $m_{L}^{s}=u^{1}_{L}-v^{1}_{L}$
($=m_{1}^{s}$ if the system is clean). Conversely if the mode $m=1$ is
occupied, one gets $m_{L}^{s}=v^{1}_{L}-u^{1}_{L}$ ($=-m_{1}^{s}$ again if the
system is clean). We therefore see that the magnitude of the surface
magnetization does not depend on the many-body energy $E_{n}^{\pm}$, but only
on the coeffiient of the lowest ($m=1$) single particle mode. Only its
relative sign between left and right boundaries oscillates across the
spectrum, see also Tab. 1.

###  B.3 Critical scaling of the surface magnetization

In the ordered regime, for $L\gg\xi_{\rm zm}$, one can use the mapping
${|\,n\,\rangle}_{+}=\psi_{a}{|\,n\,\rangle}_{-}$, that brings us to the
following relation

$\displaystyle m_{1}^{s}$ $\displaystyle=$
$\displaystyle\langle\,n_{-}\,|\sigma_{1}^{x}\psi_{a}|\,n_{-}\,\rangle=\frac{1}{N_{a}},$
(69)

where the MZM norm results from a geometrical sum of exponentials, see Eq. (5)
and Eq. (6), simply yielding

$\frac{1}{N_{a}}=\sqrt{1-\exp(-2/\xi_{\rm zm})}.$ (70)

When approaching the critical point from the ordered regime, the length scale
$\xi_{\rm zm}$ becomes very large, and therefore the surface magnetization is
given by

$m_{1}^{s}\approx\sqrt{2/\xi_{\rm zm}}.$ (71)

Figure 7: Critical scaling of the surface magnetization for the TFI chain at
$\Gamma=h/X=1$. ED data are compared with the expected algebraic behavior
$\sqrt{2/(L+\ell_{\rm eff})}$.

When $\xi_{\rm zm}\gg L$, it is rather standard to replace this length scale
by the lattice size $L$ in Eq. (71). In the main text, we have argued that the
correct length scale to use is the one which also enters in the wave vector
Eq. (22) $k_{\rm min}$, i.e. $L+\ell_{\rm eff}=L+1/2$.

Fig. 7 provides the numerical justification for this: ED data for
$m_{1}^{s}(L)$ are shown against $1/L$ up to $L=10^{4}$ for the transverse
field Ising chain at criticality. A comparison with the analytical form
$\sqrt{\frac{2}{L+\ell_{\rm eff}}}$ is shown for $\ell_{\rm eff}=0$ (open
black circles) and $\ell_{\rm eff}=1/2$ (open red squares). It is clearly
visible that using a finite shift $\ell_{\rm eff}=1/2$ gives a much better
description (see inset).

##  Appendix C MZM in the general case of the XY-Kitaev chain

###  C.1 Analytical expression in the ordered regime

####  C.1.1 Iterative procedure

Using Majoranas, the clean XY-Kitaev model reads

${\cal{H}}_{\rm
XY-K}={\rm{i}}\sum_{j}\left(Xb_{j}a_{j+1}-Ya_{j}b_{j+1}+ha_{j}b_{j}\right).$
(72)

Assuming simple linear combinations for the MZM operators

$\Psi_{a}=\frac{1}{{{N}}_{a}}\sum_{j=1}^{L}\Theta_{j}\,a_{j}\quad{\rm
and}\quad\Psi_{b}=\frac{1}{{{N}}_{b}}\sum_{j=1}^{L}\Theta_{j}\,b_{L+1-j},$
(73)

and using $\left[{\cal{H}}_{\rm XY-K},a_{j}\right]=2{\rm
i}\left(Xb_{j-1}-hb_{j}+Yb_{j+1}\right)$ and $\left[{\cal{H}}_{\rm
XY-K},b_{j}\right]=-2{\rm i}\left(Ya_{j-1}-ha_{j}+Xa_{j+1}\right)$, we
iteratively arrive at the simple recursion relation for
$\Theta_{j}=\Theta_{j}^{a,b}$:

$\Theta_{j+1}=\frac{h}{X}\Theta_{j}-\frac{Y}{X}\Theta_{j-1},$ (74)

such that

$\displaystyle\left[{\cal{H}}_{\rm XY-K},\Psi_{a}\right]$ $\displaystyle=$
$\displaystyle
2{\rm{i}}\frac{1}{{{N}}_{a}}\left(Y\Theta_{L-1}-h\Theta_{L}\right)b_{L}$ (75)
$\displaystyle\left[{\cal{H}}_{\rm XY-K},\Psi_{b}\right]$ $\displaystyle=$
$\displaystyle-2{\rm{i}}\frac{1}{{{N}}_{b}}\left(Y\Theta_{L-1}-h\Theta_{L}\right)a_{1}.$
(76)

One solve the recursion Eq. (74) with initial conditions $\Theta_{0}=0$ and
$\Theta_{1}=1$, restricting to positive couplings $h,\,X,\,Y\geq 0$, and
$X\geq Y$ (other cases can be easily derived).

####  C.1.2 MZM

We note in passing that the phase diagram of the XY-Kitaev chain model, shown
in Fig. 1 (a), can simply be inferred from the existence of normalizable MZMs,
which requires that the largest eigenvalue of the Eq. (74) has its modulus
less that one, i.e. if $X+Y>h$. Contrary to the TFIM case, here the
topological regime is richer as one can distinguish two types of MZM decays,
incommensurate and commensurate.

(i) Incommensurate regime ($h^{2}<4XY$).

In this case,

$\Theta_{j}=\frac{2X}{\sqrt{4XY-h^{2}}}\,\sin\left(\varphi
j\right)\,{\rm{e}}^{-j/\xi_{\rm zm}}$ (77)

displays oscillations and exponential decay, controlled by

$\cos\varphi=\frac{h}{2\sqrt{XY}}\quad{\rm{and}}\quad\frac{1}{\xi_{\rm
zm}}={\ln\sqrt{\frac{X}{Y}}}.$ (78)

The MZM normalization factor $N_{a}=N_{b}\equiv{\mathsf{N}}_{0}$ can be
evaluated in the large $L$ limit

$\frac{1}{{{\mathsf{N}}}_{0}}\xrightarrow[L\to\infty]{}\sqrt{2}\sin\varphi\left[\frac{1}{1-\frac{Y}{X}}+\frac{\frac{Y}{X}-\cos\left(2\varphi\right)}{1-2\frac{Y}{X}\cos\left(2\varphi\right)+\left(\frac{Y}{X}\right)^{2}}\right]^{-\frac{1}{2}}.$
(79)

(ii) Commensurate regime ($h^{2}\geq 4XY$). In this case

$\displaystyle\Theta_{j}$ $\displaystyle=$ $\displaystyle\frac{X}{\alpha
h}\left(\frac{h}{2X}\right)^{j}\times\Bigl{[}\left(1+\alpha\right)^{j}-\left(1-\alpha\right)^{j}\Bigr{]}$
(80) $\displaystyle\xrightarrow[j\gg 1]{}\begin{cases}\frac{X}{\alpha
h}{\rm{e}}^{-j/\xi_{\rm zm}}&\mbox{if }2\sqrt{XY}>h>X+Y\\\ \infty&\mbox{if
}h>X+Y,\end{cases}$

where the edge mode localization length $\xi_{\rm zm}$ is given by

$\frac{1}{\xi_{\rm zm}}=\ln\left[\frac{2X}{(1+\alpha)h}\right],$ (81)

and $\alpha=\sqrt{1-4XY/h^{2}}$. The MZM normalization factor can be expressed
in the large chain $L$ limit:

$\displaystyle\frac{1}{{{\mathsf{N}}}_{0}}$
$\displaystyle\xrightarrow[L\to\infty]{}$ $\displaystyle
2\alpha\Bigl{[}\frac{1}{(1+\alpha)^{-2}-\left(\frac{h}{2X}\right)^{2}}+\frac{1}{(1-\alpha)^{-2}-\left(\frac{h}{2X}\right)^{2}}$
(82) $\displaystyle-$
$\displaystyle\frac{2}{(1-\alpha^{2})^{-1}-\left(\frac{h}{2X}\right)^{2}}\Bigr{]}^{-\frac{1}{2}},\quad{\rm
if}~{}\alpha\neq 0.$

When $Y=0$ ($\alpha=0$), we simply recover the TFIM result Eq. (16), i.e.
$1/{\mathsf{N}}_{0}\xrightarrow[L\to\infty]{}\sqrt{1-h^{2}/X^{2}}$.

###  C.2 Critical behavior

(i) Critical commutators. At criticality when $h=X+Y$, the MZM coefficients do
not decay anymore and are given by

$\Theta_{j}^{\rm
critical}=\frac{X}{X-Y}\left[1-\left(\frac{Y}{X}\right)^{j}\right].$ (83)

Therefore, the critical commutators are easy to compute: Eq. (76) thus becomes

$\displaystyle\left[{\cal{H}}_{\rm XY-K},\Psi_{a}\right]$ $\displaystyle=$
$\displaystyle-2{\rm{i}}\frac{X^{2}}{{{N}}_{a}(X-Y)}\left[1-\left(\frac{Y}{X}\right)^{L+1}\right]b_{L}$
$\displaystyle\left[{\cal{H}}_{\rm XY-K},\Psi_{b}\right]$ $\displaystyle=$
$\displaystyle
2{\rm{i}}\frac{X^{2}}{{{N}}_{b}(X-Y)}\left[1-\left(\frac{Y}{X}\right)^{L+1}\right]a_{1},$
(84)

which indeed rapidly converges ($Y<X$) at large enough $L$ to the results Eq.
(24) that we rewrite here

$\displaystyle\left[{\cal{H}}_{\rm XY-K},\psi_{a}\right]$ $\displaystyle=$
$\displaystyle\frac{-2{\rm{i}}X^{2}}{{{N}}_{a}(X-Y)}b_{L}$ (85)
$\displaystyle\left[{\cal{H}}_{\rm XY-K},\psi_{b}\right]$ $\displaystyle=$
$\displaystyle\frac{2{\rm{i}}X^{2}}{{{N}}_{b}(X-Y)}a_{1}.$ (86)

(ii) Critical MZM norm. Building on the critical coefficients Eq. (83), one
easily gets the norm (using $y=Y/X$)

$\displaystyle N_{a,b}={\mathsf{N}}_{0}$ $\displaystyle=$
$\displaystyle\sqrt{\sum_{j=1}^{L}\left(\Theta_{j}^{\rm critical}\right)^{2}}$
(87) $\displaystyle=$
$\displaystyle\frac{X}{X-Y}\sqrt{L-\frac{y(2+y)}{1-y^{2}}}$

thus yielding Eq. (25) with the effective length shift $\ell^{\prime}_{\rm
eff}$ given by Eq. (26).

## References

  * Lutchyn  _et al._ (2010) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Sarma  _et al._ (2015) S. D. Sarma, M. Freedman, and C. Nayak, Quantum Inf 1, 15001 (2015).
  * Marra (2022) P. Marra, Journal of Applied Physics 132, 231101 (2022).
  * Mi (2022) X. et al. Mi, Science 378, 785 (2022).
  * Dvir (2023) T. et al. Dvir, Nature 614, 445 (2023).
  * Yazdani  _et al._ (2023) A. Yazdani, F. von Oppen, B. I. Halperin, and A. Yacoby, Science 380, eade0850 (2023).
  * Lieb  _et al._ (1961) E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Pfeuty (1970) P. Pfeuty, Annals of Physics 57, 79 (1970).
  * Kitaev (2001) A. Y. Kitaev, Phys.-Usp. 44, 131 (2001).
  * Fendley (2012) P. Fendley,  2012, P11020 (2012).
  * Fendley (2016) P. Fendley, J. Phys. A: Math. Theor. 49, 30LT01 (2016).
  * Kemp  _et al._ (2017) J. Kemp, N. Y. Yao, C. R. Laumann, and P. Fendley, Journal of Statistical Mechanics: Theory and Experiment 2017, 063105 (2017).
  * Else  _et al._ (2017) D. V. Else, P. Fendley, J. Kemp, and C. Nayak, Phys. Rev. X 7 (2017).
  * Olund  _et al._ (2023) C. T. Olund, N. Y. Yao, and J. Kemp, “Boundary Strong Zero Modes,”  (2023), arXiv:2305.16382.
  * Alexandradinata  _et al._ (2016) A. Alexandradinata, N. Regnault, C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. B 94, 125103 (2016).
  * Wada  _et al._ (2021) K. Wada, T. Sugimoto, and T. Tohyama, Phys. Rev. B 104 (2021).
  * Chepiga and Laflorencie (2023a) N. Chepiga and N. Laflorencie, SciPost Physics 14, 152 (2023a).
  * not (a) This becomes evident when expressing the parity operator as ${\mathbb{P}}=\prod_{j}{\rm{i}}a_{j}b_{j}$.
  * Catalan (1865) E. Catalan, Memoir on the transformation of series and some definite integrals 33 (1865).
  * not (b) If ${\cal S}=+1$ (ferromagnetic), among the two partners $p=\pm 1$ the lowest-energy is in the even sector $p=+1$, while for ${\cal S}=-1$ (antiferromagnetic), the odd state $p=-1$ has the lowest energy.
  * Suzuki  _et al._ (2013) S. Suzuki, J.-i. Inoue, and B. K. Chakrabarti,  _Quantum Ising Phases and Transitions in Transverse Ising Models_, Lecture Notes in Physics, Vol. 862 (Springer, Berlin, Heidelberg, 2013).
  * Peschel (1984) I. Peschel, Phys. Rev. B 30, 6783 (1984).
  * Karevski (2000) D. Karevski, J. Phys. A: Math. Gen. 33, L313 (2000).
  * Cabrera and Jullien (1987) G. G. Cabrera and R. Jullien, Phys. Rev. B 35, 7062 (1987).
  * Campostrini  _et al._ (2015) M. Campostrini, A. Pelissetto, and E. Vicari, J. Stat. Mech. 2015, P11015 (2015).
  * Campostrini  _et al._ (2014) M. Campostrini, A. Pelissetto, and E. Vicari, Phys. Rev. B 89, 094516 (2014).
  * Fisher (1994) D. S. Fisher, Phys. Rev. B 50, 3799 (1994).
  * Iglói and Juhász (2008) F. Iglói and R. Juhász, EPL 81, 57003 (2008).
  * Kawabata  _et al._ (2019) K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Nat Commun 10, 297 (2019).
  * Bergholtz  _et al._ (2021) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Motrunich  _et al._ (2001) O. Motrunich, K. Damle, and D. A. Huse, Phys. Rev. B 63, 224204 (2001).
  * Bravyi and König (2012) S. Bravyi and R. König, Commun. Math. Phys. 316, 641 (2012).
  * DeGottardi  _et al._ (2013) W. DeGottardi, D. Sen, and S. Vishveshwara, Phys. Rev. Lett. 110, 146404 (2013).
  * Hegde and Vishveshwara (2016) S. S. Hegde and S. Vishveshwara, Phys. Rev. B 94, 115166 (2016).
  * Gergs  _et al._ (2016) N. M. Gergs, L. Fritz, and D. Schuricht, Phys. Rev. B 93, 075129 (2016).
  * Monthus (2018a) C. Monthus, J. Phys. A: Math. Theor. 51, 465301 (2018a), publisher: IOP Publishing.
  * Lieu  _et al._ (2018) S. Lieu, D. K. K. Lee, and J. Knolle, Phys. Rev. B 98, 134507 (2018).
  * Levy and Goldstein (2019) L. Levy and M. Goldstein, Universe 5, 33 (2019).
  * Francica  _et al._ (2023) G. Francica, E. M. Tiburzi, and L. Dell’Anna, Phys. Rev. B 107, 165137 (2023).
  * Fisher (1992) D. S. Fisher, Phys. Rev. Lett. 69, 534 (1992).
  * Fisher (1995) D. S. Fisher, Phys. Rev. B 51, 6411 (1995).
  * Young and Rieger (1996) A. P. Young and H. Rieger, Phys. Rev. B 53, 8486 (1996).
  * Fisher and Young (1998) D. S. Fisher and A. P. Young, Phys. Rev. B 58, 9131 (1998).
  * Iglói and Rieger (1998) F. Iglói and H. Rieger, Phys. Rev. B 57, 11404 (1998).
  * Refael and Moore (2004) G. Refael and J. E. Moore, Phys. Rev. Lett. 93, 260602 (2004).
  * Iglói and Monthus (2005) F. Iglói and C. Monthus, Phys. Rep. 412, 277 (2005).
  * Bonesteel and Yang (2007) N. E. Bonesteel and K. Yang, Phys. Rev. Lett. 99, 140405 (2007).
  * Fidkowski  _et al._ (2008) L. Fidkowski, G. Refael, N. E. Bonesteel, and J. E. Moore, Phys. Rev. B 78, 224204 (2008).
  * Vojta  _et al._ (2009) T. Vojta, C. Kotabage, and J. A. Hoyos, Phys. Rev. B 79 (2009), 10.1103/PhysRevB.79.024401.
  * Laflorencie (2022) N. Laflorencie, in  _Entanglement in Spin Chains: From Theory to Quantum Technology Applications_, Quantum Science and Technology, edited by A. Bayat, S. Bose, and H. Johannesson (Springer International Publishing, Cham, 2022) pp. 61–87.
  * Zabalo  _et al._ (2023) A. Zabalo, J. H. Wilson, M. J. Gullans, R. Vasseur, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, Phys. Rev. B 107, L220204 (2023).
  * Fidkowski and Kitaev (2010) L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010).
  * Fidkowski and Kitaev (2011) L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011).
  * Hassler and Schuricht (2012) F. Hassler and D. Schuricht, New J. Phys. 14, 125018 (2012).
  * Thomale  _et al._ (2013) R. Thomale, S. Rachel, and P. Schmitteckert, Phys. Rev. B 88, 161103 (2013).
  * Katsura  _et al._ (2015) H. Katsura, D. Schuricht, and M. Takahashi, Phys. Rev. B 92, 115137 (2015).
  * Goldstein and Chamon (2012) G. Goldstein and C. Chamon, Phys. Rev. B 86, 115122 (2012).
  * Kells (2015) G. Kells, Phys. Rev. B 92, 155434 (2015).
  * Więckowski  _et al._ (2018) A. Więckowski, M. M. Maśka, and M. Mierzejewski, Phys. Rev. Lett. 120, 040504 (2018).
  * Yates  _et al._ (2020a) D. J. Yates, A. G. Abanov, and A. Mitra, Phys. Rev. Lett. 124, 206803 (2020a).
  * Kemp (2019) J. Kemp,  _Strong edge zero modes in interacting systems_ , Ph.D. thesis, University of Oxford (2019).
  * Więckowski and Ptok (2019) A. Więckowski and A. Ptok, Phys. Rev. B 100, 144510 (2019).
  * Kemp  _et al._ (2020) J. Kemp, N. Y. Yao, and C. R. Laumann, Phys. Rev. Lett. 125, 200506 (2020).
  * Pandey  _et al._ (2023) B. Pandey, N. Mohanta, and E. Dagotto, Phys. Rev. B 107, L060304 (2023).
  * Yates  _et al._ (2020b) D. J. Yates, A. G. Abanov, and A. Mitra, Phys. Rev. B 102, 195419 (2020b).
  * Rahmani  _et al._ (2015) A. Rahmani, X. Zhu, M. Franz, and I. Affleck, Phys. Rev. B 92, 235123 (2015).
  * Alet and Laflorencie (2018) F. Alet and N. Laflorencie, Comptes Rendus Physique Quantum simulation / Simulation quantique, 19, 498 (2018).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Chepiga and Laflorencie (2023b) N. Chepiga and N. Laflorencie, “Resilient infinite randomness criticality for a disordered chain of interacting Majorana fermions,”  (2023b), arXiv:2305.10363.
  * Monthus (2018b) C. Monthus, J. Phys. A: Math. Theor. 51, 115304 (2018b).
  * Roberts and Motrunich (2021) B. Roberts and O. I. Motrunich, Phys. Rev. B 104, 214208 (2021).
  * Slagle  _et al._ (2016) K. Slagle, Y.-Z. You, and C. Xu, Phys. Rev. B 94, 014205 (2016).
  * You  _et al._ (2016) Y.-Z. You, X.-L. Qi, and C. Xu, Phys. Rev. B 93, 104205 (2016).
  * Kells  _et al._ (2018) G. Kells, N. Moran, and D. Meidan, Phys. Rev. B 97, 085425 (2018).
  * Sahay  _et al._ (2021) R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y. Yao, Phys. Rev. Lett. 126, 100604 (2021).
  * Wahl  _et al._ (2022) T. B. Wahl, F. Venn, and B. Béri, Phys. Rev. B 105, 144205 (2022).
  * Laflorencie  _et al._ (2022) N. Laflorencie, G. Lemarié, and N. Macé, Phys. Rev. Research 4, L032016 (2022).
  * Hannukainen  _et al._ (2023) J. D. Hannukainen, M. F. Martínez, J. H. Bardarson, and T. K. Kvorning, “Interacting Local Topological Markers: A one-particle density matrix approach for characterizing the topology of interacting and disordered states,”  (2023), arXiv:2307.06447.
  * Moudgalya  _et al._ (2020) S. Moudgalya, D. A. Huse, and V. Khemani, arXiv:2008.09113 (2020).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Yates  _et al._ (2018) D. Yates, Y. Lemonik, and A. Mitra, Phys. Rev. Lett. 121, 076802 (2018).
  * Bauer  _et al._ (2019) B. Bauer, T. Pereg-Barnea, T. Karzig, M.-T. Rieder, G. Refael, E. Berg, and Y. Oreg, Phys. Rev. B 100, 041102 (2019).
  * Yates and Mitra (2021) D. J. Yates and A. Mitra, Phys. Rev. B 104, 195121 (2021).
  * Harle  _et al._ (2023) N. Harle, O. Shtanko, and R. Movassagh, Nat Commun 14, 2286 (2023).
  * Ling  _et al._ (2023) H. Ling, P. Richard, S. R. Koshkaki, M. Kolodrubetz, D. Meidan, A. Mitra, and T. Pereg-Barnea, “Disorder induced topological phase transition in a driven Majorana chain,”  (2023), arXiv:2310.17088.