About the Definition of the Local Equilibrium Lattice Temperature in Suspended Monolayer Graphene

Abstract

The definition of temperature in non-equilibrium situations is among the most controversial questions in thermodynamics and statistical physics. In this paper, by considering two numerical experiments simulating charge and phonon transport in graphene, two different definitions of local lattice temperature are investigated: one based on the properties of the phonon–phonon collision operator, and the other based on energy Lagrange multipliers. The results indicate that the first one can be interpreted as a measure of how fast the system is trying to approach the local equilibrium, while the second one as the local equilibrium lattice temperature. We also provide the explicit expression of the macroscopic entropy density for the system of phonons, by which we theoretically explain the approach of the system toward equilibrium and characterize the nature of the equilibria, in the spatially homogeneous case.

1. Introduction

The definition of temperature in non-equilibrium situations is among the most controversial questions in thermodynamics and statistical physics (for a comprehensive review, the interested reader is referred to [1,2,3,4]). Several approaches have been adopted, trying to generalize what is already known in thermostatics, e.g., the non-equilibrium temperature is defined as the partial derivative of the entropy with respect to the internal energy, even if this translates the problem to the formulation of non-equilibrium entropy, which can be introduced on the basis of considerations of extended thermodynamics or non-equilibrium statistical mechanics. In [5], a contact temperature is introduced in which a non-equilibrium system is put in contact with a control equilibrium environment and the global heat flux through the boundaries is zero, the contact temperature being that of the environment. In [6], by extending Rugh’s ideas [7] to non-equilibrium situations, a configurational temperature is introduced, related to the forces between the particles of the system under consideration.

The attempt to validate an approach by comparison with experiments is intrinsically ambiguous since it is not well understood as to what it is indeed measured by the apparatus.

Here, we focus on the issue of the definition of the local equilibrium temperature in the case of the lattice temperature of graphene. Recently, this material has received much attention, in particular because of its thermal and electrical properties [8]. At the kinetic level, the thermal effects are described with a phonon gas, which obeys a suitable Boltzmann equation. Likewise, charge transport is accurately described by a semiclassical transport equation. Electron–phonon scatterings couple the two subsystems. In [9,10] electro-thermal effects were analyzed by a Monte Carlo simulation (for results based on a finite difference scheme, see [11]; for results based on a discontinuous Galerkin scheme, see [12]), adapting the method in [13,14,15] with a description of the phonon–phonon scatterings based on a Bhatnagar–Gross–Krook (BGK) approximation, where a proper definition of a local temperature of the graphene crystal lattice is needed. This latter is deduced from the conservation properties of the phonon–phonon collision terms.

A different approach for the analysis of the electro-thermal effects in graphene was adopted in [16,17], where a hydrodynamical model based on the moment method with closure relations obtained by exploiting the maximum entropy principle (MEP) [18,19,20,21], was formulated (for similar results with phonons considered to be a thermal bath, see [22,23]). Here, in agreement with rational extended thermodynamics [24], a definition of a different local temperature was introduced, which is directly related to the Lagrange multiplier associated to the phonon energy (a hybrid MC–hydrodynamical approach can be found in [25]). We remark that the problem is indeed very complex because, even if we restrict the analysis to the phonons, there are several degrees of freedom represented by the different phonon branches, each of them having its own temperature. The phonon–phonon collision operator creates a mechanism of energy exchange, which pushes the several degrees of freedom to a common temperature in the absence of an external field. Therefore, there are two important questions to face: to give, out of equilibrium, a suitable definition of temperature both for each phonon species and for the crystal lattice.

On account of the difficulties with the experimental results, we face the question of analyzing the meaning of the above-mentioned two local lattice temperatures with two numerical experiments. In the first one, coupled electron–phonon simulations are performed, while in the second one, only phonon transport is considered. The hydrodynamic model in [17] is adopted, with a more accurate approximation of the phonon–phonon collision operator. The obtained results, particularly the second set of simulations, strongly indicate that the appropriate definition of the local equilibrium lattice temperature should be that based on the energy Lagrange multipliers, also in accordance with the current non-equilibrium theory [24,26]. In fact, in the second numerical experiment, it remains constant, as it should be, when the total energy is conserved at variance with the other local temperature, which, instead, presents an overshoot during the transient before asymptotically tending to the other temperature. Nevertheless, the definition based on the conservation property of the phonon collision operator has also a physical meaning and can be regarded as a measure of how fast the system is trying to approach the equilibrium.

The plan of the paper is the following. In Section 2, the kinetic model and the definition of the first local temperature are presented. In Section 3, a macroscopic model and the definition of the second local temperature are introduced. In Section 4, the results of the numerical simulations are shown and the difference between the two local temperatures is critically highlighted. Eventually, in Section 5, we furnish the expression for the macroscopic entropy density for the phonon system by which we prove the asymptotic stability of the phonon equilibrium states and theoretically justify the behavior of the phonon temperatures in the second numerical experiment.

2. The Kinetic Model and the Definition of the First Local Temperature

Graphene consists of carbon atoms arranged in a honeycomb hexagonal lattice. The charge transport is essentially due to the electrons, which are located around the Dirac points, K and $K^{\prime }$, which are the vertices of the hexagonal primitive cell of the reciprocal lattice. At the Dirac points, the valence and conduction bands touch each other, which makes graphene a gapless semimetal. Moreover, having the energy bands in an approximately conical shape, electrons behave as massless Dirac fermions [8]. If high enough, Fermi levels are considered. It is possible to neglect the dynamics of the electrons in the valence bands, being that the latter ones are fully occupied in this case [27]. Such a situation is similar to n-type doping for traditional semiconductors. As said, around the equivalent Dirac points, the band energy ${\epsilon }_{\ell }$ is approximately linear.

$${\epsilon }_{\ell }=\hslash v_F\left|\mathbf{k}-{\mathbf{k}}_{\ell }\right|,$$

(1)

and the group velocity is given by the following:

$${\mathbf{v}}_{\ell }={\displaystyle \frac{1}{\hslash }}{\nabla }_{\mathbf{k}}{\epsilon }_{\ell }=v_F\frac{\mathbf{k}-{\mathbf{k}}_{\ell }}{|\mathbf{k}-{\mathbf{k}}_{\ell }|}.$$

where $\mathbf{k}$ is the electron wave-vector, $v_F$ is the (constant) Fermi velocity, ℏ is the reduced Planck constant, and ${\mathbf{k}}_{\ell }$ is the position of the Dirac point $\ell =K,K^{\prime }$.

In the framework of the semiclassical kinetic theory, relative to the electrons in the conduction band, the charge transport is described by two Boltzmann equations for the K and $K^{\prime }$ valleys:

$$\begin{array}{c}\hfill \frac{\partial f_{\ell }(t,\mathbf{x},\mathbf{k})}{\partial t}+{\mathbf{v}}_{\ell }·{\nabla }_{\mathbf{x}}{f}_{\ell }(t,\mathbf{x},\mathbf{k}){\left.-\frac{e}{\hslash }\mathbf{E}·{\nabla }_{\mathbf{k}}{f}_{\ell }(t,\mathbf{x},\mathbf{k})={\displaystyle \frac{d{f}_{\ell }}{dt}}(t,\mathbf{x},\mathbf{k})\right|}_{e-ph},\end{array}$$

(2)

where $f_{\ell }(t,\mathbf{x},\mathbf{k})$ represents the distribution function of the electrons in the valley ℓ (K or $K^{\prime }$), at position $\mathbf{x}$, time t, and with wave-vector $\mathbf{k}$. By ${\nabla }_{\mathbf{x}}$ and ${\nabla }_{\mathbf{k}}$, we denote the gradients with respect to the position and the wave-vector, respectively; e is the elementary (positive) charge; and $\mathbf{E}$ is the external applied electric field.

The collision term at the right-hand side of (2) describes the scatterings occurring between electrons and phonons. They can be with longitudinal, transversal, acoustic or optical phonons, which are labeled by LA, TA, LO and TO, respectively. Both the acoustic and the optical phonon scatterings are intra-valley and intra-band. One also has to take into account the electron scattering with K phonons, which is inter-valley, pushing electrons from a valley to the nearby one. The general form of the collision term is as follows:

$$\begin{array}{ccc}\hfill {\left.{\displaystyle \frac{d{f}_{\ell }}{dt}}(t,\mathbf{x},\mathbf{k})\right|}_{e-ph}=\sum _{{\ell }^{\prime }}\left[{\int }_{{\mathcal{I}}_{{\ell }^{\prime }}}{S}_{{\ell }^{\prime },\ell }({\mathbf{k}}^{\prime },\mathbf{k})f_{{\ell }^{\prime }}(t,\mathbf{x},{\mathbf{k}}^{\prime })\left(1-f_{\ell }(t,\mathbf{x},\mathbf{k})\right)d{\mathbf{k}}^{\prime }\right.& & \\ \hfill \left.-{\int }_{{\mathcal{I}}_{\ell }}{S}_{\ell ,{\ell }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime })f_{\ell }(t,\mathbf{x},\mathbf{k})\left(1-f_{{\ell }^{\prime }}(t,\mathbf{x},{\mathbf{k}}^{\prime })\right)d{\mathbf{k}}^{\prime }\right],& \end{array}$$

where the total transition rate $S_{{\ell }^{\prime },\ell }({\mathbf{k}}^{\prime },\mathbf{k})$ is given by the sum of the contributions of the several types of scatterings [17]:

$$\begin{array}{ccc}\hfill S_{{\ell }^{\prime },\ell }({\mathbf{k}}^{\prime },\mathbf{k})& =& \sum _{\mu }{\left|G_{{\ell }^{\prime },\ell }^{\left(\mu \right)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2\left[\left(g_{\mu }^{-}+1\right)\delta \left({\epsilon }_{\ell }\left(\mathbf{k}\right)-{\epsilon }_{{\ell }^{\prime }}\left({\mathbf{k}}^{\prime }\right)+\hslash {\omega }_{\mu }\right)\right.\hfill \\ & & \left.+g_{\mu }^{+}\delta \left({\epsilon }_{\ell }\left(\mathbf{k}\right)-{\epsilon }_{{\ell }^{\prime }}\left({\mathbf{k}}^{\prime }\right)-\hslash {\omega }_{\mu }\right)\right].\hfill \end{array}$$

(3)

The index $\mu$ labels the $\mu$th phonon mode and ${\omega }_{\mu }$ is its angular frequency; ${\mathcal{I}}_{\ell }$ labels the region in the Brillouin zone corresponding to the valley ℓ. The ${\left|G_{{\ell }^{\prime },\ell }^{\left(\mu \right)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2$’s are the electron–phonon coupling matrix elements, which describe the interaction mechanism of a $\mu$th phonon with an electron, going from the state of wave-vector ${\mathbf{k}}^{\prime }$ belonging to the valley ${\ell }^{\prime }$ to the state with wave-vector $\mathbf{k}$, belonging to the valley ℓ. The symbol $\delta$ denotes the Dirac distribution; $g_{\mu }(t,\mathbf{x},\mathbf{q})$ is the phonon distribution for the $\mu$-type phonons; and $\mathbf{q}$ is the phonon wave-vector belonging to the first Brillouin zone $\mathcal{B}$, measured from $\Gamma$ or K, respectively, for $\Gamma$ and K phonons. In Equation (3), $g_{\mu }^{\pm }=g_{\mu }\left({\mathbf{q}}^{\pm }\right)$, where ${\mathbf{q}}^{\pm }=\pm \left({\mathbf{k}}^{\prime }-\mathbf{k}\right)$, stemming from the momentum conservation. The K and $K^{\prime }$ valleys can be treated as equivalent; therefore, in the following, we will take into account a unique electron population.

Similarly, the evolution of the phonon populations is determined by the following Boltzmann equations for the phonon distributions as follows:

$$\begin{array}{ccc}& & \frac{\partial g_{op}}{\partial t}=C_{op},op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & \frac{\partial g_{ac}}{\partial t}+{\mathbf{c}}_{ac}·{\nabla }_{\mathbf{x}}{g}_{ac}=C_{ac},ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA},\hfill \end{array}$$

(5)

where ${\mathbf{c}}_{ac}={\nabla }_{\mathbf{q}}$${\omega }_{ac}$ is the acoustic phonon group velocity.

The optical phonon group velocity can be neglected because of the Einstein approximation ${ϵ}_{op}=\hslash {\omega }_{op}\approx const$, which can be used for their dispersion relation, while, regarding the LA and TA phonons, the Debye approximation can be employed: ${\omega }_{ac}={ϵ}_{ac}/\hslash =c_{ac}\left|\mathbf{q}\right|$, with $c_{ac}$ the sound speed of the branch ac = LA, TA. Eventually, the dispersion relation of the ZA phonons is approximately quadratic: ${ϵ}_{\mathrm{ZA}}=\hslash {\alpha }_{\mathrm{ZA}}{\left|\mathbf{q}\right|}^2$, with ${\alpha }_{ZA}$= 0.62 nm${}^2$/ps.

The phonon collision term splits into two parts as follows:

$$\begin{array}{ccc}& & C_{\mu }=C_{\mu }^{p-e}+C_{\mu }^{p-p},\mu =\mathrm{LA},\mathrm{TA},\mathrm{ZA},\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K}.\hfill \end{array}$$

(6)

where the term $C_{\mu }^{p-e}$ is the phonon–electron collision operator, while $C_{\mu }^{p-p}$ represents the phonon–phonon interactions, which are very difficult to treat from a numerical point of view. For this reason, a BGK approximation is commonly used for them [16]:

$$C_{\mu }^{p-p}=-\frac{g_{\mu }-g_{\mu }^{LE}}{{\tau }_{\mu }},$$

which describes the relaxation of each phonon branch toward the equilibrium condition that corresponds to the local equilibrium Bose–Einstein distributions.

$$\begin{array}{ccc}& & g_{\mu }^{LE}={\left[{\displaystyle e^{\hslash {\omega }_{\mu }/k_B{T}_L}-1}\right]}^{-1}.\hfill \end{array}$$

(7)

The functions ${\tau }_{\mu }={\tau }_{\mu }\left(T_{\mu }\right)$, which are reported in Figure 1, are the relaxation times, each of them depending only on the temperature $T_{\mu }$ of the branch under consideration. The temperature $T_L$ is the same for each phonon population and is defined as follows.

If the phonon distributions $g_{\mu }$ are known, the average phonon energy densities can be calculated as follows:

$$\begin{array}{ccc}& & W_{\mu }=\frac{1}{{\left(2\pi \right)}^2}{\int }_{\mathcal{B}}\hslash {\omega }_{\mu }{g}_{\mu }d\mathbf{q},\hfill \end{array}$$

(8)

where $\frac{1}{{\left(2\pi \right)}^2}$ represents the $\mu$th phonon density of the states, and the temperature $T_{\mu }$ of each phonon branch is determined from the following condition:

$${\int }_{\mathcal{B}}\hslash {\omega }_{\mu }{g}_{\mu }\left(\mathbf{q}\right)d\mathbf{q}={\int }_{\mathcal{B}}\hslash {\omega }_{\mu }{\left[{\displaystyle e^{\hslash {\omega }_{\mu }/k_B{T}_{\mu }}-1}\right]}^{-1}d\mathbf{q},$$

(9)

where $\mu$ varies over the phonon branches.

From the conservation of the total phonon energy, under processes involving only phonons, the following relation has to hold:

$$\begin{array}{c}\hfill \sum _{\mu }\frac{W_{\mu }-W_{\mu }^{LE}}{{\tau }_{\mu }}=0\end{array}$$

(10)

where $W_{\mu }^{LE}$ is calculated by means of (7) and (8). The above-written condition allows one to state the following definition of local temperature.

Definition 1.

We define $T_L$ as the solution of the non-linear relation (10), which can be obtained numerically; see e.g., [28,29] where this local temperature is called scattered phonon pseudo-temperature.

It is possible to prove that (10) admits a unique solution. For further details, we refer to [10], where the previous approach was adopted to devise a simulation scheme for the electron–phonon transport in graphene.

3. A Macroscopic Model and the Definition of the Second Local Temperature

Macroscopic models can be derived from the kinetic one [17,19,30] by taking suitable moments of the distribution functions as state variables. Here, we present in some detail only the evolution equations of the phonon variables, and refer the interested reader to [17] for a complete treatment of the problem.

If one chooses a certain number of moments of the electron and phonon distributions as state variables, the extra fluxes and the production terms, which are present in the corresponding balance equations, are additional unknown quantities, which require constitutive relations in terms of the state variables. By exploiting the maximum entropy principle (MEP), the electron and phonon distributions can be estimated by the so-called maximum entropy distributions $f_{MEP}$ and $g_{\mu ,MEP}$, $\mu$ = LO, TO, ZO, K, LA, TA, ZA, which solve the following maximization problem:

$$\underset{f,g_{\mu }}{max}S[f,g_{\mu }],$$

under the constraint that the moments chosen as fundamental variables are known. $S[f,g_{\mu }]$ is the total entropy of the physical system, which depends on the electron and phonon distribution functions f and $g_{\mu }$, and whose expression is reported in [17].

In particular, for the phonons, the following moments can be chosen:

$$\begin{array}{ccc}& & W_{op}=\frac{1}{{\left(2\pi \right)}^2}{\int }_{\mathcal{B}}\hslash {\omega }_{op}{g}_{op}d\mathbf{q},{\mathbf{P}}_{op}=\frac{1}{{\left(2\pi \right)}^2}{\int }_{\mathcal{B}}\hslash \mathbf{q}{g}_{op}d\mathbf{q},op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & W_{ac}=\frac{1}{{\left(2\pi \right)}^2}{\int }_{{\mathbb{R}}^2}\hslash {\omega }_{ac}{g}_{ac}d\mathbf{q},{\mathbf{Q}}_{ac}=\frac{1}{{\left(2\pi \right)}^2}{\int }_{{\mathbb{R}}^2}\hslash {\omega }_{ac}{\mathbf{c}}_{ac}{g}_{ac}d\mathbf{q},ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA},\hfill \end{array}$$

(12)

which respectively represent the energy and momentum densities of the optical phonons and the energy and energy flux densities of the acoustical phonons. Solving the above constrained maximization problem, one obtains the following:

$$\begin{array}{ccc}\hfill g_{op,MEP}& =& \frac{1}{exp\left({\lambda }_{W_{op}}{ϵ}_{op}+\hslash \mathbf{q}·{\lambda }_{{\mathbf{P}}_{op}}\right)-1},op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \\ \hfill g_{ac,MEP}& =& \frac{1}{exp\left({\lambda }_{W_{ac}}{ϵ}_{ac}+{ϵ}_{ac}{\mathbf{c}}_{ac}·{\lambda }_{{\mathbf{Q}}_{ac}}\right)-1},ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA},\hfill \end{array}$$

where the $\lambda$s are the Lagrange multipliers arising from the presence of the constraints.

In order to tackle the problem of the inversion of the constraints, the distribution functions are linearized around their isotropic part, obtaining the following:

$$\begin{array}{ccc}\hfill g_{op,MEP}& \approx & \frac{1}{e^{{\lambda }_{W_{op}}{ϵ}_{op}}-1}\left[1-\frac{e^{{\lambda }_{W_{op}}{ϵ}_{op}}}{e^{{\lambda }_{W_{op}}{ϵ}_{op}}-1}\hslash \mathbf{q}·{\lambda }_{{\mathbf{P}}_{op}}\right],op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \\ \hfill g_{ac,MEP}& \approx & \frac{1}{e^{{\lambda }_{W_{ac}}{ϵ}_{ac}}-1}\left[1-\frac{e^{{\lambda }_{W_{ac}}{ϵ}_{ac}}}{e^{{\lambda }_{W_{ac}}{ϵ}_{ac}}-1}{ϵ}_{ac}{\mathbf{c}}_{ac}·{\lambda }_{{\mathbf{Q}}_{ac}}\right],ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA}.\hfill \end{array}$$

By substituting the latter expressions into the constraints (11) and (12) and by solving them with respect to the Lagrange multipliers, one finds the following:

$$\begin{array}{ccc}& & {\lambda }_{W_{op}}=\frac{1}{{ϵ}_{op}}ln\left(1+\frac{yA{ϵ}_{op}}{W_{o}p}\right),op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & {\lambda }_{W_{ac}}={\left(\frac{4\pi y\zeta \left(3\right)}{{\hslash }^2{c}_{ac}^2}\right)}^{\frac{1}{3}}{W}_{ac}^{-\frac{1}{3}},ac=\mathrm{LA},\mathrm{TA},{\lambda }_{W_{\mathrm{ZA}}}=\sqrt{\frac{\pi }{24\hslash {\alpha }_{\mathrm{ZA}}{W}_{\mathrm{ZA}}}},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\lambda }_{{\mathbf{P}}_{op}}=-\frac{A^2{ϵ}_{op}^2}{4{\hslash }^2{A}_1}\frac{y}{W_{op}(W_{op}+Ay{ϵ}_{op})}{\mathbf{P}}_{op},op=\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & {\lambda }_{{\mathbf{Q}}_{ac}}=-\frac{2}{3}{\left(\frac{4\pi y\zeta \left(3\right)}{{\hslash }^2{c}_{ac}^8}{W}_{ac}^{-4}\right)}^{\frac{1}{3}}{\mathbf{Q}}_{ac},ac=\mathrm{LA},\mathrm{TA},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & {\lambda }_{{\mathbf{Q}}_{\mathrm{ZA}}}=-\frac{16}{15}\frac{\sqrt{\pi }}{{\alpha }_{\mathrm{ZA}}^3\zeta (5/2)}{\left(\frac{24\hslash {\alpha }_{\mathrm{ZA}}}{\pi W_{\mathrm{ZA}}}\right)}^{\frac{7}{4}}{\mathbf{Q}}_{\mathrm{ZA}},\hfill \end{array}$$

(17)

where $y=\frac{1}{{\left(2\pi \right)}^2}$, $\zeta (·)$ is the zeta function, $A=\frac{8\sqrt{3}}{9}\frac{{\pi }^2}{a_0^2},$$A_1=\frac{20\sqrt{3}}{729}\frac{{\pi }^4}{a_0^4},$ with $a_0=0.142$ nm being the nearest neighbor distance between the atoms in graphene.

Eventually, using the MEP distribution functions and the relations expressing the Lagrange multipliers as functions of the fundamental variables, it is possible to close the moment equations as follows:

$$\begin{array}{ccc}\hfill \frac{\partial W_{op}}{\partial t}& =& C_{W_{op}},\frac{\partial {\mathbf{P}}_{op}}{\partial t}=C_{{\mathbf{P}}_{op}},op=\mathrm{LO},\mathrm{TO},\mathrm{K},\mathrm{ZO},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \frac{\partial W_{ac}}{\partial t}+{\nabla }_{\mathbf{x}}·{\mathbf{Q}}_{ac}& =& C_{W_{ac}},\frac{\partial {\mathbf{Q}}_{ac}}{\partial t}+{\nabla }_{\mathbf{x}}·{\mathbf{T}}_{ac}=C_{{\mathbf{Q}}_{ac}},ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA},\hfill \end{array}$$

(19)

by means of the following closure relations for the production terms $C_{W_{op}},C_{{\mathbf{P}}_{op}},C_{W_{ac}}$, $C_{{\mathbf{Q}}_{ac}}$, and the acoustical phonon fluxes of energy fluxes ${\mathbf{T}}_{ac}:=y{\int }_{\mathcal{B}}{ϵ}_{ac}{\mathbf{c}}_{ac}\otimes {\mathbf{c}}_{ac}{g}_{ac}\mathbf{d}q$:

$$\begin{array}{ccc}& & C_{W_{\mu }}=-\frac{W_{\mu }-W_{\mu }^{LE}}{{\tau }_{\mu }},C_{{\mathbf{P}}_{op}}=-{\mathbf{P}}_{op}/{\tau }_{op},\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & C_{{\mathbf{Q}}_{ac}}=-{\mathbf{Q}}_{ac}/{\tau }_{ac},\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& {\mathbf{T}}_{LA/TA}& =\frac{c_{LA/TA}}{2}{W}_{LA/ZA}\mathbf{I},{\mathbf{T}}_{ZA}=\frac{\zeta \left(3\right)}{\pi {\hslash }^2}{\left(\frac{\pi }{24\hslash {\alpha }_{ZA}}{W}_{ZA}\right)}^{\frac{3}{2}}\mathbf{I},\hfill \end{array}$$

(22)

where $\mathbf{I}$ is the identity matrix, the indices vary over the above-specified phonon branches, and only the productions due to the phonon interactions among themselves are considered.

In the previous section, we introduced a local lattice temperature by the relation (10), which stems from the properties of the phonon–phonon collision operator. However, the concept of temperature out of equilibrium is a subtle topic and still a matter of debate [24,26]. The rationale of the previous definition is that the collision operator pushes the system, in a characteristic time-related manner to the relaxation times toward a local equilibrium state with a single temperature for all the phonons. However, 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 to the energy constraint.

For the phonon transport in graphene, the approach based on the Lagrange multipliers was followed in [17] (which the interested reader is referred to for the details) within the application of the MEP (see [18,21] for a review of MEP in semiconductors). Let us recall here the main steps.

At equilibrium, the phonon temperatures are related to the corresponding Lagrange multipliers by means of the following:

$$\begin{array}{c}\hfill T_{\mu }=\frac{1}{k_B{\lambda }_{W_{\mu }}},\mu =\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\mathrm{LA},\mathrm{TA},\mathrm{ZA}.\end{array}$$

If we assume that such relations hold, even out of equilibrium, the definition of a second local temperature can be given in terms of the Lagrangian multipliers as follows.

Definition 2.

The local temperature of a system of two or more branches of phonons is $T_L^{*}:=\frac{1}{k_B{\lambda }_{W_L}}$, where ${\lambda }_{W_L}$ 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:

$$W\left({\lambda }_{W_L}\right):=\sum _{\mu }{W}_{\mu }\left({\lambda }_{W_L}\right)=\sum _{\mu }{W}_{\mu }\left({\lambda }_{W_{\mu }}\right),$$

where the sum is extended to the considered branches and the functions $W_{\mu }\left({\lambda }_{W_{\mu }}\right)$ are found from expressions (13) and (14).

In other words, we require that $T_L^{*}$ is such that, by evaluating all the average phonon energy densities with the Lagrange multiplier given by $1/k_B{T}_L^{*}$ and by summing them up, one obtains the value of the total average energy density.

The two definitions of local temperature are equivalent if, and only if, all the relaxation times are equal, that is, the following:

$$\begin{array}{c}\hfill {\tau }_{\mu }=\tau ,\mu =\mathrm{LO},\mathrm{TO},\mathrm{ZO},\mathrm{K},\mathrm{LA},\mathrm{TA},\mathrm{ZA},\end{array}$$

(23)

but this assumption is not compatible with the experimental data as clearly indicated in Figure 1, where the relaxation times of all the phonon branches are reported [31]. As a consequence, the two definitions of temperature do not coincide unless all the phonons are in local equilibrium among them, that is, they all have the same local temperature. We notice that $T_L^{*}$ is related only to the energy of the system and does not take into account any scattering mechanism. However, the collision terms are now expressed in terms of the Lagrange multipliers associated to the energies. In fact, the strength of the energy production terms is proportional, according to the relaxation times, to the differences between the energy densities relative to the different Lagrange multipliers and $W\left({\lambda }_{W_L}\right)$.

$$\begin{array}{c}\hfill \frac{W_{\mu }\left({\lambda }_{W_{\mu }}\right)-W\left({\lambda }_{W_L}\right)}{{\tau }_{\mu }\left(T_{\mu }\right)}.\end{array}$$

4. Two Numerical Experiments Investigating the Previous Definitions of the Local Temperature

In this section, we present two numerical experiments to get some insight into the behavior of the phonon system and the phonon temperatures. The first one refers to a suspended graphene sheet subjected to an external electric field with the aim of investigating how the temperatures of the various phonon branches increase due to the interaction with the electrons and among them. In the second one, the electric field is turned off, and we analyze how the phonons go to the equilibrium. Since we are interested in the behavior of the phonon system, we neglect the phonon interaction with the electrons for the sake of simplicity and also because it is weak as confirmed by the significant difference among the phonon temperatures and the electron one when a current flows through the graphene [32]. The influence of the electrons on the final temperature, on the basis of their much smaller surface heat capacity, can be estimated to be around 10%.

4.1. First Numerical Experiment

The graphene layer is considered infinite in one dimension, while the other one, which we call longitudinal, is taken as 0.1 $\mu$m long. A constant bias is applied to the transversal boundaries. For the simulation, we exploit the macroscopic model for electrons and phonons presented in [17] but using the phonon–phonon relaxation times adopted in [9,10], with ${\tau }_{\mathrm{K}}={(1/{\tau }_{\mathrm{LO}}+1/{\tau }_{\mathrm{TO}})}^{-1}$ and the temperature $T_L$ as the temperature of the local equilibrium distribution functions appearing in the phonon–phonon collision operators; see formula (7). There being no spatial dependence—in this case, $T_L$ and $T_L^{*}$ only depend on time—one has to solve a system of ordinary differential equations. As said, we focus our attention on the temperatures of the phonons of the various branches and in particular on the temperatures $T_L$ and $T_L^{*}$.

Unfortunately, a comparison with experiments is not easy because it is not clear what exactly is measured by the instruments.

Two different values of the electron Fermi energy are considered: ${\epsilon }_F=$ 0.4 and 0.6 eV, and for each of them two different electric fields are applied: $E=2$ and 5 kV/cm. A time window of 50 picoseconds is simulated starting from an equilibrium state at room temperature (300 K). In all the cases, see the left-hand sides of Figure 2, Figure 3, Figure 4 and Figure 5; it can be seen that the phonons that have the highest temperature are the LO ones since they have a greater energy exchange with the electrons. The temperature of the ZA phonons is significantly the lowest one, which means that they exchange less energy with the other phonons, while the temperatures of the other branches are near one another.

We also notice that initially, the LO, TO, and K phonon temperatures, in particular that of the LO phonons, increase more rapidly since the phonon–phonon decay acts on a longer time scale. As regards the temperatures $T_L$ and $T_L^{*}$, see the right-hand sides of Figure 2, Figure 3, Figure 4 and Figure 5; they begin to be substantially different on a time scale comparable with that of the phonon–phonon decay and after some time, their difference remains nearly constant.

4.2. Second Numerical Experiment

Guided by the results of the first experiment, we consider a sheet of graphene in which the initial temperatures of the phonons are the final ones of the first experiment. The electric field is turned off and only the phonon–phonon interactions are taken into account; therefore, the physical system can be described by the system of ordinary differential equations, which is obtained from (18) and (19) when no spatial dependence is present. Two cases are considered. In the first one, the initial phonon temperatures are as follows: $T_{\mathrm{LO}}=$ 408.6 K, $T_{TO}=$ 373.0 K, $T_K=$ 361.2 K, $T_{ZO}=$ 357.2 K, $T_{LA}=$ 359.7.2 K, $T_{TA}=$ 360.0 K and $T_{ZA}=$ 321.0 K. In the second one, they are as follows: $T_{\mathrm{LO}}=$ 528.4 K, $T_{\mathrm{TO}}=$ 468 K, $T_{\mathrm{K}}=$ 450.6 K, $T_{\mathrm{ZO}}=$ 444.8 K, $T_{\mathrm{LA}}=$ 448.7 K, $T_{\mathrm{TA}}=$ 449.2 K, and $T_{\mathrm{ZA}}=$ 362.4.8 K. The first thing that can be seen in the left-hand sides of Figure 6 and Figure 7 is that the ZA phonons (and, though less, also the ZO ones) have a great influence on the final temperature that is reached by the system; moreover, they have a greater inertia. The equilibrium is reached in about 60 ps and the final temperature is $T_L^{*}$, which is constant during the process (see the right-hand sides of Figure 6 and Figure 7), in agreement with the energy conservation and with the fact that the right-hand sides of the evolution equations are null if, and only if, all the phonon temperatures are equal to $T_L^{*}$.

As regards $T_L$, we can see in the right-hand sides of Figure 6 and Figure 7 that after an initial overshoot, which occurs in correspondence with the most rapid variation in the temperatures of the phonons, it decays toward $T_L^{*}$. Therefore, the numerical results suggest that $T_L^{*}$ has to be intended as the local equilibrium lattice temperature, while the difference between $T_L$ and $T_L^{*}$ can be interpreted as a measure of how fast the system is trying to approach the equilibrium.

5. The Entropy Density of the Phonon System in the Homogeneous Case

The phonon behavior observed in Section 4.2 can be theoretically proved starting from the entropy density of the phonon system, which is given by the following:

$$\begin{array}{c}\hfill S=-\frac{k_B}{{\left(2\pi \right)}^2}\sum _{\mu }{\int }_{\mathcal{B}}\left[g_{\eta }ln{g}_{\mu }-\left(1+g_{\mu }\right)ln\left(1+g_{\mu }\right)\right]d\mathbf{q}.\end{array}$$

Substituting in this formula the maximum entropy distributions and neglecting, in agreement with the small anisotropy hypothesis, the quadratic terms in the vector Lagrange multipliers, one obtains the following:

$$\begin{array}{ccc}\hfill S& \approx & k_B\sum _{\mu }\int {\lambda }_{W_{\mu }}d{W}_{\mu }=k_B\left\{\sqrt{\frac{\pi }{6\hslash {\alpha }_{ZA}}}{W}_{ZA}^{\frac{1}{2}}+\sum _{ac\ne ZA}{\left(\frac{\zeta \left(3\right)}{\pi {\hslash }^2{c}_{ac}^2}\right)}^{\frac{1}{3}}\frac{3}{2}{W}_{ac}^{\frac{2}{3}}\right.\hfill \\ & +& \left.\sum _{op}\left[\frac{W_{op}}{{ϵ}_{op}}ln\left(1+\frac{A{ϵ}_{op}}{{\left(2\pi \right)}^2{W}_{op}}\right)+\frac{A}{{\left(2\pi \right)}^2}ln\left(1+\frac{{\left(2\pi \right)}^2{W}_{op}}{A{ϵ}_{op}}\right)\right]\right\}.\hfill \end{array}$$

(24)

Remark 1.

From the previous expression for each type of phonon, the local temperature can be written as follows:

$$\frac{1}{k_B{T}_{\mu }}={\lambda }_{W_{\mu }}=\frac{\partial S}{\partial W_{\mu }}.$$

(25)

This strongly reminds of the entropic definition of temperature:

$$\frac{\partial S}{\partial u}=\frac{1}{k_{B}T}.$$

where u is the internal energy and T is the (entropic) temperature of the system. Note that such an approach cannot be used for the whole system because of the lack of a one-to-one correspondence between the total energy and the energies of the single species. Instead, for the whole system, definition 3 can be viewed as a sort of caloric temperature (see [33] for a more complete discussion on the subject).

The following theorem holds.

Theorem 1.

For any $T\in {\mathbb{R}}^{+}$, the entropy (24) with reverse sign, restricted to the hypersurface $\{(T_{\mathrm{LO}},T_{\mathrm{TO}},T_{\mathrm{K}},T_{\mathrm{ZO}},T_{\mathrm{LA}},T_{\mathrm{TA}},T_{\mathrm{ZA}})\in {\left(R^{+}\right)}^7|{\sum }_{\mu }{W}_{\mu }\left(T_{\mu }\right)={\sum }_{\mu }{W}_{\mu }\left(T\right)\}$, where μ runs over all the phonon types, is a Lyapunov function [34] for the dynamical system as follows:

$$\begin{array}{ccc}\hfill \frac{\partial W_{op}}{\partial t}& =& C_{W_{op}},op=\mathrm{LO},\mathrm{TO},\mathrm{K},\mathrm{ZO},\hfill \\ \hfill \frac{\partial W_{ac}}{\partial t}& =& C_{W_{ac}},ac=\mathrm{LA},\mathrm{TA},\mathrm{ZA}.\hfill \end{array}$$

Proof. 

One has the following:

$$\begin{array}{ccc}\hfill -\frac{\partial S}{\partial t}& =& k_B\sum _{\mu }{\lambda }_{W_{\mu }}\left(T_{\mu }\right)\frac{W_{\mu }\left(T_{\mu }\right)-W_{\mu }^{LE}\left(T_L\right)}{{\tau }_{\mu }\left(T_{\mu }\right)}\hfill \\ & =& k_B\sum _{\mu }\left({\lambda }_{W_{\mu }}\left(T_{\mu }\right)-{\lambda }_{W_{\mu }}\left(T_L\right)\right)\frac{W_{\mu }\left(T_{\mu }\right)-W_{\mu }\left(T_L\right)}{{\tau }_{\mu }\left(T_{\mu }\right)}\le 0,\hfill \end{array}$$

where the second equality is due to the energy conservation property of the collision operator (we remind that ${\lambda }_{W_{\mu }}\left(T_L\right)=\frac{1}{k_B{T}_L},\forall \mu$) and the inequality is due to the fact that the $W_{\mu }$ are increasing functions of their argument, while the ${\lambda }_{W_{\mu }}$ are decreasing functions. Furthermore, $-\frac{\partial S}{\partial t}=0$ if, and only if, $T_{\mathrm{LO}}=T_{\mathrm{TO}}=T_{\mathrm{K}}=T_{\mathrm{ZO}}=T_{\mathrm{LA}}=T_{\mathrm{TA}}=T_{\mathrm{ZA}}$. This implies that any $(T,T,T,T,T,T,T)\in {\left(R^{+}\right)}^7$ is a stable equilibrium state having, as its basin of attraction, the hypersurface $\{(T_{\mathrm{LO}},T_{\mathrm{TO}},T_{\mathrm{K}},T_{\mathrm{ZO}},T_{\mathrm{LA}},T_{\mathrm{TA}},T_{\mathrm{ZA}})\in {\left(R^{+}\right)}^7|{\sum }_{\mu }{W}_{\mu }\left(T_{\mu }\right)={\sum }_{\mu }{W}_{\mu }\left(T\right)\}.$ □

Eventually, in Figure 8, we report the time behavior of the entropy density for the two cases considered in the second numerical experiment, which show that the entropy is a convex function of time.

6. Conclusions

Two different definitions of local temperature were compared: one, $T_L$, based on the properties of the phonon–phonon collision operator, and the other, $T_L^{*}$, based on the energy Lagrange multipliers. Since it is difficult to have a clear meaning of the experimental results, the analysis was based on two numerical experiments represented by simulations of charge and phonon transport in suspended monolayer graphene. The results indicate that the temperature $T_L$ has to be intended as a measure of how fast the phonon system is trying to converge to the local equilibrium, while $T_L^{*}$ is the local equilibrium lattice temperature. We have also provided an explicit expression of the phonon entropy density in the homogeneous case.

The case of phonon transport in graphene is particularly challenging because, apart from its own interest due to the increasing importance of thermal effects in nano devices, it has several degrees of freedom—the several branches—and each of them possesses a different temperature. Therefore, one is faced with several important questions of non-equilibrium thermodynamics: the definition of temperature, the interaction among the several degrees of freedom, and the identification of a local temperature.

What was found contributes to the complex and controversial debate about the concept of temperature out of equilibrium. The definition based on the collision operator (Definition 1 of the paper) is employed within a kinetic context but it is little explored in a general framework of non-equilibrium thermodynamics; see the quoted review articles [2,3,4] where it is not explicitly mentioned. However, being that this temperature is related to the phonon–phonon interactions, that is, to the anharmonic interaction terms, we can say that it, in some sense, recalls the configurational temperature [6,7].

The definition based on the energy Lagrange multipliers is consistent with the statistical approach and the application of the maximum entropy principle [20]. Moreover, as shown in the present paper, for each species, it is related to the entropic definition of non-equilibrium temperature, while for the whole system, such an approach leads to a sort of caloric definition of local temperature out of equilibrium.