Zeeman and spin-orbit effects in the Andreev spectra of nanowire junctions

TL;DR

This paper investigates how Zeeman and spin-orbit interactions influence the energy spectra and electromagnetic responses of Andreev bound states in nanowire-based Josephson junctions, providing insights for spectroscopic experiments.

Contribution

It introduces a model accounting for magnetic field penetration in proximitized nanowires, revealing effects on Andreev spectra and fermion parity switches.

Findings

Zeeman and spin-orbit effects significantly alter Andreev spectra.

Magnetic field penetration modifies the energy levels and transition strengths.

Results aid in interpreting microwave and tunneling spectroscopy data.

Abstract

We study the energy spectrum and the electromagnetic response of Andreev bound states in short Josephson junctions made of semiconducting nanowires. We focus on the joint effect of Zeeman and spin-orbit coupling on the Andreev level spectra. Our model incorporates the penetration of the magnetic field in the proximitized wires, which substantially modifies the spectra. We pay special attention to the occurrence of fermion parity switches at increasing values of the field and to the magnetic field dependence of the absorption strength of microwave-induced transitions. Our calculations can be used to extract quantitative information from microwave and tunneling spectroscopy experiments, such as the recently reported measurements in Van Woerkom et al., arXiv:1609.00333.

Equations

$$H$$

$$-\frac{1}{2}g{\mu }_{B}B{\sigma }_x-{\Delta }_0{\tau }_x{e}^{-i\varphi \mathrm{sgn}(x){\tau }_z/2}+V\delta (x){\tau }_z]\psi (x).$$

$$\psi (0^{+})$$

$${\partial }_x\psi (0^{+})-{\partial }_x\psi (0^{-})$$

$$\psi (x)=e^{-i{k}_{F}x}{\psi }_L(x)+e^{i{k}_{F}x}{\psi }_R(x).$$

$$H\approx \frac{1}{2}$$

$$-\frac{1}{2}g{\mu }_{B}B{\sigma }_x-{\Delta }_0{\tau }_x{e}^{-i\varphi \mathrm{sgn}(x){\tau }_z/2}]\Psi (x).$$

$$\theta =\arccos \frac{\alpha k_F}{[{\alpha }^2{k}_F^2+(\frac{1}{2}g{\mu }_{B}B{)}^2{]}^{1/2}}.$$

$${\Psi }_S(x)=S\Psi (x),S=\exp [-i(\theta /2){\tau }_z{s}_z{\sigma }_y]$$

$$H$$

$$-{\Delta }_0(\cos \theta {\tau }_x+\sin \theta {\tau }_y{s}_z{\sigma }_y)e^{-i\varphi \mathrm{sgn}(x){\tau }_z/2}]{\Psi }_S(x).$$

$$\Psi (0^{+})=T\Psi (0^{-}),$$

$$T=1-i(V/v_F)s_z+(V/v_F)s_y.$$

$$\tau =\frac{1}{1+(V/v_F{)}^2}.$$

$$T^{†}=s_z{T}^{-1}{s}_z,$$

$$T_S$$

$$-(V/v_F)\sin \theta {\tau }_z{s}_z{\sigma }_y.$$

$$\psi (x)=e^{-im\alpha x(1+{\sigma }_z)}{\psi }_L(x)+e^{im\alpha x(1-{\sigma }_z)}{\psi }_R(x).$$

$$H\approx \frac{1}{2}\int dx$$

$$-\mu {\tau }_z-{\Delta }_0{\tau }_x{e}^{-i\varphi \mathrm{sgn}(x){\tau }_z/2}]\Psi (x).$$

$$\left[-iv{\tau }_z{s}_z{\partial }_x+O_N-{\Delta }_0{\tau }_x{e}^{-i\varphi \mathrm{sgn}(x){\tau }_z/2}\right]\Phi =E\Phi ,$$

$$\Psi (x)=\sum _n{\Gamma }_n{\Phi }_n(x)+{\Gamma }_n^{†}[\mathcal{P}{\Phi }_n(x)].$$

$$H=\sum _n{E}_n\left({\Gamma }_n^{†}{\Gamma }_n-\frac{1}{2}\right).$$

$$\Phi (x)=e^{i\varphi \mathrm{sgn}(x){\tau }_z/4}\tilde{\Phi }(x).$$

$$\tilde{\Phi }(0^{+})=e^{-i\varphi {\tau }_z/2}T\tilde{\Phi }(0^{-}).$$

$$[E-H_{\text{BdG}}({\partial }_x)]G(x,E)=iv{\tau }_z{s}_z\delta (x),$$

$$H_{\text{BdG}}({\partial }_x)=-iv{\tau }_z{s}_z{\partial }_x+O_N-{\Delta }_0{\tau }_x.$$

$$\tilde{\Phi }(x)=G(x,E)M\tilde{\Phi }(0^{-}),$$

$$M=\left(e^{-i\varphi {\tau }_z/2}T-1\right).$$

$$G(x,E)=-{\int }_{-\infty }^{\infty }\frac{vdq}{2\pi i}\frac{e^{iqx}}{E-H_{\text{BdG}}(q)}{\tau }_z{s}_z,$$

$$\det \left[1-G(0^{-},E)M\right]=0.$$

$$T=1-i\sqrt{\frac{1-\tau }{\tau }}{s}_z+\sqrt{\frac{1-\tau }{\tau }}{s}_y$$

$$E-H_{\text{BdG}}(q)$$

Full text

I Introduction

The Josephson current flowing across a weak link between two superconductors is mediated by Andreev bound states Golubov et al. (2004), sub-gap states localized at the position of the weak link. Recent years have witnessed the direct observation of Andreev bound states in different types of weak links Pillet et al. (2010); Deacon et al. (2010); Zgirski et al. (2011); Bretheau et al. (2013a, b); Chang et al. (2013); Lee et al. (2014), via either tunneling or microwave spectroscopy, as well as their coherent manipulation Janvier et al. (2015). Aside from increasing our understanding of mesoscopic superconductivity, these results pave the way to the realization of novel types of qubits Zazunov et al. (2003); Chtchelkatchev and Nazarov (2003); Padurariu and Nazarov (2010, 2012) and superconducting circuits. Among these results, of particular interest is the very recent microwave detection of Andreev bound states in InAs/Al nanowires van Woerkom et al. . Such hybrid semiconducting/superconducting systems are under intense investigation Mourik et al. (2012); Higginbotham et al. (2015); Albrecht et al. (2016); Deng et al. (2016); Zhang et al. (2017) as platforms for Majorana zero modes Alicea (2012); Beenakker (2013); Leijnse and Flensberg (2012); Stanescu and Tewari (2013) and, eventually, topological quantum computation Nayak et al. (2008); Das Sarma et al. (2015). In these devices, the study of Andreev bound states may be a prelu,de for the study of Majoranas in microwave circuits and the realization of topological qubits.

In view of these exciting applications, the measurement of the magnetic field dependence of the Andreev spectra was among the most interesting aspects of the experiment of Ref. van Woerkom et al. . InAs (or InSb) nanowires are characterized by strong spin-orbit coupling and large g-factors: both are necessary ingredients to reach the topological phase with Majorana bound states which is predicted to occur at high magnetic field and low electron density Lutchyn et al. (2010); Oreg et al. (2010). Spectroscopic studies of nanowire Josephson junctions, even if performed in the topologically trivial phase, can bring quantitative understanding about the interplay of Zeeman and spin-orbit couplings needed for the topological applications. For these purposes, an important merit of such experiments is their high degree of tunability. For instance, in the experiment of Ref. van Woerkom et al. three separate knobs could be tuned to study the behavior of Andreev bound states: the phase difference $\phi$ across the Josephson junction, the transparency of the junction (controlled by a local gate underneath the weak link), and the magnetic field $B$ (which was applied parallel to the wire). Thanks to this high tunability, the measurement of the Andreev spectra can allow one to obtain a great wealth of information about the properties of the device.

In order to understand existing experiments and design future ones, it would be beneficial to have a detailed theory of the Andreev bound states, describing their behavior as the magnetic field and other system parameters are continuously varied. This work aims at providing such theory, focusing on the simple yet experimentally relevant situation of a short, single-channel nanowire junction placed in a magnetic field parallel to the wire (see Fig. 1). Our theory covers all the important regions of the phase diagram, as depicted in Fig. 2. We pay particular attention to the behavior of the Andreev bound states in the topologically trivial phase, since such knowledge may be important to assign experimental data to the correct place of the phase diagram. Aside from the Andreev energy spectrum, we also study in detail the magnetic field dependence of the matrix elements which determine the absorption strength of microwave-induced Andreev transitions.

The study of Andreev bound states in Josephson junctions with broken time-reversal and/or spin-rotation symmetries is a very rich topic of research, covered by a large and diverse body of previous works Altland and Zirnbauer (1997); Bezuglyi et al. (2002); Kwon et al. (2004); Krive et al. (2004, 2005); Buzdin (2005); Dell’Anna et al. (2007); Béri et al. (2008); Michelsen et al. (2008); Reynoso et al. (2008); Fu and Kane (2009); Cheng and Lutchyn (2012); Reynoso et al. (2012); San-Jose et al. (2013); Beenakker et al. (2013); Yokoyama et al. (2013, 2014); van Heck et al. (2014); Mironov et al. (2015); Cayao et al. (2015); Dolcini et al. (2015); Konschelle et al. (2015); Väyrynen et al. (2015); Konschelle et al. (2016); Nesterov et al. (2016); Rasmussen et al. (2016); Peng et al. (2016); Sticlet et al. (2017). In many of the existing studies it is assumed that the time-reversal symmetry breaking is only operative in the weak link, while the effect of the magnetic field in the superconducting parts of the device is disregarded. At the technical level, this means that the effect of time-reversal symmetry breaking is incorporated in the scattering matrix of the junction but not in the description of the superconducting electrodes. In the present context, however, it is crucial to include the effect of the magnetic field on the proximity-induced pairing occurring in the nanowire segments which are in direct contact with the superconductors. Indeed, in experiments aimed at reaching the topological phase, the purpose of the magnetic field is not to influence the local properties of the weak link, but to change the nature of the superconducting pairing induced in the nanowire (whether or not the topological phase is actually reached). The theory of Andreev bound states developed here, therefore, removes the aforementioned assumption and incorporates the effect of the magnetic field in the entire semiconducting nanowire.

Let us summarize the main results presented in this work, and at the same time outline the layout of the paper. In Sec. II, we discuss the nanowire model and the different approximations used in this work. We then derive a determinant equation, Eq. (25), which allows us to solve for the discrete part of the spectrum, i.e. to determine the Andreev bound state energies and their wave functions. The bound state equation (25) makes use of the transfer matrix of the junction, unlike the commonly adopted approach based on the scattering matrix Beenakker (1991, 1992), but akin to previous examples appearing in the literature Shumeiko et al. (1993, 1997); Zazunov et al. (2005, 2014).

In a short single-channel junction, the sub-gap spectrum consists of a doublet of Andreev bound states. In Sec. III, we study the energies of this Andreev doublet by solving the bound state equation both analytically and numerically 111The numerical code accompanying this work is available online at https://github.com/bernardvanheck/andreev_spectra . The Section begins with a review of basic concepts regarding the excitation spectrum of Josephson junctions (Sec. III.1) and of known results about the Andreev bound states at zero magnetic field (Sec. III.2). We then discuss the important features of magnetic field dependence of the Andreev bound state energies at both low and high electron density, and in both the topological and trivial phases (Sec. III.3 and Fig. 6). In particular, we present analytical results for the effective g-factor which determines the linear energy splitting of the Andreev doublet in a small magnetic field, see Sec. III.4 and specifically Eqs. (45)-(47). We show that the g-factor of the Andreev bound states can be strongly suppressed by spin-orbit coupling and/or high electron density. The resulting g-factor can be much smaller than the g-factor of the conduction electrons in the normal state. Equations (45)-(47) also indicate that a measurement of the Andreev bound state g-factor, for instance by means of tunneling spectroscopy, can provide information about the other relevant system parameters. In Sec. III.5 we discuss the appearance of Fermi level crossings in the Andreev spectrum. The presence of Fermi level crossings is significant because it signals a change of the ground state fermion parity of the junction. These “fermion parity switches” can be used as a signature of the topological phase. Namely, if the nanowire is in the topological (trivial) phase, the number of fermion parity switches occurring as the phase difference $\phi$ is advanced by $2\pi$ must be odd (even). In the topological phase, this leads to the well-known $4\pi$ periodicity of the phase dependence of the Andreev bound state energies. The occurrence of fermion parity switches in the trivial phase has so far attracted less attention: here we show that they can appear once the magnetic field crosses a threshold value $B_{\textrm{sw}}$, which depends sensitively on the transparency of the junction and on the strength of the spin-orbit coupling (see Fig. 8).

In Sec. IV we turn our attention to the Josephson current carried by the Andreev bound states, introducing the current operator and briefly discussing the magnetic field dependence of the current-phase relation (Fig. 9). The matrix elements of the current operator between the Andreev bound states determine not only the equilibrium properties of the junction, but also its response to a microwave field. The microwave irradiation of the junction can induce two types of transitions within the Andreev bound state doublet: both are discussed in Sec. V within the linear response regime, appropriate if the applied microwave field is weak. In the first and most notable type of transition, microwaves resonantly excite a Cooper pair from the superconducting condensate to the Andreev doublet. In the second type of transition, instead, microwaves excite one quasiparticle from the first to the second Andreev bound state. The two transitions are distinguished by the parity of the number of quasiparticles involved and so, for brevity, we will refer to them as the “even” and “odd” transition respectively. The even transition is present already at zero magnetic field and, being very bright, is the most easily observed in experiment. The odd transition requires a quasiparticle to be present in the initial state of the junction, either due to a non-equilibrium population or as a consequence of a fermion parity switch. At zero magnetic field, the odd transition is not observable, since in this case the microwave field cannot induce a transition within the degenerate doublet of Andreev levels. However, it may become visible in the presence of both Zeeman and spin-orbit couplings. The magnetic field dependence of the current matrix elements for the even and odd transitions is studied in Secs. V.1 and V.2 respectively. The study reveals that the odd transition, while characterized by a non-zero current matrix element, remains much weaker than the even transition over a wide range of system parameters, including for magnetic fields $B>B_{\textrm{sw}}$. An important consequence of this fact is that, at low temperatures, the absorption spectrum of the junction should exhibit a sudden drop in visibility if the junction is driven through a fermion parity switch by varying the magnetic field or the phase difference (see Fig. 12).

II Model and Andreev bound state equation

Our investigation is based on the well-studied model of a one-dimensional (1D) quantum wire with Rashba spin-orbit coupling, a Zeeman field applied parallel to the wire, and a proximity-induced $s$-wave pairing Lutchyn et al. (2010); Oreg et al. (2010). We consider the Josephson junction geometry shown in Fig. 1. In the limit $L/\xi\to 0$, we can treat the junction as a point-like defect situated at a position $x=0$. Furthermore, provided that the length of the entire nanowire is much larger than $\xi$, we can ignore complications arising from its finite size and treat it as an infinite system in the $x$ direction. The effective BCS Hamiltonian of this system is (we set $\hbar=1$)

[TABLE]

In this Hamiltonian, the field operator $\psi(x)$ is the usual four-component Nambu spinor, $\psi=(\psi_{\uparrow}\,,\,\psi_{\downarrow}\,,\,\psi^{\dagger}_{\downarrow}\,,\,-\psi^{\dagger}_{\uparrow})^{T}$. The two sets of Pauli matrices $\sigma_{x,y,z}$ and $\tau_{x,y,z}$ act in spin and Nambu space, respectively. Furthermore, $m$ is the effective mass, $\alpha$ is the strength of the Rashba spin-orbit coupling, $B$ is the applied magnetic field, $g$ is the effective g-factor, $\mu_{B}$ is the Bohr magneton, $\Delta_{0}$ is the proximity-induced pairing gap, and $\phi$ is the gauge-invariant phase difference across the Josephson junction, and $\mu$ is the chemical potential measured from the middle of the Zeeman gap at $k=0$ (see Fig. 2). Finally, the phenomenological parameter $V$ is the strength of a point-like scatterer which models a potential barrier; later on $V$ will be related to measurable properties of the junction. In practice, the scattering term $V\delta(x)$ enters purely in the boundary condition for $\psi(x)$ at the position of the junction,

[TABLE]

Being a purely one-dimensional effective model, the Hamiltonian in Eq. (1) does not incorporate all the complexity of real devices. For instance, the orbital effect of the magnetic field is not included in our analysis: this is well justified if the cross-section $A$ of the nanowire is small, so that at a given field $B$ the flux piercing the cross-section is much smaller than a flux quantum ($BA\ll h/e$) 222Nanowires currently studied in experiments have a nominal area of about $(100\;\text{nm})^{2}$, but the lowest subbands wave functions are most likely extended over a much smaller area $\sim 1000$ nm2 due to gate confinement and band bending. Thus, they are only sensitive to a portion of the total applied flux. One may estimate that the orbital effect of the magnetic field can not be dominant for $B\lesssim 1$ T, while the experiment reported in Ref. van Woerkom et al. covered a field range smaller than $500$ mT.. The pairing strength $\Delta_{0}$, the spin-orbit strength $\alpha$ and the g-factor $g$ appearing in Eq. (1) should be viewed as phenomenological parameters. The value of $\Delta_{0}$, in particular, strongly depends on the transparency of the semiconductor-superconductor interface.

By virtue of simplicity, the Hamiltonian (1) has become a paradigmatic model in the study of Majorana physics in hybrid semiconducting-superconducting system Stanescu and Tewari (2013). As is well known, it exhibits two distinct topological phases (see Fig. 2). At high chemical potential and/or low magnetic fields, the system is in a trivial superconducting phase with a conventional $2\pi$-periodicity of the ground state energy with respect to the phase difference $\phi$. At low chemical potential, and provided that the condition $\tfrac{1}{2}g\mu_{B}B>(\mu^{2}+\Delta_{0}^{2})^{1/2}$ is satisfied, the system is instead in a topological superconducting phase. In the geometry of Fig. (1), the hallmark of the topological phase is the $4\pi$-periodicity of the ground state energy (for a fixed global fermion parity) with respect to $\phi$, which is associated with the presence of two coupled Majorana zero modes at the junction. The two phases are separated by a critical line $B_{c}(\mu)=2(\mu^{2}+\Delta_{0}^{2})^{1/2}/(g\mu_{B})$ at which the energy gap in the nanowire vanishes, marking a topological phase transition. Note that this criterion is appropriate if the transparency of the semiconductor-superconductor interface is low, so that the coupling between the two materials is weak: in the opposite limit of strong coupling, the critical field $B_{c}$ may depend only weakly on $\mu$ Sticlet et al. (2017). A more complicated phase diagram in the $(\mu,B)$ plane emerges in nanowires with multiple transport channels Stanescu et al. (2011), but we do not consider this situation in this paper.

For both regions of the phase diagram, we are interested in the Andreev level spectrum. That is, we want to find the discrete spectrum of sub-gap states with energy $E<\Delta(B)$, which are localized at the junction via the mechanism of Andreev reflection at the two superconducting interfaces. The energy $\Delta(B)$ is the spectral gap of the continuous part of the spectrum [at zero field, $\Delta(0)\equiv\Delta_{0}$]. In what follows, we will often omit the field argument, $\Delta(B)\equiv\Delta$. In the short junction limit, one expects the number of sub-gap states to be less than or equal to the number of pairs of left/right propagating modes at the Fermi level in the normal state of the nanowire. Thus, the discrete spectrum of the Hamiltonian in Eq. (1) should consist of at most either one or two Andreev levels, depending on whether the system is in the topological or in the trivial phase.

An established way to compute the Andreev level spectrum is the scattering approach Beenakker (1991, 1992); Lesovik and Sadovskyy (2011); van Heck et al. (2014); Sticlet et al. (2017). In this approach, the usual wave function matching problem for bound states is cast in terms of a scattering matrix $S_{N}(E)$ which characterizes the junction, and a second scattering matrix $S_{A}(E)$ which describes the Andreev reflection from the superconducting leads. The two matrices can be combined into a determinant equation for the bound state energies, $\det[1-S_{A}(E)S_{N}(E)]=0$. This approach is particularly advantageous if the following conditions are satisfied. First, the effect of a magnetic field can be neglected in the superconducting leads, so that the gap and the matrix structure of $S_{A}$ are independent of magnetic field. Second, the normal reflection at the superconducting interface can be neglected – this is the so-called Andreev approximation Andreev (1964): it requires $\Delta_{0}\ll E_{F}$ (where $E_{F}$ is the Fermi energy measured from the bottom of the conduction band) and amounts to linearizing the electron dispersion in the normal state. Third, the junction is short, so that the energy dependence of $S_{N}(E)$ can be neglected as long as $E<\Delta_{0}$. When combined, the first and second conditions guarantee that $S_{A}(E)$ is a simple sparse matrix whose energy dependence enters only as a prefactor, $S_{A}(E)\propto\exp[i\arccos(E/\Delta_{0})]$. This circumstance greatly simplifies the solution of the problem, as exemplified by the fact that the determinant equation can be transformed into a finite-dimensional linear eigenvalue problem for $E$ van Heck et al. (2014).

However, as mentioned in the introduction, for our purposes it is crucial to include the effect of magnetic field in the entire system, rather than in the junction alone. The motivation for doing so is two-fold. To begin with, many recent experiments focused on InAs nanowires with epitaxial Al: in this geometry, a parallel field penetrates uniformly the thin aluminum shell. Furthermore, the study of the evolution of the Andreev levels in the different regions of the phase diagram — and in particular across the topological transition — requires that we solve for the Andreev energies taking into account the magnetic field dependence of the spectral gap of the continuous spectrum. Unfortunately, once a magnetic field is included, the Andreev reflection amplitude is not unique anymore but may depend on the initial and final spin and/or orbital states. As a consequence, the calculation of $S_{A}(E)$ becomes non-trivial and strongly dependent on the particular Hamiltonian describing the leads. To overcome this complication, we take a slightly different route and derive an equivalent bound state equation for the Andreev spectrum which generalizes to the more complicated cases in a transparent fashion.

The first step in the derivation is the linearization of the model in Eq. (1), which we perform in two different limits allowing us to cover all relevant regimes of the phase diagram (see Fig. 2). The first limit is that of the high density, $\mu\gg m\alpha^{2},\Delta_{0},g\mu_{B}B$. For such high values of the chemical potential the nanowire will not enter the topological phase in a realistic range of magnetic fields, thus we will use this limit to model a topologically trivial nanowire. The second limit is that of low density, when the Fermi level is inside the Zeeman gap . This is the “helical” regime of the Rashba nanowire: in the normal state, the low-energy theory contains only a pair of counter-propagating modes at finite wave vectors, as well as a gapped pair of modes close to $k=0$. The line $\mu=0$ in the phase diagram coincides with the optimal point at which the critical field is minimal, $B_{c}=\Delta_{0}$ (see Fig. 2), so this limit will allow us to study the Andreev spectrum in the topological phase and around the phase transition. In both limits we will require the Andreev approximation to hold. The Andreev approximation is automatically satisfied in the high density regime, when $\Delta_{0}\ll\mu$. In the low density regime, the chemical potential is low and the Fermi energy is set by the spin-orbit energy, $E_{F}\sim m\alpha^{2}$. Thus, in this limit we must assume the spin-orbit energy to be the dominating energy scale: $m\alpha^{2}\gg\Delta_{0},g\mu_{B}B,\mu$.

In the two following subsections, we carry out the linearization procedure in these two limits, which will then allow us to derive the bound state equation that we seek.

II.1 Linearization for $\mu\gg\Delta_{0},g\mu_{B},m\alpha^{2}$

In the limit of a high chemical potential, we may linearize the normal state dispersion around the Fermi wave vectors $\pm k_{F}=\pm(2m\mu)^{1/2}$. That is, we write the field $\psi(x)$ as a linear superposition of left- and right-moving components,

[TABLE]

Since we are interested in the energy spectrum in a range of energies of order $\Delta$ around the Fermi level, we can assume that $\psi_{L}(x)$ and $\psi_{R}(x)$ vary over length scales much larger than $k_{F}^{-1}$. We may therefore use Eq. (3) in the Hamiltonian (1), organize the resulting expression as an expansion in powers of $k_{F}^{-1}$, and only keep the largest terms. The last step also includes neglecting quickly oscillating terms $\propto\operatorname{e}^{\pm ik_{F}x}$. The result of this procedure can be concisely presented by introducing an eight-component field vector $\Psi=(\psi_{R}\,,\,\psi_{L})^{T}$. In terms of the slowly-varying field $\Psi(x)$, the low-energy Hamiltonian of the nanowire is

[TABLE]

with $v_{F}=k_{F}/m$. Here, we have introduced a new set of Pauli matrices $s_{x,y,z}$ which act in the space of left- and right-movers. Let us now describe the low-energy modes described in this linearized Hamiltonian.

As illustrated in Fig. 3a, around each Fermi point there are two branches in the spectrum of the normal state. The two branches are separated in energy by an amount $2[\alpha^{2}k_{F}^{2}+(\tfrac{1}{2}g\mu_{B}B)^{2}]^{1/2}$ due to the combined effect of spin-orbit and Zeeman coupling. At finite $B$, the spin of each propagating mode is rotated with respect to its orientation at $B=0$ (see arrows in Fig. 3a). The rotation angle is

[TABLE]

and the rotation plane is defined by the Rashba and Zeeman fields. The spin rotation is clockwise (counter-clockwise) for left (right) movers and it can be incorporated in the definition of the field $\Psi$ via a unitary transformation $S$ (see for instance Ref. Nesterov et al. (2016)),

[TABLE]

This rotated basis diagonalizes the homogeneous Hamiltonian of the wire in the normal state. When we express the Hamiltonian in terms of the rotated field $\Psi_{S}$, we find

[TABLE]

This form of the Hamiltonian reveals how the tilting of the modes’ spin affects the pairing. At $B=0$, the $s$-wave pairing does not mix the inner and outer branches of the spectrum since it requires the spins of the two paired electrons to be anti-parallel. At finite field, however, a pairing coupling with strength $\Delta_{0}\sin\theta\approx\Delta_{0}\,(g\mu_{B}B)/(\alpha k_{F})$ emerges between the inner and outer branches, due to the fact that the spin tilts in opposite directions for left and right movers.

To complete the linearization procedure, we must provide the boundary conditions for the field $\Psi$ which are due to the scattering term $V\delta(x)\,\tau_{z}$ in the original model of Eq. (1). The boundary conditions for $\Psi$ can be derived by using Eq. (3) in Eqs. (2) and neglecting terms $\propto\partial_{x}\psi_{L,R}$ with respect to terms $\propto k_{F}$. The resulting boundary conditions can be arranged in the following form,

[TABLE]

with

[TABLE]

The matrix $T$ is, in fact, the transfer matrix associated with the point-like scatterer $V\delta(x)$ in the original model, computed at the Fermi level. The term $\propto s_{y}$ is a backscattering term, while the term $\propto s_{z}$ corresponds to forward scattering. The transmission probability $\tau$ through the junction in the normal state is related to the dimensionless parameter $V/v_{F}$,

[TABLE]

The transfer matrix obeys a “pseudo-unitarity” property,

[TABLE]

which is the equivalent of the most universally known unitarity of the scattering matrix.

For the rotated field $\Psi_{S}$, we must use a rotated transfer matrix $T_{S}=STS^{\dagger}$,

[TABLE]

We see that, similar to the pairing, the backscattering terms are changed when projected to the basis of momentum eigenstates of the nanowire. At zero field, only a single backscattering channel is open for each mode, because scattering preserves spin. At a finite field $B$, two backscattering terms appear, due to the fact that the spin of each left-moving mode has non-zero projection on the spin of both right-moving modes.

II.2 Linearization for $m\alpha^{2}\gg\Delta_{0},g\mu_{B}B,\mu$

When the chemical potential is low, the linearization procedure must take into account that the position of the Fermi points strongly depends on the spin orientation, since the Fermi points are shifted by the Rashba spin-orbit coupling. Specifically, the Fermi points for modes with spin down (up) are situated at $k=2m\alpha$ ($k=0$) for right-movers and at $k=0$ ($k=-2m\alpha$) for left-movers; see Fig. 3b. The linearization of the model therefore begins by writing the field in the following form Klinovaja and Loss (2012),

[TABLE]

Note the presence of the spin-dependent factors in the exponentials, which take into account the dependence of the Fermi points on spin. From here, we proceed as in the previous subsection: assuming that the left- and right- moving fields $\psi_{L}(x)$ and $\psi_{R}(x)$ vary over length scales much larger than $(m\alpha)^{-1}$, we replace Eq. (13) in Eq. (1) with $\mu=0$, and neglect all quickly oscillating terms $\propto\operatorname{e}^{\pm 2im\alpha x}$. Note that, in doing this, it is essential to assume that the spin-orbit energy dominates over the other energy scales. In other words, the spin-orbit length $(m\alpha)^{-1}$ takes the role of the Fermi wavelength as the microscopic length scale of the model.

As a result we obtain the following linearized Hamiltonian of freely propagating modes,

[TABLE]

Here, as in Eq. (4), $\Psi=(\psi_{R},\psi_{L})^{T}$ and the set of Pauli matrices $s_{x,y,z}$ acts in the grading of left- and right-movers.

When written in terms of the components of the vector $\Psi$, the Zeeman term in Eq. (14) is proportional to $\psi^{\dagger}_{R\uparrow}\psi_{L\downarrow}$. It is a mass term which gaps out the two counter-propagating modes with opposite spin crossing at $k=0$ (see Fig. 3). Note that, once this Zeeman gap is formed at the Fermi level, the presence of a scattering impurity may lead to Fano resonances Cayao et al. (2015). The Fano resonances are due to the formation — in the normal state — of quasi-bound states with a characteristic decay length $\alpha/g\mu_{B}B$. The quasi-bound states originate from the inverted part of the parabolic spectrum close to $k=0$, and in principle they can lead to a strong dependence of the transmission of the junction in the normal state on energy Cayao et al. (2015); Nesterov et al. (2016). We may neglect complications associated with their presence by assuming that the junction is short enough so that $L\ll\alpha/g\mu_{B}B$. With this assumption, boundary conditions for $\Psi$ can also be derived as in the previous subsection. We obtain the same transfer matrix $T$ of Eq. (9), except with the velocity $v_{F}$ replaced by $\alpha$; the transmission probability is now $\tau=1/(1+V^{2}/\alpha^{2})$.

II.3 Bogoliubov-de Gennes equations and bound state determinant condition

At this point, in either of the two limits $\mu\gg\Delta_{0}$ and $\mu=0$, our task is reduced to the solution of a system of Bogoliubov-de Gennes (BdG) equations

[TABLE]

for an eight-component Nambu wave function $\Phi(x)$ 333Note that we use the letter $\Phi$ for the BdG wave functions, and $\Psi$ for the corresponding second-quantized fields, to be solved with the boundary condition $\Phi(0^{+})=T\Phi(0^{-})$. The pseudo-unitarity of the transfer matrix $T$, Eq. (11), guarantees that the kinetic energy in the BdG equations remains a Hermitian operator when removing the point $x=0$ from its domain. In Eq. (15), $O_{N}=\alpha k_{F}\tau_{z}s_{z}\sigma_{z}-\tfrac{1}{2}g\mu_{B}\sigma_{x}$ for $\mu\gg\Delta_{0},g\mu_{B}B,m\alpha^{2}$ while $O_{N}=-\tfrac{1}{2}g\mu_{B}B\,\tfrac{1}{2}(s_{x}\sigma_{x}-s_{y}\sigma_{y})-\mu\tau_{z}$ for $m\alpha^{2}\gg\Delta_{0},g\mu_{B}B,\mu$. The velocity $v$ is a placeholder for $v_{F}$ in the former case, and for $\alpha$ in the latter.

As is well known De Gennes (1999), the BdG equations are inherently equipped with a particle-hole symmetry represented by an anti-unitary operator $\mathcal{P}$. The particle-hole symmetry dictates that for each solution $\Phi$ of Eq. (15) at energy $E$ there must be an orthogonal solution $\mathcal{P}\Phi$ at energy $-E$. In our case, $\mathcal{P}=\tau_{y}s_{x}\sigma_{y}\,\mathcal{K}$, with $\mathcal{K}$ the complex conjugation operator. The presence of particle-hole symmetry, and the corresponding doubling of the spectrum, is a consequence of the unphysical doubling of the Hilbert space coming from the introduction of Nambu indices; more fundamentally, it is a consequence of the mean-field approximation which allowed us to express the Hamiltonian (1) as a quadratic form of $\psi$ and $\psi^{\dagger}$ Altland and Zirnbauer (1997).

Once the complete spectrum $\{\pm\,E_{n}\}$ of the BdG equations is known, the field operator $\Psi(x)$ can be written in the eigenmode expansion

[TABLE]

Here, $\Gamma_{n}$ and $\Gamma_{n}^{\dagger}$ are Bogoliubov annihilation and creation operators, obeying fermionic anticommutation relations. They diagonalize the Hamiltonian,

[TABLE]

Our goal is to find the bound state solutions of Eq. (15), which have $\left|E\right|<\Delta$. In order to do so, we first bring the BdG equations to a more convenient form by a change of variable,

[TABLE]

The role of this transformation is to make the spatial dependence of the superconducting phase more easily tractable. The wave function $\tilde{\Phi}(x)$ satisfies a modified boundary condition at the origin,

[TABLE]

Next, we define the Green’s function $G(x,E)$ by

[TABLE]

where the operator $H_{\textrm{BdG}}(\partial_{x})$ is the linearized BdG Hamiltonian of the translationally invariant superconducting wire,

[TABLE]

Note that by definition $G(0^{+},E)-G(0^{-},E)=1$. Now, using the boundary condition (19), we may write

[TABLE]

with

[TABLE]

Equation (22) holds for any $x\neq 0$; the wave function is discontinuous at $x=0$. The Green’s function can be computed as

[TABLE]

with $H_{\textrm{BdG}}(q)$ the Fourier transform of Eq. (21). When $\left|E\right|<\Delta(B)$, the poles in the integrand of $G(x,E)$ lie away from the real axis and Eq. (24) can be computed via a contour integral closing on the upper (lower) half of the complex plane for $x>0$ ($x<0$). Requiring that a non-trivial solution $\tilde{\Phi}(0^{-})$ exists, we obtain from Eq. (22) the following determinant equation for the bound state spectrum:

[TABLE]

This bound state equation for the Andreev levels is cast in terms of a transfer matrix $T$ and a Green’s function $G(0^{-},E)$ for the superconducting leads, rather than in terms of scattering matrices. As a consequence of the short junction limit considered in this work, the energy dependence of Eq. (25) is entirely contained in $G(0^{-},E)$, while $T$ is independent of energy. Furthermore, as mentioned at the end of Section II.2, the matrix $T$ contained in Eq. (25) is the same in both linearization limits when expressed in terms of the transmission probability $\tau$:

[TABLE]

Thus, the differences in the subgap spectrum between the two regimes all arise from $G(x,E)$. In the following, we compute $G(x,E)$ for the two regimes of interest. In doing so, we also derive the magnetic field dependence $\Delta(B)$ of the continuum gap.

II.4 Green’s functions and magnetic field dependence of the induced gap

II.4.1 Green’s function for $\mu\gg m\alpha^{2},g\mu_{B}B,\Delta_{0}$

In order to obtain $G(x,E)$, we must first invert the $8\times 8$ matrix

[TABLE]

This task is simplified by the fact that $E-H_{\textrm{BdG}}(q)$ is a real matrix and thus its inverse must also be real. The result is

[TABLE]

Here, $A_{0},A_{1},A_{2}$ are $8\times 8$ matrices which do not depend on $q$. Their detailed expressions are:

[TABLE]

There are four simple poles $\pm q_{0},\pm q_{1}$ appearing on the right side of Eq. (28), given by

[TABLE]

In order to complete the calculation of the Green’s function, we must insert Eq. (28) in Eq. (24) and perform the integral over $q$. Let us choose $q_{0}$ and $q_{1}$ to be the two poles with negative imaginary part. Then, using the residue theorem and some simple algebra, we obtain the following expression for the Green’s function:

[TABLE]

From Eq. (30) we can easily extract the magnetic field dependence of the continuum gap. The gap $\Delta(B)$ is determined by the smallest value of $E$ such that the poles $q_{0,1}$ have zero imaginary part. A few lines of algebra give the following answer:

[TABLE]

The behavior of $\Delta(B)$ is discussed in detail in Fig. 4. Here we only note that $\Delta(B)$ is a smooth function of $B$, and never reaches zero provided that spin-orbit is present (so that $\alpha\neq 0$). These results are true if we assume no suppression of the gap in the parent superconductor which induces the proximity effect in the nanowire. In the case of InAs nanowires with epitaxial Al, this is justified by the smallness of Al g-factor and shell thickness.

II.4.2 Green’s function for $m\alpha^{2}\gg g\mu_{B},\Delta_{0},\mu$

At low chemical potentials, we must repeat the same calculation but starting from the BdG Hamiltonian contained in Eq. (14). We must first invert the matrix

[TABLE]

In this case, we may simplify the calculation by noting the presence of the unitary symmetry $[E-H_{\textrm{BdG}}(q),s_{z}\sigma_{z}]=0$. This symmetry is a consequence of the fact that the inner ($k\approx 0$, $s_{z}\sigma_{z}=1$) and outer ($k\approx\pm 2m\alpha$, $s_{z}\sigma_{z}=-1$) branches of the linearized spectrum are decoupled in the homogeneous wire (although they are coupled by scattering at the junction). Furthermore, the outer branches are not coupled to the magnetic field in the linearized Hamiltonian of Eq. (14), and so for these modes $s_{z}$ and $\sigma_{z}$ are also separately conserved operators. These facts allow use to separate the inverse of Eq. (33) as a sum of two parts,

[TABLE]

This time, the poles $q_{0},q_{1}$ appearing in Eq. (34) are given by

[TABLE]

while the matrices $B_{0},B_{1},B_{2}$ are

[TABLE]

From these expressions, we may compute the Green’s function in this regime:

[TABLE]

We can again extract the magnetic field dependence of the proximity-induced gap looking at the energy dependence of the poles in Eq. (35). In general, the minimal gap is dictated by the competition between that of the inner and outer modes. The spectral gap for the inner modes, which we denote $\Delta^{(k=0)}(B)$, is given by

[TABLE]

The gap of the outer modes is not influenced by the magnetic field to the leading order in the ratio $B^{2}/m\alpha^{2}$, thus in our effective model it is equal to $\Delta_{0}$ at all fields. The spectral gap is thus given by

[TABLE]

At a fixed value of $\mu$, after a slow initial decrease the proximity-induced gap decreases linearly with field and, as already mentioned, closes at $B=B_{c}(\mu)$, at which point the topological transition takes place (see Fig. 4). Increasing $B$ further, the gap $\Delta(B)$ reopens, growing linearly in field until the gap at $k=0$ becomes larger than that at finite momentum. The gap at finite momentum is equal to $\Delta_{0}$, while it is well known that this gap has in fact a weak field dependence: it is quadratically suppressed with increasing $B$ if corrections of the order $(g\mu_{B}B/m\alpha^{2})^{2}$, not included in our approximation, are taken into account. This limitation is inconsequential for our purposes, since we are mainly interested in the Andreev spectrum in the range of magnetic fields for which the relevant gap is the one at $k=0$.

III Properties of the Andreev spectrum

In this Section we discuss in detail the magnetic field and phase dependence of the Andreev bound state energies. We begin with a review of the basic notions underpinning the understanding of the excitation spectrum of a Josephson junction.

III.1 Andreev levels, excitation spectrum, and fermion parity switches

Solving the determinant equation derived in the previous Section, Eq. (25), allows us to determine the subgap spectrum of the BdG equations (15). Since we are dealing with a purely 1D model in the short junction limit, we expect that the subgap spectrum consists of (at most) two distinct Andreev levels. That is, taking into account the doubling of the spectrum enforced by the particle-hole symmetry, the subgap spectrum of the BdG equations consists of (at most) four solutions $\{\pm E_{1},\pm E_{2}\}$. Without loss of generality, we fix a hierarchy $0\leq\left|E_{1}\right|\leq\left|E_{2}\right|\leq\Delta(B)$.

Once the Andreev levels are determined, the many-body Hamiltonian can be expanded as

[TABLE]

where the dots represent the omission of states coming from the continuous part of the spectrum, with energies higher than $\Delta(B)$. Neglecting the presence of these states, we can limit ourselves to considering just four many-body eigenstates: the vacuum state $|V\rangle$, which is annihilated by both $\Gamma_{1}$ and $\Gamma_{2}$; two single-particle states $|1\rangle=\Gamma^{\dagger}_{1}|V\rangle$ and $|2\rangle=\Gamma^{\dagger}_{2}|V\rangle$; and finally the state with a pair of quasiparticles, $|P\rangle=\Gamma^{\dagger}_{1}\Gamma^{\dagger}_{2}|V\rangle$. The fermion parity of the junction, which is a global symmetry of the Hamiltonian, is even in the states $|V\rangle$ and $|P\rangle$, and odd in the states $|1\rangle$ and $|2\rangle$. Up to a common constant, the energies of these four many-body eigenstates are simply related to the Andreev levels $E_{1}$ and $E_{2}$ via Eq. (40), see the table in Fig. 5.

Note that, so far, we have not specified the sign of the energies $E_{1}$ and $E_{2}$ appearing in Eq. (40). In fact, this choice is arbitrary: as can be seen in Fig. 5 the many-body spectrum is invariant under a change of sign of $E_{1}$ and $E_{2}$. This is, again, a consequence of the particle-hole symmetry of the model. Conventionally, one chooses $E_{1}$ and $E_{2}$ to be positive in Eq. (40). In this case, the ground state of the system is identified with the even parity state $|V\rangle$. The states $|1\rangle$, $|2\rangle$ and $|P\rangle$ are excited states with excitation energies $E_{1}$, $E_{2}$ and $E_{1}+E_{2}$ respectively.

Although the initial choice of the sign of $E_{1}$ and $E_{2}$ in Eq. (40) is conventional and does not have measurable consequences, a change in the sign of $E_{1}$ is physical, and it has measurable and important consequences. Such a change in sign can occur as some of the parameters of the system are varied, typically the magnetic field $B$ or the phase $\phi$. To fix the ideas, let us assume that $E_{1}$ is initially positive and that it can be tuned through the point $E_{1}=0$ by changing a parameter — a so-called Fermi level crossing (see green arrow on the left panel of Fig. 5). When $E_{1}=0$, the states $|V\rangle$ and $|1\rangle$ are degenerate in energy: the energy cost to add a quasiparticle to the junction vanishes (see green arrow on the right panel of Fig. 5). Furthermore, when $E_{1}$ becomes negative, the odd-parity state $|1\rangle$ becomes the ground state of the junction. This ground state transition driven by a Fermi level crossing is commonly referred to as a fermion parity switch.

Fermion parity switches can be generically expected in Josephson junctions with broken time-reversal symmetry Beenakker et al. (2013), and can drastically affect the thermodynamic and transport properties of the junction. The Yu-Shiba-Rusinov states associated with magnetic impurities in $s$-wave superconductors Yu (1965); Shiba (1968); Rusinov (1969); Balatsky et al. (2006) provide an early example of this type of phenomenon. A fermion parity switch is also at the basis of the $4\pi$-periodic Josephson effect associated with Majoranas Kitaev (2001); Kwon et al. (2004); Fu and Kane (2009). In this case, the peculiarity is that there is an odd number of fermion-parity switches in a $2\pi$ phase interval, a signature of the presence of a fermion-parity anomaly in the low-energy theory of the junction (only an even number of fermion parity switches in a $2\pi$ phase interval is allowed in a topologically trivial phase). Later in this Section, we will investigate the occurrence of fermion parity switches in the model under study, both in the trivial and topological phases. Before doing so, we provide an overview of the features of the Andreev spectrum of the model, starting from the well-known case in which $B=0$.

III.2 Solution at zero magnetic field

At zero magnetic field, an analytic solution leads to a well-known universal result for the Andreev levels Beenakker (1991); Furusaki and Tsukada (1991); Furusaki (1999). The Andreev levels form a degenerate doublet, $E_{1}=E_{2}\equiv E_{A}$ with

[TABLE]

This result is valid independently on the values of chemical potential $\mu$ and spin-orbit coupling $\alpha$, provided that the Andreev approximation is applicable.

While the solution (41) is already well-known, it is instructive to reproduce this result from Eq. (25). At $B=0$, the Green’s function $G(0^{-},E)$, which can be deduced from Eqs. (31) or (37), takes a particularly simple form:

[TABLE]

with $\beta(E)=\arccos(E/\Delta)$. There are three meaningful facts about the above expression. First, it is valid for both limits $\mu\gg m\alpha^{2},\Delta_{0}$ and $m\alpha^{2}\gg\mu,\Delta_{0}$, so we already see that the solutions of the determinant equation (25) will be common to the two cases. Second, in both limits the right hand side of Eq. (42) is independent of the spin-orbit coupling strength $\alpha$. This is a consequence of the fact that spin-orbit coupling can be removed from the Hamiltonian via a local gauge transformation, and so the Green’s function evaluated at a single point can be made independent of $\alpha$. Thus, the independence of $E_{A}$ on $\alpha$ can also be explained as a consequence of the short junction limit. Third, the right hand side of Eq. (42) is proportional to the unit matrix in the spin grading, which leads to the anticipated double degeneracy of the solutions. Plugging the Green’s function from Eq. (42) and the transfer matrix from Eq. (26) into the determinant equation (25), we obtain the solution (41).

It is also possible to write down the bound state wave functions explicitly. In order to do this, the first step is to solve the system of linear equations $G(0^{-},E_{A})\,M\,\tilde{\Phi}(0^{-})=\tilde{\Phi}(0^{-})$ to find the wave functions at the position $x=0^{-}$. Then, using the knowledge of $G(x,E_{A})$ at arbitrary $x$, one can reconstruct the entire wave function using Eq. (22). Carrying out this procedure, one finds two solutions $\Phi_{1}(x)$ and $\Phi_{2}(x)$, which are written out in detail in Appendix A. In our model, spin along the $z$ direction is a good quantum number at $B=0$, and $\Phi_{1}(x)$ and $\Phi_{2}(x)$ are identical except for the fact that they carry opposite spin. As anticipated in the previous paragraph, the spin-orbit interaction is not effective in separating the two Andreev levels with opposite spins in energy. This would be true even in a model where the spin-orbit interaction takes a more general form and breaks the spin rotation symmetry completely. The degeneracy can not be explained by invoking Kramers’ theorem either: the Kramers partner of $\Phi_{1}(x,\phi)$ is $\Phi_{2}(x,-\phi)$, so that the two wave functions form a true Kramers doublet only at the time-reversal invariant points $\phi=0,\pi$. Rather, the degeneracy of the Andreev levels is a consequence of the short-junction limit. It is removed by spin-orbit coupling if corrections of order $(L/\xi)$ are taken into account Chtchelkatchev and Nazarov (2003); Béri et al. (2008), or even in the short junction limit in the case of a multi-terminal junction van Heck et al. (2014).

III.3 Magnetic field dependence of the spectrum: qualitative features

When a finite magnetic field is present, in general we find that the Andreev level spectrum cannot be found analytically. Thus, away from simple limits, we resort to a numerical search of the roots of the determinant Eq. (25). In total, once one of the two linearization limits is taken, there are four parameters which determine the spectrum: the magnetic field $B$, the phase $\phi$, the transparency of the junction $\tau$, and either the spin-orbit coupling $\alpha$ (when $\mu\gg\Delta_{0},m\alpha^{2},\tfrac{1}{2}g\mu_{B}B$) or the chemical potential $\mu$ (when $m\alpha^{2}\gg\Delta_{0},\mu,\tfrac{1}{2}g\mu_{B}B$). We focus in particular on the field and phase dependence of $E_{1}$ and $E_{2}$, since these are the two parameters which are varied systematically in experiment.

Let us first discuss the simple situation in which spin-orbit coupling is absent, $\alpha=0$. In this case, spin is a good quantum number and the Zeeman interaction is separable from the rest of the Hamiltonian. One simply obtains a linear Zeeman splitting, $E_{1}=E_{A}-\tfrac{1}{2}g\mu_{B}B$ and $E_{2}=E_{A}+\tfrac{1}{2}g\mu_{B}B$, with the same g-factor as that of the continuum states (see inset in the right panel of Fig. 6). Note that by increasing the magnetic field one reaches a field value $B_{\textrm{sw}}(\phi)=2E_{A}(\phi)/g\mu_{B}$ at which $E_{1}$ changes sign: a Fermi level crossing occurs. Because $E_{A}(\phi)$ has a minimum at $\phi=\pi$, this is the value of the phase at which the Fermi level crossing occurs first upon increasing the magnetic field. After this point, i.e. for $B>B_{\textrm{sw}}(\pi)$, a pair of fermion parity switches is nucleated symmetrically around $\phi=\pi$ (see for instance the top right panel of Fig. 6). This behavior is consistent with the fact that, in a topologically trivial phase, the number of fermion parity switches in a $2\pi$ phase interval must be even. While $E_{1}$ decreases with field, the other Andreev level $E_{2}$ increases and merges with the continuum of states with opposite spins at a field value $B_{\textrm{cross}}=2(\Delta_{0}-E_{A})/g\mu_{B}$. This crossing of the Andreev level with the continuum is protected by spin conservation.

The magnetic field dependence of Andreev level spectrum is qualitatively different in the presence of spin-orbit coupling. The typical behavior of the Andreev spectra at fixed $\tau,\phi$ and $\alpha$ is shown in Fig. 6 in both linearization limits. Before entering into the quantitative details of the features of the Andreev levels, let us discuss the important qualitative features.

We begin by discussing the case $\mu\gg m\alpha^{2},\Delta_{0},g\mu_{B}B$, illustrated in the left panel of Fig. 6. For small magnetic fields, the two Andreev levels $E_{1}$ and $E_{2}$ split linearly. The lowest-lying level $E_{1}$ maintains its approximately linear behavior in $B$ up to the occurrence of a Fermi level crossing. Similarly to the zero spin-orbit coupling case discussed earlier, Fermi level crossings first appear at $\phi=\pi$ upon increasing the magnetic field and are then nucleated in pairs around this point. The field $B_{\textrm{sw}}(\pi)$ at which the Fermi level crossing first occurs depends on $\alpha$ and $\tau$: this dependence is investigated in detail later. The energy $E_{2}$ of the second Andreev level increases with $B$, but bends down at $B\gtrsim B_{\textrm{cross}}$, when $E_{2}$ becomes close in energy to the continuum gap $\Delta(B)$, which is decreasing in field. This is due to the fact that, in the presence of both Zeeman and spin-orbit couplings, there are no symmetries in the model which protect the crossing of the Andreev level with the continuum. This avoided crossing between the Andreev level and the continuum leads to a non-monotonic dependence of $E_{2}$ on $B$. Such a non-monotonic dependence is the cause of the suppression in $B$ of the transition frequency $\omega_{\textrm{even}}=E_{1}+E_{2}$ between the two junction states with even parity, a fact which we used to explain the observed absorption spectra of an InAs/Al Josephson junction in Ref. van Woerkom et al. (see also Sec. V).

In the low chemical potential regime, shown in the middle panel of Fig. 6, the two Andreev levels also split linearly for small magnetic fields. However, their behavior at large fields is drastically different from that at high chemical potential, due to the different behavior of the gap $\Delta(B)$. The two Andreev levels merge in rapid sequence with the continuum of states – whose gap is linearly decreasing – right before the topological transition at $B=B_{c}$. In the topological phase at $B>B_{c}$, we find that the subgap spectrum consists of a single pair of Andreev levels $\pm E_{1}$. We may see the energy $E_{1}$ as the result of the coupling between two Majorana zero modes located at the two interfaces of the junction. This notion is accurate in particular for $\tau\ll 1$, when the two interfaces are weakly coupled.

The phase dependence of $E_{1}$ in the topological phase is shown in the bottom right panel of Fig. 6 for two different values of the junction transparency $\tau$. In both cases, and for any $B>B_{c}$, the energy spectrum displays a single Fermi level-crossing at $\phi=\pi$. (This behavior should be contrasted with that of the topologically trivial phase, where, as discussed earlier, Fermi level crossings appear in pairs, see top right panel of Fig. 6). The pinning of the position of the Fermi level crossing at $\phi=\pi$ for $B>B_{c}$ is due to a symmetry of our particular model. Under the combined operation $\mathcal{S}=\sigma_{x}\mathcal{R}$, where $\mathcal{R}$ is the operator of spatial inversion $x\mapsto-x$, the Hamiltonian in Eq. (1) is mapped to itself up to the change $\phi\mapsto-\phi$ Cheng and Lutchyn (2012). This dictates that the Andreev spectrum must be symmetric around $\phi=0$, i.e. that $E_{1}(\phi)=E_{1}(-\phi)$. If, additionally, we recall that the entire spectrum must be $2\pi$ periodic in $\phi$, the only allowed point where $E_{1}$ can vanish is indeed $\phi=\pi$ (this consideration holds in the case that only a single Fermi level crossing is present in a $2\pi$ period.) A Josephson junction with more transport channels or a denser Andreev spectrum may exhibit a higher number of Fermi level crossings Law and Lee (2011), and in a model where there are no constraints coming from spatial inversion the position of the Fermi level crossing may be in general different from $\pi$.

In the rest of this Section, we investigate in more detail the different qualitative features of the Andreev level spectrum described so far.

III.4 Behavior at small field: Zeeman splitting of the Andreev levels

We have seen that in both linearization limits the Andreev levels split starting from infinitesimally small magnetic fields. The linear-in-$B$ splitting can be captured by standard degenerate perturbation theory applied to the zero-field wave functions presented in Appendix A. This procedure is valid as long as $|\tfrac{1}{2}g\mu_{B}B|\ll\Delta_{0}-E_{A}$, so that the discrete Andreev levels are distant from the continuum part of the spectrum. Thus, the results presented in this section are most relevant for $1-\tau\ll 1$ and $\left|\phi-\pi\right|\ll\pi$, i.e. when the energy $E_{A}$ is much lower than the gap $\Delta_{0}$.

It is useful to cast the result of the perturbation calculation in terms of an effective g-factor which is the linear coefficient of the expansion of $E_{1}$ and $E_{2}$ around $B=0$,

[TABLE]

We find that the Andreev level g-factor $g_{A}$ is different from the “bare” value, $g$, which determines the size of the Zeeman gap at $k=0$ in the homogeneous wire, and that $g_{A}$ can depend strongly on the system parameters.

At perfect transmission, $\tau=1$, the zero-field Andreev bound state wave functions are eigenstates of the velocity operator $s_{z}$ [see Eq. (65)]. Therefore, only the part of the Zeeman coupling which mixes co-propagating modes is effective in splitting the Andreev levels (see Appendix B for a discussion). In this case, it is possible to derive an expression for $g_{A}$ which is valid at any value of the ratio $\mu/m\alpha^{2}$, provided that $\max(\mu,m\alpha{}^{2})\gg\Delta_{0}$:

[TABLE]

The equation above can be derived by using a linearization procedure which interpolates between the two limits $\mu\gg m\alpha^{2}$ and $m\alpha^{2}\gg\mu$ used in Sec. II.1 and II.2 respectively; the derivation is contained in Appendix B. When $m\alpha^{2}=0$, Eq. (45) yields $g_{A}=g$ independently on the value of all other parameters. At any finite value of $m\alpha^{2}$ the Andreev bound state g-factor $g_{A}$ is reduced with respect to the bare value, $g$. The suppression is the strongest when $m\alpha^{2}\gg\mu,\Delta_{0}$, in which case Eq. (45) yields $g_{A}\ll g$. A finite spin-orbit coupling also makes $g_{A}$ dependent on the phase difference $\phi$, with a maximum at $\phi=\pi$.

The considerations in the previous paragraph, based on Eq. (45), remain qualitatively valid also for $\tau<1$. In the presence of scattering, the Andreev bound state wave functions are superpositions of states with opposite velocity. In this case, the magnetic field mixes the counter-propagating components (originating from modes close to $k=0$) as well as the co-propagating ones (originating from modes at finite $k$) - see the discussion in Appendix B. We may write the Andreev level g-factor as a sum of two terms, $g_{A}=g_{\rightrightarrows}+g_{\rightleftarrows}$, corresponding to these two different contributions. The co-propagating contribution is given by

[TABLE]

of which Eq. (45) is a special case. Equation (46) is the dominant contribution to the g-factor when $\mu\gg m\alpha^{2},\Delta_{0}$, in which case $\left|g_{\rightrightarrows}\right|\gg\left|g_{\rightleftarrows}\right|$ and $g_{A}\approx g_{\rightrightarrows}$ for any value of the transmission $\tau$.

The counter-propagating contribution $g_{\rightleftarrows}$ becomes relevant in the opposite regime $m\alpha^{2}\gg\mu,\Delta_{0}$. Indeed, in this regime the dominant mixing introduced by a small magnetic field is the one between the counter-propagating modes at $k=0$, which both participate in the formation of the Andreev bound states provided that $\tau<1$. In the limit $\mu/m\alpha^{2}\to 0$ and $\Delta_{0}/m\alpha^{2}\to 0$, the one treated in Sec. II.2, we find

[TABLE]

which is illustrated in Fig. 7. Equation (47) is the leading contribution to the total g-factor $g_{A}=g_{\rightrightarrows}+g_{\rightleftarrows}$ at low chemical potential, except for the vicinities of $\tau=1$ and $\mu=\mu_{0}$, with $\mu_{0}=\Delta_{0}\sqrt{1-\tau}\left|\sin\phi/2\right|$. In these narrow regions of the parameter space, Eq. (47) is vanishing and thus the g-factor is determined by the co-propagating contribution $g_{\rightrightarrows}$, in spite of its smallness. Furthermore, note that $g_{\rightrightarrows}$ and $g_{\rightleftarrows}$ have competing signs when $0<\mu<\mu_{0}$, and so in this region higher-order corrections in the parameter $\mu/m\alpha^{2}$ may be crucial to determine the g-factor (including its overall sign). However, as discussed in Appendix B, the correction to Eq. (47) due to a finite ratio $\mu/m\alpha^{2}$ cannot be easily computed within a linearized spectrum approximation, since such a calculation necessarily involves the electronic state close to the bottom of the parabolic bands of Fig. 2. Finally, we note that Eqs. (46) and (47) agree in predicting a $\sim 1/\mu$ suppression of $g_{A}$ when $\mu\gg\Delta_{0}$.

The value of $g_{A}$ is not directly accessible in microwave absorption spectroscopy, since the microwave transition frequency $\omega_{\textrm{even}}=E_{1}+E_{2}$ is insensitive to the linear splitting in $B$. However, it is observable in tunneling spectroscopy, which can access $E_{1}$ and $E_{2}$ individually. The analysis contained in the above paragraphs suggests that a systematic investigation of $g_{A}$ may be valuable to obtain information about the electron density and the strength of the spin-orbit coupling in the nanowire. This investigation can be carried out at very small values of the field and may be helpful in predicting or understanding the high-field behavior of the system.

III.5 Occurrence and position of Fermi level crossings

Earlier in the text, we have seen that, in the high chemical potential regime $\mu\gg m\alpha^{2},g\mu_{B}B,\Delta_{0}$, Fermi level crossings may occur at a field $B=B_{\textrm{sw}}$ (see the left panel of Fig. 6). In Fig. 8 we study in more detail the dependence of $B_{\textrm{sw}}$, computed at $\phi=\pi$, on spin-orbit coupling strength and transmission.

There are two notable trends. First, the switching field $B_{\textrm{sw}}$ decreases upon increasing the transmission $\tau$ at fixed spin-orbit strength. This is due to the fact that the larger $\tau$ is, the closer to zero is $E_{A}$, and thus a smaller field is required to induce the Fermi level crossing. Second, when increasing the spin-orbit strength at fixed transmission, the field $B_{\textrm{sw}}$ increases. This is due to the suppression of the Andreev level g-factor $g_{A}$ with increasing spin-orbit strength or chemical potential, see Eq. (46), which leads to a slower decrease of $E_{1}$ with $B$. Our numerical results suggest that there is a value $(\alpha k_{F})_{\textrm{max}}$ above which the Fermi level crossings are absent: that is, the curves in Fig. 8 have an asymptote at finite $\alpha k_{F}$ at which $B_{\textrm{sw}}$ diverges. Qualitatively, a strong spin-orbit coupling may prevent the Fermi level crossing to occur because of the level repulsion between the Andreev level $E_{1}$ and the negative image of the rest of the spectrum. Judging from the numerical data shown in Fig. 8, $(\alpha k_{F})_{\textrm{max}}$ depends on the transmission $\tau$, and it grows with increasing $\tau\to 1$. We attribute this behavior to the fact that, in the limit $\tau\to 1$, $E_{A}(\pi)\to 0$: thus, a Fermi level crossing appears already at an infinitesimally small field, and it becomes prohibitive to remove it. Finally, we notice that numerical calculations do not reveal the presence of Fermi level crossings in the opposite regime $m\alpha^{2}\gg\mu,g\mu_{B}B,\Delta_{0}$. We attribute this behavior to the fact that, in this regime, $g_{A}\ll g$. Therefore, the decrease in energy of the Andreev bound states is much slower than that of the continuum states (see for instance the middle panel of Fig. 6), preventing the occurrence of a Fermi level crossing at a field $B<B_{c}$.

In a tunneling spectroscopy experiment, the closing of the excitation gap of the junction at a Fermi level crossing in the regime $\mu\gg\Delta_{0}$ may be naively mistaken for a bulk topological transition. Indeed, a typical magnetic field scale for a fermion parity switch is $B_{\textrm{sw}}\sim 500$ mT van Woerkom et al. , not dissimilar from that of the critical field $B_{c}$ Albrecht et al. (2016). The strong dependence of $B_{\textrm{sw}}$ on $\tau$ (as well as $\phi$), however, should allow to discriminate easily between the two cases.

IV Current operator and the equilibrium current

In this Section, we evaluate the temperature and magnetic field dependence of the equilibrium current. It is known that, in a short Josephson junction not subject to magnetic field, the current is carried almost entirely by the Andreev bound states Beenakker (1991); Furusaki (1999); Heikkilä et al. (2002); Levchenko et al. (2006). This conclusion remains true also in the presence of magnetic field (with or without spin-orbit coupling), as we will now argue following the discussion from Ref. Levchenko et al. (2006). On one hand, the energies of the Andreev bound states vary by an amount $\sim\Delta$ upon varying the phase $\phi$, and thus they provide a finite contribution to the current in the limit $L/\xi\to 0$. On the other hand, the contribution of the continuous spectrum to the current density comes from states within the energy range $\Delta<E<E_{\textrm{Th}}$. Here, $E_{\textrm{Th}}$ is the Thouless energy, i.e. the energy scale associated with the flight time of quasiparticles across the junction; in a short quasi-ballistic junction, the Thouless energy is large, $E_{\textrm{Th}}/\Delta\sim\xi/L\gg 1$. The spectral density of the current delivered by states with energy $\sim E$ scales as $\Delta^{2}/(E_{\textrm{Th}}E)$ for energies in the interval $E_{\textrm{Th}}\gtrsim E\gg\Delta$. It yields a total contribution $\propto(\L/\xi)\,\ln(\xi/L)$ to the current, which vanishes in the limit $L/\xi\to 0$ Levchenko et al. (2006). This argument remains valid even in the presence of a magnetic field or spin-orbit coupling. Therefore, in the following we neglect the contribution of the extended states to the current.

We start by finding the current operator $j(x)$ for the junction, and then evaluate the contribution of the many-body eigenstates $|V\rangle,\,|1\rangle,\,|2\rangle,$ and $|P\rangle$, see Fig. 5. The current operator can be derived from a continuity equation for the electric charge density $\rho$, which for the original model of Eq. (1) is given by the operator

[TABLE]

The continuity equation for $\rho$ can be computed using the equation of motion of the field $\psi(x)$ under the Hamiltonian of Eq. (1). It can be cast in the form

[TABLE]

with $j(x)$ the quasiparticle current operator, which includes a contribution from the spin-orbit coupling,

[TABLE]

and $s(x)$ a charge source (or drain) term due to the presence of the superconducting condensate Blonder et al. (1982),

[TABLE]

At the position of the junction, $x=0$, the source term vanishes since there is no proximity-induced pairing $\Delta_{0}$. Thus, at the junction the equilibrium current can be computed by studying only the quasiparticle contribution coming from $j(x)$. In the superconducting leads, the quasiparticle current is converted into current carried by the condensate over a length $\sim\xi$. Correspondingly, the contribution of the $j(x)$ term to the equilibrium current decays away from $x=0$. The decay is compensated by the source term Blonder et al. (1982) to ensure the current conservation along the wire.

Using Eq. (3) and Eq. (13) we find the current operator projected to low energies [recall that $\Psi=(\psi_{R}\,,\,\psi_{L})^{T}$ encodes the left- and right-moving envelope fields],

[TABLE]

where $v=v_{F}$ in the limit of high chemical potential, and $v=\alpha$ for $m\alpha^{2}\gg\mu,g\mu_{B}B,\Delta_{0}$. Using the pseudounitarity of the transfer matrix, Eq. (11), together with the boundary condition for $\Psi$ at the origin, Eq. (8), one can check that the linearized current operator in Eq. (52) is continuous across the junction, i.e. $j(0^{-})=j(0^{+})$. We will thus evaluate the current at $x=0^{-}$ from now on and omit the position argument.

The current operator can be expanded in the eigenbasis of the linearized Hamiltonian by using Eq. (16):

[TABLE]

Here we have introduced the matrix elements of the current operator between BdG eigenstates,

[TABLE]

The diagonal matrix elements $j_{n,n}$ in Eq. (53) give the dissipationless supercurrent. Including, as already discussed, only the contribution from Andreev bound states to the sum in Eq. (53), we find

[TABLE]

with $\left\langle\cdot\right\rangle$ being the quantum expectation value and $n_{n}=\left\langle\Gamma^{\dagger}_{n}\Gamma_{n}\right\rangle$ the occupation factors for the different quasiparticle states. At thermal equilibrium with temperature $T$, $\left\langle\Gamma^{\dagger}_{n}\Gamma_{n}\right\rangle_{\textrm{eq}}=f(E_{n})$, where $f(E)=[1+\exp(E/k_{B}T)]^{-1}$ is the Fermi-Dirac distribution. At $B=0$, the diagonal matrix elements have a simple analytic expression which can be computed using the wave function in Eq. (65) in Appendix A,

[TABLE]

The result is independent of $\mu$ and $\alpha$, as long as the Andreev approximation is valid, see Sec. III.2. Plugging the expression above into Eq. (56) immediately leads to the known expression for the Josephson current in a single channel weak link

[TABLE]

In this zero field case, the fact that $j_{1,1}=j_{2,2}$ has the consequence that the Josephson current vanishes if the state of junction is one of the two odd parity states. Namely, from Eq. (56) we see that $\left\langle j\right\rangle=0$ if $j_{1,1}=j_{2,2}$ and $n_{1}+n_{2}=1$. This is the so-called “poisoned” state of the junction Zgirski et al. (2011); Olivares et al. (2014), in which one excess quasiparticle can completely block the passage of current. Note that, if the junction has more than one pair of Andreev bound states, a single excess quasiparticle will not block the current completely, as there will be more contributions to the total equilibrium current.

The typical behavior of the equilibrium current-phase relation at finite magnetic field is illustrated in Fig. 9. At small fields, the behavior is not qualitatively different from that of Eq. (58). At low temperatures the current-phase relation exhibit the skewed-sine shape typical of weak links, with the skewness being suppressed with increasing temperatures (see upper panel in Fig. 9). The behavior is more interesting at higher fields, such that fermion parity switches occur as a function of the phase $\phi$, as in the Andreev spectrum in the upper right panel of Fig. 6. In this case, at $T=0$ the current exhibits a discontinuity in correspondence with each fermion parity switch (see bottom panel in Fig. 9). At finite temperatures there is no discontinuity, but a remnant of the fermion parity switches remains in the behavior of the current phase relation close to $\phi=\pi$. Finally, we mention that, as expected, the current model does not exhibit the anomalous Josephson effect (i.e., a finite supercurrent at $\phi=0$). Indeed, for single-channel nanowire Josephson junction, it is known that the latter requires a component of the magnetic field to be aligned with the spin-orbit field Nesterov et al. (2016).

V Microwave absorption

In this section we study the microwave absorption spectrum of a short Josephson junction Bergeret et al. (2010); Kos et al. (2013); Virtanen and Recher (2013); Bretheau et al. (2014); Peng et al. (2016) (for the opposite case of a long junction, see also Refs. Dassonneville et al. (2013); Ferrier et al. (2013); Väyrynen et al. (2015)). The microwave field is modeled as a monochromatic ac voltage drop $V(t)=V_{0}\cos(\omega t)$ across the junction and is minimally coupled to the electronic field $\psi$. This leads to the addition of the following time-dependent term to the Hamiltonian of Eq. (1):

[TABLE]

where $j$ is the current operator evaluated at the junction. We assume that the perturbation $\delta H$ is small, $eV_{0}/\omega\ll 1$. The form of the perturbation $\delta H(t)$ remains valid also after the spectrum linearization, since as we have discussed in the previous Section the current matrix elements at the position of the junction remain well defined and continuous, $j(0^{-})=j(0^{+})\equiv j$. Using standard linear response theory, the expectation value of the current at time $t$ is determined by the response function $\chi(t)=-i\,\theta(t)\,\left\langle[j(t),j(0)]\right\rangle_{\textrm{eq}}$,

[TABLE]

In the frequency domain, the response function $\chi(\omega)$ determines the admittance of the junction, $Y(\omega)=i\chi(\omega)/\omega$. In turn, the real part of the admittance gives the absorption power $W$ of the microwave radiation, $W=\tfrac{1}{2}\,V_{0}^{2}\,\operatorname{Re}Y(\omega)$ with $\omega>0$. Using Eq. (53) to compute the response function, we find

[TABLE]

The first line in Eq. (61) corresponds to transitions where two quasiparticles are created by breaking a Cooper pair and occupy two energy levels with energies $E_{n}$ and $E_{m}$. The second line corresponds to transitions where a single quasiparticle with energy $E_{n}$ is excited into a higher state with energy $E_{m}$. We will refer to these two types of transition respectively as the “even” or “odd” ones, since they are distinguished by the parity of the number of quasiparticles involved. Note that only transitions in which initial and final states have the same fermion parity are allowed. Transitions between the discrete states, which are accounted for in Eq. (61), produce sharp maxima in the frequency dependence of the absorption coefficient. The omitted terms in the admittance, indicated by dots in Eq. (61), involve unbound quasiparticle states and result in an absorption continuum.

We shall consider low frequencies $\omega<2\Delta$, focusing on the transitions between the Andreev bound states. Indeed, at these low frequencies the excitation of Andreev states are the only possible resonant processes (unless the system is close to the critical point separating topological and trivial phases, a case treated in Ref. Tewari et al. (2012)). Transitions between possible above-gap non-equilibrium quasiparticles are very weak and do not result in a sharp absorption line, so we will ignore them. In the case under consideration of a single-channel short junction, with only two Andreev states with energies $E_{1}$ and $E_{2}$, there is only one relevant term in each sum in Eq. (61). These terms correspond to the two allowed transitions depicted in the bottom right panel of Fig. 5: the pair excitation $|V\rangle\to|P\rangle$ with frequency $\omega_{\textrm{even}}=E_{1}+E_{2}$ and the single particle excitation $|1\rangle\to|2\rangle$ with frequency $\omega_{\textrm{odd}}=E_{2}-E_{1}$. As mentioned in the introduction, we call these “even” and “odd” transitions, respectively. The matrix elements $j_{2,\mathcal{P}1}(0^{-})$ and $j_{2,1}(0^{-})$ determine the strengths of these transitions; their dependence on the system parameters is discussed next in detail, first for the even transition and then for the odd one.

V.1 Visibility of the even transition $|V\rangle\to|P\rangle$

We start by discussing the case $B=0$. Again by using the wave functions in Eq. (65) of Appendix A, we find the analytical expression for the relevant current matrix element

[TABLE]

Equation (62) was previously derived in Ref. Kos et al. (2013)using a tunneling Hamiltonian formalism, which is in agreement with our current method based on the transfer matrix. Just like the $B=0$ Andreev energy $E_{A}$, it is independent on the chemical potential $\mu$ and the spin-orbit coupling $\alpha$ and it generalizes to the case of multiple transport channels with different transparencies. Note that $\left|j_{2,\mathcal{P}1}\right|^{2}$ vanishes for $\tau=1$: the absence of scattering at the junction prevents the excitation of the Andreev bound states since in this case the current is a diagonal operator in the eigenbasis of Eq. (53). In the presence of scattering, the Andreev bound states are superpositions of different current eigenstates and microwave-induced transitions become possible Zazunov et al. (2005). Equation (62) has a maximum at $\phi=\pi$, corresponding to the point of greater visibility of the absorption spectral line. The visibility vanishes for small phases. This behavior is in agreement with experiment both in case of nanowire Josephson junctions van Woerkom et al. as well as other types of weak links Bretheau et al. (2013a, 2014).

The dependence of $j_{2,\mathcal{P}1}$ on magnetic field can be determined by finding the wave functions of the Andreev bound states numerically via Eq. (22). We find that a finite magnetic field suppresses the magnitude of the current matrix element, while maintaining its phase dependence qualitatively similar to that described by Eq. (62), see upper panel of Fig. (10). The decrease of $j_{2,\mathcal{P}1}$ with increasing magnetic field is slow in both regimes $m\alpha^{2}\gg\mu_{0},\Delta_{0}$ (see the inset of the bottom panel in Fig. 10) and $\mu\gg\Delta_{0},m\alpha^{2}$ (see Fig. 12). We attribute this decrease to the suppression of the proximity-induced gap $\Delta(B)$ with $B$ (see Fig. 4), which makes the Andreev bound states less tightly confined to the junction and thereby decreases the effective coupling to microwaves, $j_{2,\mathcal{P}1}\sim\Delta(B)$. Finally, at finite fields, the current matrix elements also acquires a weak dependence on the chemical potential, as illustrated in the bottom panel of Fig. (10).

V.2 Visibility of the odd transition $|1\rangle\to|2\rangle$

Without magnetic field, $B=0$, the current matrix element associated with the odd transition (which has anyway zero frequency) vanishes: $j_{2,1}=0$. This is due to the fact that the zero-field Andreev bound states have opposite spin [see Eqs. (63) and (64)], while the perturbation Hamiltonian (59) preserves spin. As the magnetic field is increased from zero, the odd transition may become visible depending on the spin-orbit coupling strength. If spin-orbit is absent (or negligible), the two Andreev bound states would develop an opposite spin polarization in the presence of a Zeeman field: therefore, again due to the spin selection rule, the odd transition would remain forbidden. In the presence of both spin-orbit coupling and magnetic field, however, this spin selection rule is no longer applicable: the two Andreev bound states would have a non-zero spin overlap and one may in general expect a non-vanishing matrix element. Indeed, we determine numerically that $j_{2,1}\neq 0$ at finite $B$ in the presence of spin-orbit coupling. Importantly, even in this case we find that $\left|j_{2,1}\right|/\left|j_{2,\mathcal{P}1}\right|\lesssim 0.1$, see the bottom right panel of Fig. 11. Hence, despite not being forbidden, the dim odd transition may be still much more difficult to observe with respect to the bright even transition. We now discuss the dependence of $j_{2,1}$ on the system parameters.

As in the case of the even transition, the current matrix element $\left|j_{2,1}\right|$ has a maximum when $\phi=\pi$ and vanishes for small phase differences; in what follows, we focus on the peak value. The dependence of $\left|j_{2,1}\right|$ on chemical potential is shown in the top panel of Fig. 11 for different values of $B$. The current matrix elements is non-zero at $\mu=0$, it grows slowly and it reaches a maximum at a small value of $\mu/\Delta_{0}$ before decreasing again. After this point, we find that $\left|j_{2,1}\right|\propto(\Delta_{0}/\mu)^{2}$ when $\mu/\Delta_{0}\gg 1$, as shown in the bottom left panel of Fig. 11. These considerations are valid when $m\alpha^{2}\gg\Delta_{0},\mu$. The suppression of $\left|j_{2,1}\right|$ for $\mu\gg\Delta_{0}$ in this regime matches the numerical results that we obtain for $\mu\gg m\alpha^{2}$, where we find that the matrix element is zero (within numerical precision) at any value of the spin-orbit coupling. Finally, the numerical data indicate that the current matrix elements grows quadratically in $B$ at small fields: $\left|j_{2,1}\right|\propto B^{2}$, see the inset in the top panel in Fig. 11.

The smallness of the current matrix element $j_{2,1}$ has important consequences for the Andreev spectroscopy of the junction at low temperatures, in case the junction undergoes a fermion parity switch. For instance, suppose that the magnetic field is sweeped from a value $B<B_{\textrm{sw}}$ to a value $B>B_{\textrm{sw}}$, as in the left panel of Fig. 6. This change of magnetic field will be accompanied by a dramatic decrease in the visibility of the absorption line corresponding to the even transition at frequency $\omega_{\textrm{even}}$. Indeed, for $B>B_{\textrm{sw}}$ the ground state of the junction is the odd parity state $|1\rangle$, and at low temperatures $k_{B}T\ll E_{1}$ the occupation probability of the even parity state $|V\rangle$ is negligible. The low occupation probability of the state $|V\rangle$ and the smallness of the matrix element $j_{\textrm{2,1}}$ combine to yield a dramatic dimming of the absorption line taking place at $B=B_{\textrm{sw}}$ (see Fig. 12).

VI Conclusions

In this work, we have investigated several consequences of the competition between Zeeman and spin-orbit couplings on the Andreev bound states in semiconducting nanowire Josephson junctions. Overall, as one may have expected, spin-orbit coupling tends to reduce the effect of the Zeeman coupling on the Andreev bound states. We have seen several examples of this general trend. First, as discussed in Sec. III.4, spin-orbit coupling tends to suppress the g-factor $g_{A}$ of the Andreev bound states, potentially resulting in very small energy splittings of the Andreev doublet for small magnetic fields. The measurement of $g_{A}$, in tunneling or supercurrent spectroscopy experiments, may allow one to estimate the strength of the spin-orbit coupling. Second, spin-orbit coupling also suppresses the occurrence of fermion parity switches in the topologically trivial phase of the nanowire (see Sec. III.5 and Fig. 8). As discussed at the end of Sec. V, fermion parity switches should be easily detectable since they are accompanied by a drastic dimming of the absorption spectrum. The knowledge of the switching field $B_{\textrm{sw}}$ at which fermion parity switches take place can also be used to infer the strength of the spin-orbit coupling. Finally, spin-orbit coupling prevents the occurrence of level crossings between the Andreev bound states and the continuum. Combined with the suppression of the proximity-induced energy gap in a magnetic field, this leads to a non-monotonic dependence of the Andreev bound state energies on $B$, see the left panel of Fig. 6. The bending of the Andreev level $E_{2}$ due to the repulsion from the continuum causes a slow decrease of the even transition frequency in magnetic field, $\omega_{\textrm{even}}(0)-\omega_{\textrm{even}}(B)\propto B^{2}$ for small $B$.

Our theoretical results are in good agreement with several aspects of the existing experimental data which motivated the development of the work presented here van Woerkom et al. . In particular, we elucidated that the quadratic suppression of $\omega_{\textrm{even}}$ with increasing magnetic field can be understood in terms of the interplay of Zeeman and spin-orbit coupling. The occurrence of fermion parity switches is also compatible with the observation that the even transition visibility vanished at a field larger than $300$ mT. This threshold can be well understood within our theory assuming reasonable values of $g$ and $\alpha$ van Woerkom et al. . At the experimental level, it would be very valuable to study directly the single-particle energy spectrum via either tunneling or supercurrent spectroscopy. This would allow a measurement of the Andreev bound state g-factor $g_{A}$ as well as a precise determination of the switching field $B_{\textrm{sw}}$, both of which can be directly compared to our theory.

Our results are all based on the one-dimensional nanowire model of Eq. (1) treated within the Andreev approximation, i.e. by linearizing the normal state dispersion. This approximation amounts to neglecting the normal reflection amplitude in favor of the Andreev reflection amplitude when considering the two interfaces of the S-N-S junction. It requires that either the chemical potential $\mu$ or the spin-orbit energy $m\alpha^{2}$ are much larger than the induced superconducting gap $\Delta_{0}$. As an extension of this work, it may be valuable to relax the Andreev approximation. In particular, it would be interesting to study the Andreev spectrum of the model in the regime $\mu\ll\ m\alpha^{2},\Delta_{0}$ and $m\alpha^{2}\sim\Delta_{0}$, which may be relevant for the Majorana applications of the semiconducting nanowires. This may give more accurate predictions for the Andreev g-factor $g_{A}$ and the current matrix elements in the regime of low chemical potential.

It will also be important to extend this work beyond the model of Eq. (1), in order to capture more accurately the complexity of real devices. Nanowire junctions may naturally host more than one transport channel, and physical effects not included in this work, such as the orbital effect of the magnetic field, may have an important influence on the Andreev bound state properties. In particular, the orbital effect of the magnetic field provides an additional contribution to the reduction of $\omega_{\textrm{even}}$. Although this contribution could be heuristically ruled out to be the dominant one in the current nanowire experiments, it would be important to have quantitative theoretical estimates.

Finally, the magnetic field dependence of the absorption spectrum in the presence of multiple transport channels stands out as a particularly interesting avenue for future research, both theoretically and experimentally. In such a situation, a new type of low-frequency transitions may become visible, in which a Cooper pair is excited to a pair of Andreev levels belonging to different transport channels. In the topological phase, we expect that these inter-channel transitions can carry a signature of the Majorana bound states in the form of a kink in the phase dependence of the absorption spectrum, similar to the effect predicated in long Josephson junctions Väyrynen et al. (2015). Notably, this type of measurement is not limited by stringent requirements on fermion parity relaxation times, as opposed to other signatures of topological Josephson junctions.

Acknowledgements.

We acknowledge stimulating discussions with A. Geresdi, D.J. van Woerkom, H. Pothier, R. Lutchyn, T. Hyart and S. Park. BvH was supported by ONR Grant Q00704, JV and LG acknowledge the support by NSF DMR Grant No. 1603243.

Appendix A Bound state wave functions at zero field

Following the procedure outlined in the main text, we find the two following bound state wave functions $\Phi_{1}(x)$ and $\Phi_{2}(x)$ at $B=0$. They are a tensor product of a position-dependent part and a position-independent spinor in spin grading:

[TABLE]

Here $\chi_{\uparrow}=(1,0)^{T}$ and $\chi_{\downarrow}=(0,1)^{T}$ are the eigenspinors of $\sigma_{z}$, and $\Phi_{A}(x)$ is a space-dependent vector in Nambu and left/right gradings:

[TABLE]

with:

[TABLE]

The expression above is valid for $\mu\gg\Delta_{0},m\alpha^{2}$ and in the phase interval $\phi\,\in\,[0,2\pi]$. The wave functions for negative phase can be determined by applying the time-reversal symmetry operator $is_{x}\sigma_{y}\mathcal{K}$. In the opposite regime $m\alpha^{2}\gg\mu,\Delta_{0}$, the wave functions are identical except that the oscillating term $\operatorname{e}^{-i\alpha k_{F}/v_{F}}$ is replaced by $\operatorname{e}^{-i\mu x/\alpha s_{z}}$ and that $\alpha$ replaces $v_{F}$ in the expression for the coherence length $\xi_{A}$ of the bound state. The wave functions above are properly normalized to unity: to see this, it is convenient to use the relation $\Delta_{0}^{2}\,(1-\tau)=(E_{A}-\Delta_{0}\sqrt{\tau}\cos\phi/2)\,(E_{A}+\Delta_{0}\sqrt{\tau}\cos\phi/2)$.

Appendix B Derivation of the different contributions to the g-factor

In the main text, Secs. II.1 and II.2, we linearized the spectrum in two limits of either large $\mu$ or large $m\alpha^{2}$. In this Appendix, we show that both limits can be obtained from a single linearization which is valid on a strip of width $\mu+m\alpha^{2}$ around the Fermi level. This linearization is achieved by a projection

[TABLE]

where the fields $\psi_{L,R}$ are slowly varying. For example, the kinetic term in the Hamiltonian density becomes

[TABLE]

We used here $\mu=\frac{1}{2}mv_{F}^{2}$. When we project the Zeeman term $-\frac{1}{2}g\mu_{B}B\psi^{\dagger}\sigma_{x}\psi$ to low energies using Eq. (67), we obtain two terms,

[TABLE]

where

[TABLE]

The first term couples co-propagating states only and it is important when spin-orbit strength is not too large, $m\alpha^{2}\lesssim\Delta_{0}$. It leads to Eq. (46) of the main text, which can be derived by evaluating the matrix elements of $O_{\rightrightarrows}$ by using the wave functions from Eq. (65).

The second term, $O_{\rightleftarrows}$, mixes counter-propagating states and therefore it only contributes to the g-factor in the presence of scattering at the junction. Furthermore, it is important only for states near $k=0$ and when $\mu\ll m\alpha^{2}$, in which case the oscillating exponent vanishes. States belonging to the outer branches at finite momentum have opposite $s_{z}$ and $\sigma_{z}$ eigenvalues, and in this case the $O_{\rightleftarrows}$ term oscillates fast and is negligible. In the limit $\mu/m\alpha^{2}\to 0$, we thus obtain

[TABLE]

After calculating the matrix elements of $O_{\rightleftarrows}$ with the wave functions from Appendix A, we find Eq. (47) of the main text. Note that correction of order $\mu/m\alpha^{2}$ to the matrix elements of $O_{\rightleftarrows}$ cannot be reliably computed within the linearized Hamiltonian.