# Valley-dependent transport through graphene quantum dots
due to proximity-induced, staggered spin-orbit couplings

A. Belayadi<EMAIL_ADDRESS>Dept. of Theoretical Physics, University of
Science and Technology Houari Boumediene, Bab-Ezzouar 16111, Algeria.  École
Supérieure des Sciences de l’Aliment et Industries Alimentaires, ESSAIA, El
Harrach 16200, Algeria.  P. Vasilopoulos<EMAIL_ADDRESS>Dept. of
Physics, Concordia University, 7141 Sherbrooke Ouest, Montréal, Québec H4B
1R6, Canada  N. Sandler<EMAIL_ADDRESS>Dept. of Physics and Astronomy, Ohio
University and Nanoscale Quantum Phenomena Institute, Athens, Ohio, USA.

###### Abstract

We study a system composed of graphene decorated with an array of islands with
$C_{3v}$ symmetry that induce quantum dot (IQD) regions via proximity effects
and gives rise to several spin-orbit couplings (SOCs). We evaluate transport
properties for an array of IQDs and analyze the conditions for realizing
isolated valley conductances and valley-state localization. The resulting
transmission shows a square-type behavior with wide gaps that can be tuned by
adjusting the strength of the staggered intrinsic SOCs. Realistic proximity
effects are characterized by weak SOC strengths, and the analysis of our
results in this regime shows that the Rashba coupling is the important
interaction controlling valley properties. As a consequence, a top gate
voltage can be used to tune the valley polarization and switch the valley
scattering for positive or negative incident energies. A proper choice of SOC
strengths leads to higher localization of valley states around the linear
array of IQDs. These systems can be implemented in heterostructures composed
of graphene and semiconducting transition-metal dichalcogenides (TMDs) such as
MoSe2, WSe2, MoS2, or WS2. In these setups, the magnitudes of induced SOCs
depend on the twist angle, and due to broken valley degeneracy, valley
polarized currents at the edges can be generated in a controllable manner as
well as localized valley states. Our findings suggest an alternative approach
for producing valley-polarized currents and propose a corresponding mechanism
for valley-dependent electron optics and optoelectronic devices.

##  I introduction

The study of valleytronic materials is of significant importance in the area
of information processing and encoding [1, 2, 3] due to the alternative degree
of freedom furnished by the carrier’s valley momentum in addition to the
conventionally used charge and/or spin properties.

The benchmark material is a two-dimensional (2D) graphene-based structure with
Dirac cones at ${\bf K_{1}=-K}$ and ${\bf K_{2}=+K}$ valleys, proposed as a
strong candidate for future valley-driven computing devices through the
manipulation of valley currents [4, 5, 6], i.e., by applying external
voltages. However, the lack of external probes or contacts that can select
individual valley currents as ferromagnetic contacts separate spins polarized
currents in spintronic devices [7, 8, 9], remains a primal obstacle to
encoding and information processing through the valley index.

In addition to transport, the material’s optoelectronic properties are also
used to access the valley degree of freedom [10, 11, 12, 13]. Several
strategies have used the optical response to control, detect, and monitor
valley polarization [14, 15, 16, 17]. Generally, a combination of gates
voltages -implemented via scanning tunneling microscopy- and suitable
substrate magnetic materials bring out the mechanism that tunes the desired
electronic, spintronic, and valleytronic properties [18, 19, 20, 21] as
proposed and later demonstrated by polarization-resolved photoluminescence
experiments [10, 22, 13, 23, 24].

An alternative approach to induce valley separation involves exploiting
confined geometries. Quantum dots (QDs) can produce valley-filtered currents
and are important ingredients in modern nanotechnology devices. Typically,
confined geometries that induce valley separation are obtained via a wide
variety of methods that include electrostatic confinement produced by a
scanning tunneling microscopy tip [20], strain fields [25], bilayer graphene
structures with spatially varying broken sublattice symmetry [26] and isolated
regions defined by local broken sublattice symmetry [23]. In all these setups,
valley separation is achieved because of the effects of external fields or due
to substrate properties that are extremely difficult to design and control
with sufficient precision. As a consequence, the potential of these geometries
to induce selective valley filtering and confinement in a controllable manner
and without external fields remains untapped. To address this issue, we
investigate the properties of a proposed heterostructure that exploits
proximity effects and periodic spin-orbit interactions.

Ideally, the most efficient way to induce uniform and large staggered spin-
orbit couplings (SOC)s on graphene is via proximity to appropriate substrates
that break the sublattice symmetry, thus allowing for a clear distinction of
the two pseudospins. Recently, the role of SOCs in valley separation has been
addressed in several works, such as graphene deposited on top of hexagonal
boron nitride [27] and graphene/TMD heterostructures [28, 29]. These setups
possess staggered onsite potentials that give rise to various SOCs via
different mechanisms [23]. Interestingly, not all types of SOCs will render
valley separation as shown by the sublattice independent intrinsic spin-orbit
coupling (ISOC) in the Kane and Mele model [30], or the Rashba SOC (RSOC) that
appears in the presence of external fields. However, other appropriately
engineered interactions can break the sublattice symmetry, rendering two main
effects that we refer to as (1) the rise of a staggered potential with a
concomitant gap opening and (2) the emergence of an ISOC in a staggered form
that is sublattice dependent (i.e., sublattice-resolved SOC). In this last
case, the spin-valley transport is due to the emergence of a valley-Zeeman
type of coupling, defined by the ISOC sign change between sublattices [31, 28,
32]. This valley Zeeman effect is of great interest because it induces a giant
spin lifetime anisotropy in proximitized graphene [33]. Furthermore, Frank et
${\it al}$ [34] showed that in narrow-width cells of zigzag-terminated
graphene with a staggered ISOC, pseudo-helical and valley-centered states
(without topological protection) are localized along the edges. These results
are consistent with the bulk system’s topological invariant $Z_{2}=0$.

In this work, we propose a heterostructure composed of graphene and TMD
islands that combines the effects of confinement and SOCs in a controllable
manner. The model is inspired by recent experiments reported in Ref. 51, with
graphene deposited on top of a TMD island that induces a local region with
various SOCs, i.e., an induced quantum dot (IDQ). Our proposal generalizes the
experimental setup to a periodic array of islands placed below or deposited on
the graphene membrane. The TMD islands preserve the underlying $C_{3v}$
symmetry of graphene and introduce SOCs in the electron dynamics via the
proximity effect. We analyze the conditions for selective valley state
confinement and the generation of valley currents under applied voltages for a
generic model that is later applied to specific material combinations.

The paper is organized as follows. In Sec. II, we briefly present the model
for a system composed of a linear array of induced quantum dots in graphene
created by proximity effects and including different emerging SOC terms. In
Sec. III, we present numerical results revealing effective mechanisms for
valley filtering and confinement. We apply these results to a series of
heterostructures composed of different materials with realistic parameters and
analyze the effect of relative twisting between the two materials. A summary
and conclusions follow in Sec. IV.

##  II Model and methods

As mentioned above, we propose to study a chain of quantum dots with $C_{3v}$
symmetry in graphene created by proximity effects due to TMD islands. The
choice of TMDs that conserve $C_{3v}$ symmetry is made to ensure the largest
values of induced spin-dependent couplings in the graphene membrane[28, 35].
Such engineered IQDs will exhibit pseudohelical and valley-centered edge
states with potential for device applications [34, 36]. The salient advantages
of such structures are: 1) longer localization lengths for valley states in
narrow ribbons, 2) valley Chern numbers and localization lengths independent
of RSOCs, and 3) gapless band structures. A schematic picture of the system is
shown in Fig. 1.

Figure 1:  Panel (a) displays the overall device, which consists of a 2D
graphene ribbon with zigzag boundaries decorated with semiconducting
transition-metal dichalcogenide (TMD) islands. (b) Zoom-in of a quantum dot
made of graphene and TMD. To visualize the emergence of different SOC terms,
we duplicate the graphene membrane to emphasize the lift of the spin
degeneracy. Each layer corresponds to a different spin component, with the
blue (red) membrane representing the spin-up (spin-down) population. The
induced SOCs combined with the underlying $C_{3v}$ symmetry give rise to
sublattice-resolved intrinsic couplings $\lambda_{I}^{(A)},\
\lambda_{I}^{(B)}$, denoted by red arrows, and pseudospin inversion-asymmetric
couplings $\lambda_{PIA}^{(A)},\ \lambda_{PIA}^{(B)}$ denoted by green arrows.
The Rashba $\lambda_{R}$ coupling is represented by purple arrows.

The deposition of adsorbates on graphene, or of graphene membranes on islands,
results in profound changes in the electronic structure that depend on the
locally conserved symmetries as described in Ref. [28]. The most important
effects are the emergence of i) an effective staggered potential due to the
broken reflection symmetry imposed by different orbital interactions
experienced by the carbon atoms in proximity to the different atomic species
of the TMD material and ii) several sublattice dependent next-nearest neighbor
hopping terms originated from the proximity-induced spin-orbit interactions.
In this case, the system can be described by an extension of the models by
Kane and Mele [30] and Haldane [37]. The Hamiltonian for the QD regions is
given by: [27, 28, 29, 38]

$\displaystyle H_{QD}$ $\displaystyle=-t\sum_{\langle
i,j\rangle}a_{is}^{\dagger}b_{js}$ (1) $\displaystyle+\sum_{\left\langle
i\right\rangle}$ $\displaystyle\Delta\left(\xi^{(A)}{\bf a}_{is}^{\dagger}{\bf
a}_{is}+\xi^{(B)}{\bf b}_{is}^{\dagger}{\bf b}_{is}\right)$
$\displaystyle+\Big{(}\frac{2i}{3}\Big{)}$ $\displaystyle\sum_{\left\langle
i,j\right\rangle\sigma,\sigma^{\prime}}\left(\lambda_{R}\ {\bf
a}_{i\sigma}^{\dagger}{\bf b}_{j\sigma}\right)\left[{\bf\hat{s}}\otimes{\bf
d}_{ij}\right]_{\sigma,\sigma^{\prime}}$
$\displaystyle+\Big{(}\frac{i}{3\sqrt{3}}\Big{)}$
$\displaystyle\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle\sigma}\nu_{ij}\left(\lambda_{I}^{(A)}{\bf
a}_{i\sigma}^{\dagger}{\bf a}_{j\sigma}+\lambda_{I}^{(B)}{\bf
b}_{i\sigma}^{\dagger}{\bf
b}_{j\sigma}\right)\left[{\bf\hat{s}}_{z}\right]_{\sigma,\sigma}$

where $t$ is the nearest neighbor hopping between sites $i$ and $j$ (note that
these are spin-preserving processes). $\Delta$ is the staggered potential
induced by the TMD islands in the dot region. This potential is sublattice-
dependent with $\xi^{(A)}=1$ ($\xi^{(B)}=-1$), rendering opposite signs for
the induced gaps, i.e., $\Delta^{(A)}=-\Delta^{(B)}=\Delta$ for sublattice A
(B). The Rashba interaction (RSOC) is expressed in terms of the coupling
$\lambda_{R}$. This coupling breaks the $z$ inversion symmetry while
exchanging the spin of different sublattices. Here, the $\textbf{d}_{ij}$
vector connects site $j$ to $i$. The terms $\lambda_{I}^{(A)}$ and
$\lambda_{I}^{(B)}$ represent the intrinsic SOCs (ISOC) between next-nearest
neighbors. These terms connect the same sublattices and spins in the clockwise
($\nu_{ij}=-1$) or anticlockwise ($\nu_{ij}=-1$) direction from site $j$ to
site $i$. Finally, the spin is denoted by the vector $\widehat{\textbf{s}}$
with components written in terms of Pauli matrices. It is worth restating that
the SOCs exist only within the QD regions; outside the system is described by
the Hamiltonian of pristine graphene.

Working with a Hamiltonian in reciprocal space is more convenient for studying
valley properties. The resulting effective Hamiltonian is obtained by
linearizing Eq. (1) around the ${\bf K_{1}}$ and ${\bf K_{2}}$ valleys labeled
below by the valley index $\kappa=-1$ and $\kappa=+1$, respectively. The final
expression is given in the form $H_{QD}=H_{k}+H_{\Delta}+H_{R}+H_{I}$ [34],
where:

$\displaystyle H_{k}$ $\displaystyle=$ $\displaystyle\hbar v_{F}\left(\kappa
k_{x}\sigma_{x}+k_{y}\sigma_{y}\right)s_{0},$ (2) $\displaystyle H_{\Delta}$
$\displaystyle=$ $\displaystyle\Delta\sigma_{z}s_{0},$ (3) $\displaystyle
H_{R}$ $\displaystyle=$
$\displaystyle\lambda_{R}\left(-\kappa\sigma_{x}s_{y}+\sigma_{y}s_{x}\right)s_{0},$
(4) $\displaystyle H_{I}$ $\displaystyle=$
$\displaystyle(\kappa/2)\big{[}\lambda_{I}^{(A)}\left(\sigma_{z}+\sigma_{0}\right)+\lambda_{I}^{(B)}\left(\sigma_{z}-\sigma_{0}\right)\big{]}s_{z}.$
(5)

The Fermi velocity $v_{F}$ is expressed in terms of the hopping $t$ as
$v_{F}=\sqrt{3}a_{0}t/2\hbar$ where $a_{0}$ is the lattice constant. The
pseudospin is denoted by the Pauli matrices $\sigma$, and $s_{0}$ denotes the
spin identity matrix.

While Eqs.(2-5) provide an intuitive picture of the effect of each SOC term on
the valleys, the results presented in the following sections are obtained by
combining the S-matrix formalism with the tight-binding model for a zigzag
terminated ribbon. This boundary condition preserves the valley quantum number
and the valley topological properties of graphene. We compute the valley-
polarized conductance with the Landauer-Büttiker approach:

$G_{\kappa}^{n,m}=(e^{2}/h)\left|S_{\kappa}^{n,m}\right|^{2},\qquad(n,\
m\equiv L,R).$ (6)

Here $S_{\kappa}^{n,m}$ is the scattering matrix element between left (L) and
right (R) leads for a given valley index ${\kappa}=\pm 1$ [39]. Thus, our
calculations exploit the formalism with valley-dependent local currents as
defined in Ref. [40]. Details regarding the computation of the valley
conductance and currents are presented in Appendices A and B.

##  III  Results

In this section, we present numerical results for the valley-polarized
conductance for a range of structures that contain from a single IQD ($n=1$)
to a chain ($n>1$) of IQDs. To emphasize the qualitative features resulting
from the competition among the different interactions in the model, we adjust
the parameters’ values accordingly to present the main findings. This
procedure is usually applied to identify the role played by the various
interactions[41, 34, 42]. We note, however, that for accurate setups, one
expects weaker values for SOCs from proximity effects, and we address this
situation in Sec.III.2.

###  III.1 Qualitative analysis of results

Following Refs.[30] and [34], we solve the model using the following ranges
for the various SOC parameters: $\lambda_{R}\leq 0.075t$ and
$\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}\leq\sqrt{27}0.06t$ in units of $t=1$
[34].

We consider $n$ symmetric quantum dots ($n=1,2,3,4$), with the same spin-orbit
parameters, arranged in a chain with the same radius $r_{0}=7$ nm and the
space in-between them set to $d=5r_{0}$. The circular geometry of the dots is
chosen to preserve the symmetries of the original model and provide maximum
localization of valley states, as discussed below. We set the incident energy
at $E_{F}=0.035t$ and consider a ribbon width $w=30$ nm, with zigzag
boundaries. For this range of ISOC values, the valley Zeeman effect is the
most important term as the topological invariant $Z_{2}$ remains trivial,
i.e., $Z_{2}=0$ and the bulk system is in the same topological phase as long
as the ISOC is staggered, i.e., $\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}$ [34,
36].

The reduced number of group symmetries of pristine graphene due to the
proximity of the TMD island is reflected by the large number of lower-symmetry
allowed SOC parameters. The effects of these couplings are observed in Fig, 2
(a, b), where the band structure near the two Dirac cones appears clearly
modified. Fig. 2 (a) shows results in the absence (a) and presence (b) of a
staggered potential. The changes include a newly open gap and additional edge
bands resulting from the ISOC. These linear bands represent edge states that
are unique to the staggered case $\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}$ while
they disappear in the uniform regime $\lambda_{I}^{(A)}=\lambda_{I}^{(B)}$.
Similar results have been found in Ref.[34] where more details about the
valley selection and its related sublattice occupation are discussed. Based on
these findings, a proximity effect that induces a staggered ISOC might lead to
the emergence of valley currents and the localization of valley-centered
sublattice polarized edge states. Indeed, these results manifest themselves in
transport, providing valley states as conducting channels, as shown in Fig.
2(b), an issue discussed in the coming subsections.

Figure 2: Energy bands for uniform $\lambda_{I}^{(A)}=\lambda_{I}^{(B)}$ (a),
and staggered $\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}$) (b) for a 30 nm wide
zigzag ribbon, $\Delta=0.1t$ and $\lambda_{R}=0.075t$. Panels (c) and (d) show
valley conductances vs intrinsic spin-orbit length $\lambda_{I}$ for several
IQDs with staggered ISOC. Narrow panels on the right correspond to zoom-ins
for each valley as a function of $\lambda_{I}$. Panels (e) and (f) correspond
to conductance vs inter-dot spacing $d$ and number $n$ of IQDs.

####  III.1.1 Valley-dependent conductance through proximity IQDs

Figure 2(c) and (d) display the conductance through a group of $n$ symmetric
quantum dots with staggered ISOC. One interesting observation is the
approximately square-type dependence revealing wide gaps that can be made more
pronounced by changing the strength of ISOC for $n\geq 3$ QDs. We observe that
$100\%$ ($0\%$) of the detected conductance results from the flow of electrons
through the $\kappa=-1$ ($\kappa=+1$) valley for positive incident energy
($E_{F}=0.035t$), while $0\%$ ($100\%$) occurs for negative incident energy
($E_{F}=-0.035t$). The opposite behavior is obtained for the complementary
valley. Interestingly, the gaps occur at different ranges of ISOC values, with
the emerging valley-polarized current being switched from one valley to the
other within the gap region by an appropriate change of $E_{F}$. Furthermore,
we observe a decrease in the widths of the gaps with increasing SOC strength,
which suggests they might vanish for high enough values of ISOC.

An analysis of Figs. 2 (c) and (d) reveal resonances in the transmittance at
around the value $\lambda_{I}=0.065t$, with the spectrum in the zoom-in panels
showing that the IQDs confine electrons with index $\kappa=+1$ ($\kappa=-1$)
at positive (negative) energy. Hence, the scattering through the IQD region
tends to zero accordingly. We observe that the number and sharpness of the
resonance weakly depend on the number of IQDs along the chain, as seen for
$n=1,2,3$ curves that present at least one resonant state each at similar
values of ISOCs. More details about the observed confinement are addressed
below in Sec.III.2.3, and we discuss these results, including realistic
parameters, in III.3.

Additionally, the conductance response shows several interesting
characteristics: (1) It exhibits an oscillating behavior that becomes more
pronounced with increasing $n$. In this case, the conductance oscillations
arise from mode mixing, and their number depends on the number $n$ of the
IQDs, as shown in the zoom-in of Fig. 2 (c) and (d). (2) The conductance
plateaux become better defined as $n$ increases ($n\geq 3$), with values
$G(\kappa=+1)=0$ and $G(\kappa=-1)=2G_{0}$ for panel (c) and
$G(\kappa=+1)=2G_{0}$ and $G(\kappa=-1)=0$ for panel (d).

These results suggest that the conductance plateaux are due to states that
become valley polarized at specific strengths of the staggered ISOC, e.g., in
the range $0.015t$ $\leq\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}\leq 0.04t$ as
shown in Fig.2 (c). Interestingly, when $E_{F}$ is negative, the valley
polarization is reversed, as shown in panel (d) within the same range of
values for $\lambda_{I}$. Consequently, the transmitted current can be made to
be valley polarized from either one of the two Dirac points ${\bf K_{1}}$ or
${\bf K_{2}}$ depending on the incident energy.

####  III.1.2 Dependence on coherent inter-dot electron transfer

Within the conductance gap regions in Fig. 2 (c) and (d), and depending on the
bias defined by the sign of the Fermi energy, the valley that appears in the
output with $T=1$ seems to be barely scattered by the IQDs irrespective of the
number of dots in the chain. Inversely, the states from the valley that appear
in the output with $T=0$ are strongly reflected by the IQDs even for the
shortest chain furnished by only one dot. To better understand these current
profiles, we plot in Fig. 3 the local valley current through a chain of three
IQDs. As the figures show, the transferred valley current is due to electron
scattering processes that involve inter-dot hoppings between neighboring dots,
as illustrated in the figure, that display uniform local densities everywhere
between the dots for both valleys (Fig. 3 (a) and (d)). Notice also that in
this case, the conductance is practically the same as for the zigzag
terminated graphene ribbon with one or more quantum dots, as depicted in Fig.
2 (f).

Figure 3: Local current mapping, within the gap regions
($\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}=0.025t$) shown in Fig. 2 (c) and (d),
for both valleys for a given bias voltage while panels (e) to (h) show results
for opposite biases. The left (right) panels are for positive (negative)
values of $E_{F}$. Black arrows indicate the direction of the bias-incident
current. Red and blue solid lines indicate the direction of corresponding
current flows.

The above observations raise the question of the inter-dot spacing $d$’s
influence on the valley-transmission output. We have proposed that the valley-
scattering processes for complete transmission occur in the presence of
effective inter-dot hopping between neighboring IQDs. To test this hypothesis,
we analyze the role of $d$ on the valley transmittance. As shown in Fig. 2
(e), the inter-dot hopping is present with an inter-dot space $d\geq 5r_{0}$.
We observe non-stable transmission signatures if $d\leq 5r_{0}$. In this
regime, the edges of the IQDs are closer to each other, and the transmittance
peaks show values less (more) than unity (zero) for the ${\bf K_{1}}$ (${\bf
K_{2}}$) valley. These might be related to mode mixing that impedes the
propagation of specific transverse modes due to edge-dot coupling.
Interestingly, when the $d$ value is large enough, the valley transmittance
becomes stable, represented by steady curves with a uniform local density
between the dots.

Finally, we analyze the role of the staggered potential $\Delta$ on the
valley-transmission output, as shown in Fig. 4. The staggered potential opens
a gap, playing the role of a potential barrier. Consequently, the valley
imbalance is maintained for all values of $E_{F}<\Delta$. However, the valley
polarization is rapidly destroyed when $E_{F}\geq\Delta$. When the confining
potential of the IQD is bigger than the incident energy, the current is
completely blocked, and only resonant states from each valley ${\bf K_{1}}$
$({\bf K_{2}})$ are allowed to transmit current. In this context, the choice
of incident energy must consider the strength of the IQD confinement
potential.

Fig. 4(b) and (d) show that the conductance is non-zero only for particular
resonances where scatterings occur at specific values of the ISOCs and
energies suited by edge-dot coupling [43, 44]. In this context, the
heterostructures of graphene and hBN would exhibit this regime since, in this
system, the confinement potential $\Delta$ is significant [45] and therefore,
its influence on the valley filtering process would have a negative impact.

The results discussed in Fig. 4 are essential since they allow us to establish
a criterion for defining the best materials for islands. Indeed, for a better
response, it is necessary to use a material that provides a weak or zero
confinement potential. If so, the IQDs will operate efficiently, and the
overall system might be used to monitor valley-driven current by either tuning
the ISOC or the RSOC. The question we might ask then is which materials are
better suited for this response. Such a case would be materialized in
heterostructures of twisted graphene and monolayers of transition-metal
dichalcogenides (TMDCs). In such scenarios, the twisting angle substantially
decreases the confinement potential, and the dominant parameters will be the
valley-Zeeman and RSOCs [35]. A concrete example of a realistic proximity
effect is discussed later on.

An alternative description of these effects is by considering the area around
the IQDs with wider staggered potentials as an electron/hole bilayer system
where the electrons are essentially required to overcome an offset barrier to
be scattered through. This action is analogous to the massless-massive
electron-hole system in a transverse electric field in graphene on a TMD
substrate. In this context, the value of the Rashba coupling has a direct
effect since it impacts the offset barrier. Due to its complexity, we postpone
the study of these effects for future work. (For more details, see Refs. [29]
and [28].)

Figure 4: Valley conductance vs intrinsic spin-orbit length $\lambda_{I}$ for
different values of staggered potential $\Delta$. The left (right) panels
correspond to $\Delta<E_{F}$ ( $\Delta>E_{F}$).

####  III.1.3 Valley-Hall and bulk conductivities

We briefly discuss the emergence of Hall and bulk conductivities in these
structures. We visualize the origin of these conductivities by mapping the
local current flow and highlighting the valley polarization with solid blue
(red) curves indicating the ${\bf K_{1}}$ (${\bf K_{2}}$) index in Fig. 5.
Analysis of the figure reveals that the generated valley currents are composed
of bulk-driven and Hall-driven currents. In this process, the local current
where both valleys are scattered shows that each valley is conducting with
either bulk or Hall currents depending on the sign of the Fermi energy. This
analysis suggests that a periodic array of dots with the best choice of SOCs
offers an alternative mechanism for generating valley-neutral Hall currents
since both valleys contribute to the current in the same direction, although
through different regions. A realistic example is discussed in Sec. III.3
where we find induced SOC terms that facilitate valley-Hall current for a
given TMDs island and twist angle.

Figure 5: Local valley polarized currents at ${\bf K_{1}}$ and ${\bf K_{2}}$
valleys, in the presence of only staggered intrinsic SOC,
$\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}=0.075t$. The left (right) panels are for
$E_{F}>0$ ($E_{F}<0$).

###  III.2 Description of the model with weak couplings

A quantitative description of the model with realistic parameters requires
discussing the above results for the weaker strengths of the various SOCs. An
important question is how the valley-conductance is affected by changes in the
RSOC values.

####  III.2.1 Valley-dependent conductance by tuning the Rashba coupling

In this context, we consider weaker ISOCs as expected in realistic settings
and an increased and controllable strength of RSOC. We notice that IQDs refer,
in reality, to any heterostructure that exhibits $C_{3v}$-symmetry. This
section considers only ISO and RSO couplings, while the pseudospin inversion
asymmetry (PIA) coupling in systems with broken inversion symmetry is
discussed in the following section. A typical example of such a setup would be
twisted graphene/transition-metal-dichalcogenides heterostructures since the
relative rotation leads to a negligible value of PIA coupling. To this end,
the IQDs with $C_{3v}$-symmetry might lead to potential applications since the
Rashba coupling can be tuned using a transverse electric field. As realized,
we consider weak staggered couplings
$\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}=0.015t$, we inject electrons with
$E_{F}=\pm 0.035t$ into several-identical QDs, and by shifting the Rashba
coupling with the aid of top gates, we compute the valley conductance through
the system.

In Fig. 6 (a) and (b), we show results for valley transmittance when the top
gate is used to obtain a valley polarized conductance for coupling within the
ranges $0.07t\leq\lambda_{R}\leq 0.12t$ and $0.1t\leq\lambda_{R}\leq 0.2t$,
for the two values of the staggered potential $\Delta=0$ anc $\Delta=0.02t$
correspondingly. The valley transmission for ${\bf K_{1}}$, $({\bf K_{2}})$
jumps almost from 2 to 0 for positive (negative) incident energy, indicating
that only electrons with either valley ${\bf K_{1}}$ or valley ${\bf K_{2}}$
go through the system.

The most important conclusion from the results is that reasonable and better
control of the valley transmission can be obtained by tuning a gate bias using
IQDs with weak couplings. In this context, a quantitative understanding might
easily be reached. By controlling the Rashba strength, the wave function may
interfere destructively or constructively depending on the valley and the sign
of the incident energy. It is seen that at $E_{F}>0$ and $\lambda_{R}\geq
0.05t$, the wave function around ${\bf K_{1}}$ is transmitted through the IQDs
region, where the wave function around ${\bf K_{2}}$ vanishes (is reflected),
resulting on a polarized valley transmission, invertible by changes in the
sign of the incident energy.

These results can be understood from the analysis of the band structure in
Fig. 6 (c) and (d). By analyzing the RSOC’s contributions to the edge state,
we observe that for weak Rashba couplings, such as $\lambda_{R}=0.025t$, the
edge state is mixed with the bulk conduction bands. However, for
$\lambda_{R}=0.075t$, the system might sustain isolated and steady helical
edge states, depending on the spin nature, as discussed in earlier works [29].

Figure 6: Panels (a) and (b): Valley conductance vs. Rashba spin-orbit
coupling $\lambda_{R}$ for $\Delta=0$ and $\Delta=0.02t$, respectively. Panel
(c) and (d) show the energy bands in the staggered case
($\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}$) for $\lambda_{R}=0.025t$ and
$\lambda_{R}=0.075t$, accordingly. Panels (e) and (f) show the conductance vs
the strength of the PIA coupling. Panel (g) and (h) show the energy bands for
$\lambda_{P}=0$ and $\lambda_{P}=0.075t$, accordingly. Insets in (a, b, e, f)
show spin-resolved valley conductance per valley.

In addition to the valley filtering, the interplay between the valley and spin
degrees of freedom might realize new transport regimes. As shown in the inset
of Fig.6 (a) and (b), the wave functions around ${\bf K_{1}}$ or ${\bf K_{2}}$
(depending on the incident energy) are spin-dependent. Indeed, the Rashba
coupling allows off-diagonal spin-flipping terms. Therefore, the spin
conductance is also directly controlled by the Rashba strength, leading to
valley-spin interplay dynamics. For a given incident energy, the spin-
conductance $G_{\uparrow}({\bf K_{1}})$ oscillations have several phase shifts
that correspond to the wavelengths of the oscillations determined by the
Rashba coupling strength. Notice that the overall conductance of spin-
polarized carriers, $G_{\uparrow}({\bf K_{1}})=G_{\uparrow\uparrow}({\bf
K_{1}})+G_{\uparrow\downarrow}({\bf K_{1}})$ remains constant. Our
calculations suggest valley processes dominate the conductance as spin-flips
between neighboring sites require wider IQDs. Details on spin-dependent
conductance in these proximity effect structures have been fully addressed in
Ref. [29].

####  III.2.2  IQDs with lower- or higher-order symmetries

$C_{3v}$ symmetry is a property of hexagonal systems described by point group
symmetry methods. To generalize the applicability of the model presented, it
might be useful to lower (or raise) the symmetry of the IQDs. This strategy
provides more options to describe novel regimes as it allows for removing (or
adding) different spin-orbit terms. (i) In the case $\lambda_{R}=0$, the point
group is increased to $D_{3h}$ where the sublattice inversion asymmetry
defines a point group of a planar system sustaining a triangular lattice with
two staggered (non-equivalent) sublattices [27]. For such a case, the presence
of the ISOC leads to the transmittance spectrum similar to Fig 2 (c) and (d).
(ii) To lower the symmetry of the IQDs, one might include an additional
staggered SOC term called PIA coupling (pseudospin inversion asymmetry for
short). The spin-orbit Hamiltonian, in this case, will have an extra term:

$\displaystyle H_{PIA}$ $\displaystyle=$
$\displaystyle(\sqrt{3}a_{0}/2)\big{[}\lambda_{PIA}^{(A)}(\sigma_{z}+\sigma_{0})$
(7) $\displaystyle+$
$\displaystyle\lambda_{PIA}^{(B)}(\sigma_{z}-\sigma_{0})\big{]}(k_{x}\sigma_{y}+k_{y}\sigma_{x}).$
(8)

where $\lambda_{PIA}^{(A)}$ and $\lambda_{PIA}^{(B)}$ are the staggered PIA
coupling and $a_{0}$ is the lattice constant.

As shown in Fig. 6 (e) and (f), the PIA coupling does not affect the valley-
conductance (even when it is a staggered spin-orbit term), neither for lower
nor higher strength values. Therefore, we can confirm that only ISOC and RSOC
control the valley process, as discussed in Fig. 2 and 6 (a, b). Hence, a QD
with $\lambda_{PIA}\neq 0$ brings forth the same behavior as
$\lambda_{PIA}=0$, and a similar valley response will be observed. Indeed, the
spectrum of the energy bands in Fig. 6 (g) and (h) supports these results,
where we observe that for either neglected ($\lambda_{PIA}=0$) or strong PIA
terms ($\lambda_{PIA}=0.075t$) the lower bulk and edge states are insensitive
to the coupling strength.

This result is exciting as it implies that one might use TMDs to create IQDs
with $C_{3v}$-symmetry but with a $\lambda_{PIA}\neq 0$ without affecting
their valley transport properties. Indeed, we might easily tune the intrinsic
and Rashba parameters in such a system using TMD quantum dots. For instance,
based on density functional theory fittings, both parameters are strongly
sensitive to the electric field, twisting, and/or vertical strain effects [28,
46, 47]. The possibility of controlling the decrease or increase in the
coupling strengths is highly desirable to monitor the valley polarization.
Additionally, the spin-conductance is weakly affected (see inset of Fig. 6 (c,
d)) by the presence of PIA as compared to Rashba terms. For more details about
spin dependence and the effect of PIA coupling on spin polarization, see [48].

####  III.2.3 Valley confinement with IQDs

This subsection addresses the possibility of attaining confinement of valley-
polarized electrons. Adopting the same sample parameters as in subsection
III.1, we see that by zooming in on the transmission spectra from Fig. 2 (c,
d), one can additionally recognize a resonance in the transmittance emerging
in the IQDS, at $E_{F}=0.035t$ for higher values of ISOC
($\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}=0.065t$). The scattering through the
induced region creates electron confinement for small-size IQDs ($r_{0}=7$
nm). The resonance might be advantageously produced for weak ISOC at specific
incident energy values. This case defines a real proximity effect where we
tune $E_{F}$ to confine the valley states as shown in Sec.III.3.

Figure 7: Local density of valley-localized states near ${\bf K_{1}}$ (left)
and ${\bf K_{2}}$ (right) at $E_{F}=0.035t$, obtained for values of
$\lambda_{I}^{(A)}=-\lambda_{I}^{(B)}=0.065t$. The red and blue arrows show
each valley state’s circulation direction. Panels (a, b) and (c, d) are for
single and 2-chain-IQDs.

A detailed analysis of the local current profiles of valley states, as
depicted in Fig. 7, shows a better representation of such a kind of valley-
dependent electron confinement. Indeed, we observe that valley state
confinement can be achieved by using either a single or a chain of IQDs.
Importantly, the incident valley-states are trapped around the IQD region for
appropriate SOC choices. This results in the confinement of the states, a
process that can be regarded as the product of multiple internal reflections
that trap the Dirac fermions by interference processes.

The system with symmetric IQDs appears to provide confinement for states at
both valleys simultaneously in the same spatial region (we do not observe
splittings of the valley confinement states). However, to achieve control over
the valleys and multiple valleytronic and optoelectronic functionalities,
using these IQDs, it will be necessary to create splittings and confining with
unique quasi-bound states (belonging to either ${\bf K_{1}}$ or ${\bf K_{2}}$)
in the local regions. The valley-bound states might be split by considering
asymmetric quantum dots with different shapes and point group symmetries. This
topic will be addressed in future work.

###  III.3 Realistic IQDs in heterostructures of twisted graphene and TMD
monolayers

As discussed in Sec. III.2.2 and III.2.3, a quantitative study of the model
has been addressed in the case of weak SOCs. To provide a comprehensive
picture of concrete setups, we perform a real case study and show how an IQD
might be realized in realistic experimental conditions.

To start, we consider IQDs produced with four different semiconducting TMDs
islands: MoSe2, WSe2, MoS2, and WS2 deposited on graphene monolayer, as shown
in Fig. 1 (b). Furthermore, we analyze the proximity-induced spin-orbit
couplings in the case of twisted TMDs on graphene. Based on the work of T.
Naimer et et al. [35], twisting leads to tuning the magnitudes of the valley-
Zeeman and Rashba spin-orbit couplings. Interestingly, the amplitude of the
staggered potential $\Delta$ in twisted G-TMDs is significantly weak and tends
to zero at some specific twisting angles, thus favoring valley filtering, as
discussed in Fig. 4. Based on the tight-binding parameters derived from first
principles calculations [35], we compute the valley conductance through four
$(n=4)$ IQDs at $E_{F}=0.035t\simeq 0.095$meV and show results in Fig. 8.

Figure 8: Panels (a) and (b) show valley conductance vs twist angle $\theta$
for IQDs made of MoSe2 on graphene at $E_{F}=0.095$meV and $-0.095$meV,
respectively. Panels (c) and (d) correspond to IQDs of WSe2 on graphene for
the same values of $E_{F}$. In panels (e)-(f) and (g)-(h) we show valley-
polarized currents for MoSe2 for positive and negative bias ($E_{F}$),
respectively

.

Valley transmittance with first propagating modes in twisted TMDs: It is
important to state that in all previous sections, we considered the
transmittance spectrum at $E_{F}=\pm 0.035t$; a value that ensures that the
valley-dependent conductance is addressed independently with the simultaneous
excitation of the first propagating modes in ${\bf K_{1}}$ and ${\bf K_{2}}$
since the Fermi energy $E_{F}$ is larger than the valley-mode spacing gap. For
more details, refer to our previous work [49].

In Fig. 8, we inject electrons with $E_{F}=\pm 0.095$meV into four identical
IQDs produced with MoSe2 and WSe2 islands. The induced SOCs are tuned by
twisting the TMDs and allow for monitoring of the valley response in the
system. At positive incident energy, the valley conductance from valley ${\bf
K_{1}},({\bf K_{2}})$ is 2 (0) in the units of $e^{2}/\hbar$. Only electrons
from either valley ${\bf K_{1}}$ or ${\bf K_{2}}$ go through the system (see
panels (a) and (c)). As previously discussed, we also confirm that changing
the Fermi energy sign reverses the valley polarity, as shown in Fig. 8 (b) and
(d). Additionally, as shown in 8 (a) and (b), twisting in the case of MoSe2
(that induces a weak staggered potential [35]) leads to control of valley
transmittances similar to what we found when discussing the model with weak
couplings in Fig. 6.

The obtained results with MoSe2 can be applied to obtain pure valley polarized
conductances with the twist angle outside the limits of
$9^{\circ}-14^{\circ}$. Hence, the valley transmission from ${\bf K_{1}}$
(${\bf K_{2}}$) jumps almost from 2 to 0 for positive (negative) incident
energy, depending on twisting angles and full valley polarized currents. To
emphasize these findings, we plot the local current for $\theta=15^{\circ}$ at
$E_{F}=0.095$meV ($-0.095$meV) for MoSe2 (WSe2). Indeed, one observes that
only one valley is scattered, which leads to conducting Hall currents
depending on the sign of the incident Fermi-energy. These results are
important compared to those shown in Fig. 5, where both valleys conduct
through bulk and Hall currents. Hence, the IQDs based on semiconducting TMDs
MoSe2 and WSe2 might be used as promising islands for generating valley Hall
signals.

Figure 9: Valley conductance vs twist angle $\theta$ for IQDs made of MoSe2,
WSe2, MoS2 and WS2 on graphene, respectively. Panels in the top (bottom) show
conductances at $E_{F}=0.2$meV (0.5)meV.

Valley transmittance beyond first propagating modes in twisted TMDs: To
observe transport involving states beyond the first propagating modes, we set
the incident energy of the injected electrons at $E_{F}=\pm 0.1,\pm 0.3$meV.
In this case, we will consider all TMDs: MoSe2, WSe2, MoS2, and WS2 where the
valley-dependent conductance is shown in Fig. 9. Similarly to what was
discussed, the valley process is sensitive to the magnitudes of the induced
valley-Zeeman and Rashba spin-orbit couplings and the incident Fermi energy.
The TMD islands, as well as the twisting, allow a variety of options (several
values of ISOCs) to monitor valley-driven currents, either in bulk or along
the edges, where the best choice is to set the incident energy at
$E_{F}=0.035t$ (within the first propagating modes). We observe that for some
of the twist angle values, the transmittance of both valleys is zero, which
might be explained either by zero current (OFF) or by the presence of the
confinement states.

Valley confinement in twisted TMDs: Previously, we have shown in Fig. 7 that
at $E_{F}=0.035t$, the resonance or confinement might be produced at higher
values of ISOC strengths. This feature is somehow challenging to attain in
available setups. To confirm the presence of resonances in realistic
conditions, we tune both the incident Fermi energy and the twist angle. In
Fig. 10, we show that the resonance condition might be obtained at weak
induced SOCs at some specific incident energy and twist values. Importantly,
we observe that at $\theta=22.7^{\circ}$ ($27^{\circ}$), the system can
confine valley states in the IQD region of MoSe2 (WSe2) at a lower energy
$E_{F}=\pm 0.035t=\pm 0.098$meV (within the range of the first propagating
mode). Indeed, for appropriate choices of the SOC values related to the chosen
TMD and respective twist angle, the incoming electron is trapped around the
IQD area, as shown in Fig. 10. More precisely, for $E_{F}>0$, the IQDs of
MoSe2 and WSe2 islands (panels Fig. 10 (a), 10 (b), 10 (e), 10 (f)) reveal
higher localization throughout the scattering regions, supported by valley-
localized states from valley ${\bf K_{1}}$. However, the valley-localized
states from ${\bf K_{2}}$ are blocked at the first IQD. The process might be
reversed for $E_{F}<0$ as we have shown in the lower panels of Fig. 10 10 (c),
and 10 (d), for MoSe2 and Fig. 10 10 (g), and 10 (h) for WSe2. This is an
important result since this class of materials allows or induces valley
confinement that might be used to process the optical responses and detect and
monitor valley polarization [14, 15, 16].

By increasing the incident energy beyond the first propagating mode, it is
possible to confine both valley states around the IQDs of WS2 islands as shown
in Fig. 10 (k) and (l). Both confined states are highly localized around all
IQDs in the scattering region, as predicted in [50]. Interestingly, shifting
the sign of the incident energy reverses the path of interference, and both
valleys exchange the propagating direction.

Figure 10: Real-space profiles for both valley currents: red (blue) lines show
the current profile from states originating in valley ${\bf K_{1}}\quad({\bf
K_{2}}$). Panels (a)-(d) show the valley conductance vs twist angle $\theta$
for IQDs made of MoSe2 on graphene at $E_{F}=\pm 100$ meV. Panels (e)-(h) are
for WSe2 on graphene at $E_{F}=\pm 95$ meV. Panels in (k) and (l) correspond
to WS2 at higher energy $E_{F}=-300$ meV.

##  IV Conclusions and summary

We have used a tight-binding approach to obtain an effective, low-energy
Hamiltonian to investigate valley-dependent transport through an array of
proximity-induced quantum dots with $C_{3v}$ symmetry in a zigzag graphene
nanoribbon. The model is a natural extension of the structure studied in
Ref.[51] that analyzed spin and charge distributions in a graphene quantum dot
induced by a WSe2 island. Because the fabrication of graphene/TMD
heterostructures with quantum dots of one or the other material is rapidly
evolving [52, 53, 54], experimental implementations for arrays of quantum dots
will likely be realized in the near future.

Our results show that the valley conductance exhibits an interesting behavior
sensitive to the model parameters and the values of the various SOCs induced
by the islands that decorate the ribbon. At the same time, this sensitivity is
a great advantage for monitoring and controlling the different valley filter
regimes observed. The results reveal the localization of valley-centered
states produced by the competition between Zeeman and Rashba couplings in
narrow IQDs. For larger IQDs, the valley states combine with bulk states, and
the valley polarization is considerably deteriorated [34].

Following Ref. [30], the qualitative description of the model shows that by
varying Zeeman couplings, the conductance through a symmetric chain of quantum
dots displays square-shaped curves with wide gaps. However, these features
tend to vanish for some large ISOC values with the subsequent vanishing of
valley polarization. Furthermore, for specific ranges of ISOC values, both
valleys are present, implying that the device does not display a perfect
valley-transistor behavior. An important conclusion from the qualitative
description is the similarity of this system to the Datta-Das transistor. In
both cases, the spin conductance is directly controlled by the strength of the
Rashba SOC. In this respect, the most important SOC for good valley
polarization is the staggered intrinsic SOC. These features could be helpful
for further exploration of actual devices.

The sensitivity to the Rashba coupling is discussed in the presence of weak
ISOC to describe realistic settings where proximity effects develop. The
valley conductance shows that tuning the valley polarization and switching the
valley scattering in the system using a top gate is possible. The possibility
of using TMDs as decorating islands to form IQDs has also been discussed. We
have shown that the presence of the PIA coupling, characteristic of these
structures, does not affect the valley-polarization. Therefore, such islands
can be used to tune the polarization by either strain or twists since the
Rashba and Zeeman coupling are sensitive to external electric or strain field
effects [28, 46, 47].

Finally, we have applied these models to solve examples of realistic material
combinations. To give a comprehensive picture, we have used graphene/TMD
heterostructures with different semiconducting materials: MoS2, MoSe2, WSe2,
and WS2. We noticed that IQDs based on semiconducting TMDs might be used as
promising islands for generating valley Hall signals. Indeed, TMDs allow
valley filtering processes and break the valley degeneracy, producing a
valley-polarized current that favors valley selection by tuning the sign of
the incident Fermi energy or the value of twist angles. Notably, the TMD
island and control of the twisting angle allow various options (as they
determine different values of ISOCs) to monitor valley-driven currents,
confine both valley states simultaneously in the same region, or split the
valley confinement states.

Achieving a nearly square-wave transmission and a valley-valve effect for the
${\bf K_{1}}$ or ${\bf K_{2}}$ valleys is highly desirable for device
applications. Moreover, the confinement of quasi-bound states from either
valley is extremely important for manipulating optoelectronic interactions
[13, 23, 55] and valley-qubit systems [17, 3, 56]. Furthermore, the mechanism
for generating valley-Hall conductivity with valley-neutral currents could be
handy to obtain pure valley-Hall signals. The system proposed in this
manuscript exhibits features close to these goals; however, asymmetric QDS
with controlled shape and several induced point group symmetries could lead to
richer results. We plan to address this issue in future work.

###### Acknowledgements.

The authors acknowledge computing time on the SHAHHEN supercomputers at KAUST
University, Saudi Arabia, and the supercomputers at the Centre for Research in
Molecular Modeling (CERMM), Richard J. Renaud Science Complex, Concordia
University. AB would like to thank Dr. Adel Abbout at KFUPM, Saudi Arabia, for
helpful discussions.

##  Appendix A Valley-dependent conductance

As stated, we use the scattering matrix formalism to calculate the individual
valley conductances. We separate the propagating modes in the leads depending
on their velocity and momentum direction using the Kwant package [40]. In the
scattering region, valley currents might be mixed due to inter-valley
scattering. However, as the two valleys are far apart in the Brillouin zone,
effective valley mixing will require short-range potentials. A graphene
membrane wraps smoothly for graphene deposited on top of islands to minimize
strain effects. Thus, the sharp atomic termination of the island (that could
give rise to inter-valley mixing) is effectively masked. These issues will
have a negligible contribution for islands on top of graphene if the islands
and their separation are large enough [57].

To this end, we consider only propagating states for which $\Phi({\bf vF}>0)$.
The applied source-drain voltage chooses the direction of the current composed
of these states, which have both spin and valley degrees of freedom. We focus
first on the valley degree and lift the valley degeneracy by defining the
propagating wave functions $\Phi_{\bf K_{1}}=\Phi({\bf k}>0),\Phi_{\bf
K_{2}}=\Phi({\bf k}<0)$.

These wave functions independently solve the scattering problem in the
reciprocal space. To implement the Green’s function formalism [39] we need to
include the scattering matrices $S^{mn}=S_{\bf K_{1}}^{mn}+S_{\bf K_{2}}^{mn}$
with $S_{\bf K_{1,2}}^{mn}$ given by

$S_{\bf K_{1,2}}^{m,n}=\text{Tr}[G_{\bf K_{1,2}}\Gamma^{m}\,G^{\dagger}_{\bf
K_{1,2}}\,\Gamma^{n}],\qquad(m,n=L,R\ or\ R,L);$ (A.1)

The Green’s function and $\Gamma$ matrices are given by

$\displaystyle G(\epsilon,{\bf K_{1,2}})$ $\displaystyle=$
$\displaystyle\left[\left(\epsilon+i\eta\right)I-H_{QD}({\bf
K_{1,2}})-\Sigma\right]^{-1}$ (A.2) $\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle i(\Sigma-\Sigma^{\dagger}).$ (A.3)

$\Gamma$ defines the self-energy of the contacts placed to the left and right
of the scattering region, and the relevant Hamiltonian is $H_{QD}$, cf. Eq.
(1). Then, for each valley mode, the valley conductance at the Dirac cones is
given by Eq. (2).

##  Appendix B Valley-resolved current

To obtain the local density of states (LDOS) and currents per valley, we
obtain the wave functions of the propagating modes $\Phi$ for a given energy
$E$ and site $i$. The propagating wave functions are stored per site depending
on their momentum $\left\\{\Phi({\bf K_{1}}),\ \Phi({\bf K_{2}})\right\\}$ and
their spin degree of freedom. The resulting LDOS, at a given site $i$ in the
sample, is defined by

$\text{LDOS}^{\bf K_{1,2}}\left(E,i\right)=\sum_{l}\left|\langle i|\Phi^{\bf
K_{1,2}}_{l}\rangle\right|^{2}\delta(E-E_{l})$ (B.1)

where the sum is over all electron eigenstates
$|\Phi_{l}\rangle=c_{l}^{\dagger}|0\rangle$ of the Hamiltonian $H_{QD}$ in Eq.
(1) with energy $E_{l}$. The valley-resolved LDOS in Eq. (B.1) is calculated
using Chebyshev polynomials [58] and damping kernels [59].

The corresponding density operator and the continuity equation are expressed
as

$\rho_{q}^{\bf K_{1,2}}=\sum_{a}[\Phi_{a}^{\bf
K_{1,2}}]^{*}\,H^{s}_{q}\,\Phi_{a}^{\bf
K_{1,2}},\qquad\frac{\partial\rho_{a}^{\bf K_{1,2}}}{\partial
t}-\sum_{b}J_{a,b}^{\bf K_{1,2}}=0.$ (B.2)

where $a$ refers to all sites in the scattering region, $J_{ab}^{\bf K_{1,2}}$
is the valley-resolved current, and $H^{s}_{q}$ is the $q$-Fourier component
of the hopping term in the scattering region.

For a given site density $\rho_{a}$, we sum over all the neighboring sites
$b$. As a result, the valley current $J^{{\bf K_{1,2}}}_{ab}$ takes the form

$J_{a,b}^{\bf K_{1,2}}=[\Phi^{\bf K_{1,2}}({\bf
v>0})]^{*}\left(i\sum_{\gamma}H^{*}_{ab\gamma}H^{s}_{a\gamma}-H^{s}_{a\gamma}H_{ab\gamma}\right)[\Phi^{\bf
K_{1,2}}({\bf v>0})]$ (B.3)

where $H_{ab\gamma}$ is a component of a rank-4 tensor that can be represented
as a vector of matrices, and $\gamma$ is an index that runs over sites in real
space. In this expression, Latin indices go over sites, and Greek indices run
over the degrees of freedom in the Hilbert space. For more details about
obtaining the current operator, see references [49, 29, 48, 40].

References

## References

  * [1] N. Rohling and G. Burkard, New Journal of Physics 14, 083008 (2012).
  * [2] J. R. Schaibley, H. Yu, G. Clark, P. Rivera, J. S. Ross, K. L. Seyler, W. Yao, and X. Xu, Nature Reviews Materials 1, 16055 (2016).
  * [3] G. Szechenyi, L. Chirolli, and A. Palyi, 2D Materials 5, 035004 (2018).
  * [4] D. Culcer, A. L. Saraiva, B. Koiller, X. Hu, and S. Das Sarma, Phys. Rev. Lett. 108, 126804 (2012).
  * [5] E. A. Laird, F. Pei, and L. P. Kouwenhoven, Nature nanotechnology 8, 565 (2013).
  * [6] N. Hadadi, A. Belayadi, A. AlRabiah, O. Ly, C. A. Akosa, M. Vogl, H. Bahlouli, A. Manchon, and A. Abbout, Phys. Rev. B 108, 195418 (2023).
  * [7] L. A. Benitez , J. F. Sierra , W. S. Torres , A. Arrighi, F. Bonell, M. V. Costache and S. O. Valenzuela. Nature Physics 14, 303 (2018).
  * [8] Y. K. Luo, J. Xu, T. Zhu, G Wu, E. J. McCormick, W. Zhan, M. R. Neupane, R. K. Kawakami, Nano Lett 6, 3877 (2017).
  * [9] V. Bonbien, F. Zhuo, A. Salimath, O. Ly, A. Abbout and A. Manchon, J. Phys. D. Appl. Phys. 55, 103002 (2022).
  * [10] D. Xiao, G. B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. 108, 196802 (2012).
  * [11] J. Qi, X. Li, Q. Niu, and J. Feng, Phys. Rev. B 92, 121403(R) (2015).
  * [12] A. Cresti, B. K. Nikolic, J. H. García, and S. Roche, Riv. del Nuovo Cim. 39, 587 (2016).
  * [13] Z. Zhang, X. Ni, H. Huang, L. Hu, and F. Liu, Phys. Rev. B 99, 115441 (2019).
  * [14] A.Krasnok, and A. Alù, Appl. Sci. 7, 1157 (2018).
  * [15] S. Sharma, P. Elliott, and S. Shallcross, Optica 9, 947 (2022).
  * [16] R. H. Godiksen, S. Wang, T. V. Raziman, J. G. Rivas and A. G. Curto, Nanoscale 14, 17761 (2022).
  * [17] Y. Wu, Q. Tong, G-B. Liu, H. Yu and W. Yao, Phys. Rev. B 93, 045313 (2016).
  * [18] P. A. Maksym and H. Aoki, Phys. Rev. B 88, 081406(R) (2013).
  * [19] A. Orlof, A. A. Shylau and I. Zozoulenko, Phys. Rev. B 92, 075431 (2015).
  * [20] W. A. Behn, Z. J. Krebs, K. J. Smith, K. Watanabe, T. Taniguchi, and V. W. Brar, Nano Lett 21, 5013 (2021).
  * [21] SY. Li and L. He, Front. Phys. 17, 33201 (2022).
  * [22] T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, Q. Niu, P. Tan, E. Wang, B. Liu, Nature communications 3, 887 (2012).
  * [23] T. Aktor, J. H. Garcia, S. Roche, A-P. Jauho, and S. R. Power, Phys. Rev. B 103, 115406 ( 2021).
  * [24] Y. Iwasa, in JSAP-OSA Joint Symposia 2015 Abstracts, (Optica Publishing Group, 2015), p. $15p_{2}D_{1}$.
  * [25] Md.T Mahmud and N. Sandler, Phys. Rev. B 102, 235410 (2020).
  * [26] F. Solomon and S. R. Power, Phys. Rev. B 103, 235435 (2021).
  * [27] D. Kochan, S. Irmer, and J. Fabian, Phys. Rev. B 95, 165415(2017).
  * [28] M. Gmitra, D. Kochan, P. Hagl, and J. Fabian, Phys. Rev. B 93, 155104 (2016).
  * [29] A. Belayadi and P. Vasilopoulos, Nanotechnology 34, 085704 (2023).
  * [30] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * [31] M. Gmitra and J. Fabian, Phys. Rev. B 92, 155403 (2015).
  * [32] Z. Wang, D.-K. Ki, H. Chen, H. Berger, A. H. MacDonald, and A. F. Morpurgo, Nat. Commun. 6, 8339 (2015).
  * [33] Aron W. Cummings, Jose H. Garcia, Jaroslav Fabian, and Stephan Roche Phys. Rev. Lett. 119, 206601 (2017).
  * [34] T. Frank, P. H$\ddot{o}$gl, M. Gmitra, D. Kochan, and J. Fabian Phys. Rev. Lett. 120, 156402 (2018).
  * [35] T. Naimer, K. Zollner, M. Gmitra, and L. Fabian, Phys. Rev. B 104, 195156 (2021).
  * [36] P. Hogl, T. Frank, K. Zollner, D. Kochan, M. Gmitra, and J. Fabian, Phys. Rev. Lett. 124, 136403 (2020)
  * [37] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * [38] T. S. Ghiasi, J. Ingla-Ays, A. A. Kaverzin, and B. J. van Wees, Nano Lett. 17, 7528 32 (2017).
  * [39] S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, 1995).
  * [40] C. W. Groth, M. Wimmer, A. R. Akhmerov, X. Waintal, New J. Phys. 16, 063065 (2014).
  * [41] J. H. Garcia and T. G. Rappoport, 2D Materials, 3, 024007 (2016).
  * [42] T. P. Cysne, J. H. Garcia, A. R. Rocha, and T. G. Rappoport, Phys. Rev. B 97, 085413 (2018).
  * [43] N. Myoung, J.-W. Ryu, H.C. Park, S.J. Lee, S. Woo, Phys. Rev. B 100, 045427 (2019).
  * [44] H. C. Park, M. Son, S. J. Lee, N. Myoung, Journal of the Korean Physical Society 78, 1208-1214 (2021).
  * [45] K. Zollner, M. Gmitra, and J. Fabian, Phys. Rev. B 99, 125151 (2019).
  * [46] Li Y and Koshino M Phys. Rev. B 99, 075438 (2019).
  * [47] David A, Rakyta P, Kormanyos A, and Burkard G Phys. Rev. B 100, 085412 (2019).
  * [48] A. Belayadi and P. Vasilopoulos, Nanotechnology 34, 365706 (2023).
  * [49] A. Belayadi, N. A. Hadadi, P. Vasilopoulos, and A. Abbout Phys. Rev. B 108, 085419 (2023).
  * [50] E. Sie, J. McIver, Y. H. Lee, F. Liang, K. Jing, and G. Nuh, Nature Mater 14, 290-294 (2015).
  * [51] H-Y. Ren, Y-N. Ren, Q. Zheng, J-Q. He and L. He, Phys. Rev. B 108, L081408 (2023).
  * [52] A. Manikandan, Y-Z. Chen, Ch-Ch. Shen, Ch-W. Sher, H-Ch. Kuo and Y-L. Chueh, Progress in Quantum Electronics 68, 100226 (2019).
  * [53] K. J. Singh, T. Ahmed, P. Gautam, A. S. Sadhu, D-H. Lien, S-Ch. Chen, Y-L. Chueh and H-Ch. Kuo, Nanomaterials 11, 1549 (2021).
  * [54] A. Alfieri, S. B. Anantharaman, H. Zhang and D. Jariwala, Advanced Materials 35, 2109621 (2022).
  * [55] S. A. Vitale, D. Nezich, J. O. Varghese, P. Kim, N. Gedik, P. Jarillo-Herrero, and M. Rothschild, Small 14, 1801483 (2018).
  * [56] J. Pawlowski, New J. Phys. 21 123029 (2019).
  * [57] J. G. Nedell, J. Spector, A. Abbout, M. Vogl, and G. A. Fiete, Phys. Rev. B 105, 075425 (2022).
  * [58] A. Ferreira and E. Mucciolo, Phys. Rev. Lett 115, 106601 (2015).
  * [59] A. Weibe, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys 78, 275 (2006).