In the last few years, higher-order topological phases have attracted tremendous interest owing to the emergence of novel boundary physics beyond the description of conventional bulk-boundary correspondence [1–27]. Roughly speaking, a mth-order topological phases in n dimensions (n ≥ m) harbors robust gapless states at boundaries with dimensions lower than the bulk by m [28]. For instance, a second-order topological insulator (SOTI) will harbor robust bound states at the sample corners in two dimensions and chiral or helical propagating states on the sample hinges in three dimensions. From a boundary perspective, the higher-order topology can be intuitively understood as the existence of certain forms of Dirac domain walls on the boundary [29, 30]. If the Dirac domain walls are symmetry-enforced, the corresponding higher-order topological phases are commonly categorized to the intrinsic class, otherwise they are categorized to the extrinsic class.

The higher-order topology greatly extends the freedom of bound states on the boundary. On one hand, the positions of bound states may be highly tunable by some external fields if they are not pinned due to symmetry constraints [31, 32]. Such tunability may benefit the design of devices based on higher-order topological materials for applications [33–35]. On the other hand, bound states with different dimensions can coexist in the same system [36], which may result in novel responses absent in conventional topological phases. The higher-order topology also enriches the possibility of topological phase transitions [37–42]. With the advent of higher-order topology, topological phase transitions can happen not only between the different topological sectors of a topological phase following Z or Z ₂ classification [43, 44], but also between phases with different orders in topology [37].

Higher-order topological insulators and semimetals have already been implemented in a great amount of platforms [45–56], while the experimental progress in higher-order topological superconductors is slow as the superconductivity is an intrinsic quantum phenomena without classical counterparts [57]. Thus far, the experimental studies on higher-order topological insulators and semimetals have also been mainly done in metamaterial systems owing to their great advantages in controllability in design and measurements [58]. In comparison, in spite of the accumulation in quantum material candidates for higher-order topological insulator and semimetals [59–72], related experiments on quantum materials remain relatively rare due to the much higher complexity in material growth and detecting the expected signatures of higher-order topology [73–77].

While higher-order topological insulators and semimetals have been intensively studied in experiments in metamaterial systems, we notice that the response of higher-order topological phases to external fields remains poorly explored [78–85]. In this work, we systematically investigate the effects of Zeeman field on time-reversal invariant SOTIs, with a particular focus on the possible topological phase transitions and the associated evolution of boundary states. It is worth noting that while the Zeeman field in a quantum material can be induced by an external magnetic field or the proximity-induced exchange field from magnetic materials in contact, the Zeeman field in metamaterials cannot be implemented in the same way as the underlying spin degrees of freedom are not real spin. Nevertheless, effective Zeeman fields in metamaterials can also be generated by appropriately designing the on-site potentials and couplings between sites [86]. In a quantum material, if the Zeeman field is generated by an external out-of-plane magnetic field, it is known that orbital effect also needs to be taken into account. However, if the Zeeman field is generated by an in-plane magnetic field or exchange interaction in a quantum material, the orbital effect is absent. Moreover, the effective Zeeman field in matematerials is not tied with an orbital effect as the spin is pseudo. With these in mind, in this work we will only focus on the Zeeman field and not consider the orbit effect.

The main findings of our work can be summarized as follows. In two dimensions, when the concerned time-reversal invariant SOTI with each corner harboring two zero-energy bound states is subjected to the Zeeman field, we find that the fate of the zero-energy corner states depends on the direction, the strength, and the type of Zeeman field. When the Zeeman field is a layer-independent ferromagnetic type and in the out-of-plane direction, we find that the two zero-energy bound states per corner are robust in the weak-field regime even though the time-reversal symmetry and chiral symmetry are broken. With the increase of Zeeman field, a bulk topological phase transition will happen at a critical value, and the system will enter a Chern insulator phase [87]. The topological phase transition is accompanied with the disappearance of corner modes and the emergence of chiral edge states. When the ferromagnetic-type Zeeman field is changed to lie in the lattice plane, we find that the number of zero-energy bound states per corner also does not change in the weak-field regime. With the increase of field strength, we find that the system can first undergo a boundary topological phase transition and enter a new SOTI phase with the number and the locations of zero-energy bound states depend on the direction of Zeeman field. Remarkably, we find that corner states can still exist even when the system is driven into a topological semimetal phase with flat-band edge states. When the Zeeman field changes to be an antiferromagnetic type, we find that the number of zero-energy bound states per corner does not change in the weak-field regime if the Zeeman field is along the out-of-plane direction. With the increase of field strength, the system will undergo a bulk topological phase transition and enter a new SOTI with each corner only hosting one zero-energy bound state. When the antiferromagnetic-type Zeeman field lies in the lattice plane, we find that the energies of the corner states will in general be shifted away from zero. Only when the Zeeman field is along certain special direction, the corner states can accidentally be pinned at zero energy. In three dimensions, we focus on the ferromagnetic-type Zeeman field as it is more realistic than the antiferromagnetic one in a bulk system. Similar to the situations in two dimensions, we find that the helical states on the hinges are robust in the weak-field regime. Interestingly, when the field is strong enough to drive the system into a Weyl semimetal phase [88], we find that the resulting Weyl semimetal can host hybrid-order topological properties. That is, first-order topological surface Fermi arcs and second-order topological helical or chiral hinge states can coexist on the boundary.

The paper is organized as follows. In Section 2, we start with a minimal model of time-reversal invariant SOTI and systematically study what kind of topological phase transitions can be induced by Zeeman fields. In Section 3, we generalize the study to three dimensions and investigate what kind of novel topological phases can be achieved. We discuss the results and conclude the paper in Section 4.

As an important consequence of the presence of Zeeman field is the breaking of time-reversal symmetry, in this paper we will start with a time-reversal invariant SOTI so that the effect of Zeeman field can be fully appreciated. While theoretically time-reversal invariant SOTIs allow a lot of model realizations, the essential physics is similar. Without loss of generality, here we consider the following minimal Hamiltonian [4], H k = m − t ⁡ cos k x − t ⁡ cos k y ρ 0 σ z s 0 + λ ⁡ sin k x ρ 0 σ x s x + λ ⁡ sin k y ρ 0 σ x s y + η cos k x − cos k y ρ y σ y s 0 + ρ 0 σ 0 B ⋅ s , (1) where ρ _i , σ _i and s _i are Pauli matrices acting on layer, orbital and spin degrees of freedom, respectively. For notational simplicity, the lattice constants are set to unity throughout this paper. If the Hamiltonian only has the first three terms, the Hamiltonian simply describes two layers of decoupled topological insulators when 0 < |m| < 2t and λ ≠ 0 [89]. Accordingly, there are two pairs of helical edge states on the boundary. Since time-reversal invariant insulators in two dimensions follow a Z ₂ classification if there is no additional crystalline symmetry [90], the two pairs of helical edge states are not topologically protected and can be gapped out by perturbations respecting time-reversal symmetry. The fourth term couples the two layers and leads to the hybridization of the edge states to open a gap. Interestingly, because of the existence of line nodes along the directions k _x = ±k _y , this term results in a time-reversal invariant SOTI with doubly-degenerate zero-energy bound states appearing at the intersecting corners of x-normal and y-normal edges [4]. The last term represents the Zeeman field which breaks time-reversal symmetry. Here we first consider that the two layers have the same Zeeman field, corresponding to a ferromagnetic type.

When the Zeeman field is along the z direction, i.e., B = (B _x , B _y , B _z ) = (0, 0, B _z ), besides the time-reversal symmetry, one can find that the chiral symmetry is also broken. As a result, one should expect that the doubly-degenerate bound states per corner of the time-reversal invariant SOTI are not stable against the Zeeman field. However, counterintuitively, this is not the case. To see this, the most evident way is to derive the boundary Hamiltonian explicitly since the corner states are a result of the formation of Dirac domain walls on the boundary [12]. Without loss of generality, let us first derive the low-energy Hamiltonian describing the left x-normal edge. Assuming that the minimum of the band gap is located at the time-reversal invariant momentum Γ = (0, 0), we expand the lattice Hamiltonian around Γ to obtain the low-energy bulk Hamiltonian. Keeping momentum up to the second order, it takes the form H k = m ̃ + t 2 k x 2 + k y 2 ρ 0 σ z s 0 + λ k x ρ 0 σ x s x + λ k y ρ 0 σ x s y − η 2 k x 2 − k y 2 ρ y σ y s 0 + B z ρ 0 σ 0 s z . (2) where m ̃ = m − 2 t is assumed to be negative so that the Hamiltonian is in the SOTI phase when B _z = 0. For the convenience of discussion, below we consider that coefficients {t, λ, η, B _z } are all positive. When focusing on the left x-normal edge, we can assume that the system occupies the half-infinity plane 0 ≤ x < + ∞. Because of the breaking of translation symmetry along the x direction, the momentum k _x needs to be replaced by − i∂ _x . Next we divide the Hamiltonian into two parts, H = H ₁ + H ₂, where H 1 − i ∂ x , k y = m ̃ + t 2 k y 2 − t 2 ∂ x 2 ρ 0 σ z s 0 − i λ ∂ x ρ 0 σ x s x , H 2 − i ∂ x , k y = λ k y ρ 0 σ x s y − η 2 k x 2 − k y 2 ρ y σ y s 0 + B z ρ 0 σ 0 s z . (3)

By solving the eigenvalue equation H ₁ (−i∂ _x , k _y )ψ(x) = Eψ(x) under the boundary conditions ψ(0) = ψ(+∞) = 0, one finds the existence of four zero-energy bound states as long as m ̃ + t k y 2 / 2 < 0 . The wave functions for the four zero-energy bound states can be compactly written in the form [12]. ψ ρ σ x = N ⁡ sin κ 1 x e − κ 2 x e i k y y χ ρ σ , (4) where χ ρ σ = | ρ z = ρ , σ y = σ , s x = σ ̄ 〉 with ρ = ±1, σ = ±1 and σ ̄ = − σ , the normalization constant N = 2 κ 2 ( κ 1 2 + κ 2 2 ) / κ 1 2 with κ 1 = − 2 m ̃ / t − k y 2 − ( λ / t ) 2 and κ ₂ = λ/t. The effect of H ₂ on these zero-energy bound states can be determined by projecting H ₂ onto the four-dimensional subspace spanned by the four orthogonal wave functions of the zero-energy bound states. Since the bound states from H ₁ have zero energy, this projection directly determines the low-energy boundary Hamiltonian. Put it explicitly, the low-energy Hamiltonian is given by H x L ; ρ σ , ρ ′ σ ′ k y = ∫ 0 + ∞ ψ ρ σ † x H 2 − i ∂ x , k y ψ ρ ′ σ ′ x d x . (5)

Simple algebra calculations reveal that under the basis ( χ 11 , χ 1 − 1 , χ − 11 , χ − 1 − 1 ) T , the low-energy boundary Hamiltonian can be expressed in terms of the Pauli matrices as H x L k y = − λ k y ρ 0 σ x − η m ̃ t ρ y σ z . (6)

The second term on the right hand side corresponds to the Dirac mass induced by the η term. Remarkably, one can find that the Zeeman field does not enter the boundary Hamiltonian. In other words, the Zeeman field does not have a direct effect on the boundary Hamiltonian, so the corner modes, even though both time-reversal symmetry and chiral symmetry are broken. Using the lattice Hamiltonian, we have numerically verified that the doubly-degenerate zero-energy bound states per corner are indeed robust against the Zeeman field in the weak-field regime, as shown in Figure 1.

Although the out-of-plane Zeeman field does not directly enter the boundary Hamiltonian, it affects the boundary through the bulk indirectly. To be more specific, the Zeeman field will affect the boundary states by reducing the bulk energy gap. To see this, let us diagonalize the Hamiltonian and obtain the energy spectrum, which reads E k = ± M k ± B z 2 + λ 2 sin 2 k x + sin 2 k y , (7) where M ( k ) = m 2 ( k ) + η 2 ( k ) with m( k ) = (m − t cos k _x − t cos k _y ) and η( k ) = η(cos k _x − cos k _y ) is introduced to shorten the expression. Apparently, the Zeeman field will directly change the bulk energy gap. In fact, it is easy to find that the bulk energy gap will vanish at one or two of the four time-reversal invariant momentums, Γ = (0, 0), X = (π, 0), Y = (0, π), M = (π, π), when B _z takes any one of the follow-up values, ± (m ± 2t) and ± m 2 + 4 η 2 . To see what kind of new topological phases can be induced by the out-of-plane Zeeman field, let us assume that m is close to 2t and η is much smaller than t, then the bulk energy gap undergoes the first closing-and-reopening transition when B _z is increased to the value B _c,Γ ≡|m − 2t|, and the last transition when B _z is increased to the value B _c,M ≡|m + 2t|. Let us first consider that B _z is approaching to B _c,Γ. Near this critical value, the band-gap minimum is located at Γ. To simplify the follow-up discussions, we keep momentum only up to the linear order in the low-energy bulk Hamiltonian. Accordingly, we have H k = m ̃ ρ 0 σ z s 0 + λ k x ρ 0 σ x s x + λ k y ρ 0 σ x s y + B z ρ 0 σ 0 s z = H ρ z = 1 k ⊕ H ρ z = − 1 k . (8)

Since H ρ z = ± 1 are the same, we can focus on just one of them, e.g., H ρ z = 1 . Interestingly, H ρ z = 1 itself takes a block diagonal form, i.e., H ρ z = 1 = H + ⊕ H − with H + k = m ̃ + B z σ z + λ k x σ x + λ k y σ y , H − k = m ̃ − B z σ z + λ k x σ x − λ k y σ y . (9)

It is worth noting that although the mathematical forms of Pauli matrices are same, the physical interpretation of the Pauli matrices σ _i in H _± is distinct to that in Eq. 1 as their underlying bases are different. However, here we are only interested in the topological property hiding in the mathematical structure of the Hamiltonian. When B _z is approaching to B _c,Γ, the energy gap for H ₊( k ) will undergo a closing-and-reopening transition while the one for H ₋( k ) remains open. Therefore, only H ₊( k ) will undergo a topological phase transition. Because neither time-reversal symmetry nor chiral symmetry is conserved, the first-order topology of H ₊( k ) is characterized by the first-class Chern number [91]. For a two-band model, the Chern number can be simply determined by the formula [92]. C = − 1 4 π ∫ d ⋅ ∂ d ∂ k x × ∂ d ∂ k y d 3 d 2 ⁡ k , (10) where d = (d _x , d _y , d _z ) with d _i being the coefficient before the Pauli matrix σ _i and the integral is over the whole two-dimensional momentum space. A simple calculation reveals that the Chern number for H ₊( k ) is C ₊ = sgn (−B _c,Γ + B _z )/2. It indicates that the Chern number has a jump ΔC = 1 when B _z is changed from B _z < B _c,Γ to B _z > B _c,Γ. As the Chern number is forced to be zero by time-reversal symmetry when B _z = 0 [93], the change of Chern number indicates that the system enters a Chern insulator phase with the total Chern number C = 2 after taking into account the contributions from two layers. In other words, an out-of-plane Zeeman field can induce a topological phase transition from a SOTI to a Chern insulator with the increase of its strength. We have numerically demonstrated this topological phase transition, as shown in Figure 2. After the topological phase transition, the zero-dimensional corner states disappear and one-dimensional chiral edge states emerge on the boundary. It is known that higher-order topology can naturally descend from symmetry-protected first-order topology by appropriately lifting the protecting symmetry [7, 8]. The above result indicates that the reverse is also true. As a final remark of this part, we point out that increasing the Zeeman field will drive the system from the Chern insulator phase back to the trivial insulator phase. This is easy to see since the system will go to the topologically trivial atomic limit when the strength of Zeeman field goes to infinity.

Let us now turn our attention to the case with the Zeeman field lying in the xy plane, i.e., B = (B _x , B _y , 0). Unlike the out-of-plane Zeeman field, the in-plane Zeeman field does not break the chiral symmetry. It is worth noting that due to the simplicity of the Hamiltonian, there are two chiral symmetry operators, ρ _y σ _x s _z and ρ ₀ σ _x s _z . Since the time-reversal symmetry is broken and the chiral symmetry is preserved, the Hamiltonian falls into the symmetry class AIII. In two dimensions, this symmetry class does not support nontrivial first-order topological insulator phases [43, 44].

To see what kind of new topological phases can be induced by the in-plane Zeeman field, we follow previous procedures and firstly study the effect of Zeeman field on the boundary Hamiltonian. Let us still focus on the left x-normal edge for illustration. Similarly, the contribution from the Zeeman field to the boundary Hamiltonian is determined by H Z ρ σ , ρ ′ σ ′ = ∫ 0 + ∞ ψ ρ σ † x ρ 0 σ 0 B x s x + B y s y ψ ρ ′ σ ′ x d x = − B x ρ 0 σ z ρ σ , ρ ′ σ ′ . (11)

Adding this contribution, the boundary Hamiltonian on the left x-normal edge becomes H x L k y = − λ k y ρ 0 σ x − η m ̃ t ρ y σ z − B x ρ 0 σ z . (12)

Since [ ρ y , H x L ] = 0 , the Hamiltonian can be decomposed as the direct sum of two parts, H x L = H ρ y = 1 ⊕ H ρ y = − 1 , where H ρ y = ± 1 k y = − λ k y σ x − B x ± η m ̃ t σ z . (13)

Apparently, the x-component of the in-plane Zeeman field directly changes the boundary Dirac mass on the x-normal edges. Similar analysis can also find that the y-component of the in-plane Zeeman field directly changes the boundary Dirac mass on the y-normal edges in the same way. According to Eq. 13, it is easy to find that the Dirac mass of one of the branches vanishes when B x = ± η m ̃ / t . This means that increasing the field strength from zero will induce a closing-and-reopening transition in the boundary energy gap, leading to the occurrence of a boundary topological phase transition [94]. Since the Zeeman field neither breaks the chiral symmetry nor changes the sign of Dirac mass when its value is smaller than the critical value ( | η m ̃ / t | ) , the doubly-degenerate corner states are also expected to be robust in the weak-field regime as long as the system size is sufficiently large so that the wave functions of corner states remain well separated in real space. When the Zeeman field is strong enough to induce a closing-and-reopening transition in the boundary energy gap, the sign change of Dirac mass will change the topological property of the domain walls at the corners. Accordingly, the number of corner states will change and the locations of corner states will depend on the direction of Zeeman field.

To show the evolution of boundary physics with respect to Zeeman field more explicitly, we numerically diagonalize the Hamiltonian under different boundary conditions. In Figures 3A–C, the energy spectra correspond to a cylindrical geometry with periodic boundary condition in the x direction and open boundary condition in the y direction. It is readily seen that the boundary energy gap undergoes a closing-and-reopening transition with the increase of Zeeman field, in agreement with the above analytical analysis. In Figures 3D–F, the number and locations of corner states are shown. The parameters in Figures 3D–F are respectively the same as Figures 3A–C, what is different is that open boundary conditions are adopted in both the x and y directions. It is readily found that the eight zero-energy corner states do remain robust in the weak-field regime. After the boundary topological phase transition, the number of zero-energy corner states is reduced by half. Interestingly, the remaining four zero-energy bound states are found to appear only at two of the four corners, with each one of the two corners hosting two zero-energy bound states. In Figures 3C,F, as B _x = − B _y is considered, the result indicates that the zero-energy bound states appear at the two corners whose connecting line is perpendicular to the direction of Zeeman field, implying that the positions of corner states can be manipulated by controlling the direction of Zeeman field.

Since the Dirac mass on the x(y)-normal edge is affected by the x(y)-component of the in-plane Zeeman field, the evolutions of the corner modes can be briefly described as follows. For the convenience of discussion, we also consider B _x and B _y to be positive and define M η = | η m ̃ / t | . When B _x < M _η and B _y < M _η , each corner of the square-lattice sample harbors two zero-energy bound states as the Zeeman field is not strong enough to induce a boundary topological phase transition. When B _x < M _η < B _y or B _y < M _η < B _x , a sign change in Dirac mass has occurred on the y-normal or x-normal edges. As a result, the number of Dirac domain walls per corner will reduce from two to one, resulting in a new SOTI with one zero-energy bound state per corner. When M _η < B _x and M _η < B _y , the sign change in Dirac mass has occurred on both x-normal and y-normal edges. As we have illustrated in Figure 3F, for this case, two of the four corners harbor zero-energy bound states, and the positions of corner states are determined by the direction of Zeeman field. In Figure 4, the number and locations of zero-energy bound states are shown for the case with B _y < M _η < B _x . One can find that each corner does harbor one zero-energy bound state, in agreement with the above analysis.

For the in-plane Zeeman field, so far we have only investigated its effect on the boundary. When the in-plane Zeeman field is strong enough, we find that the system will be driven into a topological semimetal phase. To see this, let us consider that the Zeeman field is in the x direction. Accordingly, the bulk energy spectrum is E k = ± M 2 k + λ 2 ⁡ sin 2 k x ± B x 2 + λ 2 ⁡ sin 2 k y . (14)

If we assume that the value of m is close to 2t so that the minimum of band gap is located at Γ, then the system becomes a topological semimetal when B _x becomes larger than B _c,Γ = |m − 2t|, with two nodal points appearing on the k _y = 0 axis. As the chiral symmetry is conserved, the topological properties of these nodal points are characterized by the winding number defined along a closed contour enclosing the nodal point [95]. Due to this topological protection, the nodal points are robust against perturbations respecting chiral symmetry and can only be annihilated pairwise.

If one views the system as a one-dimensional system by taking one of the momentum as a parameter, then the nodal points can also be viewed as critical points across which the winding number has a discrete change. Since a one-dimensional insulator with nontrivial winding number harbors zero-energy bound states on its two ends, this implies that the two-dimensional topological semimetal harbors flat-band edge states, as shown in Figure 5A. Remarkably, we find that even though the system becomes a semimetal, sharply-localized corner states can still appear if the Zeeman field is applied along the diagonal direction of the square geometric sample, as shown in Figures 5B–D.

Above the Zeeman field has been assumed to be same in the two layers, corresponding to a ferromagnetic type, in the following we extend the study to the antiferromagnetic case for which the Zeeman fields in the two layers are opposite. Accordingly, the Hamiltonian takes the form H k = m − t ⁡ cos k x − t ⁡ cos k y ρ 0 σ z s 0 + λ ⁡ sin k x ρ 0 σ x s x + λ ⁡ sin k y ρ 0 σ x s y + η cos k x − cos k y ρ y σ y s 0 + ρ z σ 0 B ⋅ s . (15)

Similarly, let us first focus on the out-of-plane case. It is easy to find that the chiral symmetry remains for this case, with the remaining chiral symmetry operator being ρ _y σ _x s _z . According to our previous discussion, the existence of chiral symmetry implies that the antiferromagnetic Zeeman field cannot drive the Hamiltonian into a Chern insulator phase. As the chiral symmetry still remains, the doubly-degenerate corner states are robust in the weak-field regime. By projecting the Zeeman field onto the subspace spanned by boundary states, one can find that it still does not enter the Hamiltonian. In other words, the out-of-plane antiferromagnetic Zeeman field cannot induce boundary topological phase transitions. Nevertheless, the topological property of the Hamiltonian can be changed with the increase of Zeeman field. This can be simply inferred from the bulk energy spectrum which reads E k = ± m k ± B z 2 + η 2 k + λ 2 sin 2 k x + sin 2 k y . (16) With the increase of Zeeman field, the bulk energy gap will undergo closing-and-reopening transitions. To be specific, let us again assume the value of m to be close to 2t so that the minimum of band energy gap is located at Γ. Accordingly, with the increase of Zeeman field, the first transition occurs at B _z = B _c,Γ = |m − 2t|. When B _z > B _c,Γ, the Hamiltonian enters a new SOTI phase with each corner harboring one zero-energy bound state, as shown in Figure 6.

When the antiferromagnetic-type Zeeman field lies in the xy plane, similar to the ferromagnetic-type case, the chiral symmetry also remains. However, now the Zeeman field will directly enter the boundary Hamiltonian. For instance, using the same derivation, we find that the low-energy boundary Hamiltonian for the left x-normal edge becomes H x L k y = − λ k y ρ 0 σ x − η m ̃ t ρ y σ z − B x ρ z σ z . (17)

A remarkable feature of the above Hamiltonian is that the Dirac masses induced by the η term and the Zeeman field become anticommutative. As the presence of an additional anticommutative term induces an overall mass on the boundary, the corner states are thus expected to be split away from zero energy. Through numerical calculations, we find that in general this is indeed the case, as shown in Figures 7A,C. However, when the Dirac masses induced by the Zeeman field and the η term are fine-tuned to have the same ratio on two near-neighboring edges, then their intersecting corner will harbor zero-energy bound states, as shown in Figures 7B,D. This result can be simply understood by noting that when the two Dirac masses with different origins have the same ratio on two near-neighboring edges, in fact we can redefine the Pauli matrices while maintaining their anticommutative relationship. By such a procedure, the boundary Hamiltonian reduces to the standard Dirac-Hamiltonian form with only one mass term which is known to support zero-energy bound states when the mass term has a domain-wall configuration in real space [96].

Since without fine tuning the zero-energy bound states are immediately gapped by the antiferromagnetic-type in-plane Zeeman field, the evolution of boundary physics is less interesting compared to the case with ferromagnetic-type in-plane Zeeman field. In addition, we find that the antiferromagnetic-type in-plane Zeeman field cannot drive the system into a topological semimetal phase. The increase of the field strength only drives the system to another trivial insulator phase, therefore, from the perspective of topology, the bulk physics is also less interesting compared to the case with ferromagnetic-type in-plane Zeeman fields.

As the increase of dimension enriches the allowed types of boundary states, it is natural to expect that the Zeeman field can induce a diversity of topological phases transitions and exotic patterns of boundary states in a three-dimensional SOTI. To generalize the study, we also follow the minimal-model approach as adopted in two dimensions. To be specific, we consider that the minimal Hamiltonian takes the form H k = m − t ∑ i = x , y cos k i − t z ⁡ cos k z ρ 0 σ z s 0 + λ ∑ i = x , y , z sin k i ρ 0 σ x s i + η cos k x − cos k y ρ y σ y s 0 + ρ 0 σ 0 B ⋅ s . (18)

In this part, we will focus on the ferromagnetic-type Zeeman field as the study in two dimensions tells us that it can induce more interesting physics than the antiferromagnetic type. For this Hamiltonian, it realizes a time-reversal invariant SOTI when the last term representing the time-reversal symmetry breaking Zeeman field is absent and the parameters satisfy λ ≠ 0 and η ≠ 0, and |2t − t _z | < |m| < 2t + t _z . If open boundary conditions are taken in both the x and y directions, the second-order topology is manifested through the appearance of a pair of helical states on each of the z-direction hinges [3]. It is worth noting that when the open boundary condition is only applied in the z direction, the top and bottom surfaces still harbor gapless Dirac cones since the η term cannot open a gap at the time-reversal invariant momentum at which the Dirac cones are located.

In parallel to the two-dimensional case, we first consider that the Zeeman field is along the z direction. In the weak-field regime, the main effect of the Zeeman field is to open a gap on the top and bottom surfaces. Following the same analysis as done in two dimensions, it is easy to figure out that the Zeeman field has little impact on the helical states on the z-direction hinges even though the time-reversal symmetry is broken. To be more explicit, note that the three-dimensional Hamiltonian can be viewed as the stacking of an infinite number of two-dimensional layers in momentum space. In particular, the Hamiltonians at the two high-symmetry planes with k _z = 0 and π just reduce to the form of the previously studied two-dimensional Hamiltonian. According to the results in two dimensions, an out-of-plane Zeeman field does not directly affect the corner states in the weak-field regime, for the same reason, the helical states propagating on the z-direction hinges will remain stable in the weak-field regime. Because of the opening of a gap on the top and bottom surfaces and the stability of the helical states on the z-direction hinges against the z-direction Zeeman field, the helical states will no longer flow in the whole surface when they flow towards the top and bottom surfaces as the low-energy channels on these two surfaces have been emptied by the Zeeman field, instead they will be confined to flow along the x-direction or y-direction hinges so that the low-energy transport channels can form closed paths along the hinges.

With the increase of Zeeman field, we have shown that in two dimensions a Chern insulator with chiral edge states will be realized after the bulk energy gap undergoes a closing-and-reopening transition. For the Hamiltonian in Eq. 18, because of the increase in dimension, we find that the bulk energy gap does not immediately reopen when it is reduced to zero by the Zeeman field. Instead, the system enters a Weyl semimetal phase for a broad range of parameters, with the Weyl points located at some high-symmetry k _z axes. It is known that a hallmark of the Weyl semimetal is the Fermi arcs connecting the projections of Weyl points in the surface Brillouin zone [88]. Furthermore, if one treats the momentum k _z as a parameter, the Weyl points correspond to critical points separating regions with C (k _z ) = 0 from the ones with C (k _z ) ≠ 0. Interestingly, this raises the possibility of the coexistence of helical hinge states characterizing second-order topology and surface Fermi arcs characterizing first-order topology. In other words, hybrid-order topology is possible in this system. To see this, without loss of generality, let us assume that H (k _x , k _y , k _z = 0) realizes a SOTI with zero-energy corner states and H (k _x , k _y , k _z = π) realizes a trivial insulator. If the energy gap of H (k _x , k _y , k _z = π) is smaller than that of H (k _x , k _y , k _z = 0), then H (k _x , k _y , k _z = π) can first underdo a topological phase transition and enter a Chern insulator phase. Accordingly, it is clear that there exists a window within which H (k _x , k _y , k _z = 0) remains to be a SOTI due to its stability against the Zeeman field and H (k _x , k _y , k _z = π) becomes a Chern insulator with chiral edge states. In Figure 8, we have numerically demonstrated the above picture. It is worth noting that this class of hybrid-order Weyl semimetal phase is distinct to those so-called higher-order Weyl semimetals in which the hinge states are neither helical nor chiral [97, 98].

When the Zeeman field becomes aligned in the xy plane, according to the study in two dimensions one can simply infer that the helical hinge states are stable against the Zeeman field when its strength is smaller than a critical value, and they will transform into chiral hinge states when the strength of the Zeeman field goes beyond the critical value at which the energy gaps on the x-normal or y-normal surfaces get closed. In Figures 9A–C, we have numerically demonstrated this expected evolution. Similar findings have also been obtained based on low-energy boundary theories in Ref. [78].

Apparently, because chiral hinge states can exist in the strong-field regime, Weyl semimetal phases with hybrid-order topology can also be realized when the parameters satisfy appropriate conditions. For an in-plane Zeeman field, according to the energy spectrum of the Hamiltonian one can infer that the Weyl points are located at the k _z = 0 or π plane when the Zeeman field is strong enough to drive the system into a gapless phase, as shown in Figure 9D. Accordingly, Fermi arcs will appear on the top and bottom z-normal surfaces, indicating a different pattern of the distribution of the boundary states compared to the case with a z-direction Zeeman field.

For first-order topological phases protected by time-reversal symmetry, the one-dimensional lower boundary states are generally gapped by the Zeeman field due to the breaking of time-reversal symmetry [7, 8]. However, this picture cannot be simply generalized to time-reversal invariant higher-order topological phases. The underlying reason is that the gapless bound states in a higher-order topological phases live in a subspace of the total Hilbert space of the bulk Hamiltonian. While the Zeeman field directly breaks the time-reversal symmetry of the bulk Hamiltonian, in general it cannot induce direct couplings between different sectors (e.g., for two corner states, two sectors with spin-up polarization and spin-down polarization) of the subspace [99], so the gapless bound states remain stable. Nevertheless, the Zeeman field will change the localization length of the wave functions for bound states either through the reduction in boundary energy gap or through the reduction in bulk energy gap. When the boundary or bulk energy gap vanishes, the overlap of wave functions will lead to the annihilation of some or all of the bound states, corresponding to the occurrence of boundary or bulk topological phase transitions.

In this work, we have systematically investigated the possible topological phase transitions and the associated evolution of boundary bound states when a time-reversal invariant SOTI is subjected to a Zeeman field. In both two and three dimensions, we find that the bound states characterizing the second-order topology are generally stable in the weak-field regime. With the increase in field strength, both bulk and boundary topological phase transitions can be induced. Depending on the direction, the strength, as well as the type of Zeeman field, the boundary bound states will change their total number or the type across the topological phase transitions. Remarkably, we find that corner states can exist even when the system becomes a topological semimetal in two dimensions. In addition, in three dimensions we find that the Zeeman field can induce hybrid-order Weyl semimetal phases in which surface Fermi arcs and gapless hinge states coexist. Our study reveals that the Zeeman field can induce very rich physics in higher-order topological materials, and the results are applied to both quantum materials and metamaterials.