Detection of second-order topological superconductors by Josephson junctions

TL;DR

This paper proposes a method to identify second-order topological superconductors using Josephson junctions, revealing a tunable 0-pi transition that indicates topological phase changes and enables electric control of Majorana states.

Contribution

It introduces a novel experimental signature for SOTSs via Josephson junctions and demonstrates electric control of Majorana bound states without magnetic fields.

Findings

Chemical potential tuning induces topological phase transitions.

A stable 0-pi transition is observed in Josephson junctions.

Electric control of Majorana states is achieved without magnetic manipulation.

Abstract

We study Josephson junctions based on second-order topological superconductors (SOTSs) which can be realized in quantum spin Hall insulators with large inverted gap in proximity to unconventional superconductors. We find that tuning the chemical potential in the superconductor strongly modifies the induced pairing of the helical edge states, resulting in topological phase transitions. In a corresponding Josephson junction, a $0$-$\pi$ transition is realized by tuning the chemical potentials in the superconducting leads. This striking feature is stable in junctions with respect to different sizes, doping the normal region, and the presence of disorder. Our transport results can serve as novel experimental signatures of SOTSs. Moreover, the $0$-$\pi$ transition constitutes a fully electric way to create or annihilate Majorana bound states in the junction without any magnetic manipulation.

Equations

$$H(\mathbf{k})$$

$$H_0(\mathbf{k})$$

$$E_{e(h),\uparrow /\downarrow }(k_x)=\pm \text{sgn}(m_y{v}_y)v_x{k}_x-(+)\mu .$$

$${\Psi }_{e,\uparrow ,k_x}(\mathbf{r})=\mathcal{N}{e}^{i{k}_{x}x}(e^{{\lambda }_{1}y}-e^{{\lambda }_{2}y}){\left(\text{sgn}(m_y{v}_y),1\right)}^T.$$

$$h_{\text{BdG}}^x=\text{sgn}(m_y{v}_y)v_x{k}_x{\tau }_z{s}_z-\mu {\tau }_z+{\Delta }_{\text{eff}}^x{\tau }_x{s}_z$$

$${\Delta }_{\text{eff}}^x=-{\Delta }_0+{\Delta }_2\left[m_0/m_y-(1+m_x/m_y){\mu }^2/v_x^2\right].$$

$${\Delta }_{\text{eff}}^y=-{\Delta }_0-{\Delta }_2\left[m_0/m_x-(1+m_y/m_x){\mu }^2/v_y^2\right].$$

$${\mu }_{x(y)}^c=|v_{x(y)}|\sqrt{[m_0/m_{y(x)}\mp {\Delta }_0/{\Delta }_2]/(1+m_{x(y)}/m_{y(x)})}.$$

$$2m_{0}m>v^2.$$

$$E_{\pm }(\varphi )=\pm \frac{{\Delta }_L{\Delta }_R\sin \varphi }{\sqrt{{\Delta }_L^2+{\Delta }_R^2-2{\Delta }_L{\Delta }_R\cos \varphi }}.$$

$${\gamma }_{+}(x)$$

$${\gamma }_{-}(x)$$

$$h_{e,\uparrow }(\mathbf{k})$$

$${\Psi }_{\lambda }(\mathbf{r})$$

$$\det |h_{e,\uparrow }(k_x,-i\lambda )-E-\mu |$$

$${\lambda }_{\alpha }^2$$

$${\psi }_{\alpha \beta }$$

$${\Psi }_{e,\uparrow ,k_x}(E,\mathbf{r})$$

$${\Psi }_{e,\uparrow ,k_x}(y$$

$$\left|\begin{array}{cc}({\psi }_{1+}& {\psi }_{2+})\end{array}\right|$$

$$E+\mu =$$

$${\lambda }_1{\lambda }_2$$

$$k_x^2$$

$${\lambda }_1^2+{\lambda }_2^2$$

$$E_{e,\uparrow }(k_x)$$

$${\Psi }_{e,\uparrow ,k_x}(\mathbf{r})$$

$${\lambda }_{1(2)}$$

$$\frac{1}{2{\mathcal{N}}^2}$$

$${\xi }_{\text{edge}}$$

$$E_{e,\downarrow }(k_x)$$

Full text

Introduction.–The second-order topological superconductor (SOTS) is a novel topological phase of matter and features Majorana zero-dimensional (0D) corner or 1D hinge states which are two spatial dimensions lower than the gapped bulk (Langbehn et al., 2017; Khalaf, 2018; Wang et al., 2018a; Yan et al., 2018; Liu et al., 2018; Geier et al., 2018; Zhu, 2018; Hsu et al., 2018; Wang et al., 2018; Zhang et al., 2019a; Qin et al., ; Bultinck et al., 2019; Peng et al., 2019). They may form stable qubits for topological quantum computation (Kitaev, 2001, 2003; Nayak et al., 2008; Alicea, 2012; Leijnse and Flensberg, 2012; Beenakker, 2013; Elliott and Franz, 2015; Sarma et al., 2015; Sato and Fujimoto, 2016). Recently, the SOTS has been discovered in a variety of realistic materials and triggered tremendous interest (Wang et al., 2018a; Yan et al., 2018; Liu et al., 2018; Geier et al., 2018; Qin et al., ; König and Coleman, ; Zhu, 2018; Hsu et al., 2018; Volpez et al., 2019; Wang et al., 2018; Zhang et al., 2019a; Shapourian et al., 2018; Pan et al., ; Kheirkhah et al., ; Ghorashi et al., ). One way to mimic SOTSs in 2D is given by quantum spin Hall insulators (QSHIs) in proximity to unconventional superconductors with $d_{x^{2}-y^{2}}$-wave or $s_{\pm}$-wave pairing order (Wang et al., 2018a; Yan et al., 2018; Liu et al., 2018). The proximity effect of unconventional superconductivity in 2D systems has been intensively explored in theory (Linder and Sudbø, 2008; Linder et al., 2010; Black-Schaffer and Balatsky, 2013; Zhang et al., 2013; Li et al., 2015; Zareapour et al., 2016; Li et al., 2016; Wu et al., 2016; Wang et al., 2015; Zhou et al., 2019) and experiment (Zareapour et al., 2012; Wang et al., 2013; Zhao et al., 2018; Perconte et al., 2018; Xu et al., 2014; Yilmaz et al., 2014). To date, however, the only way proposed to detect 2D SOTSs is a tunneling experiment without a concrete calculation of the observable signature. An alternative approach to probe SOTSs and manipulate the Majorana corner modes is thus needed. In QSHIs, a finite doping is typically present, and the chemical potential can be far away from the Dirac points. Therefore, it is certainly interesting and experimentally relevant to explore the influence of the chemical potential in SOTSs.

In this Letter, we investigate superconductor-normal metal-superconductor (SNS) junctions formed by a 2D SOTS. The SOTS can be realized in a QSHI with a large inverted gap in proximity to an unconventional superconductor. We introduce a minimal model which is able to capture the essential physics of the SOTS. We find that due to the nontrivial momentum-dependence of the pairing potential and mass, the chemical potential in the SOTS alters the pairing gap opened within the edge states significantly. It can even switch the sign of the pairing gap, leading to a topological phase transition. While the supercurrent across the SNS junction is insensitive to the chemical potential in the N region, it depends strongly on the filling in the superconductors. Strikingly, tuning the chemical potentials in the superconductors gives rise to a [math]-$\pi$ transition, which is absent in junctions based on conventional $s$-wave pairing. These features are robust against disorder in junctions with different sizes. They offer novel signatures to detect the SOTS with Majorana corner states. Furthermore, while Majorana bound states (MBSs) emerge in the [math]-junction when the phase difference across the junction is $\phi=\pm\pi$, they appear at vanishing $\phi$ in the $\pi$-junction. Thus, Josephson junctions with such a doping-induced [math]-$\pi$ transition provide an innovative platform to create or annihilate MBSs by electric gating in the absence of $\phi$. These predictions are applicable to a number of candidate systems including high-temperature QSHIs (Qian et al., 2014; Tang et al., 2017; Fei et al., 2017; Wu et al., 2018; Chen et al., 2018; Weng et al., 2015; Si et al., 2016; Reis et al., 2017; Hsu et al., 2015; Wrasse and Schmidt, 2014; Liu et al., 2015; Wan et al., 2017) in proximity to high-temperature cuprate or iron-based superconductors.

*Minimal model for SOTSs.–*We consider the minimal model for SOTSs realized in QSHIs in proximity to superconductors,

[TABLE]

written in the Nambu basis $(c_{a,\uparrow,{\bf k}},c_{b,\uparrow,{\bf k}},c_{a,\downarrow,{\bf k}},c_{b,\downarrow,{\bf k}},$ $c_{a,\uparrow,-{\bf k}}^{\dagger},c_{b,\uparrow,-{\bf k}}^{\dagger},c_{a,\downarrow,-{\bf k}}^{\dagger},c_{b,\downarrow,-{\bf k}}^{\dagger})$, where $c_{\sigma,s,{\bf k}}^{\dagger}$($c_{\sigma,s,{\bf k}}$) creates(annihilates) an electron with spin $s\in\{\uparrow,\downarrow\}$, orbital $\sigma\in\{a,b\}$ and the momentum ${\bf k}=(k_{x},k_{y})$ measured from the band inversion point of the QSHI. $\bm{\tau}$, $\bm{\sigma}$ and $\bm{s}$ are Pauli matrices acting on Nambu, orbital and spin spaces, respectively. $m({\bf k})=m_{0}-m_{x}k_{x}^{2}-m_{y}k_{y}^{2}$ is the mass term of the QSHI and $\mu$ is the chemical potential. The band inversion implies the conditions $m_{0}m_{x}>0$ and $m_{0}m_{y}>0$ (Bernevig et al., 2006). The pairing potential is written in general as $\Delta({\bf k})=\Delta_{0}+\Delta_{2}\left(k_{x}^{2}-k_{y}^{2}\right)$. When $\Delta_{0}\neq 0$ and $\Delta_{2}=0$, it refers to conventional $s$-wave pairing. When $\Delta_{0}=0$ and $\Delta_{2}\neq 0$, the pairing is formally $d_{x^{2}-y^{2}}$-wave. It can be induced in a QSHI with band inversion at the $\Gamma$ point via the proximity to a cuprate superconductor (Yan et al., 2018). When $0\leqslant|\Delta_{0}|<m_{0}|\Delta_{2}|/2m_{x(y)}$, the system possesses a mixture of $s$-wave and $d_{x^{2}-y^{2}}$-wave pairing. It can also describe effectively a QSHI with band inversion at the $X$ point (Wan et al., 2017; Wrasse and Schmidt, 2014; Liu et al., 2015) and $s_{\pm}$-wave pairing induced from an iron-based superconductor (Stewart, 2011; Hirschfeld et al., 2011; Zhang et al., 2018; Wang et al., 2018b; Zhang et al., 2019b).

In the absence of $\Delta({\bf k})$, the system hosts gapless helical edge states across the bulk gap, which are protected by time-reversal symmetry. The pairing term with $\Delta_{2}\neq 0$ induces a pairing gap of the edge states. The gap may switch sign at the corners, resulting in Majorana corner modes (Yan et al., 2018; Wang et al., 2018a; Liu et al., 2018). We note that although the model (1) is a low-energy effective model, it captures the essential physics of the SOTS. Based on this model, we can understand the second-order topology more intuitively from the picture of edge states and show that it can be strongly altered by changing $\mu$.

Pairing gaps of edge states and topological phase transitions.–To analyze the Majorana corner states and the influence of $\mu$ on the SOTS, we analytically derive the effective model for edge states. For illustration, we consider the edge along $x$ direction of the SOTS in the half-plane $y\leqslant 0$ and assume hard-wall boundary conditions (Note, 5). As in realistic systems, we assume weak pairing. We first calculate the edge states of $H_{0}$, following the approach of Ref. (Zhang et al., 2016). In this model, $k_{x}$ is a good quantum number. The helical electron and hole edge bands are found explicitly as

[TABLE]

The wavefunctions in the orbital basis $\{a,b\}$ read

[TABLE]

They fulfill $\Psi_{e,\downarrow,k_{x}}({\bf r})=\Psi_{e,\uparrow,-k_{x}}^{*}({\bf r})$ and $\Psi_{h,\downarrow/\uparrow,k_{x}}({\bf r})=\Psi_{e,\downarrow/\uparrow,-k_{x}}^{*}({\bf r})$, due to time-reversal and particle-hole symmetries; $\lambda_{1(2)}=|v_{y}/2m_{y}|-(+)(v_{y}^{2}/4m_{y}^{2}-m_{0}/m_{y}+m_{x}k_{x}^{2}/m_{y})^{1/2}$ and $\mathcal{N}$ is the normalization factor. The decaying length of edge states is given by $\xi_{\text{edge}}=\text{max}[1/\text{Re}(\lambda_{1}),1/\text{Re}(\lambda_{2})]$. At zero energy, the electron and hole bands touch at $k_{x}=\pm k_{c}$ with $k_{c}=\text{sgn}(m_{y}v_{y})$ $\mu/v_{x}$. For $\mu=0$, $k_{c}=0$. However, for $\mu\neq 0$, the touching points shift to finite $\pm k_{c}$. Projecting the pairing term onto the edge states, the resulting Bogoliubov de-Gennes (BdG) Hamiltonian for edge states is obtained as

[TABLE]

in the basis $(\Psi_{e,\uparrow},\Psi_{e,\downarrow},\Psi_{h,\downarrow},\Psi_{h,\uparrow})$, and the pairing gap is given by

[TABLE]

Without loss of generality, a real $\Delta({\bf k})$ has been assumed (Note, 1). We provide the derivation in detail in the Supplemental Material (Sup, ). Similarly, for an edge along $y$ direction, we find the BdG Hamiltonian of the same form but with a different pairing gap

[TABLE]

The combination of $\Delta_{\text{eff}}^{x}$ and $\Delta_{\text{eff}}^{y}$ (with opposite signs) in Eq. (4) mimics the Jackiw-Rebbi model (Jackiw and Rebbi, 1976) at corners of $x$ and $y$ axes. Thus, Majorana corner states at zero energy appear if $\Delta_{\text{eff}}^{x}\Delta_{\text{eff}}^{y}<0$.

For $s$-wave pairing, $\Delta_{\text{eff}}^{x}$ and $\Delta_{\text{eff}}^{y}$ are identical and constant. Thus, no corner state emerge. In contrast, for unconventional pairings with $|\Delta_{0}|<m_{0}|\Delta_{2}|/2m_{x(y)}$, we obtain corner states at small $\mu$. When $\mu=0$, the SOTS has two reflection symmetries. When $m_{x}=m_{y}$, $v_{x}=v_{y}$ and $\Delta_{0}=0$, it possesses a fourfold rotation symmetry. In these particular cases, the system can be characterized by a topological invariant calculated from the bulk Hamiltonian (Benalcazar et al., 2017a, b; Song et al., 2017; Sup, ). However, the corner states in our model are not restricted to any crystalline symmetries. Interestingly, $\Delta_{\text{eff}}^{x(y)}$ depends strongly on $\mu$. The $\mu$ dependence stems from the quadratic terms in the model (1), which are crucial for the topological properties of the SOTS. Moreover, $\Delta_{\text{eff}}^{x(y)}$ vanishes at $\mu=\pm\mu_{x(y)}^{c}$, where

[TABLE]

This behavior indicates that we can switch the sign of $\Delta_{\text{eff}}^{x(y)}$ by varying $\mu$. Without loss of generality, we suppose $\mu_{x}^{c}\leqslant\mu_{y}^{c}$. The system is in a SOTS phase in the parameter regions $0\leqslant|\mu|<\mu_{x}^{c}$ and $\mu_{y}^{c}<|\mu|<m_{\text{g}}$ with $m_{\text{g}}$ being the bulk gap Note (2), whereas if $\mu_{x}^{c}<|\mu|<\mu_{y}^{c}$, it becomes a trivial superconductor with no corner state. For the particular case with $v_{x(y)}=v$, $m_{x(y)}=m$, and $\Delta_{0}=0$, $\Delta_{\text{eff}}^{x}$ and $\Delta_{\text{eff}}^{y}$ are always opposite. They both close at $\mu=\pm v\sqrt{m_{0}/2m}$. Thus, there is no parameter space for the trivial phase. Nevertheless, the sign of $\Delta_{\text{eff}}^{x(y)}$ can still be changed by a finite $\mu$ inside the bulk gap (Note, 3) if

[TABLE]

This condition indeed corresponds to a QSHI phase with a large inverted gap or equivalently an indirect bulk gap. It is likely realized in the inverted InAs/GaSb bilayer (Liu et al., 2008; Knez et al., 2011; Krishtopenko and Teppe, 2018), WTe2 monolayer (Qian et al., 2014; Tang et al., 2017; Fei et al., 2017; Wu et al., 2018; Chen et al., 2018), functionalized MXene (Weng et al., 2015; Si et al., 2016), Bismuthene on SiC (Reis et al., 2017; Hsu et al., 2015), and PbS monolayer (Wan et al., 2017; Wrasse and Schmidt, 2014; Liu et al., 2015).

To test our analytical results, we discretize the model (1), put it on a square lattice, choose a proper set of parameters (satisfying the inequality (8)) and calculate the energy spectrum in a ribbon geometry (Note, 6). For concreteness, we consider $\Delta_{0}=0$ and set the lattice constant $a$ to unity. As shown in Fig. 1(a), the edge states for $\mu=0$ open a gap at $k_{x}=0$. As $\mu$ is increased, the gap splits to two points away from $k_{x}=0$. The magnitude of the gap first decreases, vanishes at a critical $\mu$ and then reopens, which explicitly demonstrates a topological phase transition. This behavior is in perfect agreement with Eq. (5), cf. Fig. 1(b).

*0-$\pi$ transition and its robustness.–*We now consider an SNS junction in which two SOTSs (also called S leads below) are connected by a QSHI with length $L$ in $x$ direction, as sketched in Fig. 2(a). The width of the junction ribbon is $W$. For simplicity, we assume the chemical and pairing potentials in step-like forms. $\mu_{L(R)}$ and $\mu_{N}$ denote the chemical potentials in the left(right) S lead and N (QSHI) region, respectively. $\phi$ is the phase difference across the junction. We calculate the supercurrent $J_{s}$ by the lattice Green’s function technique (Asano, 2001; Martín-Rodero et al., 1994; Furusaki, 1994) and provide the details in the Supplemental Material (Sup, ).

At low temperatures, the transport in the junction is conducted dominantly by helical edge channels. Perfect Andreev reflection occurs at the NS interfaces. Thus, the current-phase relation (CPR) takes a sawtooth shape with a sudden jump, see Fig. 2(b). The sawtooth-like CPR is insensitive to $\mu_{N}$ and stays stable in junctions of different sizes ($L$ and $W$), provided that the two edges at $y=\pm W/2$ are well separated, $W\gg\xi_{\text{edge}}$. The sudden jump can be related to the fermion parity anomaly at each edge (Fu and Kane, 2009; Crépin and Trauzettel, 2014). It indicates the formation of degenerate MBSs in the junction discussed below. $J_{s}$ decreases monotonically with increasing $L$, see Fig. 2(b). The critical current $J_{c}$ (maximal value of $J_{s}$) decays as \text{\sim}1/L in long junctions, similar to junctions based on conventional s-wave pairing. In short junctions, $J_{c}$ is of the same order of magnitude but always smaller than $e\sqrt{|\Delta_{L}\Delta_{R}|}/\hbar$, in contrast to the case of $s$-wave pairing. In this estimate, $\Delta_{L(R)}$ is the induced pairing gap of edge states in the left(right) S lead and determined by Eq. (5). We attribute this difference to the inhomogeneity of the superconducting pairing at the boundaries of our setup.

The CPRs for a fixed $\mu_{L}$ and various values of $\mu_{R}$ are displayed in Fig. 2(a). Since $J_{s}$ is even in $\mu_{L(R)}$, we present only the results for $\mu_{L(R)}>0$. While $J_{s}$ is insensitive to $\mu_{N}$, it decreases significantly when we increase $\mu_{R}$. This behavior can be understood as a result of the reduction of $|\Delta_{R}|$ by $\mu_{R}$, see Eq. (5). Strikingly, increasing $\mu_{R}$ further, we observe a clear [math]-$\pi$ transition for the parameters satisfying the inequality (8). While $J_{s}(\phi)$ in the region $0<\phi<\pi$ is positive for $\mu_{R}<\mu^{c}$, it becomes negative for $\mu_{R}>\mu^{c}$. We coin the former case a [math]-junction and the latter one a $\pi$-junction. Meanwhile, the sudden jump of the CPR is switched to $\phi=0$ in the $\pi$-junction, which is in strong contrast to the [math]-junction where the jump is at $\phi=\pm\pi$. In Fig. 2(c), we plot $J_{c}$ as a function of $\mu_{R}$. The critical value $\mu^{c}$ for the transition is approximately given by $v\sqrt{m_{0}/2m}$, in accord with our analytical result. Close to $\mu^{c}$, $J_{c}$ drops quickly and switches sign. These features are generic and apply to junctions of different lengths and widths. They are also robust with respect to nonmagnetic disorder in the N region. To illustrate this, we model the disorder as random on-site potentials in the range $[-V_{\text{dis}}/2,V_{\text{dis}}/2]$ (Li et al., 2018; Sup, ) and calculate 200 random disorder configurations in the inset of Fig. 2(c). There is no qualitative difference in the features compared to those in clean junctions. This can be expected since the helical edge channels which mediate the transport are less sensitive to backscattering. Similar effects can be observed by tuning $\mu_{L}$ and fixing $\mu_{R}$. Finally, it is important to note that the variation of $J_{s}$ and the [math]-$\pi$ transition by tuning $\mu_{L(R)}$ are directly related to the strong $\mu_{L(R)}$-dependence in $\Delta_{L(R)}$ in the SOTS, and absent in conventional junctions based on $s$-wave pairing.

*Majorana bound states.–*Next, we discuss the Andreev bound states (ABSs) formed in the junction, which can be obtained from the lattice Green’s function. In short junctions, there are two bands of ABSs with opposite energies, see Fig. 3(a,c). When the sudden jump of the CPR occurs, the positive and negative bands touch at zero energy. This degeneracy is robust and protected by time-reversal and particle-hole symmetries. It resembles Kramers pairs of MBSs. This can be best understood from the effective Hamiltonian (4) for edge states. In the short junction limit, two ABS bands at a given edge can be described by

[TABLE]

Notably, the ABSs are confined in the pairing gaps for $\phi$ satisfying (\cos\phi$$-$$\Delta_{L}/\Delta_{R})(\cos\phi$$-$$\Delta_{R}/\Delta_{L})>0, as verified in Fig. 3(a,c). Noticing $\Delta_{L}\Delta_{R}>0$ in the [math]-junction for $\mu_{R}<\mu^{c}$, whereas $\Delta_{L}\Delta_{R}<0$ in the $\pi$-junction for $\mu_{R}>\mu^{c}$, we can see that $E_{\pm}(\phi)$ touch at $\phi=\pm\pi$ and [math], respectively. Using the valid formula at zero temperature, $J_{s}(\phi)=\partial|E_{+}(\phi)|/\partial\phi$ (Beenakker, 1991), we also reproduce the sudden jump in the CPR. The wavefunctions of the zero modes can be written as

[TABLE]

where $\Psi_{\pm}=\text{sgn}(\Delta_{L})\Psi_{h/e,\downarrow}\mp i\Psi_{e/h,\uparrow}$ and the spatial dependence is \eta(x)=\exp\{[\theta(x)(i\mu_{R}$$-$$|\Delta_{R}|)+\theta(-x)(i\mu_{L}+|\Delta_{L}|)]x/v\} (Sup, ). Since $\Psi_{h,\downarrow/\uparrow}=\Psi_{e,\downarrow/\uparrow}^{*}$, the zero modes have self-adjoint wavefunctions, $\gamma_{\pm}(x)=\gamma_{\pm}^{*}(x)$. They are Majorana fermions. Under the time-reversal operation $\mathcal{T}$, $\mathcal{T}\Psi_{e,\downarrow/\uparrow}=\pm\Psi_{e\uparrow/\downarrow}^{*}$ and $\mathcal{T}\Psi_{h,\downarrow/\uparrow}=\pm\Psi_{h,\uparrow/\downarrow}^{*}$. Therefore, $\gamma_{\pm}$ are related by time-reversal symmetry, $\mathcal{T}\gamma_{+}=\gamma_{-}$. A similar analysis can be applied to the other edge where another Kramers pair of MBSs are located. In long junctions, all the features persist but with more pairs of discrete ABS bands emerging from the continuum spectrum, see Fig. 3(b,d).

At $\phi=0$, the MBSs emerge for $\mu>\mu^{c}$, whereas they disappear for $\mu<\mu^{c}$. In this sense, we are able to switch between the presence and absence of MBSs by gating the S leads in the absence of $\phi$. Our setup indeed realizes fully electrically controllable MBSs without fine tuning of magnetic field or threaded flux. This is an important advantage compared to previous proposals based on conventional $s$-wave superconductivity (Fu and Kane, 2008, 2009; Lutchyn et al., 2010; Sau et al., 2010; Oreg et al., 2010; Volpez et al., 2019). Moreover, since the localization lengths of the MBSs in the S leads are determined by $\xi_{L(R)}=|v/\Delta_{L(R)}|$, we are also able to control the spatial profiles of the MBSs by $\mu_{L(R)}$.

Experimental relevance and summary.–Now we briefly discuss the experimental relevance of our proposal. QSHIs with large inverted gaps (Qian et al., 2014; Tang et al., 2017; Fei et al., 2017; Wu et al., 2018; Chen et al., 2018; Weng et al., 2015; Si et al., 2016; Reis et al., 2017; Hsu et al., 2015; Wrasse and Schmidt, 2014; Liu et al., 2015; Wan et al., 2017; Liu et al., 2008; Knez et al., 2011; Krishtopenko and Teppe, 2018) in proximity to cuprate or iron-based superconductors (Stewart, 2011; Hirschfeld et al., 2011; Zhang et al., 2018; Wang et al., 2018b; Zhao et al., 2018; Zareapour et al., 2012; Wang et al., 2013; Zhang et al., 2019b) could provide promising platforms to verify our predictions. For concreteness, we take the inverted InAs/GaSb bilayer and WTe2 monolayer to estimate $\mu_{x(y)}^{c}$. For simplicity, we consider $\Delta_{0}=0$ such that $\mu^{c}_{x}(y)$ is independent of the magnitude of the pairing potential. For the inverted InAs/GaSb bilayer, $m_{0}=0.0055$ eV, $m_{x(y)}=81.9$ eV$\cdot\mathring{\text{A}}$2, $v_{x(y)}=0.72$ eV$\cdot\mathring{\text{A}}$ (Liu and Zhang, 2013). To realize the [math]-$\pi$ transition, one can fabricate the Josephson junction in any direction and find that $\mu_{x(y)}^{c}=0.0042$ eV which is smaller than the bulk gap $m_{\text{gap}}=0.005$ eV. For the WTe2 monolayer with $m_{0}=0.33$ eV, $m_{x}=4.6$ eV$\cdot\mathring{\text{A}}{}^{2}$, $m_{y}=16.9$ eV$\cdot\mathring{\text{A}}$2, $v_{x}=2.55$ eV$\cdot\mathring{\text{A}}$ and $v_{y}=0.3$ eV$\cdot\mathring{\text{A}}$ (Qian et al., 2014), we have $\mu_{x}^{c}=0.252$ eV, $\mu_{y}^{c}=0.057$ eV and $m_{\text{gap}}=0.08$ eV. Thus, it is better to design the junction in $y$ direction in our model (Note, 7). According to Eq. (7), the inclusion of a small $\Delta_{0}$ would suppress $\mu_{x}^{c}$ or $\mu_{y}^{c}$ and hence make it more feasible to observe the [math]-$\pi$ transition. A particle-hole symmetry breaking term, which is neglected here, breaks the symmetry with respect to $\mu$ but does not qualitatively change our main results.

We note in passing that there have been experimental efforts trying to incorporate unconventional superconductivity in topological systems (Zareapour et al., 2012; Wang et al., 2013; Zhao et al., 2018; Zhang et al., 2018; Wang et al., 2018b; Zhang et al., 2019b). Moreover, large proximity-induced pairing gaps in 2D systems from unconventional superconductors have been probed (Zareapour et al., 2012; Wang et al., 2013; Zhao et al., 2018; Perconte et al., 2018).

In summary, we have found that the chemical potentials in superconductors can be used to modulate the supercurrent and realize a [math]-$\pi$ transition in Josephson junctions based on SOTSs. These features are attributed to the dependence of the pairing gap of edge states on the chemical potential. They could serve as novel experimental signatures of the SOTS. We have predicted the [math]-$\pi$ transition as a fully electric way to create or annihilate MBSs at elevated temperatures.

Appendix A Derivation of the effective BdG model for edge states

A.1 Edges in $x$ or $y$ direction

In the absence of the pairing potential, the Bogoliubov-de Gennes (BdG) Hamiltonian [Eq. (1) in the main text] decouples into four blocks. Each block can be analyzed separately. Let us take the block for spin-up electrons for illustration, following the approach of Ref. (Zhang et al., 2016). The block for spin-up electrons reads

[TABLE]

in the basis $(c_{a,\uparrow},c_{b,\uparrow})$, where $m({\bf k})=m_{0}-m_{x}k_{x}^{2}-m_{y}k_{y}^{2}$.

Consider the edge in $x$ direction of a semi-infinite SOTS in the half-plane $y\leqslant 0$ and impose hard-wall boundary conditions. Then, $k_{x}$ is a good quantum number. We assume the trial wavefunction of the form,

[TABLE]

where $\psi_{\lambda}$ is a two-component spinor. Plugging Eq. (12) in $h_{e,\uparrow}(-i\nabla_{{\bf r}})\Psi_{\lambda}=E\Psi_{\lambda}$, we obtain the secular equation

[TABLE]

for a nontrivial solution of $\psi_{\lambda}$. Solving Eq. (13) gives four $\lambda$, denoted as $\beta\lambda_{\alpha}$ with $\beta=\pm$ and $\alpha=1,2$,

[TABLE]

where $m_{k}=(m_{0}-m_{x}k_{x}^{2})/m_{y}.$ Each $\beta\lambda_{\alpha}$ corresponds to a spinor state written as

[TABLE]

or alternatively,

[TABLE]

Then, a general wavefunction is given by

[TABLE]

where the energy $E$ and coefficients $C_{\alpha\beta}$ are found from the boundary conditions. The hard-wall boundary conditions read

[TABLE]

The condition $\Psi_{e,\uparrow,k_{x}}(y=-\infty)=0$ requires that $\Psi_{e,\uparrow,k_{x}}$ contains only the terms with $\beta=+$ and $\mathrm{Re}(\lambda_{\alpha})>0$, i.e, $C_{1-}=C_{2,-}=0$. The other condition $\Psi_{e,\uparrow,k_{x}}(y=0)=0$ then leads to

[TABLE]

Plugging Eqs. (15) and (16) into Eq. (19), respectively, and considering $\lambda_{1}\neq\lambda_{2}$, we obtain

[TABLE]

By comparing these two equations, we identify

[TABLE]

According to Eq. (14), there are two cases of $\lambda_{1,2}$, one is $\lambda_{1(2)}>0$ and the other $\lambda_{1}=\lambda_{2}^{*}$. In both cases, $\lambda_{1}\lambda_{2}>0.$ This determines the region for well-localized edge states:

[TABLE]

From Eq. (14), we derive

[TABLE]

By exploiting Eqs. (21) and (23) and $\lambda_{1}+\lambda_{2}>0$, we obtain $\lambda_{1}+\lambda_{2}=|v_{y}/m_{y}|$. With this result in Eq. (20), we find the dispersion as

[TABLE]

and consequently the wavefunction as

[TABLE]

where the two penetration depths $\lambda_{1,2}$ and the normalization factor are given, respectively, by

[TABLE]

The decaying length of the edge states is determined by

[TABLE]

Note that in contrast to previous studies (Yan et al., 2018; Liu et al., 2018), we neither neglect nor treat the quadratic terms as perturbations.

Similarly, the edge states for the other three blocks are found as

[TABLE]

Their wavefunctions in the orbital basis $\{a,b\}$ can be related to Eq. (25) by exploiting time-reversal and particle-hole symmetries, i.e.,

[TABLE]

Next, we calculate the pairing gap $\Delta_{\text{eff}}^{x}$ of edge states. At the Fermi energy ($E=0$), $E_{e,\uparrow}$ crosses $E_{h,\downarrow}$ at $k_{c}$, while $E_{e,\uparrow}$ crosses $E_{h,\downarrow}$ at $-k_{c}$, where

[TABLE]

At the crossing point $k_{c}$, $\Delta_{\text{eff}}^{x}$ is given by

[TABLE]

Using Eqs. (25) and (29), it is found explicitly as

[TABLE]

Similarly, the pairing gap between $E_{e,-}$ and $E_{h,-}$ at $-k_{c}$ is found as $-\Delta_{\text{eff}}^{x}$. Therefore, the full BdG Hamiltonian for the edge states in $x$ direction can be written as

[TABLE]

in the basis $(\Psi_{e,\uparrow},\Psi_{e,\downarrow},\Psi_{h,\downarrow},\Psi_{h,\uparrow})$, where $\bm{\tau}$ and ${\bf s}$ are Pauli matrices acting on Nambu and spin spaces, respectively. Notably, this BdG Hamiltonian is only effective for the excitation near the crossing points $k_{x}=\pm k_{c}$.

A.2 Edges in an arbitrary direction

In this subsection, we will show that the corner states are more generic and not restricted to a specific choice of directions (i.e., $x$ or $y$ direction) of the edges. To this end, we consider the edge in an arbitrary direction $x_{1}$ which has the angle $\theta$ relative to $x$ direction. For simplicity, we consider the isotropic QSHI case, $m_{x}=m_{y}=m$ and $v_{x}=v_{y}=v$. We note that the main conclusion persists in the general anisotropic case. To derive the edge states, it is convenient to use the $x_{1}$ and $x_{2}$ (normal to $x_{1}$) coordinates. The $x_{1}\text{-}x_{2}$ and $x\text{-}y$ coordinates are related by the rotations

[TABLE]

In the $x_{1}\text{-}x_{2}$ coordinates, the BdG Hamiltonian becomes

[TABLE]

in the new basis $(c_{a,\uparrow^{\prime}},c_{b,\uparrow^{\prime}},c_{a,\downarrow^{\prime}},c_{b,\downarrow^{\prime}},c_{a,\uparrow^{\prime}}^{\dagger},c_{b,\uparrow^{\prime}}^{\dagger},c_{a,\downarrow^{\prime}}^{\dagger},c_{b,\downarrow^{\prime}}^{\dagger})$, where the subscript $\prime$ implies that the direction of spin polarization is also rotated, and

[TABLE]

The full BdG Hamiltonian (35) always decouples into two blocks, $H_{0}$ and its time-reversal counterpart, similar to that before the rotation.

In the absence of the pairing potential, each BdG block further decouples into two sub-blocks, one for electrons and one for holes. The sub-blocks take exactly the same form as those without rotation. Following the same approach, the dispersion of spin-up electrons and spin-down holes are given, respectively, by

[TABLE]

Accordingly, the wavefunctions are

[TABLE]

in the orbital basis $(a,b)$, where

[TABLE]

With the wavefunctions in Eqs. (38), the pairing interaction between the edge electrons and holes is calculated as

[TABLE]

where

[TABLE]

Using the expressions (39) of $\lambda_{1}$ and $\lambda_{2}$, we derive

[TABLE]

Therefore, the pairing interaction is given by

[TABLE]

According to Eqs. (37), the crossing point between the electron and hole bands is

[TABLE]

Thus, the pairing gap at $q_{1}=q_{c}$ is

[TABLE]

Here, the phase factor $e^{i\theta}$ stems from the rotation of spin, while the $\theta$ dependence in the brackets comes from the rotation of coordinates. Note that in this derivation, the SOTS is on the right hand side while the vacuum is on the left hand side. In the spin basis for $\theta=0$, the phase factor $e^{i\theta}$ is discarded. Thus, in this basis, the pairing gap reads

[TABLE]

When $\theta=0$ and $\pi/2$, we recover the results for the edge in $x$ and $y$ directions, respectively:

[TABLE]

To form corner states, we need another edge. Let us consider the other edge in $x_{2}$ direction and the SOTS in the $x_{1}<0$ half plane. Along the lines of that we did for the $x_{1}$ edge, we can find analytically the electron and hole edge bands as

[TABLE]

and their wavefunctions as

[TABLE]

where $\lambda_{1(2)}$ and $\mathcal{N}$ are given by Eqs. (39) and (40), respectively. The pairing interaction between the electron and hole bands is found as

[TABLE]

At the crossing point $q_{2}=q_{c}=\text{sgn}(mv)\mu/v$ and in the spin basis for $\theta=0$, the pairing gap is given by

[TABLE]

When $\theta=-\pi/2$ and [math], we recover again the results for the edges in $y$ and $x$ direction, respectively.

Denote the angle between $x_{1}$ and $x$ direction by $\theta_{1}$, and the angle between $x_{2}$ and $y$ direction by $\theta_{2}$. The pairing gap of the edge states in $x_{1}$ and $x_{2}$ directions are

[TABLE]

Note that the spin basis is the same for all edge states (i.e., the spin basis in the particular $x$-$y$ coordinates). The existence of corner states yields that $\widetilde{\Delta}_{\text{gap}}^{X_{1}}(\theta_{1})$ and $\widetilde{\Delta}_{\text{gap}}^{X_{2}}(\theta_{2})$ have different signs,

[TABLE]

For $\Delta_{0}=0$ and considering, in general, $m_{0}/m-2\mu^{2}/v^{2}\neq 0$, Eq. (54) simplifies to

[TABLE]

The phase diagram for corner states is displayed in Fig. 4. One can see that the corner states exist in a wide range of the angles $\theta_{1}$ and $\theta_{2}$ (see the blue areas). This indicates that the corner states, in general, do not require a crystalline symmetry.

Appendix B Calculations of Josephson current

There are different tight-binding lattice models having the low-energy minimal Hamiltonian we consider. In the customary regularization, we can obtain a tight-binding model by replacing $k_{x(y)}\rightarrow\sin k_{x(y)}$ and $k_{x(y)}^{2}\rightarrow 2[1-\cos k_{x(y)}]$. For convenience, the lattice constant is set to unity. Then, Fourier transforming into lattice space, the BdG Hamiltonian is given by

[TABLE]

with

[TABLE]

where the spinor operators are $C_{l,j}^{\dagger}=(c_{a,\uparrow;l,j}^{\dagger},c_{b,\uparrow;l,j}^{\dagger},$ $c_{a,\downarrow;l,j}^{\dagger},c_{b,\downarrow;l,j}^{\dagger},c_{a,\uparrow;l,j},c_{b,\uparrow;l,j},c_{a,\downarrow;l,j},c_{b,\downarrow;l,j})$; $M_{0}=m_{0}-2m_{x}-2m_{y}$; $\{l,l^{\prime}\}$ and $\{j,j^{\prime}\}$ denote the lattice sites in $x$ and $y$ directions, respectively. The identity matrices for spin, Nambu and orbital spaces are omitted for ease of notation.

For an SNS junction, the BdG Hamiltonian can be written as

[TABLE]

where

[TABLE]

and $\Theta(x)$ is the Heaviside step function; $L$ and $W$ (in units of the lattice constant $a$) are the length and width of the junction, respectively. The chemical and pairing potentials are modeled as

[TABLE]

The dc Josephson current can be calculated as (Furusaki, 1994; Martín-Rodero et al., 1994; Asano, 2001)

[TABLE]

where the over-script $\check{...}$ indicates 8W\text{\text{\times}}8W matrices expanded in Nambu, spin, orbital, and lattice ($j=1,...,W$) spaces,

[TABLE]

$\check{\mathcal{G}}(l,l^{\prime},i\omega_{\nu})$ is the Matsubara Green’s function, $\omega_{\nu}=(2\nu+1)\pi k_{B}T$ with $\nu=0,\pm 1,...$, is the Matsubara frequency and $T$ is the temperature. $\mathbb{I}_{W\times W}$ denotes the W\text{\text{\times}}W identity matrix. We find $\check{\mathcal{G}}(l,l^{\prime},i\omega_{\nu})$ numerically by the recursive Green’s function technique (Lee et al., 1981). The trace is taken over Nambu, spin, orbital, and lattice degrees of freedom. The supercurrent $J_{s}$ is independent of $l$ (Asano, 2001). Thus, it is convenient to calculate $J_{s}$ at $l=L$.

By performing the analytical continuation $i\omega_{\nu}\rightarrow E+i\Gamma$ with a positive infinitesimal $\Gamma$, we obtain the retarded Green’s function $\check{G}^{R}(l,l,E+i\Gamma)$. The density of states is then calculated as

[TABLE]

where $1\leqslant l\leqslant L$. The energy of Andreev bound states (ABSs) can be found as the peaks of $\rho(E)$. Note that $\rho(E)$ is the same for $l\in\{1,...,L\}$. A small but finite bandwidth $\Gamma$ is employed for the calculation of $\rho(E)$. In this work, we use $\Gamma=10^{-5}m_{0}$ throughout.

To show the $L$-dependence in the supercurrent $J_{s}(\phi)$, we plot in Fig. 5 the current-phase relation for different values of $L$, and $J_{c}$ and $J_{c}L$ as functions of $L$ in the insets (I) and (II), respectively. We can observe that $J_{s}(\phi)$ decays monotonically as we increase $L$. For long junctions $L\gg\text{max}\{v_{x}/|\Delta_{L}|,v_{x}/|\Delta_{R}|\}$, $J_{c}L$ saturates to a constant. This indicates that $J_{c}$ scales as $\sim 1/L$. Moreover, $J_{s}(\phi)$ is insensitive to the width $W$ as long as $W\gg\xi_{\text{edge}}$.

Appendix C Supercurrent from the edge BdG model

In this section, we look at the edge BdG Hamiltonian (33) and derive analytically the ABSs and Majorana bound states (MBSs). Without loss of generality, we assume $m_{y}v_{y}>0$. Then, one edge of the Josephson junction, say the upper one, is described by

[TABLE]

where the spatially dependent chemical and pairing potentials are

[TABLE]

In the N region, the basis functions can be written as

[TABLE]

where $k_{e(h)}=(\mu_{N}\pm E)/v_{x}.$ Thus, the wavefunction in the N region is expanded as

[TABLE]

In the S leads, the basis functions are

[TABLE]

where $k_{qe(qh)}^{s}=\mu_{s}\pm\Omega_{s}$ and

[TABLE]

with $s\in\{L,R\}$ distinguishing the left and right S leads. $\phi_{L}=0$ and $\phi_{R}=\phi$. We are most interested in the ABSs whose energies satisfy $|E|<\Delta_{s}$. Thus, the wavefunction in the S leads is given by

[TABLE]

The energies $E$ of ABSs and the coefficients $A_{e(h),\pm}$, $B_{e(h)L}$ and $B_{e(h)R}$ are found from the continuity of the wavefunction, i.e.,

[TABLE]

A nontrivial solution of these equations yields

[TABLE]

This can be recast to the transcendental equations

[TABLE]

with $\phi$ satisfying

[TABLE]

The solutions of $E$ can be found self-consistently from Eq. (72). With the obtained $E$, the coefficients are also obtained. Several salient features of ABSs are obvious: (i) ABSs appear in pairs with opposite energies; (ii) the ABS spectrum is independent $\mu_{N}$; (iii) more ABS branches appear for a longer $L$.

In the short junction limit $L=0$, the ABS spectrum can be found analytically as

[TABLE]

Correspondingly, equation (73) defines the parameter range for the existence of ABSs

[TABLE]

Note that for the 0-junction, $\Delta_{L}\Delta_{R}>0$, while for the $\pi$-junction, $\Delta_{L}\Delta_{R}<0$. From Eqs. (74) and (75), it is easy to see that the zero-energy modes are at $\phi=\pm\pi$ in the 0-junction, while they switched to be at $\phi=0$ in the $\pi$-junction. Note that this result holds also in longer junctions. In both junctions, the wavefunctions of two zero-energy modes can be written in a compact form

[TABLE]

where

[TABLE]

Restoring the basis $(\Psi_{e,\uparrow},\Psi_{e,\downarrow},\Psi_{h,\downarrow},\Psi_{h,\uparrow}),$ we can write

[TABLE]

Recall Eqs. (29), $\Psi_{h,\uparrow/\downarrow}=\Psi_{e,\uparrow/\downarrow}^{*}.$ Hence, the two zero-energy modes obey

[TABLE]

This indicates that they are related by particle-hole symmetry. We can recombine them and obtain

[TABLE]

The new zero-energy modes have self-adjoint wavefunctions

[TABLE]

and behave like MBSs. Under time-reversal operation $\mathcal{T}$,

[TABLE]

This shows that the two MBSs are connected by time-reversal symmetry,

[TABLE]

Hence, they are Kramers partners.

Appendix D Symmetries and quadrupole moment

In this section, we analyze the symmetries and calculate the quadrupole moment of the SOTS.

D.1 Symmetries

The BdG Hamiltonian Eq. (1) in the main text

[TABLE]

satisfies the following symmetries:

$\bullet$ time-reversal symmetry, $\mathcal{T}=is_{y}\mathcal{K}$ with $\mathcal{K}$ the complex conjugation:

[TABLE]

$\bullet$ particle-hole symmetry, $\Xi=\tau_{x}\mathcal{K}$:

[TABLE]

$\bullet$ inversion symmetry, $\mathcal{P}=\sigma_{z}$:

[TABLE]

$\bullet$ if $\mu=0$, combined reflection symmetries, $\mathcal{M}_{x}=\tau_{x}s_{x}\sigma_{x}$ and $\mathcal{\mathcal{M}}_{y}=\tau_{x}s_{x}\sigma_{y}$:

[TABLE]

$\bullet$ if $m_{x}=m_{y}$, $v_{x}=v_{y}=v$ and $\Delta_{0}=0$, combined four-fold rotation symmetry, $\mathcal{C}_{4}=\tau_{z}e^{i\pi\sigma_{z}/4}$:

[TABLE]

D.2 Quadrupole moment

To find the quadrupole moment by the Wilson-loop approach of Benalcazar et al. (2017a, b), we need to consider a periodic lattice model. As the previous numerical calculations, we consider the lattice model by replacing $k_{x(y)}\rightarrow\sin k_{x(y)}$ and $k_{x(y)}^{2}\rightarrow 2[1-\cos k_{x(y)}]$ in Eq. (84)

The projected position operator $P^{\text{occ}}\hat{x}P^{\text{occ}}$ into the occupied bands can define a Wilson line. In $x$ direction, the Wilson line operator is given by

[TABLE]

where $[F_{k}]^{nn^{\prime}}=\langle u_{k+\delta_{k_{x}}}^{n}|u_{k}^{n^{\prime}}\rangle$, $\delta_{k_{x}}=2\pi/N_{x}$ and $|u_{k}^{n}\rangle$ are the eigenstates of the lattice Hamiltonian. $N_{x}$ is the number of lattice sites in $x$ direction. For the limit $N_{x}\rightarrow\infty$, $[F_{k}]^{nn^{\prime}}\approx e^{-iA_{nn^{\prime}}\delta_{k}}$ with $A_{nn^{\prime}}=i\langle u_{k}^{n}|\partial_{k_{x}}u_{k}^{n^{\prime}}\rangle$ the non-Abelian gauge field.

A Wilson loop $\mathcal{W}_{k_{x}+2\pi\leftarrow k_{x}}$ is defined as a Wilson line that goes across the entire Brillouin zone. It is unitary and its eigenvalues take the form

[TABLE]

where $R\in\{0,1,...,N-1\}$. The phases $\theta_{x}^{j}$ of the eigenvalues are called the Wannier centers. They correspond to the position of the electrons relative to the center of the unit cell (Yu et al., 2010). The Wilson loop can be connected to the Wannier Hamiltonian $H_{\mathcal{W}}(k)$ of the edge, $\mathcal{W}_{x,k}\equiv e^{iH_{\mathcal{W}}(k)}.$ It can be adiabatically related to the physical Hamiltonian of the $x$ edge (Fidkowski et al., 2010). Correspondingly, $\theta_{x}^{j}(k_{y})$ with $j\in\{1,...,N_{\text{occ}}\}$ are refereed to the Wannier spectrum (or bands). It depends on the $k_{y}$ coordinate. Given the normalized Wilson-loop eigenstate, the eigenstates of $\mathcal{W}_{k_{x}+2\pi\leftarrow k_{x}}$ are written as

[TABLE]

where $[v_{k}^{j}]^{n}$ is the $n$th component of the $j$th Wilson-loop eigenstate $|v_{k}^{j}\rangle$. Note that while the Wilson-loop eigenvalues $\theta_{x}^{j}$ do not depend on the base point $k_{x}$, their eigenstates $|v^{j}\rangle$ do. The electronic contribution to the dipole moment, called polarization, is proportional to

[TABLE]

Consider the SOTS with $\mu=0$ and vary $\Delta_{0}$. There are $N_{\text{occ}}=2$ Wannier bands $\theta_{x}^{\pm}(k_{y})$ corresponding to the two occupied bands, as shown in Fig. 6. These two Wannier bands, in general, do not touch at any point over $k_{y}$ except for $\Delta_{0}=2\Delta_{2}m_{0}/m$. They obey $\theta_{x}^{+}(k_{y})+\theta_{x}^{-}(k_{y})=0$ mod 1. Thus, the total polarization is always zero. We can define the two Wannier bands as $\theta_{x}^{+}\in[0,1/2)$ and $\theta_{x}^{-}\in(1/2,1]$.

Following a similar approach as that for the lattice model, we can define a nested Wilson loop $\mathcal{\widetilde{W}}_{k_{y}+2\pi\leftarrow k_{y}}$ for the Wannier bands with the Wannier functions given by Eq. (92) and calculate the associated polarization $p_{y}^{\theta_{x}^{\pm}}$. Under reflections $\mathcal{M}_{x}$, $\mathcal{M}_{y}$, and inversion $\mathcal{I}$, the Wannier sector polarization obey

[TABLE]

mod 1. Thus, $p_{y}^{\theta_{x}^{\pm}}$ must quantize (at 0 or $1/2$) in the presence of the symmetries $\mathcal{M}_{x}$, $\mathcal{M}_{y}$ and $\mathcal{I}$. The relations and result for the other Wannier polarization $p_{x}^{\theta_{y}^{\pm}}$ are the same as above but with exchanging $x\leftrightarrow y$. In reflection symmetric insulators, the Wannier polarization can be alternatively computed from the eigenvalues of symmetry operators at the reflection-invariant momenta (Benalcazar et al., 2017b). The existence of corner states can be associated with the quantized polarization $p_{y}^{\theta_{x}^{\pm}},p_{x}^{\theta_{y}^{\pm}}=1/2$.

The quadrupole moment $Q_{xy}$ can be written as

[TABLE]

The SOTS at $\mu=0$ preserves $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$ symmetries. Thus, $Q_{xy}$ is always quantized and can be identified as the topological invariant for the SOTS. When increasing $\Delta_{0}$, the quadrupole moment changes from 1/2 to 0 at $\Delta_{0}=2\Delta_{2}m_{0}/m$ where the Wannier gap closes, as shown in Fig. 7. However, when $\mu\neq 0$, both $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$ symmetries are broken. Then, $Q_{xy}$ is no longer quantized. Nevertheless, since we can smoothly vary $\mu$ to the particular limit (with $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$ symmetries) without closing either the bulk or edge gap, the general SOTS phase is topologically equivalent to the SOTSs that preserve these reflection symmetries.

Appendix E Results from another typical lattice model

In this section, we will show that the minimal model can be derived from different lattice models for SOTSs at low energies. For instance, we consider another typical lattice model for SOTSs given by (Wang et al., 2018a)

[TABLE]

in the basis $(c_{a,\uparrow},c_{a,\downarrow},c_{b,\uparrow},c_{b,\downarrow},c_{a,\uparrow}^{\dagger},c_{a,\downarrow}^{\dagger},c_{b,\uparrow}^{\dagger},c_{b,\downarrow}^{\dagger})$, where $M({\bf k})=2t(\cos k_{x}-\cos k_{y})+4t_{1}\cos k_{x}\cos k_{y}$ and $\Delta({\bf k})=\Delta_{0}+2\Delta_{2}(\cos k_{x}+\cos k_{y})$ with $|\Delta_{0}|<2|\Delta_{2}|$. The Pauli matrices $\bm{\tau}$, $\bm{\sigma}$ and ${\bf s}$ act on Nambu, orbital and spin spaces, respectively. This lattice model describes a QSHI with $s_{\pm}$-wave pairing potential but with band inversion at the $X=(\pi,0)$ point.

Around the $X$ point, we expand all terms up to quadratic order in momentum ${\bf k}$:

[TABLE]

Here, ${\bf k}$ is measured from the $X$ point. Rearranging the basis to $(c_{a,\uparrow},-c_{b,\downarrow},c_{a,\downarrow},c_{b,\uparrow},c_{a,\downarrow}^{\dagger},c_{b,\uparrow}^{\dagger},-c_{a,\uparrow}^{\dagger},c_{b,\downarrow}^{\dagger})$, we then obtain the low-energy effective model as

[TABLE]

where $k_{\pm}=k_{y}\pm ik_{x}$. It takes the some form of the minimal model (1) in the main text.

Taking a set of parameters that satisfy the large inverted gap condition, the energy spectra of the lattice model (96) of a ribbon geometry in $x$ direction are presented in Fig. 8. At $k_{x}=\pm\pi$, band inversion happens and edge states appear nearby. More importantly, with increasing the chemical potential from zero, we can again observe a gap closing and reopening at the edge states without closing the bulk gap.