Quantized transport, strain-induced perfectly conducting modes and valley filtering on shape-optimized graphene Corbino devices

TL;DR

This paper demonstrates how shape-optimized graphene Corbino devices can utilize strain-induced pseudomagnetic fields for directional electronic transmission and valley filtering, advancing valleytronics applications.

Contribution

It introduces a framework for designing graphene devices with optimal shapes to generate specific pseudomagnetic field profiles for electronic control.

Findings

Strain-engineered graphene can produce tunable pseudomagnetic fields.

Shape optimization enables directional control of electronic transmission.

Devices can filter charge carriers by valley and current direction.

Abstract

The extreme mechanical resilience of graphene and the peculiar coupling it hosts between lattice and electronic degrees of freedom have spawned a strong impetus towards strain-engineered graphene where, on the one hand, strain augments the richness of its phenomenology and makes possible new concepts for electronic devices and, on the other hand, new and extreme physics might take place. Here, we demonstrate that the shape of substrates supporting graphene sheets can be optimized for approachable experiments where strain-induced pseudomagnetic fields (PMF) can be tailored by pressure for directionally selective electronic transmission and pinching-off of current flow down to the quantum channel limit. The Corbino-type layout explored here furthermore allows filtering of charge carriers according to valley and current direction, which can be used to inject or collect valley-polarized currents, thus realizing one of the basic elements required for valleytronics. Our results are based on a framework developed to realistically determine the combination of strain, external parameters, and geometry optimally compatible with the target spatial profile of a desired physical property --- the PMF in this case. Characteristic conductance profiles are analyzed through quantum transport calculations on large graphene devices having the optimal shape.

Figures

Equations

$$\text{Minimize}\mathcal{I}=$$

$$+$$

$$R(\theta )=\sum _k({\alpha }_k\cos 3k\theta +{\beta }_k\sin 3k\theta )$$

$$R_{\text{out}}(\theta )=\frac{{\alpha }_0}{2}+\sum _{k=1}^{N/2}{\alpha }_k\cos 3k\theta +\sum _{k=1}^{N/2-1}{\beta }_k\sin 3k\theta ,$$

$$L_{\text{min}}^{\text{in}}\le R_i^{\text{in}}\le L_{\text{max}}^{\text{in}},L_{\text{min}}^{\text{out}}\le R_i^{\text{out}}\le L_{\text{max}}^{\text{out}}.$$

$$L_{\text{in}}=(L_{\text{min}}^{\text{in}}+L_{\text{max}}^{\text{in}})/2,\text{ and }{L}_{\text{out}}=(L_{\text{min}}^{\text{out}}+L_{\text{max}}^{\text{out}})/2.$$

$$\{(x,y):\log L_{\text{in}}\le x\le 0,-\pi /3\le y\le \pi /3\}.$$

$$R_{\text{in}}(\theta )+(R_{\text{out}}(\theta )-R_{\text{in}}(\theta ))\left(\frac{R-L_{\text{in}}}{L_{\text{out}}-L_{\text{in}}}\right),$$

$${\mathcal{I}}^{\text{reg}}={\int }_{-\pi /3}^{\pi /3}(R_{\text{out}}^{\prime \prime }(\theta {)}^2+R_{\text{in}}^{\prime \prime }(\theta {)}^2)\mathrm{d}\theta .$$

$$R_{\text{out}}(\theta )$$

$$R_{\text{in}}(\theta )$$

$$\bar{p}\propto \frac{p}{B_0^{3/2}{L}^{1/2}},\kappa \propto \frac{1}{L^3{B}_0}.$$

$$t_{ij}=$$

$$\displaystyle+\bigl{[}V_{pp\sigma}(d_{ij})-V_{pp\pi}(d_{ij})\bigr{]}\frac{(\widehat{n}_{i}\cdot\vec{d}_{ij})(\widehat{n}_{j}\cdot\vec{d}_{ij})}{d^{2}_{ij}},$$

$${G_{\text{p}}}({W_{\text{eff}}})=6\mathcal{T}\,\Bigl{(}\frac{2e^{2}}{h}\Bigr{)},\qquad\mathcal{T}=\frac{T}{T+(1-T)N}.$$

Full text

\altaffiliation

These authors contributed equally to this work.

\altaffiliationThese authors contributed equally to this work.

\alsoaffiliationCentre for Advanced 2D Materials, National University of Singapore, 6 Science Drive 2, Singapore 117546

\phone+65 6601 3642

\alsoaffiliationCentre for Advanced 2D Materials, National University of Singapore, 6 Science Drive 2, Singapore 117546

1 Introduction

An extremely active and promising area of research in strain-engineered two-dimensional electronic systems is the realization of graphene-based nanostructures with tailored pseudomagnetic fields (PMF). The PMF is a fruitful concept in strained graphene arising from the fact that, in the vicinity of the Dirac point, the description of the coupling between lattice deformations and electrons can be captured by two types of strain-induced fields: a displacement potential that couples to electrons as an electrostatic potential, and a strain-induced gauge field that couples similarly to a magnetic field 4. Since the former are expected to be effectively screened by the charge carriers, electrons are then expected to respond to strain mostly through this pseudomagnetic effect. Since (pseudo)magnetic fields are very effective means of guiding the motion and confining charged carriers, this naturally led to proposals and studies of many analogues of magnetic devices using graphene without actual magnetic fields, such as pseudomagnetic barriers 6, 16, 17, 18 or pseudomagnetic quantum dots 8, 19.

There are significant advantages of using pseudomagnetic fields in graphene electronic devices at the nanoscale. The first is that electrons in graphene are not easily confined by electric fields due to the Klein tunneling effect 20 arising from the relativistic-like behavior of electrons in this system. This has been a perennial difficulty and disadvantage of graphene-based electronics, while PMFs can efficiently localize electrons or act as tunneling barriers 6. The second advantage is that PMFs are only felt by the electrons in graphene, and they can be confined to regions as narrow as 5 – 10 nm 21, 12. This is unlike real magnetic fields which are felt not only by the target electronic system but also by whole the environment nearby. In addition, fabricating real magnetic barriers of nanoscale size is practically impossible and their magnetic field would also be fixed in principle, not reconfigurable or tuneable. Finally, the magnitude of PMFs that can be achieved in graphene is remarkably high and in excess of 300 T 21, 13, 12 which means that values of the order of tens of Tesla are easily obtainable, and strong enough for many electronic applications.

There is a key challenge that has persisted since the early proposals for strain-engineered graphene: the ability to tailor the spatial distribution of the PMF in graphene hinges on the ability to define specific non-uniform strain profiles, which is a difficult experimental prospect, especially at the micro and nanoscale. Moreover, whereas it is a straightforward theoretical task to start from a given strain field and determine the associated PMF and its impact in the electronic and transport properties, it is much less trivial and non-unique to request the opposite. Yet, it is precisely the opposite (i.e., specifying the PMF patterns in real space and obtaining the necessary strain profiles) that is of most direct interest for application in electronic devices because, in nearly any such case, one should be able to specify how, where, and to which extent we need electrons to be confined, which can be easily done reasoning in terms of PMF patterns alone, but is generally hard to anticipate in terms of strain alone. The central theoretical question for actual realizations and applications of strain-engineered graphene then becomes: how should one design the experimentally accessible parameters of the system (shape, external forces, constraints, substrate features, contacts, etc.) so that any given desired PMF pattern can be generated under realistic experimental circumstances?

In this work we pose and address one such question. We demonstrate that the geometry of a substrate can be optimized so that a graphene Corbino-type device can have its conductance suppressed when hydrostatically pressurized. This, on the one hand, can lead to extremely high on/off ratios due to the two contrasting transport regimes in the pressurized and relaxed states. Pressurized devices are shown to leak current only through perfect, spatially separated 1D modes arising from snake-type states that are stabilized at the inversion regions of the PMF. This creates a very strong current pinch-off by squeezing (literally in this case) the transport channel to the quantum wire limit and, in addition, establishes the possibility of having strictly quantized conductance on demand without electrostatic confinement. Moreover, since the chirality of PMF-induced snake states is uniquely tied to the valley degree of freedom, such devices can function as effective valley filters or sources of valley polarized currents. Our calculations are performed in an optimization framework that tackles the “reverse” strain-engineering problem in graphene. We begin with an overview of its main ingredients and a discussion of the optimal shapes for nearly uniform PMF throughout the Corbino ring. Subsequently, we describe how the deformation fields so obtained are mapped to an actual graphene lattice whose electronic properties are described in a non-uniform tight-binding approximation. This constitutes the basis for our study of the quantum transport characteristics and local electronic structure in devices of different dimensions, with and without disorder, whose results establish the behavior and phenomenology indicated above.

2 Corbino shape optimization

Stating the problem as an optimization question, we wish to find the shape of a non-circular annular cavity in a substrate such that the graphene flake, when deposited on it from above and subjected to hydrostatic pressure of magnitude $p$, exhibits a pseudomagnetic field (PMF) with a magnitude $B_{0}$ that is constant in as much area of the system as possible, compatible with realistic elastic constraints. The substrate shape and dimensions so determined would constitute the basic information for the design of the corresponding transport experiment. If we assume that the graphene flake is rigidly attached to the substrate, finding the shape of the cavity is equivalent to finding the shape of the graphene flake with clamped boundary conditions at its edges, with the same optimized PMF under the applied hydrostatic pressure. This problem is solved by appealing to a Partial Differential Equation (PDE) constrained optimization technique, discretized by finite elements 22. It involves varying the control variables (the shape of the cavity) and the state variables (the plate deformation) so as to minimize a measure of the distance of the underlying PMF away from a certain desired PMF pattern, subject to the constraint that the state variables must satisfy the equations governing the elastic deformation of the graphene sheet for a given shape.

2.1 Overview of the inverse optimization problem

Our control variable is the shape of the domain occupied by the graphene flake (it is this that we numerically vary in order to approach the desired PMF in the pressurized field). We thus parametrize the domain using a set of control variables: $\Omega=\Omega[c_{1},\ldots,c_{n}]$, to be specified later. The clamped boundary conditions at the two perimeters (and beyond) where the graphene sheet contacts the substrate imply that the material domain of interest is the annular region defined by the patterning of the substrate, and not that part of the graphene that adheres to the substrate itself. The central mathematical problem is as follows:

[TABLE]

subject to six constraint equations expressing the mechanical equilibrium of the graphene sheet, where the PMF $B$ is an explicit and known function of the strain field in the 2D sheet, and $B_{0}$ is the target PMF which, in general can be any predefined function of space. This formulation ensures that the $B$ sought in the optimization satisfies $B^{2}\approx B_{0}^{2}$ as closely as is feasible under the restrictions caused by the geometry and mechanics of the problem (the scheme seeks PMFs that are as close as possible to $\pm B_{0}$; in the Supporting Information we discuss the case where the objective function is chosen to minimize $B-B_{0}$, rather than $B^{2}-B^{2}_{0}$). $\Omega$ represents the domain occupied by the unadhered part of the graphene flake, defined by a characteristic lengthscale $L$. The explicit constraint equations that accompany eq. (1) are provided in the Supporting Information in order not to obscure the main discussion here. The final term is a regularization term necessary for well-posedness of the optimization problem. In reality, the optimization problem will be solved in dimensionless form but we present the dimensional version here for clarity.

2.2 Discretization and numerical solution

Each of the state variables (the three displacement components $v_{1}$, $v_{2}$, $w$ of the plate, and the three components of the bending moment $M_{11}$, $M_{12}$, $M_{22}$) is discretized using linear finite elements. This requires the generation of a triangular mesh covering the entire domain occupied by the graphene flake. However, the shape of this domain $\Omega$ must be parametrized by a finite set of variables in order for the optimization algorithm (which works in a discrete setting) to be appropriately formulated. The simplest way to parametrize the shape of the domain $\Omega$ is to prescribe equidistant control points around the outer and inner boundaries, and interpolate them with a Fourier series in the angular coordinate $\theta$. Furthermore, due to the symmetry of the graphene lattice, the PMF is invariant if the coordinate axes are rotated by an angle of $2\pi/3$. Therefore we may restrict our domain to the region $-\pi/3\leq\theta\leq\pi/3$, and admit outlines of the form

[TABLE]

only. An example of the domain construction from control points is shown in Figure 1. Let $N$ be the (even) number of control points on the outer boundary, with the same number on the inner boundary. Then if $(R_{i}^{\textrm{out}})_{i=1,\ldots,N}$ are the radial coordinates of the outer boundary for the angular coordinates $\theta_{i}=2\pi(i-N/2)/(3N)$, we have

[TABLE]

where the coefficients are found by solving $R_{\textrm{out}}(\theta_{i})=R_{i}^{\textrm{out}}$ for $i=1,\ldots,N$. The same procedure is repeated for the inner boundary control variables $R_{i}^{\textrm{in}}$. The $2N$ control variables are allowed to vary within a user-provided interval,

[TABLE]

For each set of control variables chosen during the iterative optimization process, we will have a domain of a different shape. It would be improper to triangulate each new domain separately: not only would it be prohibitively time-consuming, but most optimization routines work best when the constraints and objective functions vary smoothly when the state and control variables are altered. Therefore, we instead choose a mesh for an annular sector and deform the mesh as the control variables accordingly alter the domain. The annular sector is defined by $L_{\textrm{in}}\leq R\leq L_{\textrm{out}}$ and $-\pi/3\leq\theta\leq\pi/3$, where

[TABLE]

To mesh this annulus in a way which ensures a good resolution at the inner boundary, we construct a uniform mesh of $N_{p}$ points for the rectangle222A mesh size giving around 8 boundary elements per control point was used for the calculations in this article. Meshes were constructed using the Matlab-based mesh routine DistMesh 23.

[TABLE]

We conformally map the mesh for the rectangular domain to the annular sector using the mapping ${X}+\mathrm{i}{Y}=\exp(x+\mathrm{i}y)$. Then, to deform this mesh to the domain defined by the control variables, we simply kept the $\theta$-component of the meshpoint fixed, and affinely displaced the radial component $R=\sqrt{{X}^{2}+{Y}^{2}}$ to

[TABLE]

where $R_{\textrm{out}}(\theta)$ and its equivalent at the inner boundary are calculated from the control points by (3). See Figures 1–1 for a visualization of this procedure. The purpose of the regularization term $\mathcal{I}^{\textrm{reg}}$ in the objective function (1) is to make the optimization problem well-posed by penalizing solutions that have unfavorable properties. In this study we set the regularization function to be

[TABLE]

in order to penalize highly-convoluted boundaries.

As discussed previously, the state variables are ${v}_{1}$, ${v}_{2}$, ${w}$, ${M}_{11}$, ${M}_{12}$, ${M}_{22}$, and are discretized such that the equations are written in terms of their values at each node in the triangulation. Additionally there are boundary conditions to be applied: the displacements are set to zero on nodes comprising the inner and outer boundaries. At the boundaries $\theta=\pm\pi/3$, we impose rotated periodicity conditions on these quantities (expressed in polar coordinates $R$, $\theta$), i.e. $w|_{\theta=\pi/3}=w|_{\theta=-\pi/3}$, and similarly for $v_{R}$, $v_{\theta}$, $M_{RR}$, $M_{R\theta}$, $M_{\theta\theta}$. The discretized PDE system becomes a discrete optimization system: to minimize an objective function $F(\bm{C},\bm{U})$ subject to $6N_{p}$ constraints $\bm{G}(\bm{C},\bm{U})=\bm{0}$, where $\bm{C}$ is the vector of control variables and $\bm{U}$ is the vector of $6N_{p}$ state variables.

3 Optimized device shapes

As the electronic and mechanical properties of graphene remain fixed, there are only three parameters that we can change in this system: the limits on the annulus shape (and hence the lengthscale $L$), the applied (constant) pressure $p$, and the target PMF $B_{0}$. In our calculations we allowed the $N=4$ outer control points to vary in a range $261\,\textup{\AA}\pm 60\,\textup{\AA}$ and the inner points in a range $61\,\textup{\AA}\pm 15\,\textup{\AA}$. We initially set a target PMF of $B_{0}=10\,\textup{T}$ and the pressure was fixed at $100\,\textup{bar}=10^{7}\,\textup{Pa}$. Writing the coefficients to the nearest $0.1\,\textup{\AA}$, the shape of the optimal outer and inner boundaries were calculated to be (in Angstroms)

[TABLE]

This outline corresponds to a striking flower-like geometry, and its deformation under hydrostatic pressure is displayed in Figure 2. To this deformation corresponds the spatial distribution of PMF shown in Figure 2. Though the target PMF was $10\,\textup{T}$, the root-mean-squared value of the calculated PMF was only $4.05\,\textup{T}$. It should be kept in mind that the target PMF is not in general attained in these optimization calculations, but is rather a mechanism for forcing the solution towards our desired direction.

Since solutions to the optimization problem in eq. (1) are obtained for a specific choice of input parameters $p$ and $B_{0}$, one immediate and important question is to assess how robust is the solution for the shape with respect to changes in those parameters. As the optimization and its constraints were solved in a dimensionless setting, the behavior of the solution is dependent on the values of dimensionless parameters rather than true physical quantities. The nondimensionalized system has only two parameters: the dimensionless bending stiffness $\kappa$, and the dimensionless pressure $\bar{p}$. In terms of the three physical quantities $p$, $L$, and $B_{0}$, these parameters behave as (details in the Supporting Information)

[TABLE]

Varying $\kappa$ has little effect on the shape obtained through the optimization procedure (see Fig. S1 and related discussion in the Supporting Information), and so only changes in $\bar{p}$ are relevant in understanding how the optimized shape depends on the choices of input parameters. For instance, increasing the hydrostatic pressure $p$ by a factor of $8$ and the target PMF $B_{0}$ by a factor of $4$ has no effect on the optimized shape, as the dimensionless constant $\bar{p}$ is kept unchanged and the effect of a changed $\kappa$ is negligible. In contrast, if $\kappa$ is kept constant and $\bar{p}$ is allowed to vary, we see a significant variation in the calculated shape as shown in Figures 3a–d. This variation in shape with changes in $\bar{p}$ deserves further explanation. A change in $\bar{p}$ corresponds to varying $p$ while keeping the lengthscale and target PMF constant. But a better way of seeing the effect is to think about an increase in $\bar{p}$ as a decrease in target PMF $B_{0}$, while keeping $p$ and $L$ constant [cf. eq. (11)]. Though this also changes the value of $\kappa$, we have noted that the effect of this is minor. To explore this line of thinking, in the lower row of Figure 3, we report a further set of calculations which kept $p=100\,\textup{bar}$ and $L=261\,\textup{\AA}$ fixed while varying the target PMF $B_{0}$ through $50\,\mathrm{T}$, $20\,\mathrm{T}$, $10\,\mathrm{T}$, and $4\,\mathrm{T}$. For high target PMFs, the optimization forces the inner boundary to be highly convoluted in order to obtain the high strain gradients needed to achieve the target PMF. For low target PMFs, an annular solution is already close to the optimal solution, and the optimization only needs to make small adjustments to the outline to make the magnitude of the field more uniform. The results for $10\,\textup{T}$ and $20\,\textup{T}$ can be seen as intermediates between these two extremes. The two extremes also explain the results of Figure 3a–d calculated using different pressures for a given $B_{0}$: for a low pressure, a convoluted interior boundary is required to hit the target PMF, while for high pressure an annular solution is already close to the target and only minor modifications of shape are needed. In the end, the two rows of Figure 3 illustrate what is implied by eq. (11): the changes in the optimal shape with increasing $p$ (and fixed $B_{0}$) follow the same trend as the changes with decreasing $B_{0}$ (and fixed $p$). Of the calculations shown in Figure 3e–h we selected the solution (g), which was shown before in Figure 2, with $B_{0}=10\,\textup{T}$ and $\bar{p}=0.337$ as giving a good middle ground between PMF smoothness and strength.

If we are now given an arbitrary lengthscale $L$ and pressure $p$, in order to find the optimal shape we need to choose an appropriate value of $B_{0}$. Since we have seen that the result for $\bar{p}=0.337$ gives an acceptable outcome, we could choose a $B_{0}$ that results in this value of $\bar{p}$. This gives us a system of equations with the same value of $\bar{p}$ but a different value of $\kappa$ as our first calculation. We have noted above that varying $\kappa$ has little effect on the shape of the outline at constant $\bar{p}$. So, as long as the new value of $\kappa$ does not fall too much beyond the range of magnitudes where we verified this to be true (supplementary Fig. S1), the solution will be similar to that in equations (9, 10), plotted in Figure 2. Since the differences are small, and would be dwarfed by manufacturing variability in a real system, we kept the same shape (9)–(10) in all our atomistic transport calculations to be discussed below for the sake of comparability. This shape is henceforth referred to as the flower device. The initial optimization was for a pressure of $p=100\,\textrm{bar}$, and we performed further pressurizations of the shape up to $5\,\textrm{kbar}$, for which the distribution of PMF was not greatly changed, though its magnitude was quite different. Each of these pressurizations is identified below by the maximum strain achieved in the graphene sheet, up to a maximum strain of $\varepsilon_{m}=6.11\%$ (note, however, that, even though our analysis is based on a set of specific parameter values, the scaling relations (11) allow the freedom to explore a large range of relevant pressures, target PMFs, or system sizes without having to run a new shape optimization procedure for each particular choice; see, for example, Fig. S6 in the Supporting Information).

4 Transport characteristics

To access the transport properties of the flower device we use the Landauer–Büttiker formalism; specifically, we are interested in the quantum conductance of the device, which is calculated by Caroli’s formula 24 $G=\frac{2e^{2}}{h}\operatorname{Tr}[\Gamma_{q}G^{r}\Gamma_{p}G^{a}]$, where $G^{r}=[G^{a}]^{\dagger}=[E+i\eta-H-\Sigma_{p}-\Sigma_{q}]^{-1}$ is the retarded [advanced] Green’s function, the coupling between the contacts and the device is represented by $\Gamma_{q}=i[\Sigma_{q}-\Sigma_{q}^{\dagger}]$, and $\Sigma_{q}$ is the self-energy of the contact $q$. The geometry of the device is easily included in a nearest-neighbors $\pi$-band tight-binding Hamiltonian by merely selecting the sites that are located between the mathematically defined inner and outer edges [cf. Figure 2]. The inner contact is essentially a circular contact with $R_{\textrm{in}}\approx 47\,\textup{\AA}$; the outer contact radius varies between $200\,\textup{\AA}\leq R_{\textrm{out}}\leq 322\,\,\textup{\AA}$ with the angular modulation prescribed by eq. (9). The calculation of the conductance requires the Green’s functions of these non traditional contacts; this problem can be resolved by recalling that, for graphene, the particular details of the contacts can be replaced by an effective self-energy term provided that the contacts inject a high number of modes (highly doped contacts) 25. Under this model, an effective self-energy term $-i|t|$ is added to the onsite energy of the atoms at the edges ($t=-2.7\,\textup{eV}$ is the graphene hopping parameter).

The mechanical deformations induced by the hydrostatic pressure are incorporated in the tight binding Hamiltonian through the modification of the hopping parameter between neighboring sites as 26

[TABLE]

where $\widehat{n}_{i}$ is the vector normal to the surface at point $i$, $\vec{d}_{ij}$ is a vector linking the two sites, and $V_{pp\sigma}$ and $V_{pp\pi}$ are the Slater–Koster parameters, modified to account for changes in the bond length through $V_{pp\pi}(d_{ij})=t\,e^{-3.37(d_{ij}/a-1)}$ and $V_{pp\sigma}(d_{ij})=-2.4V_{pp\pi}(d_{ij})$ 27, 28, 29.

To best appreciate the features that arise from the pressure-induced PMF, it is useful, first, to have a perspective over the main transport characteristics of unpressurized devices which define our scenario of reference; the results of these calculations are provided in the Supporting Information. Devices were generated from the optimized flower shape shown in Figure 2 ($R_{\textrm{in}}\approx 47\,\textup{\AA}$ and $200\,\textup{\AA}\leq R_{\textrm{out}}\leq 322\,\textup{\AA}$). Figure 4a shows the conductance of pressurized flower devices for different values of the maximal strain in the graphene flake, ${\varepsilon_{\text{m}}}$, that we use to distinguish different cases. Notably, at ${\varepsilon_{\text{m}}}\approx 6\,\%$, the conductance starts developing a plateau of $5.7(2e^{2}/h)$ with a resonant peak at $E=0.09t$; the plateau occupies the whole energy range $0\lesssim E\lesssim 0.09t$. When the corresponding density of states (DOS) is analyzed [cf. Figure 4b], one identifies two sharp peaks at $E=0$ and $E=0.09\,t$ arising from the formation of strain-induced LL with $n=0,\,1$. From these, we extract a PMF $B_{s}\approx 65\,\textup{T}$ and this estimate allows us to further confirm that local maxima occurring in the conductance of Figure 4a correlate to higher energy LL ($E_{2}=0.13t$, $E_{3}=0.16t$ and $E_{4}=0.18t$). The dependence of the conductance on device size for the same spatial distribution of PMF can be studied without having to run a new shape optimization procedure by exploring the scaling implied by eq. (11). For example, if one scales the inner and outer boundaries as $R^{\alpha}_{\textrm{in}}(\theta)=\alpha R_{\textrm{in}}(\theta)$, $R^{\alpha}_{\textrm{out}}(\theta)=\alpha R_{\textrm{out}}(\theta)$, these relations can be used to extract the pressure needed to achieve the same PMF pattern, as well as the magnitude of the quasi-uniform field (see Fig. S5 and related discussion in the Supporting Information). The conductance of a device thus scaled with $\alpha=3$ ($R_{\textrm{out}}\simeq 80\,\textup{nm}$) is shown in Figure 4c for maximal strains of ${\varepsilon_{\text{m}}}=2.48\%$ and ${\varepsilon_{\text{m}}}=5.16\%$. In the two cases both the plateau and resonant peak are again present, from which we estimate $B_{s}\approx 9\,\textup{T}$ and $B_{s}\approx 20\,\textup{T}$, respectively. However, in larger devices these features are much more sharply defined and the plateau is flatter: the conductance is therefore more perfectly quantized at $6(2e^{2}/\hbar)$ with growing device dimension. They also require less strain to emerge due to the fact that the inner and outer contacts become decoupled at lower PMF (because the condition that the magnetic length, $\ell_{B}{\,=\,}\sqrt{\hbar/eB}$, is much less than $L$ is met at smaller PMF).

4.1 Perfectly conducting channels and valley filtering

In order to understand the origin of this robust quantization induced by pressure over a large range of energies that are experimentally relevant, let us look first at the real-space distribution of the electronic wavefunctions. The local DOS (LDOS) calculated at the plateau midpoint ($E=0.04t$) and at the peak ($E=0.09t$) are shown in Figures 4d and 4e, respectively, for the device with $\alpha=1$ pressurized to ${\varepsilon_{\text{m}}}=6.11\%$ [cf. Figure 4a]. At the peak [$E=0.09t$, Figure 4e], the observation of a state entirely confined in the central portions of the device identifies it as clearly associated with one of the PMF-induced LLs 19. In sharp contrast, for energies in the plateau [$E=0.04t$, Figure 4d] the LDOS concentrates on six narrow radial channels of the flower device. Comparing the LDOS with the PMF map shown in Figure 2 reveals that the wavefunction concentrates at precisely — and only — those regions where the PMF changes polarity. Consequently, the plateau in conductance at low energies is associated with current being carried by these pseudomagnetic interfacial states, similarly to the corresponding situation in nonuniform (real) magnetic fields 30 where electrons propagate chirally along a boundary separating fields of opposite polarity.

A useful picture of the nature of these 1D modes stabilized at the polarity inversion interfaces is provided by the semi-classical limit of this problem where they become so-called snake states 31. The designation arises from the different winding sense of the electron orbits in regions where $B$ has opposite polarity, which allows them to be confined and propagate along the interface with a definite direction. This interpretation helps one to understand: (i) directionality of the states and (ii) localization in regions where $B_{s}\approx 0$ and changes sign. However, despite the usefulness of the semi-classical perspective, it is important to retain that no clear correspondence to the classical motion can be established in the case shown in Figure 4d because, as we are looking at energies between the $n=0$ and $n=1$ LL, the cyclotron radius is smaller than the magnetic length for the energies studied: $r_{c}=\ell_{B}^{2}k_{F}\sim\ell_{B}/\sqrt{2}<\ell_{B}$ 32, 33; these modes are in the quantum regime.

At low-energies, transport takes place via one-dimensional perfectly conducting modes, which is entirely consistent with the 1D-like energy dependence of the DOS plotted in Figure 4b for ${\varepsilon_{\text{m}}}=6.11\%$. Moreover, the six neighboring sectors of alternating PMF polarity in our flower device are expected to beget six snake state channels, so one may deduce the existence of a plateau in the conductance of $6(2e^{2}/h)$, which is very approximately the case seen in Figures 4a and 4c. Note that, while one might think that the difference between the observed [$5.7\,(2e^{2}/h)$] and theoretical values [$6\,(2e^{2}/h)$] of the conductance plateau is due to the inhomogeneity of PMF over each sector, a direct comparison between the conductance of the pressurized device with an unpressurized flower under a real magnetic field with artificially sharp polarity changes shows no difference. More specifically, we have calculated the conductance of the same unpressurized device with an external magnetic field of constant magnitude ($B=65\,\textup{T}$) but alternating in sign at the six places where the PMF does so [cf. Figure 2 and Figure 5a] (additional details and different cases are discussed in the Supporting Information). Figure 4f demonstrates that the conductance of the pressurized device with ${\varepsilon_{\text{m}}}=6.11\%$ is essentially the same as that of the idealized unpressurized device with real magnetic field of the same magnitude. The agreement is nearly perfect in the plateau region below the first LL, further corroborating that the flower device is the optimal solution (to emphasize this point, in the supplementary Fig. S4(e) we show that the conductance of pressurized annular devices that retain a perfectly circular shape does not develop the quantization plateau). It also validates our earlier argument regarding the origin of the chiral 1D channels in the pressurized devices and their close analogy to the semi-classical snake states. In the end, we can trace the small deviation from the ideal $6\,(2e^{2}/h)$ quantization in the pressurized devices of Figure 4a to a small coupling between the modes belonging to different radial segments: to avoid such coupling requires $R_{\textrm{in}}\gtrsim 2\ell_{B}$ so that the wave functions of neighboring modes do not overlap near the inner contact, whereas this device has $R_{\textrm{in}}\approx 4.7\,\textup{nm}$ and $\ell_{B}=3.1\,\textup{nm}$ which falls short of it. In the supplementary Fig. S4 we show that perfect quantization is indeed achieved when $R_{\textrm{in}}$ and $\ell_{B}$ fulfill that criterion and that, as expected more generally, the conductance step is quantized at $G=2n_{s}\,e^{2}/h$, where $n_{s}$ is the number of polarity changing interfaces in the disk.

Given this scenario, the way in which PMFs couple to the electronic motion in graphene has an interesting and useful implication. Since electrons belonging to distinct valleys feel PMFs with opposite sign 4, 5, they propagate in opposite directions along the polarity boundary of the field. This effect is schematically illustrated in Figures 5a for an electron from valley $K$, and 5b for valley $K^{\prime}$. Looking back at Figure 4d, each radial segment along which the LDOS is peaked supports two chiral 1D modes, each propagating in opposite radial directions, and each associated with one given valley. Suppose then that the outer contact does not cover the entire perimeter but is, instead, sectioned into six outer contacts that allow collection of current only from the vicinity of the 1D channels (for our geometries that would mean defining contacts in the vicinity of the indentations of the outer perimeter, as in Figures 5a and 5b). Under a given bias between inner and outer contacts (say, electrons flowing radially outward), the electrons arising from valley $K$ would reach three specific contacts (all equivalent, oriented $2\pi/3$ apart from each other), whereas those associated with valley $K^{\prime}$ would reach the other three contacts, as indicated in Figures 5a and 5b. In this schematic, electrons from valley $K$ are channeled by the 1D modes from the inner contact “I” to the outer contacts 1, 2 and 3, while those from valley $K^{\prime}$ can only reach contacts 4, 5 and 6. A change in the bias sign, or changing from electron to hole-doped graphene with a back gate, would exchange the valleys that “reach” a given contact. In this manner, the device can spatially separate individual contributions from each valley and deliver valley polarized current to specific contacts or, alternatively, selectively filter or probe the existence of valley polarized currents.

We highlight that the observation of these effects in our flower-like geometry is not limited to the specific range of hydrostatic pressure and device dimensions used for the calculations reported in Figure 5. For example, using the scaling relations (11), we expect that scaling up the device to $L=400\,\textup{nm}$ while preserving its shape and using a pressure of $\sim 5.5\,\textup{bar}$ shall create a PMF of $\sim 0.26\,\textup{T}$ with the same spatial distribution. The valley filtering effect will take place as long as LLs can be formed in the center of each petal guaranteeing that transport only takes place via the 1D modes. A simple estimate of the electronic mobilities needed for this can be made from the criterion that the mean free path, $\ell$, be larger than the characteristic cyclotron diameter: $\ell\gtrsim 2\ell_{B}$. In terms of the mobility defined as $\mu{\,=\,}\ell e/\hbar k_{F}$ and using $k_{F}{\,=\,}\sqrt{2}/l_{B}$ for electrons in the first LL, this becomes $\mu\gtrsim 2\ell_{B}e/\hbar k_{F}{\,=\,}\sqrt{2}\ell_{B}^{2}e/\hbar{\,\sim\,}5{\times}10^{4}\,\text{cm}^{2}/\text{V\,s}$, where we used $l_{B}{\,=\,}\sqrt{\hbar/eB}{\,\approx\,}50\,\textup{nm}$; such mobilities are rather standard in graphene devices (see Fig. S6 and related discussion in the Supporting Information for additional details).

Recent research has shown the existence of snake states in strained graphene nanoribbons PhysRevLett.112.096805, PhysRevB.94.075432, and diverse valley filtering devices have been proposed by exploiting strained graphene 34, 35, 36, 37, 38, 39. Similarly to ours, these devices require a specific geometry. However, ours does not require external electric or magnetic fields, produces valley filtered currents equally for both valleys, and whether the filtering occurs or not can be controlled by pressure or strain.

4.2 Effects of disorder and high on/off ratios

As often happens, disorder can simultaneously be detrimental or functional in solid state systems. From the first point of view, the conductance plateau associated with the snake states should be affected by disorder in real samples. As each radial channel supports one pair of opposite chirality and opposite valley quantum number, short-range scatterers will be particularly disadvantageous since they cause inter-valley scattering and, because of the valley-chirality locking, that translates into backscattering and the degradation of their perfectly conducting nature. Additionally, the snake states’ locations coincide with the sites of highest curvature in the pressurized flower where the probability of adsorbing alien atoms or molecules is higher 40, 41. In order to assess this in a specific manner, we modeled the attachment of adatoms to graphene in these regions by adding an Anderson-impurity type perturbation to the electronic Hamiltonian: each carbon can bind an adsorbate with probability $p_{ad}$; this introduces a new electronic level of energy $E_{ad}=t/16$ that the electron can hop to from the attached carbon with a hopping amplitude $V=2t$ that describes the degree of hybridization 42. The conductance averaged over an ensemble of 20 such systems with $p_{ad}=0.1$ for the device of Figure 4a is plotted in Figure 5c. As anticipated, the presence of these dopants in the regions where the snake states occur increases backscattering, thus lowering the conductance. The plateau is seen to remain, but its height depends inversely on the size of the system (not shown); in this particular case, traces of quantization are still present due to the small width of the device used (${W_{\text{eff}}}\approx 223\,\textup{\AA}$). We point out that the approximate preservation of the plateau under strong disorder (it affects 10 % of the atoms) is further evidence of the 1D character and resilience of the underlying modes: there is now backscattering (the plateau has smaller $G$), but the one-dimensionality is unaffected (a plateau largely remains).

Whereas such a system with strong short-range scattering will be less effective as a valley filter, the sensitivity of the low-energy base conductance to the radial dimension can have practical functional applications, as we now describe. Under pressure, the perfectly conducting channels limit the minimum conductance that can be achieved, irrespective of how large the magnitude of the PMF might be. From this standpoint, they are leaky electronic devices without an off state. But the residual conductance can be controlled by the type and amount of disorder, suggesting a tuneable on-off ratio. To simulate larger devices without incurring a high computational cost, and given that all of the above establishes the one-dimensional character of transport at low energies, we resort to a simple one-dimensional model 43 to scale up the conductance for larger pressurized disordered devices:

[TABLE]

This assumes that in the snake-type 1D modes, electrons transmit through one adatom with probability $T$ and $N=\left({W_{\text{eff}}}/a\right)p_{ad}$ is the number of scatterers in the channel of width ${W_{\text{eff}}}$. It follows that wide devices should have a notorious reduction of the conductance, entailing high ratios of unpressurized/pressurized conductance (${G_{\text{up}}}/{G_{\text{p}}}$) in the energy region of the snake states. This ratio can be used as a figure of merit directly related to the on-off ratio of the device. To quantify ${G_{\text{up}}}/{G_{\text{p}}}$, we employed a number of modeling assumptions. Firstly, the Fermi energy was fixed at $E_{F}=E_{1}/2$, corresponding to one half of the energy of the $n=1$ LL, the midpoint of the conductance plateau. Secondly, a perfectly round inner contact and ballistic transmission were considered in order to estimate ${G_{\text{up}}}=(4e^{2}/h)2j_{m}$, where $j_{m}\sim k_{F}R_{\textrm{in}}$ is the maximum total angular momentum 25. Third, geometrical factors such as $R_{\textrm{in}}$ and ${W_{\text{eff}}}$ were scaled from our original device ($\alpha=1$). Finally, the value $T\approx 0.9$ was extracted from the disordered conductance of Figure 5c, from where we obtain $\mathcal{T}\approx 0.6$. Figure 5d shows that, in pristine devices, ${G_{\text{up}}}/{G_{\text{p}}}$ decreases as $\propto L^{-3/2}$. This decrease is due solely to ${G_{\text{up}}}$ because, in a clean device, ${G_{\text{p}}}$ remains at the value $\approx 6(2e^{2}/h)$ for any $L$ as long as $E_{F}=E_{1}/2$. But, since the PMF scales inversely with $L$ for a fixed device shape, to keep this choice of $E_{F}$ requires its value to vary for different $L$. Hence, each $L$ in the figure corresponds to a different $k_{F}$, and these parameters vary inversely to each other; this causes the overall scaling shown by the black line in the figure. However, in the presence of the adsorbates, the 1D transport regime that occurs under pressure is expected to be much more sensitive to the disorder. The rapid decrease of the transmission probability $\mathcal{T}$ with $L$ in this case dominates over the variation of ${G_{\text{up}}}$. As a result, the figure of merit, which is governed by the overall transmission probability through ${G_{\text{up}}}/{G_{\text{p}}}\propto j_{m}/\mathcal{T}$, increases with the channel length reaching values near 12 for a device with $L=400\,\textup{nm}$ ($R_{\textrm{in}}=70\,\textup{nm}$, $\alpha=15$), a high “on/off” value considering the absence of a band gap in graphene.

5 Conclusions

In this work we set out to determine the optimal geometry of a graphene Corbino device that guarantees a PMF of nearly constant magnitude throughout most of the system when externally pressurized. Since the sign of the PMF has to inevitably alternate six times along any closed path containing the inner contact, the field cannot be strictly constant in the whole ring 8. Yet, both its magnitude and absolute value can be made satisfactorily uniform. The first case is discussed in the Supporting Information and we see that it is possible to require a specific sign and magnitude in most of the system, and the effect of the optimization process is to generate a geometry where the regions with the opposite PMF are much reduced (a more extreme spatial reduction is possible under higher pressures and higher target PMFs). If, on the other hand, the goal is only the PMF magnitude, but not its sign, the optimal geometries are the sixfold flower shapes such as in Figure 2 that we have analyzed in detail. The spatial boundaries where the PMF changes sign for the latter are sharper and better defined, which leads to robust snake states and strictly one-dimensional transport over a range of pressures. This radially one-dimensional regime is signaled by the strict quantization of the conductance at $2e^{2}/h$ per channel that is seen to survive up to high values of $E_{F}$ (e.g. 50–100 meV, cf Figure 4), the characteristically 1D behavior of the DOS as a function of energy there, and corroborated by the real space profile of the LDOS.

The strain-induced perfectly conducting channels can be exploited in two different directions for electronic applications. On the one hand, they limit from below the current pinch-off effect and make these devices leaky because of their chiral nature. Our studies of the effect of disorder show that short-range defects can activate inter-valley scattering and reduce the residual conductance . Since their spatial location is set by the geometry, we can envisage this being done in a deliberate way through adsorbates, for example, so that only the regions where snake states develop under pressure become disordered. In this way, the conductance in the unpressurized state would remain unaffected, but the perfectly conducting channels would no longer exist under pressure which would significantly boost the on/off ratio. The pressure sensitivity and its direct translation into current modulations, suggests its possible application in electromechanical sensing or transducers 44. On the other hand, we have also seen an equally interesting perspective where snake-type states and their leaky residual conductance are not detrimental but, instead, functional: clean devices can be employed as sources of valley-polarized currents in graphene. Note that the Corbino ring (the device) is defined only by the indentation of the substrate, not by an actual patterning of graphene, and the inner and outer contacts at $R_{\textrm{in}}$ and $R_{\textrm{out}}$ are still the same sheet of graphene. Hence, the sheet can extend over large distances beyond the outer radius and these structures can act as local sources of valley-polarized currents for injection into the two dimensional graphene plane beyond $R_{\textrm{out}}$.

{acknowledgement}

DAB acknowledges the support from Mackpesquisa and FAPESP under grant 2012/50259-8, VMP that of the Singapore Ministry of Education Academic Research Fund Tier 2 under grant number MOE2015-T2-2-059, and AHCN the support of the National Research Foundation of Singapore under the Mid Size Centre Grant. Numerical computations were carried out at the HPC facilities of the NUS Centre for Advanced 2D Materials.