Eigenmode analysis of the multiple temperature model: spectrum properties, hierarchical structures, and temperature inversion

Recent developments of ultrafast laser pulse techniques enable us to study the subpicosecond scale dynamics out of thermal equilibrium. Multiple temperature models (MTMs) are frequently used to describe such dynamics where the total system is divided into subsystems each of which is in local thermal equilibrium. Typical examples include the electron-lattice two temperature model and electron-spin-phonon three temperature model. We present the exact analytical solutions of linear MTM, based on the Fourier series expansion and the Green’s function method. We then discuss their properties for the case of the two and three temperature models. We show that the general solution of MTM is expressed as a linear combinations of a spatially uniform, single-temperature stationary mode and the other non-oscillatory, decaying “eigenmodes” characterized by different wave vectors and well-defined mode lifetimes. The eigenmode picture enables us to explore the hierarchical structure of models with respect to space, time and the coupling parameter. Excitation by source term is included by the Green’s function method. As an example, we derive an analytical solution for a Gaussian type source term. We report a phenomenon “temperature inversion” where the lattice temperature exceeds electron’s temperature for ns time scale. Furthermore, we show how physical requirements such as energy conservation and equilibration are realized in the general linear MTM in terms of the eigenmode picture.


Introduction
The nonequilibrium description of the condensed matter systems has remained a subject of strong interest for decades in physics. Leaving from well-defined thermodynamical equilibrium states, a possible first step towards the description of nonequilibrium dynamics is to divide the total system into subsystems, each in local thermal equilibrium. The idea of separating the total system into electronic and lattice subsystems with different temperature dates back to 1950s [1]. Early developments of this idea is detailed in a review by Kabanov [2]. A present form of the two temperature model (2TM) can be found in early 1970s [3]. A theoretical proposal to measure the electron-phonon coupling strength by pump-probe experiments [4] followed by observations in superconducting metallic systems [5,6] has paved the way to a crucial application of the 2TM.
The 2TM is now applied to extreme conditions where melting, evaporation, and material removal occur by ultrafast laser excitation [7,8], e.g., during the ultrafast laser material processing, for which higher energy efficiency and spatial precision are expected. Subpicosecond laser pulse deposits energy on the electronic subsystem in a ultrashort time scale while the lattice temperature remains relatively low. The fast thermalization process of the electronic system is considered to justify that the electronic and lattice system possess different temperatures T e and T l after the laser pulse is turned off.
The limitation of the 2TM has been recognized early on. Its failures of predicting the electron-phonon relaxation time at low temperature and its excitation intensity dependence were pointed out [9,10]. Baranov and Kabanov derived a temperature range � 2 2 D ∕E F < k B T < � D (E F ∕� D ) 1∕3 where 2TM cannot be justified [11]. D is the Debye frequency, and E F is the Fermi energy. The Boltzmann equation approach is frequently used to improve the description of nonthermal distribution function [12][13][14][15][16]. For the description of material destruction processes, a multi-scale modeling which combines the 2TM and the classical molecular dynamics is employed [17][18][19]. A recent review can be referred for this approach [20].
Yet simple but a straightforward extension of the 2TM is dividing the system into smaller subsystems. Waldecker introduced an idea to generalize the 2TM to the nonthermal lattice model where three phonon modes of Al have their own temperatures [21]. A similar approach is applied to graphene [22]. The electron-spin-phonon three temperature model has been developed to explain the ultrafast demagnetization process [23][24][25], sometimes in combination with a microscopic equation of motion [24,25]. MTMs with a finite speed limit of the heat conduction are also numerically investigated [26,27].
In this paper, we present exact solutions of the linear multiple temperature model (MTM) whose coefficients are all constant. Under the condition of vanishing heat flow of each subsystem at the boundaries, the model can be diagonalized. The system dynamics can then be described by a linear combination of damping eigenmodes. Each eigenmode is characterized by the mode lifetime which depends on the wave vector q. We first discuss the exact solution of the 2TM, which splits into two eigenmodes ± (q) on q space. We will see that the + (q) mode is smoothly connected to the solution of an effective one temperature model (1TM) in a small q limit. While in an opposite, large q limit, ± (q) modes converge to the free diffusion modes of electron and lattice temperatures where the electron-lattice coupling G becomes negligible. The mode profile indicates a counterintuitive phenomenon where the lattice temperature exceeds the electron temperature, which we call "temperature inversion". We nextly discuss the 3TM describing an electron-phonon-spin system of nickel thin film. While some features are common to the 2TM, a striking feature is a mode almost dispersionless over a wide range of q space. This is a counterintuitive behaviour since the three subsystems are mutually strongly coupled. We also investigate a "weak coupling limit" G 12 , G 13 ≪ G 23 of 3TM, where G 12 , G 13 , and G 23 are coupling parameters between subsystems. We show the 3TM can be approximated by an effective 2TM in the weak coupling limit combined with small q limit. This result clarifies a hierarchical structure of the MTM. The 3TM includes 2TM, and 2TM includes 1TM in appropriate limits of spacial, time and parameter scale.
An important application of the eigenmode picture is the Green's function method for the description of source term effects. For an arbitrary source term the general solution can be formulated by integrating a product of source term and the Green's function over space and time. We see the temperature inversion is sustained for ns scale, which can be understood from the + (q) mode's lifetime dispersion relation. We derive an analytical solution for a Gaussian type source term. This type of source term occupies a special position. From a solution approximated at small q limit, we can physically interpret the second order derivative of the eigenvalue �� (q) as an effective diffusion constant. Early attempts [28,29] to formulate analytical solutions can be compared with our work. Our formalism has wider applicability for cases like non-zero thermal diffusion constant, more-than-two spatial dimension, and N temperature model We finally derive a series of theorems which strongly restrict the eigenvalue properties of the MTM. According to the theorems the MTM eigenvalue is always non-positive real valued, which assures the mode lifetime is always well defined. A stationary solution is also guaranteed to present only in q = 0 point. The results of this paper will serve as a foundation of advanced models, e.g., with non-linearity or spatial non-uniformity. This paper is organized as follows. In Sect. 2, we introduce the 2TM and derive its exact solutions. We examine the 2TM behaviors for the case of bulk gold. In Sect. 3, we introduce the 3TM and derive its exact solutions. The mode property of an electron-photon-spin system is analyzed. A derivation of an effective 2TM from the 3TM is also discussed. In Sect. 4, we describe the inclusion of a source term by Green's function approach. In Sect. 5, we derive a series of theorems which provide strong and physically reasonable limitations on the MTM eigenvalue properties. Conclusions are given in Sect. 6 2 Linear two temperature model The linear two-temperature model is defined as, Here T e (t, r) and T l (t, r) are the electron and lattice temperature at position r and time t. C e (C l ) and e ( l ) denote the heat capacity and the thermal diffusion coefficient of (1) the electronic (lattice) subsystem, respectively, and G the electron-lattice coupling constant. Parameters C e , C l , e , l , and G are all positive real valued. Throughout this paper, we assume that the system is rectangular shaped whose side lengths are given by L i (i = x, y, z), and use a boundary condition is an eigenfunction of the diffusion term, where the wave vector is defined by q n x n y n z = (q n x , q n y , q n z ) = ( n x ∕L x , n y ∕L y , n z ∕L z ), with n i (i = x, y, z) being a non-negative integer. For simplicity, we omit the subscript from the wave vector hereafter. Thus, the general solution of Eq. (1) can be expressed by a linear combination of different wave vector components of the form as a natural extension of the Fourier series expansion common in the studies of thermal diffusion [8,30,31]. Because of the spatial uniformity, each q component is independent.
Then, the coefficients A q B q and (q) are the eigenvectors and eigenvalues, respectively, of a 2 × 2 non-symmetric matrix, where e (q) = − e C e q 2 , l (q) = − l C l q 2 , with q = |q|, and Ω e = G∕C e , Ω l = G∕C l . It is interesting to notice that the analytical form of H ′ is analogous to the Hamiltonian of other physical systems such as the quantum Rabi model and the polariton model except H ′ is not symmetric. The eigenvalue of H ′ splits into the upper and lower modes (q) = + (q) and − (q), respectively: where The corresponding right eigenmode (eigenvector) is given by where (2) r i T e = 0, r i T l = 0 for r i = 0, L i (i = x, y, z).
The general solution Eq. (3) of Eq. (1) is given by, Since the linear temperature model assumes a spatially uniform system, different wave vector components do not couple with each other. The mode amplitude A q is determined by the initial condition. It follows from e (q = 0) = l (q = 0) = 0 that, and, One can show that e ± > 0 and l ± > 0, therefore, ± monotonically decreases with increasing q (see Fig. 1a below). Furthermore, + (q = 0) = 0 [Eq. (12)], and, otherwise, ± < 0, indicating that all the modes damp except for v R + 0 , which corresponds to the final state; the larger the wave number, the faster the mode damps on each branch. It should also be noted that each individual eigenmode except for v R + 0 cannot be a physical solution alone, since it spatially oscillates around zero. The general solution must be a superposition of two or more modes to ensure non-negative temperature everywhere in the system.
Asymptotic behaviors of the solution in small and large q limit are informative to see the nature of this model. For q → 0 limit, and We see that the upper mode + reduces to an effective "one temperature model" with the effective heat capacity C eff = C e + C l and the effective thermal diffusion coefficient Up to the q 2 order, only the relative amplitude Eq. (15) between the lattice and electronic systems provides the information of the electron-lattice coupling G. Contrary to the upper branch, the lifetime of the lower mode −1∕ − (q = 0) = 1∕(Ω e + Ω l ) enables us to determine the value of G in its leading term.
Next, we examine the large q limit, which corresponds to e (q), l (q) ≫ Ω e , Ω l . For the upper branch: and for the lower branch: In this limit, the electron-lattice coupling G is negligible, and the system dynamics is dominated by "free diffusion process". Now, as a specific example, let us investigate the behaviors of the modes for the case of gold. We set C e = 7.77 × 10 5 J∕m 3 K and e ∕C e = 9.99 × 10 −3 m 3 ∕s (see Appendix A for details). We referred the literature [20,32] for the value of electronlattice coupling as G = 3.5 × 10 16 J∕m 3 K s. The lattice heat capacity C l = 2.4 × 10 6 J∕m 3 K and the lattice thermal conductivity l = 2 J∕m K s are taken from the literature [33]. In Fig. 1a, we can confirm that all eigenvalues ± (q) are negative real valued, monotonically decreasing with q, and hence the lifetime of each mode ± = −1∕ ± is well defined except for q = 0 where + = 0 . Once the initial condition is given, the system dynamics is completely described by the damping process of each mode. We find that the asymptotic solution Eq. (18) reproduces 87% of the exact value for wave length = 2 ∕q = 0.5 μm and Eq. (20) gives 101% for = 1.0 μm. The relative amplitude (Fig. 1b) shows a qualitative difference between the upper and lower branch. In the upper mode ( + ) , the electron and lattice temperatures spatially oscillate in phase, while in the lower, or − mode the oscillation is antiphase. Figure 1b also shows that in the large q limit the amplitude of electron (lattice) temperature in the upper (lower) mode vanishes, which indicates a transition to the free diffusion process. We can also see this transition in Fig. 1c, which plots the dependence of lifetime ± ; the exact solutions Eq. (5) approach to the asymptotic solutions Eqs. (18) and (20) in the small , i.e., large q, limit. Figure 1c indicates that such a transition occurs at tenth of nanometer scales in the upper mode and at sub m scale in the lower branch. The mode profile and lifetime of + (q) mode given in Fig. 1b, c lead to a consequence-the lattice temperature can be higher than the electron temperature. According to our knowledge, such a feature has been overlooked in the decades-long history of 2TM. We will see in the Sect. 4 that such a phenomenon we call "temperature inversion" can be realized with a conventional setup of an external source term.

Linear three temperature model
As a natural extension of the 2TM, the three temperature model (3TM) is defined as follows: Wave vector q dependence of a the eigenvalue of two-temperature model Eq. (5), b the relative amplitude of lattice temperature given in Eq. (8). ± in small |q| limit Eqs. (14) and (16) are shown by dotted and dashed lines, respectively. c Wave length = 2 ∕q dependence of the lifetime ± = −1∕ ± , where small limit Eqs. (20) and (18) are shown by dotted and dashed lines, respectively. We use parameters for gold given in literature [20,32,33] (see text) Here, T is a three component vector, representing the temperatures of the three subsystems. Λ and H are symmetric 3 × 3 matrices given by, We note again that the heat capacity C i , thermal diffusion coefficient i , and the coupling constant between subsystems G ij are all positive real valued. The matrices Λ and H are thus both real valued and symmetric 3 × 3 matrices.
Let us seek for the solution T of Eq. (22) expressed as a linear combination of different modes similar to Eq. (3). Then, we find three modes 1 (q), 2 (q), and 3 (q) of eigenvalues of a non-symmetric matrix H � = Λ −1 H, by using the formula for the roots of the general cubic equation, as, We show the dispersion relations of the eigenvalues of the three modes Eqs. (26)- (28) in Fig. 2a, specifically using parameters of nickel thin film deposited on a substrate at room temperature [24,34], namely, C e = 0.17 × 10 6 J∕m 3 K, C phonon = 3.50 × 10 6 J∕m 3 K, C spin = 0.20 × 10 6 J∕m 3 K for heat capacity of electron, phonon, and spin subs yst e m s , r e s p e c t i ve ly, G ep = 1.06 × 10 18 W∕m 3 K, G es = 0.10 × 10 18 W∕m 3 K for electron-phonon and electron-spin coupling, respectively, and e = 52.7 W∕m K for electron thermal conductivity. The rest of parameters, phonon , spin , spin-phonon coupling G sp are all set to be zero. We see from Fig. 2a that the eigenvalues are negative real valued or zero in this parameter range, as we demonstrate for general linear multiple temperature models below in Sect. 5. The mode lifetimes are, therefore, well-defined (Fig. 2b). Figure 2c-e show the relative weight of electron, phonon, and spin components, i.e., the eigenvector of matrix H ′ for (c) 1 , (d) 2 , and (e) 3 mode. Let us compare the mode characters of the 3TM with those of the 2TM (Fig. 1). The 1 mode has a vanishing eigenvalue or, equivalently, an infinite lifetime at q = 0, and the three components are equally weighted there (see Theorem 2 below). This observation indicates that the 1 mode can be regarded as a counterpart to the + mode of the 2TM. The 3 mode is, by contrast, dominated by the electron component (Fig. 2b) and, at large q limit, becomes purely electron temperature diffusion mode, with vanishing phonon and spin components. This feature is analogous to the − mode of 2TM. Figure 2b indicates that sub-picosecond dynamics of nickel 3TM can mainly be described by the 1 mode. Lastly, in 2 mode (Fig. 2c), the spin component dominates the 2 mode, and the weights of electron and phonon component vanish at large q limit (Fig. 2c). It is noteworthy that the lifetime of the 2 mode is almost dispersionless, nearly independent of q (Fig. 2b), suggesting that the thermal equilibration is achieved mostly by the energy exchange between subsystems while the diffusion term plays a negligible role. This is counterintuitive since the spin subsystem is coupled to the electron subsystem, in which rapid thermal diffusion can take place in principle. For this feature, the 2 mode seems to have no counterpart in 2TM. We note that the mode character of 3TM depends on the specific choice of parameters. For example, one can find a mode in which one of the three components has exactly no weight, even if the three subsystems are strongly coupled to each other. Such mode appears when the parameter set is symmetric, i.e., C 2 = C 3 , G 12 = G 31 , and 1 ≠ 0, 2 = 3 = 0. This is an analog of the dark state found in a three-level system driven by an external field [35].
We have found in the previous section that the "effective one temperature model" is embedded in the linear 2TM. Then, a question may naturally rise asking under which conditions the linear 3TM can be approximated by an effective two-temperature model. This question is of practical importance, since, for example, the electron-phonon-spin system was treated with both the 3TM [24] and approximate 2TM [36]. In what follows, we demonstrate that both of the two conditions: (i) small q limit and (ii) strong coupling limit G 23 ≫ G 12 , G 13 must be satisfied for such approximation and present its lowest order correction [Eqs. (59)-(64)].
Since the exact solution Eqs. (26)- (28) is too complicated to handle by a simple power expansion, We put a start point on a weekly coupled 1 + 2 temperature model, where matrix H ′ is decomposed to, Equation (39) is a block diagonal matrix describing a decoupled 1 + 2 temperature model, whose eigenvalues are Parameters are taken from the literature [24,34]. b Corresponding wavelength dependence of the mode lifetime. c-e Relative weight of electron, phonon, and spin temperature components for c 1 , d 2 , and e 3 mode, normalized so that the norm of the right eigenvector is unity. Insets in panels d and e: close-up of the zero-weight region simply given by (q) = 1 (= 1 ), ± , where ± is defined by Eq. (5). The two modes 1 (q) and + (q) are degenerate at q = 0 as 1 (q) = 0 and + (q) = 0. This is in contrast to the spectrum in Fig. 2a for three subsystems coupled with equal strength, where only 3 vanishes at q = 0 and 1 and 2 are degenerate there. As long as we restrict the timescale to t ≫ − = −1∕ − , we can neglect the contribution of the − branch. Then, the corresponding right eigenvectors (eigen- The amplitude A mn , or transformation matrix, is determined by solving a following eigenvalue equation: where the elements of 2 × 2 matrix K is given by or, explicitly, and m is an eigenvalue. By taking a small q limit and omitting terms smaller than O(G −1 23 ), we obtain where and the matrix elements of J(q) is given by, J(q) is the lowest order correction in large G 23 limit. Finally, we replace q by ∇ and reformulate Eq. (54) as an effective 2TM: with the effective parameters G eff , C eff , and eff given by, and the lowest order correction terms, (59) The appearance of the ∇ 2 -dependent correction terms owes to the deviation of the + (q) mode from a parabolic dispersion at large q.

Green's function approach to the 2TM with a source term
We now extend the 2TM to include the source term as follows: where Q e (t, r) and Q l (t, r) are source terms which deposit energy at position r and time t on electron and lattice subsystems, respectively. A general solution of this problem can be expressed by using the Green's function G ij (t, r;t � , r � ) as: Here, i and j are indices for subsystems, and the sum with respect to j runs over all the subsystems. The integral with respect to r ′ is taken over the system volume V. T 0 (t, r) is a solution of the homogeneous differential equation, i.e., for the source-free case: and can be expressed as a linear combination of the right eigenvectors of the 2TM, as has been discussed in Sect. 2. The Green's function G ij (t, r;t � , r � ) is constructed from the right and left eigenvectors of 2TM as, which satisfies with I ij being the 2 × 2 identity matrix. When the system is assumed to be a rectangular parallelepiped extending over a (64) finite range [0, L i ] along the i (= x, y, z) axis, the left eigenvector v L mq (r) is given by so that v R m q (r) and v L m q (r) fulfill the orthonormality condition Since we consider a finite volume, a discrete sum with respect to wave vector q is taken in Eq. (68). The index m runs over the two eigenmodes m = + and − . (t) is the step function which satisfies that (t) = 0 for t < 0 and that (t) = 1 for t > 0. In the following, as a specific example, we consider a three-dimensional metallic system irradiated by a Gaussian shaped laser field. We assume a rectangular-shaped system extending over x, y ∈ [0, L] and z ∈ [0, L z ]. When the laser field is incident on the surface at z = 0, the source term Q(t, r) can be expressed as follows: where denotes the beam radius, d the penetration depth, and Δt the temporal width. The beam axis is set to be (x c , y c ) = (L∕2, L∕2). At t = t c the intensity reaches its peak. The prefactor has been chosen so that Eq. (72) satisfies, Q e thus provides a total energy absorbed by the system when ≪ L and d ≪ L z . In this case, we can calculate Eq. (66) as, where the time-dependent mode amplitude A m q (t) induced by the source term is expressed as a function of , d, and Δt as, T e (t, r) T l (t, r) =T 0 (t, r) being the complementary error function.
In Fig. 3, we present numerical examples for a system with L = 100 μm and L z = 1.0 μm, and a source term with = 10.0 μm, d = 0.04 μm, Δt = 1.0 ps, and t c = 0.0 ps. We use the same material parameter values C e , C l , e , l , and G as in Sect. 2. Figure 3a, b display the real-space snapshot of the electron and lattice temperature distribution, respectively, 1.0 ps after the peak of the external field intensity. The corresponding distributions of A m q (t)e m (q)t in the q space are shown in Fig. 3c, d. We note that the highly spatially oscillating factor cos(q x x c ) cos(q y y c ) is excluded here to focus on the envelope structure.
The extension to multiple pulse excitation is straightforward thanks to the superposition principle. If the source term is a train of pulses of the same form, the solution is given by,   = (0, 0, 0)). The pulse interval is set to be 1.5 ns. The excitation condition of each pulse is the same as in Fig. 3. The inset shows the short time profile of the temperature dynamics where t (i) c denotes the peak time of the i-th pulse, and T(r, t) is given just by the second term of Eq. (116) with t c = 0. Thus, it suffices to shift the common solution T(r, t) by t (i) c with no need to repeat calculation for each pulse. As an example, in Fig. 4, we show the temperature dynamics of bulk gold system excited by multiple laser pulses of 1.4 ns intervals, each of which has the same temporal profile as in Fig. 3. The temperature inversion, which we anticipated in Sect. 2 that T l > T e , is evident in the ns scale dynamics in Fig. 4. Figure 3d shows the amplitude of + (q) mode extends over sub-μm −1 scale along q z direction. Figure 1c indicates the mode lifetime is nanosecond order in the corresponding wave length, which is fully consistent with our findings in Fig. 4. In most 2TM simulations the short time dynamics in picosecond time scale matters and its long time nature has been paid less attentions. In addition, the temperature gap |T e − T l | in temperature inversion is less significant than that of the early stage dynamics. Possibly for those reasons the temperature inversion has been overlooked. We, however, point out that the fully equilibration of electron and lattice temperature in 2TM requires nanosecond order time scale. This is much longer than picosecond time scale which is usually expected [20] or estimated by − (q) mode lifetime −1∕ − (0) = 1∕(G(1∕C e + 1∕C l )). Equation (15) clarifies that it originates from the difference of two parameters e ∕C e and l ∕C l . The temperature inversion emerges when e ∕C e > l ∕C l , namely the diffusivity of the electronic system is larger than that of the lattice system, which seems inevitable in metallic systems. We here highlight that the temperature inversion is inherent even in our simplest model whose parameters are all temperature-independent. For same reasons Fig. 2c indicates the emergence of the temperature inversion in electron-phonon-spin 3TM, while its Green's function based analysis is out of our scope.
The Gaussian shaped source term does not only express a typical laser intensity distribution but also provides an insight to physical interpretation of the eigenmodes. For that purpose we consider a quasi-two-dimensional system extending over x, y ∈ [0, L]. The source term can then be expressed by Further we take limit as, and impose two conditions: (78) L → ∞ With condition Eq. (78), we replace the discrete sum on the q grid by an integral on the q space. Equation (79) states that the pulse width Δt is sufficiently shorter than the mode lifetime, which is satisfied when the mode amplitude is localized around q = 0. Hence, we truncate Taylor series in terms of q up to the second order, so that we can calculate the integral analytically to obtain, assuming an isotropic system, where, with, and q = q 2 x + q 2 y . We see from Eq. (80) that the spreading of temperature distribution increases proportionally to time. Interestingly, this observation indicates that �� m (0) plays the role of an effective diffusion constant. See Appendix B for the detailed derivation. The three dimensional case is also discussed in Appendix C. Whereas we have presented a specific example for the two-temperature model here, the Green's function formalism and the physical interpretation of �� (q) can be straightforwardly extended to general multiple temperature models with three or more subsystems.
where T now denotes the N-components vector representing the subsystem temperatures, and the subscripts i, j run from 1 to N. Examples of the MTM include the nonthermal lattice model [21] or just multitemperature model [37], which assign phonon mode resolved temperatures. In the previous sections we have found that the spectra of linear 2TM is always negative real valued, or exactly zero at q = (0, 0, 0) point. The linear 3TM shows the same property within the parameter range we plot in Fig. 2a. Here we prove that this physically reasonable property holds for any N, assuring that the temperatures of all the subsystems asymptotically tend to a common, spatially uniform, final value.

Theorem 1 When the boundary condition Eq. (2) is given, the linear MTM Eq. (84) is transformed as
Once the initial condition is given, we can completely determine the MTM dynamics from the eigenvalue of a matrix H ′ can be diagonalized and its eigenvalues (q) satisfy the following two properties: 1. the eigenvalue (q) of matrix H ′ always satisfies (q) ∈ ℝ and (q) ≤ 0, 2. When (q) = 0, q always satisfies q = 0.
Proof We consider the following eigenvalue equation is a right eigenvector and is a corresponding eigenvalue. By multiplying both sides by a diagonal matrix Λ 1∕2 from the left, which satisfies (Λ 1∕2 ) 2 = Λ, we obtain Thus Λ 1∕2 v becomes an eigenvector of a symmetric matrix Λ −1∕2 HΛ −1∕2 whose eigenvalue is given by . Clearly the always satisfies ∈ ℝ.
We can further restrict the distribution of eigenvalues on the complex plain by using the Gershgorin's theorem [38], which states that the eigenvalues of N × N matrix A exist on a closed region D which is defined by where C i (i = 1, … , N) is a closed disk whose center position is given by A ii on the complex plane and its radius R i is given by In our case, the center position of the closed disc C i is given by and the radius R i of C i is since G i and C i are positive real valued parameters. The closed region D therefore extends over a semi-infinite plain whose real part is negative, and D can include the origin of complex plain only if q = 0. We therefore conclude that (q) is always non-positive real valued and can be 0 only if q = 0. ◻ This theorem strongly restricts the behavior of linear MTM. For any given initial condition, the linear MTM only provides damping solutions regardless of material parameters. Consequently, an external heat, or maybe nonlinearity is needed to excite oscillatory and amplifying behavior in its dynamics.
In addition, we can show the following property of the MTM's solution at q = 0. Theorem 2 When q = 0, at least one eigenvalue of Eq. (89) becomes 0, and at least one right eigenvector of such solutions has all its components equal.
Proof We firstly prove the first half of the theorem. The matrix H given by becomes linear dependent when q = 0, i.e., det H = 0 at q = 0. The matrix H � = Λ −1 H then becomes linear dependent: det H � = det Λ −1 det H = 0 at q = 0 limit, as well. This implies a condition Here i (q) (i = 1, … , N) is an eigenvalue of H ′ . To satisfy Eq. (98), at least one i (q) must fulfill a condition We secondly prove the latter half of the theorem. By substituting a right eigenvector with equal components: we can show the following: ◻ It should also be noted that Theorems 1 and 2 jointly assures that the temperatures of all the subsystems approach to a common, spatially uniform, finite value in the long time limit.
From a perspective of the total energy conservation, the stationary solution (q = 0) = 0 plays a special role. The total energy U total (t) of the system at time t is given by where the sum runs over all the subsystems, and the integral is taken over the system volume V. Since the temperature distribution T i (r, t) can be decomposed into the contribution of each eigenmode, we can define the eigenmode resolved energy U q (t) by where v R q,i denotes the i-th component of the right eigenvector of the MTM Eq. (89). The U total (t) is retrieved by summing up U q over the mode index and wave vector q: In the absence of the external heat source, the total energy of the system Eq. (102) must be conserved. The following theorem guarantees this requirement and, moreover, shows that all the eigenmodes but (q = 0) = 0 do not hold net energy.
Proof The case of q ≠ 0 is trivial since its spatial dependence is sinusoidal, which becomes zero when integrated over the system volume. We next show the case of (q = 0) ≠ 0.
We can then write down the MTM Eq. (84) as By using the property of H for q = 0: we can derive the following: Since (q = 0) ≠ 0, U q must be zero. The total energy Eq. (102), therefore, has a finite contribution only from the stationary mode = 0, q = 0. ◻ We finally show a subsidiary theorem about the monotonically decreasing property of the eigenvalue (q) with respect to the magnitude q of wave number. Proof It is sufficient to prove it for the symmetric matrix H = Λ −1∕2 HΛ −1∕2 as it possesses same eigenvalues with the MTM's matrix Λ −1 H. In this case the right eigenvector v coincides with left one. We can then immediately write down as follows: (107) In the third line of Eq. (110), we assumed the norm conservation of v: ◻ This theorem physically states that the larger the wave vector q of the mode, the shorter the mode lifetime (q) = −1∕ (q), or equivalently, the faster the mode damps.

Conclusions
We have presented the exact analytical solutions of the linear multiple temperature model under the insulated boundary condition. By extending the familiar Fourier series expansion of the heat equation, we have shown that the system dynamics is expressed as a linear combination of eigenmodes with different wave numbers and lifetimes (or decay constants).
The dispersion relations of the mode lifetime and the weight of subsystem components have been investigated for the 2TM and 3TM which describe bulk gold and nickel thin film, respectively. The 2TM and 3TM have some common features such as a mode whose lifetime diverges at q = 0 and asymptotic convergence to free diffusion modes at large q limit. There is, however, a mode unique to 3TM which is almost dispersionless despite a strong coupling between subsystems. Analyzing the mode characters thus helps to deepen our understanding of the temperature dynamics.
The eigenmode picture has enabled us to unveil the hierarchical structure of the MTM that an N + 1-temperature model approximates an N-temperature model in certain limits. For example, in the small wave number q and long time limits, the upper mode + (q) of the 2TM approximates the 1TM. If the weak coupling limit is additionally taken, the 3TM approximates the 2TM. We note that such approximations break in a small spacial scale. In the opposite, large wave number limit, on the other hand, the diffusion term dominates the dynamics and the electronlattice coupling is negligible. A useful application of the eigenmode picture is the Green's function method for the source term. It leads to our finding of the temperature inversion, which shows a counterintuitive behavior T e < T l and hinders fully equilibration of the system for ns time scale. We have derived the analytical solution for a Gaussian-type source typical of a laser profile. It is straightforward to extend this formulation to a case of more complex source terms such as multiple pulses, as well as to three or more temperature models.
We have shown how physical requirements are implemented in the eigenmode picture of the linear MTM. Each eigenmode neither grows nor oscillates, but exponentially decays, except for a spatially uniform, stationary mode corresponding to the final state where all the subsystems has the same temperature. The larger the mode wave vector, the shorter its lifetime. Only the stationary mode has a net finite energy, assuring the conservation of the total energy.
The present linear MTM can be further extended to better describe real materials by including non-linearity and spatial non-uniformity, etc. Even in such a model, the eigenmode picture will be a useful tool to interpret the system dynamics through, e.g., mode-mode interaction, and mode scattering processes, as is done in quantum mechanics and nonlinear optics. The temperature inversion in nonlinear system is of our great interest, which potentially serves as a useful tool to characterize the role of non-linearity.

Appendix A: Model parameters for Green's function method
Although the electron specific heat and thermal conductivity are temperature-dependent in general, we concentrate on the linear 2TM in this study and, thus, take C e and e as constant. In this Appendix we describe a method to estimate the effective (or average) constant values of C e and e .
Let us start from the temperature dependent specific heat and thermal conductivity C � e (T e ) and � e (T e , T l ) given by, respectively, with = 67.6 J∕m 3 K 2 , v F = 1.39 × 10 6 m∕s, A = 1.2 × 10 7 K −1 s −2 , and B = 1.23 × 10 11 s −1 K −1 according to literature [20,32] for bulk gold. C � e (T e ) ≡ U e Te by definition, and Eq. (112) is derived from the electron internal � e (T e , T l ) = energy U e (T e ) = T 2 e ∕2 aside from T e independent terms. On the other hand, the constant electron specific heat implicitly assumes, with being constant. In this study we determine the value of C e , so that it may reproduce the internal energy increase when T e is elevated from T (1) e to T (2) e , i.e., C e T (i) e + = T (i)2 e (i = 1, 2), where T (1) e and T (2) e are typically the room temperature and the highest temperature achieved. In Sect. 4 we use C e = 7.77 × 10 5 J∕K m 3 corresponding to T (1) e = 300 K and T (2) e = 2000 K, the latter of which is consistent with the maximum electron temperature achieved (2093 K).

Appendix B: Derivation of the analytical solution Eq. (80) on the real space
In this section we provide a detailed derivation of the analytical solution Eq. (80) on two-dimensional real space with the source term Eq. (77). We start from the expression of the solution in two-dimensional system: where the mode amplitude A m q x q y (t) is given by We ignore a term providing initial condition for simplicity. The right eigenvector v R m q in Eq. (116) is given by We expand the normalization factor f m q in Eq. (117) up to the second order in q as, where due to the system isotropy. We factorize Eq. (117) by using Eq. (120) as, We further expand the mode eigenvalue m (q) in Eq. (121) up to the second order in q as, where, (116) T e (t, x, y) T l (t, x, y) = ∑ m,q x ,q y A m q x q y (t)e m (q)t v R m q (x, y),  again due to the system isotropy. For evaluating the sum over q in Eq. (116), we firstly examine the discontinuous factor at q i = 0 (i = x, y) in Eq. (121). We transform as: and substitute it in Eq. (121). By taking the sum over q vector in Eq. (116), the contribution from the second term of Eq. (124) turns out to be proportional to 1∕L i for i = x, y. By taking a limit L i → ∞ this contribution vanishes, thus the discontinuity at q i = 0 does not make any problem. By taking a limit L i → ∞ (i = x, y) , we replace the sum in Eq. (116) by the integral about q as: where and r c,i = x c , y c for i = x, y, respectively. Δq i = 2 ∕L i is the interval of a discrete q i grid. We also use conditions − m (q)Δt

Appendix C: Analytical solution on the bulk surface
In this Appendix we consider a three-dimensional system and derive an analytical solution at z = 0, which corresponds to the temperature distribution on the bulk surface. We first consider the amplitude Eq. (75) for d ≪ L z case. As was done in Appendix B, taking a limit L z → ∞ enables us to replace a sum over discrete q z by an integral. We then need to evaluate an integral: