Nonconservative current-driven dynamics: beyond the nanoscale

Introduction

The development of electronic devices at the nanoscale is a challenging avenue of research with the aim of improving their efficiency and performance. This requires an understanding of the mechanisms for energy transfer from current carriers into atomic motion. Large current densities can generate significant additional forces on atomic nuclei [1-3], resulting in a class of phenomena known as electromigration: atomic rearrangements and mass transport driven by current flow [4,5]. Recent work has drawn attention to another aspect of these forces, anticipated in a visionary argument by Sorbello [5]: unlike equilibrium interatomic forces, they are nonconservative (NC), enabling the current to do work on atoms around closed paths [6-9].

This mechanism for energy conversion from current into atomic motion, which we refer to as the waterwheel effect, differs from Joule heating [10,11] in two key respects. First, the growth in atomic kinetic energy is exponential. Second, it is not stochastic: the energy transferred in the waterwheel effect is stored in directional motion, specifically as generalised angular momentum [12].

In the early work above, it seemed that the waterwheel effect might require rather specialised conditions. The effect operates fundamentally through the coupling of pairs of normal modes to form generalised rotors driven by the current. This requires modes that are close in frequency and are, furthermore, strongly coupled by the NC current-induced forces.

A class of systems where these requirements are met are long, low-dimensional metallic wires [13]. They have a dense frequency spectrum providing the desired degeneracies. In addition, frequency renormalisation by the current (which in general can ruin the degeneracies) is small in these quasi-ballistic systems. Finally, electrons couple strongly to extended phonon modes with the wavevectors needed for momentum conservation under backscattering [14]. Simulations under current indeed show NC dynamics in long atomic wires on a grand scale [13].

This study revisits the waterwheel effect in long wires, and reports on two further aspects of this problem. The first is the physical interpretation of the effect. Originally the effect was demonstrated for a system with just two degrees of freedom – the corner atom in an atomic wire with a bend [6]. Under the right conditions, current drives the atom around an expanding orbit in analogy with a real waterwheel, enabling an intuitive picture of how NC forces work. We will see that an intuitive analogy at the other end of the spectrum, that is, in extended systems, is also possible: it is how strong winds generate forward-travelling ripples on a lake, or the uncompensated stimulated emission of directional phonons [12,15].

However, this process hinges on momentum conservation, and for waves, this information requires a sufficiently long-ranged physical property. For atomic motion under current, this property is the nonequilibrium dynamical response matrix, whose antisymmetric part (induced by the current) describes the NC forces [8,16]. The second advance reported here is the quantitative analysis of this property of long metallic nanoconductors. We show that in long nanowires this antisymmetric part becomes very long-ranged. This nonlocality is essential for the characterisation of NC dynamics under current beyond the nanoscale.

Methods

The system investigated is illustrated in Figure 1. A central region, C, whose middle section containing 200 atoms will be treated dynamically, and two electrodes, L and R, are used to generate current flow.

We are interested in the behaviour of the ions (nuclei and core electrons) in these systems under current, and in particular, the gain in atomic kinetic energy due to work done by the current. The force we are considering is the mean force exerted by electrons on ions. It is determined by the rate of change of the expectation value of the ionic momentum. Using the Ehrenfest approximation, the force on an atomic degree of freedom n is given by

where is the one-electron density matrix and is the one-electron Hamiltonian as a parametric function of the atomic position .

We will employ two different approaches to determine : an adiabatic steady-state approach, where = (V,) is a function of the bias V and the geometry , and a nonadiabatic dynamical approach (within the mixed quantum-classical Ehrenfest method), where = (t) is obtained from an open-boundary quantum Liouville equation [17]. These approaches will be discussed in more detail later.

Electrons are described within a spin-degenerate single-orbital orthogonal nearest-neighbour tight-binding model with parameters fitted to the elastic properties of bulk gold [18]. The nearest-neighbour Hamiltonian matrix elements have the form

where R_mn is interatomic distance, ε = 0.007868 eV, a = 4.08 eV, c = 139.07, and q = 4. In addition, the model includes a repulsive pair potential of the form

with p = 11. The onsite elements of the Hamiltonian are set to zero, and the electron band filling is ν = 0.36361. Noninteracting electrons are considered throughout. As in [13], we compress the chain to a lattice spacing of R = 2.373 Å to suppress a Peierls distortion and resultant band gap that form after geometry relaxation.

Landauer steady state

In the adiabatic steady-state method for the electronic structure, we employ the Landauer picture of conduction. Here the electrodes are infinite, and electrons populate sets of stationary Lippmann–Schwinger scattering states, arriving from either side. The respective population functions, f_L(E) and f_R(E), correspond to the electrochemical potentials of the left and right source reservoirs. The steady-state electron density matrix is then given by [19]

where (E), i = L,R are the density of state operators (subsuming spin degeneracy) for the two sets of scattering states. We work at zero electronic temperature where the occupations are step functions. The density of states operators are generated by Green’s function techniques.

We can now calculate the forces on ions about a chosen reference geometry, . Under small displacements, d,

where K_nm = −∂ F_n()/∂ R_m are the elements of the dynamical response matrix.

The dynamical response matrix determines the vibrational frequencies and corresponding collective modes of motion of the ions. We ignore velocity-dependent forces in the present steady-state description (although they will be present in the nonadiabatic dynamical simulations later). These forces can be included perturbatively [8,20] and tend to dampen the atomic motion, and introduce a contribution arising from the Berry phase [8]. Force noise is also excluded here.

The vector containing the atomic displacements can then be expressed as

where {} are the eigenvectors of the dynamical response matrix, with frequencies {ω_j}, and the static contribution {} is determined by any residual forces present in the chosen reference geometry. {A_j} and {B_j} are set by initial conditions.

The dynamical response matrix can be separated into an equilibrium and a current-induced part:

The current-induced part, in turn, is separated into a symmetric and antisymmetric part [16], ΔK = S + A, with

The antisymmetric part, present only under bias, is a generalisation of the curl of the force on an ion [13,16]. The resultant nonhermiticity of the dynamical response matrix under current in general generates complex frequencies. The complex modes come in complex conjugate pairs. Via Equation 6 these modes give rise to solutions that grow or decay exponentially in time. These are the waterwheel mode pairs investigated in [6,12,13,16].

Within this steady-state approach, current is determined from the bias and the reference geometry, and is not allowed to respond to the subsequent motion of the ions. This approach is accurate for small atomic displacements and large atomic mass (suppressing the velocity-dependent forces relative to the nonconservative forces, and also the work rate due to inelastic scattering [21]), and in systems where fluctuations in the current due to deviations from ideal steady-state behaviour are not too large [13].

Dynamical simulations

To simulate departures from the above ideal conditions, we use the nonequilibrium nonadiabatic molecular dynamics method of [17]. Now the leads are finite and embedded in external electron bath supplying carriers. The leads in the present simulations will be 250 atoms long, with open-boundary parameters, Γ = 0.5 eV and Δ = 0.0005 eV [17]. The electron density matrix then evolves according to the open-boundary equation of motion with the source and sink terms present [17] in the presence of the atomic motion. Atoms obey Newtonian equations, with forces found from the time-evolving density matrix via Equation 1. This form of electron-ion dynamics is known as Ehrenfest dynamics. It captures all forces, equilibrium and nonequilibrium, with the exception of the force noise associated with spontaneous phonon emission and Joule heating [20]. By contrast with the Landauer method above, the dynamical simulations accommodate departures from steady-state conditions and allow the current to respond to changes in the vibrational amplitudes.

To compare the two methods we will further perform a short-time Fourier transform on the ion trajectories (t) from the dynamical simulations and examine the evolution of the energy distribution across the phonon band. The Fourier transform uses a Blackman window, effectively suppressing data outside a particular time interval, while ensuring the data remains continuous. The frequency spectrum of the total ionic kinetic energy (for all N ions) is then

where {A_n(ω)} is a Fourier component of {R_n(t)}. The window is then moved along in time.

Results and Discussion

Landauer steady-state calculations

First we analyse the equilibrium vibrational modes determined from the dynamical response matrix for a long metallic wire. The mode analysis is then performed under bias, where complex frequencies are present. Long-range interactions in the dynamical response matrix are investigated.

Equilibrium mode analysis

We analyse the equilibrium eigenfrequencies and eigenmodes for longitudinal phonons in a chain with 200 mobile atoms of mass M = 10 amu. The 200 eigenvalues {k_j} of the 200 × 200 dynamical response matrix give rise to 2 × 200 eigenfrequencies: positive and negative square roots of the eigenvalues divided by mass, ω_j = . Each eigenvector of the dynamical response matrix is of length 200. Its elements give the relative amplitudes of the atoms in the given mode. These real-space eigenmodes are normalised to unity and Fourier transformed into momentum space. Figure 2 presents the mode frequencies (vertical axis), together with the modulus (colour) of the Fourier components (k, horizontal axis) of the corresponding eigenvector.

Notice the dip around k = ±0.73π/R. It arises due to the long range behaviour of the dynamical response matrix. To check this, we have determined the dispersion relation by truncating the dynamical response matrix beyond a chosen cut-off range. The dip appears when the range includes at least third or fourth neighbours. The overall shape of the curve in Figure 2 is also sensitive to the truncation range, with the shape in the figure emerging at about 20 lattice spacings. The dip in Figure 2 is qualitatively similar to experiment [22], and occurs at a wavevector of about ±0.73π/R. This is twice the Fermi wavevector κ_F = νπ/R. We conclude that this dip is the result of a Kohn anomaly [23].

Mode analysis under bias

The key difference between the equilibrium and nonequilibrium dynamical response matrices is the antisymmetric part of the latter. For an infinite perfect chain it can be evaluated analytically:

where β is the hopping integral, = cos⁻¹(μ_L(R)/2β), and f() denotes the whole preceding expression in the square brackets with replaced by . The upper panel in Figure 3 shows the relative values of the antisymmetric contribution as a function of site separation and bias.

We see that A_mn is oscillatory and long-ranged and, at small bias, becomes infinitely long-ranged. The lower panel in Figure 3 shows the number of waterwheel pairs formed under bias for different truncations of the dynamical response matrix. We see (main panel) that the data deviates from the plateau by a chosen fractional amount, for longer truncations at low bias; consistent with the upper panel. The inset then shows that the sensitivity to truncation is set by the nonequilibrium, antisymmetric part of the dynamical response matrix.

We now turn to chains relaxed at zero bias as the reference geometry (which remain close to the perfect chain). Figure 4 presents the eigenvector of the dynamical response matrix, under a bias of 0.5 V, for the waterwheel mode with the largest negative imaginary part to its eigenfrequency, ω = 0.237 − 0.083 fs⁻¹ i, for a wire with 200 mobile atoms (mass 10 amu). The inset shows the real and imaginary parts of the eigenvector (now complex), together with its modulus; the main panel displays the moduli of the Fourier components of the mode in k-space. (The Fourier picture is qualitatively similar for other waterwheel modes.)

First, what motion does this mode describe? Given ω = ω’ + ω” i as the frequency (with ω’ and ω” real) and v_n the mode component at site n, the corresponding physical displacement is of the form

with an amplitude of C. The Fourier spectrum shows that the mode is dominated by a negative k-component (close in magnitude to the value 2κ_F required for momentum conservation in electron–phonon interactions, where the Kohn anomaly occurs). Since ω’ > 0 and ω” < 0, we obtain a right-travelling wave that grows in time. The small contribution in the Fourier spectrum at the corresponding positive k ≈ 2κ_F describes phonon waves travelling the other way, which are attenuated by the electron particle current (flowing to the right).

These observations can be summarised by the physical picture in the Introduction: NC forces in long metallic systems generate modes of motion, in which the current-carrying electrons close to the Fermi level emit a directed shower of forward-travelling phonons, in analogy with how a breeze creates waves on a lake.

In addition to the features at ±2κ_F, Figure 4 shows a weak background of other Fourier components. By mixing in these other wavevectors, the mode redirects some of the energy gained from the electron “wind” to phonon momenta that can no longer interact directly with the Fermi electrons, and be reabsorbed. This enables the mode energy to grow in time.

We can use Equation 6 to simulate the atomic motion. We calculate the forces under bias for the zero-bias relaxed geometry and use them with zero initial displacements and velocities to set the coefficients {A_j} and {B_j}. We do not include the friction forces here, but we cut the imaginary parts of the frequencies by a factor of 5 to stretch out the growth of the amplitudes in time. Figure 5 shows the displacements of the ions as a function of position and time for a bias of 0.2 V.

The right-travelling phonons generated by the current are evident. The velocity of the peaks and troughs is about 2.5 × 10⁴ ms⁻¹, which is approximately equal to ω/2κ_F, with ω ≈ 0.25 fs⁻¹. This representative frequency is close to: the Einstein frequency for this system; the typical real part of the (closely-clustered) waterwheel modes; the frequency where the dispersion relation (Figure 2) starts to flatten out; and as we will see, the dominant frequency observed in the dynamical simulations (below).

Nonadiabatic nonequilibrium dynamical simulations

The dynamical simulations under bias are performed for atoms initially at rest in the zero-bias relaxed geometry. The bias is then ramped up at the start. The NC forces cause the ionic kinetic energies to rapidly increase. Unlike the steady-state analysis above, the current now responds and is suppressed by the atomic motion. The system eventually settles at a dynamical steady state, where the velocity-dependent friction forces balance out (on average) the driving NC forces. This interpretation is supported by the fact that balancing these forces leads to analytical predictions that agree with the nonadiabatic simulations [13]. A further, independent verification of this balance is given below. In the simulation below the applied bias is 0.5 V; in the long-time dynamical steady state the current settles at a value corresponding to an effective reduced bias of about 0.2 V (the bias used in the adiabatic visualisation in Figure 5).

Figure 6 presents the temporal Fourier composition of the ionic kinetic energy, as in Equation 10, for the first 6 ps of the simulation. Initially, the growing energy is stored in a narrow frequency range, clustered around the representative frequency above. As the amplitudes increase, phonon–phonon scattering eventually redistributes the energy across the phonon band, resulting in approximate energy equipartioning among available frequencies [13].

The dynamical simulations can be used to give a measure of the efficiency [24] with which the NC forces convert electrical energy into atomic motion. From [13] we can estimate the average imaginary part, Φ, of mode frequencies at a given current. E₀ is the total atomic kinetic energy in the dynamical steady state (beyond about 2 ps in Figure 6). The rate of work by the NC forces is then W_NC ≈ 2E₀2Φ (a factor of 2 to give total, as opposed to just kinetic, vibrational energy, and another to convert amplitude to intensity). For the simulation in Figure 6, the current in the dynamical steady state is I ≈ 14.2 μA and E₀ ≈ 17.5 eV, giving Φ ≈ 10⁻³ fs⁻¹ and W_NC ≈ 0.07 eV fs⁻¹. This can be compared against the power, W = IV, due to the transfer of electrons between reservoirs, W ≈ 0.04 eV fs⁻¹. Thus, for the present systems, these two quantities are comparable. A more detailed investigation of this comparison, including its system-dependence, is clearly an important direction for further work.

We can also use the above estimates to independently verify the balance between friction and NC forces in the dynamical steady state. For estimation purposes, we use the analytical result for the friction coefficient, and corresponding energy relaxation time, τ_frict, for atomic Einstein oscillators [17]: 1/τ_frict = (4/Mπ)(H’/H)², for a nearest-neighbour tight-binding chain, where H and H’ are the hopping integral and its derivative with distance. For the present parameters, and for the above steady-state kinetic energy, the power lost to friction is W_frict = 2E₀/τ_frict ≈ 0.08 eV fs⁻¹, in agreement with the (independent) estimate of W_NC above.

Conclusion

Long low-dimensional metallic systems are a promising testbed for NC current-driven atomic dynamics. We have highlighted two aspects of these effects here: the physical interpretation of NC motion as “ripples” driven by the electron “wind”, and the long-ranged character of the nonequilibrium parts of the dynamical response matrix, responsible for NC dynamics. The inclusion of Joule heating (suppressed in Ehrenfest dynamics) and its interplay with the NC forces is an attractive direction for further work, as is the current-driven dynamical behaviour in the presence of the Peierls instability that occurs under compression-free conditions. We hope that this work will motivate further research into some of these questions.