Discussion and Comparison with the Literature

During these simulations, energy values appearing in pairs of opposite signs was used as a program test. That is, the resulting energy values are supposed to respect the particle-hole symmetry. Numerically, however, we observed that when N was increased, or µ chosen sufficiently close to zero, the resulting ground state energy become of equal magnitude as the floating point precision of the computer.³ Conse-quently, the lowest energy pair did not appear exactly equal in absolute value. This numerical violation the particle-hole symmetry also resulted in eigenstates rotated within the ground state manifold, preferably giving weights on one or the other side of the chain. This problem boils down to how the Armadillo function eig_sym() specifically finds eigenvectors, and it could possibly have been avoided by developing a code from scratch. Instead, we kept the system sizes small and used the appear-ance of paired eigenvalues as an indication of whether the correct eigenvectors were returned or not. Exactly this problematic effect was present when fixing ∆ =t and µ = 0. Instead of observing the weights of d0 returned, the numerics yielded the Majorana modes−c1+c^†₁ orcN+c^†_N with weights on either side of the system. The system size and parameters in the figures were chosen to avoid this problem, which is the reason for usingN = 50 and not N = 80 in Figure 4.3.

The results with constant order parameter demonstrate that there is a non-trivial and important consequence of having an open Kitaev chain in the topological phase,

|µ| < t. Namely, that one Majorana zero mode emerge on each of the open ends.

The modes are the constituents of a fermionic operator associated with (exponen-tially close to) zero energy. When the order parameter was varied spa(exponen-tially, ∆ = ∆x, we observe in Figure 4.5(a) that two new Majorana modes appear when ∆ changes

3Typically, a number withdoublefloating-point precision inC++have 15-17 significant decimal digits [27].

Section 4.4 Numerical Implementation, Results and Discussion 65

0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 Chemical potential, ^µ/t

−0.1 0.0 0.1 0.2 0.3 0.4

Energ eigenvalue, Ei

Energ spectrum with kinked order parameter,

^{∆ =∆}x

−2.0−1.5−1.0−0.5 0.0 0.5 1.0 1.5 2.0

µ/t

−1.5−1.0

−0.50.00.51.01.5

Ei

(a)

0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 Chemical potential, ^µ/t

−0.1 0.0 0.1 0.2 0.3 0.4

Energ eigenvalue, Ei

Energ spectrum with kinked order parameter,

^{∆ = |∆}x|

−2.0−1.5−1.0−0.5 0.0 0.5 1.0 1.5 2.0

µ/t

−1.5−1.0

−0.50.00.51.01.5

Ei

(b)

Figure 4.5: Energy spectrum of open Kitaev chain with (a) ∆ = ∆_x and (b)

∆ = |∆x|. The order parameter is given by equation (4.44) with s = 10 and

∆₀ = 0.5. The main plots show the zoomed energy spectra close to the phase transition µ = t. The insets show the spectra plotted for different values of chemical potential in the range µ∈[−2t,2t] with the zoom marked in red. The number of sites was set toN = 80. The hopping parameter was fixed to t= 1.0.

sign. This basic observation was seen to be robust in s. Specifically, by increasing s significantly – making the zero-crossing of ∆x more spread out in space – the only observed change was that the central peaks in Figure 4.6 were widened. In crucial contrast, new zero modes did never seem to appear with ∆x = |∆x|, still independent ofs. Moreover, the central zero modes appearing with ∆ = ∆x have a significant overlap but they nonetheless coexist, i.e. they do not split in energy and

0.0 0.2 0.4 0.6 0.8 1.0 Position mapped to unit interval, ^xj/N

−0.4

−0.3

−0.2

−0.1 0.0 0.1 0.2 0.3 0.4

Fermion weight coefficients,q

Fermion content in lowest energ state with kinked order parameter.

^Ei=1.4e−09

,

^El=2.9e−09

q_2i−1,2j q_{2i−1,2j−1} q_2l−1,2j q_{2l−1,2j−1}

Figure 4.6: Fermion content of the two lowest energy fermionic states in the open Kitaev chain with position varying order parameter. The matrix coefficients q, as in equation (4.42), are plotted against x_j/N. This plot was produced with N = 80 sites, hopping parametert= 1.0, chemical potentialµ= 0.3tand ∆ = ∆_x as in equation (4.44) with s= 10 and ∆0 = 0.5. The energies Ei andEl refer to the numerically found (positive) eigenvalues.

hybridize. This leads us to think that the sign change of ∆ has roots in topological properties of the Kitaev chain.

While working with this problem we became aware of a recent article where sim-ilar effects are studied. In [28] T. H. Hansson et al. discuss zero modes resulting from π-junctions in ∆ on more general grounds. In their treatment, phase winding junctions, where the phase of the order parameter winds by a position dependent angle, are considered. They argue that a real order parameter going through zero always gives rise to two additional zero modes in the topological phase for p-wave chains. This happens regardless of the details on how the order parameter passes zero. However, when the phase winding becomes genuinely complex, i.e. not rotat-ing from 0 to π from one site to the next, no additional zero modes than those on the ends of the chain appear.

The difference between complex and real junctions, they argue, is more deeply rooted in symmetry classes. The authors observe that the Kitaev chain with real order parameter and real hopping strength belongs to the BDI symmetry class (TRS is +1), while it belongs to class D (TRS is 0) when they are both complex. A BDI Kitaev chain has two distinct realizations of the topological phase. A key observation supporting this view can actually be made from Figure 3.3. Withφ= 0 we see that the transformation ∆→ −∆ alters the circulation ofh(k) from clockwise to counter-clockwise in the yz-plane. The two values ±∆ correspond to opposite winding numbers, ν. In other words, ±∆ are distinct topological phases, and two

Section 4.4 Numerical Implementation, Results and Discussion 67

zero energy Majorana modes appear on the boundary between the phases. The simulations shown in Figure 4.5 and 4.6 are in agreement with and underpins this conclusion. Finally, that±∆ are distinct realizations of the topological phase might explain the fact that we observe two coexisting zero modes with significant spatial overlap in Figure 4.6. Overlapping zero modes are normally expected to split in a finite ±E pair. With these remarks we conclude our discussion of the Kitaev chain.

Chapter 5

The p + ip Model and Vortices with Majorana Modes

A two-dimensional model with p-wave paired and effectively spinless electrons is studied. First, we establish the formal framework for the model. This includes de-riving the diagonalization equations. Thereafter, the idea of a space varying order parameter will be pursued. In particular, we search for Majorana zero mode solu-tions in an infinite system with a rotationally symmetric vortex. In 2007 V. Guriare and L. Radzihovsky studied this system analytically in the limit where |∆(r)| is non-constant in an infinitesimally small region [13]. We depart from their study by using a vortex solution in accordance with Ginzburg-Landau theory. Moreover, we propose an argument that implies non-Abelian exchange statistics in a system of several vortices. This property is known from Ivanov’s compact consideration [6]. However, our argument is, in contrast to Ivanov’s derivation, based on con-servation of fermion parity. The p+ip mean field model, which is valid with a non-homogeneous order parameter (see for instance [6, 12]), is given by

H = Z

d²r ψ^†(r)

− ¯h²

2m∇²−µ

ψ(r) +1 2

Z

d²rd²r⁰ ψ^†(r)D(r,r⁰)ψ^†(r⁰) + h.c.

. (5.1) Above, ψ^†(r) is the creation operator of a spinless fermion in position r. The operators satisfy fermionic anticommutation relations

{ψ(r), ψ^†(r⁰)}=δ⁽²⁾(r−r⁰) and {ψ(r), ψ(r⁰)}={ψ^†(r), ψ^†(r⁰)}= 0, (5.2) with δ⁽²⁾(r−r⁰) being the two-dimensional Dirac delta function. The pairing func-tion D(r,r⁰) is defined by

D(r,r⁰) = ∆

r+r⁰ 2

(∂x⁰ +i∂y⁰)δ⁽²⁾(r−r⁰). (5.3) 69

We open the chapter with a brief presentation of the model in an in infinite system with constant order parameter to point out the analogy to the closed Kitaev chain and to motivate the topological classification.

5.1 Homogeneous System

We summarize some essential properties of the model as discussed in [11, 24] for completeness. Assume that the order parameter is homogeneous, ∆(r) = ∆0e^iφ, with φ some general phase and ∆0 a real, positive constant. The Hamiltonian in (5.1) reduce to¹

H = Z

d²r

ψ^†(r)

−¯h²

2m∇²−µ

ψ(r) + ∆0

2 e^iφψ^†(r)(∂x+i∂y)ψ^†(r) + h.c.

. (5.4) By introducing a momentum representation of the fermionic operators,

ψ(k) = Z

d²r e^−ik·rψ(r) and ψ(r) =

Z d²k

(2π)²e^ik·rψ(k), (5.5) one may diagonalize the Hamiltonian. Inserting the transformed operators in the Hamiltonian (5.4) and using R

d²r e^ir·(k−k0) = (2π)²δ⁽²⁾(k−k⁰), results in the ex-pression (compare with equation (3.12)),

H = 1 2

Z d²k

(2π)²Ψ^†(k)H^{(k)Ψ(k), (5.6)} with

Ψ(k)≡

ψ(k) ψ^†(−k)

and H^(k)≡

(k) ∆(k)

∆^∗(k) −(k)

. (5.7)

Moreover, we introduced the functions

(k)≡ ¯h²k²

2m −µ and ∆(k)≡i∆0e^iφ(kx+iky). (5.8) The eigenvalues of H(k) define the quasiparticle spectrum and are given by

±E(k) = ±p

(k)²+|∆(k)|² =±

s¯h²k² 2m −µ

²

+ ∆²₀k². (5.9)

1This follows by making use of the Dirac delta function propertyR

dx f(x)∂x⁰δ(x−x⁰) =^df(x_dx0⁰⁾.

In document Topological Superconductivity and Majorana Fermions (sider 74-81)