In [2]:
%run ../code/init_mooc_nb.ipy


# additional imports
from matplotlib import gridspec
import matplotlib

# global variables
default_color_cycle = matplotlib.rcParams['axes.color_cycle']
sigma0 = np.array([[1., 0.], [0., 1.]])
sigmax = np.array([[0., 1.], [1., 0.]])
sigmay = np.array([[0., -1j], [1j, 0.]])
sigmaz = np.array([[1., 0.], [0., -1.]])


# simulation related functions
def spinful_kitaev_chain(L=None, t=1, delta=0.1, mu=0.1, B=0.1):
    """ 
    Return finalized kwant system for a kitaev chain. 
    
    For L=None function returns infinite (TranslationalSymmetry)
    system. If L is set to positive integer the system is finite.
    """
    lat = kwant.lattice.chain()
    
    if L is None:
        sys = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
        L = 1
    else:
        sys = kwant.Builder()
        
    # Calculate Kronecker products.
    szs0 = np.kron(sigmaz, sigma0)
    szsz = np.kron(sigmaz, sigmaz)
    sys0 = np.kron(sigmay, sigma0)
    
    # onsite terms
    for x in xrange(L):
        sys[lat(x)] = -(mu - 2*t) * szs0 + B * szsz
         
    # hopping terms
    sys[kwant.HoppingKind((1,), lat)] = - t * szs0 - 1j * delta * sys0
    
    return sys.finalized()


def nanowire_chain(L=None, t=1, mu=0.1, delta=0.1, B=0.1, alpha=0):
    """ 
    Return finalized kwant system for a nanowire chain. 
    
    For L=None function returns infinite (TranslationalSymmetry)
    system. If L is set to positive integer the system is finite.
    """
    lat = kwant.lattice.chain()
    
    if L is None:
        sys = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
        L = 1
    else:
        sys = kwant.Builder()

    # Calculate Kronecker products.
    szs0 = np.kron(sigmaz, sigma0)
    s0sz = np.kron(sigma0, sigmaz)
    sxs0 = np.kron(sigmax, sigma0)
    szsy = np.kron(sigmaz, sigmay)
    
    # onsite terms
    for x in xrange(L):
        sys[lat(x)] = (2*t - mu) * szs0 + B * s0sz + delta * sxs0
         
    # hopping terms
    sys[kwant.HoppingKind((1,), lat)] = -t * szs0 - 0.5 * 1j * alpha * szsy
    
    return sys.finalized()

    
def spinorbit_band(t, mu, delta, B, alpha, steps=101):
    """ Calculate band structure of nanowire chain with spin-orbit. """
    
    ks=np.linspace(-np.pi/3, np.pi/3, steps)
    sys = nanowire_chain(None, t, mu, delta, B, alpha)

    fig, ax = plt.subplots()
    ax.set_color_cycle(['k'])
    kwant.plotter.bands(sys, momenta=ks, ax=ax)

    ax.set_xlabel('$k$')
    ax.set_xticks([-1.0, 0.0, 1.0])
    ax.set_ylabel('$E$')
    ax.set_xlim([-1.03, 1.03])
    ax.set_ylim([-1.0, 1.0])
    ax.set_yticks([-1.0, 0.0, 1.0])
    fig.suptitle('$\mu = %1.2f $, $B = %1.2f$, ' r' $\alpha = %1.2f$, $\Delta = %1.3f$' % (mu, B, alpha, delta))
        
    return fig

def find_gap(sys, resolution=1e-4):
    """Find gap in a system by doing a binary search in energy."""
    # This tells us if there are modes at a certain energy.
    if len(sys.modes(energy=0)[0].momenta):
        return 0
    gap = step = min(abs(kwant.physics.Bands(sys)(k=0))) / 2
    while step > resolution:
        step /= 2
        if len(sys.modes(gap)[0].momenta):
            gap -= step
        else:
            gap += step

    return gap
    

def spinorbit_band_gap(mu, t, delta, B_vals):
    gaps=[]
    alphas = np.linspace(0, 0.3, 4)
    # very slow, but understandable
    for alpha in alphas:
        temp_gaps = []
        for B in B_vals:
            sys = nanowire_chain(t=t, mu=mu, delta=delta, alpha=alpha, B=B)
            temp_gaps.append(find_gap(sys))
        gaps.append(np.array(temp_gaps))
    
    fig = plt.figure()
    for gap, alpha in zip(gaps, alphas):
        plt.plot(B_vals, gap, label= (r'$\alpha = %1.1f $' % alpha))
    plt.xlabel('$B$')
    plt.ylabel('Band gap')
    plt.xticks([0, 0.1, 0.2, 0.3])
    ylim = [0.0, 0.14]
    plt.ylim(ylim)
    plt.yticks([0, 0.05, 0.1])
    plt.legend(bbox_to_anchor=(0.0, 0.75, 1.0, .09), loc=3,
               ncol=2, mode="expand", borderaxespad=0., prop={'size':16})
    #plt.subplots_adjust(top = 0.2)
    plt.title(r'$\mu = %1.3f$, $\Delta = %1.1f$' % (mu, delta))
    
    # Add a vertical line to show where the topological region begins:
    B_crit = np.sqrt(mu**2 + delta**2)
    plt.plot([B_crit, B_crit], ylim, 'r--')
    
    return fig


ERROR: File `u'../code/init_mooc_nb.ipy.py'` not found.
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-bd7a1320a9e8> in <module>()
      8 # global variables
      9 default_color_cycle = matplotlib.rcParams['axes.color_cycle']
---> 10 sigma0 = np.array([[1., 0.], [0., 1.]])
     11 sigmax = np.array([[0., 1.], [1., 0.]])
     12 sigmay = np.array([[0., -1j], [1j, 0.]])

NameError: name 'np' is not defined

Learning goal: Predict which terms in the nanowire model lead to which consequences for Majoranas.

  • Zeeman splitting allows to get rid of one spin.
  • Singlet pairing induces superconductivity, but can never open a gap on its own.
  • Spin-orbit coupling is only needed to help open the gap.
  • Size of the Majoranas depends on all the parameters, and the transition requires tuning $B_z$ and $\mu$.

NO CONTENT ABOVE THIS LINE IS VISIBLE IN EDX

From Kitaev model to an experiment

We have a special guest to begin this week's lecture, Yuval Oreg from the Weizmann Institute in Rehovot.


In [ ]:
MoocVideo("GQLfs4i22ms", src_location="2.1-intro")

Small parameters

We are now all set to make Majoranas in a real system. Or at least to invent a way to make Majoranas in a real system.

The way we approach this problem is by considering the Kitaev chain a 'skeleton', and 'dressing' it with real physics phenomena until it becomes real.

Interestingly, this is not at all how the condensed matter community came to this model. Instead, the path to it was from complex to simple. The whole story started from what we'll consider in the very end of the course, the fractional particles.

Then it was simplified to topological superconductors (that still do not exist in real nature, as far as we know). Then Majoranas were predicted to exist (week 7) in a combination of a 3D topological insulator (week 6), then that part was simplified to a two-dimensional topological insulator (week 5), and only after a few more simplification steps, the nanowire model was developed.

So once again, here is our 'skeleton', the Kitaev model Hamiltonian written in momentum space:

$$H_{Kitaev} = (-2 t \cos k -\mu) \tau_z + 2 \Delta \tau_y \sin k$$

The model seems OK for a start, because it has some superconducting pairing $\Delta$ and some normal dispersion given by terms proportional to $\mu$ and $t$.

Before we proceed further, let's understand the relation between these parameters.

First of all, we want to make a controllable system, so that we can tweak its parameters. That means that we need a semiconductor. In semiconductors the electron density is very low, so that the chemical potential is near the bottom of the band. This makes it easier to define $\mu$ with respect to the bottom of the band:

$$\mu \rightarrow \mu - 2t.$$

Now the transition between trivial and non-trivial states happens when $\mu = 0$.

Of course semiconductors are never additionally superconducting. Luckily this is easy for us to resolve. We just paste a superconductor and semiconductor together into a hybrid structure, and let the superconductor induce superconductivity in the semiconductor. Making such a hybrid is an extremely hard challenge from the material science point of view, but it's definitely not our problem for now.

The next thing we should consider is that $\mu$ will always stay small compared to the bandwidth, so $\mu \ll 2t$. The same holds for superconducting pairing: $\Delta \ll t$. This is because superconductivity is a very weak effect compared to kinetic energy of electrons. These two inequalities combined mean that we can expand the $\cos k$ term and only work with the continuum limit of the Kitaev model:

$$H = (k^2/2m - \mu) \tau_z + 2 \Delta \tau_y k.$$

The effective electron mass $m$ is just the coefficient of the expansion. Let's take a look at the band structure in this regime, both in the topological regime and in the trivial regime:


In [ ]:
# matplotlib.rcParams['axes.color_cycle'] = 'k'

f = plt.figure(figsize=[9, 3.5])
ax0 = f.add_subplot(1, 2, 1)
ax0.set_xlabel('$k$')
ax0.set_xticks([-1.0, 0.0, 1.0])
ax0.set_ylabel('$E$')
ax0.set_xlim([-1.03, 1.03])
ax0.set_ylim([-1.0, 1.0])
ax0.set_yticks([-1.0, 0.0, 1.0])
ax0.set_title('Trivial bandstructure')
ax1 = f.add_subplot(1, 2, 2)
ax1.set_xlabel('$k$')
ax1.set_xticks([-1.0, 0.0, 1.0])
ax1.set_ylabel('$E$')
ax1.set_xlim([-1.03, 1.03])
ax1.set_ylim([-1.0, 1.0])
ax1.set_yticks([-1.0, 0.0, 1.0])
ax1.set_title('Topological bandstructure')
sys0 = spinful_kitaev_chain(L=None, t=1.0, delta=0.1, mu=-0.3, B=0)
sys1 = spinful_kitaev_chain(L=None, t=1.0, delta=0.1, mu=0.3, B=0)

ax0.set_color_cycle(['k'])
ax1.set_color_cycle(['k'])

kwant.plotter.bands(sys0, momenta=201, ax=ax0)
kwant.plotter.bands(sys1, momenta=201, ax=ax1)

The need for spin

Still, there is one obvious thing missing from the model, and it is the electron spin. This model works with some hypothetical spinless fermions, that are nowhere to be found. So to make the model physical, we need to remember that every single particle has spin, and the Hamiltonian has some action in spin space, described by Pauli matrices $\sigma$.

The simplest thing which we can do is to just add the spin as an extra degeneracy, so to multiply every term in the Hamiltonian by $\sigma_0$. Obviously this doesn't change the spectrum, and a zero energy solution stays a zero energy solution.

Just kidding, this would be very bad! The problem about adding spin is that the whole point of a Kitaev chain is to create unpaired Majorana modes. Now if we add an extra spin degeneracy to these Majoranas, the edge of our chain now hosts two Majoranas, or in other words one regular fermion fine-tuned to zero energy.

What's the correct way of introducing spin then? We still need to add the spin, but let's make the Kitaev chain corresponding to one spin species topologically trivial, and another one non-trivial. We know that the chemical potential $\mu$ controls if the chain is topological, so if say spin up has $\mu > 0$ and spin down $\mu < 0$, we're back in business.

We achieve this by adding Zeeman coupling of the spin to an external magnetic field:

$$H = (k^2/2m - \mu - B \sigma_z) \tau_z + 2 \Delta \tau_y k.$$

Whenever the Zeeman energy $|B|$ is larger than $\mu$ we have one Majorana fermion at the end of the chain.

Let's look at what happens with the dispersion as we increase the magnetic field from zero to a value larger than $\mu$.


In [ ]:
# Plot of the Kitaev dispersion with Zeeman as a parameter. 
# mu should be a few times larger than the gap, so that there is a mexican hat dispersion.
def kitaev_bands(B):
    sys = spinful_kitaev_chain(None, 1.0, 0.1, 0.3, B)
    momenta   = np.linspace(-1.03, 1.03, 201)
    bands = kwant.physics.Bands(sys)
    energies = [bands(k) for k in momenta]
    
    fig = plt.figure()
    plt.plot(momenta, energies, 'k-')
    
    plt.xlabel('$k$')
    plt.xticks([-1.0, 0.0, 1.0])
    plt.ylabel('$E$')
    plt.xlim([-1.03, 1.03])
    plt.ylim([-1.0, 1.0])
    plt.yticks([-1.0, 0.0, 1.0])
    plt.title('Band structure')
    return fig

StaticInteract(lambda B: kitaev_bands(B*0.05),
               B = RangeWidget(0, 8))

We now see that we resolved the first problem:

A high enough Zeeman splitting allows to separate the different spins. Then we can make one spin species trivial, while another one is topological and hosts Majoranas.

Realistic superconducting pairing

The next part for us to worry about is the superconductor.

Something that you probably saw in the Kitaev chain Hamiltonian is that the superconducting pairing $\Delta$ has a peculiar form. It pairs electrons from neighboring sites, and not those from the same site. In a momentum space this means that the superconducting pairing is proportional to $\Delta k$.

Of course, in a Kitaev chain the superconducting pairing cannot couple two electrons from the same site since there is just one particle per site!

The real world superconductors are different. Most of them, and specifically all the common superconductors like $Al$, $Nb$, $Pb$, $Sn$ have the $s$-wave pairing. This means that the pairing has no momentum dependence, and it is local in real space. The Kitaev chain pairing is proportional to the first power of momentum and so it is a $p$-wave pairing.

High temperature superconductors like cuprates or pnictides do have a momentum-dependent pairing, but it's a yet another type ($d$-wave, or a more exotic $s\pm$-wave).

So if we want to invent a way to make Majoranas, we will need to use the $s$-wave pairing. And then, as you should remember from the previous week, due to the fermionic statistics the pairing function should be antisymmetric. In a Kitaev chain the antisymmetry was due to the real space structure of the pairing, but in an $s$-wave superconductor, the antisymmetry of the pairing should arise due to its spin structure.

This leaves only one option. All the $s$-wave superconductors are spin-singlet:

$$H_{pair} = \Delta(c_\uparrow c_\downarrow - c_\downarrow c_\uparrow) + \text{h.c.}$$

This means that now we need to modify the pairing, but before that we'll need to do one other important thing.

Important and useful basis change.

When you see Bogoliubov-de-Gennes Hamiltonians in the literature, you will find them written in two different bases. One variant is the one which we introduced last week:

$$ H_\textrm{BdG} = \begin{pmatrix} H & \Delta \\ -\Delta^* & -H^* \end{pmatrix}. $$

It has the particle-hole symmetry $H_\textrm{BdG} = - \tau_x H^*_\textrm{BdG} \tau_x$. In this basis, the $s$-wave pairing is proportional to $\sigma_y$.

However for systems with complicated spin and orbital structure, there is a different basis which makes the book-keeping much easier.

If we have a time-reversal symmetry operator $\mathcal{T} = U \mathcal{K}$, we can apply the unitary transformation $U$ to the holes, so that in the new basis we get the Bogoliubov-de-Gennes Hamiltonian

$$ H_\textrm{BdG} = \begin{pmatrix} H & \Delta' \\ \Delta'^\dagger & -\mathcal{T} H \mathcal{T}^{-1}\end{pmatrix}, $$

with $\Delta' = \Delta U^\dagger$.

Why is this basis useful?

  • First of all, because in this new basis the $s$-wave pairing is a unit matrix regardless of system we consider.
  • Second, because it's easy to get the Hamiltonian of holes. We take the Hamiltonian for electrons, changing the signs of all terms that respect time-reversal symmetry, but not for those that break it, such as the term proportional to the magnetic field $B$. So if the electrons have a Hamiltonian $H(B)$, the Hamiltonian of the holes just becomes $-H(-B)$.

There is one disadvantage. The particle-hole symmetry now becomes more complicated. For our system with only one orbital and spin it is $\mathcal{P} = \sigma_y \tau_y \mathcal{K}$. But, let us tell you, the advantages are worth it.

s-wave superconductor with magnetic field

Let's look at how our chain looks once we change the superconducting coupling to be $s$-wave. The Zeeman field (or anything of magnetic origin) changes sign under time-reversal symmetry.

This means that the Zeeman field has the same form for electrons and for holes in the new basis, and the full Hamiltonian is now:

$$ H_\textrm{BdG} = (k^2/2m - \mu)\tau_z + B \sigma_z + \Delta \tau_x. $$

This Hamiltonian is easy to diagonalize since every term either has only a $\tau$ matrix or a $\sigma$ matrix. At $k=0$ it has 4 levels with energies $E = \pm B \pm \sqrt{\mu^2 + \Delta^2}$.

We can use this expression to track the crossings. We also know that when $B=0$ the system is trivial due to spin degeneracy. Together this means that we expect the system to be non-trivial (and will have negative Pfaffian invariant) when

$$ B^2 > \Delta^2 + \mu^2.$$

Are we now done? Not quite.

Problem with singlets

A singlet superconductor has an important property: Since electrons are created in singlets, the total spin of every excitation is conserved. Zeeman field conserves the spin in $z$-direction, so together every single state of our system has to have a definite spin, including the Majoranas.

And that is a big problem. Majoranas are their own particle-hole partners, and that means that they cannot have any spin (energy, charge, or any other observable property at all).

So does this now mean that we "broke" the bulk-edge correspondence? Let's look at the band structure (tweak the Zeeman energy):


In [ ]:
StaticInteract(lambda B: spinorbit_band(1.0, 0.0, 0.1, 0.04*B, 0.0),  
        B = RangeWidget(0, 5, 1))

Of course we didn't break bulk-edge correspondence. Majoranas in our system would have to have a spin, which isn't possible. That in turn means that they cannot appear, and that means that the system cannot be gapped.

We can also approach this differently. From all the spin Pauli matrices, only $\sigma_z$ appears in the Hamiltonian, and there's a conservation law. The two bands that cross at zero energy in the band structure above belong to opposite spin bands, and so cannot be coupled.

Now we need to solve this final problem before we are done.

How to open the gap?

The final stretch is straightforward.

We know that there is no gap because of conservation of one of the spin projections, so we need to break the spin conservation.

If we don't want to create an inhomogeneous magnetic field, we have to use a different term that couples to spin. That term is spin-orbit interaction. In it's simplest form this interaction appears in our wire as

$$H_{SO} = \alpha \sigma_y k,$$

So it is like a Zeeman field pointing in $y$-direction with a strength proportional to the particle momentum. Note that this term is invariant under time reversal symmetry (both $\sigma_y$ and $k$ change sign). So now we have our final Hamiltonian:

$$ H_\textrm{wire} = (k^2/2m + \alpha \sigma_y k - \mu)\tau_z + B \sigma_z + \Delta \tau_x. $$

At $k = 0$, spin-orbit coupling vanishes, so it has no effect on the system being topologically trivial or non-trivial.

Let's now check that it does what we want, namely open the gap at a finite momentum:


In [ ]:
StaticInteract(lambda alpha: spinorbit_band(1., 0.1, 0.1, 0.2, 0.05*alpha), 
               alpha=RangeWidget(0, 8, 1, default=0))

Yep, it does :)

An important remark: You might now think that since spin-orbit depends on spin, spin-orbit we would make the magnetic field unnecessary. This is not true: Since spin-orbit preserves time-reversal symmetry, in the absence of a magnetic field the energy spectrum of the model would have a Kramers degeneracy, as you learned last week. To get one unpaired Majorana mode per edge and not two, we need to break Kramers degeneracy and therefore break time-reversal symmetry. So the combination of both Zeeman field and spin-orbit coupling is needed.

Putting everything together

Let's now rest for a moment and reflect on what we have done.

We started from a toy model, which has a very special feature. Then one by one we fixed the parts of the model that we found unrealistic and we ended up with a new system, which still has a relatively simple Hamiltonian, but already gives a hope of being realizable in a lab.

Now try to guess: how many papers were written studying this exact model? The exact number is hard to obtain, but the count is in hundreds!

Despite the model being very simple and the fact that it can be written in one line, it has four independent parameters already in our simplest formulation. Let's enumerate the parameters once again:

  • The chemical potential $\mu$, which sets the overall electron density in the wire.
  • The induced superconducting gap $\Delta$, which is required to make particle-hole symmetry play a role.
  • The spin-orbit coupling $\alpha$, which breaks spin conservation.
  • The Zeeman field $B$, which breaks Kramers degeneracy.

We need to control every single parameter out of these 4 to create Majoranas (and there are even more). This is why the task of creating Majoranas is extremely challenging.

As a final point in our story, let's see how the four parameters work together in determining how large the gap in our system is.

Obviously, this is the key parameter that we care about when creating Majoranas. The smaller the gap, the worse the protection of Majoranas, and the more we need to worry about the effects of finite temperature.

Let's calculate the gap as a function of all of the relevant parameters.


In [ ]:
B_vals = np.linspace(0, 0.3, 71)
StaticInteract(lambda mu: spinorbit_band_gap(mu*0.025, 1.0, 0.1, B_vals), 
        mu=RangeWidget(-2, 6, 1))

Here the vertical line denotes the critical value of Zeeman field at which the wire becomes topological.

Let's summarize our observations:

  • So we see that the closer $\mu$ is to 0, the lower $B$ is required to reach the topological regime.
  • After reaching the topologically nontrivial regime, the gap slowly grows as we go away from the transition region, and after reaching its peak value, it starts dropping.
  • Finally, we see that the higher the spin-orbit coupling, the larger the optimal gap in the topological regime.

We finish our investigation of this model for now with a final simple picture of the band structure of our system.


In [ ]:
StaticInteract(lambda mu: spinorbit_band(t=1.0, mu=0.025*mu + 0.02, delta=0.025, B=0.07, alpha=0.8), 
        mu=RangeWidget(-8, 8, 0.5, default=-8))

When $\mu$ is very negative we see two split electron bands at positive energy corresponding to two spin orientations.

The lower of these two bands has a characteristic double minimum due to spin-orbit coupling.

As we increase $\mu$, the bands move down in energy, until they couple with the hole bands at $E=0$. This only happens due to the combination of superconductivity and spin-orbit.

At $k=0$ the spin-orbit is ineffective, so the electron and hole bands pass through each other, changing the system first from trivial to topological and then back.

The non-monotonous behavior of the gap versus $B$ that we saw earlier is a consequence of this complicated band structure: There are different values of momenta where the dispersion has local minima. When we are close to the phase transition, $k=0$ defines the gap, while for large $B$, it is the gap at finite momentum that becomes smallest.


In [ ]:
question = ("What happens if we align the magnetic field $B$ along the $y$-direction instead of the $z$-direction?")
answers = ["Then we do not need spin-orbit coupling anymore in order to get Majoranas.",
           "Then spin projection along $y$ direction is conserved, so we can't get Majoranas.",
           "It's impossible, because a magnetic field can only be applied along $z$.",
           "Then the spin-orbit term is automatically modified to point along the $z$ direction, so nothing really changes."]
explanation = ("If both the magnetic field and the spin orbit coupling point in the $y$ direction, " +
               "then the Hamiltonian commutes with $\sigma_y$, and spin projection along $y$ is a good quantum number. " +
               "So we are back to the problem that a gap at finite momentum does not open, " +
               "and we do not get a topological phase supporting Majoranas.")
MoocMultipleChoiceAssessment(question, answers, correct_answer=1, explanation=explanation)

Outlook


In [ ]:
MoocVideo("MsFyJBAMFLI", src_location="2.1-summary")

In [ ]:
MoocDiscussion("Questions", "Majoranas in nanowires")