Wigner Equations for Phonons Transport and Quantum Heat Flux

Starting from the quantum Liouville equation for the density operator and applying the Weyl quantization, Wigner equations for the acoustic, optical and Z phonons are deduced. The equations are valid for any solid, including 2D crystals like graphene. With the use of Moyal’s calculus and its properties, the pseudo-differential operators are expanded up to the second order in ħ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbar $$\end{document}. An energy transport model is obtained by using the moment method with closure relations based on a quantum version of the Maximum Entropy Principle by employing a relaxation time approximation for the production terms of energy and energy flux. An explicit form of the thermal conductivity with quantum correction up to ħ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbar ^2$$\end{document} order is obtained under a long-time scaling for the most relevant phonon branches.


Introduction
The use of the Wigner function is one of the most promising ways to study quantum transport.Its main advantage is that a description similar to the classical or semiclassical transport is obtained in a suitable phase-space.The mean values are expectation values with respect to the Wigner function as if the latter were a probability density and the semiclassical limit of the Wigner transport equation recovers, at least formally, the Boltzmann transport one.There is a huge body of literature regarding the Wigner equation and the way to numerically solve it (see for example [2,3,4] and references therein).However, the most of the works on the subject consider a quadratic dispersion relation for the energy.Instead, for several materials like semiconductors or semimetal, e.g.graphene, other dispersion relations must be considered [5,6,7].From the Wigner transport equation quantum hydrodynamical models have been obtained in [8] for charge transport in silicon in the case of parabolic bands, while in [9] the same has been devised for electrons moving in graphene.
The enhanced miniaturization of electron and mechanical devices makes the thermal effects increasingly relevant [10,11] requiring the use of physically accurate models.At kinetic level a good description is that based on the semiclassical Peierls-Boltzmann equation for each phonon branch.However, for typical lengths smaller than the phonon mean-free path also quantum effects must be considered (see for example [11]).The Wigner equation is a natural approach that better reveal the wave nature of phonons in such circumstances, gives the Peierls-Boltzmann equation as semiclassical limit and still keeps the structure of a kinetic formulation.In this work, the focus is on the acoustic and optical phonons dynamics with a general dispersion relation.
In order to get insights into the quantum corrections, moment equations are deduced from the corresponding Wigner equation.As in the classical case, one is led to a system of balance equations that are not closed.So, the well-known problem of getting closure relations arises, that is the issue to express the additional fields appearing in the moment equations in terms of a set of fundamental variables, e.g. the phonon energy density and energy flux.A sound way to accomplish this task is resorting to a quantum formulation of the maximum entropy principle [12] (hereafter QMEP), formulated for the first time by Jaynes [13].Recently, a more formal theory has been developed in a series of papers [14,15] with several applications, for example for charge transport in semiconductors [8,16,17,18].The interested reader is also referred to [19].
We apply QMEP to the Wigner equations assuming the energy density and the energy flux for each species of phonons as basic fields.By expanding up to the second order in , quantum corrections to the semiclassical case [1] are deduced.In particular, in a long time scaling an asymptotic expression for the heat flux is obtained.The latter consists of a Fourier-like part with a highly nonlinear second order correction in the temperature gradient.Explicit formulas for acoustic phonons in the Debye approximation are written.
The plan of the paper is as follows.In section 2, the semiclassical phonon transport is summarized while in section 3 we write down the Wigner equations for acoustic and optical phonons.Section 4 is dedicated to deducing the moment equations whose closure relations are achieved by QMEP in section 5.In the last section a definition of local temperature is introduced and an asymptotic expression of the quantum correction to the heat flux is drawn.

Semiclassical phonon transport
In a crystal lattice the transport of energy is quantized in terms of quasi-particles named phonons which are present with several branches and propagation modes.The latters vary from a material to another but in any case they are grouped in acoustic and optical phonon branches which, in turn, can oscillate in the longitudinal or transversal direction.The complete dispersion relations can be usually obtained by a numerical approach in the first Brillouin zone (FBZ) B. However, in the applications some standard approximations are often adopted.
For the acoustic phonons, the Debye approximation for the dispersion relation ε s (q) is usually assumed, ε s (q) = ω s (q) = c s |q|, s = LA, T A. LA stands for longitudinal acoustic while T A for transversal acoustic.c s is the sound speed of the s-branch and denotes the reduced Planck constant.Consistently, the first Brillouin zone is extended to R d .Here d is the dimension of the space; d = 3 for bulk crystal while d = 2 for graphene or similar 2D material like dichalcogenides.
For the longitudinal optical (LO) and the transversal optical (TO) phonon, the Einstein dispersion relation, ω s ≈ const, with s = LO, T O, the phonon angular frequency, is usually adopted.
Note that under such an assumption, the group velocity of the optical phonons is negligible.
In some peculiar materials like graphene, it is customary to introduce also a fictitious branch called K-phonons constituted by the phonons having wave vectors close to the Dirac points, K or K ′ , in the first Brillouin zone (taking the origin in the center Γ of FBZ).Also in this case the Einstein approximation is used on account of the limited variability of the phonon energy near those points.Moreover, in graphene the phonons are classified as in-plane, representing vibration parallel to the material, and out of plane, representing vibrational mode orthogonal to the material.The LA, T A, LO, T O and K phonons are in plane.The out of plane phonons belong to the acoustic branch and are named ZA phonons.For them a quadratic dispersion relation is a good approximation ε ZA (q) = ω ZA (q) = α|q| 2 , where α = 6.2 × 10 7 m 2 /s (see [20]) In the following, instead of the wave vector q we will use the phonon moment q, but we retain the same character for sake of simplifying the notation.So that ε s (q) = c s |q|, s = LA, T A and ε ZA (q) = α|q| 2 , where α = α/ .
The thermal transport is usually described by macroscopic models, e.g. the Fourier one, those based on the Maximum Entropy methods [19] or on phenomenological description [10].A more accurate way to tackle the question is to resort to semiclassical transport equations, the so-called Peierls-Boltzmann equations, for each phonon branch for the phonon distributions f µ (t, x, q) where c µ = ∇ q ( ω µ ) is the group velocity of the µth phonon specie.The phonon collision term C µ splits into two terms C µ µ describes the phonon interaction within the same branch while C ν µ describes the phonon-phonon interaction between different species.To deal with the complete expressions of the C µ 's is a very complicated task even from a numerical point of view [21].So, they are usually simplified by the Bhatnagar-Gross-Krook (BGK) approximation which mimics the relaxation of each phonon branch towards a common local equilibrium condition, characterised by a local equilibrium temperature T L that is the same for each phonon population.The local equilibrium phonon distributions are given by the Bose-Einstein distributions and the functions τ µ are the phonon relaxation times.Additional BGK terms can be added to include the interaction between pairs of different branches.
If we know the phonon distributions f µ 's, we can calculate the average phonon energy densities and the expectation value of any function ψ(q) for example the energy flux if one takes ψ(q) = ω µ v µ .
The modern devices, e.g. the electron ones like double gate MOSFETs (see [19]), are undergoing more and more miniaturization.This implies that the characteristic scales are of the same order as the typical lengths where quantum effects become more and more relevant.Therefore, quantum effects must be included and the semiclassical phonon transport equations must be replaced by a more accurate model.Among the possible approaches, that one based on the Wigner equation has the advantage to be formulated in a phase-space, allowing us to guess the feature of the solutions in analogy with the semiclassical counterpart.
A huge literature has been devoted to the application of the Wigner equations to charge transport (see [2,3,4]) but a limited use has been made for phonon transport.In the next sections a transport model, based on the Wigner quasi distribution, will be devised for phonon transport in nano-structures.

Phonon Wigner functions
The main point of our derivation is the kinetic description of a one-particle quantum statistical state, given in terms of one-particle Wigner functions.Let us now briefly recall the basic definitions and properties.A mixed (statistical) one-particle quantum state for an ensemble of scalar particles in R d is described by a density operator ρ, i.e. a bounded non-negative operator with unit trace, acting on L 2 (R d , C).Given the density operator ρ on L 2 (R d , C), the associated Wigner function, w = w(x, q), (x, q) ∈ R 2d , is the inverse Weyl quantization of ρ, We recall that the Weyl quantization of a phase-space function (a symbol) a = a(x, q) is the (Hermitian) operator Op (a) formally defined by [22] Op (a)ψ( x + y 2 , q ψ(y)e i(x−y)•q/ dy dq (6) for any ψ ∈ L 2 (R d , C).The inverse quantization of ρ can be written as the Wigner transform of the kernel ρ(x, y) of the density operator.
The dynamics of the time-dependent phonon Wigner functions g µ (x, q, t), µ = LA, T A, . . .steams directly from the dynamics of the corresponding density operator ρµ (t), i.e. from the Von Neumann or quantum Liouville equation where Ĥµ denotes the Hamiltonian operators of the µth phonons and [•, •] the commutator.
If h µ = Op −1 ( Ĥµ ) is the symbol associated with Ĥµ , then, from Eq.s (8), we obtain the Wigner equation for each phonon species With the symbol # we have denoted the Moyal (or twisted) product which translates the product of operators at the level of symbols according to for any pair of symbols a and b.Here, we do not tackle the analytical issues which guarantee the existence of the previous relations but limit ourselves to the remark that if two operators are in the Hilbert-Schmidt class, that is the trace there exists and it is not negative and bounded, then the product is still Hilbert-Schmidt and the Moyal calculus is well defined.
In the sequel, we will suppose that such conditions are valid.The Moyal product, under suitable regularity assumptions (see [23]), possesses the following formal semiclassical expansion where x i and similarly for β.
The expansion (11) can be rewritten as where The first terms of (13) read where ∇ 2 denotes the Hessian matrix and : the contracted product of tensors.It is easy to see that a# n b(x, q) = (−1) n b# n a(x, q), that is the operation # n is commutative (respectively anticommutative) when n is even (respectively odd).
If we neglect, for the moment, the phonon-phonon interactions, the Hamiltonian symbol for each phonon branch is given by By using the Moyal calculus, one can expand the second members of the previous Wigner equations.Up to first order in 2 , we have where1 The previous equations describe only ballistic transport and include only the harmonic contribution to the Hamiltonian.In order to describe also intra and inter-branch phononphonon interactions, an additional anharmonic term Ĥint encompassing the high order correction to the Hamiltonian operator must be added.So doing, one gets the so-called Wigner-Boltzmann equation In the quantum case the expression of C µ is rather cumbersome.For electron transport in semiconductors the interested reader can see [24].In certain regimes it is justified to retain the same form of the semiclassical collision operator as the semiclassical case [4].Here, we adopt a quantum BGK approach and model the collision terms as where g LE µ are now Wigner functions of local equilibrium which will be defined later.The equation (20) along with the expression (21) for the collision operator represents our starting point for the phonon transport.Note that for the optical phonons under the Einstein approximation for the energy bands one has formally the same transport equation as the semiclassical case because the group velocity vanishes.
An alternative derivation of (20) can be obtained by explicitly writing the von Neumann equation (see [18,19] for the details).One obtains whose expansion is of course in agreement with the Moyal calculus.

Phonon Moment equations
Getting analytical solutions to equations ( 20)-( 21) is a daunting task.Therefore, viable approaches are numerical solutions based on finite differences or finite elements [2] or stochastic solutions, e.g.those obtained with a suitable modification of the Monte Carlo methods for the semiclassical Boltzmann equation [3].However, it is possible to have simpler, even if approximate, models resorting to the moment method for the expectation values of interest.In fact, it is well known that, although not positive definite, the Wigner function is real and the expectation values of an operator can be formally obtained as an average of the corresponding symbol with respect to g µ (x, q, t).So, for any regular enough weight function ψ(q), let us introduce the short notation which represents a partial average with respect to the phonon moment q.More in general, if a = a(x, q) is a smooth symbol then it is possible to prove that the expectation of the (hermitian) operator A = Op (a) satisfies2 where k A (x, y) is the kernel of A.
We want to consider a minimum set of moments whose physical meaning is well clear.In particular, we shall consider the phonon energy and energy flux of each branch Note that the latter is directly related to the heat flux.
The evolution equations for W µ (x, t) and Q µ (x, t) are obtained by multiplying the relative Wigner equation by h µ (q), and h µ (q)c µ and integrating with respect to q We implicitly assume that the resulting integrals there exist, at least in the principal value sense.In order to get some global insight from eq.s (25), we formally assume the following expansions for each phonon branch3 g µ (x, q, t) = g (0) µ (x, q, t) It is possible to prove, at least formally [6], that the semiclassical Boltzmann equation is recovered from the Wigner equation as → 0 + .Therefore, g µ (x, q, t) can be considered as the solution f µ of the semiclassical transport equation.Accordingly, we write where (2)  µ (x, q, t)dq, (2)  µ (x, q, t)dq.
Regarding the moments of the collision terms, only with drastic simplifications analytical expressions can be deduced.In analogy with the BGK approximation, if an average relaxation time independent on q is considered, one can expand the r.h.s. of eq.s (25) up to first order in 2 as a relaxation time terms where Note that in the evaluation of the production term of the equations for the energy-fluxes the isotropy of the equilibrium Wigner function has been invoked and therefore Q LE µ vanishes.
The energy and energy-flux relaxation times, τ W µ and τ Q µ respectively, are assumed to depends on the temperature, which will be definite in the next section, of the relative branch.
Altogether, the resulting model is made of the following fluid-type equations where µ with µ (x, q, t)dq, (2)  µ (x, q, t)dq, and the complete symmetric tensors T µ and U µ have components ∂ 3 h µ (q) ∂q i ∂q j ∂q k g (0) µ (x, q, t)dq.
If we split into zero and first order in 2 , the evolution equations read The zero order equations are the model already investigated in several papers [1,7] where is proved that it is a hyperbolic system of conservation law.So, finite propagation speed of disturbances in energy is guaranteed to overcome the well-known paradox of the classical Fourier law for the heat flux.However, the first order corrections in 2 introduce dispersive terms and it seems that in a quantum regime the requirement of finite propagation speed of the thermal effects cannot be fulfilled.On the other hand, this is not surprising since the nonlocal character of the quantum evolution equations can lead to energy propagation without a bounded speed.

QMEP for the closure relations
The evolution equations ( 29)-(32) do not form a closed system of balance laws.If we assume the energies W µ and the energy-fluxes Q µ as the main fields, in order to get a set of closed equations we need to express the additional fields J µ , T µ and U µ as functions of W µ and Q µ .A successful approach in a semiclassical setting is that based on the Maximum Entropy Principle (MEP) (see also [19] for a complete review) which is based on a pioneering paper of Jaynes [12,13] who also proposed a way to extend the approach to the quantum case.The MEP in a quantum setting has been the subject of several papers [8,14,15,16,17] with several applications, e.g. to charge transport in graphene [1,18].Here we will use such an approach for phonon transport.
The starting point is the entropy for the quantum system under consideration.In [18] the authors have employed the Von-Neumann entropy which, however, does not take into account the statistical aspects.Therefore, we take as entropy a generalization of the classical one for bosons.Let us introduce the operator which must be intended in the sense of the functional calculus.Here k B is the Boltzmann constant.The entropy of the µ-th phonon branch reads which can be viewed as a quantum Bose-Einstein entropy.
According to MEP, we estimate ρµ with ρMEP µ which is obtained by maximizing S(ρ µ ) under the constraints that some expectation values have to be preserved.In the semiclassical point case, one maximizes the entropy preserving the values of the moments we have taken as basic field variables where is the vector of the weight functions and g M EP µ is the Wigner function associated with ρMEP µ .In the previous relations the time t and position x must be considered as fixed.
The quantum formulation of MEP is given in terms of expectation values E 1 (t) = tr {ρ µ Op (h µ (q))} (t), E 2 (t) = tr {ρ µ Op (c µ h µ (q))} (t), as follows: for fixed t in the space of the Hilbert-Schmidt operators on L 2 (R d , C) which are positive, with trace one and such that the previous expectation values there exist.Note that we are applying the maximization of the entropy for each phonon branch separately.In other words, we are requiring the additivity of the entropy.If we introduce the vector of the Lagrange multipliers the vector of the moments and the vector of the moments which must be considered as known the constrained optimization problem (36)-(37) can be rephrased as a saddle-point problem for the Lagrangian in the space of the admissible ρµ and smooth function η µ .
If the Lagrangian L µ (ρ µ , η µ ) is Gâteaux-differentiable with respect to ρµ , the first order optimality conditions require δL µ (ρ µ , η µ )(δ ρ) = 0 for each Hilbert-Schmidt operators δ ρ on L 2 (R d , C) which is positive, with trace one and such that the previous expectation values there exist.The existence of the first order Gâteaux derivative is a consequence of the following Lemma (for the proof see [25]; an elementary proof in the case of discrete spectrum is given in [14]).Lemma 1.If r(x) is a continuously differentiable increasing function on R + then tr{r(ρ)} is Gâteaux-differentiable in the class of the Hermitian Hilbert-Schmidt positive operators on L 2 (R d , C).The Gâteaux derivative along δρ is given by δtr{r(ρ)}(δ ρ) = tr {r ′ (ρ)δ ρ} . (42) The extremality conditions for the unconstrained minimization problem (36)-( 37) are similar to that of the semiclassical case, as expressed by the following lemma (see [14]).
Lemma 2. The first order optimality condition for the minimization problem (36)-( 37) is equivalent to where (s ′ ) −1 is the inverse function of the first derivative of s.
Proof.By applying Lemma 1, the Gâteaux derivative of the Lagrangian is given by ∀δ ρ perturbation in the class of the Hermitian Hilbert-Schmidt positive operators on Since the function s(x) is concave, s ′ (x) is invertible.Explicitly we have and the operator solving the first order optimality condition reads Moreover, such an operator is a point of maximum for the Lagrangian.Now, to complete the program we have to determine, among the smooth functions, the Lagrange multipliers η µ by solving the constraint If such an equation has a solution η * µ , altogether the MEP density operator reads where we have rescaled the Lagrange multipliers including the factor 1/k B .
To determine conditions under which the equation ( 45) admits solutions is a very difficult task.Even in the semiclassical case there are examples (see [26]) of sets of moments that cannot be moments of a MEP distribution.We will directly find out the solution at least up to first order in 2 .
Once the MEP density function has been determined, the MEP Wigner function is given by which can be used to get the necessary closure relations by evaluating the additional fields with g µ replaced by g M EP µ .We remark that the constraints (45) can be more conveniently expressed as and indeed we will require, in analogy with the semiclassical case, the stronger conditions where the Lagrange multipliers enter through g M EP µ (x, q, t).

Determination of the Lagrange Multipliers
For the sake of making lighter the notation, let us consider a single branch and drop the index µ in the Wigner function in this section.We look formally for a solution in powers of g M EP = g M EP 0 firstly without taking into account the dependence of the Lagrange multipliers on .Of course, on account of the properties of the Weyl quantization, g M EP 0 is equal to the semiclassical counterpart [22] In order to determine the higher order terms g M EP k , k ≥ 1, given a symbol a(x, q) let us introduce the so-called quantum exponential Exp(a) defined as which can be expanded as Proposition Let a(x, p) be a smooth symbol.Then the following expansion is valid where Einstein's convention has been used.The proof can be found for example in [14].
By using what is proved in [16], we have where # 2l are the even terms of the Moyal product expansion and In particular e ξ (e ξ − 1) 3 (e ξ + 1) Therefore, up to first order in 2 we have and the constraints for each phonon branch read The previous equations form a nonlinear system of PDEs for the Lagrange multipliers whose analytical solution seems very difficult to get.Indeed, the situation is even more cumbersome because in a numerical scheme the inversion of the constraints should be performed at each time step.A viable strategy is to use the Lagrange multipliers as field variables by rewriting the evolution equations (28) in the form getting a highly nonlinear system of PDEs.Note that both ∇ η W and ∇ η Q i contain space derivatives of η.
A further simplification can be obtained by expanding the Lagrange multipliers as Therefore, the basic fields are also expanded with respect to 2 where The balance equations become We observe that the equations (55)-(56) decouple.Once they are solved, one can get the second order term of the Lagrange multipliers from (57)-(58) which form a linear system for η (2) .This is rather beneficial from a computational point of view Proposition 1.At zero order in 2 the map η → M(η) is (locally) invertible.
The proofs can be found in [19].

Local equilibrium temperature and heat conductivity
The concept of temperature out of equilibrium is a subtle topic and still a matter of debate.In the case of charge transport in semiconductors often the phonons are considered as a thermal bath and under some reasonable assumptions one can hypothesize that the electrons are in thermal equilibrium with the bath.In general if the dynamics of the phonons must be included, a thermal bath for these does not exist, unless a thermostated system is considered.Therefore, we need to introduce a local equilibrium temperature for the overall phonon system.
In statistical mechanics, one of the most reasonable and adopted ways to generalize the concept of temperature in a non-equilibrium state is that of relating it to the Lagrange multipliers associated with the energy constraint.For the phonon transport in graphene, an approach based on the Lagrange multipliers was followed in [1] (which the interested reader is referred to for the details).Let us recall here the main features.At equilibrium, the phonon temperatures and the corresponding Lagrange multipliers are related by If we assume that such relations hold, even out of equilibrium, the definition of the local temperature can be given in terms of the Lagrangian multipliers as follows.
Definition 1.The local temperature of a system of two or more branches of phonons is , where η LE 0 (x) is the common Lagrange multiplier that the occupation numbers of the branches, taken into account, would have if they were in the local thermodynamic equilibrium corresponding to their total energy density, that is, the following: where the sum runs over all the phonon branches.
At global equilibrium the temperature is constant T = T and the Wigner function reduced to the Bose -Einstein distribution with the same temperature for each phonon branch.Let us consider a small perturbation δT (x) of the temperature in the sense that δ(x)/ T ≪ 1.We can expand g M EP µ in powers of δ(x)/ T We observe that typically the relaxation energy relaxation time is much longer than the energy-flux relaxation times, that is τ Q ≪ τ W .In a long time scaling, much longer than τ Q , we get The relation between the Lagrange multipliers and the basic fields, as seen, can hardly be inverted analytically but a numerical procedure is necessary.However, if we consider a situation where the system is not too far from the equilibrium an expansion of the Lagrange multipliers around the equilibrium state can be performed.At equilibrium g M EP is isotropic and therefore η equil 1 = 0 and in a neighborhood of the equilibrium η 1 remains small.Therefore, for small deviations from the thermodynamic equilibrium the expansion is valid.By substituting in (61) one gets up to first order in 2 In particular, at the zero order we have which can be written in the Fourier form with the thermal conductivity tensor given by e hµ(q)/k B T (e hµ(q)/k B T − 1) 2 dq.
It is evident that K µ is positive definite.
Observe that ∀n ∈ S d S d being the unit sphere in R d .Therefore, if h µ (q) is isotropic K µ is isotropic as well with I identity matrix of order d and k (0) the zero order trace e hµ(q)/k B T (e hµ(q)/k B T − 1) 2 dq.
Indeed the last term in the previous relation is of order 2 δT T and can be considered negligible for small deviations from local equilibrium.The remaining part gives a highly nonlinear correction which cannot be put in a Fourier form.
As an example we consider the case of the longitudinal and transversal acoustic phonons in the Debye approximation for a single branch.In such a case the corresponding symbol of the phonon hamiltonian reads c|q| and therefore B |q|(e ξ − 1)T (x, t) 5 (δ ij |q| 2 − q i q j ) ∂T ∂x i ∂T ∂x j − q i q j T ∂ 2 T ∂x i ∂x j = − 1 8 c 2 e ξ (e ξ − 1) 3 (e ξ + 1) with now ξ = c|q|/k B T .Therefore, the second order correction to the heat flux is given by µ with