Planar Josephson Hall effect in topological Josephson junctions

Josephson junctions based on three-dimensional topological insulators offer intriguing possibilities to realize unconventional $p$-wave pairing and Majorana modes. Here, we provide a detailed study of the effect of a uniform magnetization in the normal region: We show how the interplay between the spin-momentum locking of the topological insulator and an in-plane magnetization parallel to the direction of phase bias leads to an asymmetry of the Andreev spectrum with respect to transverse momenta. If sufficiently large, this asymmetry induces a transition from a regime of gapless, counterpropagating Majorana modes to a regime with unprotected modes that are unidirectional at small transverse momenta. Intriguingly, the magnetization-induced asymmetry of the Andreev spectrum also gives rise to a Josephson Hall effect, that is, the appearance of a transverse Josephson current. The amplitude and current phase relation of the Josephson Hall current are studied in detail. In particular, we show how magnetic control and gating of the normal region can enable sizable Josephson Hall currents compared to the longitudinal Josephson current. Finally, we also propose in-plane magnetic fields as an alternative to the magnetization in the normal region and discuss how the planar Josephson Hall effect could be observed in experiments.


I. INTRODUCTION
The helical spin structure of the surface states of threedimensional topological insulators (3D TIs) offers intriguing possibilities of tailoring the surface-state properties by various proximity effects. A conventional s-wave superconductor can, for example, be used to proximityinduce superconductivity in the TI surface. The interplay between the helical spin-momentum locking of the TI surface state and the superconducting pairing then mediates an effective pairing between electrons at the Fermi level. This effective pairing features a mixture of singlet s-wave and triplet p-wave pair correlations 1-3 and turns the TI surface into a topological superconductor 2,4,6? -10 with Majorana zero modes 1 and odd-frequency pairing 11 .
Topological Josephson junctions are particularly intriguing if they are based on 3D TIs, as depicted in Fig. 1(a): Because of the 2D nature of the surface, the system supports modes that propagate along the direction parallel to the superconductor/normal TI interface, that is, the y direction in Fig. 1(a). Due to the protected zero-energy crossing occurring at zero transverse momentum and phase difference φ = π, a π-junction exhibits two counterpropagating, gapless states, so-called nonchiral Majorana modes 1 [see Fig. 1 Besides proximity-induced superconductivity, one can also envision other proximity effects whose interplay with the spin texture of the TI surface state leads to novel phenomena: In non-superconducting setups, for example, the interplay between the helical surface states and proximity-induced magnetism provides a versatile platform for studying fundamental effects and spintronic applications 4,34? ? , 35 . Ferromagnetic tunnel junctions based on 3D TIs 37,38,40,42,44? ? ? , in particular, show some promise for potential spintronic devices 45? . The combination of 3D TIs with both proximity-induced superconductivity and magnetism can prove even more interesting 12,[47][48][49] , however, and could point to novel possibilities for superconducting spintronics 50,51 .
Motivated by this prospect 86 as well as by phenomena found in non-superconducting TI tunneling junctions, such as the tunneling planar Hall effect 44 , we study 3D TI-based Josephson junctions with a ferromagnetic tunneling barrier [see Fig. 1 (c)]. In contrast to previous studies on this system 12,47,48 , we focus not only on the longitudinal response, but also on the transverse response to an applied phase bias. We find that especially the configuration with an in-plane magnetization parallel to the direction of the phase bias exhibits striking features: Such a magnetization leads to an asymmetric Andreev spectrum for a fixed finite transverse momentum. If sufficiently large, this asymmetry even induces a transition from the regime of counterpropagating, nonchiral Majorana modes to a regime with unprotected unidirectional modes at small transverse momenta [compare The Andreev spectrum is asymmetric and has been tilted in such a way that the two low-energy modes are unidirectional for small py. Note that these ABS are no Majorana modes protected against backscattering because there are additional zero-energy states for py close to the Fermi momentum (not shown). The asymmetry in the Andreev spectrum gives rise to a finite Josephson Hall current flowing in y direction.
Hall current [see Fig. 1 Below, we will discuss the origin of the Josephson Hall current, its properties and how it could be experimentally verified. The manuscript is organized as follows: After introducing the effective model used to describe the Josephson junction in Sec. II, we study its ABS in Sec. III. In Secs. IV and V, the procedure to compute the different Josephson currents is presented. These currents are then discussed in Sec. VI. A brief summary with a discussion of potential experimental measuring schemes concludes the manuscript in Sec. VII.

A. Hamiltonian and unitary transformation
In our model, we consider a Josephson junction based on the 2D surface state of a 3D TI, as depicted in Fig. 1(c), where the pairing in the superconducting (S) regions is induced from a nearby s-wave superconductor. The ferromagnetic (F) region is subject to an exchange splitting/Zeeman term proximity-induced from a nearby ferromagnet 55 . If one is only interested in an in-plane Zeeman term, an alternative way to realize such a Zeeman term is by applying an in-plane magnetic field as discussed in Sec. VII below. The surface state lies in the xy plane, with the direction of the superconducting phase bias denoted as the x direction. We take the system to be infinite in both the x and y directions. Here, we study the regime where the Fermi level is situated inside the bulk gap and where only surface states exist. Moreover, we assume that the surface considered is far enough away from the opposite surface so that there is no overlap between their states. Then, the Josephson junction based on a single surface is described by the Bogoliubov-de Gennes (BdG) Hamiltonian hole space, respectively. Moreover, σ = (σ x , σ y , σ z ) and unit matrices are not written explicitly in Eq. (1).
In this manuscript, we study two models for a Josephson junction with a F region of width d: (a) a model with a δ-like F region described by h(x) = dδ(x) and ∆(x) = ∆ and (b) a model with a finite F region where h(x) = Θ(d/2 − |x|) and ∆(x) = ∆Θ(|x| − d/2). In both cases, the phase convention is Φ(x) = φΘ(x), where φ is the superconducting phase difference between the two S regions. Furthermore, ∆ ≥ 0 is the proximity-induced superconducting pairing amplitude, v F the Fermi velocity of the surface state, and V 0 the potential in the F region, which can also be viewed as describing the difference between the chemical potentials in the S and F regions, µ and µ F = µ − V 0 . The Zeeman term due the proximityinduced ferromagnetic exchange splitting is described by the effective magnetization M = (M x , M y , M z ) 55 . Note that the direction of M is set by the magnetization in the ferromagnet.
For our calculations, it is more convenient to introduce the unitary rotation transformation in spin space U = (1 − iσ z )/ √ 2 and bring the Dirac Hamiltonian (1) into the form . Because of the rotated spin axes used in Eq. (2) M is a rotated effective magnetization, which is related to the components of the real magnetization M induced in the F region via M = (−M y , M x , M z ). From now on, we use the Hamiltonian (2) because it proves more convenient mathematically.

B. General form of the solutions
To solveĤ BdG Ψ(r) = EΨ(r) and obtain the eigenspectrum of Eq. (2), we first make use of translational invariance along the y direction, [Ĥ BdG ,p y ] = 0. Hence, we choose the ansatz Ψ(r) = e ipyy ψ(x)/ √ W , where p y is the momentum quantum number, ψ(x) is a spinor in Nambu space, and W is a unit width of the system in y direction. Here and in the remainder of this manuscript, we set = 1. The eigenenergies and ψ(x) can then be obtained from the 1D BdG equation whereĤ BdG (p y ) is given by Eq. (2) with the operatorp y replaced by the quantum number p y . The energy-momentum relation in the S regions is given by We find the following solutions in the S leads: where ξ = ±1 corresponds to particle-like and hole-like solutions and α = ±1 selects the direction of motion. Here, In the F region, the electron and hole states are given by

A. General equations
In order to understand the Josephson currents and the emergence of a Josephson Hall current, it is instructive to first look at the ABS of Eq. (3), that is, bound states decaying for |x| → ∞ and hence with energies |E| < ∆. We focus on the ABS of a junction with finite F region and refer to Appendix A for the Andreev spectrum of the δ-model, where relatively compact, analytical solutions are possible in certain limiting cases. The eigenenergies of the ABS and their corresponding eigenstates can be determined from the ansatz for a junction with a finite F region and s µ = sgn(µ). Now, the coefficients A 1 , A 2 , D e± , D h± , B 1 , B 2 have to be calculated from the boundary conditions at the S/F interfaces, The boundary conditions (8) lead to systems of linear equations for the coefficients A 1 to B 2 . By requiring a nontrivial solution of this system of linear equations, that is, by requiring its determinant to vanish, we find the ABS energies E = E(φ, p y ). in Fig. 2 for a short junction with a F region of length d = 330 nm, |µ| ∆, and different configurations of M . For these parameters, there are two ABS with energies E ± (φ, p y ) at a given momentum p y , where the subscript ± denotes which state lies higher (lower) in energy, that is, E + (φ, p y ) ≥ E − (φ, p y ). We can compare the φ and p y dependence of these ABS with the case of no magnetization, that is, M = 0 (not shown): For M = 0, the Andreev spectrum E ± (φ, p y ) exhibits a zero-energy crossing protected by fermion parity at odd integer multiples of φ = π and p y = 0, as also discussed in Appendix C. This protected zero-energy crossing is accompanied by two gapless, nonchiral Majorana modes that counterpropagate along the y direction and are localized mostly in the normal region 1,3 .
If we include a finite M , its effects on the ABS are the following: i) A component M y (not shown) shifts the entire Andreev spectrum as a function of φ, that is, E ± (φ, p y ) → E ± (φ + 2Z y , p y ), where Z y = M y d/v F , but leaves the spectrum otherwise unchanged. In particular, the protected zero-energy crossing for p y = 0 and the nonchiral Majorana modes are now shifted to φ = (2n + 1)π − 2Z y , where n ∈ Z is an integer.
ii) Finite components M x and M z , shown in Figs. 2(c,d) and (a,b) respectively, also do not remove this zero-energy crossing for p y = 0 and φ = (2n + 1)π − 2Z y . This crossing remains protected by the fermion parity and cannot be removed by a finite M x or M z 32,33 (see also Appendix C). The main effect of a finite out-of-plane magnetization M z in the F region is to detach the ABS from the continuum states with |E| > ∆ [see Fig. 2(a)], consistent with the results found in Refs. 47,48 .
iii) Intriguingly, we find that a finite M x = 0 introduces an asymmetry in the Andreev spectrum at finite p y as shown in Figs. 2(c) and (d): It does no longer satisfy E ± (φ, p y ) = −E ∓ (φ, p y ), but only the weaker condition E ± (φ, p y ) = −E ∓ (φ, −p y ), dictated by the particle-hole symmetry of the BdG formalism. This asymmetry of the Andreev spectrum emerges from the interplay between the spin-orbit coupling of the TI and M x , as we discuss in Appendix B with an effective low-energy model. In particular, Fig. 2(d), which shows the p y dependence of the Andreev spectrum, illustrates that the asymmetry E ± (φ, p y ) = −E ∓ (φ, p y ) manifests itself in a 'tilting' of the spectrum. If M x is large enough, it can even lead to a situation where the group velocities in y direction, v g ∝ ∂E ± (φ, p y )/∂p y , for ABS in the vicinity of p y = 0 and φ ≈ π − 2Z y have the same sign. Such a situation is shown in Fig. 2(d). In this regime, the ABS change from nonchiral, counterpropagating Majorana modes to modes propagating in the same direction for small p y . At small p y , the dispersion of these ABS is reminiscent of the unidirectional modes found in noncentrosymmetric superconductors [56][57][58] or in Rashba sandwiches 59 . It is important to note that the unidirectional ABS close to p y = 0 are, however, not protected against backscattering: As can be seen in Fig. 2(d), these states are accompanied by other zero-energy states with p y close to the Fermi momentum and with opposite group velocities.
In Fig. 2, we also compare the numerically obtained ABS with the analytical expressions one can derive for the ABS of a model with a δ-like F region in the Andreev approximation, as discussed in Appendix A. For short junctions and momenta close to p y = 0, these analytical expressions provide an excellent description of the ABS. In particular, these expressions also capture the asymmetry and 'tilting' of the Andreev spectrum induced by M x . Figure 3 shows the spatial dependence of the quasi- As can be discerned from Figs. 3(a-c), M z = 0 leads to ABS that are increasingly localized at the S/F interfaces as p y or M z are increased. One can understand this behavior by recalling that a magnetization component in z direction acts as a mass term that increasingly isolates the left and right S regions. If the two S regions are completely isolated from each other, that is, for M z → ∞, each S region separately corresponds to a topological superconductor that hosts one chiral Majorana mode at its boundary 1 . Hence, the results in Figs. 3(a-c) can be interpreted as the intermediate regime between M = 0 with nonchiral Majorana modes that are completely delocalized inside the F region and M z → ∞ with one chiral Majorana mode at each of the S/F interfaces. Comparing |ψ(x)| 2 for finite M z with |ψ(x)| 2 for a finite M x of the same strength, we find that |ψ(x)| 2 is not as localized at the S/F interfaces for M x , but rather constant in the whole F region, as shown in Figs. 3(d-f).
We also depict the expectation values of the spin densi- Fig. 3. The spin densities σ x (x) and σ y (x) are related to the currents in x and y directions, respectively (see Sec. IV below). By comparing right and left columns of Fig. 3, we find that for in-plane magnetization σ y (x) is delocalized within the F region. This is in contrast to the out-of-plane case, where the σ y (x) spin density and the wave functions are peaked near the S/F interfaces. Another important observation is that for finite M x the spin polarization amplitudes of the two Andreev levels are no longer equal. Together with the asymmetry of the Andreev spectrum for M x = 0 discussed above, a finite σ y (x) such as in Fig. 3 gives rise to a finite net Josephson Hall current, even for small M x . The emergence of this Josephson Hall current will be discussed next.

IV. CURRENT OPERATORS AND CONTINUITY EQUATIONS
Having found ABS with a peculiar behavior for M x = 0, we next study whether this gives characteristic signatures in observable quantities, such as for example the Josephson current. In order to derive current density operators, we consider the continuity equation for the charge density defined by the operator or equivalently by the matrix 1 2 eτ z σ 0 in the Nambu basis with e denoting the electron charge. The time evolution of the density operator is given by the equation of motion ∂ρ/∂t = i Ĥ BdG ,ρ(r, t) . After using the fermionic commutation relations for field operators, this equation of motion can be written in the form of the continuity equation Here, the quasiparticle part of the current density is proportional to the spin operator, analogous to the nonsuperconducting case for Dirac materials 44 The source term corresponding to the conversion of quasiparticles to Cooper pairs in the S leads is given bŷ The expectation values of these one-body operators can be expressed as traces of the Green's function which will be derived in the next section.

V. GREEN'S FUNCTION FORMALISM
In this section, we briefly describe the procedure for constructing the Green's function of the junction Hamiltonian (2). We follow the McMillan approach 60 and derive it from the wave-function solutions of the system. Thus, we choose the energies |E| > |∆|, where Eq. (4) describes propagating states, and solve the scattering problem where the multiindex n = (α, ξ) ∈ {n > } corresponds to an incident state from the left with fixed p y . Using the boundary conditions defined in Eq. (8) [or in Eq. (A1) for a δ-barrier], we find the reflection and transmission coefficients r nn (p y ) and t nn (p y ) correspondingly. Analogously, we can obtain states ψ n< (x) corresponding to processes when there is a quasiparticle incident from the right part of the junction. Furthermore, we employ the same procedure for the transposed HamiltonianĤ T BdG to find the conjugate stateŝ where the transpose operation acts on the Pauli matrices (Nambu space) and on the coordinate space (by replacinĝ p with −p). Afterward, we can write the Green's function for a fixed p y as an outer product of these solutions where the position-independent coefficients C nn should be determined from the boundary condition at x = x , Having determined the Green's function of the system in this way, we can express a given single-particle operator in terms of this Green's function and obtain the expectation value of the operator by evaluating a sum over fermionic Matsubara frequencies ω n . In order to do so, we perform an analytical continuation of all expressions in Eq. (15) to the complex plane with E → iω n . Furthermore, we use the fact that the retarded (advanced) Green's function is analytical in the upper (lower) half of the complex plane. To access negative Matsubara frequencies, we calculate the advanced Green's function in the same manner as the retarded one.
Finally, the expectation value of the quasiparticle part of the current density operator is given by with l = x, y 87 . If M x = 0, the summation in frequency space for j y (x) does not converge due to an oscillating behavior at high energies. This is similar to the behavior of j y (x) in the normal state, where the contributions arising from the oscillating wave functions for M x = 0 vanish only after integration over x, that is, when computing the transverse current from the transverse current density. Such a behavior can also be understood as an artifact of the continuum Dirac model. In fact, this model is only valid close to the Dirac point within the band gap of the TI. To account for this, we separate the current contributions into superconducting and normal parts, j = j SC + j N , where we define j SC = j − j N . Here, j N is evaluated for a normal system where we set ∆ = 0 and captures all divergent terms that we treat in more details in Appendix D. In the remaining expression j SC , which is also the part that does not vanish after integration over x, the sum converges fast and is performed numerically up to a cutoff. Since it can be proven that the normal part goes to zero in equilibrium, we focus only on the regular part j SC which describes the actual Josephson current in the junction.
Note that Eq. (17) only contains the spatial dependence of the quasiparticle part of the current density. In order to compute the spatial dependence of the full current density, one also needs to include contributions due to the source termŜ(x) from Eq. (12) in the S leads 62 . As a consequence the full current density in x direction, consisting of j x (x) from Eq. (17) and a term originating fromŜ(x) in the S regions, has a constant value and is independent of the position x. For the transverse current, there is no contribution due toŜ(x). Finally, we remark that the current densities computed from Eqs. (15) and (17) are the current densities for a situation where all states have equilibrium occupations without any external constraints. Therefore, Eqs. (15) and (17) describe the current densities without conservation of the fermion parity.

VI. JOSEPHSON HALL CURRENT AND CURRENT-PHASE RELATION
We are now in a position to discuss the emergence of the transverse Josephson Hall current, which is the main result of this manuscript. Without a barrier magnetization, M = 0, or if there is only an M y component of M , the transverse current density j y (x) = ĵ y (x) vanishes. On the other hand, the asymmetry in the Andreev spectrum due to a finite M x or the separation of Majorana modes localized at the S/F interfaces due to a finite M z induce a finite j y (x). This is illustrated by Fig. 4, where we present the spatial dependence of j y (x) in the presence of a finite magnetization in the barrier. For a magnetization M z [ Fig. 4(a)], we observe two transverse current densities of opposite sign localized at the S/F interfaces. At each interface, this localized current density corresponds mainly to the chiral Majorana mode that emerges at an S/F interface for large M z as discussed above in Sec. III B. The magnitude of j y (x) increases proportional to M z . In contrast to the constant longitudinal Josephson current density, j y (x) oscillates with k F and decays exponentially into the S regions. As shown from a symmetry argument in Appendix E, j y (x) is odd with respect x and consequently the total Josephson Hall current through the F region, is zero for finite M z . For a magnetization M x , there is a finite transverse Josephson current density flowing in the same direction inside the whole F region, as shown in Fig. 4(b). In this case, the current density profile j y (x) is an even function of x, which clearly allows for a finite Josephson Hall current I y as given by Eq. (18) flowing in the F region. To increase I y , one can apply an additional gate voltage V 0 inside the barrier, which reduces the effective Fermi momentum in the F region and hence suppresses the oscillating behavior inside the barrier. In Fig. 4(c), one can see that by tuning V 0 close to µ we can achieve an almost flat profile of j y (x) within the junction, thereby increasing I y .
For the case of M x = 0 we, moreover, compare I y with the corresponding longitudinal Josephson current I x in Fig. 5. The latter one is normalized to the same cross section, by multiplying with the factor d/W . Figure 5(a) shows the current-phase relation of I x and I y for several different values of V 0 . Both I x and I y are 2π-periodic in the superconducting phase difference φ since fermion parity is not conserved if all states have equilibrium occupations (see Sec. V). There is, however, a marked difference in the current-phase relation between the non-sinusoidal I x , which is an odd function of φ, and I y , which is an even function of φ. Unlike I x , I y does typically not exhibit zeros at integer multiples of φ = π. Remarkably, we see that for φ close to π the direction of the current can be controlled not only by the sign of M x , but also by modifying the gate voltage V 0 , which can be appealing for practical applications. In Fig. 5(b), we show I x and I y at φ = 0.7π as a function of M x . This illustrates that for a large enough magnetization the Josephson Hall current can exceed the longitudinal Josephson current. Such ratios I y /I x > 1 are comparable to the ratios found in normal TI-based ferromagnetic tunneling junctions and are a result of the strong SOC in 3D TI surface state. This makes the planar Josephson Hall effect in TI-based Josephson junctions a promising candidate to observe sizable transverse currents with ratios I y /I x exceeding the corresponding ratios of other Josephson Hall effects 53,54,6388 .

VII. CONCLUSIONS AND EXPERIMENTAL PROPOSALS
In this manuscript, we have studied Josephson junctions realized on three-dimensional topological insulators which are subject to a Zeeman term in the normal topological insulator region. Most importantly, we have found that the interplay between the spin-momentum locking of the topological insulator surface state, superconductivity and an in-plane Zeeman field in the normal region gives rise to a net transverse Josephson Hall current. For this Josephson Hall current to emerge, the in-plane Zeeman field has to have a component parallel to the superconducting phase bias direction [see Fig. 1(d)]. Since the effect is caused by an in-plane Zeeman term, we refer to it as the planar Josephson Hall effect to also distinguish it from other Josephson Hall effects 53,54,63 .
The emergence of the Josephson Hall current is reflected in an asymmetry and 'tilting' of the Andreev spectrum with respect to the transverse momenta p y . If sufficiently large, this asymmetry even induces a transition in the Andreev spectrum from a regime with gapless, counterpropagating Majorana modes to a regime with unprotected modes that are unidirectional at small p y . Due to strong spin-orbit coupling, the planar Josephson Hall effect in topological-insulator-based junctions enables sizable Josephson Hall currents, whose amplitudes can be further modulated by electrostatic and/or magnetic control of the normal region.
Until now, we have mainly discussed Zeeman terms induced into the normal topological insulator region by magnetic proximity effects from a nearby ferromagnet, such as in YIG/(Bi,Sb) 2 Te 3 65? , EuS/Bi 2 Se 3 67 or (Bi,Mn)Te with thin Fe overlayers 68 . Since the planar Josephson Hall effect requires in-plane Zeeman terms, an alternative realization could be by applying an in-plane magnetic field along the phase bias direction in the normal region 89 . Assuming, for example, an in-plane g factor of g = 10, an in-plane magnetic field of around B = 0.35 T corresponds to a Zeeman splitting of 0.1 meV 90 , which can already yield sizable Josephson Hall currents flowing through the normal region, as illustrated by Fig. 5(b). Indeed, in Josephson junctions composed of thin-film aluminium and HgTe quantum wells, which can also act as three-dimensional topological insulators 71 , in-plane magnetic fields of more than 1 T have been achieved [72][73][74] .
The planar Josephson Hall effect could then be experimentally verified by attaching transverse leads to the normal F region of the junction, as depicted in Fig. 6. If these leads are normal leads and act as voltage probes, that is, under open circuit conditions in the y direction, an experimentally detectable Hall voltage between the leads emerges instead of the Josephson Hall current. Assuming that the resistance arises only due to the contact resistances at the interfaces between the F region and the leads, we can obtain a rough estimate of this Hall voltage as V H = R c I y , where R c is the total contact resistance due to both leads. This expression also contains the Josephson Hall current I y that would flow through the F region if the leads did not act as voltage probes. Alternatively, the planar Josephson Hall effect could be experimentally verified by replacing the transverse normal leads with superconducting leads S L1 and S L2 75 . Then,  1) and (2). The advantage of this model is that it allows for a transparent analytical treatment of the ABS with relatively compact expressions.

3D TI surface state
For a δ-junction, the ansatz to obtain the ABS is similar to Eq. (7), with the states ψ(x < 0) given by the first line of Eq. (7) and ψ(x > 0) given by the third line of Eq. (7). Now, the coefficients A 1 , A 2 , B 1 , and B 2 have to be calculated from the boundary conditions at x = 0. This boundary condition can be obtained by integrating Eq. (3) from x = −η to x = η with η → 0 + . The corresponding procedure 17,44,81 yields wherê 2. δ-model at py = 0 First, we look at the case of p y = 0, where v F (αp ξ + ip y )/(µ + ξΩ) = α. We invoke the boundary condition (A1) on the first and third lines of Eq. (7) and require a nontrivial solution of the resulting system of linear equations. This then yields the two ABS energies E = PE 0 (φ), where P = ±1 denotes the two fermionparity branches and E 0 (φ) = ∆ cos (φ/2 + Z y ) The two ABS given by E = ±E 0 (φ) exhibit a nondegenerate zero-energy crossing at φ = π if M = 0 92 . At this crossing, the ground-state fermion parity changes, and the two branches in Eq. (A3) have been chosen such that each branch preserves its fermion parity 32,33 . As such a non-degenerate zero-energy crossing is protected by the fermion parity, it cannot be removed even for finite M = 0 (see Refs. 32,33 and Appendix C). The crossing can only be shifted, which is what happens for a finite M y = 0, where E 0 (φ) = 0 for φ = (2n + 1)π − 2Z y with n ∈ Z. This protected crossing is a hallmark of the topological Josephson junction and can also be found in models with finite F region. At p y = 0, the main effect of magnetization components M x,z = 0 is thus to reduce the bandwidth of the ABS and detach them from the continuum states. We also remark that the case of p y = 0 is equivalent to a Josephson junction based on a single quantum spin Hall edge if M x → M y , M y → −M z and M z → M x . With these replacements, Eq. (A3) describes the ABS spectrum of such Josephson junctions in the short junction regime 93 .

δ-model in Andreev approximation
Another limit that allows for closed analytical solutions is the case of |µ| ∆, where we can make use of the Andreev approximation. If we introduce the angle −π/2 < θ < π/2 via v F p y = µ sin θ, the eigenstates (4) are simplified within the Andreev approximation in so far that With these approximations, the condition for a nontrivial solution to Eq. (3) can be written as where X = √ ∆ 2 − E 2 /E and From the two solutions of Eq. (A5), X = A(θ) ± A 2 (θ) + B(θ), one can see that at a fixed angle θ the two solutions for the energy E ± (φ, θ) do not come as E ± (φ, θ) = −E ∓ (φ, θ) if A(θ) = 0. This is the case for finite Z x and finite θ. Instead, the two solutions can be obtained as which only satisfies the weaker condition E ± (φ, θ) = −E ∓ (φ, −θ) originating from the particle-hole symmetry of the formalism.
If θ = 0, Eq. (A7) reduces simply to E ± (φ, θ = 0) = ±|E 0 (φ)| with E 0 (φ) given by Eq. (A3). Note that now the sign ± does not refer to the parity branch, but instead to positive and negative energies. Another point worth mentioning with regard to Eq. (A7) is that for Z x = 0 it reduces to 3,84 is the transmission of a normal/ferromagnet/normal junction with Z x = 0 44 .
For M x = 0, Eq. (A7) exhibits several salient features: A finite M x = 0 introduces not only an asymmetry in the ABS spectrum at finite p y , but can even lead to a situation where the group velocities in y direction, v g ∝ ∂E ± (φ, θ)/∂θ, have the same sign for ABS in the vicinity of p y = 0 and φ ≈ π − 2Z y . At these momenta, the two ABS propagate in the same direction. This change from nonchiral, counterpropagating ABS to unidirectional ABS propagating in the same direction occurs for B(θ) < 0. Close to φ ≈ π − 2Z y , B(θ) < 0 is satisfied if |M x | > |V 0 |. Hence, if the Zeeman term in the direction of the phase bias φ exceeds the mismatch between the chemical potentials of the S and F regions, the ABS close to p y = 0 and φ + 2Z y ≈ π propagate in the same direction in short junctions.
At this point, it is important to remark that the Andreev approximation (A4) breaks down at large transverse momenta, that is, at momenta close to the Fermi momentum p F . Because of this, Eq. (A7) does not describe the ABS for p y ≈ p F well. This is also illustrated by Fig. 2(d), which shows a comparison between Eq. (A7) and the results for a finite barrier without any further approximations. For small p y = µ sin θ, Eq. (A7) is in good agreement with the results of the finite barrier. Equation (A7) cannot, however, capture the appearance of zero-energy ABS that occur at large momenta once the modes close to p y = 0 become unidirectional.

Appendix B: Effective low-energy model
The asymmetry of the ABS spectrum as well as the emergence of unidirectional modes for large M x and small p y can be understood from the interplay between the effective spin degree of freedom and M x . To elucidate the origin of these modes, we employ a simple effective low-energy Hamiltonian. For a δ-like F region, the BdG Hamiltonian (2) always supports two ABS. Following the procedure in Ref. 1 , we derive an effective low-energy Hamiltonian describing these two ABS in the vicinity of the protected crossing at φ = π − 2Z y for small p y .
To do so, we first note that the BdG Hamiltonian (2) can be written asĤ BdG (p y ) =Ĥ BdG (p y = 0)+v F p y σ y τ z , where we treat the term v F p y σ y τ z as a perturbation. Then, we can take the two parity-conserving ABS |± for p y = 0 discussed in Sec. A 2 and project the full Hamiltonian (2) onto these two states. This procedure yields the effective Hamiltonian where E 0 (φ) is given by Eq. (A3) and originates from H BdG (p y = 0). In Eq. (B1), the two-level system formed by the two ABS at p y = 0 is described by the Pauli matricesσ l (l = x, y, z) and the corresponding 2 × 2 unit matrixσ 0 . Moreover, we have introduced the velocities which arise from the matrix elements of the perturbation v F p y σ y τ z . Since we are mainly interested in the ABS close to the crossing at φ = π − 2Z y , we have ap- The spectrum of Eq. (B1) is given by E ± eff (φ) = v y p y ± E 2 0 (φ) + (v 0 p y ) 2 . At the crossing point, E 0 (π − 2Z y ) = 0 and E ± eff (π − 2Z y ) = (v y ± v 0 )p y . If v y = 0, that is, if M x = 0, the spectrum at φ = π − 2Z y is simply E ± eff (π − 2Z y ) = ±v 0 p y and describes two counterpropagating Majorana modes along the y direction, similar to Ref. 1 . For finite v y , on the other hand, the group velocities (v y ± v 0 ) of the two modes point into the same direction if |v y | > |v 0 |.
The appearance of a term v y p yσ0 in Eq. (B1) is thus the origin of the unidirectional modes at small p y . While there is always a finite v 0 in TI-based Josephson junctions, v y = 0 only arises for finite M x = 0. This can be understood in the following way: The terms containing v 0 and v y originate from the matrix elements v F p y P|σ y τ z |P with P, P = ±1 denoting the two parity branches of p y = 0. If M x = 0, the effective spin orientation of the eigenstates |± ofĤ BdG (p y = 0) lie in the xz plane and thus the expectation values ±|σ y τ z |± vanish and v y = 0. Only off-diagonal matrix elements ∓|σ y τ z |± are finite and give rise to v 0 = 0.
For finite M x = 0, however, the effective spin expectation values of |± now also acquire a component in the y direction and ±|σ y τ z |± = 0. The eigenstates |± satisfy the relation |± =K|∓ , whereK denotes complex conjugation. Because of this property, +|σ y τ z |+ = −|σ y τ z |− and consequently the diagonal matrix elements of the perturbation is proportional toσ 0 (and not toσ z or a linear combination ofσ 0 and σ z ). The spectrum of Eq. (B1) makes it clear that the ABS spectrum for small p y and close to the protected crossing point φ + 2Z y = π (or, more generally, close to φ = (2n + 1)π − 2Z y with n ∈ Z) can support unidirectional modes around p y ≈ 0 for finite M x .
Appendix C: Protected zero-energy crossing for py = 0 A peculiar feature of the ABS spectrum of a Josephson junction based on a single surface of a 3D TI is its protected zero-energy crossing for p y = 0, even in the presence of a Zeeman term in the F region, as discussed in Sec. A 2. Following Ref. 33 , we can understand this protection arising from the particle-hole symmetry of the BdG Hamiltonian, which allows one to define a Pfaffian, Pf Ĥ BdG (p y = 0) for any φ. The existence of a Pfaffian then implies that two-fold degenerate zero-energy states are generically protected against perturbations as long as particle-hole symmetry is preserved.
For the system studied here, particle-hole symmetry is described by the operatorĈ = σ y τ yK , whereK denotes complex conjugation and σ y and τ y are the respective Pauli matrices in spin and electron/hole space. Any BdG Hamiltonian, including Eq. (2), anticommutes withĈ, If we introduce the momentum quantum number p y , this becomesĈĤ Thus, only for p y = 0, does particle-hole symmetry imply Ĉ ,Ĥ BdG (p y = 0) = 0, while in general particle-hole symmetry connects states with p y to states with −p y . From now on, we focus only on the two ABS |± at p y = 0 and with M = 0. For a δ-like F region and M = 0, the energies are simply given by E = ±∆ cos(φ/2), that is, they possess two-fold degenerate zero-energy states at φ = π. Similar to Sec. B, the corresponding low-energy Hamiltonian is the 2×2 matrix with respect to the ABS |± , which can in turn be transformed tỗ The Pfaffian of Eq. (C3) is then given by Pf Ĥ 0 eff = iPf Â 0 eff = i∆ cos(φ/2) and can be related to the ground-state fermion parity F 0 via (−1) F0 = sgn Pf Â 0 eff . Since Pf Ĥ 0 eff exhibits only a single zero, a perturbation that preserves particle-hole symmetry cannot remove the two zero-energy states, but only shift them to other values of φ 33 . Now, we are in a position to understand why the crossing at φ = π is protected against finite M in the F region. For finite M and p y , we can write the BdG Hamiltonian asĤ We remind the reader that because of the rotated spin axes used in Eq. On the other hand, Ĉ ,Ĥ py = 0 and thus a gap is opened at finite p y because in this case particle-hole symmetry does not protectĤ py , but connectsĤ py andĤ −py (see above). While we have employed this analysis to the case of a δ-barrier for illustration, we note that this is valid for all single-energy crossings that are only double degenerate, including the case of a finite barrier also studied in this manuscript 94 . Hence, as a final remark we note that the analysis from Eqs. (C3) and (C4) applies also to the case of finite M , where ∆ cos(φ/2) should simply be replaced by E 0 (φ) from Eq. (A3). This also makes it clear that the groundstate fermion parity F 0 given by (−1) F0 = sgn [E 0 (φ)] = sgn [cos(φ/2 + Z y )] for finite M is only shifted in its φ dependence by Z y ∝ M y , but remains unaltered otherwise.
To obtain the Green's function, we proceed analogously to the main text. For example, the lead solutions are where v F q 0 = (µ + E) 2 − (v F p y ) 2 and α = ±1 gives the direction of propagation. In this case, we have two helical counterpropagating sates for each p y . The F-barrier solution is given by Eq. (6) with ξ = +1. We omit the details of solving the transposed Hamiltonian and deriving the scattering states.
To analyze the current expectation value, we stay within the real-energy picture because it allows for a discussion of high-energy contributions and the continuation to the Matsubara frequencies when the function does not decay fast for |E| → ∞. The current operator simplifies to j N = ev F σ and we obtain where n(E) is the Fermi-Dirac distribution. In the case |E + µ| ≥ |p y |, we can conduct a variable substitution in the integral, namely v F q 0 = |E + µ| cos θ and v F p y = |E + µ| sin θ. This allows us to rewrite the integration limits over p y , which yields for the transverse current where r(E, θ) is the reflection coefficient of the mode incident from the left lead. The case of |E + µ| < |p y | is treated analogously employing hyperbolic functions. In the non-superconducting case, it is possible to obtain a relatively compact form for the reflection coefficient 44 . a. δ-barrier solutions We consider the stationary states similar to Eq. (13), with the superconducting lead wave functions replaced by ψ (0) n (x) and in the limit d → 0. Then, using the boundary condition (A1) for the electron block, we obtain the reflection and transmission coefficients r = e isθ (Z x + iZ z cos θ − sZ 0 sin θ) sin Z Z cos θ cos Z + i(Z 0 − sZ x sin θ) sin Z , t = e −iZy Z cos θ Z cos θ cos Z + i(Z 0 − sZ x sin θ) sin Z , where we have defined s = sgn(E + µ). We notice the property r(E, θ) = r(−E, −θ). The δ-barrier does not introduce an energy scale. Therefore, the reflection amplitude is energy independent. This means that all states in the system are affected by the introduction of the barrier, which has significant consequences for Eq. (D4) because the spectrum is not bounded from below.
b. Finite-barrier solutions By using the boundary conditions from Eq. (8) and lead wave functions defined in Eq. (D2), we obtain the reflection and transmission coefficients as r(E) = e isθ−id(E+µ) cos θ/v F a(E) sin(dk 0 ) v F k 0 cos θ cos(dk 0 ) + ib(E) sin(dk 0 ) , (D7) with s = sgn(E + µ) and a(E) =M x + iM z cos θ − sV 0 sin θ, In this case, r(E) exhibits a behavior ∼ 1/ |E| for large |E|, but this is still not enough to make the energy integral in Eq. (D4) finite. The divergent behavior comes from the fact that the energy spectrum of the Dirac cone is not bound from below, so formally we have to include all contributions down to E = −∞ in the expectation values. At the same time, the presence of the magnetic barrier introduces spin polarization into all states, making them contribute to the integral (D4). In the real system, the low-energy model is invalid for energies far from the Dirac cone located in the gap of the topological insulator. On the other hand, the high-energy solutions become highly oscillating with wave vector E/v F . These oscillations cannot be resolved in the real system, which provides another argument why we should drop high-energy terms. Thus, we choose to use the regularization e −λ|E| in the integral. Then, we can perform the energy integration analytically which yields a prefactor λ in front of the expression for the current. Hence, after taking the limit λ → 0, j N y vanishes. When computing the current density in the superconducting case (∆ > 0), we subtract j N expression before performing integration. After that, the integral can be performed numerically or more conveniently by going to the complex plane and mapping it to the Matsubara sum.
which act in the spin space such that M y,z → −M y,z and M x,y → −M x,y , respectively and both reverse the sign of the kinetic term. If M x(z) = 0, we find that S x(z) = M yz(xy) IT is a symmetry of the Hamiltonian.
Next, we derive how the current operator transforms under given symmetries S x j y S −1 x = −j y and S z j y S −1 z = j y .
Using that Sψ(r) can be presented as U ψ * (V r), where U is a unitary matrix in spinor space and V is an orthogonal transformation in coordinate space, we obtain a relation for the contribution of the operator expectation value from a single state Sψ(r)| j y |Sψ(r) = ψ(V r)| Sj y S −1 |ψ(V r) , (E2) where the scalar product is performed only in spinor space. The expectation value of the total current density is the sum of contributions from all states weighted with the occupation number, which is a function of energy. If S is the symmetry of H BdG , states |ψ n (x, y) and S |ψ n (x, y) either have the same energy or coincide. Hence, we obtain Since the current is independent of y due to translation invariance, these symmetry relation are generalized to the whole junction. We also mention that in case of the δ-barrier we may have a discontinuity at x = 0 and the value of the current would depend on the direction from which we approach the barrier. puted as gµBB/2 = 0.1 meV, where µB is the Bohr magneton. 91 The Zeeman field needed to induce the planar Josephson Hall effect acts as an in-plane Zeeman field that is perpendicular to direction of phase bias in the transverse SL1/F/SL2 Josephson junction. 92 Note that in contrast to the main text and to Appendix A 3, the sign ± refers to the parity branch, not to the ordering of the energies. 93 Compare to Ref. 32 discussing Jospehson junctions based on a single quantum spin Hall edge. 94 In contrast, such zero-energy crossings are not protected for a two-dimensional electron gas without magnetic field because the crossings are four-fold degenerate due to the spin degeneracy.