Adiabatic pumping driven by a moving kink in a buckled graphene nanoribbon with implications for a quantum standard for the ampere

A quantum pump in a buckled graphene ribbon with armchair edges is discussed numerically. By solving the Su-Schrieffer-Heeger model and performing the computer simulation of quantum transport we ﬁnd that a kink adiabatically moving along the metallic ribbon results in highly efﬁcient pumping, with a charge per kink transition close to the maximal value determined by the Fermi velocity in graphene. Remarkably, insulating nanoribbon show the quantized value of a charge per kink (2 e ) in a relatively wide range of the system parameters, providing a candidate for the quantum standard ampere. We attribute it to the presence of a localized electronic state, moving together with a kink, whose energy lies within the ribbon energy gap.

As a generic quantum pump transfers electric charge between two reservoirs at zero external bias, solely due to periodic modulation of the device connecting the reservoirs [16], new fundamental and practical aspects of any particular pumping mechanism may be unveiled with the charge quantization at nanoscale. Various single-electron pumps were considered as candidates for the quantum standard ampere [17][18][19]: In case the charge pumped per cycle is perfectly quantized (i.e., equal to Q = ne, with n integer) in a considerably wide range of driving parameters, the output current delivered by the device is I P = ne f P , with f P being the external frequency, and the SI unit of current can be rede- 10 −19 A s, with the second defined via the ground-state hyperfine transition frequency of the cesium 133 atom, ν Cs = 9 192 631 770 Hz [20,21].
So far, single-electron pumps with potential to operate as standard ampere are predominantly based on gate-driven quantum dot systems [19]. We argue here, presenting the results of computer simulation of quantum transport, that the electromechanical pump based on a buckled graphene ribbon, which has recently attracted some attention as a physical realization of the classical φ 4 model and its topological solutions (kinks) connecting two distinct ground states [22,23], may also be considered as a counterpart to the above-mentioned single-electron pumps.
Earlier [15] we have shown that the system similar to the one presented in Fig. 1, consisting of metallic graphene nanoribbon with armchair edges coupled to heavily doped graphene leads, may operate as an efficient quantum pump, but the charge per cycle is not quantized. Here the discussion is supplemented by (i) taking the case of insulating nanoribbon into account, and (ii) by optimizing atomic bond lengths in a framework of the Su-Schrieffer-Heeger (SSH) model [24] including electron-phonon coupling of the Peierls type [25]. As a result, we find that for an insulating ribbon a single electronic state localized at the kink is well separated from extended states, and the charge pumped is quantized. Topological aspects of the system are crucial to understand the charge pump operation, since the electron-phonon coupling leads to the peculiar arrangement of shortened (lengthened) bonds being perpendicular (parallel) to the main ribbon axis in the kink area, resulting in the electron localization.
We also show in this paper that, although the kink shape is well described within the standard molecular dynamics potentials for graphite-based systems (as implemented in the LAMMPS package [26]), for accurate modeling of quantum transport phenomena one needs to include small corrections to the bond length (up to a few percent of the kink area), following from electron-phonon coupling. A minimal quantum-mechanical Hamiltonian, of the SSH type, allowing Buckled graphene ribbon as a quantum pump. Top: Nanoribbon buckled by changing the distance between fixed armchair edges from the equilibrium width of W = 11 a (with a = 0.246 nm the lattice spacing) to W = 0.9W , attached to heavily doped graphene leads (red), each of width W ∞ , separated by distance L 1 . The total ribbon length is L = L 1 + 2W ∞ + 2L s , where L s denotes the distance between the free ribbon edge and the lead edge. The kink is formed near the ribbon center. The ammeter detects the current driven by a moving kink. The gate electrode (not shown) is placed underneath to tune the chemical potential μ 0 in the ribbon area. The schematic potential profile U (x) (bottom left) and the coordinate system (top left) are also shown. Inset: Band structure of the infinite flat ribbon with armchair edges for W = 11 a (solid blue lines) and W = 10 a (dashed red lines). one to model both the kink shape and the transport, is proposed.
The remaining part of the paper is organized as follows. In Sec. II we present the model Hamiltonian and our method of approach. In Sec. III we discuss quantum states of a finite section of buckled ribbon (i.e., closed system) with a kink. The conductance and the adiabatic pumping in the ribbon coupled to the leads (open system) is analyzed numerically in Sec. IV. The conclusions are given in Sec. V.

A. The Hamiltonian
Our analysis starts from the Su-Schrieffer-Heeger (SSH) Hamiltonian for graphene nanostructures [24,27,28], with potential energy describing the covalent bonds [4] where with a constrain i j The kinetic-energy operator for π electrons T [Eq. (2)] includes the hopping-matrix elements (t i j ) corresponding to the nearest neighbors on a honeycomb lattice (denoted by using brackets i j ), with the equilibrium hopping integral t 0 = 2.7 eV. The change in bond length, δd i j = d i j − d 0 , is calculated with respect to the equilibrium bond length, d 0 = a/ √ 3, with a = 2.46 Å being the lattice spacing. The operator c † i,s (or c i,s ) creates (or annihilates) a π electron at the ith lattice site with spin s. The electron-phonon coupling, quantified by the dimensionless parameter β = −∂ ln t i j /∂ ln d i j | d i j =d 0 (to be specified later), is represented by the exponential factor in Eq. (2) replacing standard Peierls form (1 − βδd i j /d 0 ) in order to prevent t i j from changing the sign upon strong lattice deformation.
The next two terms in Eq. (1), V bond (3) and V angles (4), approximates the potential energies for the bond stretching and bond angle bending (respectively); see Fig. 2. The parameters K d = 40.67 eV/Å 2 , K θ = 5.46 eV/rad 2 , and θ 0 = π/3 are taken from Ref. [4] and restore the actual in-plane elastic coefficients of bulk graphene in the case of β = 0. Otherwise (for β = 0), a correction to the potential energy per bond can be estimated as with the number of C-C bonds N b . The second approximate equality in the above is obtained by substituting s c † i,s c j,s + c † j,s c i,s ≈ 1.050 (with i and j the nearest neighbors), being the value for a perfect, bulk graphene sheet at the half-electronic filling [29].
The expression for V angles (4) consists of two terms, each involving summation over the three angles ( j) having a common vertex at a given lattice site j (see Fig. 2). First term, ∝ K θ , represents the harmonic approximation for in-plane bond angle bending. For out-of plane deformations, quantified by the height h j of a tetrahedron formed by jth site and its three nearest neighbors, this term represents a fourth-order correction to the potential energy. A realistic description of out-of-plane deformations requires a correction of the ∼h 2 j order (for h j d 0 ). Here we propose a term proportional to the excess angle, with the coefficient V δ ≈ t 0 = 2.7 eV roughly approximating the bending rigidity of graphene [5,6]. The main advantage of such an approach is that it requires no computationally expensive operations since the four-body term (∝ V δ ) depends only on angles (θ ( j) ) earlier determined for the three-body term (∝ K θ ). The validity of our approach, in comparison with standard molecular dynamics treatments [22,23], is discussed later in this paper.

B. The optimization procedure
Throughout the paper we compare the results obtained in the absence of electron phonon coupling β = 0 and for the dimensionless parameter β = 3; the two values chosen to bound the possible range of β (see Ref. [6]). It is worth to stress here that in the forthcoming analysis physical properties are discussed as functions of the deformation applied, and thus the other choice of β = 0 (being the proportionality coefficient between the local deformation and corresponding correction to the Hamiltonian) may rather shift the characteristic features observed than change the picture in a qualitative manner.
In the β = 0, electronic and lattice degrees of freedom are decoupled, and one simply needs to solve a purely classical minimization problem for the potential energy part of the Hamiltonian H SSH [Eq. (1)], given by V bonds + V angles [see Eqs. (3) and (4)]. For the β = 0 case, the average kinetic energy can be calculated as where the factor 2 in the last two expressions follows from a spin degeneracy t i j = −t 0 exp(−βδd i j /d 0 ) if i and j are the nearest neighbors (otherwise, t i j = 0), and ψ ( j) k denotes the probability amplitude for the kth eigenstate of the kinetic energy operator T [Eq. (2)] at the jth lattice site. We further suppose that the eigenstates are ordered such that the energies E 1 E 2 · · · E N at , with the number of atoms N at , and that the number of electrons N el , is even for simplicity.
In both cases (β = 0 and β = 3), the numerical minimization of the ground-state energy with respect to atomic positions {R j }, is performed employing the modified periodic boundary conditions in y direction (see Fig. 1). Namely, the system is invariant upon y → y + L and z → − z, forcing the kink formation in a buckled ribbon. The outermost two rows of atoms near each armchair edge are fixed during the minimization (see Fig. 2), and buckling of the ribbon is realized by changing the distance between the fixed edges [see Fig. 2(a)] from W to W < W . Furthermore, the number of electrons is fixed at N el = N at , with N at = 3600 for metallic armchair ribbon (W = 10 a, L = 90 √ 3 a) or N at = 3960 for insulating armchair ribbon (W = 11 a, L = 90 √ 3 a). Although in open system, coupled to the leads, the average N el varies with the chemical potential, such fluctuations (typically, limited to μ < 0.1 eV or, equivalently, N el /N at ≈ 0.18 ( μ/t 0 ) 2 < 2 × 10 −4 ; see Ref. [30]) are insignificant when determining the optimal bond lengths. Alternatively, one can interpret the β = 0 case as a hypothetical N el = 0 situation, leading to bond length modifications not exceeding a few percent (see below).
For a fixed W /W ratio and the kink position y = y 0 , the computations proceed as follows.
The initial arrangement of carbon atoms is given by where with (x j ,ỹ j ) being the coordinates of the jth atom on a flat honeycomb lattice, and the scaling function The buckle height H in Eq. (11) is adjusted such that C = N b d 0 in Eq. (5). The kink size is fixed at = 5 a, roughly approximating the kink profiles reported in Refs. [22,23].
At first step, we minimize the potential energy term V bonds + V angles , ignoring a constrain given by Eq. (5). This gives us the solution for β = 0 in the Hamiltonian H SSH (1).
Next step, performed only if β > 0, involves a further adjustment of atomic positions {R j } such that full ground-state energy E G (9) reaches a minimum. In practice, we determine hopping parameters for given {R j }s, and then find (within the gradient descent method) a conditional minimum of E G with respect to {R j } at fixed values of { c † is c js }s, satisfying a constrain given by Eq. (5). The procedure is iterated until the numerical convergence is reached. Typically, after 3-4 iterations the atomic positions {R j }s are determined with the accuracy better then 10 −5 a.

C. Comparison with LAMMPS results
A brief comparison of the kink shape following from the numerical procedure described above with the corresponding output produced by the LAMMPS Molecular Dynamics Simulator [26,31] is presented in Fig. 3.
In order to quantify the difference in atomic arrangements obtained within different approaches, we choose the maximal absolute displacement of an atom along the x (and y) axis, max |x − x (0) | (and max |y − y (0) |), where the maximum is taken for a subset of atoms with equal initial y (0) coordinates; see Figs. 3(a) and 3(b). It is sufficient to display the data corresponding to a vicinity of the kink, 30 y (0) /a 90, since far away from the kink position (being fixed at y 0 ≈ 58.5 a) both quantities considered become y (0) independent. (As free boundary conditions are applied in case the LAMMPS package is utilized, some y (0) dependencies reappear near the free zigzag edges, but they are much smaller in magnitude than dependencies in the kink area).
It is clear from Figs. 3(a) and 3(b) that the LAMMPS results (see red dashed lines) are closer to that obtained with our optimization procedure in the presence of electron-phonon coupling, β = 3 (solid symbols), then for β = 0 (open symbols). Also, x-z views of the system, presented in Figs. 3(c), 3(d), and 3(e), show that approximate mirror symmetry of the kink appears for β = 3 [see Fig. 3(d)] and for the LAMMPS results [ Fig. 3(e)], but is absent for β = 0 [ Fig. 3(c)].
The above observations can be rationalized by taking into account that four-body (dihedral) and long-range Lennard-Jones potential energy terms are included in the LAMMPS package but absent in our model Hamiltonian H SSH (1). In the presence of electron-phonon coupling (β > 0), however, the average kinetic energy T [Eq. (8)] can be interpreted as an effective long-range (and "infinite-body") attractive interaction between atoms, restoring some features related to the Lennard-Jones forces in molecular dynamics (including an approximate mirror symmetry of the kink).

A. Bond-length modulation
Before discussing the electronic structure of the system, we briefly describe small corrections to the bond lengths appearing in the kink area due to electron-phonon coupling (see Fig. 4), which are essential to understand the results presented in the remaining parts of the paper.
In Figs. 4(a) and 4(b) we visualize the spatial arrangements of shortened and lengthened bonds; namely, d i j < d i j i, j=1,...,N at (thick black lines) and d i j > d i j (thin red lines), where the average bond length d i j = 0.998 d 0 for β = 0 and W /W = 0.9, or d i j ≡ d 0 for β = 3 due to a constrain imposed [see Eq. (5)]. Apparently, in the presence of electron-phonon coupling (β = 3), a large rectangular block is formed in the kink area [i.e., for |y − y 0 | 7 a; see Fig. 4(b)], in which almost all bonds oriented in the zigzag direction are shortened (resulting in the hopping element |t i j | > t 0 ) and almost all remaining bonds are lengthened (|t i j | < t 0 ). In the absence of electron-phonon coupling (β = 0) the situation is less clear [ Fig. 4(a)], with a few smaller blocks of shortened or lengthened bonds forming more complex patterns, some of which are isotropic, and some show various crystallographic orientations.
The qualitative finding presented above is further supported with statistical distributions of the relative bond length (d i j /d 0 ), determined using all N b = 5749 bonds in the system (for L = 90 √ 3 a and W = 11 a) and displayed in Figs. 4(c) and 4(d). In particular, the distribution for β = 3 [ Fig. 4(d)] is significantly wider than for β = 0 [ Fig. 4(c)]. Also, bimodal structure of the distribution is visible in the presence

B. The current blocking
In order to understand how the bond-length modulation may affect the transport properties, we focus now on the Dirac points (K and K ) and changes in their positions in the first Brillouin zone due to strain-induced fields (see Fig. 5).
Revisiting the derivation of an effective Dirac equation for graphene, one finds that weak deformations introduce peculiar gauge fields, with the vector potential for K valley [6] where c is a dimensionless coefficient of the order of unity, u i j = 1 2 (∂ i u j + ∂ j u i ) (with i, j = x, y) is the symmetrized strain tensor for in-plane deformations [32], and the coordinate system is chosen as in Fig. 1 (i.e., such that the x axis corresponds to a zigzag direction of a honeycomb lattice). For the K valley, the strain-induced field has an opposite sign (namely, For an approximately uniform compression along the x direction occurring in the kink area, we have u x ≈ (W /W )x, u y = y, and the K point is shift by δA ∝ − (1 − W /W )k x withk x being a unit vector in the k x direction, while the K point is shift by −δA, as visualized in the top panels of Fig. 5. Away from the kink area, buckling without changing bond lengths does not create strain-induced fields (δA ≈ 0) [33].
Additionally, a finite size along the x direction introduces the well-known geometric quantization, with the discrete values of quasimomentum k x , separated by k x ∼ π/W (see bottom panels in Fig. 5). In principle, for a particular combination of δA and k x , a nanoribbon may locally change its character from metallic to insulating (or vice versa). In a more general situation, if δA and k x are not precisely adjusted to alter the system properties at E = 0, one can find some finite energies (E > 0 for electrons or E < 0 for holes), for which quantum states are available only away from the kink area (or only in the kink area). A direct illustration is provided with the density of states discussed next.

C. Density of states
We consider here two nanoribbons with armchair edges, one of width W = 10 a (the metallic case) and the other of W = 11 a (the insulating case). The system length is L = 90 √ 3 a in both cases, with modified periodic boundary conditions (see Sec. II B) applied for both lattice and electronic degrees of freedom. The two values of β = 0 and β = 3 in the Hamiltonian H SSH [Eq. (1)] are considered; the buckling magnitude is fixed at W /W = 0.9. The above parameters allow us to define the two energy scales: The subband splitting and the longitudinal quantization In Fig. 6 we display the electronic density of states with E n denoting the nth eigenvalue of the kinetic-energy operator T given by Eq. (2), for all four combinations of β and W . For plotting purposes, the δ function is smeared by a finite ; namely, we put finite section: in the former case, ρ(E ) is elevated for any E , whereas in the later case, we have ρ(E ) ≈ 0 in a vicinity of E = 0.
The effects of electron-phonon coupling can be summarized as follows. In the metallic case, bond length modulation results in small splittings of the electronic levels [see Fig. 6(b) for β = 3], originally showing approximate degeneracy [see Fig. 6(a) for β = 0], due to amplified scattering between the k y and −k y states occurring in the kink area. In the insulating case, we have two energy levels, appearing for β = 3 [see Fig. 6(d)] but absent for β = 0 [see Fig. 6(c)], one for electrons (marked with red arrow) and one for holes, which occur in the gap range and are well separated from other levels, suggesting that they are associated with localized states.
The above expectation is further supported with local density of states (presented Fig. 7) where the δ function is represented via Eq. (17) and the remaining symbols are the same as in Eq. (8). Adjusting the energy to the isolated electronic level appearing in the insulating case (W = 11 a) at E = 0.04 t 0 , we immediately find that the corresponding quantum state is strongly localized in the kink area (see right panel in Fig. 7). In the metallic case (W = 10 a), the value of E = 0.04 t 0 belongs to a continuum of extended states in the lowest subband, but the corresponding ρ loc (R j , E ) profile shows a clear suppression in the kink area (see left panel in Fig. 7), allowing one to expect that the current propagation in y direction may be blocked, in the presence of a kink, for a whole energy window corresponding to the lowest (or highest) subband for an electron (or holes).

IV. CONDUCTANCE AND ADIABATIC PUMPING IN OPEN SYSTEM
In this section we present the central results of the paper concerning transport properties of the open system (finite section of a nanoribbon attached to the leads) presented in Fig. 1.

A. Simulation details
So far we have discussed several characteristics of the closed system with modified periodic boundary conditions in the y direction (see Sec. III), making the kink position (y 0 ) irrelevant for global characteristics, such as the density of states. Now we use the atomic positions {R j } = {(x j , y j , z j )} obtained with the optimization procedure described in Sec. II B (again, we consider the cases without and with the electron-phonon coupling, β = 0 and β = 3) for y 0 = 3 8 L. Next, the kink is placed at the desired position (say, y 0 + y) by applying a shift to all y coordinates, y → y + y. A series of consecutive shifts, such as visualized in Fig. 8, emulates the kink motion (including full kink and antikink transitions) in a real system. In case the shift is commensurate with the longitudinal ribbon periodicity, y = √ 3na with n integer, we simply apply modified periodic boundary conditions for all atoms, for which y j + y < 0 or y j + y 0. Otherwise (i.e., if y = na), atomic positions after a shift {R j } y are determined via third-order spline interpolation , and x is the floor function of x. The hopping-matrix elements (t i j ) in Eq. (8) are then determined using atomic positions after a shift {R j } y , but we set t i j = 0 in case i and j are terminal atoms from the opposite zigzag edges (i.e., periodic boundary conditions are no longer applied for electronic degrees of freedom).
The leads, positioned at the areas of x < 0 and x > W in Fig. 1, are modeled as perfectly flat (i.e., t i j = −t 0 for the nearest neighbors i and j) and heavily doped graphene areas, with the electrostatic potential energy U ∞ = −0.5 t 0 (compared to U 0 = 0 in the ribbon area, 0 < x < W ), each of the width W ∞ = 17.5 √ 3 a (corresponding to 11 propagating modes for E = 0). What is more, both leads are offset from the free ribbon edges by a distance of L s = 7.5 √ 3 a, suppressing the boundary effects. The scattering problem is solved numerically, for each value of the chemical potential μ = E − U 0 and the kink position y 0 , using the KWANT package [34] in order to determine the scattering matrix which contains the transmission t (t ) and reflection r (r ) amplitudes for charge carriers incident from the left (right) lead.

B. Landauer-Büttiker conductance
The linear-response conductance is determined from the S matrix via the Landauer-Büttiker formula [35,36], namely where G 0 = 2e 2 /h is the conductance quantum and T n is the transmission probability for the nth normal mode. In Fig. 9 we compare the conductance spectra for the same four combinations of parameters W and β as earlier used when discussing the density of states (see Fig. 6). This time, results for a buckled ribbon, with W /W = 0.9 and a kink placed at y 0 = 3 8 L, are compared with the corresponding results for a flat ribbon (solid blue and dashed red lines in Fig. 9, respectively). In the metallic case, electron-phonon coupling strongly suppresses the transport in the presence of a kink [ Fig. 9(b)]; the effect of a kink is much weaker in the absence of electron-phonon coupling [ Fig. 9(a)]. Similar effects can be noticed in the insulating case, provided that the chemical potential is adjusted to the first conductance step above (or below) the gap range [Figs. 9(c) and 9(d)].

C. The pumping spectra
In the absence of a voltage bias between the leads, the charge transferred solely due to adiabatic kink motion (i.e., 165408-7 by varying the parameter y 0 ) can be written as [16] where the summation runs over the modes in a selected (output) lead. We further notice that molecular dynamics simulations of Refs. [22,23] allow us to estimate typical kink velocity (up to the order of magnitude) as v kink ∼ 1 km/s v F , where v F = √ 3 t 0 a/(2h) ≈ 10 6 m/s is the Fermi velocity in graphene, justifying the adiabatic approximation [37].
Numerical results for Q(μ), obtained by shifting the kink from y 0 = 0 to y 0 = L, are presented in Fig. 10. Although the current blocking in the metallic case is far from being perfect [see Fig. 9(b)], the related pumping mechanism for W = 10 a appears to be rather effective (see top panel in Fig. 10), with Q(μ) approaching the total charge available for transfer in a section of the length L eff , a value of which can be approximated by [38] where we put L 1 L eff L 1 + W ∞ estimating the effective length of a ribbon section between the leads (see shaded area in Fig. 10). Significant changes to the Q(μ) spectra are observed in the insulating case of W = 11 a (see bottom panel in Fig. 10). Namely, there is an abrupt switching between Q ≈ 0 near the center of a gap (at μ = 0) and Q ≈ 2e appearing for μ exceeding the energy level localized in the kink area [see Fig. 6(d)]. The value of Q ≈ 2e remains unaffected until μ approaches the bottom of the lowest electronic subband [corresponding to the first conductance step in Fig. 9(d)]. For higher μ, the picture becomes qualitatively similar to this for a metallic case, with Q(μ) systematically growing with μ and decreasing with W /W . Noticeably, the plateau with Q ≈ 2e is well developed starting from moderate bucklings, W /W ≈ 0.95. For W /W ≈ 0.9, deviation from the quantum value in the plateau range is of the order of | Q − 2e| ∼ 10 −4 e, and can be attributed to the finite-size effects. Some stronger deviations may appear in a more realistic situation due to the finite-temperature and nonadiabatic effects, which are beyond the scope of this work.
In both (metallic and insulating) cases, the stability of numerical integration in Eq. (21) substantially improves for the lead offsets L s 5 a (being comparable with the kink size), for which parts of the ribbon attached to the leads, together with a section between the leads, are (almost) uniformly buckled for either y 0 ≈ 0 or y 0 ≈ L.
In Fig. 10 we also display maximal bond distortions for different bucklings (see the inset), showing that local deformations |δd i j | < 0.1 d 0 for all 0.9 W /W < 1.

V. CONCLUSIONS
We have demonstrated, by means of computer simulations of electron transport, that buckled graphene nanoribbon with a topological defect (the kink) moving along the system may operate as an adiabatic quantum pump. The pump characteristic depends on whether the ribbon is metallic or insulating. In the former case, even for moderate bucklings (with relative bond distortions below 10%) the kink strongly suppresses the current flow, and shifts the electric charge when moving between the leads attached to the system sides. In turn, the charge pumped per cycle is not quantized. For insulating ribbon, there are electronic states localized near the kink (with energies lying within the energy gap) which can be utilized to transport a quantized charge of 2e per kink transition (with the factor 2 following from spin degeneracy), providing a candidate for the quantum standard ampere.
Remarkably, the current suppression, and subsequent effects we have described, are visible after the bond-lengths optimization for the Su-Schrieffer-Heeger model is performed, introducing significantly stronger bond distortions than the classical (a molecular-dynamics-like) model optimization. Therefore, electron-phonon coupling appears to be a crucial factor for utilizing the moving kink for adiabatic quantum pumping in buckled graphene ribbons.