Volkov-Pankratov states in topological graphene nanoribbons

In topological systems, a modulation in the gap onset near interfaces can lead to the appearance of massive edge states, as were first described by Volkov and Pankratov. In this work, we study graphene nanoribbons in the presence of intrinsic spin-orbit coupling smoothly modulated near the system edges. We show that this space modulation leads to the appearance of Volkov-Pankratov states, in addition to the topologically protected ones. We obtain this result by means of two complementary methods, one based on the effective low-energy Dirac equation description and the other on a fully numerical tight-binding approach, finding excellent agreement between the two. We then show how transport measurements might reveal the presence of Volkov-Pankratov states, and discuss possible graphene-like structures in which such states might be observed.


I. INTRODUCTION
Graphene was the first material theoretically predicted to be a Quantum Spin Hall (QSH) insulator. In the proposal by Kane and Mele [1,2], the intrinsic spinorbit coupling (SOC) opens a topological gap in the energy dispersion of the bulk system, and edge states appear in nanoribbons due to the bulk-boundary correspondence. While possible signatures of a topological gap in graphene have been recently reported [3], the minute size of the spin-orbit gap in pristine graphene (≈ 25 µeV) complicates the observation and application of the promising electronic and spin properties of the QSH edge states. Two different approaches have mainly been followed to overcome this limitation: a) Find ways to induce a stronger SOC in graphene, for example by depositing heavy adatoms on the graphene surface [4][5][6][7][8] or by proximity to materials with much stronger SOC than carbon such as transition metal dichalcogenides (TMDs) [9][10][11][12][13][14][15]. b) To grow graphene-like honeycomb structures made of heavier elements in groups IV and V [16][17][18][19][20][21][22][23][24]. While the experimental realization of QSH physics in these systems seems challenging, and even though there is so far limited evidence for the existence of protected edge states [14], the advances in the artificially induced SOC in graphene are promising [13,25].
The experimental approaches described above are likely to result into an inhomogeneous distribution of the strength of the induced SOC, with larger inhomogeneity near the sample edges. We therefore investigate, both analytically and numerically, the effects of a reduction of the SOC near the edges in wide graphene nanoribbons, both of zigzag and armchair type. Intrinsic SOC leads to band inversion at the K and K points, and due to the smooth SOC modulation at the edges, new massive edge states appear in addition to the topologically * These authors contributed equally; tineke.vandenberg@dipc.org † These authors contributed equally; ademarti@city.ac.uk ‡ dario.bercioux@dipc.org protected ones. Such massive edge states were first described by Volkov and Pankratov [26][27][28], and are therefore refered to as Volkov-Pankratov (VP) states. Recently, VP states have attracted attention again in the context of topological insulators (TIs) [29][30][31][32][33] and topological superconductors [34]. Although VP states in TIs are not topologically protected, they are of topological origin, because they result from the band inversion between a topological and a trivial material [30]. Threedimensional TIs with band inversion at the Γ-point were investigated both theoretically, within an effective linear in momentum model [30,32], as well as experimentally, in HgTe/CdTe heterojunctions [29,31]. Other studies focused on two-dimensional quantum wells within the Bernevig-Hughes-Zhang model [33]. In these works the VP states appeared due to the smooth modulations in the band structure near the edges. Recently, the coexistence of topological and trivial modes was also reported in 2D TMD, specifically in the 1T' phase of WSe 2 [24]. Within our numerical approach, in addition to the spectral properties, we investigate the transport properties of clean and disordered nanoribbons. In general, the conductance of the system increases by 4e 2 /h every time a new VP band opens to conduction. In disordered ribbons, the opening of a new conduction channel via a VP state is accompanied by the appearance of a dip in conductance. These dips resemble the ones observed in quasi-one dimensional quantum wires in the presence of an attractive impurity [35]. In this context, the decrease in the conductance is due to the coupling of propagating modes to quasi-bound states in the scattering region. Therefore, the presence of these dips in the conductance indicates that a new subband is opened, thereby demonstrating the existence of such bands within the topological gap. This hints that the presence of disorder is not detrimental to the detection of VP states in transport experiments.
This article is organized in the following way. In Sec. II, we discuss the analytical solution for the edge mode spectrum of a semi-infinite graphene flake with a space modulation of the intrinsic SOC close to the boundary, within In the analytical calculation the system is a semi-infinte half plane. In the numerical calculations it has a ribbon geometry, attached to source and drain leads. Depending on the orientation, it is a system with either zigzag, or armchair edges. (b) Sketch of the profile of the space modulation of the intrinsic SOC given by Eq. (3). The system occupies the unshaded region.
the long-wavelength approximation. In Sec. III, we solve for the full spectrum of zigzag and armchair nanoribbons with the modulated intrinsic SOC using the tight-binding model. In Sec. IV, we analyse the transport properties within the full tight-binding description, including the effects of disorder. Finally, in Sec. V, we present conclusions and outlook. Some technical details are relegated to the three appendices at the end of the paper.

II. ANALYTICAL LOW-ENERGY APPROACH
In this section we investigate the spectrum of edge states of graphene nanoribbons by using the Dirac equation description [36,37]. As it is well-known, this emerges in the long-wavelength approximation (LWA) of the tight-binding Hamiltonian around the Dirac points K and K . Within this continuum model, we account for the presence of non-uniform intrinsic SOC. We focus on the part of the spectrum corresponding to states localized at the edges. These states decay exponentially away from the edges in the bulk, so if the width of the nanoribbon is much larger than their decay length, the two edges can be treated separately. We assume that this is the case, and study a semi-infinite system with a single edge of either zigzag or armchair type. Our results below, therefore, apply to wide nanoribbons, and do not account for finite-size effects. Within the LWA approximation, graphene's effective Hamiltonian is given by where v F denotes graphene's Fermi velocity, p = −i ∇ is the momentum operator, and {τ x , τ y , τ z }/{σ x , σ y , σ z }/{s x , s y , s z } are the Pauli matrices associated to the valley/sublattice/spin degree of freedom, respectively. Since the Hamiltonian (1) is diagonal in valley and spin space, we will work in a basis of given valley and spin projection: and we shall omit the valley and spin indices wherever there is no risk of confusion. In Eq. (1), we include a nonuniform intrinsic SOC with a smooth spatial modulation transverse to the boundary: Here, ζ represents the coordinate in the direction perpendicular to the boundary, located at ζ = 0, and the system extends on the side of positive ζ. Equation (3) describes a domain-wall profile centered at ζ = ζ 0 , with characteristic modulation length , and with asymptotic values given by Then ∆ i represents the SOC deep in the interior of the nanoribbon. We assume that the modulation occurs close to the boundary, where the SOC reduces to∆, with ζ 0 of the order of graphene's lattice constant a 0 , see Fig. 1 The choice of the hyperbolic tangent profile is convenient because it allows for an exact solution of the corresponding Dirac equation [30,38]. However, we expect that the qualitative features of the spectrum do not depend on the detailed shape of the profile, as long as the typical length scale of the SOC modulation is large on the scale of the lattice constant a 0 . Since this is also the condition for the validity of the LWA employed here, throughout this paper we assume a 0 .

II.1. Spectral properties of zigzag ribbons
We start by considering a semi-infinite graphene flake extending in the region y > 0, with a zigzag edge along the x-axis [see Fig. 1(a)], and SOC profile given by Eq. (3) with ζ = y. By exploiting the translational invariance along the x-direction, the wave function can be expressed in the form with ψ T (y) = (v A , v B ). Then, the Dirac equation reduces to where we set = 1, and measure energies in units of v F / , lengths in units of , and wave vectors in units of −1 . Equation (6) admits an exact solution [30,38], whose derivation is summarized in App. A. Here, we just present the result. We introduce the notation Then, in terms of the new variable u given by the sublattice amplitudes can be expressed as with α = A, B = ±. The functions w ± (u) satisfy a hypergeometric equation (see App. A). Selecting the solution that leads to normalizable states for y → +∞ (i.e., u → 0), we obtain where F (a, b; c; z) is the ordinary hypergeometric function [39]. Normalizability requires κ i > 0, but imposes no constraint on κ e . Therefore, at fixed k x , the edge modes exist in energy window |E| < ∆ 2 i + k 2 x . Since for y → +∞ we have u ∼ e −2y , κ −1 i actually represents the decay length of the corresponding edge mode into the bulk. The relative factor between the two spinor components is fixed by the Dirac equation. We find We now need to impose the appropriate boundary condition. In the case of a zigzag edge, the boundary condition requires that one sublattice component of the wave function vanishes at the boundary y = 0, separately for each valley and spin [36,37]: where α = A (respectively α = B) if the boundary sites belong to the B (respectively A) sublattice. Using Eqs. (8), (9), and (10), from Eq. (11) we obtain where u 0 = 1 2 (1 + tanh y 0 ). The shift y 0 is important for the comparison to the numerical tight-binding analysis, in which the SOC profile is centered exactly at the edge of the system, i.e., on the first line of carbon atoms [40]. In the continuum approach, this corresponds to the choice y 0 = a 0 / √ 3, because y = 0 corresponds to a line of auxiliary sites, where one imposes the vanishing of the wave function.
Before discussing the solutions of Eq. (12), we notice that this equation is invariant under the following transformations: The first invariance is just the consequence of the timereversal symmetry of the Hamiltonian (1): the edge Band structure of a zigzag nanoribbon of width Ly = 37.2 nm (150 rows), with modulation of the SOC at the boundary. We set λi = 0.01t, λe = −0.05t, and = 12a0. The gray, black, and red lines are the tight-binding results, describing bulk states, VP states, and topological states, respectively. The analytical results from Eq. (12) are represented as blue dots. The right inset shows the case of homogeneous SOC, with λ = 0.1t, corresponding to a gap ∆ = 3 √ 3λ ≈ 0.52t. In both cases, two topological modes (red lines) cross within the gap. The modulation of the SOC results in additional massive edge states, popping up under the conduction band and above the valence band (black lines).
states occur in pairs of counterpropagating modes with opposite spins, residing on opposite valleys. The second invariance implies that the edge state spectrum at a B-type edge can be obtained from the spectrum at an A-type edge by simply reversing energy and wave vector. In order to reproduce the full spectrum of edge states (including degeneracies) of a wide zigzag nanoribbon, which has one edge of A sites and one edge of B sites, we need to take the solutions of Eq. (12) with α = +1 and with α = −1. These solutions are shown in Fig. 2 as blue dots on top of the numerical results discussed in the next sections (gray, black, and red lines). For clarity, we only include the states at one valley, the states at the other valley follow by symmetry.
We observe that the continuum approach faithfully reproduces all the main features of the spectrum close to the Dirac points. Within the gap, there exist two topological bands with approximately linear dispersion, and ten massive VP modes (for the given values of the parameters), in agreement with the results of the numerical approach discussed in Sec. III below. All these levels are doubly degenerate if one considers a system with two edges. Figure 2 shows the agreement between analytical and numerical results. For wave vectors in the interval between the two Dirac points, the agreement is less satisfactory, which we attribute to the fact that the coupling between the valleys, neglected in the continuum approach, plays an important role at these wave vectors. This is especially evident for the topological bands. Another source of discrepancy stems from neglecting higher order terms in momentum in the LWA.

II.2. Spectral properties of armchair ribbons
Let us now turn to the case of a semi-infinite system with an armchair edge along the y-direction [see Fig. 1(b)]. The wave function can be written as Then the Dirac equation reduces to and its solutions can be written as The functionsw ± are, again, solutions of a hypergeometric equation, and are given bỹ (16) with the relative factor, fixed by the Dirac equation, given by Notice that the dependence on the valley index only appears in the prefactor.
In the armchair case, the boundary condition involves both valleys [36,37] and reads with α = A, B. In terms of the w ± we find From Eq. (17) we see that eitherw −,τ =−1 = −w −,τ =1 or w +,τ =−1 = −w +,τ =1 , therefore the boundary condition takes the simple form Here, u 0 = 1 2 (1 + tanh x 0 ). In the armchair system, the shift required to have the correct value of the SOC on the first line of carbon atoms is x 0 = a 0 /2. edges. The obvious solutions E = ∓k y describe two counterpropagating linearly dispersing topological bands. These modes only exist if s = ∓1, respectively, otherwise one gets the trivial solution w + = w − = 0. The allowed value of s guarantees that the corresponding wave function is normalizable. This is a manifestation of the spin-momentum locking. Remarkably, in contrast to the zigzag case, here the group velocity of the topological modes is equal to v F and does not depend on SOC. Moreover, we observe that, for E = ±k y , the boundary condition depends on k y and E only through the combination k 2 y − E 2 . Therefore we find solutions corresponding to VP states, whose dispersion has the form where the effective masses M n and the number of massive states N max depend on the parameter values. In Fig. 3, we compare the analytical results with the numerical ones, finding an excellent agreement. Since in this case the continuum approach incorporates the coupling between valleys in the boundary condition, it is not surprising that the agreement is better than in the zigzag case.

III. NUMERICAL TIGHT-BINDING MODEL
In order to be able to go beyond the low-energy approximation, as well as to study the transport properties, we now move on to a fully numerical approach within a tight-binding formalism, which we implement using KWANT [41]. Contrary to the analytical calculation, we will here consider a finite size system, comprising of a scattering region with two edges, along which edge states can propagate, and a source and a drain lead. We will take rather large ribbons, but still of experimentally relevant sizes, with length L ≈ 60 nm and width W ≈ 35 nm. For this width, which largely exceeds the modulation length = 12 a 0 ≈ 3 nm, the electronic states located on opposite edges do not overlap. Hence, the edges are independent, and the numerical results can readily be compared to the analytical solution for a semi-infinite system.

III.1. The model and its parameters
Within a tight-binding formalism, the Hamiltonian for graphene with intrinsic SOC reads [1,2]: where c † ns (c ns ) creates (annihilates) an electron with spin s on the site n, and the symbol . . . ( . . . ) indicates sum over nearest (next nearest) neighbour sites. In Eq. (22), the sign ν nm = ±1 depends on the orientation of the next nearest neighbour hopping: it is positive (negative) for electron making a left (right) turn to the next nearest neighbour carbon atom. The hopping parameter is t, and λ is the intrinsic SOC parameter, which is related to the gap size as ∆ = 3 √ 3λ. We consider the following space modulation of the intrinsic SOC along the coordinate corresponding to the lateral width of the graphene nanoribbon: where L ζ is the width of the ribbon in the ζ-direction, and λ i(e) is the value of the SOC in the internal (external) region of the ribbon, respectively. Throughout this paper, we use λ i = 0.1t and λ e = −0.05t. The length scale characterizes the size of the spatial region over which the variation of the intrinsic SOC takes place. This has to be compared with the three natural length scales present in the system: the lattice constant a 0 , the length scale associated to the SOC gap ξ = v F /∆, and the width of the ribbon L ζ . In order to get VP states, one has to assure ξ. Moreover, to resolve the smoothness of the SOC modulation one needs a 0 , and the two edges are independent for L ζ . In App. B we provide some details on the relation between the parameter values in the LWA and in the tight-binding model.

III.2. Spectral properties
The spectral properties for graphene nanoribbons with homogeneous SOC are shown in the insets of Figs. 2 and 3.
Without SOC modulation, each edge of the system hosts two topological states with linear dispersion, as shown in the insets in red. The bulk states, in gray, form the conductance band (CB) and the valence band (VB).

III.2.1. Zigzag ribbons
In the zigzag case, we consider a ribbon of width W = L y = 37.2 nm, corresponding to 150 rows. Figure 2 shows that, due to the suppression of the SOC near the edges, VP bands are pulled out of the CB and the VB, in symmetric fashion. Near the K and K points, these VP modes push away the topological modes, whose group velocity is thus strongly affected. This feature can be rationalized by observing that, from Eq. (12), one can see that the group velocity of the topological modes close to the Dirac points depends strongly on the gap parameters ∆ i/e . Near the boundary, as the SOC gets weaker, the effective gap becomes smaller. Looking at the VP modes under the CB, we observe that the first (lowest in energy) VP mode is the one which lies closest to the edge, and has the smallest decay length. Each consecutive VP mode has a longer decay length and extends deeper in the bulk, and therefore "sees" a slightly larger effective gap. In Fig. 4 the local density of states (LDOS) of the edge states is plotted at different energies. In Fig. 4(a) there is only the topological mode, in 4(b) there is additionally one VP mode, in 4(c) two VP modes etc. In the zigzag case we observe that there are zones on the lattice with predominant A or B contributions to the LDOS. Moreover, for each VP mode that is added, the predominance changes sublattice.

III.2.2. Armchair ribbons
In the armchair case, we consider a ribbon of width W = L x = 32.5 nm, corresponding to 128 unit cells. The band structure is shown in Fig. 3. The two Dirac points are projected onto the same point k = 0. In the absence of SOC, an armchair ribbon is metallic or semiconducting, depending on its exact width [37]. In the case studied here, the ribbon is wide enough that the semiconducting gap, which is of order of v F /L x , is exceedingly small compared to all other energy scales, and can be ignored. In the presence of SOC, the gap size is then determined by the SOC strength, and there are always two topological modes crossing the gap. With the SOC modulation, new massive bands appear under the CB and above the VB. However, in contrast to the zigzag case, the dispersion of the topological bands is not affected by the appearance of these new levels, and their group velocity is not modified. This is consistent with the analytical solution of Eq. (20), which gives a linear dispersion with slope v F , independent of the values of the gap parameters. As shown in Fig. 3, the agreement between the analytical and numerical results for all edge states is excellent.
In the LDOS for an armchair ribbon we see there are no privileged A or B sublattice zones - Fig. 5. It is harder in this case to say where on the lattice the different VP states lie. Also in this case, we observe the expansion of the area containing the edge states, as one consecutively adds subbands.

IV. TRANSPORT PROPERTIES
In this section we will investigate if two-terminal conductance measurements on nanoribbons with SOC modulation near the edges could reveal the presence of VP states. Being massive, VP states are not protected against backscattering due to disorder. We will therefore add Anderson disorder to the model, namely random spin-independent onsite energies with uniform probability distribution in the interval ε n ∈ [−U 0 /2, U 0 /2]. The disorder Hamiltonian can be written as where the sum runs over the entire system. The topological edge states will not be sensitive to this disorder, but the VP states will. How much the disorder weakens the VP states depends not only on the disorder strength, but also on the system size. Here, we consider a sample of size L × W ≈ 60 × 35 nm 2 .

IV.1. In-gap conductance
The numerical results for the two-terminal conductance are shown in Fig. 6. As one can expect for a system presenting in-gap states, the conductance never decreases to zero. Here, we have two counterpropagating doubly-degenerate topological edge states in the gap, so the minimal in-gap conductance is 2e 2 /h. In the absence of SOC modulation, this is the only in-gap contribution to the conductance. With SOC modulation near the edges, we observe clear steps at the energies at which new VP modes open, see upper panels of Fig. 6. Both in the zigzag as well as in the armchair cases, these steps are symmetric around the gap center E = 0.
In the presence of disorder, the conductance due to the VP modes is suppressed, because those modes are sensitive to backscattering. For strong disorder, the conductance reduces to 2e 2 /h, below which it can not descend because the topological modes are not sensitive to Anderson disorder, which does not break time-reversal symmetry. However, in this case the conductance in the CB and the VB is also significantly reduced due to the disorder in the ribbon. Remarkably, within a certain range of disorder strength, we observe dips in the conductance at the step edges. This reminds of the physics of a quasi-1D quantum wire containing an impurity with an attractive potential, where the conducting modes couple to a quasibound state at the impurity [35].

IV.2. Dip behavior near the steps
In the lower panels of Fig. 6, we observe a minimum at each step in the conductance curve of the disordered system. Notice that each line represents the conductance averaged over 100 disorder configurations. Plotting single disorder configurations [c.f. App. C], one observes random fluctuations at the plateaus, which average out over many disorder configurations, as well as a dip at the step edges, which most configurations have in common, and which therefore remains in the averaged conductance. Such dips in the conductance at the opening of new subbands were discussed in detail in a paper by Bagwell [35,42], in which he shows how propagating modes in a narrow wire with parabolic confinement couple to the zero energy quasi-bound states of delta shaped, negative potential impurities. In the case presented here, in the energy gap we have to deal with a quasi 1D subsystem near the edge, with triangular confinement potential. The quasi-bound states can be hosted in local energy minima, that appear due to the random energy landscape. In order to verify that we are indeed dealing with this physical phenomenon, we have simulated a clean sample, containing one Gaussian-shaped impurity at each edge. We investigated the minimal requirements for observing dips in the conductance curves right at the onset of the steps. The simulations with single impurities clearly demonstrate that the dips come from the coupling of the propagating states to quasi-bound states lying at randomly distributed energy minima on the lattice. This can be observed in strongly confined quasi-1D systems with an attractive impurity (negative potential). The precise properties of the dip in the conductance depend on the shape and strength of the impurity, and other system details. Additional information is given in App. C.
Although the situation near the edges in the system presented here differs from that investigated in earlier works, many of the physical arguments still hold [35,43]. Because of the random energy landscape, attractive sites or regions on the lattice may result in the appearance of quasi-bound states lying just under each VP subband. As we can observe, and as is normally the case, the evanescent state under the lowest VP subband is a true bound state, as it does not couple to the topological state. We therefore observe no dip before the first step, i.e. conductance never decreases below 2e 2 /h. We also observe the usual renormalization of the energy gap, which slightly shifts the energies of the onsets of the steps, and the value of the CB opening [44][45][46]. The binding energy, defined as the difference in energy between the step edge and the dip minimum, is therefore hard to quantify. However, we do observe that, as disorder gets stronger, the dips become deeper, which is in line with earlier observations. In wires with a parabolic confinement potential and single attractive impurities, the interaction strength between the impurity and the various available subbands depends on the lateral position of the impurity in the wire. For our triangular confinement near the edge there is no symmetry around any center, however, we know where at the edge each of the modes lies, and if it can interact with the impurity. In the implementation of our system with one single impurity at each edge, we clearly observe how the impurity interacts with consecutive edge states as we move it away from the edge towards the bulk of the sample, as discussed above. A more in depth investigation is beyond the scope of this work, but will be the subject of a future investigation.

V. CONCLUSION AND OUTLOOK
In this work, we have investigated the appearance of Volkov-Pankratov edge states in topological graphene nanoribbons of zigzag and armchair type. In the presence of intrinsic SOC which is smoothly suppressed near the edges, the well-known QSH edge states are accompanied by multiple massive VP states. We have demonstrated their existence by means of two complementary methods, the exact analytical solution of the low-energy effective Dirac equation, and the numerical tight-binding approach, finding good agreement between the two. Transport simulations show how the VP modes contribute to transport, also in the presence of disorder. We observe dips in the conductance at the onset of each VP band, which are due to the coupling of the propagating VP states to evanescent modes present in the random energy landscape. At sufficiently strong disorder, the VP states are entirely suppressed, and cease to contribute to transport.
Our results can be relevant to experimental systems in which the intrinsic SOC in graphene is enhanced by one of the methods mentioned in the Introduction. Both the deposition of adatoms, as well as the proximitization with a TMD layer, would likely give rise to a inhomogeneous intrinsic SOC, especially at the edges of the system.
In addition to a possible implementation in graphene or post-graphene materials, the presence of these VP states accompanying the topological modes could be achieved in systems of ultracold atoms in optical lattices. The advantage of a realization within this platform is related to the flexibility to control the parameters separately across a large range compared to condensedmatter systems, where the system parameters are generally fixed by the material properties and by the sample geometry [47].
The time-reversal symmetric Kane-Mele model for the QSH effect [1,2] can be thought of as a double copy of the Haldane model [48] in which time-reversal symmetry is broken. The implementation of the Kane-Mele model could be achieved within the state-of-the-art technology for ultracold atoms; the honeycomb lattice and the Haldane model have been already realized in this context [49,50]. In practice, the Kane-Mele model could be implemented by using an internal atomic state as a spin degree of freedom. For each spin, the same scheme as for the Haldane model could be used to implement the second-next-neighbour hopping, see Ref. [51][52][53]. This system would then correspond to two copies of the Hal-dane model [50]. Contrary to the implementation in Ref. [54], we propose to realize a system with a homogeneous intrinsic SOC and soft-boundary conditions, corresponding to an inhomogeneous onsite energy profile due to the confining potential of the atomic trap. Similar to the results presented in this work, this scheme gives rise to a set of VP states accompanying the topologically protected one, but with energy symmetry breaking. A similar approach was proposed for the case of the quantum Hall edge states [55]. Other aspects of the implementation of this model for ultracold atoms require further investigation. In order to make this paper self-consistent, in this appendix we provide the details of the exact solution of the Dirac equation in the presence of an inhomogeneous SOC with hyperbolic tangent profile. The analysis follows Refs. [30,38].

Zigzag case
We consider first the case of a semi-infinite system with a zigzag edge along the x-direction, see Fig. 1(a). After factorization of a plane wave in the x-direction with wave vector k x , the Dirac equation reads where we set = 1, and measure energies in units of v F / , lengths in units of , and wave vectors in units of −1 . We look for solution on the half-line y ≥ 0. Squaring Eq. (A1), we obtain Next, we express ψ as where |± are the eigenvectors of σ x , with respective eigenvalue ±1. By performing the change of variable Eq. (A3) can be rewritten as where we use the notation Substituting in Eq. (A3) the ansatz we find the hypergeometric equation We need to select the solution which leads to normalizable states for y → +∞, i.e., u → 0. We find where F [a, b; c; z] is the ordinary hypergeometric function [39]. Notice that F [a, b; c; z] = F [b, a; c; z]. The other solution to the hypergeometric equation does not lead to normalizable states and we omit it. Since we consider a semi-infinite system, in contrast to Ref. [30], we do not require normalizability for y → −∞, but we need to impose the appropriate boundary condition at y = 0, as discussed in the main text. The relative factor between the two components is fixed by the Dirac equation, and we find Then, to summarize, up to an overall normalization factor, we have We observe that, since for y → +∞ we have u ∼ e −2(y−y0) → 0, and F [a, b; c; u] → 1, the decay of the wave functions in Eq. (A6) is controlled by the parameter κ i , which can then be identified with the inverse decay length of the corresponding edge state.

Armchair case
In the case of a semi-infinite system with an armchair edge along the y-axis, after factorization of a plane wave in the y-direction with wave vector k y , the Dirac equation reads Squaring Eq. (A7), we obtain In this case, we express ψ as where |± denote now the eigenvectors of σ y , with respective eigenvalue ±1. Following the same steps as in the previous subsection, we find with the prefactors given by We mention in passing that the solution for the arm-chair case can also be obtained by an appropriate π/2rotation of the solution for the zigzag case. Notice that while the overall phase of the wave function is immaterial, the relative phase between φ + and φ − is important. The equation (A10) implies that either φ + or φ − has opposite signs at the two valleys, while the other has the same sign. This observation turns out to be important when we impose the boundary condition, as discussed in Sec. II.2.
Appendix B: On the units conversion between tight-binding model and LWA results As already mentioned in the main text, the SOC parameters in the continuum Dirac equation and in the numerical tight-binding model are related as In the tight-binding model, we measure energy in units of the hopping amplitude t, length in units of the lattice constants a 0 , and the Fermi velocity is given by a0t . In the continuum, we measure energy in units of v F / , and wave vectors in units of −1 . Then, the conversion formulas are where we have inserted the value = 12a 0 used throughout this paper. The linear relation E = v F k in continuum units becomes E c = k c , and in tight-binding units becomes

Zigzag case
With the zigzag boundary along x, we have where the continuum coordinate y is given by and y 0 = a 0 / √ 3. Since x = na 0 (n ∈ Z), the onedimensional Brillouin zone in the transport direction is 0 < k x < 2π/a 0 .

Armchair case
With the armchair boundary along y, we have with the continuum coordinate x = na 0 , n ∈ Z , and x 0 = a 0 /2. In this case, the continuum coordinate in the transport direction is y = √ 3na 0 , thus the corresponding Brillouin zone is |k y | < π √ 3a0 . In the plots showing the tight-binding band structure, however, the wave vectors are rescaled in such a way that the onedimensional Brillouin zone appears to be |k tb | < π, see Fig. 3. Therefore, when comparing continuum and tightbinding results, it is important to take into account this additional √ 3 rescaling factor. In particular, the linear dispersion E = v F k appears in the numerical results as Appendix C: Conductance minima due to bound states Conductance curves for single disorder configurations in the armchair case for U0 = 0.1t. At the conductance plateaus, random fluctuations eventually cancel each other out, whereas the dip near the conductance step is in common to most disorder configurations, therefore it does not average out.
In Fig. 7 we show conductance curves of single disorder configurations. It can be seen that the dip at the conductance step is in common to most realizations. This is because most disordered energy landscapes have room for quasi-bound states within certain local energy wells. These states, lying at energies just under the opening of the next VP mode, couple to the propagating mode, which therefore localises at that energy, thereby decreasing the conductance.
To have a qualitatively better understanding of what causes the dips that we systematically observe in the conductance curves for a disordered system, we have performed additional numerical simulations. Suspecting these minima are due to the coupling of the VP modes to quasi-bound states in the system, we run simulations in the more conventional setting of a clean nanoribbon with a single attractive impurity. In this case, we put a single impurity on each edge, shifted away from the edge Conductance for a clean system containing one Gaussian shaped impurity on each edge. In (a) σ = 3a0 and in (b) σ = a0, with u0 = −0.4t and positions xL = xR = Lx/2, and yL = 8.848a0 and yR = Ly − yL. In (a), due to the width of the impurity, all propagating VP modes couple to the evanescent mode at the impurity. In (b), because the impurity potential is much narrower, certain propagating modes do not interact with the evanescent mode at the impurity.
by the same amount. This makes the impurity on each edge interact with the same VP mode. The impurity has a Gaussian shape U 0 (ζ) = u 0 e − (ζ−ζ L ) 2 +(ζ−ζ R ) 2 2σ , centered at ζ L = (x L , y L ) on one edge and ζ R = (x R , y R ) on the other edge. The impurity is therefore fully characterised by its position, strength u 0 < 0, and width σ (or variance σ 2 ). As discussed in the main text, we indeed observe dips in the conductance near the steps. How many dips we see, at which steps, and their shape, depend on the properties of the impurities. For impurities that are very narrow, such as σ = a 0 , we observe dips only at certain steps, and not at others, depending on where the impurity lies -c.f. Fig. 8(b). For wider impurities, such as σ = 3a 0 , we observe dips at all steps, as well as other dips in the plateaus -c.f. Fig. 8(a). Additional dips can result from having multiple evanescent modes at the impurity due to its finite width, or geometrical resonance effects. For the case of quasi-one dimensional quantum wires, it is also known that in the presence of Rashba SOC [56] the dips never go all the way down to the level of the previous conductance step, but are lifted proportional to the forth order in the Rashba SOC parameter [57]. A systematic study of all these features is left to future investigations.