Planar Josephson junctions in germanium: Effect of cubic spin-orbit interaction

Planar Josephson junctions comprising semiconductors with strong spin-orbit interaction (SOI) are promising platforms to host Majorana bound states (MBSs). Previous works on MBSs in planar Josephson junctions have focused on electron gases, where SOI is linear in momentum. In contrast, a two-dimensional hole gas in planar germanium (Ge) exhibits SOI that is cubic in momentum. Nevertheless, we show here that due to the particularly large SOI, Ge is a favorable material. Using a discretized model, we numerically simulate a Ge planar Josephson junction and demonstrate that also cubic SOI can lead to the emergence of MBSs. Interestingly, we find that the cubic SOI yields an asymmetric phase diagram as a function of the superconducting phase difference across the junction. We also find that trivial Andreev bound states can imitate the signatures of MBSs in a Ge planar Josephson junction, therefore making the experimental detection of MBSs difficult. We use experimentally realistic parameters to assess if the topological phase is accessible within experimental limitations. Our analysis shows that two-dimensional Ge is an auspicious candidate for topological phases.

Figure 1

A germanium 2DHG (gray box) of size $L_x\times L_y$ is proximitized by two superconductors (blue boxes) of width $S$. The superconductors induce superconducting gaps in germanium of strength ${\mathrm{\Delta }}_{\text{SC}}^l$ and ${\mathrm{\Delta }}_{\text{SC}}^r$ with a phase difference $\varphi$ between them. They are separated by a normal section of width $W$. A magnetic field $\vec{B}$ is applied along the junction. In the topological phase, MBSs are localized at the ends of the normal region (indicated schematically by the red regions) (see Sec. 3).

Figure 2

Sketch of the germanium heterostructure. A thin germanium layer (gray) of thickness $l$ is sandwiched between two layers of ${\mathrm{Si}}_{1-x}{\mathrm{Ge}}_x$ (green) and an electric field of strength $\mathcal{E}$ is applied along the $z$ axis.

Figure 3

(a) Spectrum $E(k_x,k_y)$ of the 2DHG. (b) Fermi surfaces, i.e., the momenta where $E(k_x,k_y)=0$. The solid black lines in (a) and (b) correspond to the two lowest-energy bands of the Luttinger-Kohn (LK) Hamiltonian confined to two dimensions. The red dashed lines correspond to the energy eigenvalues of the effective Hamiltonian $H_{\text{eff}}+H_4$ defined in Eqs. (6) and (7). The two models agree best at small momenta. (c) Spectrum $E(k_x,k_y=0)$ of the effective model with $\beta =0$ (solid blue curve) and $\beta =0.02E_{\mathcal{E}}{l}_{\mathcal{E}}^4$ (dashed red curve). Without the quartic term, the lower-energy band bends downwards at large momenta, causing an unphysical Fermi surface. By introducing the quartic term defined in Eq. (7), this is avoided. The quartic term does not change the low-energy physics noticeably, as is shown in the inset of (c). The parameters for the Luttinger-Kohn Hamiltonian are $E_s=0$, ${\gamma }_1=13.38$, and ${\gamma }_s=4.965$ [103]. Due to the normalization, the spectrum is independent of the electric field $\mathcal{E}$. The parameters for the effective model are $\mu =0.3E_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.096E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\alpha =0.133E_{\mathcal{E}}{l}_{\mathcal{E}}^3$. These parameters were determined using perturbation theory [103].

Figure 4

(a)–(c) Topological invariant $Q_{\mathbb{Z}}$ of the BDI symmetry class, the bulk gap ${\mathrm{\Delta }}_g$ [defined in Eq. (23)], and the lowest-energy level $E_0$ (calculated in the finite geometry) as functions of the Zeeman energy ${\mathrm{\Delta }}_Z$ and the superconducting phase difference $\varphi$. The black cross in (b) indicates the maximum bulk gap ${\mathrm{\Delta }}_m$ in the topological phase. The black curves in (b) and (c) indicate where the bulk gap closes at $k_x=0$. (d) Probability distribution ${\left|\mathrm{\Psi }\right|}^2$ of a MBS (in arbitrary units). Boundaries between regions with different values of $Q_{\mathbb{Z}}$ in (a) coincide with a closing of the bulk gap, i.e., ${\mathrm{\Delta }}_g=0$, in (b). If the parity of $Q_{\mathbb{Z}}$ changes, the bulk gap closes at $k_x=0$, whereas otherwise it closes at some finite $k_x\ne 0$. Regions where $Q_{\mathbb{Z}}\ne 0$ coincide with regions where the lowest-energy level $E_0$ is close to zero in (c). These bound states are MBSs. Used parameters: $L_x=199.5l_{\mathcal{E}}$ [in (c) and (d)], $L_y=49.5l_{\mathcal{E}}$, $W=4.5l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\alpha =0.39E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, $\mu =0.3E_{\mathcal{E}}$, ${\mathrm{\Delta }}_{\text{SC}}=\mu /3$. In (d), $\varphi =1.3\pi$ and ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=1.4$ [indicated by the white cross in (c)].

Figure 5

Fitted slope of the ${\mathrm{\Delta }}_Z$ vs $\varphi$ phase transition curve for different strengths $\alpha$ of CSOI. An example of the fit is shown in the inset, where the phase transition curve for $\alpha =0.34E_{\mathcal{E}}{l}_{\mathcal{E}}^3$ (indicated by the vertical dashed black line) and the fitted lines are shown. The larger the CSOI strength, the more the left and right slopes differ, therefore, the more asymmetric the phase diagram becomes. The parameters are $L_y=49.5l_{\mathcal{E}}$, $W=4.5l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\mu =0.3E_{\mathcal{E}}$, ${\mathrm{\Delta }}_{\text{SC}}=\mu /3$. The strength $\beta$ of the quartic term is varied together with the CSOI strength $\alpha$ as follows. For several values of $\alpha$, the corresponding value of $\beta$ is determined following Appendix pp2. Using a quadratic fit, we then set $\beta =c_0+c_1\alpha +c_2{\alpha }^2$ with $c_0=0.0017E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, $c_1=0.0236l_{\mathcal{E}}$, and $c_2=0.8507E_{\mathcal{E}}^{-1}{l}_{\mathcal{E}}^{-2}$.

Figure 6

Largest bulk gap in the topological phase ${\mathrm{\Delta }}_m$ as a function of the normal section width $W$. As the normal section becomes wider, ${\mathrm{\Delta }}_m$ decreases. The parameters are $N_y=49.5l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.14E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\alpha =0.39E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, $\mu =0.3E_{\mathcal{E}}$, ${\mathrm{\Delta }}_{\text{SC}}=\mu /3$.

Figure 7

(a), (e) Energy spectrum of a planar JJ as a function of the Zeeman energy ${\mathrm{\Delta }}_Z$. The black dots are the energy levels of the system calculated numerically in the finite geometry. The colored lines indicate the bulk gap ${\mathrm{\Delta }}_g$, defined in Eq. (23). The line color indicates the value of $k_x$ for which $E\left(k_x\right)={\mathrm{\Delta }}_g$. In (a), the system has the effective time-reversal symmetry defined in Eq. (21). The background color in (a) signifies the corresponding BDI topological invariant $Q_{\mathbb{Z}}$. In the white region $Q_{\mathbb{Z}}=0$ (trivial phase, no zero-energy states), in the green region $Q_{\mathbb{Z}}=-1$ (two zero-energy states, i.e., one pair of MBSs), and in the blue region $Q_{\mathbb{Z}}=-2$ (four zero-energy states, i.e., two pairs of MBSs). In (e) the effective time-reversal symmetry is broken by having ${\mathrm{\Delta }}_{\text{SC}}^l\ne {\mathrm{\Delta }}_{\text{SC}}^r$. In the finite geometry we also assume disorder in the chemical potential of the system (see Appendix pp4-s2). The background color signifies the corresponding topological invariant $Q_{{\mathbb{Z}}_2}$ of symmetry class D. In the white region $Q_{{\mathbb{Z}}_2}=1$ (trivial phase, no MBSs) and in the gray region $Q_{{\mathbb{Z}}_2}=-1$ (topological phase, one pair of MBSs). (b)–(d) and (f)–(h) Probability distributions ${\left|\mathrm{\Psi }\right|}^2$ of the MBSs [(b), (c), (d), and (h)] and ABSs [(f) and (g)] (arbitrary units). Although the two pairs of MBSs that appear in the blue region of (a) hybridize when the symmetry is broken, the resulting ABSs are still localized in both $x$ and $y$ directions and their energies are close to zero. The parameters used are $L_x=349.5l_{\mathcal{E}}$, $L_y=49.5l_{\mathcal{E}}$, $W=4.5l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\mu =0.3E_{\mathcal{E}}$, $\alpha =0.39E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, $\varphi =0$. In (a)–(d) ${\mathrm{\Delta }}_{\text{SC}}^l={\mathrm{\Delta }}_{\text{SC}}^r=\mu /3$, in (e)–(h) ${\mathrm{\Delta }}_{\text{SC}}^l=\mu /3$ and ${\mathrm{\Delta }}_{\text{SC}}^r=0.6{\mathrm{\Delta }}_{\text{SC}}^l$. In (b), (c), (f), and (g) ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}^l=1.7$ and in (d) and (h) ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}^l=2.08$ [these values are indicated by the red crosses in (a) and (e)].

Figure 8

The normal section width $W\left(x\right)$ of a planar JJ is varied according to Eq. (24). (a) Energy spectrum (calculated in the finite geometry) as a function of the Zeeman energy ${\mathrm{\Delta }}_Z$. The blue vertical dashed lines indicate ${\mathrm{\Delta }}_{Z,l}$ and ${\mathrm{\Delta }}_{Z,r}$. (b), (c) Probability distribution ${\left|\mathrm{\Psi }\right|}^2$ (arbitrary units) of the lowest-energy states. The horizontal dashed lines indicate the boundaries of the superconductors. (d) Lowest-energy level $E_0$ (calculated in the finite geometry) as a function of the Zeeman energy ${\mathrm{\Delta }}_Z$ and the superconducting phase difference $\varphi$. The black dashed (solid) curve indicates ${\mathrm{\Delta }}_{Z,l}$ (${\mathrm{\Delta }}_{Z,r}$). In the region of parameter space where ${\mathrm{\Delta }}_Z\in [{\mathrm{\Delta }}_{Z,l},{\mathrm{\Delta }}_{Z,r}]$, near-zero-energy ABSs appear [in (b)]. The ABSs are localized to the left side of the junction only. Once the bulk of the system is in the topological phase, MBSs appear [in (c)] and they are localized on both ends of the junction. The parameters used are $L_x=149.5l_{\mathcal{E}}$, $L_y=99.5l_{\mathcal{E}}$ in (a)–(c) and $L_y=49.5l_{\mathcal{E}}$ in (d), $W_{-\infty }=4.5l_{\mathcal{E}}$, $W_{\infty }=1.5l_{\mathcal{E}}$, $\sigma =17.5l_{\mathcal{E}}$, $p=21l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\mu =0.3E_{\mathcal{E}}$, ${\mathrm{\Delta }}_{\text{SC}}=\mu /3$, $\varphi =0.8\pi$, $\alpha =0.39E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=1.1$ [in (b), indicated by the red cross in (a)], ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=1.65$ [in (c), indicated by the red cross in (a)]. Note that $W_{\pm \infty }$ and the widths $W$ used in (d) are related via Eqs. (24) and (D9) (see Appendix pp4-s1).

Figure 9

Phase diagram showing the closing of the bulk gap at $k_x=0$, i.e., the line $E(k_x=0)=0$ (calculated in the semi-infinite geometry) as a function of the Zeeman field ${\mathrm{\Delta }}_Z$ and the superconducting phase difference $\varphi$ for cases 1 and 2 discussed in Sec. 4. In case 1 the confinement axis of the 2DHG is the crystallographic $\langle 001\rangle$ axis, in case 2 it is the $\langle 110\rangle$ axis. The horizontal lines indicate the critical Zeeman energies ${\mathrm{\Delta }}_{Z,c}$ for each case. The shaded regions indicate the experimentally accessible topological regions. The parameters are given in Appendix pp5.

Figure 10

Setup for case 3 of Sec. 4. The Zeeman energy ${\mathrm{\Delta }}_Z$ is induced by proximitizing the Ge (gray) with a ferromagnetic insulator (FI, red), e.g., EuS [122]. The strength of the induced Zeeman coupling can be adapted by inserting a thin layer of insulating material (green) between the ferromagnetic insulator and the Ge.

Figure 11

(a) Spectrum $E(k_x,k_y=0)$ of the 2DHG in anisotropic Ge [see Eqs. (A1) and (A2)]. (b) Fermi surface. The black solid lines correspond to the two lowest-energy bands of the Luttinger-Kohn (LK) Hamiltonian of Eq. (A1) confined to two dimensions. The red dashed lines are the energy eigenvalues of the effective Hamiltonian $H_{\text{eff,a}}+H_4$ defined in Eqs. (A2) and (7). The blue dotted line corresponds to the effective Hamiltonian expanded up to fourth order in perturbation theory, which corresponds to $H_{\text{eff},a}^{\left(4\right)}$ defined in Eq. (A3). The parameters for the Luttinger-Kohn Hamiltonian are $E_s=0$, ${\gamma }_1=13.38$, ${\gamma }_2=4.24$, and ${\gamma }_3=5.69$ [103]. Due to the normalization, the spectrum is independent of the electric field $\mathcal{E}$. The parameters for the effective model are determined using perturbation theory [103] and they are ${\hslash }^2/2m^{*}=0.906E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\alpha =0.189E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, ${\alpha }_a/\alpha =0.146$, $\mu =0.3E_{\mathcal{E}}$, $\beta =0.06E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, ${\eta }_1=-0.480E_{\mathcal{E}}{l}_{\mathcal{E}}^4$, and ${\eta }_2=0.046E_{\mathcal{E}}{l}_{\mathcal{E}}^4$.

Figure 12

Top row: Spectrum $E(k_x,k_y=0)$, calculated using $H_{\text{eff}}+H_4$ from Eqs. (6) and (7). Bottom row: Fermi surface, i.e., the momenta where $E(k_x,k_y)=0$. The red curves indicate the physical Fermi surfaces, the black curves indicate the unphysical ones. In each column, the strength $\beta$ of the quartic term is different: $\beta =0$ in (a), $\beta =0.10E_{\mathcal{E}}{l}_{\mathcal{E}}^4$ in (b), $\beta =0.12E_{\mathcal{E}}{l}_{\mathcal{E}}^4$ in (c), and $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$ in (d). The other parameters are ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\alpha =0.39E_{\mathcal{E}}{l}_{\mathcal{E}}^3$, $\mu =0.3E_{\mathcal{E}}$. The Fermi momenta on the two physical (red) Fermi surfaces are $k_{F,1}{l}_{\mathcal{E}}\approx 0.44$ and $k_{F,2}{l}_{\mathcal{E}}\approx 0.70$. Therefore, $\beta =0.14E_{\mathcal{E}}{l}_{\mathcal{E}}^4$ gives $\alpha k_{F,2}^3=0.13E_{\mathcal{E}}\approx 4.3\beta k_{F,2}^4$, which satisfies the upper boundary set in Eq. (8).

Figure 13

Energy levels of the Hamiltonian given in Eq. (19) that includes the fictitious quartic term (black solid curve) compared to the energy levels calculated using perturbation theory to third order in the CSOI strength $\alpha$ (red dashed curve). Both calculations are done in the semi-infinite geometry at $k_x=0$. The parameters are $L_y=49.5l_{\mathcal{E}}$, $W=4.5l_{\mathcal{E}}$, ${\hslash }^2/2m^{*}=1.141E_{\mathcal{E}}{l}_{\mathcal{E}}^2$, $\mu =0.3E_{\mathcal{E}}$, and ${\mathrm{\Delta }}_{SC}=\mu /3$. In the top row $\varphi =\pi /4$, in the bottom row $\varphi =3\pi /4$. In column (a) ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=1$, in column (b) ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=2$, and in column (c) ${\mathrm{\Delta }}_Z/{\mathrm{\Delta }}_{\text{SC}}=3$. The values for $\beta$ change together with $\alpha$, as explained in the caption of Fig. 5.

Figure 14

Spectrum $E\left(k_x\right)$ of the 2DHG with parameters used for case 1 in Sec. 4. The momentum $k_y$ is fixed to $k_y=0$ in (a) and $k_y{l}_{\mathcal{E}}=0.18$ in (b). The solid black line corresponds to the two lowest-energy bands of the Luttinger-Kohn (LK) Hamiltonian confined to two dimensions [see Eq. (2)]. The red dashed lines correspond to the energy eigenvalues of the effective model [see Eqs. (6) and (7)]. Due to the fact that the chemical potential is rather small, the two models agree well at the Fermi level. The Luttinger parameters used for the LK Hamiltonian are the same as in Fig. 11, the strain is $E_s/E_{\mathcal{E}}=0.26$, and the chemical potential is $\mu =0.04E_{\mathcal{E}}$. The parameters used for the effective model are given in the column “Case 1” in Table 1.

Figure 15

Phase diagram analogous to Fig. 9, with the only difference that a potential barrier [see Eq. (F1)] is included here. Due to the potential barrier, the topological phase can no longer be reached within the experimental limitations of case 1. Case 2 still has an experimentally accessible topological phase. The barrier height is ${\mu }_{\text{b}}=0.06E_{\mathcal{E}}$ and the barrier width is $W_{\text{b}}=l_{\mathcal{E}}$. All other parameters are as in Fig. 9