## Interplay between boundary conditions and Wilson's mass in Dirac-like Hamiltonians

A. L. Ara´jo, u 1 R. P. Maciel, 1 R. G. F. Dornelas, 1 D. Varjas, 2 and G. J. Ferreira 1

1 Instituto de F´ ısica, Universidade Federal de Uberlˆndia, Uberlˆndia, MG 38400-902, Brazil a a

2 QuTech and Kavli Institute of Nanoscience, Delft University of Technology,

P.O. Box 4056, 2600 GA Delft, The Netherlands

(Dated: December 6, 2021)

Dirac-like Hamiltonians, linear in momentum k , describe the low-energy physics of a large set of novel materials, including graphene, topological insulators, and Weyl fermions. We show here that the inclusion of a minimal k 2 Wilson's mass correction improves the models and allows for systematic derivations of appropriate boundary conditions for the envelope functions on finite systems. Considering only Wilson's masses allowed by symmetry, we show that the k 2 corrections are equivalent to Berry-Mondragon's discontinuous boundary conditions. This allows for simple numerical implementations of regularized Dirac models on a lattice, while properly accounting for the desired boundary condition. We apply our results on graphene nanoribbons (zigzag and armchair), and on a PbSe monolayer (topological crystalline insulator). For graphene, we find generalized Brey-Fertig boundary conditions, which correctly describe the small gap seen on ab initio data for the metallic armchair nanoribbon. On PbSe, we show how our approach can be used to find spinorbital-coupled boundary conditions. Overall, our discussions are set on a generic model that can be easily generalized for any Dirac-like Hamiltonian.

## I. INTRODUCTION

The Dirac-like Hamiltonians play an ubiquitous role in novel materials, ranging from graphene [1-3] to topological insulators (TI) [4-9], its crystalline [10-15] and higher-order [16-19] TI counterparts, and Weyl semimetals [20, 21]. The Dirac cone structure of their low-energy band dispersion leads to great interest for possible optoelectronic and spintronic devices [22-25]. Since the Dirac cone itself is well described by linear in momentum k Hamiltonians, bulk models could be limited to this leading-order contribution. However, the k 2 Wilson's terms [26-29] are required to regularize the models for the calculation of topological invariants [8, 9, 30]. Moreover, numerical (finite-differences) implementations of k -linear models face the fermion-doubling problem [31-33]. For finite systems ( e.g. , nanoribbons), the Dirac models allow for a variety of possible non-trivial boundary conditions [34-36], depending on the broken symmetry that imposes the confinement [37, 38], as initially discussed by Berry and Mondragon [39]. However, if the k 2 corrections are included in the model, one expects that the only allowed boundary condition is that of a vanishing envelope function at the edges. Therefore, we ask the following: How can the different k -linear boundary conditions be translated to models that account for the k 2 corrections?

we apply our method to the well-known graphene zigzag and armchair nanoribbons, and to the PbSe monolayer topological crystalline insulator. Indeed, graphene is an ideal material to test our findings due to (i) the formation of edge-state bands connecting the K and K ′ valleys on the zigzag case; and (ii) the contrast between gapped and gapless dispersions of armchair nanoribbons with different widths. We show that our systematic approach allows us to derive and generalize the Brey and Fertig boundary conditions [40] from symmetry constraints, which is more general than the usual analysis of the atomistic terminations, thus extending the derivation of boundary conditions to naturally include spinful systems ( e.g. , PbSe and SnTe TCIs).

In this paper, we investigate this question to show that the proper choice of the k 2 Wilson's correction induces the desired boundary conditions on Dirac-like ( k -linear) Hamiltonians. In the first part of the paper we establish our results on a generic model that applies for all Dirac-like materials. We formulate this discussion using group-theory arguments, thus emphasizing its generality, while providing a recipe on how to apply our ideas to different materials. For the Bernevig-Hughes-Zhang (BHZ) model [4], the relation between the hard-wall boundary condition and the k 2 terms was shown in Ref. 30. Here,

In a previous paper [41], our group has shown that the k 2 Wilson's mass term [26] allows for a simple elimination of the numerical fermion-doubling problem on finitedifferences implementations [31-33]. A similar proposal is established in Ref. [33], however, their choice of Wilson's term undesirably breaks time-reversal symmetry. In Ref. [41], it was suggested, and shown as a conjecture, that the proper choice of the k 2 term avoids this undesired broken symmetry. Here, we prove this conjecture and extend it to show how it can be used to either (i) derive non trivial boundary conditions for Dirac-like materials on k -linear models, or (ii) properly regularize the Dirac-models on a lattice by choosing the appropriate k 2 Wilson's term that accounts for the desired type of boundary condition.

The band structures from the effective Hamiltonians are compared with ab initio results obtained from density functional theory (DFT). The effective Hamiltonians are obtained with support from QSYMM python's package [42], and the tight-binding models are implemented with KWANT python's package [43]. All codes, input, and data files are available as Supplemental Material [44]. For the DFT simulations, we use the generalized gradient

approximation (GGA) for the exchange and correlation functional [45]. Fully relativistic j -dependent pseudopotential, within the projector augmented wave method [46], was used in the noncollinear spin-DFT formalism self-consistently. We use the Vienna ab initio simulation package (VASP) [47, 48], with plane wave basis set with a cut-off energy of 400-500 eV. The Brillouin zone is sampled using a number of k -points such that the total energy converges within the meV scale. The optimized force criteria for convergence was less than 0.01 eV/ ˚ . A

## II. GENERIC MODEL

Effective models can be obtained from symmetry constraints imposed by method of invariants [49], which is equivalent to a k · p envelope function approach, yielding a matrix expansion of the Hamiltonian H ≡ H ( k ) in powers of the momentum k . For Dirac-like materials, one might truncate the expansion on the leading order ( k -linear terms), for which the confinement is set by non-trivial boundary conditions [37-41] (see Sec. II B). However, the k 2 corrections play a significant role in numerical simulations, allowing for a simple elimination of the fermion-doubling problem [31, 32, 41]. Moreover, a conjecture introduced in Ref. 41 states that the matrix form of the k 2 term is related to the hard-wall boundary conditions [39]. In this section, we prove this conjecture on a generic, yet complete, formulation in terms of a minimal model.

To guide our discussions, let us consider a onedimensional system given by the generic Hamiltonian

$$H = \hbar { v } _ { F } U _ { k } k + m U _ { w } k ^ { 2 } + U _ { c } V ( x ), \quad \ ( 1 ) \quad ^ { \text{anti} } _ { \text{repr} }$$

which is defined along x , and k = -i∂ x . The k -linear term gives the Dirac-like dispersion at low energies with Fermi velocity v F . The k 2 correction introduces the Wilson's mass m [26, 41]. The last term is a softwall confining potential given by a symmetric profile V ( x ) = V 0 [1 -Θ( x + L ) + Θ( x -L )], where Θ( x ) is the Heaviside step function, and V 0 is the intensity. This profile defines the physical system within the inner region | x | &lt; L , while on the outer region ( | x | ≥ L ) it opens a gap 2 | V 0 | . Later, we will consider the hard-wall limit | V 0 | → ∞ , which excludes the outer region from the physical domain. The U , U k w , U c are Hermitian matrices defined by the symmetry constraints imposed on H , which will be discussed throughout the next sections. Namely, U k defines the kinetic energy term, U w sets the type Wilson's mass, and U c sets the form of confinement.

glyph[negationslash]

glyph[negationslash]

Hereafter, we assume that the U k , U c , and U w matrices are nonsingular. Neglecting the k 2 correction by setting m = 0, U k defines the unbounded Dirac-like spectrum of the Hamiltonian, i.e. , its eigenvalues give the positive and negative velocities of the Dirac cone. Thus, we require det( U k ) = 0, otherwise, one would have a flat branch in the energy dispersion. In turn, the U c matrix defines the confinement gap on the outer region where V ( x ) = 0.

glyph[negationslash]

glyph[negationslash]

Within this region, we demand that Eq. (1) does not have any propagating modes at zero energy. To guarantee this, we demand that U -1 k U c has no purely real eigenvalues, which implies det( U c ) = 0. The anticommutation condition { U , U k c } = 0 is sufficient, but not necessary, to fulfill this. A similar argument applies to the case with a Wilson's mass term in the interior of the sample, where V ( x ) = 0 and m = 0. We demand that the Wilson's mass term does not introduce any new propagating modes at zero energy beyond the ones in the original Dirac model. This results in the identical conditions for U w as for U c (for details, see Appendix A).

## A. Symmetry constraints: U w ≡ U c

Let us consider that our generic system, i.e., H from Eq. (1), is invariant under a symmetry group G , which is composed by two types of symmetry operations: S + and S -. The S + operators leave x invariant, while S -takes x →-x . Therefore, it follows the transformations S ± x S -1 ± = ± x , S ± k S -1 ± = ± k , S ± V ( x ) S -1 ± = V ( x ). The last one is a consequence of the symmetric form of V ( x ) introduced previously. Imposing that H is invariant over the full group G , ( i.e. , S ± H S -1 ± = H ), one obtains the symmetry constraints for the matrices of H :

$$[ U _ { w }, D ^ { \psi } ( \mathcal { S } _ { \pm } ) ] = [ U _ { c }, D ^ { \psi } ( \mathcal { S } _ { \pm } ) ] = 0, \text{ \quad \ } ( 2 )$$

$$[ U _ { k }, D ^ { \psi } ( \mathcal { S } _ { + } ) ] = \{ U _ { k }, D ^ { \psi } ( \mathcal { S } _ { - } ) \} = 0, \quad \ ( 3 )$$

which defines the symmetry-allowed matrices U k , U w , and U c . Here [ · , · ] and {· , ·} are the commutator and anti-commutator operations, and D ψ ( S ± ) are the matrix representations of S ± in the Hilbert space. Since U w and U c satisfy the same constraints, it follows that they are equivalent ( U w ≡ U c ), i.e. , both are in the same linear space of allowed matrices. This equivalence between U w and U c was assumed truthful, but not rigorously proven in Ref. 41.

## B. Hard-wall boundary conditions

glyph[negationslash]

The appropriate hard-wall boundary condition depends on the order of the differential equation. For our generic effective model H in Eq. (1), the Schr¨dinger o equation HF x ( ) = EF x ( ) has order 2 if m = 0, or order 1 if m = 0. Here, F x ( ) is an envelope spinor function [49]. In all cases, the energy E = 〈 F H F | | 〉 must be bounded and well defined. Consider H from Eq. (1), with a simplified single boundary profile V ( x ) = V 0 Θ( ). x On the outer region x &gt; 0 the gap | 2 V 0 | yields evanescent solutions at low energies, i.e. , F x &gt; ( 0) ∼ F e 0 -x/λ , where F 0 is the spinor at x = 0. If the k -linear term dominates the low-energy band structure, the penetration length is λ ∝ glyph[planckover2pi1] v F / V | 0 | . In the hard-wall limit | V 0 | →∞ and λ → 0, thus near the interface x ≈ 0 we can write F x ( ) ≈ F 0 [1 -Θ( )]. x Considering only the x ≈ 0 range

on the integrals in E = 〈 F H F | | 〉 , it can be shown that the contributions from the k -linear and potential V ( x ) terms are finite, while m F U k 〈 | w 2 | F 〉 ≈ -mF U F δ † 0 w 0 (0) is ill defined due to the δ (0). Therefore, either m = 0 or F 0 = 0. In the first case, one gets a k -linear model with discontinuous F x ( ), while the second case yields a k 2 model with a continuous F x ( ) that vanishes at the hard-wall interface.

The discontinuous behavior of F x ( ) in k -linear Hamiltonians ( m = 0) was first introduced in the neutrino billiards by Berry and Mondragon (BM) [39], and further discussed in Refs. 50 and 51. Later, it was applied to graphene [37, 38] and topological insulators [30, 41, 52]. Here, we cast the BM hard-wall boundary condition for H in Eq. (1) in a form that explicitly shows U k and U c as

$$\left ( i U _ { k } ^ { - 1 } U _ { c } \right ) F _ { 0 } = \frac { \hbar { v } _ { F } } { \lambda V _ { 0 } } F _ { 0 }. \text{ \quad \ \ } \left ( 4 \right )$$

This is an eigenvalue equation for the matrix iU -1 k U c with eigenvalue α = glyph[planckover2pi1] v F / ( λV 0 ). The solutions with Re λ &gt; 0 (Re α &gt; 0) describe states decaying outside the boundary. Defining M = α -1 ( iU -1 k U c ) , the boundary condition becomes

$$M F _ { 0 } = F _ { 0 }. \quad \ \ ( 5 ) \quad \underset { \iota \iota } { \circ }$$

If U c and U k are chosen such that M is traceless, unitary, and Hermitian, this reproduces the most general currentconserving boundary conditions [37-39]. We derive formally the same boundary condition replacing U c with U w for the case of a hard wall with Wilson's mass term (for details, see Appendix A). In the remaining of the paper we choose the Hamiltonian terms such that ( iU -1 k U c ) has at least one pair of eigenvalues α = ± 1, equivalently det ( U k ∓ iU c ) = 0. Using this, we cast the boundary conditions Eq. (5) for a finite system with | x | &lt; L as ( U k ∓ iU c ) F ( ± L ) = 0.

## C. Summary of the models

From the considerations above, we conclude that one can choose to work with either the k -linear formulation with nontrivial boundary conditions, or the k 2 model with trivial hard-walls:

$$H = \hbar { v } _ { F } U _ { k } k, \quad \text{ with } ( U _ { k } \mp i U _ { c } ) F ( \pm L ) = 0, \ ( 6 ) \quad \ v b ( r )$$

$$H = \hbar { v } _ { F } U _ { k } k + m U _ { w } k ^ { 2 }, \quad \text{ with } F ( \pm L ) = 0. \ ( 7 )$$

These two approaches are equivalent, due to the analytical properties of the boundary conditions discussed above, and the equivalence U w ≡ U c shown in Sec. II A. Thus, for the k -linear model of Eq. (6), the characteristics of the confinement appear on the BM boundary condition, while on the k 2 model of Eq. (7), it enters through Wilson's k 2 term. For the k 2 model, the Wilson's mass can be set within

$$\frac { \delta _ { x } ^ { 2 } \delta _ { \varepsilon } } { 2 } \leq | m | \leq \frac { ( \hbar { v } _ { F } ) ^ { 2 } } { \delta _ { \varepsilon } },$$

which ensures that the k -linear term dominates the lowenergy spectrum within the δ ε energy range and eliminates the fermion doubling problem on a lattice [41]. In the analytical limit, the pure Dirac-like model is restored as m → 0 [53].

We explicitly show in Appendix A that the envelope functions of the two models match far from the boundaries. Near the boundary, the k 2 model has an additional localized contribution that decays quickly away from the boundary in the limit of small m [53].

## III. GRAPHENE NANORIBBONS

Graphene nanoribbons (Fig. 1) are ideal cases to present our findings on a concrete system, since its band structure and boundary conditions are well known [3, 40]. In Ref. 40, Brey and Fertig (BF) have shown that the graphene k -linear boundary conditions depend on the nanoribbon atomic termination. Particularly, their boundary condition predicts a metallic armchair nanoribbon with an identically zero gap for N A = 3 p +2 (with integer p ), while the ab initio data from Ref. 54 shows a vanishing, but finite gap. Indeed, the null gap is not consistent with its symmetry group ( ∼ D 2 h , see Appendix B) of the nanoribbons.

In this section we systematically revise the BF boundary conditions. We find that it is equivalent to the BM boundary condition, given by a proper choice of U c , which is imposed by symmetry [37, 38]. However, for the metallic armchair case, we find that U c diverges, which is a consequence of the zero gap inconsistency mentioned above. Therefore, we propose a generalization of the BF boundary condition that fixes this inconsistency in both k -linear and k 2 approaches.

Initially, let us consider the usual k -linear Dirac model of a full monolayer graphene. Later in Section III C, we introduce the k 2 model for the nanoribbons. The k · p expansion for graphene considers basis functions given by its solutions at K and K ′ valleys, i.e., ϕ A ( r ), ϕ B ( r ), ϕ A ′ ( r ) and ϕ B ′ ( r ), where A and B label the sublattices. These are illustrated in Fig. 1(a). Within the envelope function approximation [55-57], the expansion reads as

$$\psi ( r ) = \sum _ { \mu } \left [ f _ { \mu } ( r ) e ^ { i q \cdot r } \varphi _ { \mu } ( r ) + f _ { \mu ^ { \prime } } ( r ) e ^ { i q ^ { \prime } \cdot r } \varphi _ { \mu ^ { \prime } } ( r ) \right ], \ ( 9 )$$

where f µ ( r ) and f µ ′ ( r ) are the envelope functions, µ = { A,B } label the sublattices, q = k -K and q ′ = k -K ′ are the deviations from the K and K ′ valleys in k -space. Bloch theorem requires ψ ( r + R ) = e i k R · ψ ( r ), where R = n 1 a 1 + n 2 a 2 (with n 1 , n 2 integers) is a Bravais translation of the monolayer. Since the Bloch phase in ϕ µ ( ) ′ ( r + R ) = e i K ( ) ′ · R ϕ µ ( ) ′ ( r ) cancels out the opposite

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-6038f4d5-d231-4a95-9004-6c3de18db36e.jpg)

A

Figure 1. (a) Hexagonal cell of monolayer graphene with the A (blue) and B (red) sublattices emphasized. The K -point basis functions ϕ µ ( r ) (with µ = A,B ) for each sublattice is composed by p z orbitals multiplied by the Bloch phases indicated at each site, with φ = 2 π/ 3. For K ′ , ϕ µ ′ ( r ) are composed replacing φ → -φ . (b) Projections of the K , K ′ , and Γ points of the Brillouin zone into its one-dimensional counterparts for zigzag ( k y ) and armchair ( k x ) nanoribbons. Lattices for (c) zigzag and (d) armchair nanoribbons with N D ( N A ) dimers (atoms) from edge to edge. Their primitive vectors a Z/A and unit cells are highlighted. The carbon-carbon distance is a/ √ 3 ≈ 0 142 nm. . The empty circles at the edges indicate the lattice sites that were removed to form each nanoribbon, defining the effective dimensions W = W 0 + 2 a/ √ 3 and L = L 0 + a , where W 0 = (3 N / D 2 -1) a/ √ 3 and L 0 = ( N A -1) a/ 2.

phase in the q -exponentials, Bloch theorem is satisfied for a periodic f µ ( ) ′ ( r + R ) = f µ ( ) ′ ( r ). Up to leading order in k , the usual effective Dirac-like Hamiltonian is

$$H _ { G } ( k ) = h ( q ) \oplus h ^ { * } ( q ^ { \prime } ), \quad \ \ ( 1 0 ) \quad \ \psi ( j )$$

$$h ( q ) = \hbar { v } _ { F } \sigma \cdot q. \quad \ \ ( 1 1 ) \quad _ { \ f o c } ^ { J A }$$

The 4 × 4 H G ( k ) Hamiltonian acts on the envelope spinor F ( r ) = [ f A ( r ) , f B ( r ) , f A ′ ( r ) , f B ′ ( r )]. Notice that we write H G ( k ) in terms of the deviations q and q ′ , such that the Dirac cones occur at k ∼ K and K ′ . This notation will be useful to keep track of the nanoribbon confinement projections onto the k x (armchair) or k y (zigzag) axis in the next sections. There, the projections will retain the overall form of the ψ ( r ) expansion above, but they will change the definitions of q and q ′ .

## A. Revised Brey and Fertig boundary conditions

An elegant approach to the boundary conditions for graphene nanoribbons was introduced in Ref. [40] by Brey and Fertig. There, they propose that the envelope function must vanish at the sites that were removed to form the nanoribbons. Consequently, it depends on the atomic terminations, rendering different boundary conditions for the zigzag and armchair cases. We have recently used this approach to obtain boundary conditions for topological crystalline insulators [58].

Next, we revise and generalize these boundary conditions for zigzag and armchair nanoribbons, and in the next section we show their equivalence to the U k and U c matrices on the BM approach. Complementary, the boundary conditions for confinement in arbitrary directions (beyond the zigzag and armchair) and different atomistic terminations were studied in Ref. [38]. Their results can be used to define general U k and U c matrices.

## 1. Zigzag nanoribbons

To model the zigzag nanoribbons, we start from the bulk basis functions ϕ µ,µ ′ ( r ), but with a modified ψ ( r ) expansion. Namely, replace q ( ) ′ → k -¯ K ( ) ′ in Eq. (9), where ¯ K ( ) ′ are the K and K ′ projections into the zigzag k y axis, as shown in Fig. 1(b). These k projections also apply to q ( ) ′ in H G ( k ) from Eq. (10). For simplicity, here we consider only this pristine form of H G ( k ). However, for narrow ribbons, symmetry-allowed corrections due to the finite size of the nanoribbons are relevant for an improved fit with ab initio data (see Appendix B).

The zigzag nanoribbon lattice is illustrated in Fig. 1(c). The right (left) edge ( x = ± W / 0 2) is composed only of atoms from the A ( B ) sublattice. The length W 0 = ( 3 2 N D -1) a/ √ 3, where N D is the number of dimers. The next line of atoms, removed to form the ribbon, would have been located in x = ± √ W/ 2, which defines the effective length W = W 0 + 2 a/ 3. Since the p z orbitals are highly localized at each carbon atom, the absence of B atoms at the left edge x = -W/ 2 yields ϕ B ( ) ′ ( R -) ≈ 0, where R ± = ( ± W/ ,y 2 ). Imposing ψ ( R -) = 0 in Eq. (9), we get f A ( R -) e -i ¯ K R · -ϕ A ( R -)+ f A ′ ( R -) e i ¯ K R ′ · -ϕ A ′ ( R -) = 0. Moreover, the phase factors [see Fig. 1(a)] and orbitals cancel out, i.e. , e -i ¯ K R · -ϕ A ( R -) = e i ¯ K R ′ · -ϕ A ′ ( R -). Similar considerations follow for the B sublattice on the right edge, r = R + . Due to translational invariance along y , we can simplify f µ ( r ) → e ik y y f µ ( x ). Therefore, the boundary conditions for the zigzag nanoribbon envelope functions are

$$f _ { A } \Big ( - \frac { W } { 2 } \Big ) & = - f _ { A ^ { \prime } } \Big ( - \frac { W } { 2 } \Big ), \\ f _ { B } \Big ( + \frac { W } { 2 } \Big ) & = - f _ { B ^ { \prime } } \Big ( + \frac { W } { 2 } \Big ). \quad \quad ( 1 2 )$$

In the original discussion by BF [40], they consider only a single valley ( K or K ′ ) on the ψ ( r ) expansion. Indeed, near the K valley f µ ′ ( r ) ≈ 0 (with µ = A,B ). Neglecting these contributions in Eq. (12), one immediately recovers their well known result f A ( -W/ 2) = f B (+ W/ 2) = 0. Similarly, for the K ′ valley, one gets

f A ′ ( -W/ 2) = f B ′ (+ W/ 2) = 0. Therefore, Eq. (12) is a generalization of the original BF boundary condition for the zigzag nanoribbon.

## 2. Armchair nanoribbons

For the armchair nanoribbons, confinement along y projects both K and K ′ valleys into Γ on the ¯ k x armchair axis [Fig. 1(b)]. Therefore, the ψ ( r ) expansion in Eq. (9) is defined by q = q ′ → k . As in the zigzag case above, these projections also apply to H G ( k ) in Eq. (10).

glyph[negationslash]

The armchair lattice is shown in Fig. 1(d). Confinement along y defines the length L 0 = ( N A -1) a/ 2, where N A is the number of atoms along the ribbon. The line of atoms removed to form the ribbon defines the effective length L = L 0 + a . Differently from the previous case, here both edge terminations contain atoms from the A and B sublattices. Consequently, all orbitals ϕ µ ( x, ± L/ 2) = 0. In this case the system has translational invariance along y , thus f µ ( r ) → e ik x x f µ ( y ). Imposing ψ x, ( ± L/ 2) = 0 for all x in Eq. (9), we get

$$f _ { \mu } \Big ( \pm \frac { L } { 2 } \Big ) = - e ^ { \pm i \theta } f _ { \mu } ^ { \prime } \Big ( \pm \frac { L } { 2 } \Big ), \quad \ \ ( 1 3 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where e ± iθ = ϕ µ ′ ( x, ± L/ 2) /ϕ µ ( x, ± L/ 2) is the phase difference between the K and K ′ solutions at the edges. On the BF approach [40], they consider the bulk Bloch phase difference, yielding θ → θ BF = ∆ K L · = ( N A +1)2 π/ 3. However, since confinement breaks the Bloch periodicity along y , deviations from θ BF should be expected for narrow ribbons. Therefore, hereafter we consider Eq. (13) with an arbitrary θ as a generalization of the original BF boundary condition.

## B. Equivalence between Brey and Fertig and Berry and Mondragon

The Brey-Fertig boundary condition discussed above, can be equally understood via the Berry-Mondragon formalism summarized in Eq. (5). To verify this equivalence, let us compare the BF and BM boundary conditions for the zigzag and armchair nanoribbons. Within the BM approach, the matrix U k multiplies the momentum k along the confinement direction, while U c must satisfy the symmetry constraints presented previously. For the zigzag nanoribbons, these are

$$U _ { k } ^ { Z } & = \begin{pmatrix} \sigma _ { x } & 0 \\ 0 & \sigma _ { x } \end{pmatrix}, \quad & \text{A} \\ U _ { c } ^ { Z } ( \eta ) & = \begin{pmatrix} \sigma _ { y } \eta & \sigma _ { y } ( 1 - \eta ) \\ \sigma _ { y } ( 1 - \eta ) & \sigma _ { y } \eta \end{pmatrix}, \quad & ( 1 4 ) \quad & \text{corr}$$

where σ ν (with ν = 0 , x, y, z ) are Pauli matrices acting on the sublattice A B / subspace, and σ 0 is the 2 × 2 identity matrix. The parameter η is restricted to η = { 0 1 , } , which identifies the particular boundary condition, as it is discussed below. In turn, for the armchair confinement,

$$U _ { k } ^ { A } & = \begin{pmatrix} \sigma _ { y } & 0 \\ 0 & - \sigma _ { y } \end{pmatrix}, \\ U _ { c } ^ { A } ( \theta ) & = \begin{pmatrix} \sigma _ { y } \cot \theta & \sigma _ { y } \csc \theta \\ \sigma _ { y } \csc \theta & \sigma _ { y } \cot \theta \end{pmatrix}. \quad \quad ( 1 5 )$$

Armchair nanoribbons . Substituting U A k and U A c ( θ ) from Eq. (15) into Eq. (5), one obtains Eq. (13) after straightforward manipulations. This establishes the equivalence between the BF and BM approaches. Interestingly, these are two drastically distinct approaches for the boundary condition. On the BF approach, one uses the atomistic terminations of the lattice to motivate the boundary condition. On the other hand, the BM approach is based solely on the symmetries of the lattice.

Zigzag nanoribbons . First, for η = 0, replacing U Z k and U Z c (0) from Eq. (14) into Eq. (5), we reproduce our boundary conditions shown in Eq. (12), plus an additional pair of equations f B ( R -) = f B ′ ( R -), and f A ( R + ) = f A ′ ( R + ), which are trivially satisfied since, ϕ A ( ) ′ ( R + ) ≈ 0 and ϕ B ( ) ′ ( R -) ≈ 0. Second, for η = 1 the same procedure gives us the original BF boundary condition for the zigzag confinement. Therefore, the two possible values of η = 0 1 { , } label our generalized boundary condition, and the original BF result. A comparison between the band structures of these two cases will be shown in the next section.

## C. The k 2 model for graphene nanoribbons

To construct the k 2 model for the nanoribbons, we consider the symmetries of the armchair and zigzag lattices [Figs. 1(c) and 1(d)]. They are both invariant under the D 2 h group. Additionally, we consider that the system is time-reversal symmetric and chiral. The matrix representations of these symmetry operations are built from the ψ ( r ) expansion in Eq. (9). These, and the corresponding most general symmetry-allowed Hamiltonian, are shown in Appendix B. Next, we first discuss a minimal k 2 model, which matches the usual k -linear model, but allows for simpler numerical implementations. Later, we use the full k 2 model to fit and compare the results with the DFT data.

## 1. Minimal k 2 models

A minimal model for the nanoribbons must contain only the bulk-like k -linear terms and the necessary k 2 corrections, which reads as

$$H _ { Z } = h ( \mathbf k ^ { + } ) \oplus h ^ { * } ( \mathbf k ^ { - } ) + \frac { m _ { Z } } { 2 } U _ { c } ^ { Z } ( \eta ) k _ { x } ^ { 2 }, \quad ( 1 6 )$$

$$H _ { A } = h ( { \boldsymbol k } ) \oplus h ^ { * } ( { \boldsymbol k } ) + \frac { m _ { A } } { 2 } \overset { 2 } { U } _ { c } ^ { A } ( \theta ) k _ { y } ^ { 2 }, \text{ \quad \ } ( 1 7 )$$

for zigzag and armchair nanoribbons, respectively. The Dirac-like term h ( k ) is given in Eq. (10), and k ± = ( k , k x y ± ¯ ) K centers the Dirac cones into the K and K ′ projections [Fig. 1(b)]. The k 2 terms are defined by the Wilson's matrices U Z A ( ) w ≡ U Z A ( ) c and masses m Z A ( ) . Since the k 2 terms allow for a trivial boundary condition, i.e., ψ = 0 at the edges, it can be numerically implemented via simple finite-differences schemes [41]. Hereafter, all results were obtained discretizing the coordinates with ∼ 100 points, and the Wilson's masses are chosen on the mid-range of Eq. (8).

For zigzag nanoribbons, Fig. 2 compares the generalized ( η = 0) and original ( η = 1) BF conditions. Indeed, the ideal results are given by η = 0 in Fig. 2(a). Here, the absence of fermion doublers is due to their hybridization near ¯ X [Fig. 2(b)], which drives the doublers towards high energies as m Z increases. On the other hand, for η = 1 the Hamiltonian H Z splits into uncoupled ¯ K and ¯ K ′ blocks, each showing an independent doubler [Fig. 2(c)]. Since there is no hybridization, the doublers will always occur at k doubler y ≈ ± [ ¯ K -2 m / Z ( glyph[planckover2pi1] v F δ 2 x )]. Not

Figure 2. Zigzag bands ( N D = 24) from the k 2 model for different boundary parameters η and Wilson's masses m Z . (a) For η = 0 and finite m Z , the generalized BF boundary condition returns the expected zigzag band structure. (b) For smaller m Z → 0 15 . m Z the doublers appear near ¯ X as hybridized extra cones. (c) For η = 1 (original BF) the doublers at k doubler y ≈ ± [ ¯ K -2 m / Z ( glyph[planckover2pi1] v F δ 2 x )] cannot be eliminated (dotted lines). The arrows in (a) and (c) point to the state ( E,k y ) used to plot the envelope functions F ( x ). Near ¯ , K F ( x ) for both (d) η = 0 and (e) η = 1 are similar and smooth (solid lines). (e) The F ( x ) for the doublers (dotted lines) show nonphysical phase oscillations on the scale of the numerical discretization.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-66fbf307-69ed-405e-8673-0c62754fcd66.jpg)

even staggered-lattice implementations are able to fully eliminate the doublers in this case [31, 32]. Notice that in the analytical k -linear limit δ x → 0 and the doublers are eliminated as k doubler y →±∞ , thus justifying the use of the original BF boundary conditions in analytical calculations. For k y ≈ ¯ , K the envelope functions F r ( ) for η = { 0 1 , } match well, as shown in Fig. 2(d)-(e). For η = 0, a small ¯ K -¯ K ′ hybridization occurs, which is absent for η = 1 due to the block form of H Z . In contrast, the F r ( ) of the doublers show nonphysical oscillations with a period given by the numerical step size δ x [dotted lines in Fig. 2(e)].

For the armchair nanoribbons, the boundary condition is set by θ . The original BF θ → θ BF = ( N A +1)2 π/ 3 gives a qualitatively correct picture for the armchair gap ∝ 1 /L . However, it predicts that the bands for N A = 3 p and 3 p +1 (for integer p ) are degenerated, while for N A = 3 p +2 it is identically gapless [54]. However, as discussed above, the confinement breaks Bloch periodicity and a deviation from θ BF is expected. Therefore, in Fig. 3 we consider θ = θ BF +∆ . Within the θ k -linear model, the armchair gaps are

$$E _ { 3 p } \approx \frac { \hbar { v } _ { F } } { L } \left [ \frac { 2 \pi } { 3 } - 2 \Delta \theta \right ], \text{ \quad \ \ } \text{ (1 8)}$$

$$E _ { 3 p + 1 } \approx \frac { \hbar { v } _ { F } } { L } \left [ \frac { \overset { \cdot } { 2 \pi } } { 3 } + 2 \Delta \theta \right ], \text{ \quad \ \ } ( 1 9 )$$

$$E _ { 3 p + 2 } \approx \frac { \hbar { v } _ { F } } { L } \left [ 2 \Delta \theta \right ].$$

Indeed, the DFT data for the gaps obey these expres-

Figure 3. (a) The DFT gaps (circles) as a function of L (log scale) fall into straight lines on a wide range 3 ≤ N A ≤ 50. The dashed lines are fits used to obtain the linear coefficients that yields the constraints between (b) the Fermi velocity and the correction ∆ θ to the BF boundary conditions. (c) Comparison between the gaps in DFT (circles), k -linear (dotted lines) and k 2 models (solid lines). The dotted match the solid ones, but are shifted for clarity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-3c346e09-62a1-4283-a4e7-6b09fbbadf7f.jpg)

sions, i.e. , E N A ∝ 1 /L , as shown in log-log scale in Fig. 3(a). The linear coefficient of these lines gives us constraints between v F and ∆ , which we use to establish θ the correction ∆ θ to the BF boundary condition shown in Fig. 3(b). Considering v F = 0 8 . × 10 3 nm/ps (see fits in the next section), we solve the k -linear and k 2 models with the corresponding ∆ θ ≈ 20 , ◦ and both match well the DFT gaps in Fig. 3(c). Here, we have used a constant v F for all N A for simplicity. However, v F may change as a function of N A due to finite-size effects. Consequently, ∆ θ must also be N A dependent, while obeying the constraints from Figs. 3(a) and 3(b).

The comparisons above show that our k 2 models [Eqs. (16) and (17)] generalize the usual graphene Dirac model and the BF boundary conditions. On the zigzag case, the η = 0 model couples the ¯ K and ¯ K ′ valleys, such that the edge-state branch is restricted to its correct interval in Fig. 2(a), while in the usual BF case they are uncoupled, yielding edge-state branches that extend toward k y → ±∞ . Indeed, both models would match identically if the valley projections are driven far apart in Eq. (16). For the armchair case, the BF model is recovered for ∆ θ → 0, yielding the zero gap for the N A = 3 +2 p metallic case, and degenerate 3 p and 3 p +1 gaps.

## 2. Full k 2 model: Fitting the DFT data

While the minimal k 2 model above provides a sufficient approach to regularize the Dirac models on a lattice, it is also insightful to investigate the full k 2 model in comparison with DFT results. The derivation of the most general symmetry-allowed Hamiltonian for graphene up to k 2 is shown in Appendix B. In a general compact notation, it reads as

$$H = h ( \mathbf k ^ { + } ) & \oplus h ^ { * } ( \mathbf k ^ { - } ) \\ + & \left [ \frac { m _ { A 1 } } { 2 } U _ { 0 y } + \frac { m _ { A 2 } } { 2 } U _ { x y } \right ] k _ { y } ^ { 2 } + \left [ \frac { m _ { Z 1 } } { 2 } U _ { 0 y } + \frac { m _ { Z 2 } } { 2 } U _ { x y } \right ] k _ { x } ^ { 2 } \\ & \quad + \Delta U _ { x y } + \hbar { \mu } U _ { x x } k _ { x } + m _ { x y } U _ { z x } k _ { x } k _ { y }. \quad ( 2 1 ) \\ & \quad + \Delta U _ { x y } + \hbar { \mu } U _ { x x } k _ { x } + m _ { x y } U _ { z x } k _ { x } k _ { y }. \quad ( 2 1 )$$

The matrices U ij = τ i ⊗ σ j are set in terms of the Pauli matrices τ ν (with ν = 0 , x, y, z ) acting on the K K / ′ valley subspace, and σ ν acting on the sublattices A B / . The first line in Eq. (21) contains the Dirac-like terms. For armchair k ± = k , while for zigzag it sets the valley projections as discussed in the previous section. The second line shows the most general form of the k 2 terms. For armchair ribbons, the BM boundary conditions constrain m A 1 = m A cot θ and m A 2 = m A csc θ . Similarly, for zigzag m Z 1 = m η Z and m Z 2 = m Z (1 -η ). The third line shows the extra terms that allow for a fine tuning of the band structure. The ∆ term couples the projected cones from K and K ′ valleys, the velocity µ couples the dispersions at finite k , and m xy is a trigonal correction for the masses. To illustrate the results from the full k 2 model, we have considered a medium-sized armchair nanoribbon with N A = 48 (type 3 p ), as shown in Fig. 4.

The parameters used to obtain the figures, and equivalent results for N A = 49 (3 p +1), N A = 50 (3 p +2), and for the zigzag case are shown in Appendix B.

In Fig. 4(a)-(c) we compare the DFT band structure with the k -linear BF model, tight-binding (implemented with the Kwant code [43]), and our full k 2 model, respectively. As expected, at low energies all models agree reasonably well with the DFT data. However, for | E | glyph[greaterorsimilar] 0 5 eV discrepancies are visible in all cases. . The DFT data show two sets of quantized cone dispersions with different parabolicities. This is not captured by the k -linear model. The tight binding model captures these features, but the band edges are shifted. The k 2 model provides a better fit up to | E | ∼ 0 75 eV. .

The densities of the first and second conduction subbands are shown in Figs. 4(c) and 4(d). The DFT data show peaks at atomic positions. For the models, the envelope functions are extracted from Eq. (9), | ψ y ( ) | 2 ∝ | f µ ( y e ) iθy/L + f µ ′ ( y e ) -iθy/L | 2 , where the phase factors arise from the ϕ µ ( ) ′ ( r ) phases in Fig. 1(a). These are highly oscillating envelopes, thus in Figs. 4(c) and 4(d) we plot them only at the atomic positions, showing an excellent agreement with the DFT data. These also agree with the tight binding densities [59] for the low-energy subbands.

Figure 4. Comparison between the DFT band structures and (a) the k -linear, (b) tight-binding, and (c) k 2 models for N A = 48. At low energy all models agree, while for | E | glyph[greaterorsimilar] 0 5 eV . parabolic corrections become relevant and broken chirality starts to develop. The envelope densities (colored symbols) match well the DFT data for (d) the first and (e) second conduction bands at k x = 0.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-bceca9b4-d2cf-431b-aeb9-3895c08609fc.jpg)

## IV. SPINFUL CASE: TOPOLOGICAL CRYSTALLINE INSULATORS

Graphene has a very weak spin-orbit coupling [1-3]. Therefore, to illustrate our results on a spinful system, let us instead consider a monolayer of PbSe, which is a topological crystalline insulator (TCI) [13, 60, 61]. The effective model for this material was derived in Ref. [58] up to k 2 . Hereafter, we follow the notation from this reference. For simplicity, we restrict the discussion to the PbSe nanoribbons of types A, B and C (Fig. 5), while generalizations for ribbons D and E (defined in Ref. [58] and not shown here) are straightforward.

The lattices from ribbons A and B are invariant under the point group D 2 h , while for ribbon C the symmetry is reduced to C 2 v . Therefore, the sole difference between ribbons A and B is their atomic terminations [see Fig. 5], which shall reflect on their boundary conditions. For ribbon C, one edge is equivalent to that of ribbon A, and the other is of the B type. Consequently, due to the reduced symmetry, ribbon C admits extra terms in the Hamiltonian.

The model of PbSe monolayers [58] around the X point of the Brillouin zone is defined on the basis functions { ϕ xz, ↑ ( r ) , ϕ xz, ↓ ( r ) , ϕ x, ↑ ( r ) , ϕ x, ↓ ( r ) } , where the xz and x indices refer to the symmetries of the orbitals, and ↑ ↓ , label the spin states along z . Thus, similarly to Eq. (9), the wave-function ψ ( r ) expansion is given by these basis functions multiplied by the k -phase e i ( k -X r ) · and the envelope spinor F ( x, y ) = e ik x x [ f xz, ↑ ( y , f ) xz, ↓ ( y , f ) x, ↑ ( y , f ) x, ↓ ( y )] T . Here we already assume a plane-wave along x , since the confinement in along y in Figs. 5(a)-(c).

Considering an isotropic limit for simplicity, the minimal k 2 effective model for ribbons A, B and C, confined along y and extended along x , is

$$H = \Delta U _ { z 0 } + \alpha ( U _ { x x } k _ { y } - U _ { x y } k _ { x } ) + \frac { m } { 2 } U _ { c } ( \rho, \theta ) k _ { y } ^ { 2 } + \Delta _ { C } U _ { y x }, \quad \ \ \ ( 2 2 ) \quad \ \ ^ { - 2 }$$

$$U _ { c } ( \rho, \theta ) = \cos \theta \left [ \cosh ( \rho ) U _ { z 0 } + \sinh ( \rho ) U _ { 0 0 } \right ] + \sin ( \theta ) U _ { y x }. \quad \begin{array} { c } \infty \\ ( 2 3 ) & \stackrel { c } { \underbrace { 1 } } \end{array} \right ] _ { c }$$

The matrices U ij = τ i ⊗ σ j are set by Pauli matrices acting on the orbital ( τ ν ) and spin ( σ ν ) subspaces, ∆ is the gap at k = 0, α defines the Fermi velocity, and m is Wilson's mass. The coupling ∆ C is only allowed for ribbon C. Since the confinement is along y the kinetic matrix U k = U xx , yielding det[ U k ± iU c ( ρ, θ )] = 0, as expected. The boundary conditions are defined by the continuum parameters ρ and θ . For ribbons A and B, θ ≡ 0.

The effects of the boundary condition parameter ρ and θ on the band structure are shown in Figs. 5(d)-5(f). Here, we consider ∆ and α/ ∆ as the energy and distance units. For m &gt; 0 the system is trivial, thus there are no states within the gap | E/ ∆ | &lt; 1 in Fig. 5(d). For glyph[negationslash]

m&lt; 0 the system becomes topologically non-trivial with a mirror Chern number n M = -2 [58], yielding two Dirac cones (degenerate for θ = 0). In this case, the states seen within the gap in Fig. 5(e) refer to the crossing point of the Dirac dispersion seen in Figs. 5(e1) and 5(e2). The system is chiral for ρ = 0, while for finite ρ the broken chirality is a consequence of the distinct atomic terminations of ribbons A and B [58]. For ribbon C we can consider ρ = 0 for simplicity and allow θ to vary. This is shown in Fig. 5(f). In this case, both θ = 0 or ∆ C = 0 break the degeneracies between the Dirac crossings, as seen in Fig. 5(f1). These results are equivalent to those from Ref. [58], where a BF-type boundary condition was proposed.

glyph[negationslash]

Complementarily, within the k -linear model ( m = 0), the BM approach for the boundary conditions [ U k ±

Figure 5. Lattices of PbSe ribbons of types (a) A, (b) B, and (c) C, as defined in Ref. [58]. (d)-(f) PbSe band edges at k x = 0 varying the boundary condition parameters ρ and θ . (d) For m &gt; 0, the system is trivial, there are no states within the gap | E/ ∆ | &lt; 1. (e) For m&lt; 0, a pair of degenerate topological Dirac crossings appear, and its crossing point at k x = 0 is controlled by ρ : (e1) for ρ &gt; 0 the crossing is downshifted, and (e2) for ρ &lt; 0 it shifts up in energy. (f) A finite θ or ∆ C splits the Dirac crossings as shown in (f1).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-d7480b6c-dd65-4ec1-a9c7-195df010de5e.jpg)

iU c ( ρ, θ )] · F ( x, ± L/ 2) = 0 yields,

$$f _ { x, + \sigma } \Big ( \pm \frac { L } { 2 } \Big ) = \mp \frac { i e ^ { \rho } \cos \theta } { 1 \pm \sin \theta } f _ { x z, - \sigma } \Big ( \pm \frac { L } { 2 } \Big ). \quad ( 2 4 ) \quad \text{with}$$

Interestingly, this boundary condition implies a spinorbital admixture, as it couples opposite spins ± σ and orbitals x/xz [62]. In Ref. [52], this type of constraint leads to a spin texture across the ribbon.

## V. CONCLUSIONS

We have shown how the k 2 Wilson corrections not only regularize the Dirac-like models on a lattice, but are directly related to the boundary conditions of finite systems. Considering the symmetries of the finite size system ( e.g., nanoribbons), the choice of Wilson's corrections are not arbitrary. Indeed, we show that the symmetry-allowed k 2 terms are equivalent to the nontrivial boundary conditions from Berry and Mondragon [39], thus providing a recipe to regularize the Dirac model by including the k 2 term compatible with the desired boundary condition. This hidden connection between Wilson's k 2 term and the boundary conditions were taken as a conjecture in Ref. [41] to propose a simple method to eliminate the fermion doubling problem. Here, our systematic derivation now proves this conjecture.

Applying this methodology for graphene, we have found a generalization of the Brey-Fertig boundary conditions [40]. For the zigzag nanoribbons, the K K -′ coupling induced by the boundary condition restricts the edge-state bands to lie within these valleys. More interestingly, for the armchair case, it introduces ∆ θ as a deviation from the bulk Bloch phases. Particularly, for the 'metallic' armchair case a finite ∆ θ eliminates the nonphysical gapless band structure [54]. Additionally, for the spinful systems ( e.g., PbSe TCI) our approach allows for simple derivation of the spinful boundary conditions.

## VI. ACKNOWLEDGEMENTS

We thank M. Novaes and L. Chico for useful discussions. A.L.A., R.P.M, R.G.F.D., and G.J.F. acknowledge financial support from the Brazilian funding agencies CNPq, CAPES, and FAPEMIG, and D.V. acknowledges funding from the Dutch national science organization (NWO) VIDI grant 680-47-53.

A.L.A. and R.P.M. have contributed equally to this work.

## Appendix A: Rigorous proofs of the hard-wall boundary conditions

Section § II B of the main text presents an illustrative picture on how the boundary conditions change from the k -linear to the k 2 model. Here, we provide more rigorous proofs of the equivalence of these boundary conditions for the general case of N × N (for even N ) Hamiltonians with Dirac-like spectrum.

## 1. Constraints on U c and U w

In the linear model, to define a hard wall, U c needs to open a gap at zero energy in the outside region. Consequently, the Hamiltonian ( glyph[planckover2pi1] v F U k k + V U 0 c ) should have no propagating modes at zero energy, meaning the determinant cannot have any zeros for real k . Using that U k is nonsingular, this is equivalent to U -1 k U c having no real eigenvalues, which shows that U c cannot be singular. Demanding that U c anticommutes with U k is sufficient to fulfill this condition, as this ensures that U -1 k U c is skewHermitian with purely imaginary eigenvalues.

Similarly, in the quadratic case, the Wilson's mass term cannot introduce any new modes, i.e., ( glyph[planckover2pi1] v F U k k + mU k w 2 ) should have no propagating modes at zero energy beyond the N original k = 0 modes of the linear model. Furthermore, the Wilson's mass term should ensure that in the discretized model [ k → sin( δ x k /δ ) x ] has no fermion doubling at k = π/δ x , requiring nonsingular U w . Dividing the Hamiltonian by k we get ( glyph[planckover2pi1] v F U k + mU k w ) whose determinant should not vanish for any real k . Using that U k and U w are non-singular, this is equivalent to U -1 k U w having no real eigenvalues.

Hence, we showed that U c and U w obey the same constraint, and demanding anticommutation with U k is sufficient to fulfill it.

## 2. Boundary conditions for the k -linear model

Consider the general model Hamiltonian from Eq. (1) in the k -linear case ( m = 0), i.e. , H = glyph[planckover2pi1] v F U k k + U V c ( x ), and its corresponding Schr¨dinger o equation HF x ( ) = EF x ( ). Let us assume that V ( x ) = V 0 Θ( ) defines a x single wall at x = 0, opening a gap | 2 V 0 | for x &gt; 0. On the outer region, the hard-wall limit ( V 0 → ∞ ) allow us to neglect E glyph[lessmuch] V 0 , and considering the ansatz F x &gt; ( 0) ∼ F e 0 -x/λ (Re λ &gt; 0), we get

$$\left ( \frac { i \hbar { v } _ { F } } { V _ { 0 } \lambda } U _ { k } + U _ { c } \right ) F _ { 0 } = 0. \text{ \quad \ \ } ( \text{A1} )$$

Non-trivial solutions will only exist if the determinant of the matrix in parentheses vanishes, which can be cast as

$$\det \left ( U _ { k } ^ { - 1 } U _ { c } - \omega \right ) = 0, \text{ \quad \ \ } ( A 2 )$$

since U k is non-singular. Here, ω = -i glyph[planckover2pi1] v F /V λ 0 plays the role of the eigenvalues of U -1 k U c . The restriction Re λ &gt; 0, imposed by the evanescent solution for x &gt; 0, implies Im ω &lt; 0. Next, one must check that U -1 k U c admits eigenvalues with Im ω &lt; 0.

First, notice that Eq. (A1) is formally identical to the constraint on U c derived in the previous section, requiring

that there are no solutions with purely real ω . Second, using that U -1 k U c is a product of two Hermitian matrices, it is easy to prove that its eigenvalues come in complexconjugate pairs, i.e., simply take the complex conjugate of Eq. (A2) to find

$$0 & = \det \left ( U _ { c } - \omega U _ { k } \right ) = \det \left ( U _ { c } - \omega U _ { k } \right ) ^ { * } \quad \text{than} \\ & = \det \left [ ( U _ { c } - \omega U _ { k } ) ^ { \dagger } \right ] = \det \left ( U _ { c } - \omega ^ { * } U _ { k } \right ). \quad ( A 3 )$$

Therefore, both ω and ω ∗ are eigenvalues of U -1 k U c , and ∗ labels complex conjugation. This guarantees that there are always exactly N/ 2 solutions with Im ω &lt; 0.

Finally, we check that all solutions of this boundary condition obey current conservation, meaning vanishing normal current for any solution of a hard-wall boundary condition. The current operator normal to the boundary is given by

$$J = \frac { \partial H } { \partial ( \hbar { k } ) } = v _ { F } U _ { k }, \quad \quad ( \text{A} 4 ) \quad \text{a} _ { \text{tho} } ^ { \quad ] }$$

and the expectation value of the current is

$$F _ { 0 } ^ { \dagger } J F _ { 0 } & = v _ { F } F _ { 0 } ^ { \dagger } U _ { k } F _ { 0 } & \text{sum} \\ & = \omega ^ { - 1 } v _ { F } F _ { 0 } ^ { \dagger } U _ { k } ( U _ { k } ^ { - 1 } U _ { c } ) F _ { 0 } & \text{pen} \\ & = \omega ^ { - 1 } v _ { F } F _ { 0 } ^ { \dagger } U _ { c } F _ { 0 }, & \text{(A5)} & \text{ary}$$

where we used that F 0 is a solution of the boundary condition. The left-hand side is the expectation value of a Hermitian operator and is real, while the right hand side is a nonzero, not real number times another real expectation value. This is only possible if F U F † 0 k 0 = F U F † 0 c 0 = 0. This shows that the boundary condition is current conserving.

## 3. Boundary conditions for the k 2 model

For the full model Hamiltonian from Eq. (1), we expect that the hard-wall boundary condition becomes F (0) = 0 at the hard-wall set at x = 0. This is only possible if the bulk solutions in the inner region [ x &lt; 0, thus V ( x ) = 0] are linearly dependent, thus allowing for a linear combination that vanishes at the boundary.

To verify this condition, consider a plane-wave ansatz F x &lt; ( 0) ∼ F e k ikx , with Im k ≤ 0 to avoid divergent states for x →-∞ . The Schr¨dinger equation for o x &lt; 0 becomes

$$\left ( \hbar { v } _ { F } U _ { k } k + m U _ { w } k ^ { 2 } - E \right ) F _ { k } = 0. \quad \ ( \text{A6} ) \quad \text{the}$$

Nontrivial solutions require det[ S k, E ( )] = 0, with S k, E ( ) = glyph[planckover2pi1] v F U k k + mU k w 2 -E , yielding a 2 N th -order polynomial equation for k , where N is the size of the matrix S k, E ( ). For a fixed E , this provides 2 N complex roots k n , which are either purely real, or come in complex-conjugate pairs, as shown by a similar argument to Eq. (A3). Together with the constraint on U w this guarantees that there are N bulk states at a given E with purely real k n , and a total of N + N/ 2 roots with Im k n ≤ 0. Since the nullity of S k, E ( ) cannot be larger than N , this implies that the set of zero eigenvectors F k n with Im k n &lt; 0 is linearly dependent, thus validating the criteria stated above.

## 4. Relation between the wave functions of the k -linear and k 2 models

To further validate the equivalence between the k -linear and k 2 models, we show that their envelope functions F x ( ) are closely related in the limit of small energies E ≈ 0 and momentum k ≈ 0, for which the linear spectrum dominates.

First, consider the linear model H from Eq. (6), and a plane-wave ansatz F x ( ) = F e n ikx . Since [ H,U k ] = 0, they share a common set of eigenmodes F n , which does not depend on k . The set F n spans the full internal Hilbert space, thus allowing for a linear combination that satisfies the boundary condition from Eq. (6) (see Appendix A 2).

Second, for the k 2 model from Eq. (7), the boundary condition changes to F (0) = 0 (see Appendix A 3). However, in the smallk limit we can neglect the k 2 term, hence, the smallk eigenmodes are approximately the same set F n from the k -linear case. Since this set is a complete basis, the boundary condition F (0) = 0 can only be set by the trivially zero linear combination. Therefore, this boundary condition requires the modes from higher momentum states. To obtain these, consider the small energy limit ( E → 0) on the Schr¨dinger equao tion for the k 2 model, and simplifying k in the equation to eliminate the smallk solutions, we get

$$( \hbar { v } _ { F } U _ { k } + m U _ { w } k ) F _ { k } = 0. \text{ \quad \ \ } ( A 7 )$$

As shown in Appendix A 3, this equation is guaranteed to provide evanescent solutions that decay toward the interior of the sample, where k e = glyph[planckover2pi1] v F / mω ( ) with Im k e &gt; 0, and ω is a complex eigenvalue of U -1 k U w with Im ω &gt; 0. This equation is formally identical to the boundary condition of the k -linear model (4). Since the smallk modes F n form a complete set, the addition of these F k e evanescent modes makes the set linearly dependent, allowing a vanishing combination at the boundary. For small m k , e is large and the evanescent solutions decay quickly, while the propagating part of the solution with small k satisfies the same boundary condition as the linear model.

The analysis above shows that the envelope functions of the two models match far from the boundaries. Near the boundary, the k 2 model has an additional localized contribution that decays quickly away from the boundary in the limit of small m [53].

## Appendix B: Full k 2 model for graphene

The full monolayer graphene honeycomb lattice is invariant under the ( P /mmm 6 ) symmorphic space group. Eliminating the trivial Bloch translations group T B , its factor point group is ( P /mmm 6 )/ T B ∼ D 6 h . On the other hand, once the monolayer is cut to form the nanoribbons, the sixfold rotation symmetry is broken, such that the zigzag and armchair nanoribbons may transform as either the Pmmm or Pmma space groups, depending on their widths.

Armchair nanoribbons with odd N A transform as the symmorphic Pmmm space group, while for even N A , it transforms as the nonsymmorphic Pmma . Nevertheless, since the monolayer Dirac cones are projected into ¯ Γ for the armchair confinement, both have the same factor group under trivial Bloch translations, thus, Pmmm/T B ≡ Pmma/T B ∼ D 2 h . Similarly, zigzag nanoribbons with even (odd) N D belong to the Pmmm ( Pmma ) space group. In this case, the K and K ′ monolayer valleys fall into the ¯ K and ¯ K ′ under the confinement projection. Consequently, the Pmma nonsymmorphic symmetries yield an extra phase into the representation matrices. Fortunately, this phase matches that of a single Bloch translation, thus, it also follows Pmmm /T B ≡ Pmma /T B ∼ D 2 h . Therefore, hereafter it is sufficient to analyze the symmorphic lattices and the point group D 2 h . This allows us to build a single model for both armchair and zigzag nanoribbons, considering a basis set that contains both K and K ′ basis functions.

The D 2 h point group can be generated by its mirror operations M x , M y , and M z . Here, M x reflects x →-x , and similarly for M y and M z . Under the basis vector r = ( x, y, z ), the coordinates representation from the O (3) group is given by the matrices

$$D ^ { r } ( M _ { x } ) = \text{diag} ( - 1, + 1, + 1 ), \quad \text{ (B1)} \quad \text{of $i$}$$

$$D ^ { r } ( M _ { y } ) = \text{diag} ( + 1, - 1, + 1 ), \quad \ ( B 2 ) \quad \frac { \delta \ s } { \text{mein} }$$

$$D ^ { r } ( M _ { z } ) = \text{diag} ( + 1, + 1, - 1 ), \quad \ ( B 3 ) \quad \underset { \text{res} } { \ u s m }$$

where diag( · · · ) labels a diagonal matrix with elements given by its arguments.

To obtain the Hilbert space representation H we consider the basis functions { ϕ A ( r ) , ϕ B ( r ) , ϕ A ′ ( r ) , ϕ B ′ ( r ) } shown in Fig. 1(a). These are built from p z orbitals of the carbon atoms centered at the A or B lattice sites and Bloch phases related to the K or K ′ valleys. Namely, the representation matrices for the D 2 h generators are

$$D ^ { \mathcal { H } } ( M _ { x } ) = - \tau _ { 0 } \otimes \sigma _ { y }, \text{ \quad \ \ } \text{(B4)}$$

$$D ^ { \mathcal { H } } ( M _ { y } ) = + \tau _ { x } \otimes \sigma _ { 0 }, \text{ \quad \ \ } \text{(B5)}$$

$$D ^ { \mathcal { H } } ( M _ { z } ) = - \tau _ { 0 } \otimes \sigma _ { 0 }, \text{ \quad \ \ } \text{(B6)}$$

where σ = ( σ , σ 0 x , σ y , σ z ) are Pauli matrices acting on the A B / lattices subspace, and τ = ( τ 0 , τ x , τ y , τ z ) acts on the K/K ′ valley subspace.

Additionally, we consider that the system is chiral C and time-reversal T invariant. Since our graphene model is spinless, T = K is simply the complex conjugation. The C symmetry labels the sublattices. Under the r -representation D r ( T ) = D r ( C ) = 1 , while on k space D k ( T ) = -1 and D k ( C ) = 1 . Within the H representation,

$$D ^ { \mathcal { H } } ( \mathcal { T } ) = \tau _ { x } \otimes \sigma _ { z } \mathcal { K }, \text{ \quad \ \ } \text{(B7)}$$

$$D ^ { \mathcal { H } } ( \mathcal { C } ) = \tau _ { 0 } \otimes \sigma _ { z }.$$

To obtain the effective Hamiltonian H , we consider the method of invariants [49]. Thus, we seek the most general form of H ≡ H ( k ) as an expansion in powers of k = ( k , k x y ) that is invariant, i.e. , [ H, S ] = 0 for all S symmetries above. This can be easily implemented in using the QSYMM Python's package [42]. Splitting the resulting terms as H = H 0 + H A + H Z + H ft , we obtain

$$H _ { 0 } = \hbar { v } _ { x } \begin{pmatrix} \sigma _ { x } & 0 \\ 0 & \sigma _ { x } \end{pmatrix} k _ { x } + \hbar { v } _ { y } [ \begin{pmatrix} \sigma _ { y } & 0 \\ 0 & - \sigma _ { y } \end{pmatrix} k _ { y } + \begin{pmatrix} \sigma _ { y } & 0 \\ 0 & \sigma _ { y } \end{pmatrix} \Delta _ { K } ],$$

$$H _ { A } = \frac { m _ { A A } } { 2 } \begin{pmatrix} \sigma _ { y } & 0 \\ 0 & \sigma _ { y } \end{pmatrix} k _ { y } ^ { 2 } + \frac { m _ { A A } } { 2 } \begin{pmatrix} 0 & \sigma _ { y } \\ \sigma _ { y } & 0 \end{pmatrix} k _ { y } ^ { 2 } \longrightarrow H _ { A } = \frac { m _ { A } } { 2 } \begin{pmatrix} \sigma _ { y } \cot \theta & \sigma _ { y } \csc \theta \\ \sigma _ { y } \csc \theta & \sigma _ { y } \cot \theta \end{pmatrix} k _ { y } ^ { 2 },$$

$$H _ { Z } = \frac { m _ { Z 1 } } { 2 } \begin{pmatrix} \sigma _ { y } & 0 \\ 0 & \sigma _ { y } \end{pmatrix} k _ { x } ^ { 2 } + \frac { m _ { Z 2 } } { 2 } \begin{pmatrix} 0 & \sigma _ { y } \\ \sigma _ { y } & 0 \end{pmatrix} k _ { x } ^ { 2 } \longrightarrow H _ { Z } = \frac { m _ { Z } } { 2 } \begin{pmatrix} \eta \sigma _ { y } & ( 1 - \eta ) \sigma _ { y } \\ ( 1 - \eta ) \sigma _ { y } & \eta \sigma _ { y } \end{pmatrix} k _ { x } ^ { 2 },$$

$$H _ { \mathrm f l } = \Delta \begin{pmatrix} 0 & \sigma _ { y } \\ \sigma _ { y } & 0 \end{pmatrix} + \hbar { \mu } \begin{pmatrix} 0 & \sigma _ { x } \\ \sigma _ { x } & 0 \end{pmatrix} k _ { x } + m _ { x y } \begin{pmatrix} \sigma _ { x } & 0 \\ 0 & - \sigma _ { x } \end{pmatrix} k _ { x } k _ { y }.$$

Here, H 0 represents a minimal Dirac-like model with anisotropic Fermi velocities v x and v y , with the mono- layer projected cones at k y = ± ∆ K = ± ¯ K for the zigzag confinement, and ∆ K = 0 for armchair. The H A and H Z

show the most general k 2 x and k 2 y terms on the left hand side. On the right hand side, the second form of H A and H Z are written as the U c ≡ U w confinement or Wilson's matrices for the minimal armchair and zigzag models, respectively. Additionally, H ft contains fine-tuning terms that were neglected on the minimal models [Eqs. (16) and (17)].

We have fit the full k 2 model to the DFT data, trying to use a minimal set of finite parameters (see Table I). The results for the N A = 48 armchair ribbon was shown in the main text, Fig. 4. Here we show equivalent results for N A = 49 in Figs. 6(a1)-(e1), for the metallic case N A = 50 in Figs. 6(a2)-(e2), and for a zigzag ribbon with N D = 48 in Figs. 6(a3)-(e3). In all cases, the agreement between the models and the DFT data at low energies is satisfactory. At higher energies, the k 2 model provides slightly better results. Since all models are based on low-energy expansions, the agreement with DFT must improve for wider ribbons. However, the main conclusion to be extracted from the comparison is that the k 2 model provides a better approach for numerical simulations in comparison to the k -linear model, which

- [1] G. Dresselhaus and M. S. Dresselhaus, Spin-Orbit Interaction in Graphite, Phys. Rev. 140 , A401 (1965).
- [2] C. L. Kane and E. J. Mele, Quantum Spin Hall Effect in Graphene, Phys. Rev. Lett. 95 , 226801 (2005).
- [3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81 , 109 (2009).
- [4] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Quantum spin Hall effect and topological phase transition in HgTe quantum wells, Science 314 , 1757 (2006).
- [5] B. A. Bernevig and S.-C. Zhang, Quantum spin hall effect, Phys. Rev. Lett. 96 , 106802 (2006).
- [6] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82 , 3045 (2010).
- [7] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83 , 1057 (2011).
- [8] S.-Q. Shen, Topological insulators: Dirac Equation in Condensed Matter , Springer Series in Solid-State Sciences, Vol. 174 (Springer-Verlag, Berlin, Heidelberg, 2012).
- [9] B. A. Bernevig and T. L. Hughes, Topological insulators and topological superconductors (Princeton University Press, 2013).
- [10] L. Fu, Topological crystalline insulators, Phys. Rev. Lett. 106 , 106802 (2011).
- [11] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Topological crystalline insulators in the SnTe material class, Nat. Commun. 3 , 982 (2012).
- [12] R.-J. Slager, A. Mesaros, V. Juriˇi´, and J. Zaanen, The c c space group classification of topological band-insulators, Nat. Phys. 9 , 98 (2013).
- [13] E. O. Wrasse and T. M. Schmidt, Prediction of TwoDimensional Topological Crystalline Insulator in PbSe Monolayer, Nano Letters 14 , 5717-5720 (2014).

requires complex handling [31, 32].

Table I. Parameters used for the k -linear, k 2 , and tightbinding models for the armchair nanoribbons.

|                           | Armchair        | Armchair        | Armchair        | Zigzag          | Zigzag          |
|---------------------------|-----------------|-----------------|-----------------|-----------------|-----------------|
| N A                       | 48              | 49              | 50              | N D             | 48              |
| k -linear model           | k -linear model | k -linear model | k -linear model | k -linear model | k -linear model |
| v F [10 3 nm/ps]          | 0.8             | 0.8             | 0.8             | v F             | 0.8             |
| ∆ θ [deg]                 | 21.7            | 23.1            | 19.8            | -               | -               |
| k 2 model                 | k 2 model       | k 2 model       | k 2 model       | k 2 model       | k 2 model       |
| v x = v y [10 3 nm/ps]    | 0.8             | 0.8             | 0.8             | v x = v y       | 0.8             |
| m A [meV nm 2 ]           | 50              | 50              | 50              | m Z             | 50              |
| ∆ θ [deg]                 | 21.7            | 23.1            | 19.8            | η               | 0               |
| m Z 1 = m Z 2 [meV nm 2 ] | 0               | 0               | 0               | m A 1 = m A 2   | 0               |
| µ [10 3 nm/ps]            | 0.1             | -0.05           | 0               | µ               | 0               |
| ∆ [meV]                   | 0               | 0               | 0               | ∆               | 0               |
| m xy [meV nm 2 ]          | -50             | -50             | -50             | m xy            | 0               |
| Tight binding             | Tight binding   | Tight binding   | Tight binding   | Tight binding   | Tight binding   |
| Hopping [eV]              | 2.4             | 2.4             | 2.4             | Hopping         | 2.4             |

- [14] Y. Ando and L. Fu, Topological Crystalline Insulators and Topological Superconductors: From Concepts to Materials, Annu. Rev. Condens. Matter Phys. 6 , 361 (2015).
- [15] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Topological classification of crystalline insulators through band structure combinatorics, Phys. Rev. X 7 , 041069 (2017).
- [16] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Phys. Rev. B 96 , 245115 (2017).
- [17] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Quantized electric multipole insulators, Science 357 , 61 (2017).
- [18] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Reflection-symmetric second-order topological insulators and superconductors, Phys. Rev. Lett. 119 , 246401 (2017).
- [19] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Higherorder topological insulators, Sci. Adv. 4 , eaat0346 (2018).
- [20] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and fermi-arc surface states in the electronic structure of pyrochlore iridates, Phys. Rev. B 83 , 205101 (2011).
- [21] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Discovery of a Weyl fermion semimetal and topological Fermi arcs, Science 349 , 613 (2015).
- [22] Z. Yue, G. Xue, J. Liu, Y. Wang, and M. Gu, Nanometric holograms based on a topological insulator material, Nat.

Figure 6. Comparison between the DFT band structures and (a) the k -linear, (b) tight-binding, and (c) k 2 models for armchair nanoribbons with (a1)-(c1) N A = 49, (a2)-(c2) N A = 50, and (a3)-(c3) zigzag nanoribbons with N D = 48. The envelope densities (colored symbols) match well the DFT data for (d) the first and (e) second conduction bands at k x = 0 in all cases.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00145v2-figure-79d87e7a-3250-4bf9-a795-9ebc1f46f045.jpg)

## Commun. 8 , 15354 (2017).

- [23] E. Goi, Z. Yue, B. P. Cumming, and M. Gu, Observation of Type I Photonic Weyl Points in Optical Frequencies, Laser Photon. Rev. 12 , 1700271 (2018).
- Energy levels of triangular and hexagonal graphene quantum dots: A comparative study between the tightbinding and Dirac equation approach, Phys. Rev. B 84 , 245403 (2011).
- [24] C. Yue, S. Jiang, H. Zhu, L. Chen, Q. Sun, and D. W. Zhang, Device applications of synthetic topological insulator nanostructures, Electronics 7 , 225 (2018).
- [25] Y. Fan and K. L. Wang, Spintronics based on topological insulators, SPIN 6 , 1640001 (2016).
- [26] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10 , 2445 (1974).
- [27] J. B. Kogut and L. Susskind, Hamiltonian formulation of Wilson's lattice gauge theories, Phys. Rev. D 11 , 395 (1975).
- [28] H. B. Nielsen and M. Ninomiya, A no-go theorem for regularizing chiral fermions, Phys. Lett. B 105 , 219 (1981).
- [29] H. B. Nielsen and M. Ninomiya, Absence of neutrinos on a lattice: (I). Proof by homotopy theory, Nucl. Phys. B 185 , 20 (1981).
- [30] D. R. Candido, M. Kharitonov, J. C. Egues, and E. M. Hankiewicz, Paradoxical extension of the edge states across the topological phase transition due to emergent approximate chiral symmetry in a quantum anomalous Hall system, Phys. Rev. B 98 , 161111(R) (2018).
- [31] J. Tworzydglyph[suppress]o, l C. W. Groth, and C. W. J. Beenakker, Finite difference method for transport properties of massless dirac fermions, Phys. Rev. B 78 , 235438 (2008).
- [32] A. R. Hern´ndez and C. H. Lewenkopf, Finite-difference a method for transport of two-dimensional massless Dirac fermions in a ribbon geometry, Phys. Rev. B 86 , 155439 (2012).
- [33] Y.-F. Zhou, H. Jiang, X. C. Xie, and Q.-F. Sun, Twodimensional lattice model for the surface states of topological insulators, Phys. Rev. B 95 , 245137 (2017).
- [34] M. Gruji´, M. Zarenia, A. Chaves, M. Tadi´, G. A. Farias, c c and F. M. Peeters, Electronic and optical properties of a circular graphene quantum dot in a magnetic field: Influence of the boundary conditions, Phys. Rev. B 84 , 205441 (2011).
- [35] M. Zarenia, A. Chaves, G. A. Farias, and F. M. Peeters,
- [36] M. Mirzakhani, M. Zarenia, D. R. da Costa, S. A. Ketabi, and F. M. Peeters, Energy levels of ABC-stacked trilayer graphene quantum dots with infinite-mass boundary conditions, Phys. Rev. B 94 , 165423 (2016).
- [37] E. McCann and V. I. Fal'ko, Symmetry of boundary conditions of the Dirac equation for electrons in carbon nanotubes, J. Phys. Condens. Matter 16 , 2371 (2004).
- [38] A. R. Akhmerov and C. W. J. Beenakker, Boundary conditions for Dirac fermions on a terminated honeycomb lattice, Phys. Rev. B 77 , 085423 (2008).
- [39] M. V. Berry and R. J. Mondragon, Neutrino billiards: time-reversal symmetry-breaking without magnetic fields, Proc. Royal Soc. Lond. 412 , 53 (1987).
- [40] L. Brey and H. A. Fertig, Electronic states of graphene nanoribbons studied with the Dirac equation, Phys. Rev. B 73 , 235411 (2006).
- [41] B. Messias de Resende, F. C. de Lima, R. H. Miwa, E. Vernek, and G. J. Ferreira, Confinement and fermion doubling problem in Dirac-like Hamiltonians, Phys. Rev. B 96 , 161113(R) (2017).
- [42] D. Varjas, T. O. Rosdahl, and A. R. Akhmerov, Qsymm: algorithmic symmetry finding and symmetric Hamiltonian generation, New J. of Phys. 20 , 093026 (2018).
- [43] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, Kwant: a software package for quantum transport, New J. Phys. 16 , 063065 (2014).
- [44] Supplemental material is available at the arXiv repository, which includes (i) the DFT/VASP raw data; (ii) python/qsymm [42] implementation of the symmetry constraints to obtain graphene's effective models; (iii) python/kwant [43] implementation of graphene's tightbinding model; and (iv) numerical (python) implementation of the effective models using finite differences, and the codes used to generate the figures.
- [45] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett.

- 77 , 3865 (1996).
- [46] P. E. Bl¨chl, Projector augmented-wave method, Phys. o Rev. B 50 , 17953 (1994).
- [47] G. Kresse and J. Furthm¨ller, Efficiency of ab-initio total u energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 , 15 (1996).
- [48] G. Kresse and J. Furthm¨ller, Efficient iterative schemes u for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 , 11169 (1996).
- [49] R. Winkler, Spin-orbit coupling effects in twodimensional electron and hole systems , Springer Tracts in Modern Physics, Vol. 191 (Springer-Verlag, Berlin, Heidelberg, 2003).
- [50] V. Alonso, S. De Vincenzo, and L. Mondino, On the boundary conditions for the dirac equation, Eur. J. Phys. 18 , 315 (1997).
- [51] V. Alonso and S. De Vincenzo, General boundary conditions for a dirac particle in a box and their nonrelativistic limits, J. Phys. A: : Math. Gen. 30 , 8573 (1997).
- [52] G. J. Ferreira and D. Loss, Magnetically Defined Qubits on 3D Topological Insulators, Phys. Rev. Lett. 111 , 106802 (2013).
- [53] In the continuum limit δ x → 0. Additionally, requiring that the linear spectrum to dominate over the full energy range δ ε →∞ , the constraints in Eq. (8) yields m ≡ 0.
- [54] Y.-W. Son, M. L. Cohen, and S. G. Louie, Energy gaps in graphene nanoribbons, Phys. Rev. Lett. 97 , 216803 (2006).
- [55] D. P. DiVincenzo and E. J. Mele, Self-consistent effective-
- mass theory for intralayer screening in graphite intercalation compounds, Phys. Rev. B 29 , 1685 (1984).
- [56] G. Bastard, Wave mechanics applied to semiconductor heterostructures , Monographies de physique (Les ´ ditions de Physique, 1988). E
- [57] T. Ando, Theory of electronic states and transport in carbon nanotubes, J. Phys. Soc. Jpn. 74 , 777 (2005).
- [58] A. L. Ara´jo, E. O. Wrasse, G. J. Ferreira, and T. M. u Schmidt, Topological nonsymmorphic ribbons out of symmorphic bulk, Phys. Rev. B 93 , 161101(R) (2016).
- [59] K. Wakabayashi, K. i. Sasaki, T. Nakanishi, and T. Enoki, Electronic states of graphene nanoribbons and analytical solutions, Sci. Technol. Adv. Mater 11 , 054504 (2010).
- [60] J. Liu, X. Qian, and L. Fu, Crystal Field Effect Induced Topological Crystalline Insulators In Monolayer IV-VI Semiconductors, Nano Letters 15 , 2657 (2015).
- [61] C. Niu, P. M. Buhl, G. Bihlmayer, D. Wortmann, S. Bl¨gel, and Y. Mokrousov, Topological crystalline inu sulator and quantum anomalous Hall states in IV-VIbased monolayers and their quantum wells, Phys. Rev. B 91 , 201401(R) (2015).
- [62] In Ref. [58], we have considered a Brey-Fertig type of boundary condition applied at the different atomic terminations of each sublattice (Pb/Se). This gives the correct results, as it breaks chiral symmetry properly for ribbons A and B. Within our approach, the same result is much simpler to obtain with a finite ρ .
A. L. Araújo R. P. Maciel R. G. F. Dornelas D. Varjas G. J. Ferreira
Interplay between boundary conditions and Wilson's mass in Dirac-like Hamiltonians
Dirac-like Hamiltonians, linear in momentum $k$, describe the low-energy
physics of a large set of novel materials, including graphene, topological
insulators, and Weyl fermions. We show here that the inclusion of a minimal
$k^2$ Wilson's mass correction improves the models and allows for systematic
derivations of appropriate boundary conditions for the envelope functions on
finite systems. Considering only Wilson's masses allowed by symmetry, we show
that the $k^2$ corrections are equivalent to Berry-Mondragon's discontinuous
boundary conditions. This allows for simple numerical implementations of
regularized Dirac models on a lattice, while properly accounting for the
desired boundary condition. We apply our results on graphene nanoribbons
(zigzag and armchair), and on a PbSe monolayer (topological crystalline
insulator). For graphene, we find generalized Brey-Fertig boundary conditions,
which correctly describe the small gap seen on \textit{ab initio} data for the
metallic armchair nanoribbon. On PbSe, we show how our approach can be used to
find spin-orbital-coupled boundary conditions. Overall, our discussions are set
on a generic model that can be easily generalized for any Dirac-like
Hamiltonian.