Phononic heat transport in the transient regime: An analytic solution

We investigate the time-resolved quantum transport properties of phonons in arbitrary harmonic systems connected to phonon baths at different temperatures. We obtain a closed analytic expression of the time-dependent one-particle reduced density matrix by explicitly solving the equations of motion for the nonequilibrium Green's function. This is achieved through a well-controlled approximation of the frequency-dependent bath self-energy. Our result allows for exploring transient oscillations and relaxation times of local heat currents, and correctly reduces to an earlier known result in the steady-state limit. We apply the formalism to atomic chains, and benchmark the validity of the approximation against full numerical solutions of the bosonic Kadanoff--Baym equations for the Green's function. We find good agreement between the analytic and numerical solutions for weak contacts and baths with a wide energy dispersion. We further analyze relaxation times from low to high temperature gradients.


I. INTRODUCTION
Heat in a nanoscale junction is transported by electrons and molecular vibrations known as phonons or vibrons. In traversing an electronic device, an electron may transfer a portion of its energy to the structure's molecular vibrations thus creating a heat wavefront propagating through the system 1-3 . Therefore, phonons are ubiquitous in virtually any molecular or nanoscale junction. In spite of the progress and advances in studying electron transport phenomena in nanoscale structures since the early 1970s 4 , measurements of thermal conductance at this scale presents several technical and conceptual difficulties and became available only recently. Suspended nanostructures such as monochrystalline GaAs 5 and silicon nitride layers 6 have been among the first settings where quantized thermal conductance was observed. Further thermal conductance studies have been performed in carbon nanotubes 7,8 and silicon nanowires 9,10 , and ultimately reaching the scale of molecular and atomic contacts [11][12][13] .
A considerable amount of purely electronic transport simulations has been done according to the intuitive framework by Landauer 14 and Büttiker 15 . The basic idea is that the charge dynamics in a conducting channel connected to electrodes is governed by a scattering process where charge carriers are either transmitted between the electrodes through the channel, or reflected back to an electrode from the channel. From this approach an expression for the current through the channel is obtained in terms of transmission probabilities. This framework has been particularly useful after the more rigorous work by Caroli et al. 16,17 based on the nonequilibrium Green's functions. Further generalizations to the formalism are due to Meir and Wingreen 18 , Pastawski 19 , Jauho et al. 20 , Stefanucci and Almbladh 21 , Myöhänen et al. 22 who discussed electron interactions in the conducting channel, dissipation due to many-body effects, timedependent voltages and fields, and the role of initial correlations. In the case of noninteracting particles the same methodology has been proposed up to the derivation of a time-dependent Landauer-Büttiker formula in simple systems 21,23 and arbitrary structures with timedependent bias voltages [24][25][26] .
The Landauer-Büttiker formalism is not limited to electronic transport. It also applies in the context of thermal transport for phononic systems as proposed by Segal et al. 1 , Rego and Kirczenow 27 , Rego 28 via the wave scattering method (similar to the original work of Landauer) and by Yamamoto and Watanabe 29 , Wang et al. 30,31 in terms of nonequilibrium Green's functions (similar to the work of Meir and Wingreen). In both approaches, expressions for the thermal conductance and the heat current are recovered as a Landauer-Büttiker formula. These studies have further been supplemented by, e.g., exact transmission functions 32 , heat transport models for refrigeration 33 , nonlinear interactions 34 , and thermal conductance studies in atomic junctions [35][36][37][38] , graphene ribbons [39][40][41] and carbon nanotubes 29,42 .
The nonequilibrium Green's function approach offers a natural framework to study transient effects both in electronic [43][44][45] and phononic systems [46][47][48] . Also in the case of phonon transport, it would be extremely useful to have a time-dependent Landauer-Büttiker formula for interpreting the transient oscillations and relaxation times in an intuitive fashion. In addition, as the time-propagation of the equations of motion for the phonon Green's function is computationally a rather heavy task 49,50 , bringing down the computational cost when studying larger systems would be advantageous. In this article, we present a substantial advance in this direction; we consider harmonic nanoscale systems connected to heat baths in a wide-band-like approximation, and derive a closed expression for the time-dependent phonon density matrix. In the steady-state limit our result reduces to a form of Landauer-Büttiker type discussed in, e.g., Refs. 51 and 52.
We will first introduce the model and the dynamical equations of motion in our phonon transport setup (Sec. II A), and then we will discuss the embedding scheme with its approximations and limitations (Sec. II B). The solution to the equations of motion, the main result of the present work, is derived in Sec. II C with more detailed calculations presented in the appendix. As an illustration of the derived formula we study the transient heat currents in small atomic chains (Sec. III). A summary and main conclusions are drawn in Sec. IV.

A. Transport setup and assumptions
We model heat transport in a nanomechanical device coupled to phononic baths at different temperatures. The description is for noninteracting phonons only. The Hamiltonian for this setup can be written in terms of momentum and displacement field operators ( = 1) with indices j, k = 1, . . . , N running over the basis of the studied system, m j being the mass of the j-th oscillator, and K being the positive definite force constant matrix. The momentum and displacement operators obey the canonical commutation relations [û j ,p k ] = iδ jk , [û j ,û k ] = 0 = [p j ,p k ]. We define mass-normalized operators asû j = √ m jûj andp j =p j / √ m j obeying the same commutation relations, and write further where also the force constant matrix was transformed as K jk = K jk / √ m j m k . The second equality follows by introducing the composite operatorsφ with componentŝ φ 1 j =û j andφ 2 j =p j . The indices {µ, ν} ∈ {1, 2} run over the different components of the field operatorsφ. The matrix elements of the block matrix Ω jk are given by Ω 11 jk = K jk , Ω 22 jk = δ jk and zero otherwise. The motivation behind theφ-"spinor"-representation is that the equations of motion are first order differential equations, instead of second order 49,50 . The canonical commutation relations are encoded in the field operatorsφ as We accordingly define the phononic Green's function for time arguments z and z on the Keldysh contour 53 γ as where T γ is the contour-time ordering operator and · is an ensemble average 49,50,54 . The equations of motion for the Green's function can be expressed through the time evolution of the field operatorsφ, and they read as (matrices in the 2N × 2N representation are from now on denoted with boldface symbols) Here each term is a 2 × 2 matrix. Let us look more specifically at the transport setup shown in Fig. 1. The full Hamiltonian can be expressed as a composition of three parts: the left reservoir (L), the central system (C), and the right reservoir (R) The different subsystems are coupled apart from a direct coupling between the reservoirs. In principle we have an underlying potential energy surface E( R) for all nuclear coordinates R in the composite system which, in turn, determines the force constant matrix elements by K jk = ∂ 2 E/∂R j ∂R k . For an uncoupled system we would only have the diagonal blocks nonzero in Eq. (7).
Here we consider only two reservoirs but the partitioning in Eq. (7) can directly be generalized to an arbitrary number of reservoirs to include, e.g., inner reservoirs as probes 51 . In Eq. (7) each diagonal block is a 2N C × 2N C or 2N λ × 2N λ matrix for λ ∈ {L, R}. The block structures for the diagonal elements are simply those discussed earlier The different regions, however, couple only through the displacement term, so the block structures for the offdiagonal elements are given by ( FIG. 1. (Color online) Heat transport setup where a central system of interest is connected to two reservoirs of different temperatures. The internal structures and couplings are defined by the force constant matrices K. The coupling refers to partitioned approach where the different systems are uncoupled at times t < 0 and coupled at times t ≥ 0 when the system is driven out of equilibrium.
We are mainly interested in the transport properties of the central system, so we extract the component corresponding to the central region, CC, from the equations of motion (5) and (6). This procedure leads to the following set of equations (15) and the corresponding adjoint equations. In the equations above α CC ≡ α ⊗ 1 CC and α λλ ≡ α ⊗ 1 λλ are 2N C × 2N C and 2N λ × 2N λ matrices, respectively. Here d λλ is the isolated phonon Green's function in the reservoir λ satisfying the equation of motion (15) with the reservoir Hamiltonian Ω λλ . Since we are dealing with noninteracting phonons, the self-energy Π CC in the collision integral (13) is given solely as the embedding selfenergy defined in terms of the coupling Hamiltonians Ω Cλ in Eq. (14). From Eq. (12) and its adjoint we can derive an equation of motion for the lesser Green's function D < CC by using the Langreth rules for the collision integrals 53 . Particularly, we are interested in the timedependent one-particle reduced density matrix given by the time-diagonal ρ(t) = iD < CC (t, t). This is given by Keldysh functions a and b, and we also dropped the subscripts CC as we will only refer to the central region from now on 55 . The Keldysh components lesser (<), greater (>), retarded (R), advanced (A), left ( ), and right ( ) of a function k(z, z ) on the contour are defined by 53 with the contour points z = t − on the forward branch, z = t + on the backward branch, and z = t 0 − iτ on the vertical branch.
If we did not have the central system embedded into the environment, the self-energy terms in Eq. (16) would simply be zero, and we would be left with a Liouville-type equation for the time-evolution of the reduced density matrix for an isolated central region. The self-energy terms, therefore, account for the open transport setup where a finite central region is embedded into the environment.
So far, the discussion has been rather general, and Eq. (16) also applies to many different setups beyond the present study. For instance, interactions could be included by adding a many-body self-energy contribution. Our aim is to solve (analytically) this integro-differential equation for D < (in the equal-time limit) and then extract dynamical quantities such as heat currents from the time-dependent phonon density matrix. Similar derivations have been done in earlier studies 21,23-25 for electron transport, and we will work along these previous ideas. For solving the equation we need to make some approximations, the first one being the partitioned approach, i.e., all regions in the transport setup are initially (t < 0) uncoupled and in separate thermodynamical equilibria. At t = 0 we couple the different regions thus driving the system out of equilibrium, see Fig. 1. From the perspective of an underlying potential energy surface, which we discussed earlier, this would lead to nonzero off-diagonal blocks in Eq. (7). This coupling could, however, change the matrix elements in the diagonal blocks as well but we approximate them to be the same as for the uncoupled system.
The partitioning procedure disregards the initial couplings: Ω λC = 0 in equilibrium, so the integrations along the vertical track of the Keldysh contour on the righthand side of Eq. (16) are simply left out 46,47 . The remaining two terms on the right-hand side of Eq. (16) can be interpreted as a source/drain and a damping/equilibration term. The drain (source) term is a convolution between the propagator in the central region, D A (h.c.), and the lesser embedding self-energy Π < which is related to the probability of finding a phonon in the reservoirs, i.e., it describes the extraction (insertion) of phonons out of (into) the central region. The second term is a convolution between the propagator in the reservoirs, Π R (h.c.) and the lesser Green's function in the central region D < which is related to the probability of finding a phonon in the central region, i.e., it is responsible for damping (equilibration) effects.
For any type of reservoirs we may obtain Π R in a separate calculation. However, the complicated timestructure of the embedding self-energy makes it difficult to close the equations of motion into an analytically solvable form. Thus, we introduce an approximation. In the case of electronic transport, the commonly used wideband approximation makes the embedding self-energy proportional to a delta function allowing for a closed solution 24,25 . In the next section, we will make a similar approximation in the phononic case. B. Approximating the embedding self-energy Let us start by considering the coupling of the central region to the reservoirs via the embedding self-energy. In order to evaluate the embedding self-energy from Eq. (14) we need the coupling Hamiltonians and the reservoir Green's function, and, as we are considering the retarded component of the embedding self-energy, we need to find an expression for d R . For the isolated phonon Green's function in the reservoir λ we have the following expression which can be checked to satisfy Eq. (15) by calculating the derivative iα λλ ∂ z d λλ (z, z ). In the above expression f λ = 1 + f λ and f λ (ω) = (e β λ ω − 1) −1 is the Bose-Einstein distribution for reservoir λ at inverse temperature β λ = (k B T λ ) −1 . The expression in Eq. (17) is also proportional to the density matrix for an isolated system in the limit z → z + . By using the above expression for the uncoupled Green's function we may derive different Keldysh components and calculate the retarded embedding self-energy from Eq. (14), see App. A. We find that the real and imaginary parts of the retarded embedding self-energy are, respectively, even and odd functions in frequency ω. This finding also corresponds to the form calculated explicitly for a uniform one-dimensional system of N coupled springs 56 . We take this model for our reservoirs (N λ sites coupled with equal springs) and construct the embedding self-energy accordingly. In this model, the force constant matrix K λλ has diagonal elements 2k λ and the first off-diagonal elements −k λ . In the limit N λ → ∞ the retarded embedding self-energy is given by with ζ λ (ω) = sgn(ω + 2 √ k λ ). Even though we have an explicit form for the frequency dependency of the embedding self-energy in Eq. (19), we consider an approximation so that the equation of motion (16) could be solved analytically. The full numerical solution of Eq. (12) is also possible using a time-stepping algorithm in the twotime plane 49,50 but this is restricted to only small systems due to the computational cost.
We make a wide-band-like approximation to Eq. (18) so that for small frequencies (compared to other energy scales of the studied system) In contrast to the conventional wide-band approximation in electronic transport, now the retarded embedding selfenergy is not a purely imaginary constant. Instead, the real part is constant and the imaginary part is frequency dependent (linearized approximation). In fact, the imaginary part is not bounded when |ω| → ∞ leading to some technical difficulties addressed soon. Also, since the force constant matrices K are by construction positive definite, then from Eqs. (20) and (A11) we see that Λ 0,λ is negative definite, and from Eqs. (18) and (21) we get that Γ 0,λ is positive definite. Compared to Eqs. (20) and (21) similar wide-band-like approximations for the self-energy have been proposed in Refs. 57 and 58 where the embedding self-energy is approximated as a purely imaginary sign function. The real and imaginary parts of Π R λ (ω) are plotted together with the approximations in Eqs. (20) and (21) [evaluated from Eq. (19)] in Fig. 2 versus ω. [In the figure we also employ a cut-off frequency, see Eq. (23).] Looking at the form in Eq. (19) the imaginary part of the self-energy is nonzero only for |ω| < 2 √ k λ introducing the phonon bandwidth. The motivation behind the approximations in Eqs. (20) and (21) around ω = 0 is also clearly visible.
Since the coupling matrices between different regions only have nonzero components in the uu block, we express the embedding self-energy in the following form where each element of the 2×2 matrix is an N C ×N C matrix. Also, summing over the different reservoirs λ gives the total embedding self-energy: Π R = λ=L,R Π R λ , and this also applies to the real and imaginary parts Λ 0 and Γ 0 . Compared to the conventional wide-band approximation in electronic transport, where the imaginary part of the self-energy is constant and the real part is zero, here we run into problems with our approximation for the embedding self-energy because of the unboundedness of the imaginary part ωΓ 0 /2 when |ω| → ∞. This approximation for the real and imaginary parts of the embedding self-energy function does not satisfy the Kramers-Kronig relations 53,59 , and for this reason we introduce a cut-off frequency ω c,λ above which the approximation for the embedding self-energy is simply zero A natural choice for the cut-off frequency would be the phonon bandwidth, see Fig. 2. We may tune the considered frequency range by varying the reservoir force constant so that the important processes are captured around comparatively small frequencies ω/ √ k λ . The lesser component of the embedding self-energy is then simply given by the fluctuation-dissipation relation 53 The cut-off frequency is a similar concept as the Debye temperature; here they are related by ω c,λ = k B T λ 31 . Instead of fixing the cut-off frequency at the phonon bandwidth we could also determine it by a frequency sum rule.

C. Solving the equations of motion
Based on the approximation discussed in the previous section we derive expressions for the time convolutions in Eq. (16). More detailed calculations are shown in App. B and we state here only the results. First, we define the effective (nonhermitian) Hamiltonian as which reduces to the uncoupled Hamiltonian, αΩ, in the limit Λ 0 , Γ 0 → 0. For cases relevant to us 60 , this object has its eigenvalue spectrum in the lower-half plane, i.e., the imaginary part, given essentially by Γ 0 , of the eigenvalues are negative. The real part of the eigenvalues occur at the eigenfrequencies of the uncoupled system shifted by Λ 0 . For the adjoint Ω † eff , on the other hand, the eigenvalues lie in the upper-half plane.
For the convolution between the retarded Green's function and lesser embedding self-energy we have As the frequency integral is cut off, we evaluate D R in the limit ω c,λ → ∞ [see Eq. (B5)]: By conjugating Eq. (26), we also get where D A is found by conjugating Eq. (27). The convolution between the lesser Green's function and the retarded/advanced embedding self-energy is wellbehaved for all cut-off frequencies. We therefore take the limit ω c,λ → ∞ in Eq. (23), see App. B for more details, and obtain and by conjugating we also get The approximated form of the convolutions in Eqs. (26), (28), (29), and (30) are expected to be in good agreement with a full Kadanoff-Baym propagation provided that the vibrational frequencies of the central region fall in the frequency window where Im Π R (ω) is linear, see Fig. 2.
In Sec. III B we assess the accuracy of the approximations by evaluating Eqs. (26) and (28) with the exact D R/A , and Eqs. (29) and (30) with the finite cut-off Π R/A . An important thing to notice is that the the convolutions in Eqs. (29) and (30) depend not only directly on D < but also on the time-derivative of D < at equal-time limit. This means that inserting these expressions back into the equation of motion does not immediately close the equation for D < . However, we may insert iteratively the explicit time-evolution from Eq. (12) and its adjoint for the derivative terms. As Eq. (12) includes a similar structure of time-convolutions, we will gain similar terms by the iteration procedure. Interestingly, it turns out that the iteration is truncated already after the first cycle due to the structure of the transport setup, and we get a closed equation for D < . Now we are ready to insert Eqs. (26), (28), (29), and (30) into Eq. (16): Then, we insert ∂ t D(t, t ) and ∂ t D(t, t ) from Eq. (12) and its adjoint, and accordingly insert the consequent convolutions from Eqs. (26), (28), (29), and (30). This step generates a plethora of terms but we still expand all the parentheses to better see what is the overall structure of the various terms: Next step would be to again insert the derivatives from Eq. (12) and its adjoint, but already in the above step of the iteration, all terms involving derivatives of the lesser Green's function vanish. We notice this truncation (and cancellation of other terms as well) by evaluating simple matrix products All higher order derivatives are also truncated based on these matrix structures. By combining accordingly and simplifying we end up with In Eq. (34) we now have a linear, first order, nonhomogeneous differential equation for D < which can be solved uniquely with an initial condition. The solution is (see with the initial condition stemming from the uncoupled lesser Green's function in the central region as in Eq. (A2) where the distribution f C is defined via an equilibrium temperature for the central region before coupling. The spectral function Eq. (35) is our main result for the time-dependent oneparticle density matrix. Remarkably, it is a closed expression, i.e., no time propagation is needed for evaluating the time-dependent density matrix; this also is a general feature in previous studies including similar derivations 21,23-26 . The transient behavior is encoded in the exponentials: We find oscillations ω jk = | Re ω j,eff − Re ω k,eff |, where ω eff are the complex eigenvalues of the effective Hamiltonian Ω eff , as transitions between the vibrational modes in the central region . Finally, f λ (ω)B λ (ω) is well-behaving at ω = 0 (although f λ (ω) diverges at zero), and the cut-off frequency ω c,λ regulates the nonintegrable behavior at ω → −∞.
It is instructive to investigate few limiting cases for Eq. (35). At t = 0 the square brackets vanish and we are left with the uncoupled result, as should be the case due to the initial condition. This also happens if the systems remain uncoupled, i.e., Λ 0 = 0 = Γ 0 ; then we are left with the free evolution of the initial state as Ω eff → αΩ and B λ (ω) → 0. The steady-state result comes from the limit t → ∞ when the exponentials vanish due to the nonhermitian structure of Ω eff [see the discussion after Eq. (25)]: Within our self-energy approximation,  [65][66][67] , the resulting frequency integrals may be written in terms of complex logarithms and exponential integral functions.

A. Setup and quantities of interest
Here we apply the derived result in Eq. (35) to study the transient behavior of the heat current in simple lattice models when coupled to reservoirs at different temperatures. We benchmark the validity of the approximation by comparing to full numerical solution to the equation of motion (12) with the embedding self-energy in Eqs. (18) and (19).
As the derived result provides information about the one-particle density matrix in the central region, we are interested in local quantities in this region. The local energy in the central region may be calculated as a sum over the uu and pp blocks of the product of the Hamiltonian and the density matrix The local heat current between the sites of the central region may be derived by considering the temporal change in local energy in a given site; this should amount to the sum of heat currents flowing in and out of that site [68][69][70] . The local energy for site j can be written as an expectation value j = Ĥ j of the local Hamiltonian H j = [(p j ) 2 + kû j K jkû k ]/2. (This is chosen so that H = jĤ j .) Then, forĤ being the total Hamiltonian for the central region [see Eq. (2)], we get from the Heisenberg equation where we used the commutation algebra of the momentum and displacement operators. This motivates to define the local (net) heat current between sites j and k as the up component of the density matrix where the two terms can be regarded as "in-coming" (from k to j) and "out-going" (from j to k) heat current. This definition in Eq. (42) deviates a little from the conventional definitions for the heat current between a reservoir and the central region 33 since in our case there may be multiple (arbitrary) contacts between the sites contributing to the heat current in a given site.
We also note here that the above definition for a local energy (density) j is ambiguous since we can add any local contribution c j , for which j c j = 0, which leaves the total energy unchanged. Further, the definition for the local heat current, obtained from the temporal change in local energy, is also not unique although widely used [68][69][70][71][72][73] . The appropriate definition for the local energy density has recently been discussed in, e.g., Ref. 74. Our motivation here, however, is to compare the derived analytical result in Eq. (35) with a full numerical solution, and this issue does not affect the comparison.
In our transport setup, we have uniform onedimensional (semi-infinite) systems of coupled springs as reservoirs, and we fix the force constant in the reservoirs as k λ = 1 for all leads and then relate the remaining parameters to this energy scale. Also, the central region, through which we study the heat current transients, is similarly a uniform one-dimensional (but finite) system of coupled springs for which we have the inter-atom force constant k C : The terminal sites of the central region are coupled to the terminal sites of the reservoirs by the coupling force constant k λC . As discussed in Sec. II B, we can tune the embedding by choosing the coupling strength k λC . Also, for the vibrational frequencies of the central region to be inside the bandwidth (given by ±2 √ k λ ) we choose the force constant in the central region k C small enough. As the force constants in the reservoirs are equal, the cut-off frequency also becomes λ-independent ω c,L = ω c,R ≡ ω c . We also set the Boltzmann constant k B = 1. We consider the temperature scale between the subsystems as a difference ∆T = T L − T R and set the temperature for the central region as T C = (T L + T R )/2. We fix the temperature in the right reservoir T R = 1 and relate the remaining ones to this. After the systems are coupled the temperature in the central region loses its meaning during transient due to nonequilibrium conditions. Even if the system reaches a steady-state, it is not a thermal equilibrium but a nonequilibrium steady-state. Therefore, instead of "thermalization" we say that the system relaxes, and we use a measure for defining the relaxation time τ (time from t = 0 to reach the steady state) as κ = 10% 75 . Evaluation of Eq. (44) also assumes the relative error to stay under the given tolerance at times t > τ . As the structure of the central region and its couplings to the reservoirs are arbitrary, it is feasible to study also more complex systems, such as graphene ribbons 76 , with the same methodology. Also, the nature of the heat transport (normal or anomalous 77,78 ) may be analyzed by different models with impurities or inner reservoirs. These type of simulations will, however, be postponed to future work since the focus of the present work is to benchmark the derived formula against exact results from the full propagation of the Kadanoff-Baym equations.

B. Validity of the approximations
Compared to the full numerical solution of Eq. (12), we have made a number of approximations when deriving the analytic solution in Eq. (35). We study how much error each approximation causes when compared to the full numerical solution. First, we make the wide-bandlike approximation for the self-energy in Eq. (23)  We benchmark the validity of the approximations stated above by studying the heat current through a dimer molecule as the central region; the force constant matrix in Eq. (43) is then 2 × 2. The parameters are chosen so that: (a) the vibrational frequencies of the dimer molecule are comparable with the cut-off frequency, i.e., the wide-band-like approximation is expected to fail; and (b) the vibrational frequencies of the dimer are in a narrow range compared to the reservoir bandwidth, i.e., the wide-band-like approximation should not neglect the detail of the spectrum. In case (a) we take k λC = 1/2 and k C = 1, and in case (b) we take k λC = 1/4 and k C = 1/3. In both cases we set the temperature profile by T L = 5, ∆T = 4.
In Fig. 3 we show the heat current, J Q 12 , between the dimer's atoms evaluated using Eq. (42) for the two abovementioned cases. In Fig. 3(a) we see more deviation from the full solution (red, thick solid line) compared to Fig. 3(b) as was expected from the parameter choice. By neglecting the frequency dependency of the phonon band in Eq. (19) (see Fig. 2) we see in both panels (green, long-dashed curve, WB) that the transient behavior is overestimated due to too crude of an approximation for the band edges. (The approximated embedding self-energy abruptly drops to zero.) This also affects the long-time behavior as the current saturates towards a steady-state value faster than the full solution because the coupling strength (dissipation) is overestimated in WB. Also, the steady-state value for the heat current is overestimated. When we also add the approximation for the time-convolution [D R ·Π < ] (blue, dash-dotted curve, WB-1) we do not, at least in this case, see any qualitative difference to WB. This means that approximating the retarded Green's function in Eqs. (B4) and (B5) as the embedded one only slightly modifies the initial transient and the steady-state values. When we consider the approximation only for the time-convolution [D < · Π A ] (magenta, short-dashed curve, WB-2) we see a relatively good match with the full solution. However, individual density matrix elements may still differ considerably between the approximated and full solutions, see Fig. 4 and Sec. III E. When deriving [D R · Π < ] in Eq. (B10) and [D < · Π A ] in Eq. (B14) we implicitly assumed the limit ω c,λ → ∞. Now, when the cut-off frequency is set to the phonon bandwidth, in Fig. 3(a) the approximations do not fully take into account the broader spectrum of the central region. In Fig. 3(b) the spectrum of the central region is more narrow, and all approximations give a reasonable agreement. Even if the general trend of the transient is qualitatively captured in both cases, quantitative differences can still be considerable, see the insets in Fig. 3. When combining all the approximations, we get the analytic result in Eq. (35) (cyan, thin solid line). It is worth pointing out that this result, without needing any numerical evaluation of the Green's function, can give a comparatively good description.
To elucidate the transient behavior further we consider also the time-dependence of the commutator of the displacement and momentum operators [û 1 ,p 1 ] at site 1; this is shown in Fig. 4. We see that this quantity in 'Full' and 'WB' equals the canonical commutation relation, i, but when the time-convolutions are approximated in 'WB-1' and 'WB-2' we see a deviation. For these elements of the density matrix (from which the commutator is computed) this deviation is more significant than for the heat current in Fig. 3. Further, this deviation can partly be tracked down to the approximation made for the retarded Green's function in Eq. (B5). From Eq. (4) and using withφ the composite field operators in Eq. (2), and the commutator is understood component-wise. (Each block in the above 2 × 2 matrix is an N C × N C matrix.) However, if we directly Fourier transform Eq. (B5) and take the equal-time limit, we end up with This means that the commutation relation, in case of the retarded Green's function, would only be satisfied in the limit of weak coupling, Γ 0 → 0. Fig. 4(b) shows that this is related to the commutation relation obtained from the lesser Green's function (density matrix elements), i.e., we observe a deviation of only few percent when the coupling is weak.

C. Dependency on the cut-off frequency
We further compare the obtained analytic solution to the full numerical solution of the equation of motion (12). We perform numerical integrations of Eq. (35) with varying cut-off frequency in terms of the bandwidth 2 √ k λ . In Fig. 5 we show the time-dependent energy in a single site connected to two reservoirs calculated from Eq. (40). We set T L = 5, ∆T = 4, and study in Fig. 5(a) strong coupling k λC = 1/2 with broad spectrum k C = 1, and in Fig. 5(b) weak coupling k λC = 1/4 with narrow spectrum k C = 1/3. Both transients (the full solution) show similar features as the energy in the site first grows and oscillates and then saturates; reaching the steady-state takes longer in the weak coupling case, as is to be expected. In fact, since we are dealing with one site only, this quantity E(t) can by the equipartition theorem be related to the time-dependent local temperature in the site 69,70,79,80 . At t = 0 the temperature is as prepared for the uncoupled system T C , and after the transient the system saturates to a value within the temperature "bias" window [T L , T R ].
The insets in Fig. 5 show the frequency dependency of the uu component of the steady-state integrand f L (ω)B L (ω) in logarithmic scale and the integration limits corresponding to the cut-off values. As can be seen from the insets, the value of the integral depends on the cut-off frequency, and more specifically, on how well the spectral peaks in the central region fit in this frequency window. Also, as discussed earlier, the integral would not converge due to the self-energy approximation but it is regulated when ω → −∞ by the cut-off frequency. We see that in case (a) the integrand is broader and the spectral peaks are further away from ω = 0 leading to a mismatch between the numerical solution with the full embedding self-energy and the derived time-dependent density matrix with the self-energy approximation. Increasing the cut-off frequency brings the time-dependent data closer to the full solution as less spectral weight is left outside of the integration region. It is still important to remember that if we do not limit the integration by the cut-off frequency the approximate solutions will grow (slowly) without a bound; in principle they do not approach the full solution. Due to the unsuitable parameters in case (a) the curves agree only close to the steady- state with higher cut-off values. In case (b) the match is better since the integrand is more narrowly peaked and the peaks are closer to ω = 0. The oscillation frequencies correspond to the the peak separation in the spectral function. It can be seen that the transient frequency in the narrow spectrum case (b) is indeed smaller than in the broad spectrum case (a).

D. Heat current in a four-site system
Next, we study the heat current in the middle of a 4-site chain given by the local heat current in Eq. (42) between sites 2 and 3, see Fig. 6. Also here, we set the temperature difference as T L = 5, ∆T = 4. In Fig. 6(a) we have strong coupling k λC = 1/2 with broad spectrum k C = 1, in Fig. 6(b) intermediate coupling and spectrum k λC = 1/4, k C = 2/3, and in Fig. 6(c) weak coupling k λC = 1/8 and narrow spectrum k C = 1/3. Like in the previous example, we solve numerically the equation of motion (12) (solid red lines) and compare this to numerical integration of Eq. (35) (dashed green lines) by setting the cut-off frequency to the phonon bandwidth, ω c,λ = 2 √ k λ . Due to the strongest coupling in Fig. 6(a) the steadystate is reached the fastest, but, as the spectrum is broad (reaching the edge of the band with tails outside the band) the transient is not captured very accurately. In Fig. 6(b) the spectral peaks are inside the band although the broadening still causes differences with the full numerical result due to the cut-off frequency. Even though in Fig. 6(b) most transient features are already captured, in Fig. 6(c) the correspondence is even better since the spectrum is well inside the band with very narrow peaks. Also, in a longer time scale [see the insets in Figs. 6(b) and 6(c)], the transient oscillations can reach several hundred units of time before saturation. The approximate steady-state values overestimate the full solution by roughly 10%. More generally, we observe in steadystate the heat current to be positive, meaning a steady flow from the hot to the cold reservoir. However, if the temperature differences were smaller, the current could also momentarily (during the transient) flow from the cold to the hot reservoir (not shown). This finding is not related to the approximations (the effect was also seen in the full solution) but to the partitioned model with contacts being suddenly switched on 46,76 ; the effect would fade away in the case of adiabatic switching.

E. Transient time shifts in individual density matrix elements
Even though the comparisons between the full numerical solution and the derived analytical result in Figs. 3 and 6 matched reasonably well, we may still wonder about the individual density matrix elements. This is due to the heat current in Eq. (42) being an "averaged-out" quantity as it is evaluated as a sum of two contributions. In Fig. 7 we consider the same numerical simulation as in Fig. 6(c) but now we plot separately the different terms in the definition of the heat current (and also the net current for reference). We notice that, even if the net current is described very accurately, the "in-coming" and "outgoing" component of the current differ considerably (time shift) at short times between the full numerical solution (solid lines) and the analytic result in Eq. (35) (markers in corresponding color). This deviation relates to the approximation made for the time-convolution between the lesser Green's function and the advanced embedding selfenergy; the short-time failure is considered more in detail in App. B. If the cut-off frequency is not large enough compared to the relevant frequencies in the central region, the time shift becomes relatively more severe, see also Fig. 4.
This effect also relates to the fulfilment of the commutation relation investigated in Fig. 4. The approximations for the time-convolutions in deriving Eq. (35) result in up and pu blocks of the density matrix possibly having complex elements when j = k although they should be real-valued. This means that the heat current components, when considered separately as in Fig. 7, gain a nonzero imaginary part. This imaginary part, however, cancels when evaluating the net heat current as a differ-ence in Eq. (42).
If an accurate description is desired also for the individual matrix elements, it is possible as a "rule of thumb" to take the time shift into account by transforming the time axis for the approximate solution by t → t(1 + ae −bt ), where a and b are positive real numbers depending on the system parameters (a = 2 √ Γ , b = ω c,λ √ Γ /4), to better represent the initial transient. We have tested (not shown) that during the initial transient the curves will be mostly in phase due to this shifting. Also, at larger times when the shifting goes to zero, the curves will stay in phase as argued in App. B.

F. Relaxation time scales
We may also use the time-dependent formalism to estimate relaxation time scales using Eq. (44). We consider similar atomic chains as in the previous sections but now we vary the length of the chain and the strength of the coupling between the chain and the reservoirs. We expect to see longer relaxation times for longer chains as it simply takes longer time for the wavefront to propagate over longer systems. Also, by increasing the coupling strength between the central region and the reservoirs the relaxation times should decrease due to the stronger dissipation. We choose the model parameters in physical units as k λ = 1.0 eV/(Å 2 u), k C = 0.625 eV/(Å 2 u) where u is the atomic mass unit. (Notice that these are the mass-normalized force constants, i.e., in SI units [k] = 1/s 2 .) In Fig. 8 we display the relaxation times (colormap) for chains of varying length (horizontal axis) and with varying coupling strength (vertical axis) having a 10% temperature difference between the reservoirs. In Fig. 8(a) the baseline temperature is 10 K (leading to T L = 11 K and T R = 9 K), and in Fig. 8(b) we have the baseline temperature 300 K (leading to T L = 330 K and T R = 270 K).
Indeed, we observe that, for both low and high temperature regimes, the relaxation is fastest in the shortest chain (N C = 2) when the coupling strength between the central region and the reservoirs is the highest (k Cλ = 0.9k C ). Expectedly, the slowest relaxation, which can take up to picosecond scale, occurs in the longest chain (N C = 16) when the coupling is the weakest (k Cλ = 0.2k C ). The overall dependency on the studied parameters is roughly similar in both low and high temperature regimes. However, we see in high temperature and small coupling that even the mid-range chains can take a comparatively long time for relaxation. On the other hand, increasing the coupling strength in the high-temperature regime leads to comparatively faster relaxation. This is also partly due to the larger absolute difference in temperature in the high-temperature case.

IV. CONCLUSION
We discussed a quantum transport model for noninteracting phonons in an arbitrary harmonic lattice setup. Using the nonequilibrium Green's function approach we derived an analytic result for the time-dependent oneparticle reduced density matrix. The derivation involved a wide-band-like approximation for the embedding selfenergy of the reservoirs; compared to the conventional wide-band approximation in electronic transport, now the real part of the embedding self-energy was set to be a nonzero constant around zero frequency whereas for the imaginary part we introduced a linear approximation. In addition, we introduced a cut-off frequency above which the embedding self-energy is simply set to zero. In the steady-state limit our analytical formula reduces to a known result 51,52 .
As an application of the derived formula we analyzed transient heat currents in atomic chains. Using the numerical examples we were able to benchmark our derived result against a full numerical solution to the equations of motion for the Green's function, and furthermore test the validity of the approximations put forward when deriving the result. The approximations were found to be reasonable and the benchmark results congruent when the cut-off frequency was chosen large enough compared to the relevant energy scales of the studied systems, and when the coupling strength between the central region and the reservoirs was small enough for the wide-bandlike approximation to hold.
Due to its computational simplicity, the introduced method holds promise for studying transient heat transport especially in large spatial scale, e.g., in graphene ribbon or carbon nanotube circuitries. It remains to be investigated if and how a partition-free approach, conventionally used in electronic transport 21 , could be incorporated also in the context of (time-dependent) phonon transport. A related issue, in the case of heat transport due to electrons, has recently been discussed in Refs. 81 and 82, the latter concentrating on the transient regime.

ACKNOWLEDGMENTS
We thank the Vilho, Yrjö and Kalle Väisälä Foundation, the Academy of Finland, MIUR FIRB Grant No. RBFR12SW0, and EC funding through the RISE Co-ExAN (GA644076) for financial support. We further wish to acknowledge CSC -IT Center for Science, Finland, for computational resources.
Appendix A: Self-energy calculations From Eq. (17) we may deduce the greater and lesser components for the uncoupled Green's function to be and further, the retarded component to read as Since the retarded Green's function is a function of t − t , we find by Fourier transforming where, for the second equality, we used the idempotency of α. The parameter η is a positive infinitesimal accounting for the correct causal structure. We can then insert Eq. (A4) into Eq. (14) and derive an expression for the retarded self-energy for the central region embedded in the environment. We may evaluate the retarded Green's function by performing a matrix inversion for a block matrix resulting to If we know the eigen decomposition of K λλ as K λλ X = Xω 2 λ , we can then also diagonalize the full reservoir Hamiltonian with X ≡ diag(X, X) as Further, we may write the retarded embedding selfenergy in terms of the eigenmodes when inserting Eq. (A5) into Eq. (14) where we explicitly labelled the basis elements of the central region (j C , k C ) and the reservoirs (q λ ). Also, Ω j C q λ ≡ (Ω Cλ X ) j C q λ . Since the coupling Hamiltonians only have nonzero elements in the uu block, this leads to the embedding self-energy having nonzero contribution also in the uu block only, and we have where K j C q λ ≡ (K Cλ X) j C q λ . The advanced embedding self-energy Π A λ is simply given by complex conjugating Eq. (A8). Then, we may evaluate the level-width func-tion Γ λ defined as where we used the lorentzian representation for the delta where Λ λ and Γ λ are real functions related by the Hilbert transform From Eqs. (A9) and (A11) we notice that Λ λ is an even function and Γ λ is an odd function.
We keep D R so far unspecified and we calculate the time convolution where we inserted the Fourier transform of Π < , changed the integration variable as t −t = t and inserted a step function for extending the time interval to minus infinity. For the step function we may use the expression and evaluate further where we used the Fourier transform of D R (t) and changed the integration variable to ω − ω =ω. The exponential involving both frequencies can be split up, and we may also insert the approximation for the embedding self-energy The cut-off frequency ω c,λ is now explicitly in this expression without specifying D R . As we argued earlier, this expression is only valid in the cut-off regime and we simply use the embedded retarded Green's function for all frequencies ω (in the inner integral) where we inserted the approximation for the embedding self-energy and defined the effective Hamiltonian as in Eq. (25). Now the retarded Green's function is specified, and D R and Π R satisfy the Dyson equation in the limits of |ω| < ω c,λ . In Eq. (B4), if the eigenvalues of Ω eff lie in the lower-half plane [see the discussion after Eq. (25)], then also the analytical structure for the Green's function is correct: D R (ω ) has poles in the lower-half plane, and also the denominator goes to zero when ω = ω − iη (in LHP). The key for evaluating the time-convolution was the approximation for D R in Eq. (B4). As this is done within the cut-off frequency window θ(ω c,λ −|ω|), we may analyze how large a difference does this approximation make compared to a full numerical solution of the equations of motion (12). As we approximate the retarded Green's function for all frequencies ω as in Eq. (B5), this means the limit ω c,λ → ∞. On the other hand, when we specify the cut-off frequency directly to the self-energy approximation in Eq. (23), this would amount to As we discuss only the region |ω| < ω c,λ when evaluating the time-convolution in Eq. (B4), we may compare how much the approximated Green's function deviates from that in Eq. (B6) (outside the cut-off window) where ω eff are the real parts of the eigenvalues of the effective Hamiltonian Ω eff . This limit means: (1) we choose the cut-off frequency high enough so that the physical frequencies of the central region fall well inside this window; (2) if the frequency in the retarded Green's function still was higher than ω c,λ , we would have, in the limit Γ 0 → 0 (weak coupling), that the relative difference in the retarded Green's functions approaches unity. Therefore, in order to make the approximation better is then two-fold: We can tune the force constant in the reservoirs so that the cut-off window is considerably broader than the energy scale of the central region; and/or we can also couple the central region more weakly to the reservoirs, thus, decreasing the value for Γ 0 . By using the eigenbasis of the effective Hamiltonian Ω eff , and if the eigenvalues lie in the lower-half plane [again, also see the discussion after Eq. (25)], we can evaluate the integral in Eq. (B4) over ω by closing the contour in the lower-half plane (t is a positive number): Then we can take the limit η → 0 and write Inserting this into Eq. (B4) we get It is worth noticing that this result could also be derived by In the limit ω c,λ → ∞ the sinc functions become delta functions, i.e., lim →0 1 πx sin x = δ(x) and we obtain This is naturally the same result as if we put the limits of the integration in the derivation of Eq. (B11) to ±∞. Based on the above expression for the advanced embedding self-energy we aim to calculate the time convolution in the equation of motion and the corresponding hermitian-conjugated one. The higher cut-off frequency ω c,λ we choose for the advanced embedding self-energy, the faster will the oscillations (in time) be in Π A (t, t). On the other hand, the fastest oscillations for D < (t,t) correspond to the (physical) vibrational frequencies of the central region. If we choose the cut-off frequency ω c,λ to be considerably higher than the typical energy scales in the central region, then in timedomain, Π A (t, t) appears almost as Eq. (B12) compared to D < (t,t). This allows us to calculate where we integrated by parts and noticed that the boundary term vanishes. By conjugating Eq. (B14) we also find Π R · D < . The result in Eq. (B14), however, implicitly assumes the limit ω c,λ → ∞ as we motivated its derivation by comparison of energy scales in Π A and D < exactly this way: ω c,λ ω for frequencies ω in the central region.
Let us try to justify this proposition by evaluating the time convolution of Π A and D < also by using the explicit expression in Eq. (B11) and performing an asymptotic expansion in ω c,λ in where the upper limit of the integration follows from the advanced nature of Π A (t, t) ∼ θ(t −t). Using Leibniz' rule we may write the second term of Eq. (B15) as This tells us that we only need to consider "the first line"like terms in Eq. (B15) and then we get "the second line"-term in Eq. (B15) by inserting ∂ t D < (t,t)Γ /2 and D < (t,t)Γ /2 as f (t,t) and using Eq. (B16). We then consider the behavior of the function F in Eq. (B17) when ω c,λ is large. We can get rid of the sinc function structure by taking the derivative with respect to ω c,λ . (We assume we are allowed to differentiate under the integral sign since the functions are well-behaving in our case.) We obtain for which we may perform a sequential integration by parts since we know all the anti-derivatives of a cosine function dF (t) dω c,λ = π −1 f (t,t) 1 ω c,λ sin[ω c,λ (t − t)] where the remainder will only be higher order in 1/ω c,λ . This is justified by assuming that the derivatives |∂tf (k) (t,t)| remain small compared to ω k+1 c,λ which is true as the oscillations in the function f relate to the frequencies in the central region [see Eq. (B15)] which are by construction assumed small compared to the cutoff frequency. Since we are interested in the large ω c,λ limit, we may simply take only the leading term, and integrating once over ω c,λ we obtain c,λ t 0 (B20) where Si(x) is the sine integral. Now we may consider Eq. (B15) by inserting the corresponding parts as the function f in Eq. (B20) with the help of Eq. (B16). Further, by using an asymptotic expansion for the sine integral 83 we may conclude that the terms neglected by approximating Eq. (B15) by Eq. (B14) are of the order O[(ω c,λ t) −1 ]. Therefore, for long times the approximation is reasonable but even if we choose a large cut-off frequency ω c,λ , for short times t 1/ω c,λ the approximation fails.
Appendix C: Integrating Eq. (34) It is convenient to start by making a transformation Before we start integrating over t, recall the matrix structures In Eq. (C2) we have terms such as where, conveniently, which can be checked by simply evaluating the matrix products and inverses. Then, looking at the first row of Eq. (C2), we can perform the integration over t by using the formula for arbitrary matrices A and B. The second row of Eq. (C2) is simple to integrate over t since there is only one exponential depending on time in each term. After integration we arrive at On the left-hand side we have the initial-state Green's function (at t = 0). We are working in the partitioned scheme, i.e., the systems of different temperatures are coupled at t = 0, so the initial condition should be equal to the uncoupled Green's function as in Eq. (A2) but for the indices in the central region: Here f C gives the thermal distribution according to which the central system is prepared before it is connected to the reservoirs, and Ωα is correspondingly the uncoupled Hamiltonian. Now, in Eq. (C9) we may transform back from D < (t, t) to D < (t, t) by multiplying with the exponentials from left and right and obtain our final result iD < (t, t) which can be written as Eq. (35) in main text by introducing the spectral function B λ (ω) in Eq. (37).