Twin lead ballistic conductor based on nanoribbon edge transport

TL;DR

This paper models a quantum device based on graphene nanoribbons with four terminals, analyzing how two electrical circuits interact through quantum and classical effects, providing analytical and numerical insights into their coupled behavior.

Contribution

It introduces a theoretical framework for analyzing coupled quantum-classical circuits using a nanoribbon-based device with a detailed conductance matrix and self-consistent algorithms.

Findings

Identified regimes where the device supports two nearly independent currents.

Derived analytical formulas for linear-response operation.

Developed a recursive algorithm for non-equilibrium Fermi level calculations.

Abstract

If a device like a graphene nanoribbon (GNR) has all its four corners attached to electric current leads, the device becomes a quantum junction through which two electrical circuits can interact. We study such system theoretically for stationary currents. The 4-point energy-dependent conductance matrix of the nanostructure and the classical resistors in the circuits are parameters of the model. The two bias voltages in the circuits are the control variables of the studied system while the electrochemical potentials at the device's terminals are non-trivially dependent on the voltages. For the special case of the linear-response regime analytical formulae for the operation of the coupled quantum-classical device are derived and applied. For higher bias voltages numerical solutions are obtained. The effects of non-equilibrium Fermi levels are captured using a recursive algorithm in which self-consistency between the electrochemical potentials and the currents is reached within few iterations. The developed approach allows to study scenarios ranging from independent circuits to strongly coupled ones. For the chosen model of the GNR with highly conductive zigzag edges we determine the regime in which the single device carries two almost independent currents.

Tables

Table 1. Table 1: Off-diagonal elements G α ​ β subscript 𝐺 𝛼 𝛽 G_{\alpha\beta} of four linear-response conductance matrices in the units of G 0 subscript 𝐺 0 G_{0} . The first row of the numerical values describes the nanojunction based on the intact (perfect) GNR (Fig.  1 ) in the full equilibrium, i.e. without any bias or gate voltage applied. The second row describes the case of the perturbed GNR, which is specified by Fig.  1 as well, with a detail shown by the inset to Fig.  3 , again in the full equilibrium. Rows 3 and 4 provide linear conductances for the same two GNRs with the gate voltage applied upon them so that the transmission is resonantly enhanced due to the ZZ edge states (Sect.  4.3 ). The unlisted off-diagonal elements can be deduced from symmetries.

	$G_{12}$	$G_{13}$	$G_{14}$	$G_{24}$
Intact	$0.1840$	$0.1246$	$0.1816$	$0.1246$
Perturbed	$0.1269$	$0.02575$	$0.1297$	$0.08584$
Intact+$V_g$	$1.012$	$0.07974$	$0.1567$	$0.07974$
Perturbed+$V_g$	$0.9180$	$0.07857$	$0.1859$	$0.1679$

Equations

$$I_{\alpha }=-\frac{2e}{h}\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^{N_{\mathrm{T}}}{\int }_{{\mu }_{\beta }}^{{\mu }_{\alpha }}{\mathcal{T}}_{\alpha \leftarrow \beta }(\mathcal{E})\mathrm{d}\mathcal{E},\alpha \in \{1,2,\dots ,N_{\mathrm{T}}\}.$$

$$\sum _{\beta (\ne \alpha )}{\mathcal{T}}_{\alpha \leftarrow \beta }(\mathcal{E})=\sum _{\beta (\ne \alpha )}{\mathcal{T}}_{\beta \leftarrow \alpha }(\mathcal{E}).$$

$${\mu }_{\alpha }\equiv -e{U}_{\alpha }.$$

$$I_3+I_4=0,$$

$$I_1+I_2=0.$$

$$U_2$$

$$U_4$$

$$U_1=0$$

$$I_{\alpha }=\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^{N_{\mathrm{T}}}{G}_{\alpha \beta }(U_{\alpha }-U_{\beta }).$$

$$I_{\mathrm{A}}$$

$$I_{\mathrm{B}}$$

$$\gamma ={\alpha }_{\mathrm{AA}}{\alpha }_{\mathrm{BB}}-{\alpha }_{\mathrm{BA}}{\alpha }_{\mathrm{AB}}.$$

$$\displaystyle U_{4}-U_{2}=-\frac{1}{SD(G,R)}\,\bigl{[}$$

$$\displaystyle\bigl{.}-$$

$$D(G,R)=D=1+{\alpha }_{\mathrm{AA}}{R}_{\mathrm{A}}+{\alpha }_{\mathrm{BB}}{R}_{\mathrm{B}}+\gamma R_{\mathrm{A}}{R}_{\mathrm{B}}$$

$$u\equiv \frac{\mathcal{E}}{e},$$

$$G_{\alpha \beta }(u)\equiv \frac{2e^2}{h}{\mathcal{T}}_{\alpha \leftarrow \beta }(eu).$$

$$-\sum _{\begin{array}{c}\beta =1\\ \beta \ne 1\end{array}}^4{\int }_{{\mu }_{\beta }/e}^{{\mu }_1/e}{G}_{1\beta }(u)\mathrm{d}u$$

$$-\sum _{\begin{array}{c}\beta =1\\ \beta \ne 2\end{array}}^4{\int }_{{\mu }_{\beta }/e}^{{\mu }_2/e}{G}_{2\beta }(u)\mathrm{d}u$$

$$-\sum _{\begin{array}{c}\beta =1\\ \beta \ne 3\end{array}}^4{\int }_{{\mu }_{\beta }/e}^{{\mu }_3/e}{G}_{3\beta }(u)\mathrm{d}u$$

$$-\sum _{\begin{array}{c}\beta =1\\ \beta \ne 4\end{array}}^4{\int }_{{\mu }_{\beta }/e}^{{\mu }_4/e}{G}_{4\beta }(u)\mathrm{d}u$$

$${\mu }_2={\mu }_1-e{V}_{\mathrm{A}}+e{R}_{\mathrm{A}}{I}_{\mathrm{A}},{\mu }_4={\mu }_3-e{V}_{\mathrm{B}}+e{R}_{\mathrm{B}}{I}_{\mathrm{B}},{\mu }_1\equiv E_{\mathrm{F}}=0.$$

Full text

Keywords: conductance; multiterminal; non-equilibrium; transport; edge; carbon nanostructures

1 Introduction

Electronic transport through graphene nanoribbons has attracted attention in particular because their zigzag edges exhibit magnetic properties and support specific modes, which may influence the charge transport [1, 2, 3]. Edges of larger-scale graphene ribbons have been proposed to serve similarly or analogously to optical fibres or waveguides [4]. In experimental works studying electronic transport through a nanoribbon, the flake is typically contacted to electrodes across its whole width, see for instance Ref. [5]. However, if the edges are to be used for the transport, it may be appropriate to make contacts only to the corners of such a ribbon. This poses a question if a single graphene nanoribbon or an alternative planar atomistic structure could serve as two wires carrying two independent electric currents. Obviously, there would be some interaction between the currents flowing along the two parallel edges and they would not be fully independent. Still, if the energy ranges of the high density of states, resulting from the zigzag (ZZ) edge modes, are used for the transport, the conductance along the ZZ direction may be significantly larger than the conductance along the perpendicular (armchair, AC) direction. In such case the two ZZ edges might be considered as two relatively independent spatially separated conductors. More generally, such scenario would be the case if the $4\times 4$ conductance matrix of a ribbon preferred one of the directions over the perpendicular and diagonal directions. On the other hand, a nanoribbon with significant conductance along several directions would allow to study the intriguing regime in which the two classical circuits are coupled through the quantum ballistic device. Regardless of the conductance characteristics, such setups would naturally require four electrodes contacted to the corners of the nanoribbon.

Electronic transport through graphene flakes with electrodes contacted at their corners has been studied computationally in Ref. [6], with focus on non-rectangular flakes such as trapezoids or triangles. The work has addressed the question how the two-contact conductance qualitatively depends on the magnitude of the angles in the two contacted corner areas and on the types of the edges forming the corner. It was in addition found that the usual nearest-neighbor (NN) tight-binding (TB) model was insufficient for proper description of the edge-induced transport properties (provided that the electrodes were attached to the corners).

Four-terminal phase-coherent conductance has been considered by Büttiker [7] in order to clarify the occurrence of asymmetric magnetoresistances. The author considered a sample contacted to two circuits, what is a setup similar to the one assumed in the present work. On the contrary, the model of Büttiker does not include any classical resistors and is limited to the linear regime. We will return to it in the exposition of our theoretical description. Four-terminal schemes have many times been studied both theoretically and experimentally in setups with two voltage probes implying that two of the four terminals carry zero net currents. Analysis of this specific model can be found as early as in the above cited work by Büttiker. A four-terminal electron waveguide coupler was proposed and theoretically analysed in Ref. [8], assuming ballistic transport. Based on the formulae of Ref. [7], Büttiker’s resistance tensor of the about $500\,$Å long junction has been calculated. In Ref. [9] an experimental realisation of two parallel wires coupled by a ballistic window was described. Molecular logic gates with intramolecular circuits are another example of nanoscale multiterminal devices. Considered schemes include both ballistic components as well as classical resistors [10, 11]. Theoretical analysis and computational results on multiple-lead coherent conductors in case of finite applied voltages were reported in Ref. [12]. A particularly studied case was a four-terminal conductor with two of the currents set to zero. Voltage difference between the two open leads as a non-linear function of the current was calculated.

In this paper we study a coupled quantum-classical device consisting of a ballistic transport junction and of two classical circuits with the four leads contacted at the junction. As a concrete realisation of the junction a GNR is considered. The two electric currents are driven by DC bias voltage sources in the circuits that include resistors. For the linear regime, transparent analytical formulae displaying the interplay of the quantum and classical elements are derived for both the electric currents as well as for the electrochemical or electrostatic potentials in the device’s leads. A numerical scheme is described to tackle the non-linear regime of the composite system. In order to include the effects of the non-equilibrium (NE) electrochemical potentials (EChPs) on the current-voltage characteristics and other quantities, a recursive method for the coupled device is developed. Numerical demonstrations are performed for both a perfect GNR and for a perturbed one. We show that under certain conditions the two electrical circuits of the scheme operate almost independently. Hence such a GNR could serve as two relatively independent conductors, which might be found useful in experimental setups and practical applications as it would allow to shrink the device size. The regime in which the circuits are significantly coupled through the ballistic device presents an interesting physical model and is cover by our numerical demonstrations as well.

In Sect. 2 we define the system under study. In Sect. 3 we derive main results of our work both for the linear and non-linear regimes. We also provide additional description of the studied model, relevant for the NE calculations. In Sect. 4 we apply the theoretical modelling to study the setup defined in Sect. 2. Additional details can be found in the Supplementary Material (SM) [13]. We conclude our results in Sect. 5.

2 The model and methods

The physical model.

An atomistic model of the considered GNRs is shown in Fig. 1. The red-marked carbon atoms in the corners form the contact areas where the four electrodes (not shown) are attached. If this is accomplished, the whole system in a simple electrical setup can be described by the scheme drawn in Fig. 2.

The resistances $R_{\mathrm{A}}$ and $R_{\mathrm{B}}$ of the two classical resistors are known parameters. The control variables are the static bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$ on the two DC sources, not the EChPs in the leads. The potentials as well as the two currents are unknown non-trivial functions of the bias voltages and have to be determined as it is explained in detail below in the text. The stationary regime is assumed for all elements of the scheme. We consider a low-temperature regime in which inelastic scattering of the electrons in the GNR and in its contacts can be neglected. Hence the GNR with its contacts represents a quantum subsystem of the entire scheme. The electronic transport properties of such a contacted device in the linear-response regime can be described by a $4\times 4$ conductance matrix $G$. This matrix is computed using a microscopic model described below. For sufficiently low bias voltages the linear regime with $I_{\mathrm{A}},I_{\mathrm{B}}\propto V_{\mathrm{A}},V_{\mathrm{B}}$ can be assumed and the conductance matrix is then considered as energy-independent. This will allow us to derive analytical formulae for currents $I_{\mathrm{A}}(V_{\mathrm{A}},V_{\mathrm{B}})$ and $I_{\mathrm{B}}(V_{\mathrm{A}},V_{\mathrm{B}})$ in terms of the fixed system parameters $G$, $R_{\mathrm{A}}$, and $R_{\mathrm{B}}$.

We recall that in the scheme under study (Fig. 2) not only the currents but also the EChPs $\mu_{2}$, $\mu_{3}$ and $\mu_{4}$ are unknown non-trivial functions, $\mu_{\alpha}(V_{\mathrm{A}},V_{\mathrm{B}})$, which have to be determined and must strictly be distinguished from the bias voltages. This requirement is dictated by the scheme itself in which the differences $\mu_{2}-\mu_{1}$ and $\mu_{4}-\mu_{3}$ are a priori unknown because the currents through the resistors are unknown. The chemical potential at the terminal 1 is fixed by the ground-level condition. One of the foci of the present work is the non-linear response regime with generally non-equilibrium EChPs $\mu_{\alpha}$ in individual leads. In the non-linear regime the response of the junction will be described by an energy-dependent conductance matrix $G(\mathcal{E})$ instead just a simple numerical matrix. If, in addition, the NE values of $\mu_{\alpha}$’s are assumed, the conductance matrix will acquire the additional dependences: $G=G(\mathcal{E};\{\mu_{\alpha}\})$.

The microscopic description.

We consider GNRs with the electrodes attached to their corners according to Fig. 1. The model of the electrodes and the electronic structure description are similar to those used in Ref. [6] with a notable exception of the more general model of the electrodes in the present work. We use the independent-electron approximation and an extended tight-binding (TB) hamiltonian for the GNR with interactions up to the 3rd NN included. The hopping parameters $t_{\mathrm{B}}=-2.97\,\mathrm{eV}$, $t_{\mathrm{B}}^{\prime}=-0.073\,\mathrm{eV}$ and $t_{\mathrm{B}}^{\prime\prime}=-0.33\,\mathrm{eV}$ are taken from Ref. [15]. The orbital overlaps are neglected. The extended interaction range within the GNR was found important to account for the contributions of the ZZ edge state to the electronic transport if the electrodes are attached to corners of graphene flakes [6]. We do not consider the spin-related effects; the focus of the present work is the description of a four-terminal nanoribbon device coupled to two classical branches of the circuits. A generalisation of our formalism including the spin effects at least on the level of the mean-field approximation to Hubbard hamiltonian would be relatively straightforward although computationally much more expensive [16]. Quantitatively accurate modelling of GNR’s electronic structure especially for narrow ribbons would also require an inclusion of the electron-electron interactions into the model [14].

Each electrode in our model is formed by a bunch of mutually non-interacting monoatomically thin wires [17, 6] described within the $2$nd NN approximation with the couplings $t_{1}=t_{\mathrm{B}}$ and $t_{2}=t_{\mathrm{B}}^{\prime}$. The particular model drawn in Fig. 1 and used for our quantitative analyses employs $54$ wires per electrode. Using just a single wire would allow us to study very small currents only. The introduction of the $2$nd NN couplings $t_{2}$ for the atoms of the wires is the only important modification of the used microscopic theory in comparison to work [6]. The equilibrium on-site (atomic orbital) energies are denoted as $\epsilon$ and are assumed to be the same for the entire system, including the electrodes. The extended TB model of the monoatomically thin wires yields the one-electron dispersion relation $\mathcal{E}(\mathcal{K})=\epsilon+\sum_{n=1}^{N_{\mathrm{far}}}2t_{n}\cos(n\mathcal{K})$, with $N_{\mathrm{far}}$ being farthest nearest neighbor considered in the wires; here $N_{\mathrm{far}}=2$. The Fermi energy of the wires is $E_{\mathrm{F}}=\mathcal{E}(\pi/2)=\epsilon-2t_{2}$. Introduction of the extended TB model of the electrodes thus causes the positive shift $-2t_{2}$ of the Fermi energy compared to the basic TB model. We set $\epsilon=2t_{2}=2t_{\mathrm{B}}^{\prime}$ what results in having $E_{\mathrm{F}}=0$.

The couplings between the wires and the GNR again extend up to the $2$nd NN from the vertex atom of the GNR (the atom of the GNR to which a given wire is attached). In other words, the vertex atom is coupled to two closest atoms of the semi-infinite wire. Our model assumes the same coupling parameters $t_{1}$ and $t_{2}$ as defined above for the intra-wire interactions.

In the present work we pay attention to the non-linear response regime of the coupled quantum-classical device including the transport with an account of NE Fermi levels. Techniques to tackle these regimes, given the specifics of the studied scheme, are described in Sects. 3.3.3 and 3.3.4. In this way the specification of our model from Sect. 2 will be supplemented to be usable also under the NE conditions.

3 Theory

3.1 Quantum-mechanical calculations

Standard Green’s function formalism (described for instance in Ref. [18]) is used to calculate eigenfunctions of the entire system. Knowledge of the eigenfunctions allows us to compute the transmission and reflection coefficients for electrons incoming from the wires. With the exception of the extended TB model of the wires the calculations have been accomplished according to the route described in Ref. [6]. The Green’s function formalism requires knowledge of the semi-infinite wire’s eigenstates which are provided in the SM [13].

3.2 Landauer-Büttiker formula

We will employ the well-known Landauer-Büttiker (LB) type formula for stationary currents through the leads of a multi-terminal device:

[TABLE]

$N_{\mathrm{T}}$ is the number of the terminals (also the electrodes), $\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E})$ is the transmission coefficient from the terminal $\beta$ to the terminal $\alpha$ and $\mu_{\gamma}$ are the EChPs associated with the individual wires. Formula (1) assumes (i) zero temperature, (ii) the time-reversal symmetry (TRS) and (iii) the trivial description of the spin degrees of freedom (spin degeneracy). The TRS implies that

[TABLE]

We note that in our model each electrode (attached at its corresponding terminal) is assumed to be composed of a number of the mono-atomically thin wires. We assume that $\alpha$ and $\beta$ label the electrodes, not the wires. Still, the TRS formula (2) is equally exactly valid for these composite electrodes as it would be valid if the indices $\alpha$ and $\beta$ labelled the individual wires attached on the central device.

The charge conservation for the central device implies Kirchhoff’s first law: $\sum_{\alpha}I_{\alpha}=0$. The EChP $\mu_{\alpha}$ at the terminal $\alpha$ is conveniently expressed in terms of the electrostatic potential $U_{\alpha}$ associated with the terminal (the electrode):

[TABLE]

We recall that in our model (see its description in Sect. 2 and in Fig. 2) the EChPs themselves as well as the currents are non-trivial unknown functions of the two bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$ as will be further explained below.

3.3 Inclusion of the circuits with resistors into the 4-terminal device

3.3.1 Formulation of the problem

Consider a device with 4 terminals, i.e. $N_{\mathrm{T}}=4$ like that drawn in Fig. 2. Each terminal has an electrode connected to it. In the basic statement of the LB formula (1) above the (stationary) currents $I_{\alpha}$ obey to Kirchhoff’s first law and no other constraints were supposed. Now we put in consequences the assumption that the four leads form the two circuits and there are also the two resistors included so that the entire system can be schematically depicted as on Fig. 2. The lower circuit (A) has one of its branches grounded. For such 4-terminal device we have to require the following additional constraining condition,

[TABLE]

which is more restrictive than just Kirchhoff’s first law. We assume the stationary regime for each element of the entire scheme, i.e. also for the sources. Therefore the condition (4) is necessary, otherwise an unphysical charge imbalance would be gradually built on the voltage source $V_{\mathrm{B}}$ and that would in addition contradict the assumption of the stationarity. Kirchhoff’s first law together with the Eq. (4) then imply a similar condition on the currents $I_{1}$ and $I_{2}$:

[TABLE]

The presence of the circuits requires that, in addition to the Eqs. (4) and (5), Kirchhoff’s second law has to be fulfilled:

[TABLE]

We recall the relation (3): $\mu_{\alpha}\equiv-eU_{\alpha}$. According to the scheme in Fig. 2, we can choose

[TABLE]

without any significant loss of generality. This auxiliary setting allows us to accomplish derivations using the potentials $U_{2}$, $U_{3}$ and $U_{4}$ instead of the differences $U_{2}-U_{1}$, $U_{3}-U_{1}$ and $U_{4}-U_{1}$. Hence the three electrostatic potentials $U_{2}$, $U_{3}$, $U_{4}$ and the four currents $I_{\alpha}$ are unknown quantities and have to be determined as functions of the two given bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$.

3.3.2 The linearity approximation

Under the assumption of the section’s title, formula (1) for the electric currents takes the frequently accounted form

[TABLE]

We stress again that for the scheme under study (see Fig. 2) the electrostatic potentials $U_{\gamma}$ must not be reduced to or confused with the bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$. See also the reasoning in Sect. 2. The linear conductance matrix elements are $G_{\alpha\beta}\equiv(2e^{2}/h)\,\mathcal{T}_{\alpha\leftarrow\beta}$. The four equations (8) together with the constraints (4), (5) and with the Eqs. (6) form an inhomogeneous set of eight linear equations with $7$ unknowns $U_{2}$, $U_{3}$, $U_{4}$, $I_{1}$, $\dots$, $I_{4}$. Known parameters are the two resistances $R_{\mathrm{A}}$ and $R_{\mathrm{B}}$, the matrix elements $G_{\alpha\beta}$ satisfying the TRS of the form (2) and the two bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$. I.e. the two voltages, $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$, are supposed to be the control variables in such an experiment, not the three electrostatic potentials $U_{2}$, $U_{3}$, and $U_{4}$. There is a linear dependence among the 8 equations. Once the set is properly reduced to remove the dependence, we obtain an inhomogeneous set with a unique solution. The equations are solved in the SM [13] and resulting linear relations $I_{\mathrm{A}}=I_{\mathrm{A}}(V_{\mathrm{A}},V_{\mathrm{B}})$ and $I_{\mathrm{B}}=I_{\mathrm{A}}(V_{\mathrm{A}},V_{\mathrm{B}})$ are obtained:

[TABLE]

where $\alpha_{CC^{\prime}}$ are conductances introduced by Büttiker [7] and depend solely on the conductance matrix elements. For the problem in hand they take the form $\alpha_{\mathrm{AA}}=S_{1}-(G_{13}+G_{14})(G_{31}+G_{41})/S$, $\alpha_{\mathrm{AB}}=(G_{13}G_{24}-G_{14}G_{23})/S$, $\alpha_{\mathrm{BA}}=(G_{31}G_{42}-G_{41}G_{32})/S$ and $\alpha_{\mathrm{BB}}=S_{3}-(G_{31}+G_{32})(G_{13}+G_{23})/S$, with $S=G_{13}+G_{23}+G_{14}+G_{24}$ and $S_{\alpha}=\sum_{\beta\neq\alpha}G_{\alpha\beta}$. Alternative expressions are provided in the SM [13]. For convenience we have introduced the parameter

[TABLE]

In the SM [13] we find and present formulae for the all potentials differences $U_{\alpha}-U_{\beta}$. Among the other differences, we find

[TABLE]

as well as a more compact formula $U_{4}-U_{3}=[(1+\alpha_{\mathrm{AA}}R_{\mathrm{A}})V_{\mathrm{B}}+\alpha_{\mathrm{BA}}R_{\mathrm{B}}V_{\mathrm{A}}]/D$. The shorthand notation

[TABLE]

represents the commonly occurring denominator in the above expressions (9) for the currents as well as in formula (LABEL:eq:U42_solo) for the differences of the potentials. Expressions for the remaining differences $U_{\alpha}-U_{\beta}$ can be determined from the above stated formulae using the symmetries of the considered scheme.

The Eqs. (9) and (LABEL:eq:U42_solo) are our main analytical results for the linear regime, together with similar equations for the remaining $U_{\alpha}-U_{\beta}$ differences [13]. The Eqs. (9) generalise Büttiker’s formulae [7] and transparently show how the currents in the two circuits relate both to the classical and to the quantum elements of the scheme.

The Eq. (LABEL:eq:U42_solo) describes how the potentials difference $U_{4}-U_{2}$ (and analogously $U_{3}-U_{1}$) can be varied through a continuous range of values by simply tuning the classical resistors in the circuits. For example, we can obtain zero difference $U_{4}-U_{2}$, i.e. also the equality $\mu_{4}=\mu_{2}$ of the EChPs. [See the relations (3).] With one of the batteries reversed (for instance $V_{\mathrm{B}}\to-V_{\mathrm{B}}$), one can similarly try to set to zero at least one of the potentials differences across the sample as it can be seen from the results for the other $U_{\alpha}-U_{\beta}$ differences [13]. An illustration of the effect will be shown in Sect. 4.1.3. Obviously, a particular condition like $U_{4}=U_{2}$ is only achievable for proper combinations of the numerical values of the quantities entering the relation.

3.3.3 The 4-terminal device beyond the linear response

In this section we provide a generalisation of the considerations in Sect. 3.3.2 beyond the linear regime. It will allow us to more thoroughly describe nano-devices like junctions formed by graphene flakes that are often found to operate in non-linear regimes. The quantum-mechanical part of the system is parameterised by the matrix $\mathcal{T}$ standing in the Eq. (1), which is now a set of (presumably known) functions of the energy. The energy dependence is the substantial complication in comparison to the linear regime in which the matrix elements were constants. To make the formulae free of the numerical prefactors $2e^{2}/h$, we equivalently use an energy-dependent conductance matrix instead of $\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E})$. To keep the dimensions correct, we introduce an auxiliary quantity

[TABLE]

where $e>0$ is the unit of charge. The conductance matrix elements are then written as

[TABLE]

The equations to be solved are (1) with $N_{\mathrm{T}}=4$, again with the constraints (4) and (5) and Kirchhoff’s second law Eqs. (6). For convenience we collect and write down all these equations and conditions here, making their specific to the 4-terminal system under study:

[TABLE]

Accompanying relations (6) and (7) for the classical part of the system are now more conveniently expressed in terms of the EChPs $\mu_{\alpha}\equiv-eU_{\alpha}$ rather than the electrostatic ones:

[TABLE]

The accomplish the task we have to solve the coupled set of the above equations (15) and (16), most of them being integral equations. Given parameters and known functions are $R_{\mathrm{A}}$, $R_{\mathrm{B}}$, $V_{\mathrm{A}}$, $V_{\mathrm{B}}$ and $G_{\alpha\beta}(u)$. Main unknown physical quantities to be determined are the currents $I_{\mathrm{A}}$ and $I_{\mathrm{B}}$. Of secondary interest are the EChPs $\mu_{\alpha}$ which are also unknown with the exception of $\mu_{1}$. Obviously, all this is a much more difficult task than it was in the case of its linear counterpart in Sect. 3.3.2. For the non-linear case we develop a numerical algorithm which takes $\mu_{2}$, $\mu_{3}$, $\mu_{4}$, $I_{\mathrm{A}}$ and $I_{\mathrm{B}}$ as unknown functions of the two independent (and presumably experimentally controllable) variables $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$. As the key part of the derivation, we express the coupled set of the equations (15) and (16) also at varied values $V_{\mathrm{A}}+\delta V_{\mathrm{A}}$ and $V_{\mathrm{B}}+\delta V_{\mathrm{B}}$ of the bias voltages. Resulting dependent variables are $\mu_{\alpha}(V_{\mathrm{A}}+\delta V_{\mathrm{A}},V_{\mathrm{B}}+\delta V_{\mathrm{B}})=\mu_{\alpha}(V_{\mathrm{A}},V_{\mathrm{B}})+\delta\mu_{\alpha}$, $\,\alpha\in\{2,3,4\}$, and analogously for the currents. This differentiation leads to a set of algebraic equations of the form (8) and (6), where instead of the variables stand their differentials [13]. As a result, $I_{\mathrm{A}}(V_{\mathrm{A}},V_{\mathrm{B}})$, $I_{\mathrm{B}}(V_{\mathrm{A}},V_{\mathrm{B}})$ and $\mu_{\alpha}(V_{\mathrm{A}},V_{\mathrm{B}})$ are obtained as some generally non-linear functions for a given (fixed) set of $G_{\alpha,\beta}(u)$ spectra, see the Eq. (14).

3.3.4 Non-equilibrium Fermi levels: recursive calculation of the non-linear response

We now return to the Landauer-Büttiker type formula (1). Given bias voltages $(V_{\mathrm{A}},V_{\mathrm{B}})$ on the batteries induce a set of the a priori unknown EChPs $\mu_{1},\dots,\mu_{4}$, in short $\mu$, which enter the formula (1), or, equivalently, the Eqs. (15). The transmission coefficients spectra $\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E})$ depend on the unknown EChPs: $\mathcal{T}_{\alpha\leftarrow\beta}=\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E};\mu)$. In terms of the conductance matrix (14), $G_{\alpha\beta}=G_{\alpha\beta}(u;\mu)$. The dependence on the unknown $\mu$’s makes the analysis of the coupled device much more complicated even in comparison to the non-linear regime treatment in Sect. 3.3.3. The complications were not faced there thanks to the approximation consisting in the use of the the single (fixed) set of spectra $G_{\alpha\beta}(u)$ assumed to be usable for any $\mu$, hence for any $(V_{\mathrm{A}},V_{\mathrm{B}})$ within a studied domain. Obviously, this approximative treatment in general would not yield correct results for the studied model. It is an acceptable approximation for sufficiently low voltages; numerical examples for different levels of approximation will be shown in Sect. 4.2 and in the SM [13]. We note that non-equilibrium response in a weakly non-linear regime was discussed in Ref. [19]. More general approaches to quantum transport including the non-equilibrium EChPs often employ the non-equilibrium Green’s function method, see for instance [20]. In the present study we use the independent-electron approximation as defined in Sect. 2. In order to capture the effects of the non-equilibrium potentials $\mu$ (what we have not done until the present section) we calculate the Green’s functions [18] using the non-equilibrium potential profile. Given the specifics of the system under study (in particular the explicit presence of the circuits and the non-trivial dependence of $\mu$’s on chosen bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$), we proceed in the following steps:

In a quantum-mechanical (QM) calculation we obtain $\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E};\mu)$ with all $\mu$’s set to the equilibrium Fermi energy $E_{\mathrm{F}}$. We denote these coefficients as $\mathcal{T}^{(0)}$. 2. 2.

We accomplish the circuit calculations by solving the set of the Eqs. (15) and (16) using $\mathcal{T}^{(0)}$. In this manner, which was in detail explained in Sect. 3.3.3, we obtain dependences $\mu^{(1)}(V_{\mathrm{A}},V_{\mathrm{B}})$, $I_{\mathrm{A}}^{(1)}(V_{\mathrm{A}},V_{\mathrm{B}})$, and $I_{\mathrm{B}}^{(1)}(V_{\mathrm{A}},V_{\mathrm{B}})$ as a first approximation on chosen domains $V_{\mathrm{A}}\in[0,V_{\mathrm{max}}]$ and $V_{\mathrm{B}}\in[0,V_{\mathrm{max}}]$. Later we refer to the (1)-superscripted results as the iter1 level of approximation. 3. 3.

Assuming that in further search we seek for refined $\mu^{(2)}(V_{\mathrm{A}},V_{\mathrm{B}})$, $I_{\mathrm{A}}^{(2)}(V_{\mathrm{A}},V_{\mathrm{B}})$ and $I_{\mathrm{B}}^{(2)}(V_{\mathrm{A}},V_{\mathrm{B}})$, we take the values $\mu^{(1)}(V_{\mathrm{A}},V_{\mathrm{B}})$ as new reference EChPs. There are as many such references as is the number of the sampled points $(V_{\mathrm{A}},V_{\mathrm{B}})$ within the considered domain. 4. 4.

The new reference (now non-equilibrium) potentials $\mu^{(1)}$ at a chosen particular $(V_{\mathrm{A}},V_{\mathrm{B}})$ point are used for the next QM calculation of $\mathcal{T}_{\alpha\leftarrow\beta}(\mathcal{E};\mu)$, now denoted as $\mathcal{T}^{(1)}$. 5. 5.

$\mathcal{T}^{(1)}$ is used in a new circuit calculation [again solving the Eqs. (15) and (16)] thus yielding more accurate (iter2 ) results $\mu^{(2)}$, $I_{\mathrm{A}}^{(2)}$ and $I_{\mathrm{B}}^{(2)}$, now assumed to be valid at the chosen reference $(V_{\mathrm{A}},V_{\mathrm{B}})$ point only.

Steps (4) and (5) are accomplished for a chosen subset of the grid points $(V_{\mathrm{A}},V_{\mathrm{B}})$ within the considered domain. Next we can continue to improve the results at each chosen point $(V_{\mathrm{A}},V_{\mathrm{B}})$ performing iter3 step, which takes the $\mu^{(2)}$ potentials as a refined reference quite analogously to steps (4) and (5) of the above algorithm. The number of such recursive iterations may be arbitrary. To reach self-consistency of the currents and the EChPs, only few such iterations will be needed as it will be seen in our numerical demonstrations in Sect. 4.2.

The QM calculation of the transmission coefficients requires the knowledge of the electrostatic potential profile in the entire system including the graphene nanoribbon. Within the independent-electron approximation we have to define a model spatial profile of the electrostatic potential in the nanojunction. As specified in Sect. 2, we employ the extended tight-binding hamiltonian (see also Ref. [6]) in which the spatial profile of the potential is expressed in terms of the diagonal matrix elements of the hamiltonian (the on-site energies). Assuming known values of the EChPs at the GNR’s corners (we assume the same values as in the corresponding electrodes), we employ the bilinear-interpolation formula, which defines values of the on-site energies at any $(x,y)$ point of the GNR. In more detail, the on-site energy of an atom $l$ positioned at point $(x,y)$ of the GNR is calculated as $\epsilon_{l}=\epsilon+\Delta\epsilon_{l}$, where $\epsilon$ is the equilibrium on-site energy and $\Delta\epsilon_{l}$ is the bias-induced modification computed using the bilinear interpolation between the four corners. $\Delta\epsilon_{l}$ could in addition include the effect of a gate voltage.

4 Application to graphene nanoribbons

We provide an illustrative set of results for the nanoribbon junction drawn in Fig. 1. Two variants are considered: (i) The perfect GNR and (ii) the GNR with the removed strip of the $60$ atoms marked in the figure by the dark narrow strip, with the detail visualised in the inset of Fig. 3(b).

4.1 Coupled circuits in the linear-response regime at the neutrality point

In this regime very low bias voltages and zero gate voltage are assumed. We present results for both the intact GNR and the perturbed one with the cut-out. Tab. 1, which is placed within this section, lists also data computed for a gate voltage scenario, which will be considered in Sect. 4.3.

4.1.1 The perfect GNR

Using the microscopic approach described in Sect. 2 we obtain numerical values of the linear conductance for the perfect (intact) GNR shown in Fig. 1. The quantum subsystem of the whole scheme (Fig. 2) – the GNR with its four contacts – provides a purely ballistic transport as we assume the low temperature limit. The linear-response electronic transport properties of the junction under the equilibrium conditions are then represented by a matrix $\mathcal{T}$ of the transmission coefficients, or, alternatively, in terms of the corresponding conductance matrix $G=(2e^{2}/h)\,\mathcal{T}$. The matrix is symmetric because of the absence of a magnetic flux. Additional symmetries ($G_{12}=G_{34}$, $G_{13}=G_{24}$, $G_{14}=G_{23}$ and $G_{11}=\dots=G_{44}$) follow from the $D_{2h}$ geometrical symmetry of the junction. Off-diagonal matrix elements for the intact GNR are listed in the first numerical row of Tab. 1 using the units of the conductance quantum $G_{0}=2e^{2}/h$.

The ordering of the matrix elements is given by the indices $1,\dots,4$ at the contacts of the GNR as defined by Figs. 1 and 2. Not shown in the table are the diagonal elements. They take values typically around $53$; the particular intact GNR in the full equilibrium yields $\mathcal{T}_{\alpha\alpha}=53.51$ for all four terminals due to the $D_{2h}$ symmetry. The large diagonal matrix elements can easily be understood for they mainly represent transmissions of the incoming electrons back to the wires belonging to the same electrode. In our model (Fig. 1) we assume $54$ (monoatomically thin) wires per electrode, in other words $54$ conductance channels per electrode. In terms of $\mathcal{T}$, any incoming electron close to the Fermi level is most likely to be transmitted to a wire belonging to the same electrode what is the reason that the cumulative diagonal elements are all close to $54$.

In our first computational analysis of the entire scheme we use resistances $R_{\mathrm{A}}=7.3\,\hbar/e^{2}\approx 30.0\,\mathrm{k}\Omega$ and $R_{\mathrm{B}}=4.9\,\hbar/e^{2}\approx 20.1\,\mathrm{k}\Omega$.

A typical or roughly an average value of the off-diagonal matrix elements for the equilibrium-level conductance (the first row of Tab. 1) is $G^{\mathrm{off}}\approx 0.16\,G_{0}\approx 0.05\,e^{2}/\hbar$. Our choice of the magnitudes of the resistances then corresponds to the regime with the classical conductances $1/R$ of a similar order of magnitude as those of the quantum subsystem. The results are conveniently visualised in terms of currents and are shown in Fig. 3(a). The tiny magnitudes of the bias voltages and consequently of the currents (below $60\,\mathrm{pA}$) ensure the strictly linear regime. As can be seen in the figure, the current through circuit B is weakly influenced by the current in circuit A. In other words, the two circuits influence each other to a minor although noticeable extent. The seemingly simple, almost independent operation, of the two circuits is in fact a results of a cancelation of the internal currents flowing inside the GNR between its terminals. This can be easily understood using the formulae (9) for the currents together with the expressions for Büttiker’s conductances $\alpha_{CC^{\prime}}$ [7]. In terms of current $I_{\mathrm{B}}$, the term responsible for the mutual interaction of the two circuits is proportional to $-\alpha_{\mathrm{BA}}V_{\mathrm{A}}$. For the given highly symmetrical sample the conductance $\alpha_{\mathrm{BA}}=(G_{31}G_{42}-G_{41}G_{32})/S=$ $(G_{31}^{2}-G_{41}^{2})/S\approx-0.00908\,e^{2}/\hbar$ is a very small quantity compared to $\alpha_{\mathrm{AA}}\approx 0.107\,e^{2}/\hbar$ as can be calculated from the matrix elements in the first row of Tab. 1. In other words, because $G_{31}$ and $G_{41}$ are of similar magnitudes, they cause a mutual partial cancellation of the currents flowing along the AC edges and the currents across the diagonal of the perfect GNR.

4.1.2 The perturbed GNR

Electrical conductance properties of GNRs can be affected by different kinds of defects or designed features. For example, the conductance of GNRs with the current flowing along their AC edges can be destroyed by the AC edge disorder [21]. If a defect or a designed feature is such that especially the matrix elements $G_{31}$ and $G_{42}$ (or only one of them) are suppressed, then the cancellation effect discussed above would not work and, consequently, currents in the two circuits would become influencing each other in a more significant extent. For the purpose of a quantitative illustration we consider a perturbed GNR that was formed by the removal of the 60 atoms (15 quartets) from the perfect GNR according to the pattern visualised in Figs. 1 and 3(b) and described by their captions. We note in passing that using the one-orbital-per-atom model we assume that all the edges of the structures are passivated by hydrogens [21]. The linear conductance of the perturbed GNR at the neutrality point is described by the matrix elements, most important of which are shown in the second numerical row of Tab. 1. The matrix element $G_{13}$ is most significantly suppressed in the relative terms. The electric currents resulting from this conductance matrix are plotted in pane (b) of Fig. 3 at otherwise the same conditions as were assumed for pane (a). As can be inspected, pane (b) demonstrates the regime in which the two classical circuits interact through the quantum device more significantly [compared to pane (a)]: currents in circuit B depend on voltage in circuit A more strongly and also currents $I_{\mathrm{A}}$ depend more strongly on voltage $V_{\mathrm{B}}$.

The suppression of the $G_{13}$ matrix element might look intuitively obvious, given the location of the vacant strip in the GNR. Quantitative calculations however show that such a picture is not obvious at all; a removal of, for instance, one quartet of atoms fewer from the perfect GNR (thus leaving the vacant strip correspondingly shorter) would yield a structure with $G_{13}$ about 5 times higher than that for the perfect GNR. Continued removal of more and more 4-atom groups from the GNR (making the vacant strip step by step longer) would yield a sharply oscillating pattern of the conductance characteristics as we have verified. Although such an analysis is not a subject of the present work, we note that similar behaviour is well-know for armchair GNRs which, depending on their width, exhibit three different conductance patterns [1, 22]. We chose the particular perturbed structure with its conductance matrix in order to demonstrate the emergence of the more significant correlations between the currents $I_{\mathrm{A}}$ and $I_{\mathrm{B}}$ which can occur if the above discussed mutual cancellation effect of the internal currents (Sect. 4.1.1) is reduced. A more detailed description and quantitative results for the GNRs with this type of the structural perturbation can be found in the SM [13].

4.1.3 Balancing the Fermi levels

As it can be directly observed from the results of Sect. 3.3.2, the effects of the resistors are significant and most interesting in the regime with the resistances of the magnitudes $R\approx 1/G^{\mathrm{off}}$ (semi-quantitatively).

Employing a tunable resistor (or both of them) in this regime, one can balance the quasi-Fermi levels (3) of a certain pair of the leads, provided that the remaining physical parameters permit it. One particular scenario is described by the Eq. (LABEL:eq:U42_solo). An illustration of the effect for the intact GNR is shown in Fig. 4 in which the currents and the potentials differences are plotted as given by the Eqs. (9), (LABEL:eq:U42_solo) and by the formulae for the other $U_{\alpha\beta}$’s [13]. We keep $R_{\mathrm{A}}=7.3\,\hbar/e^{2}\approx 30.0\,\mathrm{k}\Omega$. The bias voltages on the sources use the values $V_{\mathrm{A}}=2.5\cdot\!10^{-7}\,|t_{\mathrm{B}}|/e\approx 0.74\cdot\!10^{-6}\,$\mathrm{\SIUnitSymbolMicro V} and $V_{\mathrm{B}}=2.6\cdot\!10^{-7}\,|t_{\mathrm{B}}|/e\approx 0.77\cdot\!10^{-6}\,$\mathrm{\SIUnitSymbolMicro V}. The electrostatic potentials of leads 2 and 4 match at the value $R_{\mathrm{B}}^{\mathrm{match}}=7.936\,\hbar/e^{2}\approx 32.6\,\mathrm{k}\Omega$; at this value the solid dark green curve crosses zero thus indicating the equality $U_{2}=U_{4}$. Because of the high symmetry of the studied sample (the intact GNR), the potentials $U_{1}\equiv 0$ and $U_{3}$ (the solid blue plot) are balanced at the same value of $R_{\mathrm{B}}$.

4.2 Non-linear response at non-equilibrium Fermi levels

If the magnitude of one or both of the bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$ is increased above certain threshold

(around $10^{-4}\,|t_{\mathrm{B}}|/e$ in our model), the EChPs $\mu_{2}$, $\mu_{3}$, and $\mu_{4}$ may take values such that the electrons most relevant for the transport acquire energies from the range of the ZZ edge modes.111Particular values of the electrodes’ EChPs at which the ZZ-edge modes are probed depend on an applied gate voltage that can effectively shift the energy levels of the central device. Apart from Sect. 4.3, zero gate voltage (no gate field) is assumed. Conductance matrix $G_{\alpha\beta}$ is significantly energy-dependent in this energy window. Fig. 5 shows the spectra of the matrix elements computed under the full-equilibrium conditions (in particular at zero gate field). The solid green lines represent the matrix elements strongly affected by the presence of the zigzag edges of the GNR.

The sharp oscillations due to the ZZ-edge modes contribute to the currents and make them non-linear with the increasing magnitudes of the the bias voltages $V_{\mathrm{A}}$ and $V_{\mathrm{B}}$. The oscillations start at energies as low as about $8\cdot\!10^{-5}\,|t_{\mathrm{B}}|$ above the equilibrium Fermi level meaning that the transport becomes non-linear even for such relatively low bias voltages. The related theory and computational procedures used to analyse the composite quantum-classical system (Fig. 2) in the non-linear regime assuming a single set of spectra $G_{\alpha\beta}(u)$ were exposed in Sect. 3.3.3. That level of approximation (denoted as iter1 in Sect. 3.3.4) does not in general provide sufficiently accurate results for the model under study. It fails especially for more sensitive quantities like differential conductances. Still, it often provides acceptable estimates for currents $I_{\mathrm{A}}$ and $I_{\mathrm{B}}$ and for the EChPs. A set of iter1 results, in comparison to the linear response regime, can be found in the SM [13].

In Fig. 6(a,d) we present non-linear I/V characteristics, which, in addition, include the effects of the non-equilibrium EChPs using the iterative scheme described in Sect. 3.3.4. The intact (perfect) GNR at zero gate field is assumed and the $I_{\mathrm{A}}(V_{\mathrm{B}})$ and $I_{\mathrm{B}}(V_{\mathrm{B}})$ curves are obtained for the two fixed values of $V_{\mathrm{A}}$. Because the $I_{\mathrm{A,B}}$ vs. $V_{\mathrm{B}}$ plots themselves do not display their non-linear features in a pronounced way, the differential conductances $\partial I_{\mathrm{A}}/\partial V_{\mathrm{B}}$ and $\partial I_{\mathrm{B}}/\partial V_{\mathrm{B}}$ are more convenient for this purpose. We visualise them in graphs (b), (c), (e) and (f) of Fig. 6 using the units $e^{2}/\hbar$ (i.e. not $G_{0}$). Main results in these panes are outlined using the green thick solid plots (iter3 ). The additional curves demonstrate how the recursive calculation converges toward self-consistency of the currents and the EChPs: three iterations of the scheme described in Sect. 4.2 were sufficient to obtain the converged results (the green thick solid plots). This is seen from panes (e) and (f) in which we provide a subset of results at the iter4 level (the discrete black points). They perfectly agree with the iter3 plots.222iter1 differential conductances were obtained using the numerical solution and analytical formulae based on the description in Sect. 3.3.3, with additional details provided in [13]. The derivatives for the higher iterations were obtained by direct numerical differentiation of the currents. The sampling density for the higher iterations was $2.0\cdot 10^{-7}\,|t_{\mathrm{B}}|/e$. Since the solutions were accomplished on the finite-density numerical grid, certain finite-size errors are necessary; we do not employ any smoothening procedure and show the results as they come from the calculations. We have verified that our numerical procedures are stable and a doubling of the grid density would not present any practical difficulty and would yield smoother results. In the SM [13] we provide a direct comparison of several subsequent iterations also for the currents $I_{\mathrm{A}}$ and $I_{\mathrm{B}}$. As a complement we outline there also the EChPs $\mu_{\alpha}(V_{\mathrm{B}})$.

4.3 Resonantly increased conductance in the linear regime

In Sect. 4.1.2 the impact of the particular perturbation of the originally perfect GNR was examined in the low-voltage regime around the neutrality point of the studied model. [Cf. pane (b) vs. (a) of Fig. 3.] The cut-out of the 60 atoms affected especially the $G_{13}$ matrix element, magnitude of which was suppressed to about $0.207$ of the original value (Tab. 1). The effectively almost independent conductance channels [Fig. 3(a)] along the two ZZ edges then turned to the more significantly coupled ones [Fig. 3(b)].

One way how to achieve a regime of the relatively independent currents along the ZZ edges regardless of possible disruption or disorder of the AC edges, is to increase the conductance along the ZZ edges. This can be achieved using resonant transport conditions such that the ZZ edge modes are employed. Such a tuning to the resonance can be obtained using a proper gate field applied to the GNR while still working at tiny bias voltages, i.e. in the linear regime. In this way, upon applying a gate voltage such that it shifts the orbital energies of the GNR by $-0.0016\,|t_{\mathrm{B}}|$, we obtain conductance matrices333The effect of such a gate voltage for our models is largely just to shift the spectra displayed in Fig. 5; either we perform an explicit QM calculation with an applied gate voltage or we just shift the spectra, we obtain practically the same matrices, with a relative difference not larger than $10^{-4}$. with their elements shown in the third and fourth numerical rows of Tab. 1, which should be contrasted to the matrices obtained at the neutrality point (the first and second rows). The $G_{12}$ column of Tab. 1 clearly demonstrates that the junctions under the gate voltage provide greatly enhanced conductance along the ZZ edges. As for the other matrix elements in Tab. 1, variations of their particular values may not be obvious. We have briefly commented on this point at the end of Sect. 4.1.2, with the link to the SM [13]. A graphical representation of the results in terms of the currents, now for the resonantly increased conductances $G_{12}=G_{34}$, may again be convenient: we prepare the two-pane Fig. 7 in the same style as Fig. 3. For easier comparisons we employ the same values of the bias voltages $V_{\mathrm{A}}$, $V_{\mathrm{B}}$ and of the resistances $R_{\mathrm{A}}$, $R_{\mathrm{B}}$ as were assumed for Fig. 3.

The vertical axes are now scaled up to higher values resulting predominantly from the resonantly increased matrix elements $G_{12}=G_{34}$. The increase of the currents in comparison to Fig. 3 is not such as might be expected from the significantly larger matrix elements $G_{12}$. We recall that the classical resistances $R_{\mathrm{A}}$ and $R_{\mathrm{B}}$ in the circuits are assumed to be the same (and significant) as they were for Fig. 3 and they consequently limit the currents. Nevertheless, the relatively large $G_{12}$ numerical values for the junction in the gate potential imply the dominance of the transport along the ZZ edges. The operation of the two circuits is then expected to be mutually almost independent as opposed to the case shown in Fig. 3(b). Graphical representations of the results in both panels of Fig. 7 confirm the expectation: the currents in the circuit B practically do not depend on the voltage applied in the circuit A, at least for the range of the voltages assumed for the graph. Analogously, the dependence of $I_{\mathrm{A}}$’s on $V_{\mathrm{B}}$ is similarly weak for the gated junctions. In this way the ZZ edges can serve as two almost uncoupled conductors for DC low-voltage transport. In the SM [13] we provide more comprehensive results and comparisons including two additional perturbed structures.

5 Conclusions

We considered a coupled quantum-classical device consisting of an atomistic model of a nanoribbon (NR) and of two classical branches of circuits with resistances and bias voltage sources. The four electrodes of the circuits were assumed to be contacted at the corners of the NR so that the edge-induced modes could be used for the electrical transport. The two bias voltages, not electrochemical potentials (EChPs), has been assumed as the control variables. In order to analyse the device operation, we have developed a model, which allowed us to calculate the two electric currents and three non-trivially dependent EChPs in the leads as functions of the two externally controlled bias voltages, with the NR’s conductance matrix (in general energy dependent) and the classical resistances as given parameters of the model. The model is usable for a general four-terminal ballistic device and presents a generalisation of Büttiker’s model [7], which assumes zero resistances in the circuits and is limited to the linear-response regime. As a particular material for our computational analysis we assumed graphene. We calculated the conductance matrix using a microscopic theory based on the Green’s function matrix formalism. Our present implementation of the formalism allowed us to consider graphene nanoribbons (GNRs) up to a size of about $30\,000$ atoms with computational resources of a desktop personal computer.

In the low-voltage linear regime we obtained transparent analytical formulae for the currents $I_{\mathrm{A},\mathrm{B}}$ and electrostatic potentials $U_{\alpha}$ in the system as functions of the applied bias voltages $V_{\mathrm{A},\mathrm{B}}$; hence also formulae for the EChPs $\mu_{\alpha}=-eU_{\alpha}$ of the leads. The formulae for the currents generalise those found by Büttiker [7]. Our formulae show how the interplay of the quantum and classical elements of the entire scheme affects the resulting electric currents. Similarly, our formulae for the electrochemical (or electrostatic) potentials at the device’s terminals display the non-trivial dependences of the potentials on the classical resistances and on the quantum conductance matrix elements.

For the voltages in the non-linear regime we developed an efficient numerical algorithm to compute the $I/V$ and $\mu/V$ curves. The non-equilibrium effects have been included using a recursive algorithm in which self-consistency between the EChPs and the currents is reached within few iterations.

Our analysis has shown under what conditions the two edges of a nanoribbon behave as two relatively independent wires transporting stationary electric currents. On the other hand, the intriguing regime in which the two classical circuits are significantly coupled through the quantum ballistic device has also been studied: an increase of the voltage (and of the current) in one of the circuits induces a substantial change in the current in the other circuit.

Quantitative results have been demonstrated for a GNR with zigzag (ZZ) edge induced high density of states close to the Fermi level. The system was found to behave similarly to two independent conductors especially for resonantly enhanced conductance along the ZZ edges. For other regimes the mutual interaction of the two associated circuits was found noticeable. Results have been shown for the regime $1/R_{\mathrm{A,B}}\lesssim G^{\mathrm{off}}$ ($R_{\mathrm{A,B}}$ - the resistances in the classical subsystem, $G^{\mathrm{off}}$ - a typical value of an off-diagonal conductance element of the quantum subsystem – the nanoribbon). Although we have considered relatively small GNRs with numbers of atoms of the order of $10^{3}$, our model can be applied to any junction supporting ballistic transport of electrons for which its conductance matrix is know. For example, larger-scale devices such as graphene ribbons used in Ref. [4] could be described as well provided that the conductance matrix could be computed or obtained from experimental measurements.

6 Acknowledgements

This work was supported in parts by the Slovak Research and Development Agency under the contract No. APVV-0108-11 and by the Slovak Grant Agency for Science (VEGA) through grant No. 1/0372/13. Useful discussions with Peter Bokes and Richard Hlubina are acknowledged. The GNR structure was visualised using VMD [23].

Contributions of the authors: P. D. did most of the analysis of the semi-infinite chains eigenstates presented in the SM [13], partially contributed to derive the compact forms (9) of the formulae for the currents and derived the relations (38) and (39) in the SM. All other contributions have been done by M. K.