## An Electrohydrodynamics Model for Non-equilibrium Electron and

## Phonon Transport in Metal Films after Ultra-short Pulse Laser Heating

Jun Zhou, 1,2 Nianbei Li,  and Ronggui Yang 2 1, ♥

1 Department of Mechanical Engineering and Materials Science and Engineering Program,

University of Colorado, Boulder, CO 80309, USA

2 Center for Phononics and Thermal Energy Science, School of Physics Science and Engineering, Tongji University, Shanghai 200092, People's Republic of China

## ABSTRACT

The electrons and phonons in metal films after ultra-short pulse laser heating are in highly non-equilibrium states not only between the electron sub-system and the phonon sub-system  but  also  within  the  electron  sub-system.  An  electrohydrodynamics  model consisting of the balance equations of electron density, energy density of electrons, and energy  density  of  phonons  is  derived  from  the  coupled  non-equilibrium  electron  and phonon Boltzmann transport equations to study the nonlinear transport phenomena, such as the electron density fluctuation and the transient electrical current in metal films, after ultra-short pulse laser heating. The time-dependent temperature distributions is calculated by the coupled electron and phonon Boltzmann transport equations, the electrohydrodynamics  model  derived  in  this  work,  and  the  two-temperature  model  for different laser pulse durations, film thicknesses, and laser fluences. We find that the twotemperature model overestimates the electron temperature at the front surface of the film and  underestimates  the  damage  threshold  when  the  nonlinear  thermal  transport  of

♥ Email: Ronggui.Yang@Colorado.Edu

electrons is important. The electrohydrodynamics model proposed in this work could be a more  accurate  prediction  tool  to  study  the  non-equilibrium  electron  phonon  transport process  than  the  two-temperature  model  and  it  is  much  easier  to  be  solved  than  the coupled electron and phonon Boltzmann transport equations.

KEYWORDS:  electron-phonon  coupling,  laser-material  interaction,  metal  thin  film, electrohydrodynamics model, Boltzmann transport equation, two-temperature model

PACS numbers: 63.20.kd, 66.10.cd, 44.10.+i, 68.60.Dv

## I. INTRODUCTION

The study of the interaction of ultra-short pulse laser with materials has a wide range of applications in the fields of fundamental physics of electron-phonon (EP) interaction, 1,2,3 hot  electron  dynamics, 4,5,6,7,8 superconductivity, 9,10 and  laser-material processing. 11,12,13,14,15,16 The  energy  deposited  from  ultra-short  pulse  laser  could  drive electrons and phonons to a highly nonequilibrium state. For laser-metal interaction, the energy  temporarily  stored  in  the  electron  sub-system  is  conducted  away  through  two channels: i) electron transport and ii) coupling of electron lattice to the lattice which is then  conducted  away  by  phonons. 17 , 18 The  interplay  between  these  two  processes determines the evolutions of both electrons and phonons.

The research  on  this  topic  can  be  traced  back  to  as  early  as  1957.  Kaganov et  al . 17 presented  a  model  for  EP  energy  exchange  which  is  still  quite  widely  used.  This phenomenological  two-temperature  model  (TTM)  assumes  that  the  sub-systems  of electrons and phonons can be described by their own temperatures and both the electron and phonon temperatures are much higher than the Debye temperature. The EP energy exchange rate can then be linearized as )  where ( p e T T G -e T is the electron temperature, p T is the phonon temperature, and G is the EP coupling constant. The coupled electron and phonon thermal transport could thus be described by the energy balance equations 18,19

$$C _ { e } ( T _ { e } ) \, \frac { \partial T _ { e } } { \partial t } = - \nabla \cdot \mathbf J _ { e, E } - G ( T _ { e } - T _ { p } ) + Q ( \mathbf r, t ) \,,$$

$$C _ { p } ( T _ { p } ) \frac { \partial T _ { p } } { \partial t } = - \nabla \cdot \mathbf J _ { p, E } + G ( T _ { e } - T _ { p } ) \,.$$

Here )  is  the  electronic  heat  capacity, ( e e T C )  is  the  lattice  heat  capacity,  and ( p p T C )  is  the  energy  deposition  by  the  ultra-short  pulse  laser, , ( t Q r E e , J and E p , J are  the energy current densities of electrons and phonons, respectively. The TTM follows closely the  Fourier's  law  of  heat  conduction e p e e E e T T T ∇ -κ = ) , ( , J and p p p E p T T ∇ -κ = ) ( , J , where ) , ( p e e T T κ is the electronic thermal  conductivity  and ) ( p p T κ is the lattice (phononic) thermal conductivity.

Although the TTM has been widely used by many researchers, there are quite some limitations of TTM to describe the ultra-short pulse laser-material interaction, especially with  the  recent  advances  of  femtosecond  and  attosecond  lasers.  One  of  the  obvious drawback is not considering the transient electrical current flow due to non-equilibrium electrons:

- (1) e p e e E e T T T ∇ -κ = ) , ( , J used in TTM is valid only when the electrical current and the chemical potential gradient can be ignored.
- (2) The divergence of energy current density of electrons E e , J ⋅ ∇ -consists of not only Fourier's heat conduction but also the electron density accumulation, the Thomson effect, and the Joule heat generation. 20 The TTM considers the heat conduction and neglects the other mechanisms.
- (3)  The  approximation e p e e e p e e T T T T T T 2 ) , ( ] ) , ( [ ∇ ≈ ∇ ⋅ ∇ κ κ is  used  in  the  TTM.  It  is valid only when the thermal conductivity are temperature-independent. In addition, the electronic heat capacity is usually chosen to be proportional to the electron temperature which is valid only when the temperature is lower than Fermi temperature.

Currently, there have been a few works which address the non-equilibrium transport that  considers  transient  electrical  current  and  chemical  potential  gradient.  Chen et  al. 21 added an additional conservation equation of electron momentum which is derived from the BTE into the TTM to study the effect of electronic kinetic pressure on the electron and  phonon  temperature  evolutions.  As  a  result,  they  obtain  significantly  different electron and phonon temperature response from the TTM. The larger laser fluence is and shorter laser pulse duration is, the larger the difference between their model and the TTM is.

In  this  paper,  we  derive  an  electrohydrodynamics  (EHD)  model  based  on  the coupled  electron  and  phonon  BTE  to  investigate  the  electron  and  phonon  transport  in metals after the ultra-short pulse laser heating. In this model, the energy current is driven by both the temperature gradient and the chemical potential gradient. The contributions to E e , J ⋅ ∇ -from  the  heat  conduction,  the  electron  density  accumulation,  the  Thomson effect, and the Joule heat generation are considered. The paper is organized as follows. After a brief introduction of the BTE in Sec. IIA, we derive the EHD model in detail in Sec. IIB. In Sec. III, the results of the calculation of the electron density fluctuation, the transient electrical current, the evolutions of electron temperature and phonon temperature, and the laser fluence damage threshold are presented. Section IV concludes this work.

## II. THEORETICAL MODELS

Figure 1 shows the coupled non-equilibrium electron and phonon transport mechanisms in a metal film after ultra-short pulse laser heating. The surface of the film is

in x-y plane (perpendicular to z -direction) and the film thickness is L . The photon energy of the ultra-short pulse laser can be absorbed by electrons near the front surface ( x = 0) of the film. The intensity of absorption exponentially decays with the depth where δ is the optical absorption depth. Highly non-equilibrium states between the electron sub-system and the phonon sub-system are induced because the electrons are excited by ultra-short pulse  laser  to  higher  energy  states.  The  highly  non-equilibrium  electrons  will  dissipate their energy through two channels: i) electron transport and ii) coupling of electron lattice to  the  lattice  which  is  then  conducted  away  by  phonons.  Diffusive  reflection  boundary conditions  are  used  in  this  study  by  assuming  that  electrons  and  phononsare  reflected back randomly when they reach the front surface or the rear surface of the film.

## A. Coupled electron and phonon BTE model

We can write the coupled electron and phonon BTE as follows: 22

$$\frac { \partial f _ { \mathbf k } ( \mathbf r, t ) } { \partial t } + \mathbf v _ { e, \mathbf k } \cdot \nabla f _ { \mathbf k } ( \mathbf r, t ) - \frac { e } { \hbar } \, \mathbf E \cdot \nabla _ { \mathbf k } f _ { \mathbf k } ( \mathbf r, t ) = \frac { \partial f _ { \mathbf k } ( \mathbf r, t ) } { \partial t } \Big | _ { c } + s ( \mathbf r, \mathbf k, t ) \,,$$

$$\frac { \partial n _ { \mathfrak { q } } ( \mathbf r, t ) } { \partial t } + \mathbf v _ { p } \cdot \nabla n _ { \mathfrak { q } } ( \mathbf r, t ) = \frac { \partial n _ { \mathfrak { q } } ( \mathbf r, t ) } { \partial t } \Big | _ { c },$$

where )  is the electron distribution function in the phase space with wave vector , ( t f r k k and  position )  at  time , , ( z y x = r t , )  is  the  phonon  distribution  function  in  the , ( t n r q phase  space  with  wave  vector q 23 and  position r at  time t . k ∇ is  the  wave  vector gradient  operator,  and E is  the  external  electric  field. k e ε k k v ∇ = h 1 , is  the  electron velocity, where k ε is  the  kinetic  energy  of  electron  and h is  the  Planck constant. p v is

the  phonon  velocity. )  is  the  source  term  due  to  laser  heating. , , ( t s k r 24 c t t f | / ) , ( ∂ ∂ r k is the electron collision term which includes the EP scattering term ep t t f | / ) , ( ∂ ∂ r k and the electron-electron  scattering  term ee t t f | / ) , ( ∂ ∂ r k . c t t n | / ) , ( ∂ ∂ r q is  the  phonon  collision term  which  includes  the  EP  scattering  term ep t t n | / ) , ( ∂ ∂ r q and  the  phonon-phonon scattering term pp t t n | / ) , ( ∂ ∂ r q .

The EP scattering terms in Eq. (2) are calculated from a microscopic approach: 9

$$\frac { \partial f _ { \mathbf k } } { \partial t } \Big | _ { e p } = & - \frac { 2 \pi } { \hbar } \sum _ { \mathbf q } | \, M _ { \mathbf q } | ^ { 2 } \, \Big \{ \delta ( \varepsilon _ { _ { k } } - \varepsilon _ { _ { k ^ { \prime } } } - \hbar { \omega } _ { _ { q } } ) \Big [ n _ { \mathbf q } ( f _ { \mathbf k } - f _ { \mathbf k ^ { \prime } } ) + f _ { \mathbf k } ( 1 - f _ { \mathbf k } ) \Big ] \\ & - \delta ( \varepsilon _ { _ { k } } - \varepsilon _ { _ { k ^ { \prime } } } + \hbar { \omega } _ { _ { q } } ) \Big [ n _ { \mathbf q } ( f _ { \mathbf k ^ { \prime } } - f _ { \mathbf k } ) + f _ { \mathbf k ^ { \prime } } ( 1 - f _ { \mathbf k } ) \Big ] \Big ],$$

$$\frac { \partial n _ { \mathbf q } } { \partial t } \Big | _ { e p } = - \frac { 4 \pi } { \hbar } \sum _ { \mathbf k } | \, M _ { \mathbf q } \, | ^ { 2 } \ f _ { \mathbf k } ( 1 - f _ { \mathbf k } ) [ n _ { \mathbf q } \delta ( \varepsilon _ { _ { k } } - \varepsilon _ { _ { k } } + \hbar { \omega } _ { _ { q } } ) \\ - ( n _ { \mathbf q } + 1 ) \delta ( \varepsilon _ { _ { k } } - \varepsilon _ { _ { k } } - \hbar { \omega } _ { _ { q } } ) \Big ]$$

We drop ( , )  in  the  rest  of  this  paper  for  simplification.  Here t r q M is  the  EP  scattering matrix elements, q ω h is  the  energy  of  phonons  whose  wave  vector  equals  to  the  wave vector  change  of  electron  after  the  scattering ' k k q -= .  Delta  functions  denote  the energy conservation of the scattering processes. Electron spin degeneracy has been taken into account by the factor 2 in Eq. (3b).

To calculate the EP scattering terms in Eq. (3) in metal without doing the summation over three dimensional momentum  space,  the Eliashberg function 25 ) ( 2 ω α F is introduced similar as shown in Ref. [9]. We then have

$$\frac { \partial f _ { \mathbf k } } { \partial t } \Big | _ { \mathrm e p } = - 2 \pi \int d \omega \alpha _ { i } ^ { 2 } F ( \omega ) \int & d \varepsilon ^ { \left \{ } \delta ( \varepsilon _ { i } - \varepsilon ^ { \prime } - \hbar { \omega } ) \Big [ n _ { 0, \omega } ( f _ { 0, \varepsilon _ { i } } - f _ { 0, \varepsilon ^ { \prime } } ) + f _ { 0, \varepsilon _ { i } } ( 1 - f _ { 0, \varepsilon ^ { \prime } } ) \right ] \\ & - \delta ( \varepsilon _ { k } - \varepsilon ^ { \prime } + \hbar { \omega } ) \Big [ n _ { 0, \omega } ( f _ { 0, \varepsilon ^ { \prime } } - f _ { 0, \varepsilon _ { i } } ) + f _ { 0, \varepsilon ^ { \prime } } ( 1 - f _ { 0, \varepsilon _ { i } } ) \right ],$$

$$\frac { \partial n _ { \mathbf q } } { \partial t } \Big | _ { \mathrm e p } = - 2 \pi \Xi ( \mu ) \Big | \, d \varepsilon \Big | \, d \varepsilon ^ { \prime } \frac { \alpha ^ { 2 } F ( \omega ) } { F ( \omega ) } ( f _ { 0, \varepsilon } - f _ { 0, \varepsilon } ) \\ & \Big [ n _ { 0, \omega } \delta ( \varepsilon - \varepsilon ^ { + } \hbar { \omega } ) - ( n _ { 0, \omega } + 1 ) \delta ( \varepsilon - \varepsilon ^ { \prime } - \hbar { \omega } ) \Big ] \Big | _ { \omega = \omega _ { \mathbf q } },$$

where )  is  the  electron  density  states  (DOS)  at  the  chemical  potential  and ( µ Ξ )  is ( ω F the phonon DOS. k f ε , 0 is the equilibrium Fermi-Dirac distribution of electrons and ω , 0 n is the equilibrium Bose-Einstein distribution of phonons.

In  non-equilibrium systems, the electron-electron scattering and the phonon-phonon scattering drive the non-equilibrium  electron  and  phonon  sub-systems  towards  to equilibrium  ones,  respectively.  Therefore,  we  assume [ ] ) ( / | / , 0 k f f t f ee ee k τ ε --= ∂ ∂ k k 24 and [ ] )  using relaxation time approximation where ( / | / , 0 q n n t n pp pp q τ ω --= ∂ ∂ q q )  is ( k ee τ the electron-electron relaxation time and )  is the phonon-phonon relaxation time. In ( q pp τ this work, )  and ( k ee τ )  are assumed to be constants so that we could focus on the ( q pp τ EP scattering on the energy transfer pathway.

With  the  scattering  terms  explicitly  written  as  above,  the  coupled  electron  and phonon BTE model shown in Eq. (2) can be solved iteratively for both directions ( 0 &gt; z k and 0 )  by  using  finite  difference  method  with  a  numerical  scheme  similar  to  that &lt; z k described in Ref. [26]. In our calculation, the Eliashberg function is taken from Ref. [27] and the phonon DOS is taken from Ref. [28]. Both the front surface and rear surface are assumed  to  be  rough  so  that  the  equilibrium  distribution  function  due  to  diffusive reflection of electrons and phonons can be written at the boundaries: 29

$$f _ { \mathbf k } ( z = 0, t ) = f _ { 0, \varepsilon _ { \mathbf k } } ( z = 0, t ) \,, \ f _ { \mathbf k } ( z = L, t ) = f _ { 0, \varepsilon _ { \mathbf k } } ( z = L, t ) \,,$$

$$n _ { \mathfrak { q } } ( z = 0, t ) = n _ { 0, \omega _ { \mathfrak { q } } } ( z = 0, t ) \,, \ n _ { \mathfrak { q } } ( z = L, t ) = n _ { 0, \omega _ { \mathfrak { q } } } ( z = L, t ) \,.$$

Table I lists all the parameters used in the calculation of BTE.

## B. Electrohydrodynamics (EHD) Model

The coupled electron and phonon BTE is a set of multi-variable integral-differential equations  in  the  phase  space  where  the  scattering  integral  is  very  complicated  to  be solved. 30 A  convenient  approach  to  simplify  the  problem  is  to  use  balance  equations, instead of directly solving the BTE to obtain the distribution function of electrons and phonons, where the macroscopic variables such as electron density, electron temperature, and phonon temperature are calculated. The computational cost is remarkably reduced in spite of the losing of some detailed information in distribution functions. Here, we write a set  of  EHD  equations  with  the  electron  density  and  energy  density  of  electrons  as macroscopic variables by multiplying 1/ V and V k / ε in Eq. (2a) and then summing over k , where V is the unit volume. 22

$$\frac { \partial N } { \partial t } = - \nabla \cdot { \mathbf J } _ { e, n } \,,$$

$$\frac { \partial U _ { e } } { \partial t } = - \nabla \cdot \mathbf J _ { e, E } - e \mathbf E \cdot \mathbf J _ { e, n } - G ( T _ { e } - T _ { p } ) + Q ( \mathbf r, t ) \,.$$

Here ∑ = k k f N 2 is the local electron density and ∑ = k k k v J f e n e , , 2 is the particle current density of electrons (n e e , J is the electrical current density). ∑ = k k f Ue k ε 2 is the energy

density  of  electrons  and ∑ = k k k v J f k e E e ε , , 2 is  the  energy  current  density  of  electrons.

The  factor  2  comes  from  electron  spin  degeneracy. ∑ = q v J q q p p n ω h is  the  thermal

current density carried by phonons. We focus our study on noble metals, the photoelectric

effect and inter-band transition which excites d - or  -electrons to be itinerant electrons in f transition  metals  are  not  considered.  Therefore,  2 0  since  there  is  no  net ) , , ( = ∑ k k r t s photon )

electron generated or annihilated in electron intra-band transition due to absorption  as  shown  in  Fig.  1.  The  electron  energy  source  term  is ∑ = k k r , , ( 2 t s Q k ε

which is assumed to be independent on the temperature. 19 The balance equation which is the same as Eq. 1(b) can also be obtained by multiplying V q / ω h in  Eq.  (2b)  and  then summing over q .

Under the free electron approximation, the left side of Eq. (6b) can be written as:

$$\frac { \partial U _ { e } } { \partial t } = \frac { \partial U _ { e } } { \partial \mu } \frac { \partial \mu } { \partial t } + \frac { \partial U _ { e } } { \partial T _ { e } } \frac { \partial T _ { e } } { \partial t } = \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial t } + C _ { e } ( T _ { e } ) \frac { \partial T _ { e } } { \partial t } \,,$$

where  the  electronic  heat  capacity  is  defined  as e e e e e e T N N U T U T C ∂ ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂µ ∂ ∂µ ∂ -∂ ∂ = -1 ) ( which can be calculated as

$$C _ { e } ( T _ { e } ) = \gamma T _ { e } \left [ 1 - \frac { 3 } { 1 0 } \left ( \frac { k _ { B } T _ { e } } { \mu _ { 0 } } \right ) ^ { 2 } + O \left ( \frac { k _ { B } T _ { e } } { \mu _ { 0 } } \right ) ^ { 4 } \right ],$$

where 0 µ is  the  chemical  potential  at  zero  temperature. 31 See  Appendix  A  for  detailed derivations.

n e , J and E e , J in  Eqs.  (6a)  and  (6b)  can  be  obtained  by  solving  the  BTE  under  the relaxation time approximation 22 as shown in Appendix B for details:

$$\mathbf J _ { e, n } = - \frac { \sigma } { e } \mathbf E - D _ { n } \nabla N - S _ { n } \nabla T _ { e },$$

$$\mathbf J _ { e, E } = - \left ( - T _ { e } \sigma \alpha + \frac { \mu \sigma } { e } \right ) \mathbf E - D _ { E } \nabla N - S _ { E } \nabla T _ { e } \,. \quad \quad \quad ( 9 b )$$

Here σ is the electrical conductivity, 0 0 2 3 2 N e Dn µ σ = is the diffusion coefficient of electron,

$$S _ { n } = - \frac { \sigma \alpha } { e } \frac { 2 \lambda - 1 } { 2 \lambda } \quad \text{is} \quad \text{the} \quad \text{Soret} \quad \text{coefficient}, \quad \alpha \quad \text{is} \quad \text{the} \quad \text{Seebeck}$$

$$\text{coefficient}, \ \ D _ { E } = - \frac { T _ { e } \sigma \alpha } { e } \frac { 2 \mu _ { 0 } } { 3 N _ { 0 } } + \mu D _ { n } \quad, \quad \text{and} \quad S _ { E } = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \frac { 1 } { 2 \lambda } + \mu S _ { n } \quad, \quad \text{where}$$

$$\kappa _ { e } ^ { 0 } = \kappa _ { e } + T _ { e } \sigma \alpha ^ { 2 } \,. \, \text{Here } \, N _ { 0 } \text{ is the electron density without fluctuation} \, \lambda = - \alpha \left ( \frac { \pi ^ { 2 } k _ { B } ^ { 2 } T _ { e } } { 3 e \mu _ { 0 } } \right ) ^ { - 1 }$$

is almost a constant for given metal. 32

Then n e , J ⋅ ∇ -and E e , J ⋅ ∇ -in Eqs. (6a) and (6b) can be calculated as follows (see detailed derivation in Appendix C):

$$- \nabla \cdot \mathbf J _ { e, n } \approx D _ { n } \nabla ^ { 2 } N + S _ { n } \nabla ^ { 2 } T _ { e } - \frac { D _ { n } } { 3 N _ { 0 } } ( \nabla N ) ^ { 2 } - \frac { 2 S _ { n } ^ { \bullet } } { T _ { e } } ( \nabla T _ { e } ) ^ { 2 } - \frac { 2 D _ { n } } { T _ { e } } ( \nabla N \cdot \nabla T ),$$

$$- \nabla \cdot \mathbf J _ { e, E } \approx D _ { E } \nabla ^ { 2 } N + S _ { E } \nabla ^ { 2 } T _ { e } - \frac { D _ { E } } { 3 N _ { 0 } } \left ( \nabla N \right ) ^ { 2 } - \frac { 2 S _ { E } ^ { \prime } } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 } - \frac { 2 D _ { E } } { T _ { e } } \left ( \nabla N \cdot \nabla T _ { e } \right ).$$

$$\text{Here } \, S _ { \, _ { n } } ^ { \prime } = - \frac { \sigma \alpha } { \ e } \frac { 4 \lambda - 1 } { 4 \lambda } \text{ and } S _ { \, _ { E } } ^ { \prime } = \kappa _ { _ { e } } ^ { 0 } - T _ { _ { e } } \sigma \alpha ^ { 2 } \frac { 1 } { 4 \lambda } + \mu S _ { \, _ { n } } ^ { \prime } \,.$$

In the rest of this work,  we  do  not  consider  the  external electric field for simplification.  By  substituting  Eq.  (10)  into  Eq.  (6),  the  EHD  model  can  finally  be written as (see Appendix D for details):

$$\frac { \partial N } { \partial t } = D _ { n } \nabla ^ { 2 } N + S _ { n } \nabla ^ { 2 } T _ { e } - \frac { D _ { n } } { 3 N _ { 0 } } \left ( \nabla N \right ) ^ { 2 } - \frac { 2 S _ { n } ^ { \prime } } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 } - \frac { 2 D _ { n } } { T _ { e } } \left ( \nabla N \cdot \nabla T \right ),$$

$$C _ { e } ( T _ { e } ) \frac { \partial T _ { e } } { \partial t } & = D \nabla ^ { 2 } N + S \nabla ^ { 2 } T _ { e } - \frac { D } { 3 N _ { o } } \left ( \nabla N \right ) ^ { 2 } - \frac { 2 S ^ { \prime } } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 } - \frac { 2 D } { T _ { e } } \left ( \nabla N \cdot \nabla T _ { e } \right ) \\ & \quad - G ( T _ { e } - T _ { p } ) + Q,$$

$$C _ { p } ( T _ { p } ) \, \frac { \partial T _ { p } } { \partial t } = \kappa _ { p } ( T _ { p } ) \nabla ^ { 2 } T _ { p } + \frac { \partial \kappa _ { p } ( T _ { p } ) } { \partial T _ { p } } ( \nabla T _ { p } ) ^ { 2 } + G ( T _ { e } - T _ { p } ) \,,$$

$$\text{where } D = - \frac { T _ { e } \sigma \alpha } { e } \frac { 2 \mu _ { 0 } } { 3 N _ { 0 } } \frac { 2 \lambda - 1 } { 2 \lambda } \,, \, S = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \frac { 4 \lambda - 1 } { 4 \lambda ^ { 2 } } \,, \, \text{and } S ^ { 0 } = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \frac { 6 \lambda - 1 } { 8 \lambda ^ { 2 } } \,.$$

The models for the temperature dependence of the transport coefficients are chosen as

$$\text{follows in our calculation: } \sigma = \sigma _ { r n } \frac { T _ { r n } } { T _ { p } } \,, \, \alpha = \alpha _ { r n } \frac { T _ { e } } { T _ { r n } } \,, ^ { 3 2 } \, \kappa _ { e } ^ { 0 } = \kappa _ { e, r n } ^ { 0 } \frac { T _ { e } \, \vdash 1 8 } { T _ { p } } \, \text{ where } \sigma _ { r n } \,, \, \alpha _ { r n } \,, \, \text{ and}$$

0 , rt e κ are the transport coefficients when  both electron temperature and phonon temperature  are  at  room  temperature K .  The  electron  energy  source  term  is 300 = rt T written as 19

$$Q ( { \mathbf r }, t ) = 0. 9 4 ( 1 - R ) I \frac { 1 } { \tau _ { p } \delta } \exp \left [ - \frac { z } { \delta } - 2. 7 7 \frac { t ^ { 2 } } { \tau _ { p } ^ { 2 } } \right ],$$

where R is the reflectivity, I is the fluence of laser pulse, and p τ is the laser pulse time duration. The time zero is noted as the center of laser pulse, 0. = t

The thermal-insulation  boundary  conditions  for  temperature  and  zero  electrical current boundary conditions are used:

$$\nabla T _ { e } \Big | _ { z = 0 } = \nabla T _ { e } \Big | _ { z = L } = 0 \,, \quad \nabla T _ { p } \Big | _ { z = 0 } = \nabla T _ { p } \Big | _ { z = L } = 0 \,, \ \nabla N _ { e } \Big | _ { z = 0 } = \nabla N _ { e } \Big | _ { z = L } = 0 \,.$$

The parameters used in the EHD model calculation in gold thin film are listed in Table I.

## C. Comparison between TTM and EHD Model

By comparing TTM in Eq. (1) with EHD model in Eq. (11), the major differences between the proposed EHD model and the TTM are as follows:

(1). The addition of the balance equation on electron density in the EHD model (Eq. (11a)) enables  us  to  study  the  electron  density  fluctuation  and  transient  electrical  current induced  by  ultra-short  pulse  laser  heating.  The  chemical  potential  gradient µ ∇ is considered as the driving force for  current  flow  in  the  EHD  model,  which  is  not  even considered in the TTM  model.  A  complete energy current density of electron n e Q e E e , , , J J J µ + = is considered in the EHD while n e , J µ is not included in TTM.

- (2).  The  divergence  of  energy  current  density  of  electron  in  the  presence  of  electrical current is 20

$$- \nabla \cdot \mathbf J _ { e, E } = \nabla \cdot [ \mathcal { K } _ { e } ( T _ { e }, T _ { p } ) \nabla T _ { e } ] - \mu \nabla \cdot \mathbf J _ { e, n } - \xi \mathbf J _ { e, n } \cdot \nabla T _ { e } + \frac { 1 } { \sigma } [ e \mathbf J _ { e, n } ] ^ { 2 },$$

where ξ is  the  Thomson  coefficient.  Comparing  to  the  TTM  model,  the  EHD  model considers the electron accumulation, the Thomson effect, and the Joule heat generation through the last three terms on the right side of Eq. (14)

- (3). The divergence of heat current of electrons is

$$\nabla \cdot [ \kappa _ { e } ( T _ { e }, T _ { p } ) \nabla T _ { e } ] = \kappa _ { e } ( T _ { e }, T _ { p } ) \nabla ^ { 2 } T _ { e } + \left [ \frac { \partial \kappa _ { e } ( T _ { e }, T _ { p } ) } { \partial T _ { e } } \right ] \nabla T _ { e } \right ) ^ { 2 } + \left [ \frac { \partial \kappa _ { e } ( T _ { e }, T _ { p } ) } { \partial T _ { p } } \right ] \nabla T _ { p } \cdot \nabla T _ { e } \,, \ ( 1 5 )$$

where the first term on the right side in describes the linear thermal transport due to the electron temperature gradient while the last two terms are additional nonlinear thermal transport  terms.  The  TTM  model  only  explicitly  considers  the  linear  transport  term/ 19 Empirically,  the  nonlinear  terms  can  be  calculated  when  one  has  an  explicit  analytical temperature dependence of ) . Anisimov and Rathfeld , ( p e e T T κ 33 showed that p e e T T / ∝ κ

when F e T T &lt;&lt; and 2 / 5 e e T ∝ κ when F e T T ~ , where TF is  the  Fermi  temperature  of material (usually over 10 K in metals). A more complicated analytical expression of 4 e κ is given by Refs. [34,35]. The question here is whether it is correct to directly substitute the  above  temperature  dependence  of e κ into  Eq.  (15)?  A  revisit  is  needed  from  more basic  approaches  such  as  the  BTE.  We  assume p e e T T / ∝ κ as  an  example  and  set  the Seebeck coefficient 0 which means = α 0  and = D e e S S κ κ = = = 0 ' , then the divergence of heat current in Eq. (11b) becomes

$$\kappa _ { e } ( T _ { e }, T _ { p } ) \nabla ^ { 2 } T _ { e } - \frac { 2 \kappa _ { e } ( T _ { e }, T _ { p } ) } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 }.$$

In contrast, the divergence of heat current in Eq. (1a) in the TTM can be rewritten as:

$$\kappa _ { e } ( T _ { e }, T _ { p } ) \nabla ^ { 2 } T _ { e } + \frac { \kappa _ { e } ( T _ { e }, T _ { p } ) } { T _ { e } } ( \nabla T _ { e } ) ^ { 2 } - \frac { \kappa _ { e } ( T _ { e }, T _ { p } ) } { T _ { p } } \nabla T _ { p } \cdot \nabla T _ { e } \,.$$

The  coefficients  of  linear  term e T 2 ∇ in  Eq.  (16)  and  Eq.  (17)  are  the  same,  and  the coefficient  of  nonlinear  term ( ) 2 e T ∇ in  Eq.  (16)  is e e T / 2 κ -which  is  different  from e e T / κ in Eq. (17), and there is one more term e p p e T T T ∇ ⋅ ∇ -κ in Eq. (17). The difference between the EHD model and the TTM can be neglected only when ( ) e e e T T T 2 2 / ∇ &lt;&lt; ∇ and ( ) e p e p T T T T 2 / ∇ &lt;&lt; ∇ ⋅ ∇ . When (∇ ) e e e T T T 2 2 ~ / ∇ and ( ) e p e p T T T T 2 ~ / ∇ ∇ ⋅ ∇ one have to directly solve the BTE or solve the EHD model instead.

## III. RESULTS AND DISCUSSIONS

In  the  following,  we  compare  the  calculation  results  from  the  coupled  electron  and phonon BTE model, the proposed EHD model and the classical TTM for the coupled electron and phonon transport in gold thin films after ultra-short pulse heating, as show in Fig. 1.

## A. Electron density fluctuation and transient electrical current

We first compare the evolution of electron temperature at the front surface for a 50nmthick  gold  film  after  an  ultra-short  pulse  laser  heating  with  fluence 2 1mJ/cm = I and pulse duration 96fs calculated from the BTE model, the EHD model, and the TTM. = p τ The  EP  coupling  constant  is  chosen  to  be K throughout  this  paper W/m 3 10 1 . 2 16 × = G which is obtained by fitting the experimental data measured by Brorson et al. in Ref. [8], which is close to the values in literature,  2 2 K. W/m 3 10 4 -K W/m 3 10 . 16 16 × × Error! Bookmark not defined. As shown in Fig. 2(a), the electron temperature at the front surface calculated from the EHD model is in good agreement with the prediction from the BTE while the TTM overestimates the electron temperature.

The  large  temperature  gradient  due  to  the  non-equilibrium  ultra-short  pulse  laser heating can induce electron density redistribution and transient electrical current due to electron  diffusion  and  thermoelectric  effect.  Such  electron  density  fluctuation  and transient  electrical  current  after  ultra-short  pulse  laser  heating  can  be  calculated  by  the BTE model and the EHD model while the Fourier-like TTM model is not able to. Figure 2(b) shows the electron density fluctuation in gold film at different time delays,  =0.05ps, t 0.1ps, and 0.2ps. At the beginning, the amplitude of the electron density fluctuation is on the  order  of 10 3 24 /m ,  which  is  about  1%  of  the  electron  density  near  Fermi  surface

0 0 / µ T k N B .    Figure  2(c)  shows  the  particle  current  density  of  electrons  in  gold  film  at different time delays,  =0ps, 0.1ps, 0.2ps. At the beginning, a negative particle current is t noted due to flow of particles towards the front surface as a result of thermoelectric effect. However, the particle  current  density  becomes  positive,  i.e.,  flow  away  from  the  front surface at  =0.1ps, due to the diffusion of electrons driven by electron density gradient as t shown  in  Fig.  2(b).  The  maximum  amplitude  of  the  particle  current  density  is  about 2 30 /sm 10 which  results  an  electrical  current  density  on  the  order  of 10 2 11 A/m .  After t =0.2ps, the particle current density becomes negligible. In our calculation, the Seebeck coefficient of gold is positive 36 which results in a negative Soret coefficient ( 0) in &lt; n S Eq. (9a). This negative Soret coefficient results in the flow of electrons towards the front surface driven by temperature gradient. In other metals with negative Seebeck coefficient and positive Soret coefficient (for example, Al, Pb et al. ), the electrons would flow away from the front surface driven by temperature gradient.

## B. Dependence of temperature evolution on laser pulse duration, laser fluence and film thickness

The  time  scale  that  the  laser-deposited  energy  is  conducted  away  from  the  front surface is on the order of sub-picosecond. When the laser pulse time duration is smaller than this time scale, most energy absorbed from laser is stored at the front surface and the interior of the film is not heated up yet at the beginning of the time evolution. This results in  an  abrupt temperature drop near the front surface. In other words, e T ∇ is  very  large and ( ) e e e T T T 2 2 ~ / ∇ ∇ is  satisfied  which  means  that  the  nonlinear  thermal  transport could be important as we discussed in Sec. IIC. Figure 2(a), Fig. 3(a), and Fig. 3(b) show

the  evolution  of  electron  temperature  at  the  front  surface  of  the  film  with  thickness 50nm for different laser pulse time durations, = L fs , 150fs, and  200fs , when the 96 = p τ laser fluence is 2 1mJ/cm . For all these three time durations, the electron temperature = I calculated  from  the  EHD  model  is  closer  to  that  calculated  from  BTE  while  the  TTM overestimate the electron temperature. The maximum electron temperatures for the BTE and  the  EHD  model  are  almost  the  same  when 96fs .  The  calculated  maximum = p τ electron  temperatures  are  634K  from  the  BTE,  660K  from  the  EHD  model,  and  686K from the TTM when the pulse duration is 150fs. When the duration time increases = p τ to 200fs, the maximum electron temperatures become 621K from the BTE, 653K from the  EHD  model,  and  672K  from  the  TTM.  The  reason  is  that ( ) e e e T T T 2 2 ~ / ∇ ∇ for shorter time duration and the nonlinear thermal transport must be considered.

When the electrons are heated by the laser pulse, they would transport inside the film with the Fermi velocity F v (about 1 4 m/s  in  gold),  and  then  be  reflected  back  by 10 . 6 × the  rear  surface.  These  hot  electrons  return  to  the  front  surface  after  a  time  interval, F v L / 2 ~ . When  this  time  interval  is  smaller  than  the  laser pulse  duration, the temperature evolution at the front surface can be affected by the boundary reflection from the rear surface which results in highly nonlinear transport. Therefore, the thickness of film  is  another  important  parameter  in  our  study.  Figures  3(c)  and  3(d)  show  the evolution  of  electron  temperature  at  front  surface  of  the  film  for  different  thicknesses, 40nm and 60nm , when the laser fluence is = L 2 1mJ/cm  and the laser pulse duration = I is 96fs .  Figure  2(a)  also  shows  the  results  with = p τ 50nm  while  other  conditions = L are the same. For all these three film thickness, the electron temperature calculated from

the EHD model is closer to that from the BTE, and the TTM overestimates the electron temperature.  The  calculated  maximum  electron  temperatures  are  696K  from  the  BTE, 723K from the EHD model, and 750K from the TTM when 40nm. When the film = L thickness increases to 50nm , the electron temperature calculated from the EHD model is almost the same as that from the BTE. When the film thickness is 60nm , the electron temperature calculated from the EHD model is slightly smaller than that from the BTE.

Larger laser fluence results in higher electron temperature at the front surface and the electron  distribution  is  farther  away  from  the  equilibrium  distribution.  Therefore,  the electron transport is highly nonlinear for large laser fluence. Figures 4(a) and 4(b) show the  evolution  of  electron  temperature  at  front  surface  of  the  film  for  different  laser fluences, 2 5mJ/cm  and 10mJ/cm = I 2 when the film thickness is 50nm  and the laser = L pulse  time  duration  is 96fs ,  respectively.  Figure  2(a)  also  shows  the  results  with = p τ 2 1mJ/cm  when  other  conditions  are  the  same.  For = I 2 1mJ/cm  and 5mJ/cm = I 2 ,  the electron temperature calculated from the EHD model is close to that calculated from the BTE for all the time delays considered. The TTM overestimates the electron temperature. For higher laser fluence 2 10mJ/cm , both the EHD model and the TTM underestimate = I the  electron  temperature  in  comparison  with  the  BTE  at  the  beginning  of  the  time evolution  (time  delay&lt;1.5ps).  The  reason  is  that  both  the  EHD  model  and  the  TTM assume that the electron distribution is slightly deviated from the Fermi-Dirac distribution  which  does  not  hold  at  the  beginning  of  the  time  evolution  for  high  laser fluence.  When  time  delay&gt;1.5ps,  the  electron  temperatures  calculated  from  the  EHD model and the BTE become getting closer since the electron-electron scattering drives the

electron sub-system from non-equilibrium distribution to equilibrium one which is called thermalization process. Error! Bookmark not defined.

From Fig. 3 and Fig. 4, we find that the EHD model is more accurate than TTM in thinner  film  and  shorter  laser  pulse  time  duration  because  of  the  contribution  from transient  electrical  current  and  nonlinear  thermal  transport.  When  the  thermalization process  is  important  for  higher  laser  fluence,  both  the  EHD  model  and  the  TTM underestimate the electron temperature at the beginning of the time evolution.

## C. Damage thresholds

We now study the laser damage threshold calculated from the EHD model and the TTM with high laser fluences. We do not show the results from the BTE calculations due to the computational cost. The BTE model calculates the distribution functions in wave vector  and  position  phase-space.  Larger  laser  fluence  requires  larger  truncation  wave vector that significantly increases the computational cost. Figure 5(a) shows that the time evolutions of both electron temperature and phonon temperature at the front surface of the  film  calculated  from  the  EHD  model  slightly  differ  from  that  calculated  from  the TTM  when nm .  Therefore,  the  converged  temperature  (electron  temperature 25 = L becomes equal to phonon temperature) calculated from the EHD model is also slightly smaller than that from TTM. When the film thickness increases to 50nm, the maximum electron temperature calculated from the EHD model could be 800K lower than that from TTM and the converged temperature calculated from EHD model is 100K higher than that from TTM as shown in Fig. 5(b). Overall, the EHD model gives a lower converged

temperature  than  the  TTM,  which  means  that  more  energy  at  the  front  surface  is conducted away from the surface.

When the laser fluence is large enough, the phonon temperature at the front surface can reach the melting point of the metal, which is 1337K for gold. Such laser fluence is called  damage  fluence  threshold.  As  we  shown  in  Fig.  5,  the  converged  temperature depends on the film thickness. Therefore, the damage fluence threshold is also a function of film thickness. Figure 6 shows the calculated damage fluence thresholds from the EHD model and the TTM for an ultra-short pulse laser for different film thickness with a laser pulse  duration 96fs .  As  shown  in  Figs.  5  and  6,  the  TTM  overestimates  the = p τ converged temperature at front surface of the film and the difference of this overestimation  increases  with  the  film  thickness.  Therefore,  the  EHD  model  gives  a larger  damage  fluence  threshold  in  comparison  with  the  TTM.  For  example,  when  the film thickness is smaller than 25nm, the damage fluence thresholds calculated from the EHD  model  and  the  TTM  are  both  around 125mJ/cm 2 .  When  the  film  thickness  is increased  to  60nm,  the  damage  fluence  threshold  calculated  from  the  EHD  model  is 2 mJ/cm  which is much larger than that calculated from the TTM,  236mJ/cm . 375 2

## IV. CONCLUSIONS

In summary, we established an EHD model which consists of the balance equations of  electron  density,  energy  density  of  electrons,  and  energy  density  of  phonons  to investigate the coupled non-equilibrium electron and phonon transport in thin metal films after the ultra-short pulse laser heating. The non-equilibrium transport of electrons such as  the  electron  density  fluctuation  and  transient  electrical  current  are  considered  in  the

EHD model beyond  the  phenomenological  TTM.  By  comparing  the  calculated  results from the BTE for different laser pulse time durations, film thicknesses, and laser fluences, we  found  that the EHD  model  is  more  accurate  than  the  TTM  which  usually overestimates the surface electron temperature when the nonlinear thermal transport in electron transport is important. We find that shorter laser pulse duration and thinner metal film result in stronger nonlinear transport. This model gives a more accurate tool to study the non-equilibrium coupled electron and phonon transport process than the TTM and it is much easier to be solved than the BTE.

Acknowledgements :  J.  Z.  would like to thank Dr. J. L. Cheng and Dr. Y . Y . Wang for valuable discussions. R. Y. would like to acknowledge the support from AFOSR (Grant No. FA9550-11-1-0109) and NSF (Grant No. 0846561). J. Z. and B. L. are supported by the  National  Natural  Science  Foundation  of  China  grant  no.  11334007.  J.  Z.  is  also supported by the program for New Century Excellent Talents in Universities grant no. NCET-13-0431.

## APPENDIX A. ELECTRONIC HEAT CAPACITY

We calculate the temperature-dependent electronic heat capacity by using the Sommerfeld expansion method under the free electron approximation and local equilibrium approximation, the local electron density and energy density can be written as:

$$N = \chi \overset { \infty } { \int } d \varepsilon _ { k } \, \frac { \varepsilon _ { k } ^ { 1 / 2 } } { \exp [ ( \varepsilon _ { k } - \mu ) / k _ { B } T _ { e } + 1 ] } = \chi ( k _ { B } T _ { e } ) ^ { 3 / 2 } \, F _ { 1 / 2 } ( \eta ) \,, \quad \quad ( \text{A1a} )$$

$$U _ { e } = \chi \overset { \infty } { \int } d \varepsilon _ { k } \, \frac { \varepsilon _ { k } ^ { 3 / 2 } } { \exp [ ( \varepsilon _ { k } - \mu ) / k _ { B } T _ { e } + 1 ] } = \chi ( k _ { B } T _ { e } ) ^ { s / 2 } \, F _ { 3 / 2 } ( \eta ) \,, \quad \, \quad ( \text{A1b} )$$

where T k / , ) ( F is the dimensionless integral e B µ η = η j

( ) ∫ ∞ --+ -η = 0 1 ] 1 ) / /[exp( ) ( ε ε ε η e B k j k k j e B j T k d T k F , the electron density of states is

2 / 1 ) ( k k χε ε = Ξ , 2 / 3 2 2 2 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = h m V π χ , V is  the  unit  volume, m is  the  electron  mass,  and

$$j = - 1 / 2, 1 / 2, 3 / 2 \,.$$

The derivatives in Eq. (7) can be calculated as

$$\frac { \partial N } { \partial \mu } = \frac { 1 } { 2 } \, \chi ( k _ { _ { B } } T _ { e } ) ^ { 1 / 2 } \, F _ { - 1 / 2 } ( \eta )$$

$$\frac { \partial N } { \partial T _ { e } } = \frac { 1 } { 2 T _ { e } } \, \chi ( k _ { B } T _ { e } ) ^ { 3 / 2 } \left [ 3 F _ { _ { 1 / 2 } } ( \eta ) - \frac { \mu } { k _ { B } T _ { e } } \, F _ { _ { - 1 / 2 } } ( \eta ) \right ],$$

$$\frac { \partial U _ { e } } { \partial \mu } = \frac { 3 } { 2 } \, \chi ( k _ { _ { B } } T _ { e } ) ^ { 3 / 2 } \, F _ { _ { 1 / 2 } } ( \eta ) \,,$$

$$\frac { \partial U _ { e } } { \partial T _ { e } } = \frac { 1 } { 2 T _ { e } } \, \chi ( k _ { B } T _ { e } ) ^ { 5 / 2 } \left [ 5 F _ { 3 / 2 } ( \eta ) - \frac { \mu } { k _ { B } T _ { e } } 3 F _ { 1 / 2 } ( \eta ) \right ].$$

$$\text{Here we used } \frac { d F _ { 1 / 2 } ( \eta ) } { d \eta } = \frac { 1 } { 2 } F _ { - 1 / 2 } ( \eta ) \, \ \frac { d F _ { 3 / 2 } ( \eta ) } { d \eta } = \frac { 3 } { 2 } F _ { 1 / 2 } ( \eta ) \. \ \text{The Taylor expansion of}$$

)  can be found in Ref. [31] ( η j F

$$F _ { 3 / 2 } ( \eta ) = \frac { 2 } { 5 } \eta ^ { 5 / 2 } \left [ 1 + \frac { 5 \pi ^ { 2 } } { 8 } \eta ^ { - 2 } - \frac { 7 \pi ^ { 4 } } { 3 8 4 } \eta ^ { - 4 } + O ( \eta ^ { - 6 } ) \right ],$$

$$F _ { 1 / 2 } ( \eta ) = \frac { 2 } { 3 } \eta ^ { 3 / 2 } \left [ 1 + \frac { \pi ^ { 2 } } { 8 } \eta ^ { - 2 } + \frac { 7 \pi ^ { 4 } } { 6 4 0 } \eta ^ { - 4 } + O ( \eta ^ { - 6 } ) \right ],$$

$$F _ { - 1 / 2 } ( \eta ) = 2 \eta ^ { 1 / 2 } \left [ 1 - \frac { \pi ^ { 2 } } { 2 4 } \eta ^ { - 2 } - \frac { 7 \pi ^ { 2 } } { 3 8 4 } \eta ^ { - 4 } + O ( \eta ^ { - 6 } ) \right ].$$

By ignoring the high order terms, the electronic heat capacity defined for Eq. (7) can be written as

$$C _ { e } ( T _ { e } ) & = \frac { \partial U _ { e } } { \partial T _ { e } } - \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial T _ { e } } \\ & = \frac { 1 } { 2 T _ { e } } \, \chi ( k _ { b } T _ { e } ) ^ { 5 / 2 } \left [ 5 F _ { 3 / 2 } ( \eta ) - 9 F _ { 1 / 2 } ^ { 2 } ( \eta ) / F _ { - 1 / 2 } ( \eta ) \right ] \\ & \approx \frac { 1 } { 3 } \, \chi \pi ^ { 2 } k _ { b } ^ { 2 } \mu ^ { 1 / 2 } T _ { e } \left [ 1 - \frac { 3 1 \pi ^ { 2 } } { 1 2 0 } \left ( \frac { k _ { b } T _ { e } } { \mu } \right ) ^ { 2 } \right ] \\ & \approx \gamma T _ { e } \left [ 1 - \frac { 3 \pi ^ { 2 } } { 1 0 } \left ( \frac { k _ { b } T _ { e } } { \mu _ { 0 } } \right ) ^ { 2 } \right ] \\ & \quad \lceil \quad \longrightarrow \frac { \pi ^ { 2 } \left ( \left ( \left ( \right ) ^ { 2 } \right ) } { \nu ^ { 2 } } \right ] \quad \left ( \left ( \left ) ^ { 2 } \right ) \right )$$

Here,  we  use  the  relation ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ -≈ 2 0 2 0 12 1 ) ( µ π µ µ e B e T k T and  replace µ e B T k 31 in  the  square bracket by 0 µ e B T k where the higher order terms are ignored. 2 2 0 0 2 / 1 0 2 2 2 3 1 B B k N k π µ µ χπ γ = =

is the ideal heat capacity constant.

## APPENDIX  B.  PARTICLE  CURRENT  DENSITY  AND  ENERGY  CURRENT DENSITY OF ELECTRONS

The  particle  current  density  and  energy  current  density  of  electrons  are  calculated  by solving electron BTE under relaxation time approximation. The distribution function in the BTE can be written into two parts k k k , 1 , 0 f f f + = , 22 where 1 1 / ) ( , 0 + = -µ e B T k k e f ε k is the equilibrium Fermi-Dirac distribution. The non-equilibrium part driven by thermodynamic force along z-direction can be written as:

$$f _ { 1, k } = e E _ { z } \frac { \hbar { k } _ { z } \tau _ { k } } { m } \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } - \frac { \hbar { k } _ { z } \tau _ { k } } { m } \frac { \partial f _ { 0, k } } { \partial z } \,,$$

with k , 1 f under  the  relaxation  time  approximation.  The  particle  current  density  and  the energy current density of electrons alone  -direction are defined as z

$$J _ { e, n, z } = \sum _ { k } v _ { e, k, z } f _ { 1, k }$$

$$J _ { e, E, z } = \sum _ { k } v _ { e, k, z } \mathcal { E } _ { k } f _ { 1, k } \,,$$

We  replace 2 z v by  2 m k 3 / ε using  isotropic  condition,  and  substitute  Eq.  (B1)  into  Eq. (B2).

$$J _ { e, n, z } = & \frac { 2 e } { 3 m } \sum \tau _ { k } \varepsilon _ { k } \left [ \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } \, E _ { z } - \frac { 1 } { e } \frac { \partial f _ { 0, k } } { \partial z } \right ] \\ = & - \frac { N } { e } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle \Big ( e E _ { z } + \frac { \partial \mu } { \partial z } \Big ) - \frac { N } { e T _ { e } } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } ( \varepsilon _ { k } - \mu ) \rangle \frac { \partial T _ { e } } { \partial z },$$

$$J _ { e, E, z } = & \frac { 2 e } { 3 m } \sum \tau _ { k } \varepsilon _ { k } ^ { 2 } \left [ \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } E _ { z } - \frac { 1 } { e } \frac { \partial f _ { 0, k } } { \partial z } \right ] \\ = & - \frac { N } { e } \langle \tau _ { k } \varepsilon _ { k } ^ { s / 2 } \rangle \left ( e E _ { z } + \frac { \partial \mu } { \partial z } \right ) - \frac { N } { e T _ { e } } \langle \tau _ { k } \varepsilon _ { k } ^ { s / 2 } ( \varepsilon _ { k } - \mu ) \rangle \frac { \partial T _ { e } } { \partial z },$$

where the relation ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ -µ + ∂ ∂µ ∂ε ∂ -= ∂ ∂ z T e T e z f z f k k k ε , 0 , 0 is used for the derivation. The bracket

) ( k A ε denotes the integration of arbitrary function ) ( k A ε

$$\left \langle A ( \varepsilon _ { k } ) \right \rangle = - \frac { 2 e } { 3 m N } \chi \overset { \infty } { \int } d \varepsilon _ { k } A ( \varepsilon _ { k } ) \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } \,,$$

The derivative of chemical potential z ∂ ∂µ can be transformed to the derivative of electron density z N ∂ ∂ and electron temperature z T e ∂ ∂ as follows:

$$\frac { \partial \mu } { \partial z } & = \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial z } - \frac { \partial N } { \partial T _ { e } } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial T _ { e } } { \partial z } \\ & = \frac { 2 } { \chi ( k _ { B } T _ { e } ) ^ { 1 / 2 } F _ { - 1 / 2 } ( \eta ) } \frac { \partial N } { \partial z } - \left [ \frac { 3 k _ { B } F _ { 1 / 2 } ( \eta ) } { F _ { - 1 / 2 } ( \eta ) } - \frac { \mu } { T _ { e } } \right ] \frac { \partial T _ { e } } { \partial z }. \\ & \approx \frac { 1 } { \chi \mu ^ { 1 / 2 } } \frac { \partial N } { \partial z } - \frac { \pi ^ { 2 } k _ { B } ^ { 2 } T _ { e } } { 6 \mu } \frac { \partial T _ { e } } { \partial z } \\ \longrightarrow \lambda \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot$$

Then Eqs. (B3a) and (B3b) can be rewritten by substituting Eq. (B5) into them

$$\int n \in \psi \colon \infty \, \infty \, \frac { 1 } { 2 } \, \infty \, \infty \, \frac { 1 } { 3 } \, \infty \, \infty \, \frac { 1 } { 4 } \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \incolon \\ J _ { e, n, z } = & - \frac { N } { e T _ { e } } \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \right \rangle e E _ { z } - \frac { N } { e \chi \mu ^ { 1 / 2 } } \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \right \rangle + \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \left ( \varepsilon _ { k } - \mu \right ) \right \rangle \right ] \frac { \partial T _ { e } } { \partial z } \\ = & - \frac { \sigma } { e } E _ { z } - D _ { n } \frac { \partial N } { \partial z } - S _ { n } \, \frac { \partial T _ { e } } { \partial z }, \\ = & - \frac { \sigma } { e } E _ { z } - D _ { n } \frac { \partial N } { \partial z } - S _ { n } \, \frac { \partial T _ { e } } { \partial z }, \\ = & \quad \varepsilon \quad \varepsilon \quad \varepsilon \quad \varepsilon \, \quad \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \. \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \: \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \numpy \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \inosity \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \inocity \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \infty \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \invisible \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inulative \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inexclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inlike \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inless \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusively \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \begin{matrix} \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inclusive \, \inwikipedia \, \inclusive \, \$$

$$& \varepsilon \quad \varnothing _ { \sigma } \quad \varnothing _ { \sigma } \\ J _ { _ { e, E, \tau } } = & - \frac { N } { e } \langle \tau _ { _ { k } } \varepsilon _ { _ { k } } ^ { 5 / 2 } \rangle e E _ { \tau } - \frac { N } { e \chi \mu ^ { 1 / 2 } } \langle \tau _ { _ { k } } \varepsilon _ { _ { k } } ^ { 5 / 2 } \rangle \frac { \partial N } { \partial z } \\ & - \frac { N } { e T _ { e } } \left [ - \frac { \pi ^ { 2 } ( k _ { _ { B } } T _ { e } ) ^ { 2 } } { 6 \mu } \langle \tau _ { _ { k } } \varepsilon _ { _ { k } } ^ { 5 / 2 } \rangle + \langle \tau _ { _ { k } } \varepsilon _ { _ { k } } ^ { 5 / 2 } ( \varepsilon _ { _ { k } } - \mu ) \rangle \right ] \frac { \partial T _ { e } } { \partial z } \\ = & - \left ( - T _ { e } \sigma \alpha + \frac { \mu \sigma } { e } \right ) E _ { \tau } - D _ { E } \frac { \partial N } { \partial z } - S _ { E } \frac { \partial T _ { e } } { \partial z }.$$

The transport coefficients are defined as follows: the electrical conductivity is 2 / 3 k k Ne ε τ σ = ,  the  diffusion  constant  of  electron 0 0 2 2 / 1 2 3 2 1 N e e Dn µ σ χµ σ ≈ = ,  the  Soret

$$\text{coefficient is } S _ { n } = - \frac { \sigma \alpha } { e } \frac { 2 \lambda - 1 } { 2 \lambda } \,, \text{ the Seebeck coefficient is } \alpha = - \frac { 1 } { e T _ { e } } \frac { \left \langle \tau _ { t } \varepsilon _ { k } ^ { 3 / 2 } \left ( \varepsilon _ { k } - \mu \right ) \right \rangle } { \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \right \rangle } \,,$$

$$D _ { E } = - \frac { T _ { e } \sigma \alpha } { e } \frac { 2 \mu _ { 0 } } { 3 N _ { 0 } } + \mu D _ { n } \ \, \ \text{ and } \ S _ { E } = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \frac { 1 } { 2 \lambda } + \mu S _ { n } \ \. \ \text{ The                                                                                                                                                                                                         thermal }$$

conductivity is 2 0 σα κ κ e e e T -= where ( ) 2 2 / 3 0 µ ε ε τ κ -= k k k e e eT N . The relations

λ µ π α = 0 2 2 3 e T e k B , 2 / 3 0 0 3 2 χµ = N , and 0 µ µ ≈ are  also  used 31 where λ is a  material parameter. 32 Mathematical tricks

$$\text{ics} \quad \left \langle \tau _ { k } \varepsilon _ { k } ^ { 5 / 2 } \right \rangle = \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \left ( \varepsilon _ { k } - \mu \right ) \right \rangle + \mu \left \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \right \rangle \quad \text{and}$$

$$\left \langle \tau _ { k } \mathcal { E } _ { k } ^ { 5 / 2 } ( \mathcal { E } _ { k } - \mu ) \right \rangle = \left \langle \tau _ { k } \mathcal { E } _ { k } ^ { 3 / 2 } ( \mathcal { E } _ { k } - \mu ) ^ { 2 } \right \rangle + \mu \left \langle \tau _ { k } \mathcal { E } _ { k } ^ { 3 / 2 } ( \mathcal { E } _ { k } - \mu ) \right \rangle \text{ are use} \right \rangle$$

are used in the calculation.

The current components in x and y direction can be calculated in the same way. We can then extend Eqs. (B6a) and (B6b) to three dimension equations

$$\mathbf J _ { e, n } = - \frac { \sigma } { e } \mathbf E - D _ { n } \nabla N - S _ { n } \nabla T _ { e } \,,$$

$$\mathbf J _ { e, E } = - \left ( - T _ { e } \sigma \alpha + \frac { \mu \sigma } { e } \right ) \mathbf E - D _ { E } \nabla N - S _ { E } \nabla T _ { e } \,.$$

## APPENDIX  C.  DIVERGENCES  OF  ELECTRON  DENSITY  AND  ENERGY DENSITY OF ELECTRONS

The divergences of electron density and energy density of electrons can be calculated as:

$$- \nabla \cdot \mathbf J _ { e, n } = - \frac$$

$$- \nabla \cdot \mathbf J _ { e, E } = - \frac$$

when  the  thermodynamic  forces  such  as  temperature  gradient  and  chemical  potential gradient are alone z -direction, i.e. 0 , , , , = ∂ ∂ -= ∂ ∂ -y n e x n e J y J x and

$$- \frac { \partial } { \partial x } J _ { e, E, x } =$$

We assume that only the equilibrium distribution function in k , 1 f vary with position and the external electric field is absent. Eqs. (C1a) and (C1b) become

$$\sum _ { \substack { w a t i n d e s t i n d e s$$

$$\begin{array} {$$

Here the second order derivative of distribution function is

$$\frac { \partial ^ { 2 } f _ { 0, k } } { \partial z ^ { 2 } } = & - \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } \left [ \frac { 1 } { k _ { 2 } T _ { e } } \frac { e ^ { ( e _ { k } + \mu ) / k _ { 2 } T _ { e } } } { e ^ { ( e _ { k } - \mu ) / k _ { 2 } T _ { e } } } + 1 } \right ] \left ( \frac { \partial \mu } { \partial z } + \frac { \varepsilon _ { k } - \mu } { T _ { e } } \frac { \partial T _ { e } } { \partial z } \right ) ^ { 2 } \\ & + \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } \frac { 2 } { T _ { e } } \frac { \partial T _ { e } } { \partial z } \left ( \frac { \partial \mu } { \partial z } + \frac { \varepsilon _ { k } - \mu } { T _ { e } } \frac { \partial T _ { e } } { \partial z } \right ) - \frac { \partial f _ { 0, k } } { \partial \varepsilon _ { k } } \left ( \frac { \partial ^ { 2 } \mu } { \partial z ^ { 2 } } + \frac { \varepsilon _ { k } - \mu } { T _ { e } } \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } \right )$$

The  first  term  on  the  right  in  Eq.  (C3)  is  symmetric  on  the  chemical  potential,  its contribution to z n e J z , , ∂ ∂ -and z E e J z , , ∂ ∂ -is  negligible.  That  is  why  we  use  ' ≈ '  in  Eqs.

(C2a) and (C2b). 2 2 z ∂ ∂ µ in Eqs. (C2a) and (C2b) can be transformed to

$$\frac { \partial ^ { 2 } \mu } { \partial z ^ { 2 } } = \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \left [ \frac { \partial ^ { 2 } N } { \partial z ^ { 2 } } - \frac { \partial N } { \partial T _ { e } } \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } - \frac { \partial } { \partial z } \left ( \frac { \partial N } { \partial \mu } \right ) \frac { \partial \mu } { \partial z } - \frac { \partial } { \partial z } \left ( \frac { \partial N } { \partial T _ { e } } \right ) \frac { \partial T _ { e } } { \partial z } \right ].$$

And z ∂ ∂µ in  Eqs.  (C2a),  (C2b),  and  (C4)  can  be  transformed  to z N ∂ ∂ and z T e ∂ ∂ using  Eq.

(B5). Then Eqs. (C2a) and (C2b) are rewritten as

$$( B 5 ). \text{ Then Eqs. } ( C 2 a ) \text{ and } ( C 2 b ) \text{ are rewritten as } \\ - \frac { \partial } { \partial z } J _ { e, n, z } \left ( \frac { N } { \exp } \frac { N } { \exp } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle \frac { \partial ^ { 2 } N } { \partial z ^ { 2 } } + \frac { N } { \exp } \left [ - \frac { \pi ^ { 2 } ( k _ { b } T _ { e } ) ^ { 2 } } { 6 \mu } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle + \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } ( \varepsilon _ { - } \mu ) \rangle \right ] \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } \\ & - \frac { N } { 2 e \chi ^ { 2 } \mu ^ { 2 } } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle \left ( \frac { \partial N } { \partial z } \right ) ^ { 2 } - \left [ \frac { 2 N } { \exp } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } ( \varepsilon _ { - } \mu ) \rangle - \frac { \pi ^ { 2 } k _ { b } } { 6 \mu } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle \right ] \left ( \frac { \partial T _ { e } } { \partial z } \right ) ^ { 2 } \\ & - \frac { 2 N } { e T _ { e } \chi \mu ^ { 1 / 2 } } \langle \tau _ { k } \varepsilon _ { k } ^ { 3 / 2 } \rangle \frac { \partial N } { \partial z } \frac { \partial T _ { e } } { \partial z } \\ & = D _ { n } \frac { \partial ^ { 2 } N } { \partial z ^ { 2 } } + S _ { n } ^ { \ } \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } - \frac { D _ { n } } { 3 N _ { 0 } } \left ( \frac { \partial N } { \partial z } \right ) ^ { 2 } - \frac { 2 S ^ { \prime } _ { n } } { T _ { e } } \left ( \frac { \partial T _ { e } } { \partial z } \right ) ^ { 2 } - \frac { 2 D _ { n } } { T _ { e } } \frac { \partial N } { \partial z } \frac { \partial T _ { e } } { \partial z }, \\ \Delta \quad \lambda ^ { \prime } \quad \Delta ^ { 2 } \lambda ^ { \prime } \quad \lambda ^ { \prime } \Gamma \quad \mu ^ { 2 } ( \L _ { k } \ \ T \ \L ^ { 2 } \. \.$$

$$\Lambda ^ { \prime } \, \dot { \ } d z ^ { \prime } \quad \Lambda ^ { \prime } \, \dot { \ } d z ^ { \prime } \quad \Lambda ^ { \prime } \quad T _ { e } \, \left \{ \, \dot { \ } d z \, \right \} \quad T _ { e } \, \left \{ \, \dot { \ } d z \, \right \} \quad T _ { e } \, \left \{ \, \dot { \ } d z \, \right \} \, \\ - \frac { \partial } { \partial } J _ { e, E, \varepsilon } \approx \frac { N } { e \chi \mu ^ { \prime 2 } } \left \langle \tau _ { e } \varepsilon _ { k } ^ { s / 2 } \right \rangle \frac { \partial ^ { 2 } N } { \partial z ^ { 2 } } + \frac { N } { e T _ { e } } \left [ - \frac { \pi ^ { 2 } ( k _ { B } T _ { e } ) ^ { 2 } } { 6 \mu } \left \langle \tau _ { e } \varepsilon _ { k } ^ { s / 2 } \right \rangle + \left \langle \tau _ { e } ^ { s / 2 } _ { k } ( e _ { k } - \mu ) \right \rangle \right ] \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } \\ - \frac { N } { 2 e \chi ^ { 2 } \mu ^ { 2 } } \left \langle \tau _ { e } ^ { s / 2 } _ { k } \right \rangle \left ( \frac { \partial N } { \partial z } \right ) ^ { 2 } - \left [ \frac { 2 N } { e T _ { e } ^ { 2 } } \left \langle \tau _ { e } ^ { s / 2 } _ { k } ( e _ { k } - \mu ) \right ) - \frac { \pi ^ { 2 } k _ { B } ^ { 2 } } { 6 \mu } \left \langle \tau _ { e } ^ { s / 2 } _ { k } \right \rangle \right ) \left ( \frac { \partial T _ { e } } { \partial z } \right ) ^ { 2 } \\ - \frac { 2 N } { e T _ { \chi } \mu ^ { \prime 2 } } \left \langle \tau _ { e } ^ { s / 2 } _ { k } \right \rangle \frac { \partial N } { \partial z } \frac { \partial T _ { e } } { \partial z } \\ = & \, D _ { E } \frac { \partial ^ { 2 } N } { \partial z ^ { 2 } } + S _ { E } \frac { \partial ^ { 2 } T _ { e } } { \partial z ^ { 2 } } \frac { D _ { E } } { 3 N _ { 0 } } \left ( \frac { \partial N } { \partial z } \right ) ^ { 2 } - \frac { 2 S _ { E } } { T _ { e } } \left ( \frac { \partial T _ { e } } { \partial z } \right ) ^ { 2 } - \frac { 2 D _ { E } } { T _ { e } } \frac { \partial N } { \partial z } \frac { \partial T _ { e } } { \partial z }, \\ \\ \tau \alpha \alpha \, 4 \lambda - 1 \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime } \quad \Lambda ^ { \prime }$$

$$\text{where } S _ { \, _ { n } } ^ { \prime } = - \frac { \sigma \alpha } { e } \frac { 4 \lambda - 1 } { 4 \lambda } \text{ and } S _ { \, _ { E } } ^ { \prime } = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \frac { 1 } { 4 \lambda } + \mu S _ { \, _ { n } } ^ { \prime } \,.$$

Eqs. (C5a) and (C5b) can be extended to three dimension form

$$- \nabla \cdot \mathbf J _ { e, n } \approx D _ { n } \nabla ^ { 2 } N + S _ { n } \nabla ^ { 2 } T _ { e } - \frac { D _ { n } } { 3 N _ { 0 } } ( \nabla N ) ^ { 2 } - \frac { 2 S _ { n } ^ { \prime } } { T _ { e } } ( \nabla T _ { e } ) ^ { 2 } - \frac { 2 D _ { n } } { T _ { e } } ( \nabla N \cdot \nabla T ) \,,$$

$$- \nabla \cdot \mathbf J _ { e, E } \approx D _ { E } \nabla ^ { 2 } N + S _ { E } \nabla ^ { 2 } T _ { e } - \frac { D _ { E } } { 3 N _ { 0 } } ( \nabla N ) ^ { 2 } - \frac { 2 S _ { E } } { T _ { e } } ( \nabla T _ { e } ) ^ { 2 } - \frac { 2 D _ { E } } { T _ { e } } ( \nabla N \cdot \nabla T _ { e } ).$$

## APPENDIX D. BALANCE EQUATION OF ENERGY DENSITY OF ELECTRONS

$$\text{Using the relation in } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text
 \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{  the balance}$$

equation Eq. (9b) in the absence of external electric field becomes

$$\frac { \partial U _ { e } } { \partial t } = \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial t } + C _ { e } ( T _ { e } ) \, \frac { \partial T _ { e } } { \partial t } = - \nabla \cdot \mathbf J _ { e, E } \, - \, G ( T _ { e } \, - \, T _ { p } ) \, + \, Q ( \mathbf r, t ) \,.$$

Eq. (D1) can be rewritten as

$$C _ { e } ( T _ { e } ) \frac { \partial T _ { e } } { \partial t } = - \nabla \cdot \mathbf J _ { e, E } - \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial t } - G ( T _ { e } - T _ { p } ) + Q ( \mathbf r, t ) \,.$$

Using the divergences of electron density and energy density of electrons obtained in Eqs.

(C6a) and (C6b), the first two terms on the right side of Eq. (D2) are

$$- \nabla \cdot \mathbf J _ { e, E } - \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } \frac { \partial N } { \partial t } \approx - \nabla \cdot \mathbf J _ { e, E } + \mu \left [ 1 + \frac { \pi ^ { 2 } \left ( k _ { b } T _ { e } \right ) ^ { 2 } } { 6 \mu ^ { 2 } } \right ] \nabla \cdot \mathbf J _ { e, n } \\ = D \nabla ^ { 2 } N + S \nabla ^ { 2 } T _ { e } - \frac { D } { 3 N _ { 0 } } \left ( \nabla N \right ) ^ { 2 } - \frac { 2 S ^ { \prime } } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 } - \frac { 2 D } { T _ { e } } \left ( \nabla N \cdot \nabla T _ { e } \right ),$$

$$\text{where } \ D = - \frac { T _ { e } \sigma \alpha } { e } \frac { 2 \mu _ { 0 } } { 3 N _ { 0 } } \frac { 2 \lambda - 1 } { 2 \lambda } \,, \, S = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \, \frac { 4 \lambda - 1 } { 4 \lambda ^ { 2 } } \,, \, \text{and } S ^ { \prime } = \kappa _ { e } ^ { 0 } - T _ { e } \sigma \alpha ^ { 2 } \, \frac { 6 \lambda - 1 } { 8 \lambda ^ { 2 } } \,. \text{ The}$$

$$\text{relations} \ \frac { \partial N } { \partial t } = - \nabla \cdot \mathbf J _ { e, \eta } \text{ and } \ \frac { \partial U _ { e } } { \partial \mu } \left ( \frac { \partial N } { \partial \mu } \right ) ^ { - 1 } = 3 k _ { _ { B } } T _ { e } \frac { F _ { 1 / 2 } ( \eta ) } { F _ { - 1 / 2 } ( \eta ) } \approx \mu \left [ 1 + \frac { \pi ^ { 2 } \left ( k _ { _ { B } } T _ { e } \right ) ^ { 2 } } { 6 \mu ^ { 2 } } \right ] \text{ are } \text{ used.}$$

Combining Eqs. (D2) and (D3), the balance equation of energy density of electrons is finally written as:

$$C _ { e } ( T _ { e } ) \frac { \partial T _ { e } } { \partial t } & = D \nabla ^ { 2 } N + S \nabla ^ { 2 } T _ { e } - \frac { D } { 3 N _ { 0 } } \left ( \nabla N \right ) ^ { 2 } - \frac { 2 S } { T _ { e } } \left ( \nabla T _ { e } \right ) ^ { 2 } - \frac { 2 D } { T _ { e } } \left ( \nabla N \cdot \nabla T _ { e } \right ) \\ & \quad - G ( T _ { e } - T _ { p } ) + Q ( { \mathbf r }, t ).$$

Table I. Parameters used in the calculation of the BTE model and the EHD model.

| R a          | 0.93                | δ a    | 15.3 nm              |
|--------------|---------------------|--------|----------------------|
| rt σ b       | m / 10 52 . 4 7 Ω × | rt α c | μ V/K 94 . 1         |
| 0 , rt e κ a | W/mK 315            | λ c    | -1.48                |
| 0 µ b        | 5.5 eV              | 0 N b  | 3 28 m / 10 86 . 5 × |
| γ a          | 2 3 K J/m 9 . 62    | p C a  | K J/m 70 3           |
| p v b        | 3240m/s             | ee τ d | 0.04ps               |
| pp τ d       | 0.08ps              |        |                      |

a Ref. [19];   [28];   [32];    Ref. [37] b c d

## Figure Captions

FIG.  1  (Color  online)  Schematic  diagram  of  the  coupled  non-equilibrium  electron  and phonon transport in metal film after ultra-short pulse laser heating. The film thickness is L and the surface of the metal film is in x-y plane. The low energy electrons (black dot) near  the  front  surface  are  excited  to  be  higher  energy  ones  (red  dot)  after  absorbing  a photon.  The  absorption  intensity  exponentially  decays  with  the  radiation  penetration depth δ .  The  hot  electrons  release  energy  through  two  channels:  transport  to  the  rear surface  and  colliding  with  the  lattice  to  emit  phonons,  i.e.  EP  scattering.  When  the electrons  reach  both  the  front  and  rear  surfaces,  they  are  scattered  back  randomly  if  a rough  surface  is  assumed.  In  this  work,  diffusive  reflection  boundary  conditions  of electrons and phonons are used.

FIG. 2 (Color online) (a) Time-dependent electron temperature at the front surface of gold film calculated from the BTE, the EHD model, and theTTM. Evolution of (b) electron density fluctuation and (c) particle current density of electrons in the film calculated from the BTE and the EHD models at different delay time, respectively. The calculations use a laser fluence 2 mJ/cm  and film thickness 1 = I nm  as input parameters. 50 = L

FIG. 3 (Color online) Time-dependent electron temperature at the front surface of the film calculated  from  the  BTE,  the  EHD  model,  and  the  TTM  after  ultra-short  pulse  laser heating with different thicknesses and laser durations of (a) nm  and 50 = L fs , (b) 150 = p τ nm  and 50 = L fs ,  (c) 200 = p τ nm and 40 = L fs ,  (d) 96 = p τ nm and 60 = L fs . 96 = p τ The calculations use laser fluence 2 mJ/cm 1 = I as input parameter.

FIG. 4 (Color online) Time-dependent electron temperature at front surface of the film calculated from the BTE, the EHD model, and the TTM for different laser fluences of 2 mJ/cm 5 and  10mJ/cm 2 . The  calculations  use  laser  duration fs  and  film 96 = p τ thickness nm as input parameters. 50 = L

FIG.  5  (Color  online)  Time-dependent  electron  and  phonon  temperature  at  the  front surface calculated from the EHD model and the TTM for different film thicknesses of 25nm and 50nm for a high laser fluence 2 mJ/cm  with a laser duration 100 = I fs . 96 = p τ

FIG. 6 (Color online) Comparison of the damage threshold for a fs  laser pulse as 96 = p τ a function of film thickness, calculated from the EHD model and the TTM.

Figure 1

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-048e27dd-56ae-4f67-8798-d559bd7ccbdb.jpg)

FIG.  1  (Color  online)  Schematic  diagram  of  the  coupled  non-equilibrium  electron  and phonon transport in metal film after ultra-short pulse laser heating. The film thickness is L and the surface of the metal film is in x-y plane. The low energy electrons (black dot) near  the  front  surface  are  excited  to  be  higher  energy  ones  (red  dot)  after  absorbing  a photon.  The  absorption  intensity  exponentially  decays  with  the  radiation  penetration depth δ .  The  hot  electrons  release  energy  through  two  channels:  transport  to  the  rear surface  and  colliding  with  the  lattice  to  emit  phonons,  i.e.  EP  scattering.  When  the electrons  reach  both  the  front  and  rear  surfaces,  they  are  scattered  back  randomly  if  a rough  surface  is  assumed.  In  this  work,  diffusive  reflection  boundary  conditions  of electrons and phonons are used.

Figure 2

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-2b0b9ff3-e5a3-495a-b677-becb3ec621a7.jpg)

FIG. 2 (Color online) (a) Time-dependent electron temperature at the front surface of gold film calculated from the BTE, the EHD model, and the TTM. Evolution of (b) electron density fluctuation and (c) particle current density of electrons in the film calculated from the BTE and the EHD models at different delay time, respectively. The calculations use a laser fluence 2 mJ/cm  and film thickness 1 = I nm  as input parameters. 50 = L

Figure 3

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-76269a7f-b4be-47d1-8528-ca23303c907a.jpg)

FIG. 3 (Color online) Time-dependent electron temperature at the front surface of the film calculated  from  the  BTE,  the  EHD  model,  and  the  TTM  after  ultra-short  pulse  laser heating with different thicknesses and laser durations of (a) nm  and 50 = L fs , (b) 150 = p τ nm  and 50 = L fs ,  (c) 200 = p τ nm and 40 = L fs ,  (d) 96 = p τ nm and 60 = L fs . 96 = p τ The calculations use laser fluence 2 mJ/cm 1 = I as input parameter.

Figure 4

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-8c1762c9-12aa-4e45-aef9-943c50d530a9.jpg)

FIG. 4 (Color online) Time-dependent electron temperature at front surface of the film calculated from the BTE, the EHD model, and the TTM for different laser fluences of 2 mJ/cm 5 and  10mJ/cm 2 . The  calculations  use  laser  duration fs  and  film 96 = p τ thickness nm as input parameters. 50 = L

Figure 5

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-23c43451-da9d-416a-9347-d2b4007e87f4.jpg)

FIG. 5 (Color online) Comparison of time-dependent electron and phonon temperature at the  front  surface  calculated  from  the  EHD  model  and  the  TTM  for  different  film thicknesses  of  25nm  and  50nm  for  a  high  laser  fluence 2 mJ/cm 100 = I with  a  laser duration fs . 96 = p τ

Figure 6

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0426v1-figure-78af9d86-a750-42ca-9700-d3fff6117718.jpg)

FIG. 6 (Color online) Comparison of the damage threshold for a fs  laser pulse as 96 = p τ a function of film thickness, calculated from the EHD model and the TTM.

## REFERENCES

| 1 P. B. Corkum, F. Brunel, N. K. Sherman, and T. Srinivasan-Rao, Phys. Rev. Lett. 61 ,      |
|---------------------------------------------------------------------------------------------|
| 2886 (1988).                                                                                |
| 2 J. G. Fujimoto, J. M. Liu, E. P. Ippen, and N. Bloembergen, Phys. Rev. Lett. 53 , 1837    |
| (1984).                                                                                     |
| 3 H. E. Elsayed-Ali, T. B. Norris, M.A. Pessot, and G. A. Mourou, Phys. Rev. Lett. 58 ,     |
| 1212 (1986).                                                                                |
| 4 G. L. Eesley, Phys. Rev. Lett. 51 , 2140 (1983).                                          |
| 5 C. Suarez, W. E. Bron, and T. Juhasz, Phys. Rev. Lett. 75 , 4536 (1995).                  |
| 6 W. S. Fann, R. Storz, H. W. K. Tom, and J. Bokor, Phys. Rev. Lett. 68 , 2834 (1992).      |
| 7 R. W. Schoenlein, W. Z. Lin, J. G. Fujimoto, and G. L. Eesley, Phys. Rev. Lett. 58 , 1680 |
| (1987).                                                                                     |
| 8 S. D. Brorson, J. G. Fujimoto, and E. P. Ippen, Phys. Rev. Lett. 59 , 1962 (1987).        |
| 9 P. B. Allen, Phys. Rev. Lett. 59 , 1460 (1987).                                           |
| 10 S. D. Brorson, A. Kazeroonian, J. S. Moodera, D. W. Face, T. K. Cheng, E. P. Ippen, M.   |
| S. Dresselhaus, and G. Dresselhaus, Phys. Rev. Lett. 64 , 2172 (1990).                      |
| 11 M. Bass, Laser Materials Processing (North-Holland, Amsterdam, 1993).                    |
| 12 H. Häkkinen and U. Landman, Phys. Rev. Lett. 71 , 1023 (1993).                           |
| 13 C. L. Cleveland, U. Landman, and R. N. Barnett, Phys. Rev. Lett. 49 , 790 (1982).        |
| 14 C. Schäfer, H. M. Urbassek, and L. V. Zhigilei, Phys. Rev. B 66 , 115404 (2002).         |
| 15 K. Sokolowski-Tinten, J. Bialkowski, A. Cavalleri, D. von der Linder, A. Opatin, J.      |
| Meyer-ter-Vehn, and S. I. Anisimov, Phys. Rev. Lett. 81 , 224 (1998).                       |

| 16 V. Schmidt, W. Husinsky, and G. Zetz, Phys. Rev. Lett. 85 , 3516 (2000).                 |
|---------------------------------------------------------------------------------------------|
| 17 M. I. Kaganov, I. M. Lifshits, and L. V. Tanatarov, Sov. Phys. JETP 4 , 173 (1957).      |
| 18 S. I. Anisimov, B. L. Kapeliovich, and T. L. Perel'man, Sov. Phys. JETP 39 , 375 (1974). |
| 19 T. Q. Qiu and C. L. Tien, Int. J. Heat Mass Transfer 35 , 719 (1992).                    |
| 20 Z. C. Wang, Thermodynamics and Statistical Physics , (Higher Education Press, Beijing,   |
| 1993)                                                                                       |
| 21 J. K. Chen, D. Y. Tzou, and J. E. Beraun, Int. J. Heat Mass Transfer, 49, 307 (2006).    |
| 22 M. Lundstrom, Fundmentals of Carrier Transport (Cambridge University Press,              |
| Cambridge, 2000).                                                                           |
| 23 The phonon branch index is absorbed into the wave vector for simplification.             |
| 24 L. D. Pietanza, G. Colonna, S. Longo, and M. Capitelli, Eur. Phys. J. D 45 , 369 (2007). |
| 25 W. L. McMillan, Phys. Rev. 167 ,331 (1968).                                              |
| 26 J. L. Cheng and M. W. Wu, J. Appl. Phys. 101 , 073702 (2007).                            |
| 27 R. Bauer, A. Schmid, P. Pavone, and D. Strauch, Phys. Rev. B 57 , 11276 (1998).          |
| 28 H. R. Schober, P. H. Dederichs, Au , Landolt-Börnstein Vol.III/13A, edited by K.-H.      |
| Hellwege, J. L. Olsen (Springer-Verlag, Berlin, 1981).                                      |
| 29 J. M. Ziman, Electrons and Phonons , (Oxford University Press, London,1962, p. 456)      |
| 30 L. D. Pietanza, G. Colonna, and M. Capitelli, AIP Conf. Proc. 762 , 1241 (2004).         |
| 31 N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston,        |
| New York, 1976).                                                                            |
| 32 N. F. Mott and H. Jones, The Theory of the Properties of Metals and Alloys , (Dover      |
| Publications, Inc., New York, 1958).                                                        |

- 33 S. I. Anisimov and B. Rethfeld, Proc. SPIE 3093 , Nonresonant Laser-Matter Interaction (NLMI-9) , 192 (1997).
- 34 D. G. Yakovlev and V. A. Urpin, Sov. Astron. J., 24 , 303 (1980).
- 35 Yu. V. Martynenko and Yu. N. Yavlinskii, Sov. Phys. Doklady, 28 , 391 (1983).
- 36 D. K. C. MacDonald, Thermoelectricity: an Introduction to the Principles , (Dover Publication, Inc., New York, 2006)
- 37 A. Pattamata and C. K. Madnia, J. Heat Trans. 131 , 082401 (2009).
Jun Zhou Nianbei Li Ronggui Yang
An Electrohydrodynamics Model for Non-equilibrium Electron and Phonon Transport in Metal Films after Ultra-short Pulse Laser Heating
The electrons and phonons in metal films after ultra-short pulse laser
heating are in highly non-equilibrium states not only between the electron
sub-system and the phonon sub-system but also within the electron sub-system.
An electrohydrodynamics model consisting of the balance equations of electron
density, energy density of electrons, and energy density of phonons is derived
from the coupled non-equilibrium electron and phonon Boltzmann transport
equations to study the nonlinear transport phenomena, such as the electron
density fluctuation and the transient electrical current in metal films, after
ultra-short pulse laser heating. The time-dependent temperature distributions
is calculated by the coupled electron and phonon Boltzmann transport equations,
the electrohydrodynamics model derived in this work, and the two-temperature
model for different laser pulse durations, film thicknesses, and laser
fluences. We find that the two-temperature model overestimates the electron
temperature at the frontsurface of the film and underestimates the damage
threshold when the nonlinear thermal transport of electrons is important. The
electrohydrodynamics model proposedin this work could be a more accurate
prediction tool to study the non-equilibrium electron phonon transport process
than the two-temperature model and it is much easier to be solved than the
coupled electron and phonon Boltzmann transport equations.