Conductance of a superconducting Coulomb blockaded Majorana nanowire

TL;DR

This paper develops a comprehensive theoretical model for conductance in superconducting Coulomb blockaded nanowires hosting Majorana bound states, explaining experimental observations and the effects of temperature, subbands, and Andreev states.

Contribution

It introduces a generalized finite temperature master equation including multiple subbands, Andreev bound states, and Majorana states, providing a complete understanding of tunneling conductance.

Findings

Conductance peak spacings oscillate with magnetic field, matching some experimental trends.

Inclusion of all effects shows oscillation amplitudes decrease with magnetic field at high temperatures.

The theory explains experimental results only under specific conditions involving temperature and state coexistence.

Abstract

In the presence of an applied magnetic field introducing Zeeman spin splitting, a superconducting (SC) proximitized one-dimensional (1D) nanowire with spin-orbit coupling can pass through a topological quantum phase transition developing zero-energy topological Majorana bound states (MBSs) on the wire ends. One of the promising experimental platforms in this context is a Coulomb blockaded island, where by measuring the two-terminal conductance one can in principle investigate the MBS properties. We theoretically study the tunneling transport of a single electron across the superconducting Coulomb blockaded nanowire at finite temperature to obtain the generic conductance equation. By considering all possible scenarios where only MBSs are present at the ends of the nanowire, we compute the nanowire conductance as a function of the magnetic field, the temperature, and the gate voltage. In the simplest 1D topological SC model, the oscillations of the conductance peak spacings (OCPSs) arising from the MBSs overlap from the two wire ends manifest an increasing oscillation amplitude with increasing magnetic field (in agreement with theories without Coulomb blockade and in disagreement with a recent experimental observation). We develop a generalized finite temperature master equation theory including not only multiple subbands in the nanowire, but also the possibility of ordinary Andreev bound states in the non-topological regime. Inclusion of all four effects (temperature, multiple subbands, Andreev bound states, and MBSs) provides a complete picture of the tunneling transport properties. Based on this complete theory, we indeed obtain OCPSs whose amplitudes decrease with increasing magnetic field in qualitative agreement with recent experimental results, but this happens only for rather high temperatures with multisubband occupancy and the presence of both Andreev bound states and MBSs.

Equations

$$U(N)=$$

$$I=$$

$$\displaystyle\Big{\{}\Gamma_{p}^{l}\big{[}P_{N-1}(\{n_{i}\})\delta_{n_{p},0}f(\epsilon^{l}_{p})-P_{N}(\{n_{i}\})\delta_{n_{p},1}(1-f(\epsilon^{l}_{p}))\big{]}$$

$$\displaystyle+\Lambda_{p}^{l}\big{[}P_{N-1}(\{n_{i}\})\delta_{n_{p},1}f(\tilde{\epsilon}^{l}_{p})-P_{N}(\{n_{i}\})\delta_{n_{p},0}(1-f(\tilde{\epsilon}^{l}_{p}))\big{]}\Big{\}}$$

$$0=$$

$$=$$

$$\displaystyle-\sum_{p}P_{N}(\{n_{i}\})\delta_{n_{p},0}\Big{\{}\Lambda_{p}^{l}(1-f(\tilde{\epsilon}_{p}^{l}))+\Lambda_{p}^{r}(1-f(\tilde{\epsilon}_{p}^{r}))\Big{\}}$$

$$\displaystyle+\sum_{p}P_{N-1}(n_{p}=0)\delta_{n_{p},1}\Big{\{}\Gamma_{p}^{l}f(\epsilon_{p}^{l})+\Gamma_{p}^{r}f(\epsilon_{p}^{r})\Big{\}}$$

$$\displaystyle+\sum_{p}P_{N-1}(n_{p}=1)\delta_{n_{p},0}\Big{\{}\Lambda_{p}^{l}f(\tilde{\epsilon}_{p}^{l})+\Lambda_{p}^{r}f(\tilde{\epsilon}_{p}^{r})\Big{\}},$$

$$0=$$

$$=$$

$$\displaystyle-\sum_{p}P_{N-1}(\{n_{i}\})\delta_{n_{p},1}\Big{\{}\Lambda_{p}^{l}f(\tilde{\epsilon}_{p}^{l})+\Lambda_{p}^{r}f(\tilde{\epsilon}_{p}^{r})\Big{\}}$$

$$\displaystyle+\sum_{p}P_{N}(n_{p}=1)\delta_{n_{p},0}\Big{\{}\Gamma_{p}^{l}(1-f(\epsilon_{p}^{l}))+\Gamma_{p}^{r}(1-f(\epsilon_{p}^{r}))\Big{\}}$$

$$\displaystyle+\sum_{p}P_{N}(n_{p}=0)\delta_{n_{p},1}\Big{\{}\Lambda_{p}^{l}(1-f(\tilde{\epsilon}_{p}^{l}))+\Lambda_{p}^{r}(1-f(\tilde{\epsilon}_{p}^{r}))\Big{\}},$$

$$P_N^{\mathrm{eq}}(n_p=1)e^{\beta {ϵ}_p}=$$

$$P_N^{\mathrm{eq}}(n_p=0)e^{\beta {\widetilde{ϵ}}_p}=$$

$$\displaystyle P_{N}^{\rm{eq}}(\{n_{i}\})=Z^{-1}\exp\Big{[}-\beta(\sum_{i}E_{i}n_{i}+U(N))\Big{]}$$

$$P_N(\{n_i\})=P_N^{\mathrm{eq}}(\{n_i\})(1+\beta eV{\Phi }_N(\{n_i\})),$$

$$\frac{dI}{dV}=-\beta e^2\sum _p\sum _{\{n_i\}}$$

$$\displaystyle\bigg{\{}\delta_{n_{p},0}P_{N-1}^{\rm{eq}}(\{n_{i}\})f(\epsilon_{p})\Gamma_{p}^{l}\Big{[}\Phi_{N-1}(\{n_{i}\})-\eta f(-\epsilon_{p})\Big{]}-\delta_{n_{p},1}P_{N}^{\rm{eq}}(\{n_{i}\})\Gamma_{p}^{l}\Big{[}\Phi_{N}(\{n_{i}\})\big{(}1-f(\epsilon_{p})\big{)}+\eta f(\epsilon_{p})f(-\epsilon_{p})\Big{]}$$

$$\displaystyle+\delta_{n_{p},1}P_{N-1}^{\rm{eq}}(\{n_{i}\})f(\tilde{\epsilon}_{p})\Lambda_{p}^{l}\Big{[}\Phi_{N-1}(\{n_{i}\})-\eta f(-\tilde{\epsilon}_{p})\Big{]}-\delta_{n_{p},0}P_{N}^{\rm{eq}}(\{n_{i}\})\Lambda_{p}^{l}\Big{[}\Phi_{N}(\{n_{i}\})\big{(}1-f(\tilde{\epsilon}_{p})\big{)}+\eta f(\tilde{\epsilon}_{p})f(-\tilde{\epsilon}_{p})\Big{]}\bigg{\}}$$

$$\frac{dI}{dV}=$$

$$\displaystyle+\delta_{n_{p},1}P_{N-1}^{\rm{eq}}(\{n_{i}\})f(\tilde{\epsilon}_{p})\Lambda_{p}^{l}\Big{[}\Phi_{N}(n_{p}=0,\{n_{i\neq p}\})-\Phi_{N-1}(n_{p}=1,\{n_{i\neq p}\})+\eta\Big{]}\bigg{\}}$$

$$0$$

$$\displaystyle+\sum_{p}P_{N-1}^{\rm{eq}}(n_{p}=1,\{n_{i\neq p}\})f(\epsilon_{p})(\Lambda_{p}^{l}+\Lambda_{p}^{r})\Big{(}\Phi_{N}(n_{p}=0,\{n_{i\neq p}\})-\Phi_{N-1}(n_{p}=1,\{n_{i\neq p}\})+\eta-\frac{\Lambda^{r}_{p}}{\Lambda_{p}^{l}+\Lambda_{p}^{r}}\big{)}\delta_{n_{p},0}$$

$$0$$

$$\displaystyle+\sum_{p}P_{N-1}^{\rm{eq}}(n_{p}=1,\{n_{i\neq p}\})f(\epsilon_{p})(\Lambda_{p}^{l}+\Lambda_{p}^{r})\Big{(}\Phi_{N}(n_{p}=0,\{n_{i\neq p}\})-\Phi_{N-1}(n_{p}=1,\{n_{i\neq p}\})+\eta-\frac{\Lambda^{r}_{p}}{\Lambda_{p}^{l}+\Lambda_{p}^{r}}\big{)}\delta_{n_{p},1}$$

$$\displaystyle\frac{\partial}{\partial t}\Big{(}\sum_{\{n_{i}\}}P_{N}{\{n_{i}\}}+\sum_{\{n_{i}\}}P_{N-1}{\{n_{i}\}}\Big{)}=0$$

$$\displaystyle P_{N}(n_{p}=1)\Big{[}\Gamma^{l}_{p}(1-f(\epsilon^{l}_{p}))+\Gamma_{p}^{r}(1-f(\epsilon^{r}_{p}))\Big{]}$$

$$=$$

$$\displaystyle P_{N}(n_{p}=0)\Big{[}\Lambda^{l}_{p}(1-f(\tilde{\epsilon}^{l}_{p}))+\Lambda_{p}^{r}(1-f(\tilde{\epsilon}^{r}_{p}))\Big{]}$$

$$=$$

$$\displaystyle\beta eP_{N}^{\rm{eq}}(n_{p}=1)\bigg{\{}\Phi_{N}(n_{p}=1,\{n_{i\neq p}\})(\Gamma_{p}^{l}+\Gamma_{p}^{r})(1-f(\epsilon_{p}))+f(\epsilon_{p})f(-\epsilon_{p})\Big{[}\eta\Gamma_{p}^{l}-(1-\eta)\Gamma_{p}^{r}\Big{]}\bigg{\}}$$

$$=$$

$$\displaystyle\beta eP_{N}^{\rm{eq}}(n_{p}=0)\bigg{\{}\Phi_{N}(n_{p}=0,\{n_{i\neq p}\})(\Lambda_{p}^{l}+\Lambda_{p}^{r})(1-f(\tilde{\epsilon}_{p}))+f(\tilde{\epsilon}_{p})f(-\tilde{\epsilon}_{p})\Big{[}\eta\Lambda_{p}^{l}-(1-\eta)\Lambda_{p}^{r}\Big{]}\bigg{\}}$$

$$=$$

$$0={\Phi }_N(n_p=1,\{n_{i\ne p}\})-{\Phi }_{N-1}(n_p=0,\{n_{i\ne p}\})+\eta -\frac{{\Gamma }_p^r}{{\Gamma }_p^l+{\Gamma }_p^r},$$

$$0={\Phi }_N(n_p=0,\{n_{i\ne p}\})-{\Phi }_{N-1}(n_p=1,\{n_{i\ne p}\})+\eta -\frac{{\Lambda }_p^r}{{\Lambda }_p^l+{\Lambda }_p^r},$$

$$\frac{dI}{dV}=$$

$$=$$

$$\frac{dI}{dV}=\frac{e^{2}\alpha}{kT}\frac{\sum_{p}\big{[}A_{\rm{odd}}(\epsilon_{p})+B_{\rm{even}}(\tilde{\epsilon}_{p})\big{]}}{F_{\rm{even}}e^{-\beta\Delta U}+F_{\rm{odd}}}$$

$$A_{even/odd}({ϵ}_p)$$

Full text

I Introduction

The integer quantum Hall effect Klitzing et al. (1980) ushered in the era of topological systems and phenomena in condensed matter physics, although the fact that the precise quantization of the Hall conductance in two dimensions (2D) is indeed a direct manifestation of topological robustness took several years to be appreciated. Of course, the 2016 Nobel Prize Thouless et al. (1982); Haldane (1988) in physics has made this fact rather universally celebrated. The classification of topological insulators and superconductors Kitaev (2009); Schnyder et al. (2008) provides guidance to look for topological systems and materials. The essential signature of the topological phase is the presence of stable gapless (zero energy) states on the boundary of the system with the bulk having a robust energy gap. In fact, this can be construed as an equivalent definition of a topological phase for the quantum Hall effect, the bulk gap corresponds to the cyclotron gap imposed by the external magnetic field whereas the boundary gapless states are the edge states confined to the 1D edge of the 2D layer. Indeed , the symmetry-based topological classification scheme Hasan and Kane (2010); Qi and Zhang (2010); Chiu et al. (2016) is restricted to insulators (e.g. 2D quantum Hall states, 2D quantum spin Hall insulators, 3D topological insulators) and superconductors simply because these are the systems with bulk energy gaps separating ground states from excited states. But, obviously, very special constraints are necessary to ensure that an insulator or a superconductor would have gapless or zero energy boundary states, and this is why the subject of topological systems has become active only in recent years because of deep theoretical advances in spite of insulators and superconductors having been known for more than hundred years. Insulators and superconductors all have bulk gaps, but only the ones having robust boundary gapless states are called topological (whereas those not having such special boundary states are called nontopological or trivial).

The topological quantum phase transition (TQPT) between topological and trivial phases is also a subject of great current interest, and such a transition can only happen through the vanishing of the bulk gap at the TQPT. Recently, there has been an enormous inter-disciplinary interest in one particular type of zero energy boundary states associated with both 1D and 2D topological superconductor systems. These are the so-called Majorana (zero-energy) bound states (MBSs) Elliott and Franz (2015); Kitaev (2001), which form the central theme of the current work. In particular, these MBSs are strange quantum objects obeying anyonic non-Abelian braiding statistics Sarma et al. (2015); Nayak et al. (2008); Stanescu (2016), which can be used for fault-tolerant quantum computation, thus bringing together physicists, mathematicians, computer scientists, electrical engineers, and materials scientists in an effort to build such a quantum computer. Much popular interest has recently focused on the possibility of Majorana-based topological quantum computation since Microsoft Corporation has just announced a large commercial effort to build such a computer. An MBS, which is its own anti-particle, has zero energy protected by particle-hole symmetry stemming from superconductivity. Furthermore, these MBSs can be realized on the ends of 1D topological superconductors Kitaev (2001) and at point defects of 2D topological superconductors Fu and Kane (2008). The most promising experimental platform for the realization of MBSs is SC-proximitized semiconductor nanowire with strong spin-orbit coupling Lutchyn et al. (2010); Oreg et al. (2010); Sau et al. (2010a) in an applied magnetic field to create a Zeeman spin splitting. The precise theoretical prediction Lutchyn et al. (2010); Oreg et al. (2010) is that a spin-orbit-coupled semiconductor nanowire (e.g. InSb, InAs) with superconducting proximity effect induced by a nearby regular superconductor (e.g. Nb, Al) would become a topological superconductor with zero-energy (i.e. mid-gap) MBSs localized at the wire ends provided the Zeeman spin splitting induced by the applied magnetic field is large enough to overcome the induced SC gap. Here, the applied magnetic field is the tuning parameter inducing the TQPT with MBSs appearing as localized zero energy bound states at the wire ends at high enough magnetic field values. One technique to experimentally look for MBSs in the nanowires is to study the tunneling conductance Sau et al. (2010a); Sengupta et al. (2001), which should manifest a quantized zero bias peak associated with the MBS in the topological phase. The recent experimental observation of such a predicted zero-bias conductance peak in nanowires above a critical applied magnetic field has been a milestone to hint at the possible existence of the MBSs in the nanowire-superconductor hybrid structure Mourik et al. (2012); Das et al. (2012); Zhang et al. (2016); Deng et al. (2012); Churchill et al. (2013); Finck et al. (2013); Chen et al. (2016). These experimental observations of zero bias conductance peaks in nanowires in a finite external magnetic field by multiple independent laboratories using different semiconductor (InSb or InAs)-superconductor (NbTiN or Al) combinations have created a great deal of excitement because of the implication that the conductance peak is providing strong evidence for the existence of topological Majorana modes in these nanowires. But the non-Abelian nature of these possible MBSs localized in nanowires still remains to be demonstrated experimentally, and this area is one of the most active current research areas in all of physics.

A key development in this subject is the recent experimental paper by Albrecht et al. Albrecht et al. (2016) reporting the apparent observation of exponential “topological protection” in Majorana nanowires, which, if validated and understood, would be a singular landmark in the field. In particular, the most straightforward interpretation of the observation of Albrecht et al. Albrecht et al. (2016) is that the two MBSs localized at the two ends of the nanowire in their system are sufficiently far apart so that their wavefunction overlap is exponentially small. This exponential weakness is reflected in the MBS splitting oscillation showing an exponential decrease with increasing wire length as predicted theoretically Cheng et al. (2009, 2012); Das Sarma et al. (2012); Hützen et al. (2012). This would imply that each MBS can now be thought of as an independent topological entity obeying non-Abelian statistics and hence suitable for use in topological quantum computation. Unfortunately, however, such a straightforward interpretation seems inapplicable to the experiment of Ref. Albrecht et al. (2016) since the magnetic field dependence of the MBS overlap seems to disagree with the theoretical predictions Das Sarma et al. (2012); Stanescu et al. (2013) in spite of the length dependence manifesting the predicted theoretical exponential behavior. The current work is aimed at an understanding of the Albrecht experiment Albrecht et al. (2016), which because of its singular importance (i.e. “topological protection”) must be taken extremely seriously.

The conundrum here is the following. The theoretical exponential behavior Das Sarma et al. (2012), as manifesting in the oscillations of the conductance peak spacings (OCPSs) as a function of the applied magnetic field or the wire length, reflects an $e^{-L/\xi}$ dependence in the overlap between the MBSs localized at the two wire ends, where $L$ and $\xi$ are respectively the wire length (or more precisely, the separation between the two MBSs) and the Majorana localization length (or more precisely, the nanowire superconducting coherence length). Since the coherence length increases with increasing magnetic field B Das Sarma et al. (2012), the Majorana oscillations in the wire length and magnetic field are intimately coupled. An observation of the exponential decrease of oscillations in the wire length (as reported in Ref. Albrecht et al. (2016)) must therefore automatically come with an increase in the oscillation amplitude as a function of magnetic field since $L/\xi$ decreases with increasing magnetic field at fixed $L$. Seeing one without the other does not make any sense from the perspective of the minimal theory Das Sarma et al. (2012). The experimental situation in Ref. Albrecht et al. (2016) is actually worse since the oscillation amplitude seems to decrease (instead of increasing) with increasing magnetic field (while at the same time, the amplitude decreases with increasing L), which is completely the opposite of the predicted theoretical behavior. A resolution of this conundrum in Ref. Albrecht et al. (2016) is obviously of key importance in the context of the “exponential protection” claim in terms of wire length. The hope for understanding the puzzling results reported in Ref. Albrecht et al. (2016) (i.e. apparent exponential decrease in MBS overlap as a function of wire length along with a decreasing overlap as a function of increasing magnetic field) lies in the fact that the experimental system in Ref. Albrecht et al. (2016) is not consistent with the minimal model of a nanowire coupled to a SC as considered in most MBS theories, but presents a more complex situation (discussed below) involving a Coulomb blockade in the nanowire. Our goal in the current work is to generalize Majorana theories to include Coulomb blockade to see if the results of Ref. Albrecht et al. (2016) can be explained and understood.

Coulomb blockade at the basic level means that the system is small enough so that the Coulomb charging energy for putting electrons into the system is significant in affecting the experimental behavior. We must therefore incorporate Coulomb blockade in the Majorana nanowire theory to develop an appropriate model for the situation studied Ref. Albrecht et al. (2016). This new model inspired by Ref. Albrecht et al. (2016) includes a spin-orbit coupled proximitized SC nanowire in the presence of non-SC quantum dots at wire ends under an applied magnetic field. The study of Coulomb blockade in quantum dots (without any superconductivity) has a long history Kastner (1992). Coulomb blockade was first discovered experimentally by Dolan and Fulton in small metallic tunnel junctions manifesting charging oscillations associated with the finite Coulomb energy in small systems Fulton and Dolan (1987). Subsequently, Fulton et al. observed the interplay of superconductivity and charging energy in small SC tunnel junctions in 1989 Fulton et al. (1989). The periodic Coulomb oscillation of conductance peaks in semiconductor quantum dots was first observed also in 1989 Scott-Thomas et al. (1989) and supported by followup experiments Meirav et al. (1989, 1990). The phenomenon was explained by single-electron tunneling in the presence of Coulomb blockade physics Glazman and Shekhter (1989); Beenakker (1991); van Houten and Beenakker (1989). The single-electron Coulomb blockaded normal tunneling through small dots (with small capacitance) exhibits $1e$ periodicity of the conductance peak oscillation as the gate voltage varies. As the dot becomes superconducting and the SC gap is greater than the charging energy, $1e$ tunneling becomes blocked since single electrons cannot tunnel through SC gaps. It is, however, possible for small SC dots to manifest $2e$ Cooper pair tunneling Hekking et al. (1993) leading to $2e$ periodicity in the tunneling conductance oscillation as a function of gate voltage as observed in Refs. Eiles et al. (1993); Tuominen et al. (1992). It was also observed there can be a transition from $2e$ periodicity to $1e$ periodicity by increasing the temperature to suppress superconductivity Amar et al. (1994). In principle, the possibility exists in superconducting tunnel junctions at finite temperatures for both $2e$ and $1e$ charge oscillations to occur as a function of gate voltage since transport could take place through sub-gap and above-gap states. On the other hand Coulomb blockaded transport through zero energy MBSs at the wire ends in a 1D topological SC should manifest 1e charge oscillations since nonlocal resonant tunneling of single electrons are allowed through the MBSs Zazunov et al. (2011); Sau et al. (2015); Fu (2010); van Heck et al. (2016).

In addition to the experimental work by Albrecht et al. Albrecht et al. (2016), which motivates our work, there have been a few recent theoretical papers Zazunov et al. (2011); Sau et al. (2015); Fu (2010); van Heck et al. (2016) on Coulomb blockaded Majorana nanowires, but our work extends the theory to the realistic situation of finite temperature and multi-subband occupancy along with the inclusion of both ordinary Andreev bound states and Majorana states in contributing to transport so that a meaningful comparison to the important results of Ref. Albrecht et al. (2016) could be carried out. We first briefly discuss the qualitative expectations for tunneling transport in a Coulomb blockaded Majorana nanowire with increasing magnetic field so that the system evolves from being a non-topological ordinary (proximitized) SC at low field to being a topological SC at high fields, and then perhaps becoming a gapless SC at very high magnetic fields. At low magnetic field, the oscillation has $2e$ periodicity since the superconducting gap is larger than the charging energy. For higher magnetic field, while still being in the non-topological phase, the superconducting gap decreases eventually becoming smaller than the charging energy, and therefore, each conductance peak starts to split into two. Eventually, at high enough magnetic field the system goes through the TQPT, and in the presence of zero energy MBSs the conductance peak is expected to exhibit $1e$ periodicity since transport is now occurring through resonant tunneling through both MBSs. The possibility that at very high magnetic field, the 1e metallic periodicity may arise simply because the nanowire has developed a zero SC gap cannot be ruled out. In Ref. Albrecht et al. (2016), modifications to this strictly 1e resonant conduction due to MBS are claimed to be arising from the presence of MBS overlap from the two ends, and this modification leads to a measurement of the Majorana energy splitting, leading to the topological protection observed in Ref. Albrecht et al. (2016). Thus, it is crucial to understand how MBS overlap affects the resonant tunneling in the topological regime so that the topological protection claim of Ref. Albrecht et al. (2016) can be validated. This is particularly critical given that Ref. Albrecht et al. (2016) finds an exponential MBS splitting in the wire length, but inconsistent results as a function of magnetic field in the same wire as discussed above and also below.

The major problem in the Majorana interpretation is the inconsistency in the magnetic field behavior between the experimental observation and the theoretical prediction in the $1e$ periodic region. The oscillation amplitude of the even and odd conductance peak spacings observed in the experiment Albrecht et al. (2016) always decreases as the magnetic field increases. At low enough temperatures the conductance peak spacing should reflect the positive and negative energies of the two superconducting states close to zero energy van Heck et al. (2016) (i.e. the conductance should be an approximate map of the underlying mid-gap energy spectra near zero energy). The spectrum of a simple 1D SC proximitized semiconductor nanowire Stanescu et al. (2013) shows the oscillation amplitude of the two energy levels close to zero energy becoming larger as the magnetic field increases corresponding to an increasing MBS overlap due to the SC gap suppression with increasing field. Thus, if the Coulomb blockaded nanowire in the experiment possessing the MBSs is described by the 1D semiconductor model, the theoretical prediction and the experimental observation have a serious inconsistency in the oscillation amplitude. The goal of this manuscript is to resolve this inconsistency by considering possible scenarios in the presence or absence of the MBSs (i.e. by probing both the topological and the non-topological regimes of the nanowire). Several mechanisms could possibly lead to this observed damped oscillation with increasing field, such as finite temperature, small bulk gap, multi-subband contribution, and the presence of Andreev bound states in the topologically trivial region. To understand the experimental results, in particular the $1e$ periodicity region, which arises presumably from the MBSs localized in the nanowire, we develop a transport theory for a superconducting Coluomb blockaded nanowire in the presence of spin-orbit coupling and Zeeman splitting including effects of finite temperature, multi-subband occupancy, and Andreev bound states. Surprisingly, such a generic transport equation has not been derived in the literature in spite of the fact that Beenakker developed the non-equilibrium transport theory for a non-SC metallic Coulomb blockaded quantum dot 25 years ago Beenakker (1991). We find that such a transport theory description enables us to obtain the non-local conductance contributions resulting from both conventional Andreev bound state Sau et al. (2015) and topological MBSs Sau et al. (2015); Fu (2010); van Heck et al. (2016). Such contributions lead to qualitative distinctions between the conductance peak structure and the low-lying energy spectra even in the tunneling limit. This will lead us to reevaluate the simple interpretation of Majorana oscillations measured from Coulomb-blockaded resonances allowing an understanding of the experimental results, which are much more nuanced and subtle than the simple exponential protection in length scenario envisioned in Ref. Albrecht et al. (2016). Our theory obviously thus has wide-ranging consequences in the current search for non-Abelian Majorana modes in Coulomb blockaded nanowires in the quest toward building a topological quantum computer. We emphasize that the Coulomb blockade may arise simply from the superconductor-semiconductor hybrid structure being small in size without there being any explicit quantum dots being present in the system.

The rest of this manuscript is organized as follows. In section II, we derive the generic tunneling transport equation for a superconducting Coulomb blockaded nanowire at finite temperature in the $1e$ periodicity region with the tunneling leads weakly coupled to the nanowire. In section III, we describe the nanowire Hamiltonian along with the description for the leads, which are relevant for studying the system. In section IV, we numerically compute the OCPSs at different temperatures by using a simplified version of the general transport equations. We also present in section IV our main results in the context of the experimental findings in Ref. Albrecht et al. (2016), carefully searching for situations where the OCPSs could have decreasing amplitude with increasing magnetic field as observed experimentally. By considering several scenarios leading to the bulk gap shrinking as the magnetic field increases (as observed in Ref. Albrecht et al. (2016)), we compare the conductance peak spacings with the experimental observation. By introducing Andreev bound states in the non-topological (or trivial) region in the 1D superconducting nanowire model, we also look for an alternative explanation for the observed OCPSs in the absence of the MBSs. In the last part of section IV, we combine all of the aforementioned physics to find the best scenario that explains the conductance peak spacing oscillation in the experimental observation. We consider the length dependence of the conductance peak spacings in sec. V, concluding in sec. VI with a summary of our results and a discussion of the open questions.

II Transport formalism for Coulomb blockaded superconducting nanowires

In this section we develop a superconducting generalization of the Coulomb blockaded transport theory originally considered for metallic islands Beenakker (1991). Consider the two ends of a SC quantum dot or nanowire that are weakly coupled with the left lead with zero voltage and the right lead with a small applied voltage $V$. The voltage on the nanowire is described by $\eta V$ with the indeterminate factor $\eta$ dropping out of our final equations. We label $E_{p}$ in ascending order as the energy levels of the quasi-particle of the nanowire with respect to the BCS ground state. The possibility of adding hole excitations, which are electron excitations combined with the loss of a Cooper pair to the SC (i.e. the Andreev process) is what distinguishes the SC case from the metallic dot considered by Beenakker Beenakker (1991). We define $P_{N}(\{n_{i}\})$ to indicate the probability of the fermion configuration $\{n_{i}=0,1\}$ in the nanowire with the total particle number $N$. Since moving in and out by the Cooper pairs does not cost any energy in a SC, $N=\sum_{i}n_{i}+2N_{c}$, where $N_{c}$ is the additional number of the Cooper pairs away from the charge neutral point as the gate voltage of the nanowire is zero. To simplify the problem, we assume that the charging energy $E_{c}$ is large enough so that only the two lowest energy levels $U(N)$ and $U(N-1)$ of the electrostatic energy

[TABLE]

are involved with the energy levels of the other electron numbers being too high to be included, where $n_{g}$ is proportional to the gate voltage of the nanowire. As the gate voltage is zero ($n_{g}=0$), at the charge neutral point $N=0$ the electrostatic energy reaches its minimum. More energy levels are straightforward to include theoretically, but the results become nontransparent and complicated with no new qualitative insight. We further assume that one-electron transfer process at the left and right junctions are incoherent and independent (the conductance measured in Ref. Albrecht et al. (2016), which is less than $e^{2}/h$, indicates the process is incoherent) and the tunneling rates $\Gamma_{p}^{l,r}$ and $\Lambda_{p}^{l,r}$ are much less than temperature ($\Gamma_{p}^{l,r},\Lambda_{p}^{l,r}\ll T$), where labels $\Gamma_{p}^{l,r}$ and $\Lambda_{p}^{l,r}$ indicate the tunneling rates of the left($l$) and right($r$) for quasiparticle with energy $E_{p}$ and quasihole with energy $-E_{p}$ respectively (see section III.1). The current flowing from the left lead to the nanowire is given by (see Fig. 1 for the actual tunneling processes)

[TABLE]

where $\epsilon^{l}_{p}=E_{p}+\Delta U+\eta eV$, $\tilde{\epsilon}^{l}_{p}=-E_{p}+\Delta U+\eta eV$, the Fermi-Dirac distribution $f(\epsilon)=1/(1+e^{\beta\epsilon})$ and the difference between the electrostatic energies with the different particle numbers $\Delta U=U(N)-U(N-1)$. The first term describes that an electron with energy $\epsilon^{l}_{p}$ at the left lead moves to the energy level $E_{p}$ of the nanowire and its electrostatic energy is changed to $U(N)$ from $U(N-1)$ (see Fig. 1 (a)); the second term describes that an electron in the energy level $E_{p}$ of the nanowire moves to the energy level $\epsilon^{l}_{p}$ at the left lead and the electrostatic energy of the nanowire changes to $U(N-1)$ from $U(N)$ (see Fig. 1 (b)). Since the nanowire is superconducting, the last two terms, which are not included in Beenakker’s paper [Beenakker, 1991], represent the movement of a hole. The third term represents the hole ($+2e$ Cooper pair) with energy $\tilde{\epsilon}^{l}_{p}$ moving from the left lead to the energy level $-E_{p}$ of the SC nanowire (see Fig. 1 (c)); likewise, the fourth term represents the hole ($+2e$ Cooper pair) moving from the energy level $-E_{p}$ of the SC nanowire to the energy level $\tilde{\epsilon}^{l}_{p}$ of the left lead (see Fig. 1 (d)).

Since the nanowire is in a stationary state, the probability of the electron configuration should be time-independent. That is, since the flow in and out of the two leads are equal, the stationary equations for particle number $N$ and $N-1$ are given by

[TABLE]

[TABLE]

where $\epsilon_{p}^{r}=E_{p}+\Delta U-(1-\eta)eV$, and $\tilde{\epsilon}_{p}^{r}=-E_{p}+\Delta U-(1-\eta)eV$. Signs “+/-” in the front of the summations indicate an electron moving in and out in the nanowire respectively. At the applied voltage $V=0$, $P_{N}(\{n_{i}\})$ reaches equilibrium described by detailed balance as

[TABLE]

where $\epsilon_{p}=E_{p}+\Delta U$ and $\tilde{\epsilon}_{p}=-E_{p}+\Delta U$, since in the equilibrium the probability distribution is described by the Gibbs distribution

[TABLE]

where $Z$ is the partition function. We note that in the presence of Cooper pairs the total electron number $N$ is not always $\sum_{i}n_{i}$. When a small applied voltage $V$ is turned on, the probability distribution in the nanowire can be expanded to the linear order around equilibrium

[TABLE]

where $\Phi_{N}(\{n_{i}\})$ is an unknown functional to be obtained by solving the transport master equations to be described below.

By using the two equations above, we expand the current $I$ (2) to the linear order of $V$

[TABLE]

Again using the identity of the Fermi-Dirac distribution $1-f(\epsilon)=f(\epsilon)e^{\beta\epsilon}$ and Eq. 5, we have

[TABLE]

To solve for the conductance $dI/dV$, we can keep the linear terms in $V$ in the stationary equations (Eqs. 3 and 4) by using Eq. 7

[TABLE]

If the nanowire has $n_{\rm{max}}$ energy levels, then a set of $2^{n_{\rm{max}}}$ stationary equations determine $\Phi_{N}$ and $\Phi_{N-1}$. One of the stationary equations is redundant due to the charge conservation

[TABLE]

Although the total number of the $\Phi_{N}$ and $\Phi_{N-1}$ is $2^{n_{\rm{max}}}$, one of $\Phi_{N}$’s can be set to zero with no loss of generality and $\eta$ can be neglected. There are then $2^{n_{\rm{max}}}-1$ variables to be determined, where $n_{\rm{max}}$ is the number of nanowire energy levels participating in the transport process. For an arbitrary $n_{\rm{max}}$, it is difficult to analytically solve for $dI/dV$ except perhaps under some special conditions. Later, in computing the conductance of the superconducting nanowire, we will numerically solve Eqs. 10 and 11 for generic cases (e.g. the presence of Andreev bound states) and use the analytical solution for special cases, which are discussed in the next subsection.

The physical parameters affecting the conductance are finite temperature $T$, all of the energy levels $E_{p}$’s less than and close to temperature, the charging energy $E_{c}$, and the tunneling rates $\Gamma_{p}^{l,r},\Lambda_{p}^{l,r}$. In the next section, we consider the realistic 1D SC proximitized semiconductor model and use its physical parameters to compute the Coulomb blockaded conductance.

II.1 The analytic solution for fixed tunneling ratio

Obtaining the analytic expression for the current (9) for an arbitrary number of energy levels $p$ is complicated. In the special situation when the tunneling ratios are fixed $\{\Gamma_{p}^{l}/\Gamma_{p}^{r}=\Lambda_{p}^{l}/\Lambda_{p}^{r}=\rm{const}\}$ for all energy level $p$, the stationary current conditions (Eqs. 3 and 4) can be analytically solved, and the expression for the current can be explicitly written down. If the two ends of the nanowire are almost identical and localized states are absent, then $\{\Gamma_{p}^{l}=\Gamma_{p}^{r}$, $\Lambda_{p}^{l}=\Lambda_{p}^{r}\}$ fulfills this special condition. First, we assume that these stationary equations still hold as the summation over $p$ is removed (later, we will go back to check if this assumption is valid.)

[TABLE]

[TABLE]

We expand the lead voltage $V$ to the linear order by using Eq. 7

[TABLE]

After the simplification

[TABLE]

we use the condition $\{\Gamma_{p}^{l}/\Gamma_{p}^{r}=\Lambda_{p}^{l}/\Lambda_{p}^{r}=\rm{const}\}$ and let $\alpha\equiv\frac{\Gamma^{r}_{p}}{\Gamma_{p}^{l}+\Gamma_{p}^{r}}=\frac{\Lambda^{r}_{p}}{\Lambda_{p}^{l}+\Lambda_{p}^{r}}$. Eqs. 17 and 18 have the relevant solutions; hence, dropping $\sum_{p}$ in Eqs. 3 and 4 is legitimate. Again using the identity $1-f(\epsilon)=f(\epsilon)e^{\beta\epsilon}$ and Eq. 5, we obtain the expression for the tunneling conductance from Eq. 9

[TABLE]

Since the fermion parity is the only conserved quantity and the particle number is not conserved, the physics is not altered by the transformation $N\rightarrow N+2$ and $n_{g}\rightarrow n_{g}+2$ in Eq. 1; hence, only even and odd $N$’s lead to distinct conductances, which can be explicitly written down for numerical calculations. When $N$ is even, the conductance in the explicit form is given by

[TABLE]

where

[TABLE]

Likewise, when $N$ is odd, the conductance is given by

[TABLE]

We further consider the low temperature limit $T\ll|E_{i>1}|$ (still $T\gg\Gamma_{p}^{l,r},\ \Lambda_{p}^{l,r}$). The conductance can then be written in the simple form, which is consistent with [van Heck et al., 2016]. For even $N$, the low-temperature conductance can be simply written as

[TABLE]

When $N-1$ is even and $N$ odd, we have the following low-temperature limit

[TABLE]

The conductance reaches the maximum as $\tilde{\epsilon}_{1}=0$ or $\epsilon_{1}=0$, and the broadening of the conductance peak is proportional to $T$. Hence, at low temperatures the conductance peaks for even and odd $N$ are located at $n_{g}(N_{e})=N_{e}-1/2-E_{1}/2E_{c}$ and $n_{g}(N_{o})=N_{o}-1/2+E_{1}/2E_{c}$ respectively, where $N_{e}(N_{o})$ indicates even(odd) $N$. The key quantities studied in the experiment Albrecht et al. (2016) are the even and odd $1e$ oscillation peak spacings

[TABLE]

which are the differences between the two closest peaks (see fig. 2). Except for odd and even $N$, the spacings should be independent of $N$. In the experiment, the multiple conductance peaks are measured as the gate voltage of the nanowire varies in a wide region. Since the peaks might fluctuate for different $N$, the even and odd peak spacings are averaged over multiple $N$’s to suppress non-universal effects. All information about the underlying physics is extracted from these conductance peak spacings in Ref. Albrecht et al. (2016).

II.2 The analytic solution for 2-level systems (tunneling from localized states)

Now we consider only 2 energy levels in the nanowire when $N$ is even. This may apply at very low temperatures where only the lowest two energy levels near zero energy are operational in transport processes. The total number of the occupied quasiparticle $n_{1}+n_{2}$ is [math] or $2$. The four stationary equations can be explicitly written as

Eq. 10: $\partial P_{N}(1,1)/\partial t=0$

[TABLE]

Eq. 10: $\partial P_{N}(0,0)/\partial t=0$

[TABLE]

Eq. 11: $\partial P_{N-1}(0,1)/\partial t=0$

[TABLE]

Eq. 11: $\partial P_{N-1}(1,0)/\partial t=0$

[TABLE]

where $\Gamma_{i}=\Gamma^{l}_{i}+\Gamma^{r}_{i}$. By solving 3 of these 4 stationary equations, we obtain the explicit expression for the conductance

[TABLE]

where $\gamma^{l,r}_{i}=\Gamma^{l,r}_{i}/\Gamma_{i}$ and $\lambda^{l,r}_{i}=\Lambda^{l,r}_{i}/\Lambda_{i}$. Similarly, when $N-1$ is even and $N$ is odd, the expression for the conductance is given by

[TABLE]

As shown in the next section, the tunneling rates $\Gamma_{p}^{l}$, $\Gamma_{p}^{r}$, $\Lambda_{p}^{l}$, and $\Lambda_{p}^{r}$ depend on the wavefunctions of the quasiparticles and quasiholes on the ends of the nanowire. Since the tunneling ratios are the same constant $(\Gamma_{p}^{l}/\Gamma_{p}^{r}=\Lambda_{p}^{l}/\Lambda_{p}^{r}=\rm{const})$ in our approximation, the last terms vanish; the conductance equations are reduced to Eq. 19. Because the MBSs are located on the two ends, the conductance peaks stem from the first two terms of Eqs. 30 and 31.

Consider two MBSs, which possess the hybridization energy $E_{1}$ due to wavefunction overlap, located on the left and the middle of the nanowire and assume that the first excited state with energy $E_{2}$ is delocalized. Since $\Gamma_{1}^{r},\Lambda_{1}^{r}\sim 0$ in the absence of the Majorana on the right end, the first and second terms vanish in $dI/dV$. The tunneling ratios are different constants so that the last term survives. However, an electron still can propagate through the MBSs in the presence of the last term of the conductance equations. The physical meaning of the last term mixing with the two-energy-level rates is that after encountering the MBS near the left lead, an electron moves out to the right lead through the extended states in the nanowire.

These transport equations can capture $1e$ tunneling of the conventional (i.e. non-topological) Andreev bound states. Consider two Andreev bound states localized on the two ends separately and the other energy levels much higher than the temperature; the localization leads to $\Gamma_{1}^{r}=\Lambda_{1}^{r}=\Gamma_{2}^{l}=\Lambda_{2}^{l}=0$. The conductance does not vanish in this limit and is given by:

[TABLE]

This result is consistent with the observation in an earlier work Sau et al. (2015) that the so-called “teleportation” phenomenon typically associated with MBSs in the Coulomb blockade regime Fu (2010), can also occur for non-topological Andreev bound states Sau et al. (2015).

To see the clear meaning of Eq. 32 we consider a special case that the system preserves reflection symmetry and the particle and hole tunneling rates are identical. Hence, the notations can be simplified to $E\equiv E_{1}=E_{2}$, $\Gamma\equiv\Gamma_{1}=\Gamma_{2}=\Lambda_{1}=\Lambda_{2}$. The conductance in Eq. 2.32 can be rewritten as

[TABLE]

where $\epsilon=E+\Delta U$ and $\tilde{\epsilon}=-E+\Delta U$. As $\Delta U=0$, the conductance reaches to the maximum with the value

[TABLE]

The location of the peak is independent of the energy of the two states. Furthermore, the conductance peak corresponds to $n_{g}=N+1/2$. We note that this case is one of the extreme limits. It does not imply the conductance peak is always independent of the energy of the states.

III Tunneling and nanowire Hamiltonians

The energy levels $E_{p}$ and the tunneling rates $\Gamma_{p}^{l}$, $\Gamma_{p}^{r}$, $\Lambda_{p}^{l}$, and $\Lambda_{p}^{r}$ are the necessary microscopic input parameters in order to perform the conductance calculations for the SC proximitized semiconductor nanowire. We assume that the superconducting order parameter $\Delta$ in the semiconductor nanowire is proximity induced through contact with a regular metallic superconductor and study the BdG Hamiltonian of the 1D model described by reasonable physical parameters. We note that the energy scales of our results depend on our assumptions about these physical parameters, in particular the SC gap, the Coulomb blockade energy, the nanowire effective mass and g-factor, the chemical potential, the spin-orbit coupling, various hopping amplitudes, and the actual confinement potential in the nanowire. Most, if not all, of these parameters are unknown for the real experimental systems. Therefore, one should not attach special significance to our absolute numbers, particularly the precise temperature scales defining our high and low temperature regimes.

III.1 Tunneling rate

The tunneling rates $\Gamma_{p}^{l}$, $\Gamma_{p}^{r}$, $\Lambda_{p}^{l}$, and $\Lambda_{p}^{r}$ are related to the wavefunction with energy $E_{p}$ near the leads. We first write the BdG Hamiltonian of the superconducting nanowire in real space

[TABLE]

where the annihilation operator including all of the lattice sites is written as $C=(\ldots,c_{x\uparrow},c_{x\downarrow},c_{x\uparrow}^{\dagger},c_{x\downarrow}^{\dagger},\ldots)^{T}$. By diagonalizing $H_{\rm{nano}}$, the Hamiltonian can be rewritten in the diagonal form

[TABLE]

The quasiparticle and quasihole for the energy level $p$ are given by

[TABLE]

The normalization leads to $\sum_{x,\alpha}(|u_{p,x\alpha}|^{2}+|v_{p,x\alpha}|^{2})=1$. The electron creation and annihilation operators are written in terms of quasiparticles and holes

[TABLE]

Now we determine the tunneling rates $\Gamma_{p}^{l}$, $\Gamma_{p}^{r}$, $\Lambda_{p}^{l}$, and $\Lambda_{p}^{r}$ by assuming $kT,\ \Delta E\gg\Gamma,\ \Lambda$, where $\Delta E$ is the energy level separation. According to Fermi’s golden rule, the tunneling rates are proportional to

[TABLE]

where $|f\rangle$ and $|i\rangle$ are the initial and final states respectively and $H_{t}$ is the tunneling part of the Hamiltonian. First, consider the tunneling between the left lead and the superconducting nanowire; the tunneling Hamiltonian can be written as

[TABLE]

where $c_{1,\alpha},c_{1,\alpha}^{\dagger}$ are the annihilation and creation operators located at the left end of the nanowire. As the electron moves from the left lead to the energy level $p$ in the nanowire, the initial and final states are given by

[TABLE]

where $|O_{L}\rangle$ and $|BCS\rangle$ are the normalized wavefunction in the left lead and the superconducting nanowire respectively and $L^{\dagger}_{\epsilon_{p}}$ is the electron creation operator with energy $\epsilon_{p}$. Hence,

[TABLE]

We can obtain the tunneling rate

[TABLE]

where $t_{p}^{\prime}$ absorbs all of the constants from Fermi’s golden rule. Similarly, for an electron moving from the nanowire to the left end the tunneling rate is identical. On the other hand, for a hole movement, by following a similar derivation, we have the tunneling rate

[TABLE]

We assume the undetermined tunneling constants $t_{p}^{\prime}=w^{\prime}_{p}=1$ ($\Gamma_{p}^{l,r},\Lambda_{p}^{l,r}\ll T$ still holds) for convenience leading to

[TABLE]

where $L$ is the length of the nanowire (We intentionally do not provide the unit of $\Gamma_{p}^{l}$ and $\Lambda_{p}^{l}$ since the details of the experimental setups are unknown and non-universal. Our focus is on the OCPS, which is universal and is not necessarily determined by the exact value of the conductance).

III.2 The Hamiltonian for 1D superconducting nanowire

The 1D SC proximitized semiconductor nanowire with spin-orbit coupling in the presence of a field-induced Zeeman spin splitting Oreg et al. (2010); Lutchyn et al. (2010) can be described in momentum space as

[TABLE]

Using the open boundary condition, the Hamiltonian can be written in the following form suitable for numerical calculations

[TABLE]

where $C_{i}=(c_{\uparrow i},c_{\downarrow i},c_{\uparrow i}^{\dagger},c_{\downarrow i}^{\dagger})$. In the following calculations (unless specified otherwise), our choice of the physical parameters is based on Liu et al. (2017) with slight differences: hopping strength $t=6$meV, spin-orbit coupling $\alpha=1.2$meV, superconducting order parameter $\Delta=0.9$meV, the chemical potential $\mu=0.2$meV, and the Zeeman splitting $V_{z}=1.2B$meV, $V_{y}=0$, where $B$ is the magnetic field in unit of Tesla Zhang et al. (2016). We consider the number of lattice sites to be $L=80$ (we vary $L$ later in presenting our results) in the unit of the lattice constant $a=10$nm so the length of the nanowire is 0.8 $\mu$m and show the wire spectrum in fig. 3. These parameters are representative for the experimental system used in Ref. Albrecht et al. (2016), and changing these parameters to other reasonable values for currently used semiconductor-superconductor hybrid systems does not change any of our qualitative conclusions. We have made no attempt to obtain quantitative agreement with experiments since our goal here is to understand the findings of Ref. Albrecht et al. (2016) qualitatively, in particular to see if the intrinsic inconsistency between the length and field dependence of the quoted oscillation amplitude in Ref. Albrecht et al. (2016) can be explained by some mechanism. There are sufficient numbers of unknown parameters (e.g. the chemical potential and the nanowire confinement potential) in the problem (even in the clean limit without invoking any disorder) rendering quantitative comparisons meaningless.

By performing exact-diagonalization of the Hamiltonian, the eigenenergies of the nanowire are given by $\pm\varepsilon_{i}$, where $i$ is a positive integer and $\varepsilon_{i}$ is positive due to particle-hole symmetry in the BdG Hamiltonian. The relation between the quasi-particle and the quasi-hole can be described by

[TABLE]

Since varying the magnetic field through the nanowire is an adiabatic process, the fermion parity of the BCS wavefunction is conserved. In the absence of the magnetic field, we start with the BCS ground state obeying

[TABLE]

and define the energy levels as $E_{p}=\varepsilon_{p}$ for all $p$ for the transport calculation. As the magnetic field increases, the lowest positive energy of the quasiparticle $\varepsilon_{1}$ reaches zero at the TQPT. Due to the Majorana wavefunction overlap in the finite wire, the quasiparticle $a_{1}$, which is the hybridization of the two MBSs, adiabatically evolves to become a quasiparticle with negative energy $-\varepsilon_{1}$. After this energy level crossing at zero energy, the BCS ground state evolves to the first excited state obeying

[TABLE]

Due to the parity fermionic conversion, we still use this state as the basis state (not the ground state) to compute the conductance; the energy levels have to be redefined $E_{1}=-\varepsilon_{1}$, $E_{p}=\varepsilon_{p}$ for all $p>1$. As the magnetic field keeps increasing, after the next zero energy crossing, the BCS state goes back to the original definition in Eq. 49 until the third level crossing and so on. In the following conductance calculation, the basis BCS states are determined by the level crossings as the magnetic field varies. We carry out our numerical calculations following the above prescription and present our results in the next section.

IV Numerical results for conductance

We numerically compute the magnetic-field-dependent conductance of the superconducting nanowire (47) hosting MBSs on the ends after passing through the TQPT. The nanowire described by the Hamiltonian (47) accidentally preserves reflection symmetry and all of the eigenstates are extended, except for the two localized Majoranas at the ends; hence, the tunneling rate on the left and right are equal ($\Gamma_{i}^{r}=\Gamma_{i}^{l},\ \Lambda_{i}^{r}=\Lambda_{i}^{l}$ for all the energy levels). Due to the Majorana hybridization in the finite length nanowire, the tunneling ratio of the hybridized states should be $\Gamma_{1}^{r}=\Gamma_{1}^{l},\ \Lambda_{1}^{r}=\Lambda_{1}^{l}$. With these fixed tunneling ratios, the conductance can be simply computed by using Eqs. 21 and 20. The energy spectrum and the wavefunction on the ends can be obtained by solving the eigenvalue problem of the Hamiltonian (47). For computing the conductance, we still need to know one more physical quantity, which is the magnitude of the electrostatic energy difference between the total numbers $N-1$ and $N$ arising from Coulomb blockade physics

[TABLE]

When $E_{c}$ is less than the nanowire gap ($<\Delta=0.9$meV), the 2$e$ periodicity of the conductance peak, arising from pure Andreev process in the SC wire, dominates. Here our focus is only on the $1e$ periodicity region since we are interested in the MBS physics; therefore, we choose $E_{c}=2$meV greater than the superconducting nanowire gap. (We do not discuss the 2e-periodic physics van Heck et al. (2016) in this work focusing entirely on the large charging energy, $E_{c}>\Delta$, regime since this is presumably not connected with MBSs.) Furthermore, the gate voltage ($V_{g}$) of the nanowire is an experimentally controllable physical parameter, which is proportional to $n_{g}$ since $n_{g}=CV_{g}/e$, where $C$ is the capacitance of the nanowire. In the following, the conductance of the nanowire will be computed in wide regions of normalized gate voltage $(n_{g})$ and Zeeman splitting ($V_{z}$) from the magnetic field ($B$). Furthermore, the energy levels much higher than temperature can be neglected since they do not participate in the transport process. We keep the ten lowest energy levels for the conductance calculations carried out in the current work although it should be feasible to add a few more levels in the numerical work if there is good reason to do so. We believe that 10 levels should suffice for qualitative conclusions as long as the temperature is not too high (which would kill the superconductivity any way). We emphasize that at very low temperatures, the theoretical results for OCPSs disagree qualitatively with the experimental results of Ref. Albrecht et al. (2016) in terms of the magnetic field dependence as described already in the introduction (section I) of this paper.

IV.1 Finite temperature

We numerically compute the conductance of the superconducting nanowire at three different temperatures $T=0.1,\ 0.5,\ 1.0$meV as the magnetic field increases. We emphasize that these absolute temperature scales are somewhat arbitrary being totally dependent on the input parameters of our theory (the corresponding parameters for the actual wires in Ref. Albrecht et al. (2016) are not known), and the three cases should be qualitatively considered as the low, high, and very high temperatures respectively. The low temperature in this manuscript means the temperature is much smaller than the second lowest energy state and still $T\gg\Gamma_{p}^{l,r},\Lambda_{p}^{l,r}$, while the high temperature means at least more than two energy levels are close to or smaller than the temperature scale. (Note that all of our chosen temperatures are on the high side since otherwise the experimental results simply cannot be understood as arising from Majorana physics at all; obviously, if all the energy scales are much lower than our input parameter values, then the same physics could emerge at lower temperatures – we discuss later in this paper the relevance of the energy scale for the experimental nanowires; our choice of model parameters dictate these particular temperature scales given here, which are indeed higher than the quoted temperatures in Ref. Albrecht et al. (2016) perhaps because the experimental induced gap is lower than our theoretical value.) Fig. 4 shows the calculated nanowire conductance for different $V_{z}$ and $n_{g}$. The conductance peaks ($n^{o}_{g}(N),\ n^{e}_{g}(N)$) for even and odd $N$’s can be clearly seen at the low temperature as shown in panel (a) and (d) respectively. At the low temperature, keeping just the lowest energy level is enough to compute the conductance in the topological region (after the bulk gap closes) and OCPSs depend on the the hybridization energy of the MBSs. Thus, at the lowest temperature, the oscillation amplitude as a function of Zeeman splitting is indeed a reflection of the low-lying MBS energy spectrum, thus always manifesting an oscillatory amplitude increasing with increasing magnetic field in contradiction to the experimental finding in Ref. Albrecht et al. (2016). Since at higher temperature the conductance peaks, which are thermally broadened, are difficult to visualize, we present in Fig. 5, the even and odd conductance peak spacings $S_{e}=n_{g}^{e}(N+1)-n_{g}^{o}(N)$ and $S_{o}=n_{g}^{o}(N+1)-n_{g}^{e}(N)$ (see the special case in Eqs. 24 and 25). At the low temperature, as expected and as shown in Fig. 5 (a), the conductance can be described by Eqs. 22 and 23. The spacings should be identical to the two energy levels close to zero energy in Fig. 3. The important feature of these two energy levels is that after the bulk gap closes, the OCPS amplitude becomes larger as the magnetic field increases due to the Majorana hybridization. This is the same as in the simpler theory Das Sarma et al. (2012) without Coulomb blockade and in disagreement with the experimental data in Ref. Albrecht et al. (2016). Although high temperature suppresses the OCPS due to thermal damping, the amplitude of the OCPS still increases as the magnetic field increases in agreement with the low temperature results even if the actual increase is quantitatively damped by temperature. Therefore, the observed amplitude deceasing with increasing magnetic field in Ref. Albrecht et al. (2016) cannot be explained by temperature effects (and the associated multilevel occupancy) alone – some other element of physics is still missing in the theory.

To understand this thermal damping effect of the OCPS, we can consider the limit as the temperature is extremely high. Although the conductance peaks are broadened, the system becomes an effective non-SC metallic Coulomb blockade; the conductance peaks exhibit $1e$ periodicity. Hence, the conductance peak spacings are fixed as $V_{z}$ increases. In this limit, the OCPS completely vanishes. In the other limit, at low temperature the conductance peak is reflected by the lowest energy of state in the nanowire. Then, the OCPS due to the Majorana hybridization is expected. Therefore, these two limits show that high temperature suppresses the MBS-induced OPCS without changing its qualitative behavior.

In the next two sub-sections we consider two possible physical effects, shrinking bulk gap with increasing field and contributions from ordinary Andreev bound states in the trivial regime, to see if we can qualitatively reproduce the puzzling results of Ref. Albrecht et al. (2016). We note that the orbital effect Nijholt and Akhmerov (2016), stemming from parallel magnetic field through the nanowire with a finite cross-sectional diameter, does not explain the decreasing amplitude of the OCPS since the Majorana splitting energy still increases as the magnetic field increases.

IV.2 Shrinking bulk gap

Although the conductance behavior from the standard 1D Majorana Hamiltonian (47) is not consistent with the experimental observation Albrecht et al. (2016), we might suspect that the conductance peak might be affected not only by the Majorana hybridization but also by the size of the bulk SC gap. The reason is that, when the bulk gap collapses, both $S_{e}$ and $S_{o}$ should be unity without any oscillation similar to an effective Coulomb blockaded normal quantum dot with $1e$ periodicity. We mention that indeed in Ref. Albrecht et al. (2016) there is experimental evidence supporting the collapse of the bulk gap near where the oscillations are suppressed with increasing field. There are several ways or mechanisms for the bulk gap to shrink with increasing magnetic field – the most obvious one being that the gap of the parent SC producing the proximity effect itself shrinks with increasing magnetic field. Some possible mechanisms could be: a. reduce the strength of the spin-orbit coupling ($\alpha$) b. reduce the strength of the superconductor order parameter ($\Delta$) as the magnetic field increases. c. change the direction of the magnetic field. d. introduce sub-bands having bulk gap closing after the TQPT point. Since at low temperature the conductance peak spacings are almost identical to the two levels close to zero energy; we only need to compute the conductance spacing at high temperature ($T=0.5$meV unless specified) in the following. Clearly, at low temperatures, the gap closing does not affect the results of the last subsection.

$(\alpha)$

the spin-orbit coupling $\alpha$ is adjusted to $0.5$meV from $1.2$meV. As shown in Fig. 5 (f), the bulk gap near $V_{z}=8$meV is around $0.25$meV compared with $0.4$meV in Fig. 3 with $\alpha=1.2$meV. At low temperature, the OCPS amplitude still increases as the magnetic field increases based on the spectrum (Fig. 5 (d)). Thus, at any temperature, the oscillation amplitude still increases at higher magnetic field in spite of decreasing bulk SC gap.

$(\Delta)$

keeping the strong spin-orbit coupling $\alpha=1.2$meV, we change the superconducting order parameter exhibiting arbitrarily exponential decay $\Delta=\Delta_{o}e^{-V_{z}/5}$meV, where $\Delta_{o}=0.9$meV. (There is no particular significance to this particular form for the bulk gap except that it captures its decay with increasing field in a quantitative manner.) Similarly, as shown in Fig. 5 (g) the oscillation amplitude still increases as the magnetic field increases, although the last peak (near the bulk gap) becomes smaller. Although real $\Delta$ might exhibit a different decay behavior as the magnetic field grows, the qualitative trend of increasing amplitude of the OCPS with increasing magnetic field should remain unchanged.

Direction of the magnetic field

To keep the bulk gap open after the TQPT, the direction of the magnetic field has to be perpendicular to the direction ($\hat{y}$) of the spin-orbit coupling. On the contrary, when the magnetic field is parallel to the the direction of the spin orbit coupling, the bulk nanowire is gapless after the TQPT point. We define an angle $\theta$ between the direction of the magnetic field in the $yz$-plane and the $z$ direction, which is perpendicular to the spin-orbital direction. That is, $V_{z}=V\cos\theta$ and $V_{y}=V\sin\theta$.

For $\theta=0$, the direction of the spin-orbit coupling is perpendicular to the magnetic field, and the bulk spectrum (Fig. 3) is gapped hosting stable MBSs. Now we slightly tilt the direction of the magnetic field. For $\theta=\pi/8$, after the transition point, the bulk gap is open until $V_{z}\sim 3$meV as shown in Fig. 6 (c). The MBSs hybridize and are no longer zero energy modes in the high magnetic field region. When the angle $\theta$ increases, the gap region becomes smaller. At $\theta=\pi/2$, the nanowire is completely gapless after the transition point as shown in Fig. 6 (a).

Now we can consider the peak spacings at higher temperature ($T=0.5$meV) in this tilted magnetic field situation. As shown in Fig. 6 (f-h), the peak spacings exhibit very similar oscillations as the $\theta=0$ case. It is quite surprising that even without the MBSs, the peak spacing can manifest oscillations, Fig. 6 (f), which are almost identical to the nanowire possessing MBSs. Therefore, even if this oscillation of the peak spacings is observed in experiment, it is difficult to definitely conclude the existence of MBSs.

$subbands$

To have the bulk subband gap closing as the magnetic field increases, we consider two copies of the nanowire Hamiltonians (47). The chemical potential $\mu$ in one of the Hamiltonians is adjusted to $5$meV so that the bulk gap starts to decrease in the zero field and close near $V_{z}=5$meV Liu et al. (2017). This model is obviously not able to describe the spectrum beyond $V_{z}=5$meV. In reality, after the subband gap closes, the density of states should be large at zero energy. Hence, the OCPS is suppressed since $S_{\rm{e}}=S_{\rm{o}}=1$ in normal quantum dots.

To see a larger OCPS, the length of the nanowire is reduced to $L=40$ lattice sites from $L=80$. With the remaining parameters unchanged, the energy spectrum is shown in Fig. 6 (e). While the oscillation of the Majorana hybridization energy becomes larger, the bulk gap of the subbands becomes smaller. Consider the peak spacings at two different temperatures $T=0.1,\ 0.3$meV. At the lower temperature ($T=0.1$meV), the spacings are almost identical to the two energy levels close to the zero energy (i.e. the OCPS is a map of the energy spectrum). Although the last oscillation is smaller, the oscillation magnitude is similar to the case without the subbands. At the higher temperature ($T=0.3$meV), the subbands close to zero energy squeeze the OCPSs and the oscillation amplitude decreases, but the conductance is no longer a map of the energy spectrum. Fig. 6 (j) shows that for $V_{z}=3$meV, where the subbands come down at energy level $0.3$meV, the oscillation amplitude starts to become smaller as the magnetic field increases. Although this scenario seems to be partially consistent with the experimental observation, the experimental oscillation amplitude never increases with magnetic field, which disagrees with Fig. 6(j). In Fig. 6(c) the second and third oscillations are still larger than the previous ones; thus, these 1D models cannot completely explain the oscillation observed in the experiment Albrecht et al. (2016) in spite of our incorporating many mechanisms affecting the SC gap. We need something more (i.e. an additional mechanism beyond just gap collapse, finite temperature, and multi-subband occupancy) to make theory and experiment consistent. We include Andreev bound state contribution to the conductance in the next subsection as this new mechanism.

IV.3 Contribution of end Andreev bound states

The results obtained in the previous section for the ideal topological Coulomb blockaded nanowire with MBSs at the end of the nanowire do not appear to produce conductance peak spacings in agreement with the experimental observation Albrecht et al. (2016). Therefore, motivated by the experimental geometry and the conductance results, we consider a case where there are Andreev states at the wire end. Such Andreev bound states are generated by ensuring that the superconductivity vanishes near the ends of the nanowire as shown in Fig. 7 (a). This is a reasonable (perhaps even necessary) consideration since in the experiment Albrecht et al. (2016) the ends of the nanowire, which do not touch the parent superconductor, do not possess proximity-induced superconductivity in all likelihood (see Fig. 7 (a) for the experimental schematic). That is, between the leads and the superconducting nanowire there are small normal metallic regions, which could induce Andreev bound states in the system. There could also be unintentional (and therefore unknown) effective quantum dots inside the nanowire leading to Andreev bound in the system. Low energy Andreev bound states are difficult to distinguish from a split pair of MBSs, which have essentially the same energy spectrum – in fact, one central concern in the Majorana nanowire field is how to distinguish the effects of MBS from those of regular non-topological subgap low-energy Andreev bound states. However, weakly split MBSs are, in principle, distinguishable from Andreev bound states via the localization properties of their wavefunctions. Suppose a low energy Andreev bound state is represented by a creation operator $a_{\epsilon_{1}}^{\dagger}$. Such an Andreev state can be thought of as a pair of weakly coupled MBSs if there exists a phase $\theta$ such that $(e^{i\theta}a_{\epsilon_{1}}+e^{-i\theta}a_{\epsilon_{1}}^{\dagger})/\sqrt{2}$ is localized at one end and $(e^{i\theta}a_{\epsilon_{1}}-e^{-i\theta}a_{\epsilon_{1}}^{\dagger})/\sqrt{2}$ is localized at the other end. In Fig. 8 (b), after the TQPT ($V_{z}\sim 0.92$meV), the lowest energy state satisfies this Majorana criterion. On the other hand, before the transition point the low-energy states, which fail this MBS criterion, are the Andreev bound states in Fig. 8 (b,c). Hence, in this scenario the Andreev bound states and MBSs do not coexist at the same magnetic field since MBSs appear after the bulk gap closes. Thus, we are considering Andreev states and Majorana states on an equal footing because the experiment is carried out by sweeping the magnetic field, and since experimentally the topological transition point (i.e. the critical field at TQPT) is not known, the possibility that the low-field (high-field) behavior in the peak spacings arises from Andreev (Majorana) bound states cannot be ruled out. This is the possibility we are investigating in this subsection. Of course, the criterion of localization is one of degree and therefore this does not precisely define a transition unless the system size goes to infinity where the MBSs can be infinitely far apart. However, for the clean system we consider the critical value of the Zeeman field at the TQPT in agreement with the analytical result $V_{z}\sim\sqrt{\Delta^{2}+\mu^{2}}$ Sau et al. (2010b). Andreev (Majorana) states exist for magnetic field below (above) this critical field. One should think of the transition as a crossover in a finite system, and we need to take into account the possibility that the experimentally observed oscillations arise from a combination of Andreev and Majorana physics as the magnetic field is increased through the critical field near the first peak oscillation.

The localization of the Andreev bound states leads to different tunneling ratios ($\Gamma^{l}_{p}/\Gamma^{r}_{p},\ \Lambda^{l}_{p}/\Lambda^{r}_{p}$) for different energy levels. This means that we must do considerable additional work in order to include the Andreev states in our theory at low magnetic field values below the TQPT. The conductance cannot be simply calculated using Eqs. 21 and 20; we have to directly solve the master equations Eqs. 10 and 11 for the numerical values of $\Phi_{N}$’s and $\Phi_{N+1}$’s and use the conductance equation (9) to obtain the conductance. Since this new calculation including Andreev states is much more involved and is computationally much more demanding than the MBS calculations just using Eqs. 21 and 20, we consider only nine energy levels, instead of ten. To set up the nanowire hosting Andreev bound states for the simulation, we let $\Delta=0$ on the first four left and first three right lattice sites. The reason to choose the different numbers of the left and right lattice sites is to avoid any accidental energy level degeneracy, which leads to the incorrect and nongeneric conductance results from master equation numerical code.

First, we assume that only Andreev bound states are present in the nanowire and MBSs are absent by adjusting the spin-orbit coupling to a low value $\alpha=0.1$meV. The spectrum (Fig. 7 (b)) shows that it is hard to distinguish the MBSs and the bulk states. By comparing with Fig. 3 without Andreev bound states and examining the wavefunctions of the low energy states, the energy levels of the Andreev bound states are the states closest to zero energy in Fig. 7 (b,c) before the bulk gap closing and these localized states become extended after the bulk gap closing point. Unfortunately, the random oscillation of the peak spacings in Fig. 7 (d) appears different from results reported in the experiment. Therefore, only Andreev bound states without any MBS do not appear to explain the experimental results in Ref. Albrecht et al. (2016), although it gives a clue that the presence of Andreev bound states may indeed lead to a situation where peak spacings could sometime reflect a decreasing oscillatory amplitude with increasing magnetic field as observed experimentally. A word of caution is, however, in order here. If chosen appropriately over different samples, it is certainly possible that some of these random Andreev peak spacings could lead to the observed experimental behavior except that other samples would manifest a different behavior. We cannot therefore decisively rule out the possibility that the physics described in Ref. Albrecht et al. (2016) arises from Andreev bound states. We believe that a recent paper Liu et al. (2017) gives a clue to why the Andreev bound states appear to produce random oscillatory behavior in our calculated OCPS results. It is shown Liu et al. (2017) that the Andreev bound states could actually come very close to zero energy accidentally and stay there for finite regimes of magnetic field (“zero-sticking property”) because of the presence of SO coupling and Zeeman splitting. Such accidental “zero-sticking” of Andreev bound states is nonuniversal as a function of the magnetic field and depends on all the parameters of the system. The corresponding OCPS arising from such zero-sticking of Andreev bound states would not manifest any systematic magnetic field dependence and would appear random in a sample to sample measurement. We believe that this zero-sticking property of Andreev bound states is responsible for the OCPS behavior in Fig. 7(d), which may very well be what is being observed in Ref. Albrecht et al. (2016).

Now we introduce MBSs in the nanowire by tuning the spin-orbit coupling back to $\alpha=1.2$meV. Fig. 7 (b) shows the energy of the MBSs to be close to zero and the Majoranas are protected by the large bulk gap. The Andreev bound states are the lowest energy states before the TQPT point. For $V_{z}$ less than $3$meV, as shown in Fig. 7 (e), the peak spacings are similar to the experimental observation since the oscillation becomes smaller as the magnetic field increases. However, in the high magnetic field region the peak spacings are affected by the Majorana hybridization, and the oscillation becomes larger. Of course, the magnetic field at which this change occurs is non-universal, and it is possible that Ref. Albrecht et al. (2016) probes only the low field regime where the Andreev states dominate leading to suppressed oscillations, but then at much higher field eventually oscillations with increasing amplitude as appropriate for MBSs should return. Furthermore, comparing Fig. 7 (c,e), we find the crossing of the conductance peak spacings ($V_{z}\sim 0.6$meV) to occur before the bulk closing point ($V_{z}\sim 0.92$meV). This example shows that the conductance peak spacings cannot be simply described by the lowest energy level (Eqs. 24 and 25). The localized states play a subtle role changing the conductance and the possibility that the Andreev bound states is playing a role cannot be ruled out.

Finally, we add the following ingredient to the physics of the nanowire keeping the Andreev states in the analysis: the bulk gap closing subbands. Since at higher temperature the subbands lead to smaller oscillations at higher magnetic field as shown in Fig. 6 (j), inclusion of all the ingredients (the higher temperature, the Majoranas in high magnetic field, the Andreev bound states in low magnetic field, and the gap closing subbands) leads to a decreasing oscillation amplitude with increasing magnetic field except perhaps at very high magnetic field (which may be outside the experimental regime) as shown in Fig. 8 (d,e). These peak spacings are qualitatively similar to the experimental finding in Ref. Albrecht et al. (2016). Thus, higher temperature and the presence of the closing subband gap as well as the MBSs in the topological region and the Andreev bound states in the trivial region (all of these mechanisms taken together) are currently our best explanation for the experimental observations. Although we do not find a situation where MBSs without Andreev bound states produce results in agreement with experiment, the reverse is not, strictly speaking, true – just Andreev bound states without MBSs give random peak spacings, which, in some situations, may mimic decreasing oscillations with increasing field in a narrow field region. We note that Figs. 7 and 8 show that the same physics applies for long (L=60; Fig. 7) and short (L=40; Fig. 8) wires with the inclusion of Andreev bound states in the calculation always producing decreasing oscillations with increasing magnetic field as observed experimentally. Our findings are therefore generic.

The conductance including all effects at $T=0.2$meV is shown in Fig. 9. The conductance before the TQPT point is much smaller than the one after the TQPT. This result is in agreement with the experiment Albrecht et al. (2016) and has been explained in Ref. van Heck et al. (2016).

The discontinuity appears in the conductance plot in Fig. 9 (a) near $V_{z}=3.6$meV. The conductance discontinuity occurs when the two lowest energy levels are degenerate $E_{1}=E_{2}$ as shown in the spectrum plot Fig. 8 (a). That is, the qausiparticles $a_{\epsilon_{1}}^{\dagger}$ and $a_{\epsilon_{2}}^{\dagger}$ exchange so that the BCS ground state in the odd parity |\rm{BCS_{o}}\rangle=$$a_{\epsilon_{1}}^{\dagger}|\rm{BCS_{e}}\rangle changes sharply as $V_{z}$ passes through the degenerate points; this sharp change leads to the conductance discontinuity. In the reality, since $a_{\epsilon_{1}}^{\dagger}$ and $a_{\epsilon_{2}}^{\dagger}$ usually weakly couple, $E_{1}$ and $E_{2}$ are close but are never identical. The change of the $|\rm{BCS_{o}}\rangle$ is smoothed out and then the conductance discontinuity should not be expected.

To suppress the oscillation as the magnetic field increases (purely in the high-field topological regime where Andreev states cannot play any role) we estimate the lower bound of the temperature to be around $0.2$meV, which corresponds roughly to one fifth of the SC gap, (for the specific parameters used in our study) in Fig. 8 (e) since in high magnetic field the oscillation grows at $T=0.1$meV as shown in Fig. 6 (i). We emphasize that the Andreev bound states can only explain the decrease in the oscillation amplitude in the low-field region (near the first oscillation), and once the MBS oscillations set in, the theory predicts unequivocally that the oscillation amplitude must increase with increasing magnetic field, which is not seen experimentally. We therefore need a mechanism to suppress the oscillations at higher field, and we use temperature as this damping mechanism.

One may wonder what happens at much lower temperature in the presence of both Andreev and Majorana bound states as the magnetic field sweeps through the TQPT. We show our calculated low-temperature results for this general situation in Fig. 8 (f) to be compared with the high-temperature results shown in Fig. 8 (e). As shown in Fig. 8 (f), although near the TQPT the presence of the Andreev bound states does slightly change the oscillation at the low temperature, the oscillation amplitude always grows with increasing magnetic field in contrast to the experimental results of Ref. Albrecht et al. (2016). This is, in fact, expected since at low enough temperatures, the conductance should be a direct map of the low energy spectrum, and the MBS splitting always increases with increasing magnetic field. Thus, even in the presence of Andreev (and Majorana) bound states, the OCPS amplitude increases with increasing magnetic field at low temperatures. The damped oscillation should not be apparent in the experiment unless the electron temperature is pretty high. Of course, as mentioned earlier, what this “pretty high” temperature scale is in absolute units must depend crucially on the experimental parameters which are unknown. All we can claim is that experiments should see distinct behaviors at low and high temperatures with the high-temperature result manifesting decreasing oscillation amplitude with increasing magnetic field.

A serious note of caution is in order about the absolute temperatures used in our simulations, which are much higher than the quoted temperature ($T=50\sim 100mK$) in Ref. Albrecht et al. (2016). We emphasize that our temperature scale is determined entirely by the parameters used in our model. For example, if we use an increased effective mass (i.e. a lower nanowire hopping matrix element), the energy scale goes down, and such an adjustment can induce arbitrary lowering of the temperature in our model. Another (perhaps even more important) point is that the precise chemical potential is not known in the experiment. If the Fermi level is somehow near the bottom of the 1D subbands in the nanowire (not an unlikely scenario given the large density of states near the 1D subband bottom), again the effective energy scale is suppressed lowering the temperatures used in our simulations. Similarly, if the typical subband energy spacing in the nanowire is small (our confinement model for the nanowire corresponds to a hard-wall infinite square well-confinement), that will again make our effective temperatures much smaller. Since these energy scales in the experimental nanowires are unknown, not much significance should be attached to our absolute energy scales. In addition, we assume that the nanowire possesses a hard gap Chang et al. (2015); Cole et al. (2015) after the TQPT point. In reality, the experimental gap is very soft after the TQPT where a zero bias peak emerges, and therefore some (unknown) low energy states are present in the gap. (These low-energy subgap states producing the soft gap most probably arise from disorder in the underlying parent superconductor due to the strong superconductor-semiconductor coupling as discussed in Cole et al. (2016); Takei et al. (2013); Hui et al. (2015), an effect beyond the scope of the current work.) Therefore, even if the temperature is low, these low energy states in the soft gap might be able to suppress the oscillation spacings by acting like effective Andreev states as considered in our theory. Since the nature of these low energy states is unknown, they cannot be included in the theory, but it is quite possible that these low energy fermionic subgap states at high magnetic field act similar to the low magnetic field Andreev subgap states included in our theory, leading to suppressed oscillations with increasing magnetic field even at low temperatures. Since our current conductance numerical program can compute at most 10 energy levels, the numerical result in the manuscript is only for the hard gap. We leave the soft gap simulation, which must include many more energy levels as well as a physical model for the subgap states causing the soft gap, for future work. But, other than suppressing the energy scale at which oscillation amplitude dampens thermally (i.e. by providing many available low energy states already at low temperatures), such low energy states should not change any of our results qualitatively. More experimental and theoretical work will be necessary to decisively settle the question of which subgap states may actually be contributing to the finite temperature conductance causing the increasing oscillation with increasing field. What we have accomplished is to show that a mechanism combining finite temperature and additional subgap states provides a possible explanation for the observed oscillation amplitude decreasing with increasing magnetic field.

Finally, we show our calculated numerical results in Fig. 10 as a function of different number of energy levels included in the theory in order to establish that the 9-level calculation is adequate to obtain convergent results for the parameters used in our theory in the high-temperature regime. We emphasize that the necessary number of levels is obviously a non-universal quantity depending crucially on the input parameters (and temperature) used in our simulation, but the point we make is that multilevel effects must play an important role in understanding the observations in Ref. Albrecht et al. (2016). It is possible (actually likely) that the renormalization by the superconductor makes the higher nanowire energy levels almost degenerate Stanescu , leading to a participation by several energy levels even at much lower temperatures. As discussed already, ‘low’ and ‘high’ temperatures in our calculation are nonuniversal quantities (in absolute terms) as they depend crucially on the unknown microscopic details of the hybrid semiconductor-superconductor structures. The numbers used in our theory are for demonstrative purposes only.

V Length dependence

So far all our results and discussions focused entirely on the magnetic field dependence of $1e$-tunneling OCPSs in the Coulomb blockaded nanowire experiment of Ref. Albrecht et al. (2016), leaving out all considerations of the length dependence. This is appropriate since the measurements are always carried out on a wire of fixed length ($L$) as a function of an applied magnetic field. The experiment is NOT done as a function of length keeping all other system parameters (e.g. magnetic field, chemical potential, superconducting gap) fixed. Thus, any conclusion about an estimated length dependence of splitting is subject to criticism since the other relevant system parameters certainly vary along with the wire length. In fact, in Ref. Albrecht et al. (2016), the length dependence is extracted from the measurement of peak splitting in five different wires with each splitting in each wire measured at different magnetic fields. In fact, not only the magnitude, but even the direction of the magnetic field, is different in some of the wires used to obtain the length dependence. It is expected that these wires of different lengths have different chemical potentials, confinement potentials, superconducting gaps, Coulomb energies, and disorder as well since these are not controllable parameters in the experiment. What was kept common is that the splitting measurement was always done at the first oscillation maxima of the peak spacing, which was often the only oscillation presenting a measurable amplitude since the oscillation amplitude typically decayed rapidly with increasing magnetic field (which is the main point of our theoretical analyses). Some additional limitations of the experiment are that the measured splitting is around $1K$, which is much larger than even the parent SC gap in Al at that magnetic field, for two of the shorter wires ($<0.5\mu m$), and the measured splitting is around $\sim 10mK$ for the longest wire ($L\sim 1.6\mu m$), which is almost an order of magnitude lower than the experimental temperature, leaving only two (or even three) data points in the experiment where the measured spacings can be reasonably construed to be meaningful for extracting an exponential dependence. These two splittings are both $\sim 100mK$ for two wires of $L\sim 1\mu m$. (Here by assuming the low temperature limit, the maximum oscillatory amplitude in the unit of eV is identical to the energy splitting. We refer the readers to Fig. 2 in Ref. Albrecht et al. (2016) for the details.) Whether a decisive $\exp(-L)$ behavior can be meaningfully extracted on the basis of two data points of nearby $L$ values remains unclear. In addition, the fact that the experimental L-dependence in Ref. Albrecht et al. (2016) is extracted by using different samples of different lengths at different magnetic fields (which are sometimes oriented in different directions) casts some doubt in the accuracy of the exponential dependence conclusion reached in Ref. Albrecht et al. (2016). Obviously, many more results are necessary for compellingly establishing an exponential dependence of the splitting in wire length. Nevertheless, the question arises about our theoretical $L$-dependence, specifically for those results in our simulations (including both Andreev and Majorana bound states) which approximately mimic the experimental magnetic field dependence in the OCPSs (i.e. oscillations decreasing with increasing magnetic field).

Although we do not expect any well-defined generic length dependence in our simulated peak spacings since the results include effects of both Andreev and Majorana states as well as finite temperature and multilevel contributions, we show our calculated length dependence in Fig. 11 for the peak spacings obtained operationally in the same way as in Ref. Albrecht et al. (2016). We simply plot the amplitude of the first peak spacing oscillation as a function of $L$ in Fig. 11 making sure that for each value of $L$, the corresponding magnetic field dependence manifests qualitatively similar behavior as observed in Ref. Albrecht et al. (2016) (i.e. amplitude deceasing with increasing field). We keep all other parameters fixed except for the spin splitting, which must change since wires of different lengths manifest the first oscillations at different magnetic field values (precisely as happens in the experiment also). As can be seen in Fig. 11, the theoretically extracted $L$-dependence is nonmonotonic – it shows a strong decrease with increasing length for shorter wires and then a very slow increase with increasing length for longer wires. We have no precise explanation for the nonmonotonic length dependence in Fig. 11 except that this is what we get from the simulations (and the behavior is most certainly non-universal here as there is no reason to expect a universal behavior). Obviously, the decrease in the oscillation for shorter wires, mimicking the findings of Ref. Albrecht et al. (2016), is not a manifestation of the exponential $\exp(-L/\xi)$ behavior of Majorana splitting, since such a behavior should be more prominent for $L\gg\xi$, i.e. for longer wires. In fact, our estimated $\xi$ values for our results are around $100\ (10)$ for high magnetic field (TQPT) values, and since the first oscillation always occurs just below the TQPT in our simulations, it is unclear that the $L$-dependence can have much to do with MBS properties. It is possible, but by no means certain, that the experimentally observed $L$-dependence in Ref. Albrecht et al. (2016) is related to the $L$-dependence in our Fig. 11 for smaller values of $L$, but we emphasize that this short-length behavior has nothing whatsoever to do with topological protection as should be obvious from the larger-$L$ behavior in Fig. 11. While we cannot comment on the significance of the experimentally observed length dependence in Ref. Albrecht et al. (2016), our theoretical length dependence in Fig. 11 is likely to be nongeneric arising from the complex interplay between finite temperature and Andreev bound states and not much significance should be attached to our theoretical $L$-dependence in Fig. 11.

To understand our calculated $L$-dependence in Fig. 11 better, we show in Fig. 12 two sets of numerical results (in panels (a) and (b), respectively) for the theoretical energy splittings and low-lying Andreev bound state wave functions (actually, the squared amplitude, i.e., the probability density) for different values of the wire length $L$ (keeping the quantum dot size always 3 and 4 $\ll L$). The basic picture that emerges is that for small $L$, the first oscillation peak is always close to the TQPT point, and the Andreev bound state is very extended at this low magnetic field, covering much of the wire whereas at higher values of $L$, the first oscillation is associated with Andreev states which are mostly localized near the wire ends. Thus, the physics of the first oscillation is different for small and large values of $L$ with the crossover occurring at some nonuniversal characteristic $L$-value ($\sim 60$ in Fig. 11 for our parameters) which is determined by the system parameters including the size of the quantum dot regimes (see Fig. 7(a)). Typically this characteristic $L$-value defines the approximate crossover for the Andreev bound states being mostly extended over the wire (small-$L$) to being mostly localized near the wire ends (large-$L$). Thus, the decrease in the energy splitting with increasing $L$ in our Fig. 11 has nothing to do with Majorana protection. It arises from the overlap of the Andreev bound states at the two ends rather than the Majorana bound states at the two ends. At this stage, our results, therefore, disagree with an exponential protection conclusion for the nanowire OCPS.

VI Conclusion

Our focus in this work is on the transport properties of the superconducting Coulomb blockaded nanowire in the weak tunneling situation with small applied voltage between the two leads. The generic transport equations in the $1e$ periodicity Coulomb blockaded superconducting region are derived including Coulomb blockade and superconductivity on an equal footing. These transport master equations capture various physical effects, such as multi-band tunneling, finite temperature, tunneling from localized Andreev states, and possible gap closing. The goal is to qualitatively understand in depth the recent experimental report of topological protection in Coulomb blockaded nanowires Albrecht et al. (2016), which is warranted since the direct interpretation of the experiment in terms of the existence of Majorana zero modes shows inconsistency with the theoretical expectation that in a given wire the Majorana oscillations must increase in amplitude with increasing magnetic field (because the SC gap decreases with increasing magnetic field leading to an increasing coherence length). We obtain analytical solutions to our equations to calculate the low temperature conductance arising from the resonant tunneling through MBSs in several simple situations. These low temperature Majorana conductance results cannot explain the experimental observation Albrecht et al. (2016) of decreasing oscillation amplitude with increasing magnetic field. In the presence of Andreev bound states (at low fields) and/or at higher temperatures where many levels contribute to transport, we solve our master equations numerically, and find that the presence of both finite temperature and Andreev bound states are essential in obtaining results which are in qualitative agreement with experiment. We also find that eventually with increasing field, the oscillation amplitude must always increase, but it is possible that the experiment would not work in this high field region because of the complete collapse of the SC gap. We do not find any way to explain the experimental behavior in the low-temperature regime, but we emphasize that the temperature scale in the problem is nonuniversal and is affected by numerous unknown system parameters (effective mass, hopping amplitude, SC gap, chemical potential, SO coupling, confinement, soft gap, subgap states, etc.). We cannot rule out the possibility that Ref. Albrecht et al. (2016) indeed observes our “high-temperature” behavior reflecting the anomalous decrease in the oscillation amplitude with increasing field because of the combined contributions from Andreev states and Majorana states.

When the temperature is less than or close to the lowest energy level of the superconducting nanowire and the next excited energy level is much higher, the OCPSs (as studied experimentally in Ref. Albrecht et al. (2016) and theoretically in Ref. van Heck et al. (2016)) is proportional to the lowest energy level. Therefore, for the simplest Majorana nanowire model Stanescu et al. (2013), the OCPS always grows as the magnetic field increases and the thermal broadening of the conductance peak is proportional to $T$. To explain the damped oscillation (i.e. decreasing with increasing field) observed in the puzzling experiment, we consider several physical mechanisms. Unfortunately, higher temperature, small spin-orbit coupling, and changing of the magnetic field direction always lead to growing oscillations with increasing magnetic field in conflict with the experimental report in Ref. Albrecht et al. (2016). We further introduce the SC gap closing of the subbands after the TQPT point in order to see if gap closing could explain the anomalous behavior. At high temperature in high magnetic field (away from the TQPT point) the oscillation becomes smaller in the presence of gap closing (and in qualitative agreement with experiment); however, the oscillation still grows with increasing magnetic field right after the TQPT point in the low field region. Finally, when we introduce Andreev bound states in the trivial region, the theory becomes much more complicated and the conductance has to be computed by numerically solving the master equations. The calculation time grows exponentially with the number of the included energy levels in the theory. Although the oscillation is suppressed near the first transition point, only the presence of the Andreev bound states by itself does not lead to the consistency of the oscillation behavior with the experiment. In particular, the oscillations in the trivial phase arising purely from Andreev bound states (and no Majorana states at any magnetic field) appear to have random amplitudes although some samples in some limited regimes of magnetic field (for fine-tuned gate voltage values) may very well manifest OCPSs in qualitative agreement with the results in Ref. Albrecht et al. (2016), but this agreement is nongeneric because of the random nature of the theoretical results. Including the high temperature, the Majorana bound states (in the topological regime), the subband gap closing, and the Andreev bound states (in the trivial regime), we find that the OCPS amplitude always decreases as the magnetic field increases. Obviously in this scenario, which is in excellent qualitative agreement with the experimental data presented in Ref. Albrecht et al. (2016), the first oscillations in peak spacing as a function of magnetic field are taking the system from the trivial phase (with Andreev states) to the topological phase (with Majorana states), and as such, no concept of topological protection applies to these first oscillations. Therefore, the Majorana scenario at high field values along with Andreev states at lower magnetic fields can explain the recent puzzling experimental results. We emphasize, however, that in this scenario no significance can be attached to the length dependence of the peak size of the first oscillation since the first oscillation happens below TQPT and is dominated by Andreev, and not Majorana, physics. Although our theoretical length dependence for the peak size of this first oscillation mimics the experimental finding for shorter wire lengths, this behavior is not directly connected with the Majorana splitting. We in fact establish through direct numerical simulations that the decreasing splitting with increasing length (in the shorter wire regime) arises entirely from the physics of Andreev bound states which tend to be more extended through the wire for shorter wires leading to this decrease. For longer wires, where the Andreev states from the two ends are localized and no longer overlap with each other, we find the splitting to be essentially independent of length (actually increasing slowly with increasing length).

Our work shows that it is difficult to understand the important findings in Ref. Albrecht et al. (2016) by invoking only Majorana bound states in the nanowire. At the minimum one also needs finite temperature and Andreev bound states playing crucial roles. These additional mechanisms together can not only explain the puzzling observed decrease in the OCPS amplitude with increasing magnetic field, but also implies that the concept of an exponential protection in length may not apply to the situation since the system is transitioning from the trivial to the topological phase in its oscillatory regime. The only way we see to experimentally establish the topological protection is to use wires of variable length (imposed for example by applying suitable gate voltages) and then study oscillations as a function of length at a fixed magnetic field. These oscillations at a fixed field (in the topological regime) should eventually fall off exponentially as the length increases, signifying topological protection Das Sarma et al. (2012). Our work showing that the first oscillation may arise from the physics of Andreev bound states makes the situation very complex with the conclusion that great care is necessary to establish the purely isolated Majorana topological protection regime experimentally.

If our interpretation of Ref. Albrecht et al. (2016) is qualitatively correct (i.e. both Andreev and Majorana bound states contributing to OCPSs), then there is an immediate significant (and highly encouraging) implication. Since the interpretation most likely (but not absolutely certainly) requires the presence of both Andreev and Majorana bound states in the system, it is clear that going to higher magnetic field should confine the physics entirely to MBSs since the Andreev bound states are operational only near the low magnetic field part involving the first oscillation. This means that the higher magnetic field regime of these nanowires, if accessible, should enable a study of only MBS without any complications arising from Andreev bound states. Unfortunately, this may not be possible in hybrid core-shell epitaxial structures where Al is used as the superconductor since the Al SC gap collapses precisely around where MBS physics becomes operational. But the data of Ref. Albrecht et al. (2016) may be indicating that the true topological Majorana nanowire regime may very well be very close to the parameter values used in Ref. Albrecht et al. (2016) except that one needs to somehow suppress the soft gap behavior so that the SC gap persists to somewhat higher magnetic field. We believe that at magnetic field values higher than that used in Ref. Albrecht et al. (2016) Majorana physics will indeed manifest itself if the SC gap can be kept finite in this high-field MBS regime.

Before concluding, we make two salient comments. First, our conclusion regarding the key role that Andreev bound states may be playing in the experimental results of Ref. Albrecht et al. (2016) is completely consistent with the recent work by Liu et al.Liu et al. (2017), who shows that Andreev bound states have an interesting generic zero-sticking tendency in Majorana nanowires, thus necessarily conflating the physics of Majorana bound states and Andreev bound states. More work is therefore necessary in understanding the role of Andreev bound states vis a vis Majorana properties in nanowire experiments. Second, our results, albeit including many effects, are still obtained within the minimal Majorana nanowire model, and we have no way of ruling out a more complicated model going beyond the minimal model providing a more compelling explanation for the data of Ref. Albrecht et al. (2016) although the interpretation provided in Ref. Albrecht et al. (2016) also uses a minimal model leaving out all the complications (e.g. Andreev bound states, finite temperature, several energy levels, gap collapse, etc.) considered in our work.

VII Acknowledgements

The authors are indebted to J. Alicea, W. S. Cole, A. P. Higginbotham, C.-X. Liu, F. Kuemmeth, R. M. Lutchyn, F. Setiawan, and T. D. Stanescu for discussions. This work is supported by Microsoft Q and LPS-MPO-CMTC.