Current-induced runaway vibrations in dehydrogenated graphene nanoribbons

Introduction

The electronic and transport properties of graphene has been the focus of intense study since its discovery in 2004 [1]. Due to the strong σ-bonding between carbon atoms, graphene has a very high thermal conductivity, and can potentially sustain much higher current intensities than other materials. Graphene nanoribbons (GNR) exhibit a bandgap opening due to quantum confinement in the transverse ribbon direction. This opens the possibilities of realizing various electronic devices, especially field-effect transistors, using graphene nanoribbons. Atomically precise ribbons [2], as well as more advanced ribbon-based structures [3,4], have been fabricated “bottom-up” on metal surfaces. The conductance through the ribbons has been investigated using STM [5], and the signals of electron vibrations in the current have been addressed by theory [6].

When cutting graphene into one-dimensional ribbons, dangling bonds emerge at the boundary carbon atoms. If there is an electrical current passing through the ribbon, we expect that these boundary atoms with dangling bonds are mechanically weak compared to the central atoms. Actually, it has been observed experimentally that these atoms can be removed by the passing electrical current due to current-induced local heating [7,8]. One theoretical study suggests that the carbon dimers at the armchair edge vibrate locally and interact strongly with the electrical current [8]. They can be thought as atomic scale defects at the boundary. How the current-induced forces affect the dynamics of these dimers is an interesting question to ask since it could be addressed by experiments. Employing a semi-classical Langevin approach, we have previously studied the current-induced atomic dynamics of a graphene nano-constriction [9]. However, the number of atoms involved even in such a small system makes a simple analysis difficult.

On the other hand, since the first prediction that current-induced forces are nonconservative with respect to energy [10], there has been a substantial theoretical effort aimed at exploring its consequences for the stability of current-carrying nano-systems [11-15], or the possibility of driving atomic motors [10,16]. Moreover, it has been shown that, in addition to the nonconservative force, the current-induced forces also include an effective Lorentz force or pseudo-magnetic force, originated from the Berry phase of electrons [11]. Performing similar analysis using a scattering theory approach shows that the predictions apply equally well to much larger mesoscopic coherent conductors [16,17]. A requirement of impact of the nonconservative force is that two or more vibrational modes close in frequency couple to each other via the current carrying electronic states. This can establish a generalized circular “water-wheel” motion, either in real space [10] or in mode space. Another requirement is that these modes have little damping due to the coupling to the phonon reservoir. Unfortunately, there has not been a clear experimental setup where these new theoretical findings can be put to a test proving their effect in an unambiguous way. Thus, it is of interest to be able to propose such a setup based on first principles calculations with realistic unadjustable parameters.

In this paper, we study the current-induced dynamics in a partially dehydrogenated armchair graphene ribbon. We show that, atomic motion of the dehydrogenated carbon dimer at the nanoribbon boundaries are relatively decoupled from other dimers and also from the rest carbon atoms. This results in several nearly degenerate atomic vibrations, where each of these involves mainly one dimer. However, a coupling of the dimer vibrations takes place via the flowing electrical current. All these features are favorable to observe the effect of current-induced forces, thus making armchair nanoribbon an ideal candidate to study.

In the rest of the paper, we summarize our theoretical (section “Theory”) and numerical (section “Numerical Calculation”) methods, and present our analysis of the armchair graphene ribbon (section “Results and Discussion”). We end this paper with our concluding remarks (section “Conclusion”).

Theory

We consider a standard Landauer-type transport setup described in Figure 1a. The system/device is in contact with the left and right leads. Each lead serves as both electronic and phononic bath. The bath degrees of freedom are non-interacting. We are interested in the atomic dynamics in the device region (displacements U), which can be described by the semi-classical generalized Langevin equation (SGLE) [12,18-20],

Here F(U) is the force between atoms in the device region and f(t) is a random force due to thermal and voltage-bias-induced fluctuations. The Π^r (self-energy) describes the time-delayed back action of the bath on the system due to the motion of the system. It has three contributions,

where and describe the coupling to the phonon reservoirs outside device, while the self-energy describes the coupling to the electrons. In equilibrium, this result has been obtained by, i.e., Head-Gordon and Tully [21]. In principle we may apply the SGLE including the non-linear part of F [22], but here we will restrict ourselves to the harmonic approximation, F(U) = −KU.

Non-equilibrium

The semi-classical Langevin equation can be extended to include the non-equilibrium effects in the electronic system due to the current [11,12,23]. In accordance with intuition the “traditional” Joule-heating is present in the fluctuating force, f, while the current-induced forces show up in . Here we focus on the latter. The contribution from the electron degrees of freedom including the non-equilibrium effects can be expressed in terms of the coupling-weighted electron–hole pair density of states, Λ,

with Λ (including spin), given by,

Here, A_α/β is the density of scattering states incoming from left and right electrodes (indices α and β), while M describes the electron–phonon couplings (k and l phonon indices). One can loosely think of the motion of phonon k excites an electron–hole pair of energy which is absorbed by phonon l.

The SGLE, in Equation 1, is given in the time domain. However, since we are considering steady state, it is convenient to work in the frequency domain. Thus, by Fourier transformation we obtain,

By applying the Sokhatsky–Weierstrass theorem Π^r(ω) can be split into four contributions giving rise to the four forces

Here, FR, NC, RN, BP represent the electronic friction, nonconservative force, renormalization of the atomic potential, and Berry-phase-induced pseudo-magnetic force, respectively [12].

Run-away modes

In order to analyze the influence of the current we define the nonequilibrium phonon density of states (DOS) as,

where D^r(ω) is the nonequilibrium phonon Greens function obtained from the SCLE,

Note that we introduced boldface to underline that these are matrices with mode-index (k,l). Contrary to the equilibrium situation, the DOS given in Equation 8 can take negative values at certain peak values, due to the electronic current. We can interpret a negative peak in the DOS at a frequency ω₀ as modes at ω₀ with a negative lifetime, i.e., with growing in amplitude as a function of time and denote these by ”run-away” modes.

Mode analysis

In order to identify the modes that can show run-away behavior, we need to find the solutions to Equation 1 by setting the driving noise force, f(ω), to zero. This is done by treating the velocity and displacement as independent variables and use the relation to obtain the double-sized eigenvalue problem. However, the self-energy Π^r is frequency-dependent. Thus, to analyze a specific runaway mode giving rise to a negative peak in Figure 2a (see below), we evaluate the self-energy at the negative peak frequency ω₀ as given below in Equation 10,

Thus, the dynamical matrix is renormalized by and the friction originates from . Solving Equation 10 gives a set of eigenmodes and complex eigenfrequencies, but only the “self-consistent” mode, which fulfills Re(ω) = ω₀, is relevant.

For a given eigenmode a corresponding positive imaginary part of the eigenfrequency designates that the mode is a run-away mode, while if the imaginary part is negative the mode is damped. The damping can be quantified by the inverse Q-factor giving the change in energy per period

Thus, the run-away modes can be identified as the modes where Im(ω) > 0. The run-away modes are a linear combination of the non-perturbed normal modes. Normally, the runaway makes closed loops in real or in abstract mode space. Thus, the NC force allows the mode to pick up energy every time a loop is completed, eventually leading to break down of the harmonic approximation, ending with, e.g., rupture or damping by anharmonic effects leading to a limit cycle motion [24].

Numerical Calculation

We have calculated the electronic and phononic structure of the graphene nanoribbon from density function theory (DFT) using the SIESTA/TranSIESTA codes [25,26]. The generalized gradient approximation is used for the exchange–correlation functional, and a single-ζ polarized basis set is used for the carbon and hydrogen atoms. A cut-off energy of 400 Ry is used for the real-space grid. The electron-vibrational coupling is calculated using the INELASTICA package, which uses a finite difference scheme [27].

Results and Discussion

The partially dehydrogenated graphene nanoribbon we considered is shown in Figure 1a, where four hydrogen atoms have been removed on each side of the ribbon. In principle, dehydrogenation could be performed at chosen positions with an STM [28]. The same structure has been considered in our recent work, focusing on the asymmetry in phonon emission and heat distribution due to the nonequilibrium for lower voltage than where run-away instability occurs [29]. This asymmetry is intrinsically linked to momentum transfer from electrons to phonons and thus to the nonconservative current-induced forces (“electron wind”) [30]. As we mentioned before, the reason we choose this structure is that the dehydrogenated carbon dimers can be considered as “defects” in the armchair ribbon. In general defects give rise to modes localized around the defect. They originate from the local change in force constants shifting the mode out of its unperturbed subband [31].

The relative motion of the two carbon atoms in each dimer only couples weakly to the motion of neighbouring dimers and to the phonons in the leads. This high-energy mode is shifted out of the entire phonon bandstructure. Meanwhile, the flowing electrical current passing through these dimers introduce a small bias-dependent coupling via the self-energy (Π_e) in Equation 10. In principle, one can tune the relative distance between different dimers by changing the ribbon width or the position of the dehydrogenation. It is an ideal and clean system to study the current-induced atomic dynamics, with some tunability since one may imagine doping or gating to shift the Fermi level, E_F, as well as changes in geometry such as varying the distance between dimers.

In this study, the nanoribbon has a width of 7 dimers corresponding to a C–C edge distance of 7.5 Å. The lateral confinement introduces a direct semi-conducting band gap, giving rise to the gap in the electronic transmission for the perfect ribbon as shown in Figure 1b. We have introduced a broadening of the electronic states as in Christensen et al. [6] to mimic coupling to metallic electrodes. This in effect smoothens the transmission curves Figure 1b (dotted lines) akin to the experimental conductance [5]. The introduction of the defects results mainly in a potential shift, but besides this does not impact the transmission dramatically, as seen in Figure 1b (solid lines).

In order to characterize the phonons we show the phononic transmission in Figure 1c. The phonon transmission shows a significant reduction of approx. 50% for above 25 meV. The fact that we see little difference between a perfect and a defect system at low phonon energy is expected when the wavelength is much larger than the defect. We note that, ideally, the translation (T) in the three spatial directions and rotation around the longitudinal direction of the GNR should lead to perfect T = 4 at , equal exactly to zero for both pristine and defected structure. The deviation is due to our numerical neglect of long-range elastic forces. However, while the low energy/long wavelength modes are important for heat transport, we are here concerned with modes with a higher frequency above 25 meV, where the calculation is expected to be accurate.

The influence of the current-induced forces on the phonon self-energy depend on the underlying electronic properties. Thus, besides the unperturbed phonon DOS exclusive the current-induced forces (black solid line in Figure 2a), we also calculate the non-equilibrium DOS using Equation 8 shifting the Fermi energy away from the gap to E_F = −1.4 eV (green dotted line) and E_F = 1.4 eV (red dashed line), and applying a bias of V_b = 0.5 V. Comparing these results, we see that there are several run-away modes (negative peaks) for E_F = 1.4 eV, but not for E_F = −1.4 eV. In Figure 2b and Figure 2e we show two typical run-away modes marked as 1 and 2 in Figure 2a. In Figure 2c,d and Figure 2f–h, we also plot the “bare” normal modes (without coupling to electron and phonon baths), which give the largest contribution to the two selected run-away modes. These run-away modes have in common that they are spatially localized, meaning a vanishing damping due to the coupling to the phonon reservoirs, which is a prerequisite for the run-away instability here. The driving of the motion by the current has to exceed this damping. Both run-away modes involve mainly the dehydrogenated carbon dimers, but while mode 2 lies outside the bulk phonon bands and is thus localized, the frequency of mode 1 is well within the entire phonon band, illustrating that this is not a necessary requirement. In the case of 1 the localization is due to a shift out of a ribbon phonon sub-band. Most of the “bare” normal modes (Figure 2c,d and Figure 2f–h contributing to run-away motion can be considered as a in or out of phase combination of different dimer vibrations.

The two selected modes illustrate how the run-away modes can display circular motion in real-space or in abstract mode space. For the modes showing circular motion it is intuitively clear why these have been dubbed “water-wheel” modes [10], and that the direction of rotation is linked to the direction of current via (angular) momentum transfer [30]. Mode 1 is made up by two principal bare modes, while mode 2 does so in abstract mode space, and consists of mainly three bare modes. In both cases the nonconservative force pump energy into these modes when they oscillate around closed loops. We note that anharmonic coupling will lead to additional damping, which is not included in the harmonic approximation applied here, but we expect a lower anharmonic coupling to the mode well outside the bulk band (mode 2) due to the frequency mismatch. The structure we consider has mirror symmetry in the lateral direction perpendicular to the transport. The resulting motion of the run-away mode respects this symmetry. But, the motion along the current direction is asymmetric since the current breaks the symmetry. This is more obvious for mode 1 (Figure 2c). Mode 2 is strictly localized in the center, but mode 1 has weak coupling to the leads. Consequently, mode 1 shows larger asymmetry. We should mention that this asymmetry in the local heating already shows up before the run-away modes emerge [29].

Conclusion

In conclusion, we have studied the effect of current-induced forces on the dynamics of dehydrogenated carbon dimers at the edges of graphene armchair ribbon. These carbon dimers are weakly coupled to each other and the rest of the carbon atoms, but they interact with the electrical current. This induces effective coupling between them, and the nonconservative and effective magnetic force become important in describing their dynamics. Using a simplified eigenmode analysis, we analyze how the carbon dimer motion is modified by these forces at different Fermi level positions. The possibility of observing the atomic structure of the two-dimensional structures in microscopy, gating or doping, and atomic scale modification of graphene ribbon boundaries makes it an ideal candidate to study current-induced forces in nanoconductors, where interesting theoretical predictions are awaiting for experimental confirmation.

So far current-induced motion and desorption have been observed around edges in graphene sheets [7,8]. One signature of the nonconservative forces is, besides the asymmetry build into the momentum transfer, the highly non-linear heating of modes with bias [11], which in principle could be observed around edges [32].

Acknowledgements

We acknowledge computer resources from the DCSC, and support from Center for Nano-structured Graphene (Project DNRF58). J.T.L. acknowledges support from the National Natural Science Foundation of China (Grants No. 11304107 and No. 61371015).