Non-equilibrium phonon distribution caused by an electrical current

In an attempt to explain flash sintering experiments, it had been proposed that the electron–phonon coupling leads to a proliferation of short wave-length lattice vibrations in an electric field. In this paper we investigate this by solving two coupled Boltzmann equations, describing a free electron gas in an electric field scattering from a crystal lattice coupled to a heat bath. The electric field imposes cylindrical symmetry and drives the electrons and the phonons into a non-equilibrium steady state. We find that the phonon distribution shows a strong excess population at the Brillouin zone edge in the direction of the electric field. We argue analytically, that this can be traced back to the shifted Fermi sphere for the electrons. Furthermore, not only energy but also momentum is exchanged in the electron–phonon system, which defies any attempt at describing the system by a two-temperature model.


Introduction
In a flash sintering experiment, a suitable combination of electric field and furnace heating leads to a complete densification of a green body within seconds and at significantly lower temperatures compared to conventional sintering. This considerable improvement on time and energy consumption was discovered in 2010 by Cologna et al on Zirconia [1]. The phenomenon is always accompanied by a surge in electrical conductivity and emission of bright visible light, which define the flash event even for single crystalline samples [2]. The explanation of this phenomenon is still controversial. One hypothesis that is able to describe all these effects in a wide variety of materials is that the electric field somehow causes a rapid generation of Frenkel defects in the sample [3,4]. Recent molecular dynamics simulations [5,6] have shown that proliferation of short wavelength lattice vibrations can cause a Frenkel defect concentration 10 15 times higher than the thermal concentration, enough to explain flash sintering. What remains to be understood, though, is, whether the electric field can cause such phonon proliferation at the Brillouin zone edge. The purpose of this paper is to show that this is the case.
More generally, one may ask the question, whether an electric current can drive the phonon system into a non-equilibrium steady state, and what its characteristics are. This question is related to recent research on the relaxation of optical excitations of the electron system due to electron-phonon interaction [7][8][9]. However, there is a significant difference: the relaxation after a short laser pulse leads ultimately to thermal equilibrium between the electron and the phonon system, with a transient that can be described approximately by a two temperature model. A steady state in the presence of an electrical current, by contrast, cannot be expected to be anything close to a thermal distribution. It turns out that the transient is surprisingly complex and evades any attempt to be describable in terms of thermal distributions.
With novel, sophisticated experimental techniques it became feasible to measure deviations from thermal distributions in great detail, even with temporal resolution. Recent evidence for a non-equilibrium phonon distribution was reported by Hanisch-Blicharski et al [10]. This paper also provides an example for the impact on macroscopic processes (in this case the cooling of a thin metal film), which can be explained microscopically by the non-equilibrium phonon distribution.
In the present paper we investigate the effects of an electric field on a phonon distribution in a metal. For this purpose, we propose a phenomenological model of a free electron gas in an electric field in section 2. To maintain the electric field inside the metal, a constant bias voltage is applied. This causes an electrical current density j proportional to the electric field. The electrons scatter from the crystal lattice, which leads to a net transfer of momentum to the phonons. The electron-phonon system interacts with the environment. We model this interaction by coupling the phonon system to a heat bath of a fixed temperature T hb and at rest, which acts as an energy and momentum sink for the system. The electron and phonon distributions are calculated from two coupled Boltzmann equations. The key difference to previous work, where the electron system was excited at time t = 0 by a laser pulse [7-9, 11, 12], is that our excitation term models the effects of a constant electric field, hence imposing a cylindrical symmetry as opposed to the previous spherical one. We solve the computational challenges of that in section 3. Moreover, the electron system as well as the phonon system adopt a nonzero total momentum, which requires a substantial modification of the theory. The results and the summary are given in sections 4 and 5, respectively.

Modeling approach
In order to study the influence of an electric field on a coupled electron-phonon system, we describe the electrons with a distribution function f( k) and the phonons with a distribution function g( q). In thermal equilibrium, the distribution function of the electrons is a Fermi-Dirac distribution with the electron temperature T e , the chemical potential μ e and the free electron dispersion function ( k) = 2 k 2 /2m. On the other hand, the distribution of the phonons in thermal equilibrium is a Bose-Einstein distribution with the phonon temperature T ph and dispersion function ω( q). In our model we consider three acoustic branches with Debye's linear dispersion ω(q) = c S q and the speed of sound c S . For simplicity, we do not distinguish between the different sound velocities for transversal and longitudinal waves. Furthermore, we approximate the Brillouin zone by a sphere | q| π/a, where a denotes the lattice constant. The time evolution of the distribution functions is governed by two coupled Boltzmann equations For a better visualisation of our model and the five terms in the Boltzmann equations, figure 1 shows a schematic illustration of our system and the involved interactions.

Electron acceleration by the electric field
The field acts on the electrons with a force F = −e E with the elementary charge e. It gives rise to the term [13] ∂f( k) ∂t in the Boltzmann equation. The solution of equation (4) is f( k, t) = f( k + e Et). Therefore, the electric field will cause the center of the electron distribution f( k) determined by to drift in k-space along an axis parallel to E without changing the shape of the distribution. Hence, if we start in thermal equilibrium with a Fermi-Dirac distribution in k-space, the electric field will only push the Fermi sphere while keeping its shape and therefore its temperature and chemical potential unchanged. In this respect, it is helpful to consider the total energy of the electrons consisting of kinetic center of mass energy and inner energy U e associated with the shape of the distribution f , which then gives the total electron energy E e = U e + 2 k 2 S /2m k f( k).

Electron-phonon scattering
The accelerated electrons will be slowed down due to the scattering with the phonons described by the collision term [13,14] ∂f( k) ∂t where the factor 3 accounts for the three branches of the phonon dispersion. The transfer function describes transitions between k + q and k by emitting (first term) or absorbing (second term) a phonon with wave vector q, whereas describes transitions between k − q and k by absorbing (first term) or emitting (second term) a phonon with wave vector q. The electron-phonon matrix element is given by [15] with the susceptibility 0 , Debye's length 1/σ and the crystal volume V. The delta functions in equation (6) ensure the conservation of energy during a scattering process. Correspondingly, the electron-phonon scattering changes the phonon distribution by [13,14] ∂g( q) where the factor 2 accounts for spin degeneracy of the electron states. The coupling between the electron system and the phonon system, described by (6) and (10), enables both energy and momentum exchange. Obviously, these equations obey momentum conservation and energy conservation As mentioned above, the electrons acquire a non-zero total momentum due to the electric field. The net effect of their interaction with the lattice is that momentum is transferred to the phonon system. In a steady state, the electron acceleration by the electric field will be compensated by this transfer. Hence, the phonon distribution gets a non-zero total momentum q S E, as well. The electrical current drives the phonons out of equilibrium. In a steady state, the total lattice momentum is, of course, limited by the coupling to the heat bath (or environment) and by Umklapp-processes, which transfer momentum to the crystal as a whole [15].

Phonon relaxation
What happens with the non-equilibrium phonon distribution, if the electric field as well as the coupling to the heat bath are switched off? In the long run, anharmonic effects, i.e. phonon-phonon interactions, will lead toward a Bose-Einstein distribution g BE ( q; T ph ) [15]. The total momentum will be lost due to Umklapp-processes, while the total energy will be redistributed among the phonon modes, giving rise to the temperature T ph determined implicitly by energy conservation, A simple relaxation time ansatz with a single relaxation time τ ph would then read However, if the electric field remains switched on, there is continuous energy flux from the electrons into the phonon system. The heat bath is then vital for obtaining a steady state, because it acts as an energy sink by imposing a temperature T hb . As a result, the relaxation time ansatz (14) must be modified. We postulate that the target temperature of the relaxation process will no longer be T ph , but will be replaced by T = T ph − δT < T ph , accounting for the energy transfer to the heat bath: We determine T from Fourier's law of the heat current from the lattice to the heat bath, with the heat conductance G Q . We will be looking at the system not far below the Debye temperature, which makes it safe to assume a constant heat capacity C for the phonons. With this we can set with τ hb := C/G Q identified as the relaxation time of the phonons with the heat bath.

Redistribution of electron energy and momentum
The relaxation of the phonons acts back on the electrons via the g-dependence in (7) and (8). The electron-phonon scattering not only decelerates the electrons, but also leads to Joule heating, hence an increase of the electronic inner energy U e . Since the electrons interact with each other, this inner energy is redistributed among them. During this relaxation process in the free electron gas, not only the inner energy U e , but also the total momentum k S are constant. We describe this momentum conserving electron relaxation with the ansatz with the relaxation time τ e . This means that the electron distribution relaxes toward a shifted Fermi-Dirac distribution, which keeps the total momentum of the electrons conserved. The temperature T e and the chemical potential μ e are determined such that the inner energy and the density of the electrons are conserved during this process, which means that the conditions , μ e (t))), , μ e (t))) (19) must hold at all times.

Computation of the electron-phonon interaction
Examining the symmetries of the different contributions to the Boltzmann equation (3), we note that the electron-phonon interaction, as well as the thermalizing terms exhibit spherical symmetry, as they do not distinguish one particular direction in space. Therefore, a typical approach when considering the excitation of an electron-phonon system by a laser pulse is to average over all polarisation directions, thereby keeping the spherical symmetry of the problem [12]. In our case, however, we are interested in the effects of an electric field on such systems. It enters explicitly in the excitation term and breaks the spherical symmetry of the problem. Consequently, we choose cylindrical coordinates ( e ρ , e z , e φ ) such that e z E , leaving the system symmetric in the φ-direction. The two distribution functions for the electrons and the phonons then become dependent on two coordinates f( k ) = f(ρ k , z k ) and g( q ) = g(ρ q , z q ).
In this section, we present some details, how to evaluate equations (6) and (10) in cylindrical coordinates. For large volume, the sums over k and q can be replaced approximately by the integrals and We do not limit the free electron k vectors to the Brillouin zone. In the present model, the electrons interact with the lattice exclusively via inelastic phonon scattering. For the numerical evaluation of the k integrals we make use of the fact that f( k ) vanishes for large k. This allows to limit ρ k and z k to adequately chosen finite intervals ρ k ∈ [0, ρ k,max ] and z k ∈ [z k,min , z k,max ]. As mentioned above, we consider a spherical Brillouin zone for the phonons with the radius ρ 2 q + z 2 q π/a.
The condition of energy conservation (expressed by the delta functions in equations (6) and (10)), with the free electron dispersion, rewritten in cylindrical coordinates, reads with To getḟ e−ph ( k ), we integrate over the q -space, cf equation (6). Therefore, the delta function has to be solved for one of the q coordinates, to eliminate one of the q integrals. On the other hand, forġ e−ph ( q ) we integrate over the k -space and the delta function must be solved for a k coordinate. To get the same term in both cases, we solve the delta function, i.e. equation (22), for φ q and for φ k , respectively. In both cases, there are two solutions φ 1 and φ 2 = 2π − φ 1 . Starting with the electrons, we can set φ k = 0 without loss of generality, and get , the electrons are accelerated significantly more, since there are less phonons in the system and therefore less scatterers. This is a typical result for metals, as the electric resistivity increases with increasing temperature [15].
Inserting this into equation (6) using equation (21) we can execute the φ q integration and get with z q and ρ q integration ranges discussed below. Similarly for the phonons, the delta function becomes by setting φ q = 0 without loss of generality. With this, equation (10) becomes Obviously, the integrals in equation (25) and in equation (27) only exist in the regions, where Therefore, the integration boundaries ρ q for equation (25) and ρ k for equation (27) have to be found. Whereas the latter is rather easy with ρ k = m|Ω + |/(ρ q 2 ), the integration boundaries for the ρ q integration have to be found numerically with the Newton-Raphson method. Depending on z q one obtains one or more integration segments for ρ q . The integrand has poles at the boundaries of each segment, at ρ q = ρ q,L and ρ q = ρ q,R . The integrand p(ρ q ) is a product of a known function u(ρ q ) with poles of order 1/2 at ρ q,L and ρ q,R , and a tabulated transfer function F ± at discrete values ρ q,i . In order to minimize the discretization error of the integral, we do not apply the trapezoidal rule to the integrand directly, but approximate p(ρ q ) (ρ q − ρ q,L ) · (ρ q,R − ρ q ) by a piece-wise linear function, which is devoid of poles. The integration of this function divided by (ρ q − ρ q,L ) · (ρ q,R − ρ q ) from ρ q,i to ρ q,i+1 is then performed analytically. This integration scheme has an integration error of order N −2 , where N is the number of discretization points. This is smaller than the error of order N −1 , which is anyway present due to the approximation of the spherical Brillouin zone by a cubic lattice.

Remarks on the simulation and parameters
Obviously, this study is a phenomenological one, and although the model is formulated in terms of a metal, we do not claim to describe one specific material. However, we argue that the qualitative results presented in the following are of generic nature, and that our simplified model serves as a proof of principle for these phenomena.
One must be aware, though, that the model has numerical limitations. The discretization of the k -space poses limits on (a) how steep the Fermi edge and hence how low the temperature T e can be, and (b) on the minimum possible electric field, as the distribution f has to be shifted by an amount comparable to the discretization. On the other hand, the fact that the electron distribution vanishes for large k permits to restrict the k -space. This, however, limits temperature and electric field strength. For a k -space grid of N ρ × N z = 50 × 100, we limit the electric field strength to E/E 0 ∈ [0.001, 0.007] and the heat bath temperatures to T hb /T 0 ∈ [0.3, 0.9] with our units E 0 = 2 /(ema 3 ) = 1.14 · 10 7 V cm −1 , T 0 = 2 /(ma 2 k B ) = 5680 K and the time unit t 0 = ma 2 / ≈ 1.42 fs. This corresponds to the lattice constant of aluminum, a = 4.05 Å. The electron mass is m = 9.1 · 10 −31 kg.
The sound velocity is set to c S = 0.2a/t 0 . The chemical potential μ e (T e ) is calculated so that the electron density is n e = 3/a 3 . Comparing with the conditions of a flash sintering experiment [16], the temperature is around a factor three and the electric field is around five orders of magnitude higher. However, we also need to choose the characteristic time τ hb for the coupling between the lattice and the heat bath sufficiently small in order to keep the duration of the transients in a manageable range. This conspires favorably with the electric field being too large, as the influx of energy into the system and the outflux are both too large.   Furthermore, the electron-electron interaction time is set to τ e /t 0 = 0.7, which is in the range of femtoseconds and hence a typical value. The phonon-phonon interaction is typically a thousand times slower and is therefore set to τ ph /t 0 = 700.

Non-equilibrium electron distribution
We start our simulations in thermal equilibrium, meaning f( k ) = f FD ( k , T e , μ e ) and g( q ) = g BE ( q , T ph ) with T e = T ph = T hb =: T at time t = 0. The average electron momentum k S starts at zero and increases due to the electric field. According to equation (3), the time evolution of k S is: (29) Figure 2 shows the absolute value k S (t) as a function of time for different electric fields. Two regimes can be distinguished. The first regime extends from t = 0 till t ≈ τ ph , which is the time scale, on which Umklapp processes remove momentum from the phonon system. Before t ≈ τ ph , the momentum exchange between the electron and the phonon system therefore cannot decelerate the electrons efficiently. This explains the roughly linear increase of the electron momentum due to the electric field. Above t ≈ τ ph (second regime), the electron momentum decreases because-once transferred to the phonon system-a part is irreversibly lost due to Umklapp processes. Figures 3(a) and 4(a) show that the electron distribution is essentially a Fermi-Dirac distribution shifted by the momentum k S . The difference to a shifted Fermi sphere is shown in figures 3(b) and 4(b). Due to the Pauli principle and the relatively low phonon energy, the electron-phonon scattering can only change the electron distribution close to the Fermi edge. This change is relatively slow compared to τ e , so that the relaxation due to electron-electron interaction keeps the modification of the Fermi sphere small. The The north pole occupation goes through a maximum at t/t 0 = 1400 = 2τ ph /t 0 (vertical black line). In (b) the relative north pole occupation compared to a Bose-Einstein distribution of the same amount of energy is shown. In steady state, the non-equilibrium distribution is 1.8 times higher than a thermal one, independent of the electric field or temperature. The black horizontal line marks g/g BE = 1.8. difference shown in figures 3(b) and 4(b) is small and originates from the non-thermal phonon distribution, which we describe in the next paragraph.

Non-equilibrium phonon distribution
In contrast to the electrons, the phonon distribution function for t = 0 differs significantly from a Bose-Einstein distribution. The electric field leads to a strong proliferation of phonons at the 'north pole', e.g. phonons at the Brillouin zone boundary in field direction, cf figure 5. We note that the distribution is not shown for states of small q, in order to keep the rest of the distribution visible, since for acoustic phonons g(q → 0) → ∞. The electron-phonon matrix element (9) does not distinguish between positive and negative q-vectors and therefore cannot be the reason for the proliferation at the north pole. It rather seems to be a generic effect of the electric field, in contrast to what has been assumed previously [6].
In order to understand what actually causes this proliferation, we take a look at how the phonon distribution is affected by the electrons, i.e.ġ e−ph ( q ) given in equation (10) and shown in figure 6. We see thatġ e−ph ( q ) is highest at the north pole, which means that these phonons are emitted primarily. According to (10), the transfer function F + must be responsible for the generic structure ofġ e−ph ( q ), as the electron-phonon matrix element was ruled out. This is confirmed by considering F + in q -space for f( k ) = f FD ( k − k S , T e ) and g( q ) = g BE ( q , T ph = T e ). The reason for this choice is to see how a shifted Fermi-Dirac distribution drives the Bose-Einstein distribution away from thermal equilibrium. With these assumptions, equation (7) becomes (30) This function couples the component of q in field direction z q to k S , because where we have used the energy conservation condition of a scattering process ( k + q ) = ( k ) + ω( q ). Of course, to get the fullġ e−ph ( q ) we need to sum the transfer function over k -space. However, the largest contributions are from states k at the Fermi edge. In figure 6(b) we show the transfer function for a specific k at the shifted Fermi edge and k S · a = 0.348. As one can see, the structure ofġ e−ph ( q ) is essentially produced by the transfer function with a shifted Fermi-Dirac distribution. We therefore conclude that for a proliferation of phonons at the north pole it is sufficient to have a shifted electron distribution. This means that this effect is rather robust, as it does not depend on material properties. Figure 7(a) shows that the phonon population at the north pole goes through a pronounced maximum at around t/t 0 = 1400 which compares well with τ ph /t 0 = 700 and corresponds to 2 ps in SI-units. A higher temperature T leads to a lower maximum. This can be explained with the temperature dependent resistivity: at higher temperatures, there are more phonons in the system that scatter with and decelerate the electrons. The Fermi sphere is not shifted as much (as shown in figure 2), and therefore, less phonons at the north pole are emitted. Toward the steady state, the north pole occupation decreases to about half the maximal value. The north pole occupation relative to a Bose-Einstein distribution with the same amount of energy is shown in figure 7(b). In steady state, the north pole has a non-equilibrium population 1.8 times higher than the thermal population, independent of the electric field or the temperature T. In figure 8(a) the phonon distribution in the steady state is shown. Small q-values are again left out, since for acoustic phonons g(q → 0) → ∞. The relative steady state phonon distribution compared to a Bose-Einstein distribution with the same amount of energy, g( q )/g BE ( q , T ph ), is shown in figure 8(b). In steady state, there are more phonons with q · E > 0 and less phonons with q · E < 0 compared to a spherically symmetric Bose-Einstein distribution.

Failure of the two-temperature model
The two temperature model is based on the assumption that electrons excited by a laser pulse thermalise due to their interaction much faster than the time scale of their interaction with the phonons, adopting a Fermi-Dirac distribution with a temperature T e . As the phonons are assumed to be in thermal equilibrium at a temperature T ph < T e to begin with, the equilibration of the coupled electron-phonon system can be modeled by a Fourier heat conduction law [17].
In our case, the continuous momentum transfer leads to the failure of this model. The phonons represent the only momentum sink for the electrons, and therefore the electrons cannot thermalise on their own. Instead, they can be described approximately by a shifted Fermi sphere to which one can reasonably associate an electron temperature T e . By contrast, the definition (13) of T ph is rather artificial and only gives a qualitative measure of the energy in the phonon system.
As the electrons transfer momentum to the phonon system primarily by exciting phonons at the Brillouin zone edge, the energy in the phonon system increases, and with it the associated temperature T ph . The interactions that irreversibly remove momentum from the system, namely the Umklapp-processes, happen on larger time scales. Therefore, momentum builds up both in the electron and in the phonon system. For the electrons, the associated energy is primarily given by 2 k 2 S /2m and only to a smaller extent reflected by an increase of T e . However, for the phonon system, the total energy associated with the momentum contributes to the increase of T ph , which leads to T ph > T e seen in figure 9. Of course, a higher electric field means a stronger Joule heating. This explains the steeper temperature curves for higher field strengths in figure 9.

Summary
We developed a dynamical model for describing how the occupation numbers of electron and phonon states evolve in a metal under the combined influence of an electric field and a heat bath. It is based on the coupled Boltzmann equations for the electrons and the phonons. Since the electric field imposes cylindrical symmetry, we calculated the distribution functions in cylindrical coordinates of the reciprocal space. In the model presented in this paper, a free electron gas is considered. This means that the interaction with a static periodic lattice potential is neglected, but only inelastic phonon scattering is taken into account. In principle, the dispersion function ( k ) could be folded back into the Brillouin zone, but this would require to introduce a band index for the electrons, which is not necessary in the present paper. One of the interesting future extensions of the model would be to consider semiconductors, in which case at least a valence band and a conduction band with k vectors restricted to the Brillouin zone must be considered.
We have found that the electron distribution function maintains the thermal shape of a Fermi sphere very well, however gets shifted in reciprocal space by a wave vector parallel to the electric field. This is attributed to the Pauli principle, the short relaxation time τ e , and the relatively small phonon energy, which confine the modifications due to the electron-phonon interaction to the neighborhood of the Fermi surface. The phonons on the other hand deviate strongly from a thermal distribution. Their average number increases with the electric field. More precisely, we have found that phonons with wave vectors at the Brillouin zone edge in direction of the electric field, i.e. at the north pole of the Brillouin zone, proliferate maximally for a transient of about 2 ps. After that, the relative occupation compared to a Bose-Einstein distribution with the same energy decreases to a steady state value of 1.8, as shown in figure 7(b) for electric fields E/E 0 ∈ [0.001, 0.007] and temperatures T hb /T 0 ∈ [0.3, 0.9]. Whether or not this proliferation is able to stimulate the flash phenomenon, as suggested by Jongmanns et al [5,6], remains to be shown.
We argued that the phonon proliferation has nothing to do with the electron-phonon matrix element. We showed by examining the transfer function, that a shifted Fermi sphere is reason enough, which indicates that the proliferation of phonons at the north pole is universally related to a finite electron momentum. A consequence of the momentum transfer is that the system cannot be described by a two temperature model. In fact, we show that if we describe the momentum transfer to the phonon system in terms of a temperature T ph , this temperature would be higher than that of the electrons, T ph > T e . The reason is that the non-zero total momentum of the electrons does not contribute to T e , while the definition of the phonon temperature does not count the energy contribution of the total momentum separately.
Our prediction that an electrical current drives the phonon distribution away from a Bose-Einstein distribution, and in particular leads to an overpopulation of the north pole in reciprocal space by a factor 1.8 in steady state for strong fields and high temperatures, should be experimentally observable by thermal diffuse electron and x-ray scattering [18][19][20]. A further numerical evaluation of our model for a specific material and realistic fields and temperatures is desirable and feasible, but requires much more computing power than available for the present study, as we pointed out at the beginning of section 4.