Hinge states in a system of coupled Rashba layers

We consider a system of stacked tunnel-coupled two-dimensional electron- and hole-gas layers with Rashba spin-orbit interactions subjected to a staggered Zeeman field. The interplay of different intra-layer tunnel couplings results in a phase transition to a topological insulator phase in three dimensions hosting gapless surface states. The staggered Zeeman field further enriches the topological phase diagram by generating a second-order topological insulator phase hosting gapless hinge states. The emergence of the topological phases is proven analytically in the regime of small Zeeman field and confirmed by numerical simulations in the non-perturbative region of the phase diagram. The topological phases are stable against external perturbations and disorder.


I. Introduction
Over the last decade, topological insulators and superconductors (TIs) have attracted a lot of attention in the domain of condensed matter physics [1][2][3][4][5]. One of the main reasons for such popularity resides in the bulk-boundary correspondence relating the topological description of a d-dimensional insulating bulk to the presence of gapless modes on a (d − 1)-dimensional boundary of TIs. Moreover, the topological nature of the boundary modes makes them insensitive to a wide spectrum of external perturbations and disorder -the property of greatest importance in the area of quantum information and quantum computing [6][7][8]. Recently, the concept of bulk-boundary correspondence has been generalized by introducing a new class of topological systems, called higher-order TIs [9][10][11]. In contrast to conventional TIs, the (d − 1)-dimensional boundary of an n-th order TI is insulating, similarly to the bulk. Instead, it exhibits protected gapless modes on the (d − n)-dimensional boundaries. The corresponding gapless modes are called corner states in the case d ≥ 2, n = d or hinge states for d ≥ 3, n = d − 1.
Motivated by this context, we propose here a mesoscopic setup that can be controllably brought into the second-order TI (SOTI) phase in three dimensions.
The key component of our setup is an array of twodimensional electron gas (2DEG) layers with Rashba spin-orbit interaction (SOI) and electron/hole-like dispersions [38,39] (see Fig. 1). Different Rashba layers are tunnel coupled to each other, and the whole heterostructure is subjected to a staggered Zeeman field. The phase diagram of the resulting model comprises a strong threedimensional topological insulator (3DTI) phase at zero Zeeman field, which transforms into the SOTI phase when the Zeeman field is switched on. In this work we do not include the effect of interactions, but we point out the possibility to create even more exotic fractional higherorder topological states in a version of our system that contains strong electron-electron interactions [38,39].
This work is organized as follows. First, in Sec. II we introduce the model for coupled Rashba layers. Second, in Sec. III we briefly present the phase diagram of the model in the regime of a vanishing staggered Zeeman field. We show that the system exhibits a phase transition from the topologically trivial phase to the 3DTI phase hosting gapless surface states. Next, in Sec. IV we analyze the effect of a staggered Zeeman field. We show that it gaps out the surface states and leads to the emer-arXiv:1910.01655v1 [cond-mat.mes-hall] 3 Oct 2019 gence of the SOTI phase hosting gapless hinge states. We also calculate the phase diagram of the model and study the stability of the topological phases against the disorder. Finally, in Sec. V we propose two modifications of the setup, which could facilitate the experimental realization. In the first proposal, the staggered Zeeman field is generated only on the surface of the heterostructure, using an antiferromagnet with a bicollinear antiferromagnetic order or via magnetic impurities. In the second proposal, the pairs of the 2DEGs are replaced by thin TI films.

II. Model
We start by introducing a system of coupled Rashba 2DEG layers [38,39] placed in the xy plane (see Fig. 1). We assume that the corresponding SOI vector is along the z axis and defines the spin quantization axis. The unit cell is composed of two pairs of layers with electron-and hole-like dispersions described by opposite signs of the mass m and labeled by an index η ∈ {1,1}. The layers of the same dispersion type are distinguished by the value of the SOI strength, which have the same amplitude α but opposite signs indexed by τ = {1,1}. Creation operators ψ † ητ σk act on electrons with spin component σ along the z axis and momentum k = (k x , k y , k z ) ≡ (k , k z ) in the layer (η, τ ). Using this notation, the Hamiltonian describing the uncoupled layers reads where the Pauli matrices σ j act in spin space. The chemical potential µ is tuned to the crossing point of the spinsplit bands at k = 0. In the following, it will be convenient to introduce the SOI energy E so = 2 k 2 so /(2m) and the SOI momentum k so = mα/ 2 .
The nearest-neighbor Rashba layers are coupled to each other via spin-conserving tunneling processes. We distinguish two different types of tunnelings: between the layers with the same mass with the amplitude t 1 and between the layers with opposite masses with the amplitude t 2 . The corresponding processes are described by the following terms in the Hamiltonian where a denotes the size of the unit cell in z direction. Finally, the system is subjected to an effective Zeeman term Here, t Z ηn j = µ B g j B j , with B -the magnetic field vector, n -the direction of the Zeeman field, and g j -the g-factor along a given axis. We assume that the sign of the Zeeman splitting alternates between the Rashba layers with electron-and hole-like dispersions. Such a staggered effective Zeeman field can be achieved in several ways. The first option is to use a uniform magnetic field and 2DEGs with opposite signs of the g-factor in layers with different dispersions. A second possibility consists of using magnetic impurities which order ferromagnetically within the same Rashba layers but have a staggered magnetization direction in different layers [40][41][42][43][44]. Yet another option is to use thin ferromagnetically ordered layers, contained between the 2DEGs with the same dispersions. Alternatively, the effective Zeeman field can be applied only to the surface of the system, for example by placing it in the vicinity of a magnetic material, or by doping the surface with magnetic adatoms. As we will show later in Sec. V, these constructions allow us to reach the desired topological phases without affecting the bulk of the system.
After combining all previously introduced ingredients, the total Hamiltonian can be conveniently written in the basis Ψ k = (ψ 11↑k , ψ 11↓k , ψ 11↑k , ψ 11↓k , ψ1 1↑k , ψ1 1↓k , ψ11 ↑k , ψ11 ↓k ) T as where the Pauli matrices η j and τ j act in the space of four layers (η, τ ). In the following, we will present the solution of the problem described by the Hamiltonian H.

III. Strong 3D topological insulator
In order to describe the properties of our model, we first focus on the regime of vanishing Zeeman field (t Z = 0). In this section we calculate the topological phase diagram and deduce the effective degrees of freedom emerging at low energies, associated with the topological surface states.

A. Topological phase diagram
When the Zeeman field is absent, the system is described by the Hamiltonian density with the bulk spectrum given by where = k 2 /(2m). The spectrum does not depend on the direction of the vector k as a consequence of the rotational symmetry in the xy plane. We also see that at zero momentum, k = 0, the two tunneling processes compete with each other, leading to the closing of the gap at the critical point t 1 = t 2 . As a result, the phase diagram is divided into two regimes: t 1 < t 2 and t 1 > t 2 .
The Hamiltonian density satisfies the relation T H TI (k)T −1 = H TI (−k), implying that the system is time-reversal symmetric. Here the time-reversal symmetry operator is expressed as T = iσ y K, with K being the complex conjugation operator. However, we note that the particle-hole symmetry and the chiral symmetry are absent. Hence, the system belongs to the AII symmetry class characterized by a topological invariant ν 0 ∈ Z 2 . The bulk-boundary correspondence relates this topological invariant to the presence of gapless states Φ s ± localized at every surface s of the system and decaying exponentially into the bulk. We verify numerically the presence of such gapless surface states in the regime t 1 < t 2 of our model by diagonalizing the discretized version of the Hamiltonian H TI density, as shown in Figs. 2(a)-(b). This allows us to identify the regime t 1 < t 2 (t 1 > t 2 ) with the topologically trivial (topological) phase, where the topological phase corresponds to a strong 3DTI [38,39].
We study the surface states more closely by considering several non-equivalent geometries of the system with open boundary conditions (OBCs) along one direction and periodic boundary conditions (PBCs) along two other directions. In particular, we study the xy surface states in a geometry with OBC along z and PBCs along x and y [see Fig. 2(a)], and the xz surface in a geometry with OBC along y and PBCs along x and z [see Fig. 2 We observe that the surface states have a linear dispersion relation with one Dirac cone per surface. The shape of the Dirac cone and the spin texture of the surface states change depending on the surface considered. As a consequence of the rotational symmetry in the xy plane, the Dirac cone on the xy surface is isotropic, and the spin texture is anti-vortex like [see

B. Gapless surface states
In order to acquire a better theoretical understanding of the 3DTI phase, we calculate the wavefunctions of the states Φ s ± situated at different surfaces s. We solve the Schroedinger equation obtained by making the substitution k j → −i ∂ j with j ∈ {x, y, z}. For clarity, we only focus on the solutions at the Dirac point, which are obtained by taking two of three components of the momentum k to zero. The effect of the neglected terms is treated perturbatively, based on the solutions at the Dirac point.

Top and bottom surfaces
First we focus on the xy surface by taking k = 0. We rewrite the problem in a geometry with OBC along the z axis by substituting k z → −i ∂ z . We assume that the linear size L z of the system in z direction is sufficiently large, such that the surface states at opposite surfaces can be treated as independent. The resulting simplified real space Hamiltonian density reads Here, we only consider the terms linear in momentum k z . We note that, in addition to a time-reversal symmetry T already present in the full problem, the simplified Hamiltonian has a particle-hole symmetry represented by an operator P = −τ z K satisfying P H xy TI P −1 = −H xy TI and [P, T ] = 0. Moreover, the chiral symmetry described by an operator C = T P = −iτ z σ y is also present in the system. Enforced by these symmetries, the Dirac point of the surface states is located exactly at zero energy.
Solving the Schroedinger equation for H xy TI results in two zero-energy solutions localized at the top surface at z = 0: where ξ = t 1 t 2 a/(t 2 1 − t 2 2 ) and N is a normalization constant. We remind that we work in the basis Φ xy More generally, any combination of Φ xy + and Φ xy − (with its time-reversal partner) is also a solution of the problem, however, only these two states respect the rotational symmetry in the xy plane present in H TI .
The dispersion relation and the effective Hamiltonian at finite momentum k is obtained perturbatively by calculating the expectation values of the remaining terms in H TI in the basis of the surface states Φ xy ± . Here we only describe the low-energy physics close to the Dirac point, assuming that the amplitude of the momentum k is small. Hence, we only take into account the terms linear in k x and k y . We also distinguish between the diagonal and off-diagonal components of the expectation values. We find that all the diagonal components vanish exactly, i.e.
However, the off-diagonal components are non-zero: This allows us to write the effective low-energy Hamiltonian density of the top xy surface as where ρ j act in the space of the surface states Φ xy ± . We recover the Dirac cone dispersion relation of the surface states and note that at finite momentum k = 0 the spinorbit polarization is momentum-locked. Combined with the rotational symmetry in the xy plane, this results in the anti-vortex like spin texture observed numerically in Fig. 2(c).
The low-energy description of the bottom xy surface at z = L z is obtained using the symmetry transformation z → L z − z, e j → −e j , where e j are three basis unit vectors. Under this transformation, the Hamiltonian density is modified as H xy TI → τ x H xy TI τ x . Such a transformation preserves the orientation of the surface, so that the outward-pointing normal vector v xy ⊥ always coincides with the vector e z . Using this symmetry transformation, we deduce that the wavefunctions of the states localized at z = L z read τ x Φ xy ± (L z − z). The effective low-energy Hamiltonian expressed in the basis of the surfaces states remains unaffected,

Lateral surfaces
In the same way as for the top and bottom surfaces, we express the states Φ xz ± localized on the xz surface in the neighborhood of the Dirac point, by considering (k x , k z ) = 0 and making the substitution k y → −i ∂ y . We also assume that the linear size L y of the system in y direction is sufficiently large, such that the overlap between the states at opposite surfaces can be neglected. The corresponding problem reads We note that the simplified Hamiltonian density has a new particle-hole symmetry P = η x τ z K and a new chiral symmetry C = iη x τ z σ y , which both commute with T . As a consequence, the lateral surface states also have to be at zero energy and satisfy T Φ xz The zero-energy solutions of Eq. (13) describing the state localized at the y = 0 surface have the following general form: Φ xz where C j are complex coefficients and Φ xz ±j are wavefunctions satisfying T Φ xz ±j = ±Φ xz ∓j , P Φ xz ±j = Φ xz ±j , and σ x Φ xz µj = ±Φ xz µj . Each wavefunction Φ xz ±j is associated with a decay length ξ j and a wave-vector q j , which can be calculated as ξ j = 1/Re(λ j ) and q j = Im(λ j ) [45], where λ j are solutions of Solutions of the above equation can be obtained analytically. For the sake of readability, we do not provide an exact expression here, however, we point out an important difference with respect to the lateral surface case. The solutions of Eq. (15) are, in general, complex, such that both ξ j and q j are non-zero. This corresponds to an exponentially decaying solution which oscillates, as compared to the surface states on the top and bottom surfaces which are simply decaying. Due to the high complexity of the problem, we do not provide an analytical expression of the wavefunctions Φ xz ±j and the coefficients C j . Nevertheless, we are able to solve the one-dimensional problem associated with Eq. (13) numerically. This allows us to estimate the expectation values of the remaining terms in H TI in the basis of Φ xz ± . In general, we find that the diagonal components, as well as the real and the imaginary part of the off-diagonal components, are non-zero for all the terms, as a result of the strong anisotropy between the x and z directions. However, we find that one of the contributions is usually by an order of magnitude larger than the other terms. Below, we express the diagonal components of the expectation values that accounts for this largest contribution, Here, A is a coefficient that depends on the parameters of the system and the constants C j . Similarly, we calculate the off-diagonal components As a result, the effective low-energy Hamiltonian of the surface states at y = 0 can be written as where ρ j act in the space of the surface states Φ xz ± . We recover the tilted Dirac cone dispersion relation observed in Fig. 2(d).
The effective Hamiltonian of the lateral surface at y = L y is obtained by using the transformation y → L y − y, e j → −e j which preserves the surface orientation such that the normal vector v xz ⊥ pointing outwards of the surface coincides with the vector e y . Under this transformation, the Hamiltonian density transforms as H xz TI → σ y H xz TI σ y , which allows us to express the surface state wavefunctions as σ y Φ xz ± (L y − y). We find that the low-energy effective Hamiltonian density remains invariant under this transformation Finally, the low-energy description of the two remaining yz surfaces can be obtained using the rotational symmetry in the xy plane described by the transformation {L x − x, y} → {y, x}. Using this symmetry allows us to express the wavefunctions of the surface states localized at the yz surfaces at x = 0 as iσ z Φ xz ± (x). Similarly, we find that the wavefunctions of the surface states localized at . Under this transformation, the low-energy Hamiltonian describing the two yz surfaces becomes Here we see that, the momentum component k x is replaced by k y , but the structure of the low-energy Hamiltonian remains invariant.

IV. Second-order 3D topological insulator
In this section, we describe the low-energy physics of our system under applied staggered Zeeman field (t Z = 0). For this purpose, we consider the full Hamiltonian density The crucial role of the Zeeman term H Z consists in gapping out the surface states Φ s ± . As a result, the system ceases to be a 3DTI. Nevertheless, as we will explain below in this section, the topological description of the resulting system becomes even richer, supporting the emergence of chiral hinge states characteristic to threedimensional SOTI. We also note that the spatial oscillation of the staggered Zeeman field in H Z is crucial to achieve our goal. More precisely, one can show that a uniform Zeeman field fails to gap out all the surface states, and does not lead to hinge states (see Appendix 1 for more details).

A. Small Zeeman field
In a first step, we consider the regime of a small staggered Zeeman field t Z t 1 , t 2 . This allows us to use the results of Sec. III and treat the Zeeman term H Z perturbatively.

Top and bottom surfaces
Using the expression of the wavefunctions Φ xy ± (z) of the surface states located on the top xy surface at z = 0, we calculate the expectation values of the Zeeman term H Z in the neighborhood of the Dirac point. We obtain the following diagonal contributions Similarly, we calculate the off-diagonal contributions Using these results, we write the effective low-energy Hamiltonian density describing the top xy surface of the 3DTI under the applied Zeeman field as We see that x and y components of the Zeeman field simply shift the position of the Dirac cone in momentum space. However, z component of the Zeeman field induces a 'mass term' with the amplitude M xy = t Z n z (t 1 + t 2 )/t. This mass term anticommutes with all the remaining terms in H xy SOTI and gaps out the surface states Φ xy ± . These theoretical predictions are confirmed by a numerical diagonalization of H SOTI presented in Fig. 3(a).
Next, we perform the calculations for the bottom xy surface at z = L z . We note that the projection of the Zeeman field vector n onto the normal unit vector v xy ⊥ , pointing outwards of the surface, changes sign. In other words, applying the surface-interchanging transformation from Sec. III B 1 results in the inversion of the direction of the Zeeman field n → −n. As a consequence, the expectation value of the Zeeman term in the basis of the surface states of the bottom surface is exactly opposite to the case of the top xy surface, leading to In the following, we will show an important consequence of this result.

Lateral surfaces
Using the results of Sec. III B 2, we analyze the effect of the Zeeman field on the surface states localized on the lateral xz surfaces at y = 0. Below we provide the result of a numerical calculation of the expectation value of the Zeeman term H Z in the basis of states Φ xz ± . For the diagonal components of the expectation values, we obtain Similarly, for the off-diagonal components we find This allows us to write the effective low-energy Hamiltonian of the y = 0 surface as Surprisingly, we find that the Zeeman field along the x direction does not have any effect on the surface states. The z component of the Zeeman field simply shifts the Dirac cone along the x axis in the momentum space, while the y component gaps out the surface states by inducing a mass term with the amplitude M xz = t Z Cn y . We confirm these analytical predictions by performing a numerical diagonalization of H SOTI , with the results shown in Fig. 3(b). In order to describe the low-energy physics of the y = L y surface, we use the transformation y → L y − y, e j → −e j , which preserves the orientation of the lateral surface. Under this transformation, the direction of the Zeeman field changes sign such that the new low-energy Hamiltonian becomes Finally, the effective description of the two yz surfaces is obtained by using the rotational symmetry in the xy plane. We obtain and We note that the roles of the Zeeman field components along the x and y axis are interchanged. However, the mass term still gaps out the surface states and the sign of the mass term is opposite for the two surfaces.

B. Emergence of the hinge states
As we have seen previously, the effective low-energy description of the surface states Φ s ± of the 3DTI at zero Zeeman field is characterized by the Dirac cone dispersion relation. The staggered Zeeman field opens the gap at the Dirac point by inducing a mass term with the amplitude M s . The sign of M s is equal to the projection of the Zeeman field vector n onto the outward-pointing normal unit vector v s ⊥ , which changes sign between opposite surfaces. Hence, if all the surfaces of the system are gapped, a set of gapless chiral modes has to emerge at the hinge interfaces that connect the surfaces with opposite signs of M s . Such gapless chiral modes are denoted as hinge states. We use the existence of the hinge states as a criterion to identify the emergence of a SOTI topological phase.
An exact position of the hinge states is determined by calculating v s ⊥ · n on every surface. The presence of the hinge states is confirmed numerically by diagonalizing H SOTI in a geometry with OBCs along all three directions shown in Figs. 3(c)-(d) where we consider two most general orientations of the Zeeman field. In Fig. 3(c) the Zeeman field is taken to be parallel to the xy surface, but making a finite angle with both xz and yz surfaces. As a result, the two hinge states emerge at the interfaces (x = L x , y = 0) and (x = 0, y = L y ) which connect the xz and yz surfaces with opposite signs of v s ⊥ · n. The mass term on the xy surfaces is zero, such that the whole surfaces remain gapless. In Fig. 3(d) the Zeeman field makes a finite angle with all three surfaces, gaping the xy surface out and leading to the emergence of four additional hinge states at (x = L x , z = 0), (y = L y , z = 0), (x = 0, z = L z ), and (y = 0, z = L z ). In both cases the gapless chiral state splits the whole system surface into two regions with opposite values of the mass M s . Deforming the system geometry affects the exact position of the hinge state, but does not remove it. In particular, the rectangular geometry considered above can be deformed into a sphere. Then, the Zeeman field will split the sphere into two symmetric semi-spheres with opposite signs of M s , leading to the emergence of a gapless chiral mode at the equator.

C. Phase diagram
Previously, we obtained a simple low-energy description of the model using the projection onto the space of surface states Φ s ± . This description is valid only as long as t Z t 1 , t 2 . A complete phase diagram, including the regime of strong staggered Zeeman field, is obtained numerically by calculating the gap to the first excited state as a function of the ratios t 1 /t 2 and t Z /t 2 . The result of the calculation is shown in Fig. 4(a). There we see that the 3DTI phase is stable in the regime t 1 > t 2 at strictly zero t Z , characterized by a gapless surface. At finite values of t Z , the SOTI phase emerges from the 3DTI phase, leading to the apparition of a gap to the first excited surface state. Both 3DTI and SOTI phases are gapped in the bulk. We also clearly see the critical line which separates the topologically trivial phase at strong t 2 and the topological phase at strong t 1 . The Zeeman field shifts the phase boundary, so that higher values of the ratio t 1 /t 2 are required to reach the SOTI phase. Along this critical line, the system is gapless in the bulk. The numerical calculations are performed for several directions of the Zeeman field vector n, but we find that it hardly Finally, we discuss the stability of our setup with respect to external perturbations and disorder. We verify numerically that the topological states are not affected by the local on-site disorder such as a fluctuating chemical potential or a Zeeman field generated by impurities. In our calculations the disorder is implemented using a uniform distribution with zero mean and standard deviations S µ and S Z . Similarly, we check that the presence of topological phases is not affected by the offset of fine tuned parameters in different Rashba layers, such as the SOI amplitude α and the chemical potential µ. For this purpose, each Rashba layer is equipped with a set of two random variables, describing the offset of the SOI amplitude α and the chemical potential µ. The random variables are taken from a uniform distributions with standard deviations S α and S µ . These random variables are constant within each Rashba layer and only differ between different Rashba layers. The topological phase diagram showing different phases as well as the values of bulk gaps (gap to the first excited state) in a disordered system is presented in Fig. 4(b). We find that the SOTI phase is stable against the perturbations of a strength larger than the bulk gap.

V. Modified setups
In the last section, we consider two alternative approaches to generate the topological phase hosting hinge modes, which could facilitate the experimental realization.

A. Magnetic proximity setup
From the consideration above we know that the staggered Zeeman field has an effect only on the surface as the three-dimensional bulk modes are already gapped out by the tunneling between layers. Thus, we can consider a magnetic proximity setup where the staggered Zeeman field is induced only on the surface of the 3DTI and not in the bulk. Such a surface Zeeman field gaps out the states Φ s ± by locally breaking the time-reversal symmetry at the surface. It leads to the emergence of hinge states without affecting the bulk properties of the system. The effect of the surface Zeeman field is tested numerically by calculating the gap to the first excited surface state, as a function of the ratio t Z /t 2 and the penetration length of the Zeeman field l Z . The result of the calculation is shown in Fig. 5(a). We find that even a very small penetration length, of the order of a few unit cells, is enough to fully gap out the surface states. Moreover, we observe that the topological phase transition occurs as a function of the surface Zeeman field. The position of the phase transition is not affected by the value of the penetration length l Z and is independent of the bulk properties of the system.
An experimental realization of the staggered surface Zeeman field can be achieved by using magnetic adatoms [40][41][42] or by placing an antiferromagnetic material with a bicollinear antiferromagnetic ordering close to the Rashba layer heterostructure, such that the size of the magnetic unit cell along the z direction of the ferromagnet matches the distance a, as shown in Fig. 5(b). The required antiferromagnetic materials have been studied both theoretically [46] and experimentally [47]. Such materials have also been used experimentally in a proximity setup with a two-dimensional TI [48]. This setup adds flexibility in controlling the value of the mass term M s on different surfaces, which can be achieved by adjusting the position of the antiferromagnetic material or controlling the deposition of magnetic atoms. This allows one to shift the position of the hinge state or even generate several interfaces of zero M s on a single 3DTI surface.

B. Coupled TI films
Second, we propose an alternative coupled-layer construction where the 2DEGs with opposite dispersions are substituted by thin three-dimensional TI films. A thin TI film can be interpreted as a bilayer material hosting a  Similarly to the coupled 2DEG layer system, the TI layer setup requires a staggered Zeeman field to gap out the lateral surfaces. The Zeeman field has to change sign from one surface of the TI film to another, and keep the same sign at the interface between different TIs. The simplest way to realize such a Zeeman field consists in using thin ferromagnetic layers contained between the TI films. Alternatively, the effective Zeeman field can be generated using magnetic impurities. The schematic representation of the TI layer setup subjected to a staggered Zeeman field is shown in Fig. 5(d).

VI. Conclusions
To summarize, in this work we considered a system of coupled Rashba layers subjected to a staggered Zeeman field. Such a model can be realized experimentally by using 2DEGs with opposite signs of the g-factor, by depositing magnetic adatoms or by placing an antiferromagnet with a bicollinear antiferromagnetic order close to the system surface. A similar effective model can also be obtained by considering an array of coupled thin TI layers. At zero Zeeman field, the system experiences a topological phase transition from the trivial insulator to the strong 3DTI which hosts gapless modes at its surface. Focusing on the low-energy degrees of freedom associated with the surface states, we calculated perturbatively the effect of the Zeeman term. We found that it leads to the emergence of a mass term for gapless surface excitations that gaps them out. The amplitude of the mass term is determined by the direction of the Zeeman field and inevitably changes sign at opposite surfaces of the system. It leads to the emergence of a zero-mass hinge interfaces hosting chiral gapless hinge states, characteristic for the SOTI phase. We confirmed numerically the presence of the hinge state and calculated the topological phase diagram of our model. We found that the SOTI phase is stable up to relatively large values of the Zeeman field. The topological phase diagram also remains unaffected by external perturbations and disorder of strength larger than the bulk gap.

B. Lateral surfaces
Second, using the results of Sec. III B 2, we calculate the effect of the uniform Zeeman field on the surface states localized on the lateral xz surface at y = 0. A numerical calculation of the diagonal components of the expectation values provides the following result Φ xz ± t Z n x σ x Φ xz ± = Φ xz ± t Z n y σ y Φ xz ± = Φ xy ± t Z n z σ z Φ xy Similarly, for the off-diagonal components we find Φ xz − t Z n y σ y Φ xz This allows us to write down the effective low-energy Hamiltonian of the y = 0 surface Surprisingly, we find that the only component of the Zeeman field that affects the low-energy degrees of freedom is the x component. Its main effect consists in shifting the Dirac cone along the k z axis in momentum space. We confirm these analytical predictions using numerical simulations, shown in Fig. A1(b). We see that the gap remains closed for any direction of the Zeeman field. We also note that the momentum shift originating from the uniform Zeeman field is much larger than the one originating from the staggered field.

C. Numerical confirmation
As a final confirmation of the effect of the uniform Zeeman field, we perform numerical simulations of the system in a geometry with OBC along all the three axis directions. The result of such simulations is shown in Figs. A1(c)-(d). Similar to the calculation in the main text, we consider two different geometries, with the Zeeman field parallel to the xy surface, and with the Zeeman field making a finite angle with all the three surfaces. As expected from our theoretical analysis, the z component of the uniform Zeeman field gaps out the surfaces states localized at the top and bottom surfaces. No components of the Zeeman field open a gap for the surfaces state localized at the lateral surfaces.