Fraunhofer pattern in the presence of Majorana zero modes

We propose a new platform to detect signatures of the presence of Majorana bound states (MBSs) in the Fraunhofer pattern of Josephson junctions featuring quantum spin Hall edge states on the normal part and Majorana bound states at the NS interfaces. We use a tight-binding model to demonstrate a drastic change in the periodicity of the Fraunhofer pattern when comparing trivial and non-trivial regimes. We explain these results in terms of the presence of additional parallel-spin electron-hole reflections, which due to the spin-momentum locking, occur as cross Andreev reflections, accumulating a different magnetic flux and yielding a change in the Fraunhofer periodicity. We show that this detection scheme exhibits some advantages compared to previous ones as it is robust against disorder, finite temperature and works in equilibrium. Furthermore, we introduce a scattering model that captures the main results of the microscopic calculations with MBSs and extend our discussion to the main differences found using accidental zero energy ABSs.


I. INTRODUCTION
6][7] Furthermore, Majorana zero modes exhibit a fractional nature, and therefore, they always appear in pairs as the result of delocalizing the information of a single fermion c = (γ 1 + iγ 2 )/2 onto the boundaries of the system.For example, in one-dimensional (1D) systems a pair of zero-dimensional bound states becomes localized at the boundaries of the topological superconductor 4,6,7 , whereas in two-dimensions, they emerge as an even number of chiral vortices. 8Apart from their intrinsic interest, there are practical applications due to their individual charge neutrality, which protects them from the local coupling to environmental charge fluctuations, and more importantly, due to the possibility of performing computational operations by the adiabatic exchange of Majorana bound states, also known as braiding, which leads to an adiabatic state change within a degenerate ground state manifold 3 .For further information and references there is a collection of reviews that cover different aspects of these exotic particles, see for example Refs.9-15  Signatures of MBSs are found in several transport experiments.Two of the most studied experiments, the zero bias conductance in a NS junction 8,16,17 and the fractional Josephson effect 4 in the Josephson junction, report deviations from the theoretical predictions.In the case of the zero bias conductance, the signature consists of a quantized value G = 2e 2 /h at zero temperature.How-ever, most of the experimental realizations have shown a substantially suppressed value [18][19][20][21] , and only in one of them the conductance is consistently close to 2e 2 /h, exhibiting deviations below and above 22 .There are different explanations that can justify such deviations, some of them are compatible with a topological ground state, like effects of a finite temperature, or a finite coupling to the opposite MBS, while other explanations are compatible with a trivial ground state, like the scattering with quasi-Majorana bound states [23][24][25] or coupling to trivial zero-energy Andreev bound states (ZEABSs) 26,27 .The situation is similar for the fractional Josephson effect, which can be probed by means of the Shapiro experiment, where odd Shapiro steps vanish 4,5,[28][29][30][31][32] or the Josephson radiation 33,34 .Experimentally however, in most cases only few odd steps are suppressed when the driving frequency is low enough [35][36][37][38] , and only one contribution reports the lack of the first four odd steps 39 .In this occasion, the signal can also be explained in terms of a topological state that coexists with trivial ones 31,40,41 , however, it is also possible that the behavior can be explained in terms of non-adiabatic transitions between Andreev bound states [30][31][32][42][43][44][45][46] . Therefre, the need for detection schemes that are more susceptible to the triplet superconducting nature of the MBSs, such as the measurement of triplet correlations by coupling the topological superconductor to a spin-dependent current in a 3-terminal setup 47,48 or the spin susceptibility of a Josephson junction are of utmost importance 49 .
In this work, we investigate signatures of the presence of MBSs at the NS interfaces that arise in the Fraunhofer pattern of a planar Josephson junction featuring quantum spin Hall edge states on the normal part, see Fig. 1.Here, the spin-momentum locking of the helical edge states forces local and crossed Andreev reflections (LAR and CAR) to take different spin-symmetries, that is, r he ss for the LAR and r he ss for the CAR.Note that in trivial junctions, the CAR is zero for a homogeneous or effectively linear in momentum spin-orbit coupling.In contrast, in the presence of a MBS, LAR and CAR are of the same order.Interestingly, the presence of a magnetic flux can unravel both contributions since electrons accumulate a different magnetic flux when performing a LAR or a CAR process 50,51 .In this way, we investigate the critical current of a Josephson junction as a function of an applied magnetic flux, also known as Fraunhofer pattern, expecting to observe an abrupt change on the profile when passing from the trivial to the topological regime.Since this detection scheme is performed in equilibrium, it removes the complications of parity conservation required in Fu and Kane's proposal 5 and the possibility of non-adiabatic transitions that can appear in the Shapiro experiment or the Josephson radiation.Furthermore, using a scattering model we demonstrate that the change in Fraunhofer pattern profile is absent in the presence of accidental zero energy Andreev bound states if the direction of the Zeeman field is parallel to the helicity operator of the quantum spin Hall states.
In the same spirit but with a different mechanism, other proposals suggest to use quantum Hall edge states to observe a change in the Fraunhofer pattern profile when MBSs are present 52 .Alternatively, one can study the change of periodicity in the conductance of a NSN junction interferometer fabricated on a 2D-topological insulator 53 .A recent experiment shows an abrupt π shift in the Fraunhofer pattern of a Josephson junction in a SQUID loop configuration, formed by two parallel Al/InAs Josephson junctions 54 .There, the phase shift is attributed to a topological phase transition produced by an external in-plane Zeeman field, but the mechanism for this transition is not further discussed.
The paper is organized as follows.In Sec.II we introduce the proximitized BHZ model and the topological phase studied in this work.Then, in Sec.III we explore the Fraunhofer pattern of the thin superconducting junction where to expect a change in the Fraunhofer pattern periodicity.Finally, in Sec.IV we introduce an effective scattering model that captures the essential ingredients of the microscopic model and compute the resulting Fraunhofer pattern in the presence of MBSs and accidental zero energy crossings.Further information is given in several appendices.

II. MICROSCOPIC MODEL AND TRANSPORT FORMALISM
In this section, we introduce the proximitized Bernevig-Hughes-Zhang (BHZ) model 55 and the transport formalism used to calculate the critical current of the Josephson junction.

A. Tight-binding Hamiltonian
The proximitized Bogoliubov de Gennes (BdG) BHZ Hamiltonian is given by H = (1/2) dxdyΨ † HΨ with with the Fermi energy E F and H h = T H e T −1 is the Hamiltonian for holes, which is obtained from the one for electrons (H e ) by performing a time-reversal transformation T = is y σ 0 C, with C being the complex conjugate operator.Here, ∆(y) = ∆ 0 (y) exp[iφ(y)]1 4×4 with the pairing potential ∆ 0 (y), which is finite and constant for |y| > L n /2 and 0 otherwise, with φ(y) = ϕ/2 [φ(y) = −ϕ/2] for y > L n /2 (y < −L n /2), see Fig. 1.This Hamiltonian is written in the electron-hole basis Ψ(x, y) = (ψ e , ψ h ) t , with where c a,σ is the annihilation (creation) operator of an electron with orbital a = E, H and spin ↑ / ↓ at position (x, y).The Pauli matrices σ i and s i , span the orbital and spin degrees of freedom.2|M | gives approximately the energy band gap of the system and its sign determines the topological character of the Hamiltonian: M > 0 (M < 0), sets the Hamiltonian in the trivial (topological) insulating regime.Its topological character expresses through the emergence of helical quantum spin Hall edge states.In the presence of timereversal symmetry, these quantum spin Hall edge states propagate without backscattering even in the presence of disorder 55 .H R is the Rashba spin-orbit coupling with α being the Rashba spin-orbit constant, which can be tuned by an external electric field 56 .H D accounts for the bulk inversion asymmetry (BIA) contribution with δ 0 , δ e , δ h , are bulk inversion asymmetry parameters, which are specific to the material under consideration.The effects of an external magnetic field H enter via the Peierls substitution, discussed below, and also via the Zeeman Hamiltonian H Z , with in-plane H ∥ and out of plane H ⊥ components.The corresponding Zeeman energies are given by with the parallel g ∥ ≈ 2 and perpendicular electron/hole band g-factor g e/h ⊥ ≈ 22.7(−1.21)and the Bohr magneton µ B .Further, we have introduced the in-plane spin-operator s θ = cos θs x + sin θs y with θ the in-plane polar angle.Guided by HgTe, we use different in-plane and out of plane g-factors 57 .While and thus, we assume a negligible g h ⊥ ≈ 0. Nevertheless, similar results can be found using g e ⊥ = g h ⊥ , as realized in InAs/GaSb wells 58 .Following standard finite difference methods, we discretize the Hamiltonian on a 2D lattice, turning the continuum momenta kx/y Ψ(x, y) x/y (Ψ i x/y +a x/y − 2Ψ i x/y + Ψ i x/y −a x/y ).In the normal part of the junction |y| < L n /2, we add a perpendicular magnetic field, whose orbital effects enter via the Peierls substitution kx,y → kx,y − (e/ℏ)A x,y , with the Landau gauge A = |H ⊥ |xe y .Thus, ψ † e,iy ψ e,iy±ay → exp[±i (πΦ p /Φ 0 )i x ]ψ † e,iy ψ e,iy±ay , with the flux quanta Φ 0 = h/2e and the magnetic flux in a plaquette Φ p = A|H ⊥ | for |y| < L n /2 and 0 otherwise, and with A = a x a y being the area of a single plaquet.Besides, we assume a negligible Zeeman contribution generated by the magnetic field responsible for the magnetic flux.

B. Topological superconductivity
There are different ways to induce Majorana bound states on the proximitized quantum spin Hall edges.The most direct one, a proximitized helical edge state 59 , is not relevant here because without conserving the parity, its corresponding Fraunhofer pattern gives no difference in the trivial and topological regimes.Here, we explore a different configuration closer to the Rashba quasi 1D well, where the proximitized BHZ model supports a chiral Majorana edge mode for a Zeeman field above the critical field B ⊥ > B c 60,61 in the NS junction.In addition, in the quasi 1D limit a pair of proximitized quantum spin Hall edge channels can exhibit a superconducting topological phase when the width of the sample is smaller than the localization length of the quantum spin Hall edge states ξ qsh , that is, W ≲ ξ qsh 48 .To understand this scenario, we first represent the non-proximitized spectrum of a thin sample in Fig. 2(a), which develops a gap at the Dirac point.In the vicinity of the gap, the resulting spectrum exhibits the same structure as the Rashba wire, that is, parabolic bands shifted in k by the spin-orbit coupling ∆k = ±k so , see magnification in Fig. 2(a).Hence, it is not surprising to expect an inversion of the proximitized superconducting band produced by a Zeeman field applied perpendicularly to the effective spin-orbit field direction, see bottom panels in Fig. 2(c).As a result, two Majorana bound states appear at the extremes of the superconductor, which becomes visible imposing hard-wall boundary conditions.See LDOS at the NS interface in Fig. 2(b) 48 .It is important to remark that even though we need a finite overlap between opposite quantum spin Hall edge states on the superconducting part, the quantum spin Hall edge states placed on the normal part propagate without scattering to the opposite edge.This is possible by setting C in the normal region, i.e.C N , in such a way that the Dirac point is away from the Fermi energy.Alternatively, one can use the same chemical potential as in the superconducting region, i.e.C N = C S , with a different width W N > W S on the normal part and superconducting parts.Although we use the former configuration (C N ̸ = C S ), the latter (W N > W S ) might be easier to accomplish experimentally since it does not involve a spatially dependent gate voltage control.
Experimentally, there are already two scenarios where non-proximitized quantum spin Hall modes of opposite edges overlap.Firstly, the addition of an electric field can push the edge modes towards the bulk 62 .Secondly, materials like bismuthene on SiC develop topological line defects through which the quantum spin Hall edges propagating on opposite sides can overlap 63 .

Set of parameters
Unless it is specified otherwise, we use the following set of parameters: A = 373 meV nm, B = −857 meV nm 2 , D = −682 meV nm 2 , M = −10.0meV.The parameters associated to spin-orbit coupling are α = 0.0, δ 0 = 1.6 meV, δ e = −12.8meV nm, δ h = 21.1 meV nm, T = 0.1 K and ∆ 0 = 0.15 meV.In the tight-binding calculations, we use a discretization constant a ≡ a x = a y = 2.4 nm.We set the Fermi energy E F = 0 and the chemical potential on the normal and superconducting parts is tuned by C N and C S , respectively.
In this work, we use the magnetic field for two different purposes: a Zeeman field to invert the superconducting gap and the magnetic flux to study the Fraunhofer pattern.Even though both fields are applied in the zdirection, it is possible to study the critical current as a function of both contributions independently because of the different order of magnitude of each contribution.Here, the Zeeman field is of the order of B ⊥ ∼ 1T, while the orbital field is a few mT.

III. FRAUNHOFER PATTERN
The Fraunhofer pattern, i.e. the critical current as a function of the magnetic flux, contains important information about the spatial distribution of the supercurrent in a Josephson junction 64 .For example, when the normal part of the junction is dominated by bulk modes, the critical current oscillates and decays for an increasing magnetic flux.Whereas, when the supercurrent is carried along the edges (quantum spin Hall or quantum Hall states), the critical current oscillates with a period (Φ/Φ 0 or 2Φ/Φ 0 ) without decaying.Using standard equilibrium Green's function methods introduced in App.A, we calculate the critical current of a Josephson junction with quantum spin Hall edges on the normal part in the absence and presence of Majorana bound states at the NS interfaces.
In the 2D limit, where chiral Majorana modes propagate at the boundary of the superconducting slabs, the Fraunhofer pattern does not exhibit a qualitative change when comparing the trivial and topological regimes, see Fig. 11 in App.B 1.This is caused by the dominant coupling between normal and proximitized quantum spin Hall edge states, which gives rise to a supercurrent contribution based on LAR, yielding negligible traces of the presence of the chiral Majorana contribution.Therefore, in order to observe a specific signature of the MBS in the Fraunhofer pattern profile, we need to reduce the supercurrent carried by local ABS so that we can enhance the coupling between the quantum spin Hall edge states and the central part of the superconducting slab, where the MBSs are located.To this aim, we induce MBSs directly on the proximitized quantum spin Hall edge states by considering a thin superconductor junction.Alternatively, there is another way to increase the coupling between the helical modes and the MBSs, that involves a positive mass M > 0 in the superconducting leads, having no helical edge states in the normal state.Since the results are analogous to the thin superconducting junc- tion we place this setup and its results to App.B 2. In addition, we show there that it is also possible to observe the transition with an in-plane magnetic field B ∥ .
In the trivial regime (B ⊥ < B c ), the Fraunhofer pattern exhibits the usual SQUID pattern with period p ∼ Φ 0 , see Figs. 3 (a) and (b).In contrast to the wide superconductor junction, see App.B 1, we observe a reduction of one lobe every Φ ≈ 2Φ 0 .This beating pattern, also known as even-odd flux quanta effect, occurs because the quantum spin Hall edges placed on opposite sides become partially coupled at the superconductor interface 65 , resulting into normal electron-electron and hole-hole reflections involving both edges.In this situation, particles can encircle the normal part several times before an electron-hole scattering event takes place, accumulating extra magnetic phase factors, and modifying the Fraunhofer pattern.Increasing the Zeeman field above the critical field (B c ∼ 1.5 meV), we observe the suppression of the originally "odd" lobes placed around Φ ≈ (2n − 1)Φ 0 (with n being an integer), yielding a Fraunhofer pattern with a periodicity of p ≈ 2Φ 0 , see linecuts in Fig. 3(b).Numerical calculations show a similar change in periodicity for different set of parameters, such as, C S , ∆ 0 as long as B ⊥ > B c , and thus, becoming a robust signature of the topological transition.Furthermore, we can observe a constant shift of the Fraunhofer pattern as a function of the Zeeman energy B ⊥ , which can be gauged into the Peierls phase since the Zeeman field is parallel to the helicity operator of the quantum spin-Hall edges, see more details in App.B 1.
The abrupt change of the Fraunhofer periodicity can be understood from the new scattering channels introduced by the presence of MBSs, which enhance triplet electron-hole reflections with parallel spin, i.e. r he ss .Due to spin-momentum locking, the ABSs involving r he ss take place on opposite edges, while r he ss on the same edge.Due to this geometrical difference, the magnetic flux accumulated in each process is different.Crossed ABSs do not accumulate a magnetic flux, while local ABS accumulate the flux Φ.In the tunneling limit 51 , the critical current is given to lowest order in the particle tunneling by with I CAR/LAR the respective crossed and local critical currents.Using this simple equation, we can extract that the presence of MBSs at the NS boundaries generates a local and crossed supercurrents with approximately the same strength, which is in accordance with the electronhole reflection coefficients of MBS 47 and the suppression of CAR in the trivial regime 44,66,67 .
Intra and inter-edge GFs-To support this interpretation, we compute the intra and inter-edge GFs defined respectively as, as a function of the frequency ω, see Fig. 4. Here, the subscripts T, B refer to top and bottom edges, respectively.Note that these GFs are the result of tracing the orbital degrees of freedom and they are defined at the center of the junction y = 0. Besides, we have used a superconducting phase difference ϕ = 0 and the magnetic flux Φ that gives rise to the maximal critical current for the given Zeeman field.We can observe that in the trivial regime (B ⊥ < B c ), intra and inter-edge GFs contain only the singlet component  Robustness against disorder-Due to the presence of the magnetic field, time reversal symmetry is broken, and thus, the helical edge states are no longer protected against backscattering.Therefore, we check the robustness of the signature in the presence of disorder.To this aim, we set two stripes of length L trench and width W at both NS interfaces, see Fig. 10.This is consistent with an actual experiment, where the superconducting leads can introduce some disorder to the semiconductor system.In this way, we consider an on-site disorder potential that is randomly assigned in the range of ±σ.We first note that for values of disorder strength σ = 1 meV, the Fraunhofer pattern exhibits the same change of periodicity as the clean system, even when L trench extends over the whole setup (∼ 0.74 µm), not shown.For larger disorder strenghts, σ = 5 meV, we observe no change in the Fraunhofer pattern relative to the clean system for L trench ≲ 100 nm (not shown).However, increasing further L trench to 150 nm the height of the odd lobes decays in the trivial regime, see two linecuts in Fig. 5.Although in this case the periodicity is strictly not doubled as a function of the Zeeman field, it is still possible to appreciate the impact of the Zeeman field on the Fraunhofer profile by performing the Fourier transform on the critical current, that is, In Fig. 6, we represent the Fourier transform I c (f ) as a function of the frequency f and the Zeeman energy B ⊥ .
We can observe that in the trivial regime (B ⊥ < B c ), I c (f ) exhibits three maxima around f ≈ 0, 0.5 and 1.0, and smaller and narrower contributions for 0 < f < 1.
Here, the contribution around f ≈ 0 is related to the constant background observed in Fig. 3, and the contributions with frequencies at f < 0.5, correspond to the small oscillations observed on top of the odd lobes.Moreover, the contribution close to f ≈ 0.5 is related to the evenodd flux quanta effect mentioned above, whereas f ≈ 1 is the characteristic frequency of a wide quantum spin-Hall Josephson junction. 68For B ⊥ > B c , we can observe an abrupt change in I c (f ), with a dominant component at f ≈ 0.5 resulting from the reduction of the odd lobes, which serves as an indication of the presence of MBSs as a function of the Zeeman field.

IV. SCATTERING MODEL
In order to have a deeper insight into the periodicity changes introduced in the Fraunhofer pattern by the presence of MBSs, we construct a phenomenological scattering model including the main ingredients participating in the setup, that is, helical modes propagating around the normal part in the presence of a magnetic flux, and two superconducting slabs accounting for scattering events with a bulk superconductor and MBS, see Fig. 7. Furthermore, this model also allows us to consider accidental zero energy modes instead of MBSs and to study the resulting Fraunhofer pattern.
The supercurrent at temperature T can be written in terms of the S matrices S A and S N , namely, 65,69,70 where the sum runs over the Matsubara frequencies ω n = πkT (2n + 1).Here, S N is the S matrix that accounts for the propagation along the normal part, while S A accounts for the scattering events taking place at the left and right NS interfaces.
To model the processes occurring at each NS interface, we combine three S matrices represented by crossed boxes in Fig. 7, see more details in App. C. Two of them (green crossed boxes) capture the physics of the scattering events between the quantum spin Hall edges and the bulk superconductor, where we also account for finite electron-electron deflection probability (Γ d ) that connects quantum spin Hall edge states placed on opposite sides.Additionally, the scattering events produced by MBSs are modelled by the addition of a third scatterer (red crossed boxes) centered in between the bulk superconducting scatterers.Similarly as in the superconducting bulk scatterers, the intensity of the scattering processes with the MBSs scatterers is controlled by the parameter Γ M : for 1 − Γ M = 0 particles propagate along the NS interface without scattering with the MBS, while for 1 − Γ M = 1 any particle propagating along the NS junction scatters with the MBS.
As a first step, we recover a similar Fraunhofer pattern as the ones observed using the microscopic model for the trivial regime i.e.B ⊥ = 0, see blue linecuts in Figs. 3  and Fig. 8.To do so, we set Γ M = 1.0, which sets a zero coupling to the MBSs and a finite Γ d > 0, which introduces finite deflection probability, connecting opposite edges with Γ d = 0.2.Here, the resulting Fraunhofer pattern exhibits the previously observed even-odd flux quanta effect by the reduction of one of the lobe maxima every 2Φ 0 , see blue curves in Fig. 8.Note that in contrast to the microscopic model, here the reduced lobes are always placed at odd multiples of Φ 0 .In general, the parity of the reduced lobes depends on the phase factors picked up during the deflection scattering processes 65 .For simplicity, we have set all the phases picked up during the deflection to zero, and therefore, we cannot change the parity of the reduced lobes.

A. Majorana bound states
Setting a finite scattering probability with the MBSs (1 − Γ M > 0), introduces a finite electron-hole process with parallel spin 71 .Using the Mahaux-Weidenmüller formula we can obtain the reflection amplitudes of an isolated MBS, namely 47,72,73 where ) and t σ is the effective coupling amplitude between the helical edge states and the Majorana bound states, see more details in App. C. As we have discussed above, these scattering processes have a strong impact on the Fraunhofer pattern because the resulting crossed ABS do not accumulate a net magnetic flux.
Guided by the microscopic model, we assume t σ to be real 74,75 .Concretely, we analyze the tunnel amplitude configuration: t ↑,L = t ↓,L , t ↑,R = −t ↓,R , which reproduces approximately the Fraunhofer patterns observed using the microscopic model, see red dashed curves in Fig. 8.Moreover, the form of the Fraunhofer pattern can only be recovered taking into account that the MBSs appear on both sides, thus, the change in the Fraunhofer profile also provides a non-local probe of a pair of separated MBSs appearing at the left and right interfaces.

B. Zero energy Andreev bound states
Accidental zero energy Andreev bound states (ZEABS) can arise in the presence of a Zeeman field when for example the NS interfaces of the proximitized region are less doped than the interior regions 23 .These states can give rise to similar signatures in the conductance as the MBSs, like a zero bias conductance peak.Thus, to test the robustness of the signatures observed in the Fraunhofer pattern, we replace the MBS S matrix, by one resulting from the scattering with an accidental zero energy ABS (ZEABS) and study under which conditions we observe similar changes in the Fraunhofer pattern as the ones observed with the MBSs.
We model the ZEABS by means of the Hamiltonian 47 where ϵ 0 sets the doping of these regions and it is assumed to be smaller than the chemical potential of the bulk proximitized region.∆ 0 is the pairing amplitude of this region and in general it can be different from the bulk ∆ 0 .Furthermore, B ξ is the Zeeman energy of a field applied in the ξ direction.This Hamiltonian exhibits zero energy states for B ξ = ∆ 2 0 + ϵ 2 0 independently of the Zeeman field direction.This region is coupled to the quantum spin Hall edge states by means of the tunnel Hamiltonian H T = s ω s ψ † s (x 0 )d s + h.c., with ψ(x) being the electronic field operator of the quantum spin Hall edges at position x = x 0 , where the ZEABS is placed.Moreover, the ratio ϵ 0 /∆ 0 modifies effectively the tunnel couplings ω ↑/↓ = 3ℏv F / √ L n between the quantum spin Hall edges and the ZEABS.The reason for this reduction is that the larger ϵ 0 , the higher B α we need to apply and the more spin-polarized the ZEABS become.See further details of the model in App. C.
As we have seen above, the presence of the accidental ZEABS depends on the specific conditions of the NS interfaces, like a spatial dependence of the chemical potential or pairing amplitude.Thus, it is possible that only one of the NS interfaces may develop a ZEABS.
For this reason, we analyze two different scenarios depending on the number of ZEABS present in the setup, that is, whether only one or two NS interfaces develop a ZEABS 7677,78 .Besides, we consider two different orientations of the Zeeman field: in-plane and an out of plane.In the latter case, we expect to observe no change in the Fraunhofer pattern because electron-hole reflections with parallel spin are zero 47 .In contrast, the ZEABS generated by an in-plane Zeeman field (independently of the direction) exhibit a finite electron-hole reflection coefficient with parallel spin and, as we will see, it can lead to similar Fraunhofer patterns as those shown by the MBSs.
In Fig. 9, we show a representative example of the Fraunhofer pattern generated with the scattering model using ZEABS.The main idea consists on checking whether it is possible to recover the Fraunhofer pattern observed in Fig. 8, that is, the suppression of the half of the lobes by switching on the coupling to the MBSs, that is, for Γ M = 1 → Γ M = 0. To this aim, we use Γ d = 0.2, and analyze all four configurations of Zeeman fields and number of ZEABS (Γ M = 0).In panels (a) and (b) we show the results for the out of plane Zeeman field with one and two ZEABS, respectively, while in panels (c) and (d) we use an in-plane Zeeman field and again one and two ZEABS.In all panels, we represent the Fraunhofer pattern for three values of ϵ 0 /∆ 0 = 0.5, 1, 2, which effectively change the coupling between the quantum spin Hall edges and the ZEABS.We observe that when only one ZEABS is present and independently of the Zeeman orientation, the relative height of the Fraunhofer pattern lobes does not change significantly with respect to the uncoupled case (p ≈ Φ 0 ).In addition, when ZEABS are present on both sides, the relative height of the lobes remains unchanged relative to the trivial case p ≈ Φ 0 only for an out of plane Zeeman field B ⊥ , see panel (b).However, for an in-plane field [panel (d)] the Fraunhofer exhibits a periodicity of p = 2Φ 0 .Note that this in stark contrast to the case of MBSs, where the scattering matrix constructed from Eq. ( 13) is independent on the direction of the Zeeman field (as long as it is perpendicular to the spin-orbit field).

V. CONCLUSIONS
In conclusion, using tight-binding calculations based on the BHZ model, we predict that the presence of MBSs on both NS interfaces leads to a drastic change in the Fraunhofer pattern periodicity with respect to the trivial regime, where the half of the lobes become suppressed.We explain the change in periodicity in terms of the introduction of triplet electron-hole reflection amplitudes r he ss enforced by the exotic Majorana bound states, which, when coupled to the spin-momentum locked quantum spin Hall edge states, leads to the coexistence of two different kinds of Andreev bound states: local and crossed, which are mediated by antiparallel and parallel spin electron-hole reflection amplitudes, respectively.Furthermore, this signature remains visible even in the presence of disorder strengths of 5 meV as long as the disorder trench is shorter than L trench ≲ 100 nm.For larger disorder trenches L trench ≥ 150 nm, the Fraunhofer pattern exhibits a larger even-odd effect in the trivial regime, making it more difficult to observe a change in periodicity as a function of the Zeeman energy.Nevertheless, making use of the Fourier transform, we observe an abrupt increment of the fractional frequencies for B ⊥ > B c .
Finally, we have developed a phenomenological scattering model with which we recover similar results as those observed in the microscopic model with MBS, supporting our physical picture of local and crossed Andreev reflections.Furthermore, we have used this model to explore the robustness of the detection scheme when accidental zero energy Andreev bound states (ZEABS) are present.We have found that when only one NS interface holds a ZEABS, the Fraunhofer pattern does not change its profile independently of the Zeeman orientation.Also when both NS interfaces hold a ZEABS the relative height of the Fraunhofer lobes remains unchanged for a perpendicular Zeeman field B ⊥ .However, when both NS interfaces exhibit ZEABS generated by an in-plane Zeeman field B ∥ , the resulting Fraunhofer pattern exhibits similarities to the topologically non-trivial case.edge states on opposite sides are decoupled.In Fig. 11(a) we can observe that the critical current shows the usual oscillatory pattern exhibited by quantum spin Hall junctions, with periodicity p = Φ 0 39,85,86 .In contrast to the case analyzed in the main text, the Fraunhofer pattern profile remains unchanged as a function of the Zeeman field B ⊥ even though the Zeeman field exceeds largely the critical field, i.e.B ⊥ > B c = 0.3 meV, see also linecuts in panel (b).Due to the large coupling between the normal and proximitized quantum spin Hall edge states, the critical current is mainly given by LAR.Thus, the presence of MBSs located in the central regions of the NS junctions is not noticed by the edge states.
The shift observed in the Fraunhofer pattern as a function of the Zeeman field in Fig. 11 comes from the fact that the Zeeman field is parallel to the helicity of the QSH edges and thus, it shifts the momentum in the same manner as the magnetic flux introduced by the Peierls substitution.To see this, we write the effective Hamiltonian for the QSH edge states with a Zeeman field, being s z and τ z the Pauli matrices describing the spin and the top/bottom edges.The addition of a magnetic flux enters via the Peierls substitution, shifting the momentum as ky → ky − e/ℏA y , where A y changes the sign from the top to the bottom edge, and thus, A y ∝ τ z , which cancel the τ z appearing in Eq. (B1) and entering similarly as the Zeeman term.
2. Mass domain junction MN < 0, MS > 0 with an in-plane Zeeman field B ∥ An alternative way to enhance the coupling between the quantum spin Hall edge states propagating on the normal part and the MBSs placed at the NS interfaces is to consider M S > 0 on the superconducting part.In this configuration, the quantum spin Hall edge states propagate around the normal part and can scatter with MBSs when present.Experimentally, this setup requires a spatial control of the mass sign, which can be realized, for example, in InAs/GaSb wells 58 .Alternatively, one could form a Josephson junction by coupling laterally the superconductor 87 to a quantum spin Hall system.In addition, we show that the topological transition is also visible using an in-plane Zeeman field.Note that since the spin-orbit coupling depends on both momenta k x and k y , it is not possible to invert the superconducting gap 88 .However, this restriction can be overcome since the width of the junction W is small, and thus, it imposes a confinement energy given by ∼ ℏ 2 /2m * W 2 , with m * being the effective mass of the band, is much larger than the spin-orbit energy ∼ δ e,h /W, that couples different subbands.
In the trivial regime, the Fraunhofer pattern exhibits a doubled period with respect to the one obtained in the proximitized junction, i.e. p ≈ 2Φ 0 , compare the blue linecuts in Figs.12(b) with the SQUID pattern resulting from a junction with M N,S < 0 with only LAR processes in Fig. 11(b).This scenario represents an extreme version of the even-odd flux quanta effect observed in the thin superconducting junction.Here, particles can encircle the normal part due to a finite Fermi velocity mismatch between the quantum spin Hall edge states and the superconducting part, which enhances the electronelectron reflection amplitude 89 .In the topological regime (B ∥ > B c ), we again observe an abrupt change in the Fraunhofer pattern periodicity.In the particular case shown in Fig. 12, the periodicity is halfed to p ≈ Φ 0 .However, the Fraunhofer pattern in the mass domain and the corresponding changes caused by the presence of a MBS do not exhibit a universal behavior.In contrast to the thin superconducting junction, here the coupling to the bulk s-wave superconductor is more visible yielding an extra source of LAR with respect to the electronhole reflections provided by the MBS.The ratio between the two contributions is effectively controlled by C S : the larger k F , the more favored the LAR processes become.In Fig. 13, we show the Fraunhofer pattern for two values of the Zeeman field: B ∥ = 0 (blue) and B ∥ = 1.2 meV (red) and nine values of the chemical potential on the superconducting part (C S ).Panels from (a) to (g) show a blue and a red curve corresponding to the trivial and topological regimes, respectively.The shaded panels (h) and (i) correspond to the trivial regime for both curves.In these two latter panels, an even number of superconducting bands are occupied and consequently above the critical Zeeman field, an even number of Majorana bound states form on each side, which overlap, gap out and yield a trivial superconductor 61,90 .This is known as the even/odd effect in the Majorana wells.Indeed, we can observe that the Fraunhofer pattern changes its profile for the six first panels (a)-(f), whereas close to the odd-even transition and within the even sectors, the Fraunhofer pattern does not exhibit a qualitative change in periodicity.In contrast to the thin junction regime where the superconducting electrodes are composed by a single proximitized mode, here, the presence of extra modes can make the distinction between the trivial and topological regimes more difficult, see for example panels c and h, where blue and red curves exhibit a similar periodicity.To distinguish between both cases we study the critical current for Φ = 0 as a function of the Zeeman field B ∥ , see Fig. 14.We can observe that in the case of a topological transition (blue curve), the critical current exhibits a sharp drop around B ∥ ∼ B c , whereas the curve corresponding to the "even" case decreases continuously.This has been proposed as a detection scheme in Ref. 26.

Appendix C: Scattering model
In this section, we provide the explicit form of S N and S A used in Eq. ( 12).
This model is adapted from the one presented in Ref. 65, where we use helical propagation along the NS interfaces and add the possibility of scattering with MBSs or trivial ZEABSs.
We start with S N , which gives the relation between R , here T /B stand for top and bottom edges.Note that the lack of some of the incoming and outgoing states comes from the spin-momentum locking of the edge states, e.g.there is no left-top incoming down-spin state a T e/h,↓ .The form of S N is given by which is an off-diagonal matrix since S N connects the incoming of the L/R side with the outgoing scattering states of the R/L side.Here, In the following we consider the short junction limit, i.e. ∆ 0 ≪ E T,t,b = ℏv F /L n .However, the periodicity of the Fraunhofer pattern is similar in intermediate or the long junction limit.
The S matrix describing the scattering processes at the NS junctions is given by S A,L/R describes the scattering events taking place at the left/right interfaces and fulfills Guided by the microscopic model, we include three different scattering processes: • Scattering with the bulk superconductor.
• Normal deflection between the top and bottom edges.
To do so, we model each NS interface by the combination (×) of three S matrices, that is, S A,i = S T green,i × S red,i × S B green,i , with i = L/R.Here, S T /B green,i describes the scattering between the quantum spin Hall edges with a bulk superconductor that includes a deflection probability towards the opposite edge.Besides, S red,i is the MBS scattering matrix on the i-side.
We start introducing S T /B green,i , which are represented by green crossed boxes in Fig. 7 and are obtained by the combination of two S matrices, i.e. S T /B green,i = S T /B pot,i × S Bulk,i , where S T /B pot,i describes the scattering processes between the quantum spin Hall edges and a (ficticious) normal potential that introduces a deflection probability (Γ d ) towards the opposite edge.For simplicity, we consider the same deflection probability on the top, bottom, left and right, and thus, the S matrix becomes independent of the T /B and L/R labels, that is, S T /B pot,L/R ≡ S pot and consequently S T /B green,L/R ≡ S green,L/R .For the sake of simplicity, we model the potential scatterer by means of the simplest S matrix that respects time-reversal symmetry and does not mix the spin degree of freedom, namely where ψ 1 and ψ 2 are phases acquired after crossing the potential barriers.Furthermore, we set both phases to zero and also neglect the possibility of a spin-rotation due to the presence of spin-orbit coupling.The scattering matrix connects the states where "aux" refers to the auxiliary states placed between the potential and the superconductor and "edge" refers to the edge states propagating along the normal part, see Fig. 7. Their specific components depend on the considered scattering center, for example, in Fig. 7 we illustrate the modes scattering on the top-right scattering center.Here, electron and hole components of the incom-ing top-right "T" (outgoing "inter,T") states with spin up fulfill (a inter,T e/h,↑ , a aux e/h,↑ ) t = S pot (b T e/h,↑ , b aux e/h,↑ ) t .A similar expression can be written for the spin-down components.Due to their helical character, the incoming and outgoing components with spin down are exchanged, that is, (a T e,↓ , a aux e,↓ ) t = S pot (b inter,T e,↓ , b aux e,↓ ) t .
So much for S pot , we now introduce S Bulk,L/R , which describes the electron-hole reflection processes mediated by a trivial bulk superconductor and fulfills which relates the incoming a aux e/h,↑↓ and the outgoing b aux e/h,↑↓ scattering amplitudes in the "aux" region.Using the Andreev approximation, the right reflection coefficients are given by r he , with χ = −arccos(E/∆ 0 ), while the left-reflection coefficients are obtained from the right coefficients by substituting ϕ → −ϕ.An alternative way to model S green,i is to relax the conditions that lead to the Andreev approximation.Away from this limit one can naturally obtain electron-electron reflection amplitudes which lead to a finite deflection probability 67,[91][92][93] .
To proceed further, we combine S pot in Eq. (C3) with the corresponding superconducting S matrix Eq. (C4) eliminating the "aux" modes.In the trivial case, the resulting matrix is given by The S matrices S red,i , represented by red crossed boxes in Fig. 7, are obtained in a similar fashion as S green,i , that is, S red,i = S pot,M × S MBS,i .Here, we combine two S matrices, one accounting for the probability of particles propagating along the NS junction without scattering with the MBS, Γ M = 1, and the MBS S matrix.While S pot,M is obtained from Eq. C3 using the replacement Γ d → Γ M .
We now introduce the matrix S MBS,i , which accounts for the scattering with MBSs or ZEABSs.In both cases, we use the Weidenmüller formula, 72,7347 ) with the density of states of the quantum spin Hall edges ν 0 = 1/πℏv F and W accounting for the tunnel amplitudes between the MBSs/ZEABSs and the quantum spin Hall edge states.Zero energy bound states-In the same way, we can obtain the S matrix describing the scattering processes between the quantum spin Hall edges and an accidental ZEABS.To do so, we use an effective model describing a proximitized region (∆ 0 ) that can develop a zero energy state under a Zeeman term (B ∥/⊥ ), namely 47 here d ( †) are fermionic destruction (creation) operators and ϵ 0 is the energy of the particle.This region is coupled to the quantum spin Hall edge states by means of the tunnel Hamiltonian H T = s ω s ψ † s (x 0 )d s + h.c., with ψ(x) being the electronic field operator of the quantum spin Hall edges at position x = x 0 , where the ZEABS is placed.
In-plane magnetic field-For a finite in-plane field Thus, in the low energy sector, the tunnel Hamiltonian reads, and the tunnel amplitudes t 1 [2],↑ = ω * ↑ sin(α)[cos(α)] and t 1 [2],↓ = ω * ↓ cos(α)[sin(α)].Now, W ∥ enters into the Wei-denmüller formula (C6), so we obtain the reflection coefficients.In this case, the analytical expressions for r ee ss ′ and r he ss ′ are cumbersome, thus, we restrict ourselves to evaluate it numerically.In general, we observe a finite electron-hole reflection amplitude with parallel spin.Thus, we expect to observe a similar (but not equal) behavior in the Fraunhofer pattern as with MBSs.
Perpendicular magnetic field-When the Zeeman field is applied in z-direction, that is, B ∥ = 0 and B ⊥ ̸ = 0 the gap also closes at B ⊥ = ∆ 2 0 + ϵ 2 0 and the low energy sector allows us to write where in contrast to the in-plane ZEABS, here r he ↑↑ = r he ↓↓ = 0. Thus, we expect to observe a completely different behavior as compared to the in-plane fields.

Interedge transmission probability and total conductance
To have more insight into the reflection properties of the combined S matrix S A,i , we calculate numerically the inter-edge electron-hole reflection probability T ↑↑ he and the zero bias conductance at zero temperature, given by with Γ M = 0 and E = 0.
Here, we consider two different scenarios, the first one with a MBS and the second one with an accidental zero energy ABSs with an in-plane magnetic field.In both cases we restrict to the case where all tunnel amplitudes are real t σ ∈ R.
We start with the MBS case, whose results are shown in Fig. 15(a,c).For zero deflection probability Γ d = 0,

FIG. 1 .
FIG. 1. [Panel (a)] Schematic representation of the planar quantum spin Hall Josephson junction.Here, pink and white areas represent superconducting and normal parts subjected to a magnetic flux Φ. [Panel (b)] Possible Andreev bound states present in the setup.Due to spin-momentum locking, the crossed (local) Andreev bound states hold parallel (antiparallel) spin.Solid and dashed arrows represent electronand hole-like helical edge states, respectively.
FIG. 2. Panel (a), energy dispersion of the (normal) overlapped quantum spin Hall edge states in the thin limit, with (CN = 0).Inset: zoom within the gap opened by the overlap between the quantum spin Hall edge states.Panel (b), local density of states as a function of the energy ω and Zeeman energy B ⊥ .Panel (c), BdG band spectrum inversion as a function of the Zeeman energy B ⊥ , similar to the energy dispersion of the Rashba wire in Refs.6 and 7.The width of the sample is W= 144 nm, the critical field Bc ≈ 1.3 meV and CS = −9.2meV.

FIG. 4 .
FIG. 4. Intra and inter-edge GFs as a function of the frequency ω/∆ for the trivial (a, b) and topological (c, d) regimes.The set of parameters is the same as in Fig. 3, with ϕ = 0 and (a) B ⊥ = 0, Φ/Φ0 = 0 and (b) B ⊥ = 3.5 meV, Φ/Φ0 = −0.85 which lead to the maximum critical current for each Zeeman fields.

FIG. 5 .
FIG. 5. Fraunhofer pattern for the thin superconducting junction with a disorder trench of L trench = 150 nm and a disorder strength σ = ±5 meV.The rest of the parameters are the same as in Fig. 3.
, see panels a and c.In contrast, in the topological regime (B ⊥ > B c ), we can observe that the triplet components are of the same order as the singlet one, see panels b and d.

FIG. 6 .
FIG. 6. Absolute value of the Fourier transform of the critical current in the presence of disorder as shown in Fig. 3 as a function of the characteristic frequency f and the Zeeman field B ⊥ .

FIG. 7 .
FIG. 7. Panel (a): Josephson junction setup with helical modes confined in the normal region threaded by a magnetic flux Φ.We highlight the normal and superconducting regions, where the S matrices SN and S A,L/R act by light blue and red, respectively.On the superconducting sides, the crossed boxes represent trivial (green) and topological (red) superconducting scattering centers, where the helical edge states perform electron-hole reflections.In panel (b), we show the structure of SN and SA, which connect the incoming and outgoing (a/b) scattering states on the left and right sides (L/R), by means of SA(bL, bR) t = (aL, aR) t and (bL, bR) t = SN (aL, aR) t .In panels (c) and (d), we detail the top-right (trivial) and center-right (topological) scattering processes.

FIG. 9 .
FIG. 9. Fraunhofer pattern obtained using the scattering model with ZEABS scatterers for B ∥/⊥ = ∆ 2 0 + ϵ 2 0 and the same deflection probabilities as in Fig. 8, i.e.Γ d = 0.2, ΓM = 0.0.We set two different magnetic field directions B ⊥ (a and b) and B ∥ (c and d) and use either one (a and c) or two (b and d) ZEABS.We use in all panels three different values of the ratio ϵ0/∆0, which controls indirectly the coupling between the quantum spin Hall edge states and the ZEABS, see more details in App. C.

FIG. 10 .
FIG.10.Sketch of the discretized Josephson junction, with the mesh ax = ay.Here, we split the junction into two slabs to visualize the recursive inversion method, where the surface Green's functions gLL and gRR (blue areas) are coupled by V LR/RL and contain the information of the semi-infinite normal-superconducting leads, see further details in Sec. A. Besides, we introduce a magnetic flux Φp that threads each plaquette on the normal part of the junction.

FIG. 11 .
FIG. 11.Panel (a), critical current as a function of the magnetic flux and Zeeman field B ⊥ for a proximitized Josephson junction, with M = −8.0meV for the whole junction.Panel (b), critical current for two different Zeeman fields corresponding to trivial and topological superconductors.The parameters used are ∆0 = 4 meV, W = 1 µm, Ln = 0.35 µm, CN = −3 meV, CS = 10 meV and lattice constant a = 5 nm.

FIG. 12 .
FIG. 12. Fraunhofer pattern for the mass domain configuration: Panel (a) critical current as a function of the flux and parallel Zeeman field for the mass domain configuration.The parameters of the model are M S/N = ±10 meV, CS = 12.0 meV, CN = 2.0 meV, W= 0.24 µm and Ln = 0.4 µm.Panel (b): Fraunhofer pattern for two different linecuts B ∥ = 0.0 (trivial regime) and B ∥ = 1.3 meV (topological regime).

FIG. 13 .
FIG.13.Critical current for the mass domain configuration as a function of the magnetic flux, for two values of parallel Zeeman field B ∥ = 0 (blue) and 1.2 meV (red).Different panels correspond to values of CS, starting from CS = 11.0 meV and increasing in steps of ∼ 0.3 meV from (a) to (i).There are two types of panels, shaded and non-shaded, corresponding to an even and odd number of occupied bands respectively.Thus, for finite magnetic field (red curves) only the non-shaded panels can develop MBSs, whereas for the shaded panels both blue and red curves correspond to the trivial case.Besides, W= 0.24 µm and Ln = 0.4 µm.

FIG. 14 .
FIG.14.Critical current as a function of the Zeeman field B ∥ for Φ = 0. We have used the same set of parameters as in Fig.13, panels (c) and (h) for the blue and red curves, respectively.The topological transition is evidenced by a sharp decrease of the critical current, see Ref.26.

FIG. 15 .
FIG. 15.Panels (a) and (b): Interedge electron-hole reflection probability T ↑,↑ he as a function of the deflection probability Γ d for the MBS (a) and the ZEABS (b), with ΓM = 0 and E = 0. We have used three different tunnel amplitude configurations t ↑ = 2t ↓ , t ↑ = t ↓ , t ↓ = 2t ↑ in the MBS case and ϵ0 = ∆0/2, ϵ0 = ∆0 and ϵ0 = 2∆0 for the ZEABS case.Panels (c) and (d): zero bias conductance as a function of Γ d for the MBS and ZEABS with the same parameters as in panels (a) and (b).

G = e 2 h 2 −
Tr[S † ee S ee ] + Tr[S † eh S eh ] ,(C21) = 2πν 0 (|t ↑ | 2 + |t ↓ | 2) and t σ is the effective cou-pling amplitude between the helical edge states and the Majorana bound states.The analytical expressions for S red,L/R is cumbersome and therefore, we obtain it only numerically.The same expressions can be rederived without using explicitly the Weidenmüller formula following Refs.48 and 94.