Majorana bound states and subgap states in three-terminal topological superconducting nanowire-quantum dot hybrid devices
Guang-Yao Huang Xin Liu H. Q. Xu
cond-mat.mes-hall






Beijing Key Laboratory of Quantum Devices, Key Laboratory for the
Physics and Chemistry of Nanodevices, and Department of Electronics,
Peking University, Beijing 100871, China


School of Physics and Wuhan National High Magnetic Field Center,
Huazhong University of Science and Technology, Wuhan, Hubei 430074,
China


[Corresponding author: ]hqxu@pku.edu.cn; hongqi.xu@ftf.lth.se
Beijing Key Laboratory of Quantum Devices, Key Laboratory for the
Physics and Chemistry of Nanodevices, and Department of Electronics,
Peking University, Beijing 100871, China
Division of Solid State Physics, Lund University, Box 118, S-221
00 Lund, Sweden




Three-terminal topological superconducting nanowire (TSNW)-quantum dot (QD) hybrid junction devices are studied. The energy spectra and the wave functions  of the subgap states
are calculated as a function of the superconducting phase
differences between TSNWs and as a function of the QD level energy based on the Bogoliubov-de Gennes tight-binding Hamiltonians. It is shown that when the QD level is located near or inside the superconducting gap, there can exist eight subgap states. Among them, four low energy (two positive and two negative) subgap states are essentially formed by linear combinations of the six Majorana bound states (MBSs) located at the ends of the three TSNWs. The remaining four high energy subgap states are mainly built from linear combinations of the QD state and the three MBSs of the TSNWs adjacent to the QD. When there is no QD level near or inside the superconducting gap, only six subgap states built from linear combinations of the six MBSs of the three TSNWs can be present in the system. It is also shown that there exists a unique point in the parameter space of the superconducting phase differences between TSNWs, at which the energies of the six low energy subgap states move close to each other towards nearly zero energies. Simple but general effective model Hamiltonians for the three-terminal TSNW-QD hybrid devices have also been developed. Based on the effective model Hamiltonians, the subgap states of the three-terminal TSNW-QD hybrid devices in the limit of the three infinitely long TSNWs are studied.  The results of the calculations and the effective model Hamiltonians could be used as a starting point to construct and  investigate the braiding schemes of MBSs in TSNW junction devices.


Majorana bound states and subgap states in three-terminal topological superconducting nanowire-quantum dot hybrid devices
    H. Q. Xu
    December 30, 2023
=========================================================================================================================





§ INTRODUCTION


Currently, Majorana bound states (MBSs) in condensed matter physics <cit.>
attract great attention because of its non-Abelian statistics and potential applications in topological quantum computation
(TQC) <cit.>. A series
of experiments towards detecting signatures of the MBSs have been performed on nanowire devices <cit.> made from semiconductors InSb and InAs nanowires, and ferromagnetic atomic chains. These nanowire systems are believed to be the most promising systems for generating MBSs. To realize MBSs with a semiconductor nanowire, it is required
that the nanowire should possess a large g-factor and strong spin-orbit
interaction (SOI). A topological superconducting nanowire with MBSs present at its two ends can be created by placing such a semiconductor nanowire in the proximity of an
s-wave superconductor and under application of a moderate magnetic
field <cit.>.

Stimulated by the experimental progress <cit.>, many theoretical proposals
<cit.>
to implement TQC with MBSs in nanowire devices are reported. Because
braiding two MBSs, the essential operation to realize the TQC, is
not straightforward in one dimension, a solution to walk around
this difficulty is to build a nanowire network for MBS braiding or exchanging.
Before studying the actual physical process of MBS braiding in a nanowire network, the first natural step is to study the energy spectra of MBSs and subgap states in a minimal structure, a
three-terminal hybrid nanowire junction with the nanowires in the topological superconducting
phase. Various nanowire junctions <cit.>
have been realized and characterized by transport measurements. Several theoretical studies on multiterminal superconductor junctions have also been reported <cit.>.

In this paper, we report on a theoretical study of three-terminal topological superconducting nanowire (TSNW) junction devices with the nanowires connected to each other via a quantum dot (QD) in the junction region in each device. We explore the effects of the couplings among MBSs presented in the junction devices. Since, as shown in a recent work,<cit.> a QD can be an important constitution to realize
MBS braiding, here we include a QD in the devices. The QD is located at the connecting junction of the three nanowires and has an energy state which is assumed to be tunable using, e.g., a gate. Thus, the QD can play a role of high barrier (without
a QD level at the Fermi energy) or a conduction channel (with a QD
level at the Fermi energy) and can therefore be used to control the couplings among MBSs in the TSNWs. Our three-terminal TSNW junction devices are described by standard Bogoliubov-de Gennes (BdG) tight-binding Hamiltonians. The energy spectra and the probability distributions of the wave functions
of subgap states in the devices are calculated and analyzed. We have also extracted low-energy effective models that describe the results of the numerical calculations for the subgap states of the three-terminal TSNW-QD hybrid devices based on the BdG tight-binding Hamiltonians. We demonstrate that the coupling strengths between the MBSs adjacent to the QD in the three-terminal TSNW-QD hybrid devices are functions of the
azimuths of the TSNWs, the QD  level energy, and the superconducting phase differences between the TSNWs. Tuning superconducting phase differences could be achieved by employing local magnetic fields to SQUID setups connecting the TSNWs.  A similar but different all-magnetic-flux-controlled prototype device for MBS braiding is proposed in Refs. njp14.035019,prb88.035121,
in which the magnetic flux is used to control the coupling of MBSs in a single
nanowire. Here, the magnetic flux is used to tune the couplings between MBSs in
different nanowires.

The paper is organized as follows. In Sec. II, we present the BdG tight-binding model that is used to describe our three-terminal TSNW-QD hybrid junction devices.
In Sec. III, we present the results of the calculations for the three-terminal TSNW-QD hybrid junction devices based on the BdG tight-binding model. The energy spectra and wave functions of the subgap states are calculated in the parameter space of the phase differences between TSNWs for two cases, i.e., the case when there is no QD level located near or inside the superconducting gap and the case when the QD level is located near or inside the superconducting gap. Simple but general effective models extracted for describing the subgap states of the three-terminal TSNW-QD hybrid devices in the two cases will also be presented and discussed. Finally, in Sec. IV, we summarize and conclude the paper.




§ MODEL AND FORMALISM




As depicted in Fig. <ref>(a), we consider a three-terminal TSNW-QD hybrid junction device which is placed on the
x-y plane with three TSNWs connected together via a
QD at the junction. An applied magnetic field along the z direction is assumed
to drive the three semiconductor nanowires with strong SOI and a sizable g factor into the topological superconducting phase. Figure <ref>(b) displays a lattice model of the system with the TSNWs
represented by discrete sites and the QD by one site at the center of the junction. The two of the TSNWs are assumed to be connected via superconductor leads and
ϕ_ij is the phase difference between TSNWs i and j, which
can be tuned by applying a local magnetic field. The coordinate system is established by taking
the QD site as the origin.

The total Hamiltonian includes three parts and can be written as

    H_total=∑_i=1^3H_i+H_QD+H_c,

where H_i is the Hamiltonian of the i-th TSNW, H_QD is the Hamiltonian of the QD in normal phase (i.e., without including
a superconducting pairing potential), and H_c describes the couplings between
the TSNWs and the QD. In order to take into account the spatial directions of
the wires, the azimuth of the i-th TSNW q_i
enters H_i following the formula developed in Ref. prb85.144501. In the tight-binding approach, the Hamiltonian H_i takes the form of

    H_i=∑_n,m1/2Ψ_n^†(ℋ_BdG^i)_n,mΨ_m,

where ℋ_BdG^i is the BdG Hamiltonian
for the i-th TSNW in the Nambu spinor basis of Ψ_n=[ψ_n↑,ψ_n↓,ψ_n↓^†,-ψ_n↑^†]^T,

    (ℋ_BdG^i)_n,m=h_0^iδ_n,m+h_1^iδ_n,m-1+h_-1^iδ_n,m+1,

with

    h_0^i    =     (2t-μ)τ_z+V_zσ_z+Δ(τ_xcosϕ_i-τ_ysinϕ_i),  
    
    h_1^i   =     [-t+iα_0(σ_ycos q_i-σ_xsin q_i)]τ_z,
    
    h_-1^i    =    [-t-iα_0(σ_ycos q_i-σ_xsin q_i)]τ_z.

Here, m and n are the site indices, t=ħ^2/2m^*a^2
is the hopping parameter, which is related to the band width, with ħ standing for the reduced
Planck constant, m^* the effective mass, and a the lattice
spacing, σ_u and τ_u, where u=x,y or z,  are the Pauli matrices acting, respectively,
on spin and particle-hole spaces,  α_0 is the Rashba SOI strength,
μ is the chemical potential, V_z is the Zeeman energy, Δ and ϕ_i
are the magnitude and phase of the superconducting pairing potential, and q_i is the
azimuth of the i-th TSNW.

The topological phase transition point for each nanowire is at V_z=√(Δ^2+μ^2),
with V_z<√(Δ^2+μ^2) for the nanowire in trivial phase and V_z>√(Δ^2+μ^2)
for the nanowire in topological superconducting phase <cit.>.
For simplicity, in the calculations and discussions presented below, the chemical
potential μ is set to 0 and V_z=2Δ, such that all the nanowires
are in topological superconducting phase.

The Hamiltonian of the QD in the basis of Ψ_QD=[ψ_QD↑,ψ_QD↓,ψ_QD↓^†,-ψ_QD↑^†]^T
is

    H_QD=1/2Ψ_QD^†ℋ_QDΨ_QD,

with

    ℋ_QD=-eV_gτ_z+V_zσ_z,

where -e is the single electron charge and V_g is the tunable
gate voltage.

The coupling Hamiltonian is

    H_c=∑_i=1^31/2Ψ_i^†T_iΨ_QD+H.c.,

where Ψ_QD and Ψ_i are the annihilation operators
at the QD site and its adjacent site in the i-th TSNW, and T_i=-t_ciτ_z
is the spin independent coupling matrix.




§ SUBGAP STATES AND EFFECTIVE HAMILTONIANS


As can be seen in the above formulation, the energy spectra and the coupling between MBSs adjacent to the QD are functions of the phase differences, ϕ_21=ϕ_2-ϕ_1 and ϕ_32=ϕ_3-ϕ_2 between TSNWs, and the energy position of the QD level. In the following, we first in Sec. III A present the results of the calculations for the case when there is no QD level in or near the superconducting gap. Here, we will show that the MBSs adjacent to the QD and the couplings between them can be tuned through tuning phase differences ϕ_21 and ϕ_32. In Sec. III B, we will present the results of calculations for the case when a QD level is present in or near the superconducting gap and show that the tunings of the MBSs and the couplings between them can be also achieved through controlling of the energy level in the QD. Simple effective model Hamiltonians for describing the results of the numerical calculations will also be extracted for the two cases. A brief discussion about the generality of the effective models will be given in Sec. III C.



 §.§ System with no QD level in and near the superconducting gap


In a three terminal TSNW-QD hybrid system with all QD levels being located far from the superconducting gap and, thus, far from the Fermi
level, the QD just acts as a potential barrier at the connection point of
the three TSNWs. Below, we will present the results of calculations for the energy spectra and the wave functions of the subgap states of the system by fixing ϕ_21 while varying ϕ_32. If the calculations are performed by fixing ϕ_32 and varying ϕ_21, the same results will be obtained.




Figure <ref> shows the calculated energy spectra of subgap state of the system at different values of ϕ_21 and ϕ_32. In Figs. <ref> (a1) and <ref>(a2), the subgap energy spectra of the system with V_z=2Δ at a fixed value of ϕ_21=0 and different values of ϕ_32 are plotted in a linear and a logarithmic scale, respectively. There are six subgap states
inside the superconducting gap. These states appear in pairs at energies of ∼±10^-7Δ, ∼±10^-4Δ and ∼±10^-1Δ, see Fig. <ref>(a2) for the three positive energy subgap states. Note that at V_z=2Δ, the superconducting gap parameter E_g shrinks to E_g≈0.7Δ. The wave functions of the three positive energy subgap states calculated at ϕ_21=0  and ϕ_32=π  are shown in Figs. <ref>(a3) to <ref>(a5). The probability distribution of the lowest positive energy subgap state at the energy of ∼ 10^-7Δ is
shown in Fig. <ref>(a3). It is seen that this state
is formed from the three MBSs located at the outer
ends of the three TSNWs. The probability distribution of the next lowest positive energy subgap state at the energy of ∼ 10^-3Δ is shown in Fig. <ref>(a4). This state is formed from a superposition of all the six MBSs. Note that the energies of these two states are very close to zero, but are finite due to the fact that the TSNWs are all finite in length. The probability distribution of the remaining positive energy subgap state at the energy of ∼ 10^-1Δ is shown in Fig. <ref>(a5). This state is formed from the hybridization of the MBSs at the inner ends of the three TSNWs. Thus, the energy of this state would not change with increasing TSNW lengths. However, the state can be effectively tuned by
the couplings t_ci between the TSNWs and the QD, and by the phase differences between the TSNWs, as we will show below. The three corresponding negative energy subgap states show the same probability distributions as their positive energy counterparts.

Figures <ref>(b1) and <ref>(b2) show the subgap energy spectra of the system at a fixed value of ϕ_21=2π/3 and different values of ϕ_32.
Figures <ref>(b3) to <ref>(b5) show the wave functions of the subgap states of the system at ϕ_21=2π/3 and ϕ_32=3π/5. Clearly, the wave functions of the six subgap states show the similar
probability distributions as their corresponding subgap states shown in Figs. <ref>(a3) to <ref>(a5). In fact, similar characteristics in the probability distributions of the six subgap states are observed at other values of ϕ_21 and ϕ_32, except for the case of ϕ_21=4π/3 and ϕ_32=4π/3 as we will show below.

Figures <ref>(c1) and <ref>(c2) show the subgap energy spectra of the system at a fixed value of ϕ_21=4π/3. Here, it is seen that the
energies of the subgap states are varied with varying ϕ_32 and all the subgap states move together in energy when the phase difference ϕ_32 moves to ϕ_32=4π/3. We note that this particular case is found only when ϕ_21=4π/3. We note also that at this particular case, all the subgap states are close to zero in energy as shown in Figs. <ref>(c1) and <ref>(c2)
and have the same form of the probability distributions as seen in
Figs. <ref>(c3) to <ref>(c5). It can also be observed that the subgap states within
each TSNW have the same form of probability distributions as that in an isolated single TSNW, which implies that
the couplings between the three TSNWs become vanishingly small at this specific case. These results indicate that the phase differences between the TSNWs can play a role in tuning of couplings between the three TSNWs and we
can control the couplings between the TSNWs via controlling phase differences between the three TSNWs.



Based on the numerical results presented above, an effective model involving only the MBSs
in the three-terminal TSNW-QD hybrid structure can be constructed and can be used for simplification in considering employing the MBSs in the system. Similar to the procedure employed for two MBSs coupled
through a weak link in a one-dimensional topological superconducting system <cit.>, the effective Hamiltonian of the three-terminal TSNW-QD hybrid structure can be written as

    H_eff=H_eff^S+H_eff^D,

with

    H_eff^S    =     i(t_14γ_1γ_4+t_25γ_2γ_5+t_36γ_3γ_6),
    
    H_eff^D    =     i(γ_2γ_1t_21sinβ_21/2+γ_3γ_2t_32sinβ_32/2
           +γ_1γ_3t_13sinβ_13/2),

where t_ij∈ℝ denotes the coupling between MBSs γ_i
and γ_j [see the schematic in Fig. <ref>(b) for the locations of the MBSs] and satisfies t_ij=-t_ji, and phase difference β_ji is defined as β_ji≡β_j-β_i=q_j-q_i+ϕ_j-ϕ_i and can be written explicitly as

    β_21    =    ϕ_21+q_2-q_1,
    β_32    =    ϕ_32+q_3-q_2,
    β_13    =     -β_21-β_32.

For an arbitrary junction, the azimuths q_1,2,3 are in the range [0,2π]. For the specific case drawn in Fig. <ref>(b), q_1=π/6, q_2=5π/6, and q_3=3π/2, therefore

    β_21    =    ϕ_21+2π/3,
    β_32    =    ϕ_32+2π/3,
    β_13    =     -ϕ_21-ϕ_32-4π/3.


There are two types of couplings between MBSs. One contains the couplings between two MBSs in the same TSNW (H_eff^S) and the other one contains the couplings between MBSs adjacent to the QD but belonging to different TSNWs (H_eff^D). Comparing with
the couplings (t_ij) between MBSs in H_eff^S, the couplings
(t_ijsinβ_ji/2) in H_eff^D has two parts–a constant amplitude part t_ij and an oscillatory part arising from
the phase difference between TSNWs. In contrast to a straight nanowire system with the two inner MBSs coupled with a weak link, where the phase difference
is solely determined by the magnetic flux through a circuit loop, the phase differences between different TSNWs
in the three-terminal TSNW-QD hybrid system could contain both terms induced by the magnetic fluxes and terms related to the azimuths of the TSNWs. The azimuths of the TSNWs enter the Hamiltonian due to the presence of SOI in the nanowires
<cit.>. It has been shown for a continuous system that
when a one-dimensional spinful Rashba s-wave superconducting nanowire model
is mapped to a spinless p-wave superconducting nanowire model, the effective phase of the system is given by the summation of the azimuth and the phase in the
s-wave superconducting pairing potential of the nanowire. Thus, the differences in these two parts constitute the total phase differences between the TSNWs.

Figure <ref> shows the results of calculations based on the effective Hamiltonian. It is seen that the results of the full numerical calculations based on the BdG tight-binding Hamiltonian are recovered by the calculations based on the effective Hamiltonian. In particular, the energy spectra and the probability distributions seen in Figs. <ref>(a1) to <ref>(a5) are well reproduced by the calculations based on the effective model shown in Figs. <ref>(a1) to <ref>(a5). Here, we should note that the wave function of each MBS obtained based on the effective Hamiltonian is a point-like function, while the wave function of a MBS calculated based on the BdG tight-binding model exhibits a spatial but strongly  localized distribution.
Figures <ref>(b1) to <ref>(b5) show the calculated
energy spectra for the case of ϕ_21=4π/3 and the probability distributions for the case of ϕ_21=ϕ_32=4π/3 based on the effctive Hamiltonian, corresponding to the results calculated based on the BdG tight-binding Hamiltonian shown in Figs. <ref>(c1) to <ref>(c5).
Here, it is seen that at ϕ_21=ϕ_32=4π/3 where the six MBSs become nearly degenerate, the three TSNWs are decoupled in the effective Hamiltonian and each MBS resembles exactly that in
an isolated TSNW. Note that these states are still
superpositions of MBSs within each TSNW, since the coupling of the MBSs
within each TSNW remains non-zero in order to simulate the finite length effect in
our effective model. Figure <ref>(c) shows a view of the energy spectra of the three-terminal TSNW-QD hybrid system in the whole
parameter space of ϕ_21 and ϕ_32 based on the effective Hamiltonian.

If the three TSNWs are semi-infinitely long, the overlap of two MBSs within each
TSNW can be neglected, the effective Hamiltonian can be reduced to

    H_SInf    =     i(γ_2γ_1t_21sinβ_21/2+γ_3γ_2t_32sinβ_32/2
        +    γ_1γ_3t_13sinβ_13/2).

In this limit, the three-terminal TSNW-QD hybrid structure consists of four zero-energy MBSs with three of them located at the outer ends of the three TSNWs and the remaining one in the vicinity of the QD.<cit.> There are also two finite energy states in the system located in the vicinity of the QD. These two states originate from hybridization of the MBSs adjacent to the QD. Explicitly, based on the Hamiltonian H_SInf given in Eq. (<ref>), the eigenenergies of the three states localized in the vicinity of the QD can be readily obtained as

    E_0 = 0      E_±  =  ± s_0/ 2,

where s_0=√(s_12^2+s_23^2+s_13^2) with s_ji=t_jisinβ_ji/2. The corresponding eigenvectors are

    |0⟩   =   1/s_0(s_32γ_1+s_13γ_2+s_21γ_3),
    
    |E_±⟩   =   1/s_N[
    (s_13s_32± is_21s_0)γ_1
    -(s_21^2+s_32^2)γ_2 
       +   (s_21s_13∓ is_32s_0)γ_3],

where the normalized constant s_N=s_0√(2(s_21^2+s_32^2)).



 §.§ System with a QD level inside or near the superconducting gap




We now study the three-terminal TSNW-QD hybrid device in the case when a QD level is tuned into or near the superconducting gap. In the BdG particle-hole representation, a QD level contributes two states, namely the particle state and its hole partner, with opposite energies. Thus, there could totally exist eight subgap states in the system, instead of six seen in the previous subsection. However, we will show that this can only occur when the QD level is located near or inside the superconducting gap. Below, we will first discuss how the subgap states evolute with change in the energy position of the QD level, which can be tuned as

    E_QD=-eV_g+V_z+δ_Γ,

through a gate voltage V_g, where δ_Γ is the normalized energy of the QD level due to the couplings with the three TSNWs. Numerical calculation shows that δ_Γ≈0.234Δ in the system we have considered with coupling parameters t_ci=Δ in this subsection. Then, we will study the tuning of the subgap states by shifting the phase differences between the TSNWs.

Figure <ref> shows the subgap states of the three-terminal TSNW-QD hybrid system at given values of the phase differences between the TSNWs, ϕ_21 and ϕ_32, at different QD level energies. In Fig. <ref>(a), the energy spectra calculated for the system with ϕ_21=π/3 and ϕ_32=3π/4 based on the full tight-binding Hamiltonian are plotted against the QD level energy (E_QD). Figure <ref>(b) shows the energy positions in logarithmic scale of the positive energy subgap states in the region marked by a rectangle in Fig. <ref>(a).
It is seen in Fig. <ref>(a) that when the QD level is located far from the superconducting gap, only six subgap states are found in the gap. When the QD level is moved into the superconducting gap, a region, in which eight subgap states are present in the gap, exists [see Fig. <ref>(b) for the existence of the four positive energy subgap states]. When continuing moving the QD level away from the superconducting gap, the system hosts only six subgap states again and the results presented in the previous subsection are recovered.

Figures <ref>(c1) to <ref>(c4) show the wave functions of the four positive energy subgap states of the system at the QD level energy of E_QD=-0.766Δ, while Figs. <ref>(d1) to <ref>(d4) show that at E_QD=0.
It is seen in Fig. <ref>(c1) and Fig. <ref>(d1) that the lowest energy subgap state is built from a superposition of the three MBSs at the outer ends of the three TSNWs. Thus, the state is very close to zero in energy and does not show a significant change with shifting QD level. The second lowest subgap state is built from a superposition of all the six MBSs in the TSNWs and also does not show a significant change in energy with shifting QD level. The third and fourth lowest energy subgap states are built mainly from superpositions of the three MBSs adjacent to the QD and the QD state. Thus, the energies of these two states show clearly visible changes with shifting QD level position. However, the behaviors of these two subgap states are different. The highest subgap state moves towards a lower energy when the QD level moves towards the center of the superconducting gap, and is then moving back to higher energies and eventually to the outside of the superconducting gap with continuously increasing QD level energy. The wave function of this state is seen to be slightly more extended for the QD level at E_QD=-0.766Δ [see Figs. <ref>(c4-i)] than at E_QD=0 [see Figs. <ref>(d4-i)]. This is due to the fact that the subgap state is located deeper in energy in the superconducting gap at E_QD=0 than at E_QD=-0.766Δ. In fact, this subgap state becomes more extended when the QD level is moved to a higher energy and eventually develops into an extended state when the QD level  merges into continuous band of quasiparticle states. On contrast, the third lowest subgap state moves to a higher energy when the QD level moves towards the center of the superconducting gap and then moves back to a lower energy after the QD level passes through the center of the gap. The wave function of this subgap state remains localized to the vicinity of the QD and does not show a significant change in localization with change in the QD level energy.



We now discuss the evolutions of the subgap states with changes in phase differences between TSNWs at different given QD level energies. Figure <ref> shows the results of calculations based on the BdG tight-binding Hamiltonian for the system at the QD level energies of E_QD=0.234Δ [Figs. <ref>(a1) to <ref>(a6)] and E_QD=0 [Figs. <ref>(b1) to <ref>(b6)]. At both energies, the system hosts eight subgap states. it is seen in Figs. <ref>(a1) to <ref>(a2), where E_QD=0.234Δ and ϕ_21=π, that the energies of the eight subgap states only change slightly with varying ϕ_32. The four positive energy subgap states are located at ∼ 10^-7Δ, ∼ 10^-4Δ, ∼ 10^-1Δ  and ∼ 10^-1Δ in energy [see Fig. <ref>(a2)]. The wave functions of the four positive energy subgap states calculated at ϕ_21=π and ϕ_32=π/3 are shown in Figs. <ref>(a3) to <ref>(a6). It is seen in Fig. <ref>(a3) that the lowest positive energy subgap state at the energy of ∼ 10^-7Δ is formed from the three MBSs located at the outer ends of the three TSNWs. The next lowest positive energy subgap state at the energy of ∼ 10^-4Δ is formed from a superposition of all the six MBSs [see Fig. <ref>(a4)]. The reason for that the energies of these two subgap states are very close to zero but are still finite is the same as in the previous case, namely that the TSNWs are all finite in length. The probability distributions of the remaining two positive energy subgap states at the energies of ∼ 10^-1Δ are shown in Figs. <ref>(a5) and <ref>(a6). These two subgap states are formed from the hybridizations of the MBSs at the inner ends of the three TSNWs and the QD state. The probability distributions of all the eight subgap states vary little when ϕ_21 and ϕ_32 are changed, except for the case when ϕ_21 and ϕ_32 are in the vicinity of ϕ_21=4π/3 and ϕ_32=4π/3.

Figure <ref>(b1) shows the evolution of the subgap state energy spectra at E_QD=0 and ϕ_21=4π/3 with change in ϕ_32, while Fig. <ref>(b2) shows the energy evolution of the four lowest positive energy subgap states in logarithmic scale. It is seen that the two lowest positive energy subgap states at the energies of ∼ 10^-9Δ and ∼ 10^-4Δ are very close to zero and change little with change in phase difference ϕ_32. The highest positive energy subgap state at the energy of ∼ 10^-1Δ also exhibits a weak dependence on phase difference ϕ_32.  The remaining, third lowest positive energy subgap state shows, however, a different behavior–it varies slowly in energy with increasing ϕ_32 from zero and turns to move quickly towards smaller energy when ϕ_32 approaches 4π/3. The energy of this subgap state reaches the minimum of ∼ 10^-4Δ at ϕ_32=4π/3.
Figures <ref>(b3) to <ref>(b6) show the probability distributions of the four positive energy subgap states of the system at ϕ_21=ϕ_32=4π/3 and E_QD=0. It is seen in Fig. <ref>(b3) that the lowest positive energy subgap state at the energy of ∼ 10^-9Δ is built dominantly from the QD state and the three MBSs at the outer ends of the three TSNWs. Note that three small but broad peaks are visible at the outer ends of the three TSNWs in Fig. <ref>(b3).
The second and third lowest positive energy subgap states at the energies of ∼ 10^-4Δ, seen in Figs. <ref>(b4) and <ref>(b5), are formed dominantly from the MBSs in the three TSNWs. By a close examination, we also find that the probability distributions of these two subgap states in each TSNW have almost the same form of the probability distribution as in an isolated single TSNW, similar to the case discussed in the previous subsection. The  fourth lowest positive energy subgap state shown in Fig. <ref>(b6) is formed primarily from the QD state and the three MBSs adjacent to the QD, similar to the subgap state shown in the Fig. <ref>(a6).



Based on the numerical results presented above, we repeat the procedure in the previous subsection to derive an effective model for describing the subgap states of the system by involving only the MBSs and the QD state as

    H_eff=H_eff^S+H_eff^D+H_eff^QD,

where H_eff^S describes the couplings of the MBSs in the same TSNW,
H_eff^D is the couplings between the MBSs and the QD state, and H_eff^QD
is the Hamiltonian of the QD,

    H_eff^S    =     i(t_14γ_1γ_4+t_25γ_2γ_5+t_36γ_3γ_6),
    
    H_eff^D    =    ∑_i=1^3t_QDi(e^iβ_i/2d^†γ_i+H.c.),
    
    H_eff^QD    =     E_QDd^†d,

with the phases (not unique, up to a global phase background),

    β_1    =     q_1,
    β_2    =     q_2+ϕ_21,
    β_3    =     q_3+ϕ_21+ϕ_32.

Similarly, the azimuths q_1,2,3 are in the range [0,2π] for an arbitrary junction. Here, in the particular case shown in Fig. <ref>(b), q_1=π/6, q_2=5π/6, and q_3=3π/2, therefore

    β_1    =    π/6,
    β_2    =    5π/6+ϕ_21,
    β_3    =    3π/2+ϕ_21+ϕ_32.

By solving the effective model Hamiltonian, the results of the BdG tight-binding calculations
can be reproduced. Figure <ref> shows the results obtained based on the effective model Hamiltonian. When compared with Figs. <ref>(a1) to <ref>(a6), it is seen that the energy spectra and the probability distributions are well recovered in the calculations based on the effective model Hamiltonian of Eq. (<ref>).
Figure <ref>(b) shows the energy spectra of the system at phase differences ϕ_21=π/3 and ϕ_32=3π/4 calculated against the QD level energy E_QD based the effective model Hamiltonian of Eq. (<ref>). It is seen that the main features in the energy spectra of all the subgap states seen in Fig. <ref>(a) are found in the calculations based on the effective model Hamiltonian. In Fig. <ref>(c), we show the energy spectra of the system with the QD level at E_QD=0.234Δ in the whole parameter space of ϕ_21 and ϕ_32 based on the effective model Hamiltonian of Eq. (<ref>).

For the TSNW-QD hybrid system with three semi-infinitely long TSNWs, H_eff^S can be neglected and the effective model Hamiltonian reads

    H_SInf    =     H_eff^D+H_eff^QD
        =     E_QDd^†d+∑_i=1^3t_QDi(e^iβ_i/2d^†γ_i+H.c.).

Analytically solving the secular equation of the effective Hamiltonian in this infinitely long TSNW limit, we obtain eight subgap states with four zero-energy states and four finite-energy states. Similar to the case presented in the previous subsection, among the four zero-energy subgap states, three are the MBSs located at
the outer ends of the three TSNWs (γ_4, γ_5 and γ_6) and the remaining one is formed from a superposition
of the three inner MBSs adjacent to the QD as

    E_0=0,  
    |E_0⟩=1/l_N(l_32γ_1+l_13γ_2+l_21γ_3),

where l_ij≡ t_QDit_QDjsinβ_ij/2, and
the normalized constant l_N≡√(l_21^2+l_32^2+l_13^2).
The four finite-energy states are formed by hybridizations of the QD state and the three inner MBSs adjacent to the QD. Mathematically, these states are the solutions of a
quartic equation which can be simplified, by setting x≡ E^2, to a quadratic equation,

    x^2-Mx+N=0,

where

    M=   E_QD^2+2(t_QD1^2+t_QD2^2+t_QD3^2)/4,
    
       N=l_N^2/4.

Thus, the four finite-energy states have the eigenenergies of

    E_s_1,s_2=s_1√(M+s_2√(M^2-4N)/2),

where s_1,s_2=±. However, the eigenvectors |E_s_1,s_2⟩ of these four finite-energy states do not possess
a simple expression. Thus, it is more convenient to resort to numerical calculations for the eigenvectors of these states.



 §.§ Generality of the effective model Hamiltonians


In the above subsections, the calculations based on the effective model Hamiltonians for the three-terminal TSNW-QD hybrid structures with the orientations of the three TSNWs set at q_1=π/6, q_2=5π/6, and q_3=3π/2 are presented and discussed.
However, the application of the effective model Hamiltonians to a three-terminal TSNW-QD hybrid structure with arbitrary TSNW orientations is straightforward. In fact,
the effective model Hamiltonians given by Eqs. (<ref>) and (<ref>) and by Eqs. (<ref>) and (<ref>) are in their general forms and can be readily used for the study of a three-terminal TSNW-QD hybrid structure with arbitrary TSNW orientations. For example, in the energy spectra shown in Figs. <ref>(c) and <ref>(c) obtained  based on the effective model Hamiltonians in the full parameter space of phase differences ϕ_21 and ϕ_32,  a specific point at which the six low energy subgap states are very close in energy (all in the orders of ∼±10^-4Δ or smaller). This point is seen to appear at ϕ_21=ϕ_32=4π/3 in Figs. <ref>(c) and <ref>(c). For a three-terminal TSNW-QD hybrid structure with arbitrary TSNW orientations, such a specific point always exists but shifts to

    ϕ_21    =     2π+q_1-q_2,
    ϕ_32    =     2π+q_2-q_3.

For a T-shaped three-terminal TSNW-QD hybrid structure <cit.>, for which q_1=0, q_2=π, and q_3=3π/2, the specific point is found to appear at ϕ_21=π and ϕ_32=3π/2. For a parallel junction TSNW-QD hybrid system<cit.>, for which q_1=0, q_2=π, and
q_3=2π, the specific point shifts to ϕ_21=ϕ_32=π.
Thus, our effective model Hamiltonians given in Eqs. (<ref>) and (<ref>) and in Eqs. (<ref>) and (<ref>) are very general and can be applied to the calculations for a three-terminal TSNW-QD hybrid device with arbitrary nanowire orientations.



§ CONCLUSIONS


In conclusion, the subgap states in three-terminal TSNW-QD hybrid junction devices are studied.
The energy spectra and the wave functions  of the subgap states
are calculated as a function of the superconducting phase
differences between TSNWs and as a function of the QD level energy based on the BdG tight-binding Hamiltonians. In a system with no QD level located near or inside the superconducting band gap, there can exist six subgap states. Among them, two lowest energy (one positive and one negative energy) subgap states stay at nearly zero energies and are essentially formed by linear combinations of the three MBSs located at the outer ends of the three TSNWs. The next two lowest energy subgap states are built from linear combinations of all the six MBSs in the system and stay also at very small energies. However, the remaining two subgap states are mainly formed by linear combinations of the three MBSs adjacent to the QD and have significantly large finite energies. The calculations also show the existence of a specific point in the parameter space of phase differences between TSNWs, at which all the six subgap states move close to each other towards nearly zero energies. In a three-terminal TSNW-QD hybrid junction device with a QD level located near or inside the superconducting gap, there can appear eight subgap states. Among them, six low energy subgap states exhibit similar characteristics in energy positions and wave function localizations as in the system without a QD level located near or inside the superconducting gap except for the case where the QD level approaches zero energy. At this case, the two of the six subgap states at high energies contain significant contributions from the QD state. The remaining two large energy subgap states are the inherent property of the system with the QD level located near or inside the superconducting gap. These two states are constructed as linear combinations of the QD state and the three MBSs adjacent to the QD.

Simple but general effective model Hamiltonians for three-terminal TSNW-QD hybrid devices have also been extracted. The calculations based on the effective model Hamiltonians recover excellently the results of the calculations based on the full BdG tight-binding Hamiltonians. The effective model Hamiltonians can easily be generalized to the limit cases where the three TSNWs are  infinite in length. In such limit cases, some analytical solutions of the subgap states can be derived. The effective model Hamiltonians could be used as a starting point to construct and  investigate the braiding schemes of MBSs in TSNW junction devices.




§ ACKNOWLEDGMENTS


This work was supported by the Ministry of Science and Technology
of China (MOST) through the National Key Research and Development Program of
China (Grant Nos. 2016YFA0300601 and 2017YFA0303304),
the National Natural Science Foundation of China (Grant Nos. 91221202, 91421303, and 11604005), and the Swedish
Research Council (VR). GYH would also like to acknowledge financial
support from the China Postdoctoral Science Foundation (Grant No. 2016M591001)
to this work.
99

pu44.131A. Y. Kitaev, Phys. Usp. 44, 131 (2001).

np5.614F. Wilczek, Nature Phys. 5, 614 (2009).

rmp87.137Steven R. Elliott and Marcel Franz, Rev. Mod.
Phys. 87, 137 (2015).

sst27.124003Martin Leijnse and Karsten Flensberg, Semicond. Sci. Technol. 27, 124003 (2012).

nnano8.149M. Franz, Nature Nano. 8, 149 (2013).

arcmp4.113C. W. J. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).

rpp75.076501J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).

rmp80.1083C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).

ap303.2A. Kitaev, Ann. Phys. 303, 2 (2003).

science336.1003V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).

nl12.6414M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff and H. Q. Xu, Nano Lett. 12, 6414 (2012).

np8.795Leonid P. Rokhinson, Xinyu Liu and Jacek K. Furdyna, Nature Phys. 8, 795 (2012).

prb87.241401H. O. H. Churchill, V. Fatemi, K. Grove-Rasmussen, M. T. Deng, P. Caroff, H. Q. Xu, and C. M. Marcus, Phys. Rev. B 87, 241401 (2013).

sr4.7261M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Sci. Rep. 4, 7261 (2014).

arxiv1603.04069Hao Zhang et al., arXiv:1603.04069.

np8.887Anindya Das, Yuval Ronen, Yonatan Most, Yuval Oreg, Moty Heiblum, and Hadas Shtrikman, Nature Phys. 8, 887 (2012).

prl110.126406A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Phys. Rev. Lett. 110, 126406 (2013).

nature531.206S. M. Albrecht, A. P. Higginbotham, M. Madsen, F. Kuemmeth, T. S. Jespersen, J. Nygard, P. Krogstrup, and C. M. Marcus, Nature 531, 206 (2016).

science354.1557M. T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygard, P. Krogstrup, C. M. Marcus, Science 354, 1557 (2016).

science346.602Stevan Nadj-Perge, Ilya K. Drozdov, Jian Li, Hua Chen, Sangjun Jeon, Jungpil Seo, Allan H. MacDonald, B. Andrei Bernevig, and Ali Yazdani, Science 346, 602 (2014).

npjqi2.16035Remy Pawlak, Marcin Kisiel, Jelena Klinovaja, Tobias Meier, Shigeki Kawai, Thilo Glatzel, Daniel Loss, and Ernst Meyer, NPJ Quantum Inf. 2, 16035 (2016).

np13.286Benjamin E. Feldman, Mallika T. Randeria, Jian Li, Sangjun Jeon, Yonglong Xie, Zhijun Wang, Ilya K. Drozdov, B. Andrei Bernevig, and Ali Yazdani, Nature Phys. 13, 286 (2017).

prl105.077001R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).

prl106.127001R. M. Lutchyn, Tudor D. Stanescu, and S. Das Sarma, Phys. Rev. Lett. 106, 127001 (2011).

prl105.177002Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).

prb84.144522Tudor D. Stanescu, Roman M. Lutchyn, and S. Das Sarma, Phys. Rev. B 84, 144522 (2011).

pra82.052322J. D. Sau, S. Tewari, and S. Das Sarma, Phys. Rev. A 82, 052322 (2010).

prb84.035120D. J. Clarke, J. D. Sau, and S. Tewari, Phys. Rev. B 84, 035120 (2011).

prb85.144501Bertrand I. Halperin, Yuval Oreg, Ady Stern, Gil Refael, Jason Alicea, and Felix von Oppen, Phys. Rev. B 85, 144501 (2012).

np7.412J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nature Phys. 7, 412 (2011).

prb84.094505J. D. Sau, D. J. Clarke, and S. Tewari, Phys. Rev. B 84, 094505 (2011).

prb91.174305Cassio Sozinho Amorim, Kazuto Ebihara, Ai Yamakage, Yukio Tanaka, and Masatoshi Sato, Phys. Rev. B 91, 174305 (2015).

pra87.022343M. Burrello, B. van Heck, and A. R. Akhmerov, Phys. Rev. A 87, 022343 (2013).

njp14.035019B. van Heck, A. R. Akhmerov, F. Hassler, M. Burrello, and C. W. J. Beenakke, New J. Phys. 14, 035019 (2012).

prb88.035121T. Hyart, B. van Heck, I. C. Fulga, M. Burrello, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev. B 88, 035121 (2013).

prb91.201102Torsten Karzig, Falko Pientka, Gil Refael, and Felix von Oppen, Phys. Rev. B 91, 201102(R) (2015).

prb94.014511Xin Liu, Xiaopeng Li, Dong-Ling Deng, Xiong-Jun Liu, and S. Das Sarma, Phys. Rev. B 94, 014511 (2016).

prb94.035424Michael Hell, Jeroen Danon, Karsten Flensberg, and Martin Leijnse, Phys. Rev. B 94, 035424 (2016).

prx6.031016David Aasen, Michael Hell, Ryan V. Mishmash, Andrew Higginbotham, Jeroen Danon, Martin Leijnse, Thomas S. Jespersen, Joshua A. Folk, Charles M. Marcus, Karsten Flensberg, and Jason Alicea, Phys. Rev. X 6, 031016 (2016).

nl8.1100Dmitry B. Suyatin, Jie Sun, Andreas Fuhrer, Daniel Wallin, Linus E. Froberg, Lisa S. Karlsson, Ivan Maximov, L. Reine Wallenberg, Lars Samuelson and H. Q. Xu, Nano Lett. 8, 1100 (2008).

nnano8.859Sebastien R. Plissard et al., Nat. Nano. 8, 859 (2013).

nl13.5190Jung-Hyun Kang, Yonatan Cohen, Yuval Ronen, Moty Heiblum, Ryszard Buczko, Perla Kacman, Ronit Popovitz-Biro, and Hadas Shtrikman, Nano Lett. 13, 5190 (2013).

am26.4875Diana Car, Jia Wang, Marcel A. Verheijen, Erik P. A. M. Bakkers, and Sebastien R. Plissard, Adv. Mater. 26, 4875 (2014).

aem2.1500460Sebastian Heedt, Daniil Vakulov, Torsten Rieger, Daniel Rosenbach, Stefan Trellenkamp, Detlev Grutzmacher, Mihail Ion Lepsa, and Thomas Schapers, Adv. Electron. Mater 2, 1500460 (2016).

nl16.1933Torsten Rieger, Daniel Rosenbach, Daniil Vakulov, Sebastian Heedt, Thomas Schapers, Detlev Grutzmacher, and Mihail Ion Lepsa, Nano Lett. 16, 1933 (2016).

nl17.531Jung-Hyun Kang, Marta Galicka, Perla Kacman, and Hadas Shtrikman, Nano Lett. 17, 531 (2017).

nl17.2596Johannes Gooth, Mattias Borg, Heinz Schmid, Vanessa Schaller, Stephan Wirths, Kirsten Moselund, Mathieu Luisier, Siegfried Karg, and Heike Riel, Nano Lett. 17, 2596 (2017).

arXiv1703.05195Elham M. T. Fadaly et al., arXiv:1703.05195

arXiv1705.01480Sasa Gazibegovic et al., arXiv:1705.01480

prb92.155437Tomohiro Yokoyama and Yuli V. Nazarov, Phys. Rev. B 92, 155437 (2015).

prb92.205409C. Padurariu, T. Jonckheere, J. Rech, R. Melin, D. Feinberg, T. Martin, and Yu. V. Nazarov, Phys. Rev. B 92, 205409 (2015).

ncomms7.11167Roman-Pascal Riwar, Manuel Houzet, Julia S. Meyer, and Yuli V. Nazarov, Nat. Comm. 7, 11167 (2016).

prb95.235305Torsten Karzig, Christina Knapp, Roman M. Lutchyn, Parsa Bonderson, Matthew B. Hastings, Chetan Nayak, Jason Alicea, Karsten Flensberg, Stephan Plugge, Yuval Oreg, Charles M. Marcus, and Michael H. Freedman, Phys. Rev. B 95, 235305 (2017).



