Higher-order topological insulator in cubic semiconductor quantum wells

The search for exotic new topological states of matter in widely accessible materials, for which the manufacturing process is mastered, is one of the major challenges of the current topological physics. Here we predict higher order topological insulator state in quantum wells based on the most common semiconducting materials. By successively deriving the bulk and boundary Hamiltonians, we theoretically prove the existence of topological corner states due to cubic symmetry in quantum wells with double band inversion. We show that the appearance of corner states does not depend solely on the crystallographic orientation of the meeting edges, but also on the growth orientation of the quantum well. Our theoretical results significantly extend the application potential of topological quantum wells based on IV, II–VI and III–V semiconductors with diamond or zinc-blende structures.

The left-hand and right-hand solid curves correspond to the crossing between E1-H1 subbands and E2-H2 subbands, respectively. These curves divide the plane into three parts with trivial band ordering corresponding to band insulator (BI, see the left-hand white region), single band inversion (grey and blue regions) and double-band inversion (right-hand white region). The striped region defines a semimetal (SM) phase with vanishing indirect band-gap 50,51 . (E,F) Band structure calculated on the basis of effective 2D lowenergy Hamiltonian for the QWs with the layer thicknesses marked by the red symbols in (B) and (C). The blue and red curves represent band dispersion of electron-like and hole-like subbands, respectively. The wave vector is oriented along (100) crystallographic direction. The dotted curves represent the calculations based on realistic multi-band k · p Hamiltonian 49 .
Scientific Reports | (2021) 11:21060 | https://doi.org/10.1038/s41598-021-00577-z www.nature.com/scientificreports/ phase and the insulator state with a double band inversion. In this work, we identify a double-band-inversion insulator state 2D HOTI with the corner states arising due to cubic symmetry of II-VI and III-V semiconductors.
To describe double band inversion at the Ŵ point of the Brillouin zone, we derive an effective 2D low-energy Hamiltonian taking into account E1, E2, H1 and H2 subbands and preserving the cubic symmetry of the prototype QWs. Starting from a realistic multi-band k · p Hamiltonian for (0mn)-oriented cubic semiconductor QWs 49 and following expansion procedure described in the Supplementary Materials, in the basis |E1+� , |H1+� , |H2−� , |E2−� , |E1−� , |H1−� , |H2+� , |E2+� , the effective 2D Hamiltonian in the vicinity of the Ŵ point has the form: where the asterisk denotes complex conjugation. The diagonal blocks of H 2D (k x , k y ) are split into isotropic and anisotropic parts: The isotropic part H (i) 4×4 (k x , k y ) is written as 52,53 : where Here, k ± = k x + ik y , k x and k y are the momentum components in the QW plane, and C 1,2 , M 1,2 , A 1,2 , B 1,2 , D 1,2 , S 0 and R (i) 1,2 are isotropic structure parameters being defined by the QW geometry, the growth orientation and the materials. The Hamiltonian H 2D (k x , k y ) has a block-diagonal form because we keep the inversion symmetry (see Supplementary Materials) by neglecting the terms resulting from the anisotropy of chemical bonds at the QW interfaces 54 and the bulk inversion asymmetry of the unit cell of zinc-blende semiconductors 55 . The latter is absent for diamond-like semiconductors.
The most important quantities in H 2D (k x , k y ) are two mass parameters M 1 and M 2 describing the band inversion between E1-H1 subbands and E2-H2 subbands, respectively. The trivial band insulator corresponds to positive values of M 1 and M 2 . The QSHI and "bilayer graphene" states arise if M 1 < 0 and M 2 > 0 , and the difference between these states is defined by the gap between H1 and H2 subbands, which is zero in the case of "bilayer graphene" state 51 . The insulator state with double band inversion is defined by the negative values of M 1 and M 2 . We note that since the semimetal phase represented by the striped areas in the diagrams is formed by non-local overlapping of the valence and conduction subbands, it cannot be described within the low-energy Hamiltonian for the small values of k x and k y .
The isotropic term H (i) 4×4 (k x , k y ) preserves the rotational symmetry in the QW plane, therefore it is independent of the orientation of x and y axis. In contrast, the form of the anisotropic term H (a) 4×4 (k x , k y ) in Eq. (2), resulting from the cubic symmetry of diamond and zinc-blende semiconductors, depend not only on the QW growth orientation but also on the orientation of x and y axis (see Supplementary Materials). For (0mn)-oriented QWs, H (a) 1,2 and R (a) 1,2 are cubic structure parameters, which depend on the QW geometry and materials. In Eq. (5), θ = arctan(m/n) is the angle defining the QW growth orientation, while the angle ϕ is the angle between the x axis and the (001) crystallographic direction.
As clear from above, depending on the structure parameters, the effective 2D Hamiltonian in Eq. (1) describes QSHI, "bilayer-graphene" state, trivial band insulator or insulator with double band inversion. Note that thirteen parameters involved in H 2D (k x , k y ) cannot take arbitrary values. Since their calculation is based on the wavefunctions of the multi-band k · p Hamiltonian at k x = k y = 0 (for details, see Supplementary Materials), the set of parameters corresponding to specific topological state is determined by the thicknesses and materials of the QW layers, its growth orientation, and the buffer on which the QW is grown. The latter is crucial for taking into account the effect of lattice-mismatch strain on the band structure in the QW. Thus, all of the structure parameters are the functions of d, t, m, n and the QW layer and buffer materials. Further, we perform the calculations for two sets of structure parameters involved in H 2D (k x , k y ) , which correspond to the prototype QWs with the layer thicknesses marked by the red symbols in Fig. 1C,D. Figure 1E,F compare the band structure calculations based on realistic multi-band k · p Hamiltonian 49 and H 2D (k x , k y ) in Eq. (1). The isotropic and cubic structure parameters of H 2D (k x , k y ) for the insulator state with double band inversion are given in the Supplementary Materials.

Anisotropic edge states.
Let us now analyze the edge states arising in the insulator with double band inversion. Since H 2D (k x , k y ) in Eq. (1) has block-diagonal form, we further focus on the upper block H 4×4 (k x , k y ) only, while the calculations for the lower block H * 4×4 (−k x , −k y ) can be performed in the similar manner. To derive an effective 1D low-energy edge Hamiltonian, we split H 2D (k x , k y ) into two parts so that the first part represents two independent BHZ-like models 2 with M 1 < 0 and M 2 < 0 for the pairs of E1-H1 subbands and E2-H2 subbands, while the second part includes the rest isotropic and cubic terms describing the interpairs mixing. Then, assuming open-boundary conditions in a semi-infinite plane y > 0 , we solve the eigenvalue problem for the independent BHZ-like blocks to find the edge wave-functions at k x = 0 . In this case, the edge orientation represented by the x axis is defined by the angle ϕ measured from (100) crystallographic direction. Finally, we construct a low-energy edge Hamiltonian by projecting H 4×4 (k x , k y ) onto the obtained set of the basis edge functions (see Supplementary Materials).
The projection of two independent BHZ-like blocks 2 with non-zero k x leads to where η 2 n = (B n + D n )/(B n − D n ) with n = 1 and 2 corresponding to the pairs of E1-H1 and E2-H2 subbands, respectively. One can see that H where δk = k x − k c , I 2 is a 2 × 2 identity matrix, σ z is one of the Pauli matrices, and ε 0 is a constant corresponding to the energy of the crossing point at k x = k c . In Eq. (8), v 0 and v z are written as Note that the crossing of other Kramer's partners occurs at k x = −k c .
The projection of the rest terms of H 4×4 (k x , k y ) representing the mixing between the pairs |E1, +�-|H1, +� and |E2, −�-|H2, −� results in anti-diagonal mass terms describing the band-gap opening. After straightforward calculations with the details provided in the Supplementary Materials, we finally obtain the low-energy 1D edge Hamiltonian: where v x , v y , m x , m y , δ x , δ y include the angle dependence on θ and ϕ: www.nature.com/scientificreports/ Here, κ 1 and κ 2 are defined by the matrix elements ∂/∂y and ∂ 2 /∂y 2 , respectively, both calculated by using the basis edge functions at k x = 0 (see Supplementary Materials). In Eqs. (10), we have also introduced the edge isotropic ( F i , F 0 ) and cubic ( F a , F a ) parameters: The similar calculations for the block H *  Figure 2 shows the edge state picture of the prototype QWs with double band inversion based on low-energy 1D edge Hamiltonian for the edge oriented along (100) crystallographic direction, i.e. ϕ = 0 . One can see that the edge dispersion of each Kramer's partner mimics the dispersion of "tilted" 1D massive Dirac fermions. However, an important distinctive property of the edge states shown in Fig. 2 is that their gap is described simultaneously by two mass m x and m y parameters that prevent the gap vanishing at specific edge orientation (cf. Refs. [35][36][37][38]. Indeed, if one neglects the cubic terms of H 2D (k x , k y ) , resulting in F a =F a = 0 in Eq. (10), the mass parameter m y vanishes, while m x becomes independent of the edge orientation. Nevertheless, even such complex structure of the edge states yields the corner states in prototype QWs with double band inversion.
0D corner states and boundary conditions. To calculate the energy of the corner states, we apply linear approximation of the effective low-energy 1D edge Hamiltonian H 1D (δk, θ , ϕ) . Before going further, we should make a certain remark simplifying the calculations. With the parameters provided in the Supplementary Materials, one can verify that both v x and v y are significantly lower than v 0 and v z for any orientation of the edges (also (10) − F a sin 2 2θ 2k c cos 2ϕ sin 2 ϕ + κ 1 sin 2 2ϕ , v y = 2F a sin 4ϕ(κ 1 + k c ) + F a sin 2ϕ sin 2 2θ 2k c sin 2 ϕ − κ 1 cos 2ϕ , δ x = F i − F a cos 4ϕ − F a cos 2ϕ sin 2 ϕ sin 2 2θ , δ y = F a sin 4ϕ + F a sin 2ϕ sin 2 ϕ sin 2 2θ .
(11) www.nature.com/scientificreports/ see Fig. 3). Therefore, one can neglect these terms in the first approximation and take them into account by means of perturbation theory. Thus, by applying the unitary transformation, the linear part of Eq. (9) can be written as It is clear that Eq. (12) represents a 1D Dirac Hamiltonian with the mass −m y , which changes its sign with ϕ (see Fig. 3), modified by the presence of "tilted" term v 0 δkI 2 and second mass term m x σ x . Let us now define the coordinate x along the curved edge so that x = 0 corresponds to the meeting corner as shown in Fig. 3A. The latter means that m x , m y in Eq. (12) are the function of x, and δk = −i∂/∂x . Note that H 1D (−i∂/∂x) is defined in disjoint regions out of x = 0 . To fully define the 1D system, one needs to specify the boundary conditions that the wave functions must satisfy in the vicinity of x = 0.
As shown in the Supplementary Materials, the general linear boundary condition, conserving the probability current along the curved edge, can be written in the form where 1 and 2 are the wave-functions defined from different sides of the corner, while β 1 and β 2 are real parameters lying in the range from −π/2 to π/2 . A physical interpretation of these parameters will be discussed later.
On the basis of Eq. (13), it is convenient to introduce a new Hamiltonian H where α and m are real constants defined as where C is the normalization constant and Note that in Eqs. (24) and (25), the sign of σ should be chosen in accordance with normalized condition of � 0D (x) . If W (+∞) > 0 , σ = −1 , while for W (+∞) < 0 , σ = 1 . These two cases correspond to the internal and external corners at the same orientations of the two edges.
Let us make few remarks concerning the results obtained above. First, Eq. (23) can be also written in equivalent form which is reduced to the well-known condition for the existence of the bound state in 1D Dirac system if one of the mass parameters M z or M x is absent.
Second, E 0D depends on the values of β (−∞) and β (+∞) as seen from Eq. (15). Nevertheless, since Scientific Reports | (2021) 11:21060 | https://doi.org/10.1038/s41598-021-00577-z www.nature.com/scientificreports/ takes place for any values of β (−∞) and β (+∞) , the corner state energy always lies in the band-gap of the edge states as soon as Eq. (26) is fulfilled. One can show that E 0D may formally achieve the energies of the 1D band edges at certain values of β * (−∞) and β * (+∞) , corresponding to The latter however represents the moment, when the corner state becomes delocalized. Finally, we take into account the small terms proportional to v x and v y previously neglected in Eq. (12). The straightforward calculations on the basis of � 0D (x) (see Supplementary Materials) lead to the first-order energy shift: where One can verify that indeed |δE 0D | ≪ |E 0D | for both prototype QWs. The calculations for another Kramer's pair of the edge states described by H * 1D (−k x − k c , θ, ϕ) results in the same energy of the corner state.

Trivial corner states.
Let us now discuss the physical origin of the localized 0D corner state found above.
Since it is clear from Eq. (15), the existence condition of Eq. (23) (or equivalently Eq. (26)) is fulfilled for different functions m x (x) , m y (x) and β (x) . The function β (x) characterizes the corner itself, while m x and m y include the characteristic of the entire system. Further, we show that, in the most general case, arising of 0D corner state simultaneously depends on the corner boundary conditions and the cubic symmetry of the system. For a better understanding of the existence condition of Eq. (26), we first neglect the terms arising due to the cubic symmetry. In this case, one should set parameters R (a) 1,2 and R (a) 1,2 of 2D Hamiltonian in Eq. (5) to zero also resulting in the zero values of F a and F a in Eq. (11). The latter, in turn, leads to the vanishing of m y , v y and δ y in Eq. (10), while m x , v x and δ x become independent of x. Thus, in the absence of the cubic symmetry, by means of Eqs. (15) and (18), the existence condition of Eq. (26) is reduced to where β 1 =β(x = −∞) and β 2 =β(x = +∞) chosen in accordance with Eq. (13). Since it is easy to see, the localized corner state exists if β 1 = β 2 . These parameters can be given a precise physical interpretation.
Let us consider a δ-function electrostatic potential at x = 0 . Then, if we integrate the Schrödinger equation with H 1D (δk, x) + V 0 δ(0)I 2 (where V 0 is a real parameter) in the range of −η≤x≤η (where η is a positive infinitesimal quantity), we directly obtain Eq. (13) . Therefore, β 1 = β 2 can be interpreted as the presence of δ-like electrostatic potential localized at the corner. Thus, the corner states arising in this case are topologically trivial, since they are independent of the cubic symmetry of the system. We also note their independence of the sign of V 0 , as clearly seen from Eq. (30).

Topological corner states.
We now focus on the opposite case, when there is no potential barrier at the corner i.e. β 1 = β 2 . As seen from Eq. (13), these parameters can be both set to zero without loss of generality, which leads to M x (x, 0) = m x (x) and M z (x, 0) = −m y (x) . In this case, the presence of the corner states is governed by the cubic symmetry of 2D system, represented by the non-zero values of F a and F a in Eq. (11). The latter means that the picture of the symmetry-protected corner states should strongly depends not only on crystallographic orientations of the meeting edges but also on the QW growth direction. Figure 4A,B present the calculations for the prototype QWs with double band inversion grown along (001) crystallographic orientation (see the diagrams in Fig. 1). For both QWs, the angles ϕ 1 and ϕ 2 defining the corner orientation in Fig. 3A are chosen to be ϕ 1 = π/3 and ϕ 2 = 2π/3 , while θ in Eq. (10) is assumed to be zero. One can verify that the condition of existence of Eq. (26) is fulfilled for these angles, and, therefore, the corner states arise in the system. As seen from Fig. 4A,B, due to the relationship between m x (ϕ) and m y (ϕ) , the localized energies for internal σ = +1 and external σ = −1 corners are very close to the extrema of 1D dispersion of the edge states. Figure 4C,D represent the angle diagram showing the edge orientations in the (001)-oriented prototype QWs, for which the existence condition of Eq. (26) is valid. It is seen that the values of ϕ 1 and ϕ 2 yielding the corner states are qualitatively represented as the grey regions elongated along the lines defined by excluding the points lying at verticals and horizontals corresponding to ϕ 1 = m 0 π/4 and ϕ 2 = l 0 π/4 , where n 0 , m 0 and l 0 are integers. As seen from Fig. 3A, the 2D system edges at these angles, coincides with [110], [110] , [100] or [010] crystallographic planes, being the faces of (001)-oriented QWs in this case. Since we have neglected the terms resulting from possible breaking of inversion symmetry in the system, the (001)-oriented prototype QWs www.nature.com/scientificreports/ preserve the mirror symmetry about the mentioned planes. Therefore, one concludes that if one of the QW edge coincides with the mirror symmetry planes, the corner states are absent in the system. In light of the above, Eq. (31) describes the corners, whose bisector corresponds to the one of the mirror symmetry planes. The latter results in reflection symmetry of the angle diagram with respect to the lines described by Eq. (31). As expected from the four-fold rotational symmetry of (001)-oriented QWs with respect to the growth direction, the angle pattern in Fig. 4C,D also possess a π/2-periodicity in ϕ 1 and ϕ 2 .
So far, we have discussed the case of (001)-oriented QWs, by whose example we have shown a direct relationship between topological corner states and the symmetry elements of the QW. As the inversion symmetry is preserved in our model, the (001)-oriented QWs have the point group D 4h origin from the point group O h of bulk multi-band k · p Hamiltonian (see Supplementary Materials). Obviously, if one reduces the symmetry of 2D system, the picture of the corner states should change as well. Further, we consider C 2h -symmetric QWs grown along (013) crystallographic direction inspired by recent experimental investigations of double HgTe QWs of the same orientation [58][59][60] . Figure 5A,B provide the angle diagram of topological corner states in the prototypes three-layer InAs/GaInSb and double HgTe/CdHgTe QWs oriented along (013) crystallographic direction. The calculations are performed by assuming θ = arctan(m/n) in Eq. (10) with (m, n) = (1, 3). It is seen that the values of ϕ 1 and ϕ 2 , for which the existence condition is fulfilled, represent significantly different picture of the corner states than the one for (001)-oriented QWs.
We first note π-periodicity in ϕ 1 and ϕ 2 of the angle pattern due to the two-fold rotational symmetry of (0mn)oriented QWs in contrast to the π/2-periodicity of (001)-oriented 2D systems. Second, due to the symmetry lowering, (013)-oriented QWs have only two [100] and [031] mirror planes perpendicular to the QW plane. These planes are represented by dotted verticals and horizontals at ϕ 1 = m 0 π/2 and ϕ 2 = l 0 π/2 in Fig. 5A,B. It is seen that the corner states are absent if one of the QW edge coincides with the mirror symmetry planes. Finally, in contrast to (001)-oriented QWs, (0mn)-oriented QWs possess the corner states if one of the edges is at π/4 ° with respect to the (100) crystallographic axis. The latter is clearly seen in Fig. 5C, which shows the negative values of the function illustrating the existence condition of Eq. (26) for the corner states.

Discussion
To summarize, we have investigated the existence conditions for 0D corner states in cubic semiconductor QWs with double band inversion. We have demonstrated that 0D corner states in such 2D system can appear either due to the presence of electrostatic potential localized in the corner, or due to the crystalline symmetry of the QW. The former case can be physically interpreted as the impurity potential resulting in trivial corner states. The latter corresponds to the symmetry-induced corner states inherent for 2D HOTI state in cubic semiconductor QWs with double band inversion. We have shown that the corners hosting 0D topological states depend not only on crystallographic orientations of the meeting edges but also on the growth orientation and parameters of the QWs.
Let us now make a few remarks concerning the approximation used in this work. As it is clear from the above, theoretical investigations were performed on the basis of multi-band k · p Hamiltonian 49 for the envelope Bloch functions, which is indeed valid on scales much larger than the unit cell. In this case, all functions that are changed on a scale comparable to the unit cell are considered as step-like or delta-functions. In this sense, a meeting point of the corner shown in Fig. 3A represents a 2D unit cell at the corner of an effective 2D system in the QW plane. The special arrangement of atoms in the vicinity of the corner preserving C n symmetry of the bulk may lead to the fractional charge of the corner states, known as a fractional corner anomaly 22,23,[61][62][63] . We note that since fractional corner anomaly requires a theoretical description on the unit cell scale, it cannot be treated on the basis of multi-band k · p Hamiltonian.
Throughout the work, we have neglected the inversion symmetry breaking terms resulting from anisotropy of chemical bonds at the QW interfaces 54 and possible bulk inversion asymmetry of the unit cell 55 . The former is known as an interface inversion asymmetry 54 , while the latter causes the difference between O h and T d point groups of diamond and zinc-blende semiconductors. Taking into account these additional terms in bulk k · p Hamiltonian leads to the non-diagonal blocks of the effective 2D Hamiltonian H 2D (k x , k y ) in Eq. (1). Since the inversion symmetry breaking terms do not affect the very fact of double band inversion, their presence will modify only the edge state parameters v x , v y , m x , m y , δ x , δ y keeping the form of 1D edge Hamiltonian in Eq. (9).
For instance, taking into account the interface and bulk inversion asymmetry in (001)-oriented zinc-blende semiconductor QW reduces the point symmetry from D 4h to C 2v 54 , which, in their turn, also changes the periodicity of v x (ϕ) , v y (ϕ) , m x (ϕ) , m y (ϕ) , δ x (ϕ) , δ y (ϕ) in Eq. (9) from π/2 to π . In this case, the symmetry-induced corner states will arise at the corner with the meeting edges, whose orientations differ from those shown in Fig. 4. For the prototype InAs/GaInSb and HgTe/CdHgTe QWs considered in this work, the terms due to interface and bulk inversion asymmetry are small and induce therefore only a slight modification of the phase diagrams.
We emphasize that since the presence of the corner states depends on the mutual ratio between m x and m y , one may find the QW orientation and strength of inversion-symmetry breaking terms, for which the existence condition of Eq. (26) cannot be fulfilled. The study of all possible cases of violation of Eq. (26) is however out the scope of our work first considering the QWs based on IV, II-VI and III-V semiconductors, in which inversionsymmetry breaking terms are usually small.
Finally, we stress the importance of theoretical results obtained in this work in view of possible applications and their impact on further experimental investigations. After the tremendous interest in II-VI and III-V www.nature.com/scientificreports/ semiconductor QWs induced by prediction and observation of QSHI [39][40][41] , our work shows the importance of cubic semiconductor QWs for the realization of high-order topological states as well. In view of mature growth of IV, II-IV and III-V semiconductor QWs on Si-wafers 64,65 as well as device fabrication technology, our results provide an important first step to future realistic electronics operating on the basis of higher-order topological states including higher-order topological superconductors in hybrid devices 66,67 .

Methods
Band structure calculations were performed by using multi-band k · p Hamiltonian 49 , which directly takes into account the interactions between Ŵ 6 , Ŵ 8 , and Ŵ 7 bands in bulk materials. This model well describes the electronic states in a wide range of narrow-gap semiconductor QWs, particularly in the InAs/GaInSb 46,48 and HgTe/CdHgTe QWs [58][59][60] . In the multi-band k · p Hamiltonian, we also took into account the terms, describing the strain effect arising because of the mismatch of lattice constants in the buffer, QW layers, and barriers. The calculations had been performed by expanding the eight-component envelope wave functions in the basis set of plane waves and by numerical solution of the eigenvalue problem. Details of calculations and the form of the Hamiltonian can be found in the study of Krishtopenko et al. 49 . Parameters for the bulk materials and valence band offsets for the the InAs/GaInSb and HgTe/CdHgTe QWs used in the calculations are provided in Ref. 50 and Ref. 49 , respectively. To derive effective 2D Hamiltonian valid in the vicinity of the Ŵ point from the multi-band k · p Hamiltonian, we implied the procedure proposed by Bernevig et al. 2 and described in details in the Supplementary Materials.