Non-linear conductivity of metals from real-time quantum simulations

Non-linear conductivity of metals from real-time quantum simulations

Xavier Andrade    Sébastien Hamel    Alfredo A. Correa Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94550, USA

We simulate bulk materials under strong currents by following in real-time the dynamics of the electrons under an electric field. By changing the intensity of the electric field, our method can model, for the first time, non-linear effects in the conductivity from first principles. To illustrate our approach, we show calculations that predict that liquid aluminum exhibits negative-differential conductivity for current densities of the order of . We find that the change in the non-linear conductivity emerges from a competition between the accumulation of charge around the nuclei that increases the scattering of the conduction electrons, and a decreasing scattering cross-section at high currents.

In 1827, Georg Ohm established the law that carries his name by discovering that when a voltage is applied over a piece of a certain material, it will carry a current proportional to that voltage. The proportionality constant depends on the geometry of the conductor, and a property of the material: the conductivity. Since then, the conductivity has become one of the most basic quantities that characterize a material. In fact, we usually classify solids into three mayor groups: metals, semiconductors and insulators, based on their ability to conduct electricity.

As it is common with 19th century “laws”, Ohm’s law is more a first order approximation to a more complex behavior than a fundamental property: many materials do not follow Ohm’s law, diodes for example, and for those that do, we can expect that for large enough voltages the conductivity of a material will change. Probably the simplest case of non-linear conduction is the dielectric breakdown, which occurs when an insulator becomes conductive if it is exposed to a very strong electric field. In some materials it can even happen that the current decreases as the voltage is increased, leading to the phenomenon of negative differential conductivity (also known as negative differential resistance) Ridley (1963); Volkov and Kogan (1969); Pamplin (1970). Negative differential conductivity has been observed mainly in semiconductors, like GaAs or GaN, and is used for technological applications Gunn (1963); Gružinskis et al. (1999). It has also been reported in molecular electronic devices Chen (1999); Xue et al. (1999); Dalgleish and Kirczenow (2006); Perrin et al. (2014) and nano-structures Lyo and Avouris (1989); Rinkiö et al. (2010); Zheng et al. (2010); Wu et al. (2012); Du et al. (2012); Lin et al. (2015). As far as we know, negative differential conductivity has not been observed in metals due to the high fields that would be required Pamplin (1970).

For some materials the conductivity is difficult to measure, for example in the case of matter under extreme conditions that are not easy to replicate in an experimental setup. In situations like this, the computational study and prediction of transport properties is an essential tool to understand experiments, and to help in the development of technological applications like thermoelectric power generation Chen et al. (2003) or molecular electronics Lortscher (2013).

The theoretical determination and prediction of linear transport coefficients in materials has been studied from many different points of view Drude (1900); Ziman (1961); Ashcroft and Mermin (1976); Mahan (1990). In most cases, the electrons are the ones that carry electricity, so the conductivity is determined by the electronic structure of material. This means that an accurate theoretical description of conduction requires a quantum mechanical treatment of the electrons. From a first-principles perspective, the electrical conductivity of crystalline systems is usually calculated using perturbation theory, by combining density-functional theory and the Kubo-Greenwood formula Kubo (1957); Greenwood (1958). This method only yields the conductivity in the low current regime, so up to now non-linear effects in the conductivity have been out of reach for first principles methods.

In this work, we present a new method to obtain the conductivity of bulk materials that not only allows us to access the standard conductivity, even in the limit of zero frequency, but also the effective non-linear conductivities for high-intensity currents. To do it, we apply an electric field, and we simulate the response of the material by following the quantum dynamics of the electrons in real-time. By controlling the intensity of the applied field we gain access to the non-linear regime. We can also monitor the time- and space-resolved electronic and current densities to try to understand the microscopic nature of conduction and how non-linear effects appear.

We apply our method to the calculation of the conductivity of aluminum at high temperatures and under strong electric fields. Under these conditions, we find that the conductivity depends on the intensity of the current, leading to negative differential conductivity. To validate our results, we develop a simple model based on ion-in-jellium scattering theory that shows a similar dependence of the conductivity with the current density. We also find that the inhomogeneous accumulation of electronic charges in the system shows distinct patterns for the different conduction regimes, allowing us to gain insight into the cause of the changes in conductivity.

I Theory

To describe the dynamics of the electrons we use an effective Schrödinger equation in the form of time-dependent density functional theory (TDDFT) Runge and Gross (1984) by solving the time-dependent Kohn-Sham (KS) equation (Gaussian atomic units, where , are used throughout)


where are the time-dependent KS orbitals, is the KS single-particle Hamiltonian that reproduces the electron-electron interaction, and is the electronic density.

Periodic boundary conditions are the preferred way to deal with condensed-matter systems. Usually the electric field enters the Schrödinger equations through an external potential field in the Hamiltonian; however, the scalar potential associated to a uniform electric field does not satisfy the periodic boundary conditions. Instead, we use a gauge where the uniform electric field is generated using a uniform, but time-dependent, vector potential as Bertsch et al. (2000)


This additional vector potential is coupled to the electrons through the kinetic energy in the KS Hamiltonian


where is the potential generated by ions and includes the Hartree, exchange and correlation potentials (HXC), that replicate the electron-electron interaction in the KS picture.

The fundamental observable in our approach is the time-dependent macroscopic current density . In a single-particle picture, the gauge-invariant current density can be calculated as


where is the volume of the unit cell. This generalized form based on the commutator has the advantage that it is explicitly invariant under translations, and that it is correct when contains non-local pseudopotential terms.

The basic macroscopic quantity for characterizing the capacity of a medium to conduct electricity is the conductivity, . It relates, to first order, a monochromatic electric field, , of frequency with the current density it generates at that frequency, , by Ohm’s law  Landau and Lifshitz (1984)


If we consider Ohm’s law in real-time Allen (2006), we can obtain the conductivity from an electron dynamics simulation by applying a time-dependent electric field , and measuring for each time the induced current . It is convenient to perturb the system using a so-called kick Yabana and Bertsch (1996), an electric field of the form , where determines the intensity and direction of the field. By inserting the kick electric field into Ohm’s law, we find that the conductivity can be calculated as


Here we assume that the conductivity is a scalar, as the systems we are interested in this work are isotropic on average. For a more general case where the conductivity is a tensor, three different simulations are required, with pointing in different directions.

We implemented this scheme in the real-space TDDFT code Octopus Castro et al. (2006); Andrade et al. (2015). The computational cost of our implementation scales quadratically with the number of atoms and does not require unoccupied states, while linear-response approaches typically scale cubically and require the calculation of unoccupied states. This fact, combined with the high parallel efficiency of real-time TDDFT Andrade et al. (2012); Draeger et al. (2016), allows us to simulate supercells with hundreds of atoms with a reasonable computational cost. For the real-time TDDFT simulations in this article we use the adiabatic local density approximation (ALDA), a grid spacing of , and a time step of .

Figure 1: Current density vs. time after a kick of (continous blue line) for a specific ionic configuration of a supercell with 256 Aluminum atoms in the liquid phase. The inset shows in detail the current fluctuations around zero for long times. The red segment line shows an exponential fit to the current, based on the range . The dotted black line shows the spatial average of the absolute value of the local current density, which illustrates how the macroscopic current density decays due to a randomization of the microscopic local current density, which retains a finite intensity.
Figure 2: Stream lines of the current density for different times for a kick of in liquid aluminum. Initially the lines are parallel and follow the direction of the current (from left to right), as time passes the current density becomes more disordered and isotropic. This disorder in the direction of the local current density is partially responsible for making the average current density vanish. The stream lines represent the trajectories of fictitious particles under the current density as a velocity field.

Ii Results

As a test case for our formalism we consider liquid aluminum at high temperature () at solid density (). The liquid state is chosen to simplify the problem as in this state, the intrinsic disorder avoids pathologies that can arise in a perfectly crystalline system and its associated translation symmetries. Our simulation supercell contains 256 atoms and we use the gamma point for Brillouin zone sampling. To account for the ionic temperature, we perform a Born-Oppenheimer molecular dynamics (MD) simulation at the target temperature and take snapshots of the ionic positions.

We start by showing the results for a single simulation using one snapshot of the MD run. In Fig. 1 we plot the current density evolving in time after a kick of As seen in the figure, the current steadily decays towards zero. Still, small fluctuations of the current are observed after it has decayed (as seen on the inset of Fig. 1). If we fit an exponential to the initial part of the current density (), we can see that the current is very close to an exponential function.

Figure 3: Frequency-dependent conductivity in the low frequency regime (blue segmented line), calculated as a Fourier transform of the time-dependent current density (Eq. 6). We compare with the frequency-dependent conductivity calculated by the Kubo-Greenwood method (black continuous line), for the same ionic snapshot. 3000 bands are used, which allows for vertical excitation energies of up to 40 eV. Additionally, a Drude fit, obtained from the Fourier transform of the exponential fit of the real time current, is presented for comparison (red dotted line). These results are for a single snapshot for a fair comparison between the two methods.
Figure 4: Ratio of the initial current per atom with respect to the applied kick. This number is equivalent to the fraction of the electrons per atom that participate in the current at time zero. The value grows as the intensity of the kick increases and approaches 3, the valence electrons in aluminum.

The real-time simulation allows us to get some insight into the nature of the decay process. Since the current density is the spatial average of the microscopic current density, the current density might decay because the microscopic current density is decaying to zero at each point, which is to be expected of a final time-reversible state. Alternatively, it is possible that the microscopic current density becomes disordered, making the average zero while retaining a finite value at each point.

To differentiate between those two scenarios we consider the microscopic current density (defined as the integrand in Eq. 4). In Fig. 1 we include the average value of the norm of the microscopic current density, and in Fig. 2 we plot the streamlines of the current density for different times.

Initially, both the current density and its average norm decay in a similar fashion. However, after some time the average norm remains constant while the average current continues decaying. This tells us that both mechanisms are at play: initially () there is a local decay of the magnitude of the current, but later () the disorder and spatial randomness of the current is responsible for nullifying the average current. The immediate conclusion is that the final state, although is time-reversible on its macroscopic averages is still a time-dependent quantum state without time-reversal symmetry. We note, however, that given the conservative and time-reversible nature of the equations (Eq. 1) and the accuracy of our propagator Castro et al. (2004), if we invert time and propagate the system backwards, starting from the final quantum state, we observe the same dynamics in reverse, with the spontaneous regeneration of the current.

To obtain the frequency-dependent conductivity we Fourier-transform the time-dependent current density, according to Eq. 6. We can do a numerical Fourier transform of the current, or analytically transform the exponential fit, that gives us a Drude-like form for the conductivity. Both results are shown in Fig. 3, where we compare them with the results from Kubo-Greenwood Kubo (1957); Greenwood (1958) for the same single ionic configuration. It can be observed that there is a good agreement between real-time TDDFT and Kubo-Greenwood. In Kubo-Greenwood the conductivity is a sum over discrete transitions between orbitals, so a broadening of is applied to produce a continuous curve, still, the conductivity is much rougher than the real-time approach. The small differences in the results can be attributed to the additional approximations introduced in Kubo-Greenwood: the truncation of the sum over unoccupied states, and the lack of a full self-consistent response.

The results we have shown so far correspond to a single snapshot, that serve to illustrate the detailed properties of our new method. However, if we want to make predictions, we need to average the results over multiple ionic configurations obtained from MD. We took 200 ionic snapshots from a equilibrated MD run using VASP Kresse and Furthmüller (1996) at the LDA level of approximation to the exchange-correlation functional and a 2x2x2 Monkhorst-Pack grid. We used a MD time step of 2 fs. For each sample we calculated the evolution of the current for different kick intensities from to , that allows us to access linear and non-linear effects.

Figure 5: Time-evolution of the current density for a family of initial conditions with different intensities of electric kick intensities, and therefore initial current densities, averaged over ionic configurations. For initial current below we observe a decay with the same time constant, as expected from a linear response. For larger perturbations we enter into a non-linear regime that is evidenced by the changing slopes of the graphs.

We first check how the initial current depends on the kick by calculating the ratio of the initial current density over the kick intensity, which we plot in Fig. 4. Since the current corresponds to the momentum of the electrons, it can be shown that this ratio is essentially the amount of charge that participates in the conduction. In a disordered metallic system this value, when normalized per atom, should be close to the valence charge. We see that for small kicks, below , a constant fraction of around electrons per atom are excited, while for higher intensity kicks this value starts to grow and saturates towards 3 for the kick considered, which is the valence charge for aluminum.

We now consider the actual decay behavior to determine the linearity of the response with respect to the intensity of the perturbation. For each kick value we obtain the current density, averaged over configurations, as a function of time. The values are shown, in logarithmic scale, in Fig. 5. If we consider that the slope of the current decay gives us an idea of the value of the conductivity, we can clearly see that there is a non-linear behavior. As the intensity of the initial current increases, the decay is less exponential-like and the rate of decay changes, first increasing and then decreasing.

To further analyze the non-linearity, we determine the DC conductivity for each curve in Fig. 5 by an exponential fit to the initial part of the time-dependent current (). The resulting conductivity vs. initial current is shown in Fig. 6. The conductivity exhibits a strong dependency with the current density. For low currents, we obtain a conductivity that is close to the extrapolated Kubo-Greenwood value of . These values also agree with previous theoretical results Ashcroft and Guild (1965); Silvestrelli (1999); Desjarlais et al. (2002); Recoules and Crocombette (2005); Vlček et al. (2012). For higher intensities, the conductivity exhibits a minimum value of for which corresponds to a kick, or electronic velocity, of , comparable to the Fermi velocity we calculate for this system: For larger currents the conductivity increases rapidly.

Figure 6: DC-conductivity as a function of the current density. For low current the value of the conductivity agrees with Kubo-Greenwood and at large currents there is a minimum of conductivity at . At even larger currents the conductivity increases as the electrons become less affected by scattering processes. The points encased in a black circle indicate the values of the current shown in Fig. 7.
Figure 7: Regions of accumulation of electronic charge for different values of the initial current, for a single snapshot. The red surfaces represent isosurfaces of the difference of the density at time with respect to the ground state density. The value of the isosurfaces are proportional to the perturbation intensity and are , , , , and , respectively.

To explore the origin of the non-linearity in the conductivity, we plot in Fig. 7 the change in the electronic density for different initial currents for a particular snapshot and a fixed time. It can be seen that there is a clear change in regime depending on the initial currents. For low currents there is some small accumulation of charge, mainly in the interatomic regions. For larger currents the charge starts to accumulate close to the nuclei. For , which is close to the conductivity minimum, the accumulation of the charge occurs in front of the atoms. As the current is increased the charge accumulation regions are centered around the atoms, and in some cases some appendices form behind the atoms. Finally, for large currents these regions of charge accumulation start to disappear.

We theorize that the generation of these regions of charge accumulation increase the amount of scattering that the electrons are subject to, reducing the conductivity. At the same time, for larger currents we expect the conductivity to increase, as the electrons become ballistic and are less affected by scattering processes. The combination of these two effects explain the shape of the conductivity curve that we see in Fig. 6.

Based on the Linhard ion-in-jellium scattering theory Lindhard (1954), we have developed a simple model to approximate the conductivity of metals. The formula proposed is derived from electronic stopping theory Schleife et al. (2015), where the collective force on the ions is reinterpreted as a deacceleration of the electronic stream. An adjustable spatial cutoff is introduced to capture a finite radius of action, , for each individual scatterer. The conductivity, as a function of the current, in our model is given by


where is the number of participating electrons per ion, and is some appropriate model of the dielectric response. In the model, is assumed to be , based on the valence of aluminum and backed up by the results of Fig. 4. is set to one fourth of the average interatomic distance, which replicates the high current portion of the results.

The model is formally justified at large currents as it is perturbative in ; however as we can see from Fig. 6, it replicates qualitatively well the shape of the simulated conductivity, reproducing the main features of the curve: a constant value for low currents, a minimum close to the current corresponding to the Fermi velocity, and a conductivity that increases rapidly for large values. This divergence behaves as , as one would expect from the Born Born (1926) and Coulomb-logarithm scattering theory Sivukhin (1966).

It is remarkable that this duality exists between the theories of the perturbative and non-perturbative regimes: the non-linear conductivity in the high-current limit appears to be related to the linear dielectric response . Of course, both theories are linear in different parameters, the former is perturbative in and the latter is perturbative in the external field.

Figure 8: Current density vs electric field curve for liquid aluminum calculated from the current dependency of the conductivity. This curve predicts that aluminum is an S-shape negative-differential conductor for very high currents and fields.

Based on the conductivity results for each current density, we can obtain the electric field vs. current density curve shown in Fig. 8. It shows that our simulations predict a negative differential conductivity in liquid aluminum at for electric fields of the order of . The induced current density in this regime is of the order of . We have not found any indication that this effect has been directly measured, as such regimes do not seem to be experimentally reachable today. The negative differential conductivity that we find is of S-type, or current-controlled Ridley (1963). In practice, the multi-valued current density for a given field causes the system to separate in different spatial regions with different currents. We expect that the curve will shift again to the right, completing the S shape, for higher currents, as semi-core electrons will become involved in the conduction.

Iii Conclusions

We have presented a new approach to study electrical conduction, and other current-related phenomena, in metallic materials. Rather than studying response properties, we directly apply an electric field on the system and follow, in real-time, the evolution of the current density. This method has practical advantages over the standard Kubo-Greenwood approach, traditionally used to calculate the electrical conductivity of materials. Interestingly, we show a fundamentally irreversible process in the simulation of a quantum system of many electrons from microscopically reversible equations.

The most important feature of the new method, however, is that it uses finite electric fields and currents, allowing us to explore non-linear effects in the conductivity. We apply this feature to the conductivity of liquid aluminum for very-high current densities (). We find that in our simulations the conductivity of aluminum starts to change under these conditions, leading us to predict that aluminum presents negative differential conductivity. To the best of our knowledge, negative differential conductivity has not been observed in metals, probably due to the large electric fields and currents involved. However, such non-linear effects could be relevant in microscopic regions where the current density could reach much higher values than the macroscopic average. It is also possible that the non-linear regime becomes accessible when the material is under more extreme conditions of temperature or pressure; for example an analogous effect, the ‘two-stream instability’, is known to develop in classical plasmas Chen (1984). Future applications of our method will be to study the conductivity of semiconductor and insulator materials, where negative differential conductivity is well known, and routinely used in electronic devices Gunn (1963); Gružinskis et al. (1999).

Our real-time method is certainly not restricted to electrical conductivity in metals, which is one of the simplest applications of modeling currents in materials. This first application opens the path for studying other phenomena that up to now have been beyond the reach of realistic quantum simulations. For example, an interesting extension of our approach would be to study electronic thermal conductivity and mixed electrical-thermal transport properties, so that we can simulate, understand and predict non-linear effects in heat transport, including negative-differential heat conduction Hu et al. (2011); Zhu et al. (2012); Zhou and Zhang (2016). In our framework, it will require the introduction of recent developments of the description of thermal properties in TDDFT Eich et al. (2014); Tatara (2015); Eich et al. (2016). We could also study the molecular dynamics of a material under an electric field or a current. A particularly interesting application of this idea is electromigration Ho and Kwok (1989): the movement of ions in a material due to the momentum transfer from the electronic current. Electromigration is very important from a technological point of view, as it is one of the main causes of failure in integrated circuits Malone and Hummel (1997). Our ultimate goal is to build a general theoretical and computational tool that can combine different types of perturbations and observables into time-resolved simulations. By accessing combined responses and non-linear effects, this tool would give us an unprecedented capacity to study different phenomena that are of interest in condensed matter physics, material science, and plasma physics.

This work was performed under DOE Contract No. DE-AC52-07NA27344. Work performed at the Energy Dissipation to Defect Evolution Center, an Energy Frontier Research Center funded by the U.S. Department of Energy (Award Number 2014ORNL1026). Computing support for this work came from the LLNL Institutional Computing Grand Challenge program. LLNL IM release number LLNL-JRNL-715217.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description