A standard model of electric conduction for molecular dynamics simulation

Physics Procedia Physics Procedia 00 (2010) 1–3

Physics Procedia 7 (2010) 107–111 www.elsevier.com/locate/procedia

A Standard Model of Electric Conduction for Molecular Dynamics Simulation Tatsuro Yugea,1 , Akira Shimizua , Nobuyasu Itob a Department b Department

of Basic Science, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan of Applied Physics, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

Abstract We propose a simple microscopic model of electric conduction, and study transport phenomena using molecular dynamics simulation. We observe that nonequilibrium steady states are realized in the nonlinear response regime as well as in the linear response regime, and that the electric conductivity is well-defined. It is confirmed that the fluctuation-dissipation theorem and the Kramers-Kronig relation hold. Furthermore the long-time tail of current autocorrelation functions which is observed in systems without impurity and causes anomalous transport coefficient, is sufficiently suppressed. This model will provide a standard model to study a nonequilibrium statistical mechanics c 2010 Published by Elsevier Ltd. of particle transport.

1. Introduction Transport phenomena have been one of the most interesting subjects in nonequilibrium thermodynamics and statistical physics. In particular, the problem of reproducing the phenomena from microscopic level is yet to be solved. Although it was partly solved by the linear response theory [1, 2], in which it is assumed that the linear response of open system agrees with that of the lowest order perturbation to closed system, there remains a long way from the complete solution and it is necessary to construct a statistical mechanics beyond the linear response theory. One of great difficulties in solving this problem is the difficulty in calculating nonequilibrium states theoretically. A promising approach seems to be computer simulations of nonequilibrium states with the powerful methods of computational physics. In fact, for heat transport system it was recently shown using nonequilibrium molecular dynamics (MD) simulation that the Fourier-type heat conduction is reproduced in basic models [3, 4, 5]. In heat conduction, however, momentum transfer is absent and high temperature difference (ΔT ∼ 104 K) is impossible because materials melt. In electric conduction, by contrast, momentum transfer is essential and high chemical potential difference(Δμ ∼ 1eV ∼ 104 K) is easily realized, which leads to the breakdown of the local equilibrium [6]. Therefore, electric conduction is the best stage to explore transport phenomena. In the present paper, we propose a new model of electric conduction, which has all the essential elements of particle transport. Using MD simulation, we study nonequilibrium as well as equilibrium properties of the model, and show that the model successfully reproduces transport phenomena [7]. 1 Current

adress: IIAIR, Tohoku University, 6-3 Aoba-ku, Sendai, Miyagi, 980-8578, Japan

c 2010 Published by Elsevier Ltd. 1875-3892 doi:10.1016/j.phpro.2010.09.054

108

T. Yuge et al. / Physics Procedia 7 (2010) 107–111 T. Yuge et al. / Physics Procedia 00 (2010) 1–3

2

2. Model Before constructing a model, let us summarize the essential elements in nonequilibrium steady states in electric conduction. In a typical experiment on electric conduction, a uniform conductor is mounted on a large sample holder. Although the conductor is macroscopically uniform, it contains random impurities which violate the microscopic translational invariance. The energy supplied from an external electric field to the conductor, is dissipated into the sample holder through the walls of the conductor, as the Joule heat. The heat flow in the conductor is mediated by electrons and phonons, while the heat flow across the walls of the conductor is mediated not by electrons but by phonons. Hence electrons must transfer their energy to phonons in order to transfer heat outside the conductor. Therefore the electron-phonon interactions and heat contact of phonons at the walls are essential to realize a nonequilibrium steady state. Moreover, in order to construct a thermodynamics, a nonequilibrium steady state must be specified by a small number of macroscopic parameters. We now propose a model which has all these essential elements. The system is composed of three types of classical particles, which we call electrons, phonons and impurities. For simplicity, we assume a two-dimensional system, the size of which is L x × Ly . There exists an external electric field E in the x-direction, which acts on only electrons, and the boundary condition in this direction is periodic. The boundaries in the y-direction are potential walls for the electrons, and thermal walls with temperature T for the phonons, and phonons carry heat to outside of the conductor through these thermal walls. The impurities are at rest and break translational invariance in the bulk region. Moreover, we assume that short-range interactions are present among all particles. Note that this model is an almost purely mechanical except for the thermal walls for the phonons. The typical systems corresponding to this model are doped semiconductors at the room temperature. Because in semiconductor systems the Fermi energy is much smaller than the room temperature, electrons can be treated as classical particles. Around the room temperature, the number density of phonons is high, and hence it seems irrelevant whether the number of phonons conserves. Furthermore, the energy-momentum dispersion relations of phonons are complicated in real conductors. Such complications, however, seem unimportant when discussing general properties of nonequilibrium states, although they will be important when discussing properties of individual materials. We therefore model here the phonons by particles whose number is conserved, with a parabolic dispersion relation with constant mass. The radius of the interaction potential range is the same order of Debye-H¨uckel wave length in classical region, and is estimated as 1nm under the conditions of our numerical experiments. 3. Numerical analysis We study transport properties of this model by MD simulation. First, we check whether the system reaches a nonequilibrium steady state in the presence of an external electric field. We calculate the electric current in the x-direction, defined as I(t) = enel Ly velx (t) where nel is the density of the Nel x vi (t)/Nel is the average velocity in the x-direction per electron, and vix is the velocity in the electrons, velx (t) = i=1 x-direction of each electron. Setting the parameters L x = 480, Ly = 24, Nel = 300, Nph = 90, Nimp = 120, we observe that I(t) increases from I(0)  0 and that after some sufficient time it fluctuates around the time-averaged value which is almost independent of the upper limit of the time range of the averaging. Therefore, a nonequilibrium steady state is realized. Next, we calculate the time-averaged values of I(t) for five different realizations of impurity arrangement. The results indicate that they are almost the same among the five distributions for each value of E. Therefore, the electric conductivity, σ = I/ELy , is almost independent of the impurity arrangement, and steady states can be specified by the small number of macroscopic parameters, T, E, L x , Ly , Nel , Nph , Nimp . We also study the system-size dependence of σ. We calculate the average value I of the electric current for L x = 320, 480 and 640 with keeping Ly and the particle densities constant. We see that I and thus σ are almost independent of the system size.  Next, we investigate the complex admittance Y(ω) =  I(ω)/E(ω)L x , which is the response to an AC electric field. We apply E(t) = E0 sin(ωt) to the system, where E0 = 0.010 (which corresponds to the linear response regime). The

109

T. Yuge et al. / Physics Procedia 7 (2010) 107–111 T. Yuge et al. / Physics Procedia 00 (2010) 1–3

3

result shown in Fig. 1 indicates that both the real and imaginary parts are well fitted with the Drude formula, Re[Y(ω)] =

e2 nel Ly τ e2 nel Ly τ2 ω   ,  , Im[Y(ω)] = mel L x 1 + (τω)2 ) mel L x 1 + (τω)2

(1)

where τ  5.5. Therefore, the Kramers-Kronig relation holds for this model. Recall that our model has multiscatterings and many-body interactions, while the Drude model does not. It is rather surprising that in spite of these differences the Drude formula seems to be valid in our model. We also study whether the fluctuation-dissipation theorem (FDT) holds for this model. The FDT for classical electric conduction is expressed as gI (ω) = 2kB T Re[Y(ω)],

(2)

where gI (ω) is the spectral intensity of electric current of equilibrium state at E = 0, and Y(ω) is the complex admittance calculated in Fig. 1. In Fig. 2 we compare the left hand side (LHS) and the right hand side (RHS) of eq. (2). We observe that the LHS and RHS agree within the error bars. Therefore the FDT holds at all frequencies for this model. Finally, we investigate the long-time tails, i.e., power-law decay of time correlation. Two types of long-time tails are well known. One is the fluid type, which decays as t−1 in two-dimensional systems [8] and which results in logarithmic divergence (with the system size) of the transport coefficient. The other is the Lorentz-gas type. The Lorentz gas model describes the motion of a point particle in random potential. This type of tail is the negative tail decaying as −t−2 in two-dimensional systems [9] and does not lead to divergence of the transport coefficient. With changing the number of impurities, we calculate the autocorrelation function v1 (t)v1 (0) of the velocity of an electron at equilibrium [10]. The parameters are set to be L x = Ly = 480, Nel = 5000, Nph = 0, and Nimp = 0–5000. The result is shown in Fig. 3. We see that v1 (t)v1 (0) decreases exponentially from the initial time and then shows power-law decay. When Nimp = 0, the fluid-type tail exists. As Nimp increases, this type of tail disappears and instead the Lorentz-type tail appears. This observation indicates that the fluid-type long-time tails are sufficiently suppressed in the systems with appropriately large number of impurities. Reflecting this fact, the autocorrelation function I(t)I(0) of electric current does not have the fluid-type long-time tail, while we can not observe the Lorentz-type one in the precision of the present work. The reason of the absence of the fluid-type tail in the present model would be that scatterings by impurities introduce a natural cut-off of fluid modes, and therefore impurities are necessary to construct a universal theory of particle transport. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

H. Nakano, Prog. Theor. Phys. 15, 77 (1956); ibid. 17, 145 (1957). R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). T. Shimada, T. Murakami, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 69, 3150 (2000). T. Murakami, T. Shimada, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 72, 1049 (2003). F. Ogushi, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 74, 827 (2005). A. Shimizu and M. Ueda: Phys. Rev. Lett. 69, 1403 (1992). T. Yuge, N. Ito and A. Shimizu: J. Phys. Soc. Jpn. 74, 1895 (2005). B. J. Alder and T. E. Wainwright: Phys. Rev. A 1, 18 (1970). M. H. Ernst and A. Weijland: Phys. Lett. A 34, 39 (1971). T. Yuge and A. Shimizu: in preparation.

110

T. Yuge et al. / Physics Procedia 7 (2010) 107–111 T. Yuge et al. / Physics Procedia 00 (2010) 1–3

4

0.008 complex admittance Y (ω)

0.007 0.006 0.005 0.004 Re[Y(ω)] Im[Y(ω)] 0.003 0.002 0.001 0 0.001

0.01 0.1 angular frequency ω

1

Figure 1: The complex admittance, plotted against angular frequency. Triangles and squares indicate the real and imaginary parts of the admittance, respectively. Data points are averaged values over five distributions of impurities. The error bars indicate the standard deviations over them. The solid curves show the Drude formula.

0.018

LHS and RHS of FDT

0.016

0.02

0.014

LHS RHS

0.016

0.012

0.012

0.01

0.008

0.008

0.004

0.006

0 0.001

0.004

0.01

0.1

1

0.2 0.4 0.6 0.8 angular frequency ω

1

0.002 0

0

Figure 2: LHS (circles) and RHS (triangles) of the FDT, eq. (2), plotted against angular frequency. The inset shows the semi-logarithmic plot of the same data. Data points are averaged values over 20 and five distributions of impurities for LHS and RHS, respectively. The error bars indicate the standard deviations over them.

111

T. Yuge et al. / Physics Procedia 7 (2010) 107–111

autocorrelation function

T. Yuge et al. / Physics Procedia 00 (2010) 1–3

0.004 0.003 0

0.002 0.001

500 1000

0

100 200

2000

-0.001 -0.002

5000

-0.003 -0.004

0

100

200

t

300

400

500

Figure 3: The autocorrelation function v1 (t)v1 (0). The values near the curves imply the numbers of impurities.

5