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

TL;DR

This paper demonstrates how including a minimal $k^2$ Wilson's mass correction in Dirac-like Hamiltonians enables systematic derivation of boundary conditions, improving modeling of finite systems like graphene nanoribbons and topological insulators.

Contribution

It introduces a method to derive boundary conditions from Wilson's mass corrections, connecting them to Berry-Mondragon boundary conditions for Dirac models.

Findings

Wilson's mass corrections are equivalent to Berry-Mondragon boundary conditions.

The approach accurately models boundary effects in graphene nanoribbons.

The method applies to topological crystalline insulators like PbSe.

Abstract

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.

Tables

Table 1. Table 1: Parameters used for the k 𝑘 k -linear, k 2 superscript 𝑘 2 k^{2} , and tight-binding models for the armchair nanoribbons.

	Armchair			Zigzag
$N_A$	48	49	50	$N_D$	48
$k$-linear model
$v_F$ [${10}^3$ nm/ps]	0.8	0.8	0.8	$v_F$	0.8
$\Delta \theta$ [deg]	21.7	23.1	19.8	-	-
$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 nm2]	50	50	50	$m_Z$	50
$\Delta \theta$ [deg]	21.7	23.1	19.8	$\eta$	0
$m_{Z1}=m_{Z2}$ [meV nm2]	0	0	0	$m_{A1}=m_{A2}$	0
$\mu$ [${10}^3$ nm/ps]	0.1	-0.05	0	$\mu$	0
$\Delta$ [meV]	0	0	0	$\Delta$	0
$m_{xy}$ [meV nm2]	-50	-50	-50	$m_{xy}$	0
Tight binding
Hopping [eV]	2.4	2.4	2.4	Hopping	2.4

Equations

$$H=\hslash v_F{U}_{k}k+m{U}_w{k}^2+U_{c}V(x),$$

$$[U_w,D^{\psi }({\mathcal{S}}_{\pm })]$$

$$[U_k,D^{\psi }({\mathcal{S}}_{+})]$$

$$\left(i{U}_k^{-1}{U}_c\right)F_0$$

$$M{F}_0$$

$$H=\hslash v_F{U}_{k}k,\text{ with }(U_k\mp i{U}_c)F(\pm L)$$

$$H=\hslash v_F{U}_{k}k+m{U}_w{k}^2,\text{ with }F(\pm L)$$

$$\frac{{\delta }_x^2{\delta }_{\epsilon }}{2}\le |m|\le \frac{(\hslash v_F{)}^2}{{\delta }_{\epsilon }},$$

$$\displaystyle\psi(\bm{r})=\sum_{\mu}\Big{[}f_{\mu}(\bm{r})e^{i\bm{q}\cdot\bm{r}}\varphi_{\mu}(\bm{r})+f_{\mu^{\prime}}(\bm{r})e^{i\bm{q}^{\prime}\cdot\bm{r}}\varphi_{\mu^{\prime}}(\bm{r})\Big{]},$$

$$H_G(\mathbf{k})$$

$$h(\mathbf{q})$$

$$\displaystyle f_{A}\Big{(}-\frac{W}{2}\Big{)}$$

$$\displaystyle f_{B}\Big{(}+\frac{W}{2}\Big{)}$$

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

$$U_k^Z$$

$$U_c^Z(\eta )$$

$$U_k^A$$

$$U_c^A(\theta )$$

$$H_Z$$

$$H_A$$

$$E_{3p}$$

$$E_{3p+1}$$

$$E_{3p+2}$$

$$H=h(\bm{k}^{+})\oplus h^{*}(\bm{k}^{-})\\ +\Big{[}\dfrac{m_{A1}}{2}U_{0y}+\dfrac{m_{A2}}{2}U_{xy}\Big{]}k_{y}^{2}+\Big{[}\dfrac{m_{Z1}}{2}U_{0y}+\dfrac{m_{Z2}}{2}U_{xy}\Big{]}k_{x}^{2}\\ +\Delta U_{xy}+\hbar\mu U_{xx}k_{x}+m_{xy}U_{zx}k_{x}k_{y}.$$

$$H=\Delta U_{z0}+\alpha (U_{xx}{k}_y-U_{xy}{k}_x)+\frac{m}{2}{U}_c(\rho ,\theta )k_y^2+{\Delta }_C{U}_{yx},$$

$$U_{c}(\rho,\theta)=\cos\theta\Big{[}\cosh(\rho)U_{z0}+\sinh(\rho)U_{00}\Big{]}+\sin(\theta)U_{yx}.$$

$$\displaystyle f_{x,+\sigma}\Big{(}\pm\frac{L}{2}\Big{)}$$

$$\displaystyle\Big{(}\dfrac{i\hbar v_{F}}{V_{0}\lambda}U_{k}+U_{c}\Big{)}F_{0}=0.$$

$$\displaystyle\det\Big(U_{k}^{-1}U_{c}-\omega\Big{missing})=0,$$

$$0$$

$$\displaystyle=\det[(U_{c}-\omega U_{k})^{\dagger}\Big{]}=\det\Big(U_{c}-\omega^{*}U_{k}\Big{missing}).$$

$$J=\frac{\partial H}{\partial (\hslash k)}=v_F{U}_k,$$

$$F_0^{†}J{F}_0$$

$$={\omega }^{-1}{v}_F{F}_0^{†}{U}_k(U_k^{-1}{U}_c)F_0$$

$$={\omega }^{-1}{v}_F{F}_0^{†}{U}_c{F}_0,$$

$$\displaystyle\Big{(}\hbar v_{F}U_{k}k+mU_{w}k^{2}-E\Big{)}F_{k}=0.$$

Full text

I Introduction

The Dirac-like Hamiltonians play an ubiquitous role in novel materials, ranging from graphene [1, 2, 3] to topological insulators (TI) [4, 5, 6, 7, 8, 9], its crystalline [10, 11, 12, 13, 14, 15] and higher-order [16, 17, 18, 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, 23, 24, 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, 27, 28, 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, 32, 33]. For finite systems (e.g., nanoribbons), the Dirac models allow for a variety of possible non-trivial boundary conditions [34, 35, 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?

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, 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^{\prime}$ 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 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 finite-differences implementations [31, 32, 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/Å.

II Generic model

Effective models can be obtained from symmetry constraints imposed by method of invariants [49], which is equivalent to a $k\cdot p$ envelope function approach, yielding a matrix expansion of the Hamiltonian $H\equiv H(\bm{k})$ in powers of the momentum $\bm{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 [39, 37, 40, 38, 41] (see Sec. II.2). 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 one-dimensional system given by the generic Hamiltonian

[TABLE]

which is defined along $x$, and $k=-i\partial_{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 soft-wall confining potential given by a symmetric profile $V(x)=V_{0}[1-\Theta(x+L)+\Theta(x-L)]$, where $\Theta(x)$ is the Heaviside step function, and $V_{0}$ is the intensity. This profile defines the physical system within the inner region $|x|<L$, while on the outer region ($|x|\geq L$) it opens a gap $2|V_{0}|$. Later, we will consider the hard-wall limit $|V_{0}|\rightarrow\infty$, which excludes the outer region from the physical domain. The $U_{k},U_{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.

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})\neq 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)\neq 0$. Within this region, we demand that Eq. (1) does not have any propagating modes at zero energy. To guarantee this, we demand that $U_{k}^{-1}U_{c}$ has no purely real eigenvalues, which implies $\det(U_{c})\neq 0$. The anticommutation condition $\left\{U_{k},U_{c}\right\}=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\neq 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).

II.1 Symmetry constraints: $U_{w}\equiv U_{c}$

Let us consider that our generic system, i.e., $H$ from Eq. (1), is invariant under a symmetry group $\mathcal{G}$, which is composed by two types of symmetry operations: $\mathcal{S}_{+}$ and $\mathcal{S}_{-}$. The $\mathcal{S}_{+}$ operators leave $x$ invariant, while $\mathcal{S}_{-}$ takes $x\rightarrow-x$. Therefore, it follows the transformations $\mathcal{S}_{\pm}x\mathcal{S}_{\pm}^{-1}=\pm x$, $\mathcal{S}_{\pm}k\mathcal{S}_{\pm}^{-1}=\pm k$, $\mathcal{S}_{\pm}V(x)\mathcal{S}_{\pm}^{-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 $\mathcal{G}$, (i.e., $\mathcal{S}_{\pm}H\mathcal{S}_{\pm}^{-1}=H$), one obtains the symmetry constraints for the matrices of $H$:

[TABLE]

which defines the symmetry-allowed matrices $U_{k}$, $U_{w}$, and $U_{c}$. Here $[\cdot,\cdot]$ and $\{\cdot,\cdot\}$ are the commutator and anti-commutator operations, and $D^{\psi}(\mathcal{S}_{\pm})$ are the matrix representations of $\mathcal{S}_{\pm}$ in the Hilbert space. Since $U_{w}$ and $U_{c}$ satisfy the same constraints, it follows that they are equivalent ($U_{w}\equiv 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.

II.2 Hard-wall boundary conditions

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 equation $HF(x)=EF(x)$ has order 2 if $m\neq 0$, or order 1 if $m=0$. Here, $F(x)$ is an envelope spinor function [49]. In all cases, the energy $E=\bra{F}H\ket{F}$ must be bounded and well defined. Consider $H$ from Eq. (1), with a simplified single boundary profile $V(x)=V_{0}\Theta(x)$. On the outer region $x>0$ the gap $|2V_{0}|$ yields evanescent solutions at low energies, i.e., $F(x>0)\sim F_{0}e^{-x/\lambda}$, where $F_{0}$ is the spinor at $x=0$. If the $k$-linear term dominates the low-energy band structure, the penetration length is $\lambda\propto\hbar v_{F}/|V_{0}|$. In the hard-wall limit $|V_{0}|\rightarrow\infty$ and $\lambda\rightarrow 0$, thus near the interface $x\approx 0$ we can write $F(x)\approx F_{0}[1-\Theta(x)]$. Considering only the $x\approx 0$ range on the integrals in $E=\bra{F}H\ket{F}$, it can be shown that the contributions from the $k$-linear and potential $V(x)$ terms are finite, while $m\bra{F}U_{w}k^{2}\ket{F}\approx-mF_{0}^{\dagger}U_{w}F_{0}\delta(0)$ is ill defined due to the $\delta(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, 51. Later, it was applied to graphene [37, 38] and topological insulators [52, 41, 30]. 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

[TABLE]

This is an eigenvalue equation for the matrix $iU_{k}^{-1}U_{c}$ with eigenvalue $\alpha=\hbar v_{F}/\left(\lambda V_{0}\right)$. The solutions with $\real\lambda>0$ ($\real\alpha>0$) describe states decaying outside the boundary. Defining $M=\alpha^{-1}\left(iU_{k}^{-1}U_{c}\right)$, the boundary condition becomes

[TABLE]

If $U_{c}$ and $U_{k}$ are chosen such that $M$ is traceless, unitary, and Hermitian, this reproduces the most general current-conserving boundary conditions [39, 37, 38]. 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 $\left(iU_{k}^{-1}U_{c}\right)$ has at least one pair of eigenvalues $\alpha=\pm 1$, equivalently $\det\left(U_{k}\mp iU_{c}\right)=0$. Using this, we cast the boundary conditions Eq. (5) for a finite system with $|x|<L$ as $(U_{k}\mp iU_{c})F(\pm L)=0$.

II.3 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:

[TABLE]

These two approaches are equivalent, due to the analytical properties of the boundary conditions discussed above, and the equivalence $U_{w}\equiv U_{c}$ shown in Sec. II.1. 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

[TABLE]

which ensures that the $k$-linear term dominates the low-energy spectrum within the $\delta_{\varepsilon}$ 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\rightarrow 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 [40, 3]. 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}=3p+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 ($\sim D_{2h}$, 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.3, we introduce the $k^{2}$ model for the nanoribbons. The $k\cdot p$ expansion for graphene considers basis functions given by its solutions at $K$ and $K^{\prime}$ valleys, i.e., $\varphi_{A}(\bm{r})$, $\varphi_{B}(\bm{r})$, $\varphi_{A^{\prime}}(\bm{r})$ and $\varphi_{B^{\prime}}(\bm{r})$, where $A$ and $B$ label the sublattices. These are illustrated in Fig. 1(a). Within the envelope function approximation [55, 56, 57], the expansion reads as

[TABLE]

where $f_{\mu}(\bm{r})$ and $f_{\mu^{\prime}}(\bm{r})$ are the envelope functions, $\mu=\{A,B\}$ label the sublattices, $\bm{q}=\bm{k}-\bm{K}$ and $\bm{q}^{\prime}=\bm{k}-\bm{K}^{\prime}$ are the deviations from the $K$ and $K^{\prime}$ valleys in $k$-space. Bloch theorem requires $\psi(\bm{r}+\bm{R})=e^{i\bm{k}\cdot\bm{R}}\psi(\bm{r})$, where $\bm{R}=n_{1}\bm{a}_{1}+n_{2}\bm{a}_{2}$ (with $n_{1}$, $n_{2}$ integers) is a Bravais translation of the monolayer. Since the Bloch phase in $\varphi_{\mu^{(\prime)}}(\bm{r}+\bm{R})=e^{i\bm{K}^{(\prime)}\cdot\bm{R}}\varphi_{\mu^{(\prime)}}(\bm{r})$ cancels out the opposite phase in the $q$-exponentials, Bloch theorem is satisfied for a periodic $f_{\mu^{(\prime)}}(\bm{r}+\bm{R})=f_{\mu^{(\prime)}}(\bm{r})$. Up to leading order in $\bm{k}$, the usual effective Dirac-like Hamiltonian is

[TABLE]

The $4\times 4$ $H_{G}(\bm{k})$ Hamiltonian acts on the envelope spinor $F(\bm{r})=[f_{A}(\bm{r}),f_{B}(\bm{r}),f_{A^{\prime}}(\bm{r}),f_{B^{\prime}}(\bm{r})]$. Notice that we write $H_{G}(\bm{k})$ in terms of the deviations $\bm{q}$ and $\bm{q}^{\prime}$, such that the Dirac cones occur at $\bm{k}\sim\bm{K}$ and $\bm{K}^{\prime}$. 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 $\psi(\bm{r})$ expansion above, but they will change the definitions of $\bm{q}$ and $\bm{q}^{\prime}$.

III.1 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.

III.1.1 Zigzag nanoribbons

To model the zigzag nanoribbons, we start from the bulk basis functions $\varphi_{\mu,\mu^{\prime}}(\bm{r})$, but with a modified $\psi(\bm{r})$ expansion. Namely, replace $\bm{q}^{(\prime)}\rightarrow\bm{k}-\bm{\bar{K}}^{(\prime)}$ in Eq. (9), where $\bar{K}^{(\prime)}$ are the $K$ and $K^{\prime}$ projections into the zigzag $k_{y}$ axis, as shown in Fig. 1(b). These $k$ projections also apply to $\bm{q}^{(\prime)}$ in $H_{G}(\bm{k})$ from Eq. (10). For simplicity, here we consider only this pristine form of $H_{G}(\bm{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=\pm W_{0}/2$) is composed only of atoms from the $A$ ($B$) sublattice. The length $W_{0}=(\frac{3}{2}N_{D}-1)a/\sqrt{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=\pm W/2$, which defines the effective length $W=W_{0}+2a/\sqrt{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 $\varphi_{B^{(\prime)}}(\bm{R}_{-})\approx 0$, where $\bm{R}_{\pm}=(\pm W/2,y)$. Imposing $\psi(\bm{R}_{-})=0$ in Eq. (9), we get $f_{A}(\bm{R}_{-})e^{-i\bm{\bar{K}}\cdot\bm{R}_{-}}\varphi_{A}(\bm{R}_{-})+f_{A^{\prime}}(\bm{R}_{-})e^{i\bm{\bar{K}}^{\prime}\cdot\bm{R}_{-}}\varphi_{A^{\prime}}(\bm{R}_{-})=0$. Moreover, the phase factors [see Fig. 1(a)] and orbitals cancel out, i.e., $e^{-i\bm{\bar{K}}\cdot\bm{R}_{-}}\varphi_{A}(\bm{R}_{-})=e^{i\bm{\bar{K}}^{\prime}\cdot\bm{R}_{-}}\varphi_{A^{\prime}}(\bm{R}_{-})$. Similar considerations follow for the $B$ sublattice on the right edge, $\bm{r}=\bm{R}_{+}$. Due to translational invariance along $y$, we can simplify $f_{\mu}(\bm{r})\rightarrow e^{ik_{y}y}f_{\mu}(x)$. Therefore, the boundary conditions for the zigzag nanoribbon envelope functions are

[TABLE]

In the original discussion by BF [40], they consider only a single valley ($K$ or $K^{\prime}$) on the $\psi(\bm{r})$ expansion. Indeed, near the $K$ valley $f_{\mu^{\prime}}(\bm{r})\approx 0$ (with $\mu=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^{\prime}$ valley, one gets $f_{A^{\prime}}(-W/2)=f_{B^{\prime}}(+W/2)=0$. Therefore, Eq. (12) is a generalization of the original BF boundary condition for the zigzag nanoribbon.

III.1.2 Armchair nanoribbons

For the armchair nanoribbons, confinement along $y$ projects both $K$ and $K^{\prime}$ valleys into $\bar{\Gamma}$ on the $k_{x}$ armchair axis [Fig. 1(b)]. Therefore, the $\psi(\bm{r})$ expansion in Eq. (9) is defined by $\bm{q}=\bm{q}^{\prime}\rightarrow\bm{k}$. As in the zigzag case above, these projections also apply to $H_{G}(\bm{k})$ in Eq. (10).

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 $\varphi_{\mu}(x,\pm L/2)\neq 0$. In this case the system has translational invariance along $y$, thus $f_{\mu}(\bm{r})\rightarrow e^{ik_{x}x}f_{\mu}(y)$. Imposing $\psi(x,\pm L/2)=0$ for all $x$ in Eq. (9), we get

[TABLE]

where $e^{\pm i\theta}=\varphi_{\mu^{\prime}}(x,\pm L/2)/\varphi_{\mu}(x,\pm L/2)$ is the phase difference between the $K$ and $K^{\prime}$ solutions at the edges. On the BF approach [40], they consider the bulk Bloch phase difference, yielding $\theta\rightarrow\theta_{BF}=\Delta K\cdot L=(N_{A}+1)2\pi/3$. However, since confinement breaks the Bloch periodicity along $y$, deviations from $\theta_{BF}$ should be expected for narrow ribbons. Therefore, hereafter we consider Eq. (13) with an arbitrary $\theta$ as a generalization of the original BF boundary condition.

III.2 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

[TABLE]

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

[TABLE]

Armchair nanoribbons. Substituting $U_{k}^{A}$ and $U_{c}^{A}(\theta)$ 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 $\eta=0$, replacing $U_{k}^{Z}$ and $U_{c}^{Z}(0)$ from Eq. (14) into Eq. (5), we reproduce our boundary conditions shown in Eq. (12), plus an additional pair of equations $f_{B}(\bm{R}_{-})=f_{B^{\prime}}(\bm{R}_{-})$, and $f_{A}(\bm{R}_{+})=f_{A^{\prime}}(\bm{R}_{+})$, which are trivially satisfied since, $\varphi_{A^{(\prime)}}(\bm{R}_{+})\approx 0$ and $\varphi_{B^{(\prime)}}(\bm{R}_{-})\approx 0$. Second, for $\eta=1$ the same procedure gives us the original BF boundary condition for the zigzag confinement. Therefore, the two possible values of $\eta=\{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.

III.3 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_{2h}$ group. Additionally, we consider that the system is time-reversal symmetric and chiral. The matrix representations of these symmetry operations are built from the $\psi(\bm{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.

III.3.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

[TABLE]

for zigzag and armchair nanoribbons, respectively. The Dirac-like term $h(\bm{k})$ is given in Eq. (10), and $\bm{k}^{\pm}=(k_{x},k_{y}\pm\bar{K})$ centers the Dirac cones into the $K$ and $K^{\prime}$ projections [Fig. 1(b)]. The $k^{2}$ terms are defined by the Wilson’s matrices $U_{w}^{Z(A)}\equiv U_{c}^{Z(A)}$ and masses $m_{Z(A)}$. Since the $k^{2}$ terms allow for a trivial boundary condition, i.e., $\psi=0$ at the edges, it can be numerically implemented via simple finite-differences schemes [41]. Hereafter, all results were obtained discretizing the coordinates with $\sim 100$ points, and the Wilson’s masses are chosen on the mid-range of Eq. (8).

For zigzag nanoribbons, Fig. 2 compares the generalized ($\eta=0$) and original ($\eta=1$) BF conditions. Indeed, the ideal results are given by $\eta=0$ in Fig. 2(a). Here, the absence of fermion doublers is due to their hybridization near $\bar{X}$ [Fig. 2(b)], which drives the doublers towards high energies as $m_{Z}$ increases. On the other hand, for $\eta=1$ the Hamiltonian $H_{Z}$ splits into uncoupled $\bar{K}$ and $\bar{K}^{\prime}$ blocks, each showing an independent doubler [Fig. 2(c)]. Since there is no hybridization, the doublers will always occur at $k_{y}^{\text{doubler}}\approx\pm[\bar{K}-2m_{Z}/(\hbar v_{F}\delta_{x}^{2})]$. Not even staggered-lattice implementations are able to fully eliminate the doublers in this case [31, 32]. Notice that in the analytical $k$-linear limit $\delta_{x}\rightarrow 0$ and the doublers are eliminated as $k_{y}^{\text{doubler}}\rightarrow\pm\infty$, thus justifying the use of the original BF boundary conditions in analytical calculations. For $k_{y}\approx\bar{K}$, the envelope functions $\bm{F}(\bm{r})$ for $\eta=\{0,1\}$ match well, as shown in Fig. 2(d)-(e). For $\eta=0$, a small $\bar{K}-\bar{K}^{\prime}$ hybridization occurs, which is absent for $\eta=1$ due to the block form of $H_{Z}$. In contrast, the $\bm{F}(\bm{r})$ of the doublers show nonphysical oscillations with a period given by the numerical step size $\delta_{x}$ [dotted lines in Fig. 2(e)].

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

[TABLE]

Indeed, the DFT data for the gaps obey these expressions, i.e., $E_{N_{A}}\propto 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 $\Delta\theta$, which we use to establish the correction $\Delta\theta$ to the BF boundary condition shown in Fig. 3(b). Considering $v_{F}=0.8\times 10^{3}$ nm/ps (see fits in the next section), we solve the $k$-linear and $k^{2}$ models with the corresponding $\Delta\theta\approx 20^{\circ}$, 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, $\Delta\theta$ 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 $\eta=0$ model couples the $\bar{K}$ and $\bar{K}^{\prime}$ 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}\rightarrow\pm\infty$. 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 $\Delta\theta\rightarrow 0$, yielding the zero gap for the $N_{A}=3p+2$ metallic case, and degenerate $3p$ and $3p+1$ gaps.

III.3.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

[TABLE]

The matrices $U_{ij}=\tau_{i}\otimes\sigma_{j}$ are set in terms of the Pauli matrices $\tau_{\nu}$ (with $\nu=0,x,y,z$) acting on the $K$/$K^{\prime}$ valley subspace, and $\sigma_{\nu}$ acting on the sublattices $A$/$B$. The first line in Eq. (21) contains the Dirac-like terms. For armchair $\bm{k}^{\pm}=\bm{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_{A1}=m_{A}\cot\theta$ and $m_{A2}=m_{A}\csc\theta$. Similarly, for zigzag $m_{Z1}=m_{Z}\eta$ and $m_{Z2}=m_{Z}(1-\eta)$. The third line shows the extra terms that allow for a fine tuning of the band structure. The $\Delta$ term couples the projected cones from $K$ and $K^{\prime}$ valleys, the velocity $\mu$ couples the dispersions at finite $\bm{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 $3p$), as shown in Fig. 4. The parameters used to obtain the figures, and equivalent results for $N_{A}=49$ ($3p+1$), $N_{A}=50$ ($3p+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|\gtrsim 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|\sim 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), $|\psi(y)|^{2}\propto|f_{\mu}(y)e^{i\theta y/L}+f_{\mu^{\prime}}(y)e^{-i\theta y/L}|^{2}$, where the phase factors arise from the $\varphi_{\mu^{(\prime)}}(\bm{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.

IV Spinful case: topological crystalline insulators

Graphene has a very weak spin-orbit coupling [1, 2, 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_{2h}$, while for ribbon C the symmetry is reduced to $C_{2v}$. 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 $\{\varphi_{xz,\uparrow}(\bm{r}),\varphi_{xz,\downarrow}(\bm{r}),\varphi_{x,\uparrow}(\bm{r}),\varphi_{x,\downarrow}(\bm{r})\}$, where the $xz$ and $x$ indices refer to the symmetries of the orbitals, and $\uparrow,\downarrow$ label the spin states along $z$. Thus, similarly to Eq. (9), the wave-function $\psi(\bm{r})$ expansion is given by these basis functions multiplied by the $k$-phase $e^{i(\bm{k}-\bm{X})\cdot\bm{r}}$ and the envelope spinor $\bm{F}(x,y)=e^{ik_{x}x}[f_{xz,\uparrow}(y),f_{xz,\downarrow}(y),f_{x,\uparrow}(y),f_{x,\downarrow}(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

[TABLE]

[TABLE]

The matrices $U_{ij}=\tau_{i}\otimes\sigma_{j}$ are set by Pauli matrices acting on the orbital ($\tau_{\nu}$) and spin ($\sigma_{\nu}$) subspaces, $\Delta$ is the gap at $\bm{k}=0$, $\alpha$ defines the Fermi velocity, and $m$ is Wilson’s mass. The coupling $\Delta_{C}$ is only allowed for ribbon C. Since the confinement is along $y$ the kinetic matrix $U_{k}=U_{xx}$, yielding $\det[U_{k}\pm iU_{c}(\rho,\theta)]=0$, as expected. The boundary conditions are defined by the continuum parameters $\rho$ and $\theta$. For ribbons A and B, $\theta\equiv 0$.

The effects of the boundary condition parameter $\rho$ and $\theta$ on the band structure are shown in Figs. 5(d)-5(f). Here, we consider $\Delta$ and $\alpha/\Delta$ as the energy and distance units. For $m>0$ the system is trivial, thus there are no states within the gap $|E/\Delta|<1$ in Fig. 5(d). For $m<0$ the system becomes topologically non-trivial with a mirror Chern number $n_{M}=-2$ [58], yielding two Dirac cones (degenerate for $\theta=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 $\rho=0$, while for finite $\rho$ the broken chirality is a consequence of the distinct atomic terminations of ribbons A and B [58]. For ribbon C we can consider $\rho=0$ for simplicity and allow $\theta$ to vary. This is shown in Fig. 5(f). In this case, both $\theta\neq 0$ or $\Delta_{C}\neq 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.

Complementarily, within the $k$-linear model ($m=0$), the BM approach for the boundary conditions $[U_{k}\pm iU_{c}(\rho,\theta)]\cdot\bm{F}(x,\pm L/2)=0$ yields,

[TABLE]

Interestingly, this boundary condition implies a spin-orbital admixture, as it couples opposite spins $\pm\sigma$ 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 non-trivial 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^{\prime}$ coupling induced by the boundary condition restricts the edge-state bands to lie within these valleys. More interestingly, for the armchair case, it introduces $\Delta\theta$ as a deviation from the bulk Bloch phases. Particularly, for the “metallic” armchair case a finite $\Delta\theta$ 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.2 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\times N$ (for even $N$) Hamiltonians with Dirac-like spectrum.

A.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 $(\hbar v_{F}U_{k}k+V_{0}U_{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_{k}^{-1}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_{k}^{-1}U_{c}$ is skew-Hermitian with purely imaginary eigenvalues.

Similarly, in the quadratic case, the Wilson’s mass term cannot introduce any new modes, i.e., $(\hbar v_{F}U_{k}k+mU_{w}k^{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\to\sin(\delta_{x}k)/\delta_{x}$] has no fermion doubling at $k=\pi/\delta_{x}$, requiring non-singular $U_{w}$. Dividing the Hamiltonian by $k$ we get $(\hbar v_{F}U_{k}+mU_{w}k)$ whose determinant should not vanish for any real $k$. Using that $U_{k}$ and $U_{w}$ are non-singular, this is equivalent to $U_{k}^{-1}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.

A.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=\hbar v_{F}U_{k}k+U_{c}V(x)$, and its corresponding Schrödinger equation $HF(x)=EF(x)$. Let us assume that $V(x)=V_{0}\Theta(x)$ defines a single wall at $x=0$, opening a gap $|2V_{0}|$ for $x>0$. On the outer region, the hard-wall limit $(V_{0}\rightarrow\infty$) allow us to neglect $E\ll V_{0}$, and considering the ansatz $F(x>0)\sim F_{0}e^{-x/\lambda}$ ($\real\lambda>0$), we get

[TABLE]

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

[TABLE]

since $U_{k}$ is non-singular. Here, $\omega=-i\hbar v_{F}/V_{0}\lambda$ plays the role of the eigenvalues of $U_{k}^{-1}U_{c}$. The restriction $\real\lambda>0$, imposed by the evanescent solution for $x>0$, implies $\imaginary\omega<0$. Next, one must check that $U_{k}^{-1}U_{c}$ admits eigenvalues with $\imaginary\omega<0$.

First, notice that Eq. (25) is formally identical to the constraint on $U_{c}$ derived in the previous section, requiring that there are no solutions with purely real $\omega$. Second, using that $U_{k}^{-1}U_{c}$ is a product of two Hermitian matrices, it is easy to prove that its eigenvalues come in complex-conjugate pairs, i.e., simply take the complex conjugate of Eq. (26) to find

[TABLE]

Therefore, both $\omega$ and $\omega^{*}$ are eigenvalues of $U_{k}^{-1}U_{c}$, and $*$ labels complex conjugation. This guarantees that there are always exactly $N/2$ solutions with $\imaginary\omega<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

[TABLE]

and the expectation value of the current is

[TABLE]

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_{0}^{{\dagger}}U_{k}F_{0}=F_{0}^{{\dagger}}U_{c}F_{0}=0$. This shows that the boundary condition is current conserving.

A.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<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<0)\sim F_{k}e^{ikx}$, with $\imaginary k\leq 0$ to avoid divergent states for $x\rightarrow-\infty$. The Schrödinger equation for $x<0$ becomes

[TABLE]

Nontrivial solutions require $\det[S(k,E)]=0$, with $S(k,E)=\hbar v_{F}U_{k}k+mU_{w}k^{2}-E$, yielding a $2N^{\rm th}$-order polynomial equation for $k$, where $N$ is the size of the matrix $S(k,E)$. For a fixed $E$, this provides $2N$ complex roots $k_{n}$, which are either purely real, or come in complex-conjugate pairs, as shown by a similar argument to Eq. (A.2). 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 $\imaginary k_{n}\leq 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 $\imaginary k_{n}<0$ is linearly dependent, thus validating the criteria stated above.

A.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\approx 0$ and momentum $k\approx 0$, for which the linear spectrum dominates.

First, consider the linear model $H$ from Eq. (6), and a plane-wave ansatz $F(x)=F_{n}e^{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 small-$k$ limit we can neglect the $k^{2}$ term, hence, the small-$k$ 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\rightarrow 0$) on the Schrödinger equation for the $k^{2}$ model, and simplifying $k$ in the equation to eliminate the small-$k$ solutions, we get

[TABLE]

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}=\hbar v_{F}/(m\omega)$ with $\imaginary k_{e}>0$, and $\omega$ is a complex eigenvalue of $U_{k}^{-1}U_{w}$ with $\imaginary\omega>0$. This equation is formally identical to the boundary condition of the $k$-linear model (4). Since the small-$k$ 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 ($P6/mmm$) symmorphic space group. Eliminating the trivial Bloch translations group $T_{B}$, its factor point group is ($P6/mmm$)/$T_{B}\sim D_{6h}$. 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 $\bar{\Gamma}$ for the armchair confinement, both have the same factor group under trivial Bloch translations, thus, ${Pmmm}/T_{B}\equiv{Pmma}/T_{B}\sim D_{2h}$. Similarly, zigzag nanoribbons with even (odd) $N_{D}$ belong to the $Pmmm$ ($Pmma$) space group. In this case, the $K$ and $K^{\prime}$ monolayer valleys fall into the $\bar{K}$ and $\bar{K}^{\prime}$ 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 $\text{Pmmm}/T_{B}\equiv\text{Pmma}/T_{B}\sim D_{2h}$. Therefore, hereafter it is sufficient to analyze the symmorphic lattices and the point group $D_{2h}$. This allows us to build a single model for both armchair and zigzag nanoribbons, considering a basis set that contains both $K$ and $K^{\prime}$ basis functions.

The $D_{2h}$ point group can be generated by its mirror operations $M_{x}$, $M_{y}$, and $M_{z}$. Here, $M_{x}$ reflects $x\rightarrow-x$, and similarly for $M_{y}$ and $M_{z}$. Under the basis vector $\bm{r}=(x,y,z)$, the coordinates representation from the $O(3)$ group is given by the matrices

[TABLE]

where $\text{diag}(\cdots)$ labels a diagonal matrix with elements given by its arguments.

To obtain the Hilbert space representation $\mathcal{H}$ we consider the basis functions $\{\varphi_{A}(\bm{r}),\varphi_{B}(\bm{r}),\varphi_{A^{\prime}}(\bm{r}),\varphi_{B^{\prime}}(\bm{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^{\prime}$ valleys. Namely, the representation matrices for the $D_{2h}$ generators are

[TABLE]

where $\bm{\sigma}=(\sigma_{0},\sigma_{x},\sigma_{y},\sigma_{z})$ are Pauli matrices acting on the $A$/$B$ lattices subspace, and $\bm{\tau}=(\tau_{0},\tau_{x},\tau_{y},\tau_{z})$ acts on the $K/K^{\prime}$ valley subspace.

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

[TABLE]

To obtain the effective Hamiltonian $H$, we consider the method of invariants [49]. Thus, we seek the most general form of $H\equiv H(\bm{k})$ as an expansion in powers of $\bm{k}=(k_{x},k_{y})$ that is invariant, i.e., $[H,\mathcal{S}]=0$ for all $\mathcal{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_{\rm ft}$, we obtain

[TABLE]

Here, $H_{0}$ represents a minimal Dirac-like model with anisotropic Fermi velocities $v_{x}$ and $v_{y}$, with the monolayer projected cones at $k_{y}=\pm\Delta_{K}=\pm\bar{K}$ for the zigzag confinement, and $\Delta_{K}=0$ for armchair. The $H_{A}$ and $H_{Z}$ show the most general $k_{x}^{2}$ and $k_{y}^{2}$ 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}\equiv U_{w}$ confinement or Wilson’s matrices for the minimal armchair and zigzag models, respectively. Additionally, $H_{\rm 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 1). 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 requires complex handling [31, 32].