Non-conservation of the valley density and its implications for the observation of the valley Hall effect

We show that the conservation of the valley density in multi-valley and time-reversal-invariant insulators is broken in an unexpected way by the electric field that drives the valley Hall effect. This implies that fully-gapped insulators can support a valley Hall current in the bulk and yet show no valley density accumulation on the edges. Thus, the valley Hall effect cannot be observed in such systems. If the system is not fully gapped then valley density accumulation at the edges is possible and can result in a net generation of valley density. The accumulation has no contribution from undergap states and can be expressed as a Fermi surface average, for which we derive an explicit formula. We demonstrate the theory by calculating the valley density accumulations in an archetypical valley-Hall insulator: a gapped graphene nanoribbon. Surprisingly, we discover that a net valley density polarization is dynamically generated for some types of edge terminations.

Introduction-The valley Hall effect (VHE) in nontopological systems has recently stirred considerable controversy [1][2][3][4][5][6][7][8][9].When the band structure features two valleys with a non-vanishing distribution of Berry curvature, electrons skew in the direction orthogonal to the applied electric field, even in the absence of magnetic field.However, since the system is not topological, electrons originating from one valley skew in the opposite direction of those from the other valley giving rise to a zero (charge) Hall current but to a finite valley Hall current j v (r, t).This is defined as the difference between charge currents of electrons originating in opposite valleys.When this current hits the edge of the system, a valley density n v (r, t) (or, more physically, a density of orbital magnetic moment [10]), is expected to accumulate at its boundaries.This assumes that the valley density obeys a standard continuity equation [5,6].This seems a reasonable assumption: the two valleys are wellseparated in momentum space, up to the point that they could ideally be taken as completely disconnected.
Some authors [1,6] went further and claimed that even a fully-gapped non-topological insulator such as graphene aligned with hexagonal boron nitride (hBN) [3,4] can exhibit nonlocal charge transport mediated by transverse undergap valley currents flowing in the bulk of the material.The authors of Ref. [6] argued that, at finite temperature, the valley-density accumulation could drive a "squeezed edge current" (parallel to the edges) in apparent agreement with experimental observations [2].However, other authors [7][8][9] found from microscopic calculations that there is no valley density accumulation and no edge current in the simple graphene/hBN model.They proposed that the observed nonlocal resistances are caused by substrate-induced edge states crossing the Fermi level [7] or by substrate-induced valley-dependent scattering [9].In the case of a fully gapped insulator this leaves us with the following puzzle: on one hand, the electric field drives a finite dissipationless valley Hall current in the bulk; on the other hand, time reversal symmetry implies that a valley density accumulationa time-reversal-odd quantity-cannot appear in response to an electric field, unless there is dissipation, which is impossible if there are no states at the Fermi level.So where did the valley current go?
In this paper we solve the puzzle by observing that the valley density does not satisfy the conventional continuity equation when an electric field is present.This includes the field applied in order to drive the valley Hall current.The reason is that the electric field breaks the conservation of crystal momentum and therefore of valley number, which depends explicitly on crystal momentum.As a result, the bulk valley current is internally shortcircuited as electrons flow from one valley to the other (and thus switch the sign of the Berry curvature) under the action of the very same electric field that drives the valley Hall current in the first instance.The process is schematically illustrated in Fig. 1.
Our results have profound implications for the observation of the VHE [11].In a fully gapped time-reversal invariant insulator the undergap valley current is incapable of producing a valley density accumulation on the edge.This makes observing the VHE impossible in such systems, unless, e.g., the valley degeneracy is lifted (see, for example, Refs.[12][13][14]) or carriers are selectively injected into a single valley via circularly polarized light [15].In these cases, however, an anomalous Hall effect "in disguise" is measured.This also means that the non-local resistance detection in Refs.[1,[16][17][18][19] must have been caused by partially occupied bands or edge states.
In metallic systems, which support a Fermi surface, our predictions for the valley density accumulation are quite different from those of the conventional theory which treats the entire bulk valley Hall current as the source of the accumulation.In particular, the value of the accumulation depends on the form of the electronic wave functions near the edge.The length over which it occurs is not related to the carrier diffusion length as in, e.g., Ref. [5], but reflects the much shorter localization length of edge states, as observed in some experiments [20], or the Fermi wavelength of bulk states.Perhaps the most important result of this study is that the valley density in the VHE is not simply transported from one edge to the other: it can be simultaneously generated on both edges by processes that involve the electric field in the bulk of the material.
Summary of main ideas-We consider a generic system in the shape of a strip of finite width which is indefinitely extended along the x axis.As we show below, the continuity equation satisfied by the valley density is where the electric field is in the x direction, which is parallel to the edge, and the valley current is in the y direction, perpendicular to the edge [21].The system is assumed to be macroscopically homogeneous along x so that the valley density and current depend only on y.The electronic states (in the absence of the electric field) are taken to be of the form ψ k,n (x, y) = e ikx u k,n (x, y)/ √ 2π where k is the x-component of the Bloch wave vector and n is the band index.The sum over k in Eq. ( 1) stands for dk/(2π).The mixed electronic distribution f k (y) = a −1 a 0 dx n f k,n u k,n (x, y) 2 is defined in terms of the electronic wave functions and the occupations of the corresponding states f k,n , with the integral taken over one period a in the x direction.S(k) is a "valley charge" function (odd under time reversal), which is a smooth periodic function of k in the Brillouin zone.We assume that the band structure features only two valleys, thus S(k) assigns number +1 to states around one valley and −1 to states around the other valley.The valley density operator is nv (r) ≡ −(e/2) j S( kj ), δ(r − rj ) where rj and kj are the position and Bloch momentum operator (along the edge) of the j-th electron, respectively, and { Â, B} ≡ Â B + B Â [22].The valley current density is ĵv (r) ≡ −(e/4) j S( kj ), vj , δ(r − rj ) , where vj is the velocity operator of the j-th electron.
Because k is a constant of motion, nv and ĵv (r) obey a conventional continuity equation in the absence of the electric field.
As we show below, in a fully gapped time-reversal invariant insulator, in which no edge or bulk state crosses the Fermi level, and at zero temperature, the right-hand side of Eq. ( 1) completely cancels the nearly-quantized contribution due to the second term on the left hand side.In this case, therefore, the valley density accumulation vanishes, even though there is a finite valley current in the bulk.In all other cases the cancellation is not exact.The correct equation for the density accumulation in the absence of relaxation processes is then , where the source term is a Fermi surface property.Note that Q s (y) cannot be written, in general, as the divergence of a current.In fact, this is only possible if its integral over the whole strip vanishes, which implies that density accumulates at one edge and depletes at the other [23].However, if the width of the strip is macroscopically large, the source term is localized on the edges.One can then define the "effective current" I s , obtained by integrating Eq. ( 2) across a given edge, that feeds the valley number accumulation thereat.Since valley density is not conserved, the sum of the effective currents associated with the two edges does not have to be zero.I s can be split as is the contribution of edge states.Here, the sum over e is that over edge states.The calculation of the contribution of bulk states, I b s , is complicated by the fact that the integral over y cannot be extended to infinity before performing the sum over n and k: the result would diverge.Nevertheless, a closed expression can be obtained in terms of the probability amplitude for Bloch waves to scatter off the edge [see Eq. ( 6) below].Once I s is known, the valley number accumulation can be estimated as I s τ tr , where τ tr is the intra-or inter-valley momentum relaxation time for the bulk or edge states' contribution, respectively.
Anomalous continuity equation-We consider a 2D crystal periodic in the x direction with period a = 1 and with the edges positioned at y = 0 and y = −W .A uniform electric field of magnitude E oscillating at frequency ω is applied along the x direction.For the sake of conciseness, hereafter we set ℏ = 1.Thus the conductance quantum e 2 /h is equal to e 2 /(2π), where e is the electron charge.From the Kubo formula [24,25], the y component of the valley current (averaged over x) is [26] and the valley density (also averaged over x) where is the Berry connection.The Fourier transform of Eq. (1) follows directly [26] from Eqs. ( 3) and (4): The vanishing of valley density accumulation-Let us first assume that the system is a fully gapped timereversal invariant insulator.The first term on the right hand side of Eq. (4) vanishes because ∂ k f k,n = 0, since there are no bands that cross the Fermi level.In [26] we show that, due to time-reversal symmetry, the second line on the right hand side of Eq. ( 4) is proportional to ω, so the valley density accumulation vanishes in the limit of static electric field.This result implies that ∂ y j y v (y) can be different from zero-as it must necessarily be, since the valley Hall current is finite in the bulk but vanishes at the edges-yet this finite divergence does not cause any density change at the edge or anywhere else.The resolution of this apparent paradox is provided by the anomalous term on the right hand side of Eq. ( 1) which exactly matches the divergence term on the left-hand side when the system is gapped.The undergap current does not produce a density accumulation.
The source of valley density-Let us now consider the case in which some energy levels cross the Fermi level.The first term on the right hand side of Eq. (4) causes the density to grow at a constant rate, leading to a breakdown of linear response theory unless a limiting momentum relaxation mechanism, such as intra-or inter-valley scattering, is taken into account.The Fermi surface term, obtained by multiplying Eq. ( 4) by −iω and taking the ω → 0 limit, is the "source term" Q s (y) in Eq. (2).As discussed above, the integral of Q s (y) over y across a single edge can be interpreted as an effective current I s that feeds the density accumulation thereat.I s has contributions from both edge and bulk states that cross the Fermi level.The latter give (for the edge at y = 0) where momentum integration is restricted to the valley with valley number +1, p is momentum in the y direction measured from the valley bottom, v λ k,p are envelope amplitudes of propagating stationary states, labelled by index λ, R λ (k, p) is the reflection probability amplitude (|R λ (k, p)| = 1) (see [26] for details).
Example: "gapped graphene"-To demonstrate the main features of the general theory developed above, we calculate the valley Hall current and the valley density accumulation rate for a nanoribbon of "gapped graphene"-a model system that captures some aspects of monolayer graphene on a gap-inducing hBN substrate.Lattice sites are labelled with a unit cell number l and a composite index (m, σ), where m = 1, . . ., N denotes the row, while σ = A, B distinguishes the sublattice within a given row.The y coordinate will be assumed to take integer values to mark the position within a row and halfinteger values to mark the position in between the rows.The two sublattices, A and B, have different on-site potentials ±∆.Electrons are assumed to hop only between nearest neighbors.We neglect spin-orbit interaction and therefore consider spinless electrons.
For the nanoribbon we consider two terminations: a) zig-zag boundaries on both edges [Fig.2(a)] and b) a zig-zag and a bearded edge [Fig.2(c)].These lattice terminations ensure that the valley number is conserved by the unperturbed Hamiltonian.Each unit cell consists of N horizontal rows, with two atoms in each row as shown in Fig. 2(a), except the edge rows, where one atom may be missing as shown in Fig. 2(c).
The band structures for the two terminations, shown in Figs.2(b) and (d), respectively, feature two bands separated by a gap equal to 2|∆| with minima at k = ±2π/3.These points define the two valleys in the onedimensional Brillouin zone.When the lattice is terminated with zig-zag boundaries on both edges, two dispersionless bands of edge states [blue lines in Fig. 2(b)] connect the two valleys and become bulk states for |k| < 2π/3.The upper (lower) band of edge states resides on the upper (lower) edge in Fig. 2 (a).
Our main results are presented in Fig. 3.For a Fermi energy in the gap (ε F = 0) and at zero temperature, we find that n v (m, 0) = 0 for either termination, consistent with the fact that there are no states at the Fermi level.At the same time j y v (m+1/2, 0) = −Ee 2 •sign(∆)/(2π)+ O(∆/t) for 1 ≤ m ≤ N − 1 as shown in panel (a), blue line: this is the undergap current associated with the nearly-quantized Hall conductance (the actual value −0.9 deviates from the ideal quantized value −1 due to the finite bandwidth of the model) [26].
When the system is doped with electrons (ε F = 0.3t), the current distributions differ dramatically for the two terminations, as shown by the red lines in Figs. 3 (a) and (b).In the case of the double zig-zag termination the current shows a linear variation across the ribbon (red line in (a)), changing sign about the center of the ribbon.This behavior is completely at odds with our intuition, which would lead us to expect an approximately constant current in the bulk, but not entirely unexpected, because there is no scattering and electrons propagate ballistically.Of greater physical interest, however, is the valley density accumulation rate which is shown in Fig. 3

(c).
There is a significant cancellation between −∂ y j y v (green line) and the non-conservation term (black line) at the edges.The sum of the two results in a density accumulation rate which displays oscillations (red dots) on the scale of half the Fermi wavelength and two spikes of equal signs at the edges.These are the result of interference between the electronic waves incident on and reflected off the edge.The fact that the accumulation rate does not integrate to zero is the result of the anomaly on the righthand side of Eq. ( 5): valley number is pumped from one valley into the other via a partially filled band of edge states connecting the two (upper blue line in Fig. 2 (b)).This opens the way to an intriguing possibility of generating a net valley density polarization by purely electrical means, as opposed to the standard optical methods.Notice, however, that the form of valley density accumulation rate cannot be predicted from the valley Hall current alone and depends on the boundary conditions.The zig-zag+bearded termination presents us with a more familiar scenario.Panel (b) of Fig. 3 shows that the valley Hall current is approximately constant (−0.55 in units of e 2 E/(2π) at ε F = 0.3t) in the bulk.At the same time the valley density accumulation rate, presented in panel (d) has spikes of opposite signs on the two edges.This suggests a more conventional picture of valley density being transported from one edge to the other.The reason for the overall valley number conservation is, in contrast to the previous example, absence of any partially filled bands that would connect the two valleys.
Conclusion-The modified continuity equation ( 1) al-lows us to explain how a non-vanishing undergap valley current can coexist with a vanishing valley density accumulation in a fully gapped non-topological timereversal-invariant system with perfectly degenerate valleys.Any valley density accumulation requires the existence of states at the Fermi level and furthermore it is a dissipative process which requires a scattering mechanism to reach a steady state.We have provided closed expressions for calculating valley density accumulation rates on the edges of a two-dimensional material and we have applied them to the gapped graphene model: these formulas show that the connection between bulk currents and measurable edge accumulations is much more complex than previously suspected.This, in particular, leads us to surmise that any physical system in which evidence of the VHE has been found either by Kerr rotation microscopy [11]

I. THE SETUP
We consider a 2D system periodic in the x direction and of finite width W in the y direction.We will take the period in the x direction to be equal to one.We divide the system into unit cells, labelled by index l.We assume that the periodic direction is chosen such that the valley index remains a good quantum number, i.e., edges are chosen to be parallel to the separation between the valleys in momentum space (one of the edges will be taken to be at y = 0 and the other at y = −W ).The Bloch eigenstates are denoted as ψ k,n (x, y) = u k,n (x, y)e ikx / √ 2π, where k is the wavevector in the x direction and n is the band index.To take account of the spin-orbit interaction, we will assume that both ψ k,n and u k,n are two-dimensional spinors.We apply the electric field of magnitude E in the direction parallel to the edges (i.e., the x direction), i.e., the Hamiltonian is perturbed by a potential term eE x, where −e is the electron charge and x the position operator.
For graphene nanoribbon with zig-zag and bearded edges (which preserve the valley number), the role of coordinates (x, y) is played by unit cell number l and combined index (m, σ), where m = 1 . . .N denotes the two-atom horizontal row in each unit cell and σ the atom (A or B) within it (see Fig. 1 (a), (c) in the main text).Note that a row may miss an atom of either sublattice, as in Fig. 1 (c), where it misses an A atom in row 1 of each unit cell.Everywhere below when displaying the results for densities and currents in the nanoribbon, the position on the y axis will only be resolved down to the row number m.Therefore, the y coordinate will be replaced by a discrete index that will take integer values marking the position in a certain row and half-integer values marking the position (half-way) between the rows.

II. POSITION OPERATOR IN BLOCH STATE BASIS
The derivation of the position operator representation in the Bloch state basis here follows closely Ref. [1].Consider a wavepacket χ(x, y) that is a superposition of Bloch eigenstates with coefficients χ k,n The application of the position operator x to χ(x, y) We integrate by parts in the first term, while discarding the boundary contribution.Furthermore, we expand i∂ k u k,n in the series of functions u k,n so that with As a result, we obtain (lightening the notation by suppressing x and y dependence) where

III. VALLEY NUMBER AND VALLEY HALL CURRENT OPERATORS
Here everything refers to just one electron, the generalization to many electrons is trivial.The operator that gives valley density at position x = (x, y) has the form nv (x, y) = − e 2 n(x, y), S( k) , where n(x, y) = |x, y⟩⟨x, y| is the particle density operator at position (x, y) (with any spin polarization) and the curly brackets stand for an anticommutator.Note that the particle density operator satisfies the continuity equation with a particle number current where v is the velocity operator.Operator S( k), on the other hand, is a constant of motion.Therefore, taking the time derivative of Eq. ( 6) at time t, we obtain Note also that n(x, y) = δ(x− r), where r is the electron's position operator.

IV. KUBO FORMULA FOR VALLEY NUMBER AND VALLEY HALL CURRENT
The Kubo formula [2,3] allows one to calculate the change in the ensemble average of an observable caused by a perturbation to first order in its strength.If the Hamiltonian is perturbed by an operator B exp(−iωt) then, long after the perturbation was turned on, the ensemble average of the observable Â will be oscillating at the frequency ω with the amplitude Here, α and β label eigenstates of the unperturbed Hamiltonian, whose eigenvalues are ε α and ε β , respectively.Furthermore, A αβ and B βα are matrix elements of operators Â and B between eigenstates α and β, while ) −1 are the occupation numbers of such states.For the nanoribbon, the unperturbed stationary states are labelled by the value of (quasi)momentum k and by the band index n.In the case of a harmonic electric field applied parallel to the edge, B = eE x.Using Eq. ( 5), the matrix elements of the perturbation are then given by the equation The role of observable Â is played by either nv (x, y) or ĵv (x, y), depending on the response function under consideration.Matrix elements of operator nv (x, y) between stationary states are where we have used the fact that, by definition, n(x, y) = |x, y⟩⟨x, y| and ⟨x, y|k, n⟩ = ψ k,n (x, y).For k = k ′ , Eq. ( 13) becomes [n v (x, y)] kn;kn ′ = − e 2π u † k,n (x, y)u k,n ′ (x, y)S(k).( 14) Matrix elements of operator ĵv (x, y) are given by At Note that due to Eq. ( 7) the following identity is true Plugging into this equation n(x, y) = |x, y⟩⟨x, y| and Let us integrate this equation along the length of one unit cell and using the periodicity of [j(x, y)] kn,kn ′ in the x direction discard the boundary terms.The result will be Now using the fact that ĵy must vanish at the boundary of the strip at y = 0, we can integrate this equation with respect to y to obtain From Eqs. ( 16) and (20) we can now obtain that In calculating the linear response, we will average over the length of the unit cell in the x direction, so that Eq. ( 21) will turn out to be useful.The linear response formula for the valley Hall current, averaged over the length of the unit cell, has the form Integrating by parts the first term in the round brackets and then evaluating the integral with respect to k ′ will give The first term vanishes because the unit cell average 1 0 dx [j y v (x, y)] kn,kn in that term vanishes, see Eq. ( 21) at n = n ′ .Plugging Eq. ( 21) into the second line, we obtain which gives Eq. ( 3) of the main text.Next, let us calculate the valley number response.According to Eq. ( 11), it is given by the equation, averaged across the length of the unit cell, Integrating by parts the first term in the round brackets and using Eq. ( 14), we obtain which is Eq. ( 4) in the main text.In the limit of zero frequency, if there is a partially filled band such that ∂ k f k,n does not vanish, the first term gives the dominant contribution.One can demonstrate that if the system is gapped and time-reversal symmetry is not broken, the second term on the right-hand side of Eq. ( 26) is of order O(ω) and, therefore, vanishes.This fact, which is demonstrated in Sect.VIII below, was used in the main text to discard its contribution to the valley density accumulation.

V. THE VALLEY DENSITY RESPONSE IS GAUGE INVARIANT
Let us demonstrate in this section that the valley density response, Eq. ( 4) in the main text, is gauge invariant.The first line of Eq. ( 4) is obviously gauge invariant so let us turn to the second line.This is proportional to Let us change phases of all the stationary states as follows u k,n (x, y) → e −iϕ k,n u k,n (x, y). (28) After this the Berry connection A k,n ′ n changes as Using the normalization condition this can be rewritten as Let us now plug the modified states (28) into Eq.( 27).This will result in the following change The second term in the round brackets multiplies f k,n − f k,n ′ and disappears.The phases in the rest of the equation cancel each other.Therefore, P (ω + i0) will not change, which means that it is gauge invariant.

VI. PROOF OF THE CONTINUITY EQUATION
In this section we demonstrate the validity of Eq. ( 1) in the main text.Let us multiply Eq. ( 26) by −iω.By Let us work on the sum in the second line of this equation, Let us divide this sum in two across the minus sign and relabel n to n ′ and vica versa in the second sum.The result reads

Now let us plug into this equation the definition for
The result will read Let us first perform summation over n ′ by using the completeness relation where 1 2 is the identity operator in the spin space, and then perform integration over x ′ and y ′ .This will lead to Joining the two sums together we obtain Plugging this into Eq.( 32), we obtain We combine the first two lines to obtain From Eq. ( 24), the divergence of the current is obtained to be which is the term on the second line of Eq. ( 41).Combining Eq. ( 41) and (42), we obtain Eq. ( 1) in the main text.

VII. LIGHTNING SPEED DERIVATION OF THE CONTINUITY EQUATION
The continuity equation derived in the previous section can be easily derived using Heisenberg equations of motion.Everywhere below operators are given in the Heisenberg picture.The equation of motion for the valley density has the form Recall that in the absence of the electric field valley density satisfies a continuity equation, i.e., i Ĥ(t), nv (r, t) = −∇ ĵv (r, t), with the current defined as in Eq. ( 10).Now using definition ( 6) and the fact that where the prime indicates a derivative, one can obtain Thus the equation of motion for the valley density has the form Replacing all the operators in this equation with their many-body versions and then calculating the ensemble average to first order in E will give Eq.( 1) in the main text.

VIII. TIME REVERSAL SYMMETRY IMPLIES NO STATIC VALLEY POLARIZATION
Consider Eq. ( 4) in the main text and take the second line from that equation.It is proportional to Assume presence of time-reversal symmetry which acts on the electron's wavefunction as For the Bloch states this implies that or where ñ labels another stationary state which satisfies ε k,n = ε −k,ñ and α k,n is a possible phase factor.Plugging Eq. ( 51) into the connection A k,n ′ n in Eq. ( 48), we obtain where superscript t stands for 'transposed'.Interchanging the electron wavefunctions in the bilinear products and using the normalization condition, this equation can be rewritten as Analogously, Let us now substitute Eqs. ( 53) and (54) into Eq.( 48).We obtain The term in the round brackets proportional to δ n,n ′ gives a vanishing contribution to the sum and overall all the phases disappear.Let us also take into account that Let us make use of these identities to make the summand only depend on ñ and ñ′ and then also note that summation over n and n ′ is the same as summation over ñ and ñ′ .Changing the summation variables from the former to the latter and then relabelling these back to n and n ′ , we obtain Changing the integration variable from k to −k, relabelling the summation indices from n to n ′ and vice versa, and using the fact that S(−k) = −S(k), we obtain If the system is fully gapped, at small frequency the i0 prescription is irrelevant and can be neglected so that Eq. (57) implies that P (ω) = −P (−ω), which means that the second line of Eq. ( 4) is of order O(ω) and can be neglected at zero frequency.

IX. VALLEY HALL CURRENT AT VANISHING FREQUENCY
At vanishing frequency, the following relation holds (compare Eqs. ( 24) and ( 33)): Plugging the result (39) into this equation, we obtain

X. EFFECTIVE VALLEY CURRENT
Let us try to evaluate the integral in the thermodynamic limit, i.e., for the ribbon width W → ∞.Assume that the chemical potential is in the conduction or valence band.Due to the factor ∂ k f k,n the integral over k and sum over n are restricted to the Fermi surface which consists of two disjoints parts, one in each valley.Due to time reversal symmetry the two parts give equal contributions to the integral, so we can calculate just one and multiply the result by two, i.e., Everything from now on refers to the valley with valley number +1.Assuming that in the sum over n and integral over k we never go too far away from the bottom of the valley, we can use the envelope wave function description for the stationary states.Suppose energy eigenstates for a system without boundaries are described by a set of multicomponent valence (conduction) band envelope amplitudes v λ k,p with energies ε λ (k, p), where p is the component of momentum along y measured from the bottom of the valley and λ is a discrete label counting the stationary states.Now let us introduce the boundaries at y = 0 and y = −W .Assume that the scattering off the boundaries does not mix eigenstates with different values of λ and that ε λ (k, −p) = ε λ (k, +p).Then the valence (conduction) band envelope wave-functions for the system with boundaries will have the form which is nothing but a sum of an incident and a reflected wave.Here R λ (k, p) is the probability amplitude for scattering off the boundary at y = 0, momenta p m take discrete positive values (labelled by m), with the distance between them approaching π/W as W → ∞, and an overall factor N λ (k, p m ) is chosen such that In terms of the envelope wave-functions the expression for Q s (y) takes the form Plugging Eq. (62) into this equation, we obtain where we used the normalization condition |v λ k,p | 2 = 1 and the fact that |R λ (k, p)| = 1.Let us consider values of y such that |y| ≪ W , which means that we are keeping very close to the edge at y = 0.In this case each term in the sum over m is not much different from the next one and we can exchange the sum over m for an integral with respect to p.More precisely, to evaluate the sum over m we can use the Euler-Maclaurin formula, which will give us an expansion in powers of 1/W and the leading-order term is obtained by simply exchanging the sum over m for an integral over m.This in turn can be exchanged for an integral over p via p = πm/W .Taking into account Eq. ( 64), this will lead to where we neglected the terms of order O(1/W ).Contribution of the first term in the outer round brackets in the second line to the integral is equal to zero as vanishes at the values of k far enough from the bottom of the valley.Therefore we are left with By Riemann-Lebesgue lemma, Q s (y) → 0 (up to terms of order O(1/W ) that we discarded) as y → ∞.Because integration in Eq. ( 68) is restricted to the Fermi surface, it is clear that Q s (y) oscillates in space at a wavelength corresponding to double the Fermi momentum.We can now introduce an "effective current" feeding valley number accumulation at the edge XI. GRAPHENE NANORIBBON.

DESCRIPTION OF THE STATIONARY STATES
In this section and the next we prove analytically that in the completely gapped state the valley Hall current for a graphene nanoribbon with zig-zag edges is non-zero and quantized.The treatment here closely follows that of Ref. [4].Because the spin-orbit coupling in graphene is small, we will neglect it completely.In this case spin polarization is a good quantum number and its only effect is the introduction of a degeneracy factor of 2 in all linear response formulae.We will avoid that by considering a spinless electron from now on.For a spinful electron the results given here will have to be multiplied by a factor of 2.
As promised in Section I, the pair of coordinates (x, y) will be replaced by unit cell number l and combined index (m, σ), where m = 1 . . .N and σ = A, B. Integrals over x ′ and y ′ within the unit cell will be replaced as follows Furthermore, in these expressions, the current that flows between rows m and m + 1 is denoted as j y v (m + 1/2, ω).The tight binding Hamiltonian for graphene nanoribbon with zig-zag edges in second quantized form reads where t is the nearest neighbour hopping parameter and âl (m) and bl (m) destroy an electron on sublattice A or B, respectively, in a two-atom row m of unit cell l.They satisfy the usual anticommutation relations We introduce the Fourier-transformed operators αk (m) and βk (m) from the relations where x l,m,A and x l,m,B are the positions of the atoms of the A and B sublattices in the direction along the ribbon and k takes values in the interval (−π, π].The Fouriertransformed creation and annihilation operators satisfy the anticommutation relations In terms of Fourier transformed creation and annihilation operators, the Hamiltonian reads where g k = 2 cos(k/2).From now on, we will set t = 1.The one-particle eigenstates of the Hamiltonian are Bloch states ψ k,n (l, m, σ) = u k,n (m, σ) exp (ikx l,m,σ )/ √ 2π.As is customary, we will combine the amplitudes u k,n (m, σ) for σ = A, B into a sublattice pseudospinor u k,n (m) = (u k,n (m, A), u k,n (m, B)) t .For bulk states, this takes the form where ε kps = s ∆ 2 + g 2 k + 2g k cos(p) + 1 and the role of the band index n is played by the combination ps, with s = ± and p ∈ (0, π) the solution of the equation In Eqs.(77)-(78), j = 1, 2, . . ., N for g k > N/(N + 1) and j = 1, 2, . . ., N − 1 for g k < N/(N + 1).In Eq. (77) the normalization constant equals We note that, for g k < N/(N + 1), apart from the states described by Eqs. ( 77)-( 78), there is also an edge state where ε kes = s ∆ 2 + g 2 k − 2g k cosh(η k ) + 1, with s = ±, and η k > 0 is the solution of the equation The role of the band index n here is played by the combination es, where e stands for 'edge' and s is described above.The normalization factor equals XII. CURRENT IS QUANTIZED FOR THE TOTALLY GAPPED GRAPHENE NANORIBBON Let us calculate the valley Hall current for the totally gapped graphene nanoribbon with zig-zag edges in the limit of vanishing frequency.Using Eq. ( 59) (and neglecting spin) we obtain where n stands for either ps or es, as explained in Sect.XI.Note that explicit summation over σ here is replaced by matrix multiplication of the hermitian conjugate of a pseudospinor with itself.When the chemical potential is in the gap, f kps = f kes and equals 1 for s = −1 and 0 for s = +1.For the evaluation of the integral with respect to k note that for any n holds u k,n = u −k,n and we can take S(k) equal to 1 for 0 < k < π and −1 for −π < k < 0. Thus keeping only the occupied states in the sum over n and evaluating the integral with respect to k, we obtain The states u k,n form a complete set, therefore where 1 2 is the 2 × 2 identity matrix in the sublattice space.Using the explicit form of the wavefunctions [see Eq. (77)] one can find that Using Eqs. ( 85) and (86) one can rewrite the sum in the round brackets in Eq. (84) in the form Plugging this into Eq.( 84) and taking into account the fact that there are no edge states at k = 0 we obtain Let us estimate the sum over p in Eq. ( 88).Since by Eq. ( 78), the factors remain of order O(N 0 ) at k = 0 or k = π, one can see that each term in the sum in the second line of Eq. ( 88) is of order ∆/N .In making this estimate we also took into account that ε kp− is of order one [in units of t] at k = 0 and k = π.Therefore the whole sum over p in Eq. ( 88) is of order ∆.Hence the current density response equals (restoring the hopping parameter t) For edge states, as k approaches π, g k → 0 while η k approaches positive infinity as η k ∝ − ln(g k ).In this limit ε kes → s|∆| and the probability distribution for the edge states has the form Plugging this into Eq.( 91), for 1 ≤ m ≤ N − 1 we obtain Let us also, for future reference, quote here the result for the divergence of the current, div j y v (m, 0) = j y v (m − 1/2, 0) − j y v (m + 1/2, 0), which follows from this equation

XIII. QUALITATIVE EXPLANATION OF THE VALLEY HALL CURRENT PROFILE FOR DIFFERENT BOUNDARY CONDITIONS
Throughout this section we assume that the Fermi level lies in the conduction band and takes a value in the interval |∆| < ε F < t.Let us first consider the zig-zag nanoribbon.Expression for the current, Eq. ( 83), can be written in an alternative form as j y and and we set The integral with respect to k can be easily calculated to produce div T The sum in the second and third lines on the right hand side of this equation (including the factor of −2) has already been calculated in the previous section, compare Eq. ( 84).It is given by whatever multiplies e 2 E/(2π) in Eq. ( 94).The term in the first line is the contribution of the upper band of edge states, which has appeared because we raised the Fermi energy and this band became occupied.So, using Eq. ( 92) and Eq. ( 94), we obtain (ignoring terms of order O(∆/t)) Note that div T changes extremely fast on the boundaries (it goes from −1 to zero on the scale of one inter-atomic distance) and does not change at all inside the ribbon.This can be traced back to the fact that ultimately the spatial behavior of T (m) is governed by the localized edge states (see Eq. ( 91) and Eq. ( 99)).(101) Let us point out that div G(m), as opposed to div T (m), changes slowly in space, because its behavior is governed by the states on the Fermi surface.Therefore it varies significantly on the length scale defined by the inverse Fermi momentum ( ε 2 F − ∆ 2 ) −1 ≫ 1 for ε F close enough to the bottom of the conduction band.Therefore the roughest (but only the roughest) estimate of div G(m) can be given by just the spatial average ⟨div G(m)⟩ sp = (1/N ) m div G(m).This is not zero.Indeed, using m u † k,n (m)u k,n (m) = 1 we obtain which is non-zero because of the left-mover-right-mover imbalance in each valley.Indeed, using f k,n = θ(ε where the sum runs over all values of k at which the Fermi level crosses an energy band and v i is the group velocity of the band crossed at a point k i .Because there is one more right-mover than there are left-movers in the left valley and one more left-mover than there are rightmovers in the right valley, see Fig.To get j y v (m + 1/2, 0) we need to integrate (i.e., sum over m) Eq. (100) and Eq. ( 104) with the boundary conditions that the current vanishes outside the ribbon.It is not difficult to observe that Eq. (100) determines the one-sided limiting values of the current on the boundaries as approached from within the ribbon and Eq. ( 104 Now let us briefly discuss the case of the nanoribbon with a bearded edge.Because there is only one band of edge states whose occupation number does not change as we raise ε F , we see that div T does not change compared to the undoped case.This means that the values of the current on the boundaries will stay roughly the same as in the undoped case.Next, because there is no left-moverright-mover imbalance in the valleys, the overall average slope will be zero.This very rough analysis is confirmed by the numerical results, see Fig.  bearded, see Fig. 1.In an infinite system this is predicted to be fixed and quantized when the chemical potential is in the band gap and to go down to zero as −∆/|ε F | (in units of e 2 E/(2π)) for |ε F | > |∆|.For the nanoribbon the behavior of the curve is a bit different, see Fig. 1.It is indeed fixed and quantized (up to corrections of order O(∆/t)) when the Fermi energy is in the gap but outside the gap it does not conform to the ∆/|ε F | law and settles on a value of around ±1/2 for high enough hole or electron doping.
At ε F = −|∆| the valley Hall current has a discontinuity due to all the edge states suddenly changing their occupation numbers.The oscillations below ε F = −|∆| are a finite size effect due to small values of the Fermi momentum.

FIG. 1 .
FIG. 1. Cyclic flow of the Bloch wave vector under the action of an electric field in the x direction.In both panels, the valley charge is −1 in the red regions and +1 in the green regions.Also shown are the Berry curvature hot spots with positive value near K and negative value near K ′ .Panel (a) shows the flow in a one-dimensional ribbon where the valley charge is non-conserved (changes sign) when k crosses the origin of momentum space.At the same time the Berry curvature, and hence the anomalous velocity also changes sign ensuring the existence of a steady valley Hall current.In a fully gapped insulator the flow does not alter the occupation of the states (i.e., a full band remains full) so there is no change in the valley density.Panel (b) shows examples of flows in a twodimensional periodic system.In both cases the assumption of time-reversal symmetry rules out the possibility of a valley density response arising from the magneto-electric polarizability of the Bloch states.

FIG. 2 .
FIG. 2. Panel (a) and (c): Gapped graphene nanoribbon with zig-zag boundaries on both edges and a zig-zag and a bearded edge (at the top), respectively.Red (blue) discs signify atoms of the σ = A (σ = B) sublattice, l labels the unit cells, m the rows in each unit cell.Panel (b) and (d): Band structures of nanoribbons of panels (a) and (c), respectively, for N = 20 and ∆ = 0.2t.Non-topological edge states are depicted by blue lines.

FIG. 3 .
FIG. 3. Panel (a): Gapped graphene nanoribbon with zigzag edges: the blue dashed line shows the valley Hall current as a function of position at Fermi energy εF = 0 and the red solid line shows the same quantity at εF = 0.3t Panel (b): Same as in (a) for a nanoribbon with one edge zig-zag and the other bearded.Panel (c): The green line shows the valley density accumulation rate contributed by the valley Hall current in the doped nanoribbon with zig-zag edges.The black line shows the contribution of the non-conservation term, and the red dotted line is the sum of the two, i.e., the total accumulation rate.Panel (d): Same as in (c) for the nanoribbon with one zig-zag and one bearded edge.In all plots N = 100, ∆ = 0.1t.In plots (c) and (d) εF = 0.3t.
where a factor of 2 appeared because the contribution of the left valley to the integral equals the contribution of the right valley, so we left only the latter and multiplied it by two.Plugging in the occupation numbers, we obtain div T(m) = − 2 u † πe+ (m)u πe+ (m) − 2 p u † kp− (m)u kp− (m) + u † ke− (m)u ke− (m) Consider now div G(m), which is given by the equationdiv G(m) = n dkS(k)(∂ k f k,n )u † k,n (m)u k,n (m).
1(b) (the upper blue line) in the main text, we obtain ⟨div G(m)⟩ sp = ) its overall slope as a function of position.It then follows that the current j y v (m + 1/2, 0) is approximately equal to Ee 2 /(2π) at m = 1, to −Ee 2 /(2π) at m = N − 1 and those two values are connected roughly by a straight line with the slope −[e 2 E/(2π)](2/N ).

2
FIG. 1.Valley Hall current as a function of the Fermi energyin graphene nanoribbon with one edge zig-zag and the other bearded.The value is taken in the middle of the ribbon, the width N = 100, ∆ = 0.1t.