Controllable Thermal Conductivity in Twisted Homogeneous Interfaces of Graphene and Hexagonal Boron Nitride

Thermal conductivity of homogeneous twisted stacks of graphite is found to strongly depend on the misfit angle. The underlying mechanism relies on the angle dependence of phonon–phonon couplings across the twisted interface. Excellent agreement between the calculated thermal conductivity of narrow graphitic stacks and corresponding experimental results indicates the validity of the predictions. This is attributed to the accuracy of interlayer interaction descriptions obtained by the dedicated registry-dependent interlayer potential used. Similar results for h-BN stacks indicate overall higher conductivity and reduced misfit angle variation. This opens the way for the design of tunable heterogeneous junctions with controllable heat-transport properties ranging from substrate-isolation to efficient heat evacuation.


Model system
The initial intralayer carbon-carbon covalent bond length was taken as 1.420 Å (the equilibrium value obtained using the REBO intralayer potential for single layer graphene) and the initial intralayer boron-nitrogen covalent bond length was taken as 1.442 Å (the equilibrium value obtained using the Tersoff intralayer potential for single layer h-BN). The initial interlayer distance across the layered stack was set equal to 3.4 Å and 3.3 Å for graphite and bulk h-BN, respectively. Periodic boundary conditions were applied in all directions. It should be noted that the lattice structure is rigorously periodic only at some specific twist angles, the values of which are listed in Table S1 in section 3 below. While the cross-sectional area for each misfit angle, , is different, all systems considered have a contact area exceeding 12 nm 2 , which was shown to provide converged results with respect to unit-cell dimensions (see Section 2.4). 1 The intralayer interactions within each graphene and h-BN layer were modeled via the second generation REBO potential 2 and Tersoff potential, 3 respectively. The interlayer interactions between the layers of graphite and bulk h-BN were described via our dedicated interlayer potential (ILP), 4 which is implemented in the LAMMPS 5 suite of codes. 6 Figure S1. Schematic representation of the simulation setup (a) and steady-state temperature profile (b), respectively. In panel (a), two identical AB-stacked graphite slabs (gray and orange, respectively) are twisted with respect to each other to create a stacking fault of misfit angle θ. A thermal bias is induced by applying Langevin thermostats to the two layers marked by dashed red (Thot) and green (Tcold) rectangles. The arrows indicate the direction of the vertical heat flux. Since periodic boundary conditions are applied also in the vertical direction, two twisted interfaces are, shown across which heat flows in opposite directions. The steady-state temperature profiles are illustrated in panel (b), where N is the total number of layers in the model system and RAB, dAB and Rθ, and dθ mark the interfacial Kapitza resistance 7,8 and interlayer distance for contacting graphene layers with ABstacking and misfit angle θ, respectively. The red lines in panel (b) mark the temperature variation across the twisted interface, where the vertical axis corresponds to the position of the various layers along the stack and the horizontal axis marks the temperature of the various layers.

Simulation Protocol
All MD simulations were performed with the LAMMPS simulation package. 5 The velocity-Verlet algorithm with a time-step of 0.5 fs was used to propagate the equations of motion. A Nosé-Hoover thermostat with a time constant of 0.25 ps was used for constant temperature simulations. To maintain a specified hydrostatic pressure, the three translational vectors of the simulation cell were adjusted independently by a Nosé-Hoover barostat with a time constant of 1.0 ps. 9 To relax the box, we first equilibrated the systems in the NPT ensemble at a temperature of T = 300 K and zero pressure for 250 ps (see Figure S2). After equilibration, Langevin thermostats with damping coefficients 1.0 ps -1 were applied to the bottom and middle layer of the graphene stack (see Figure S1) with target temperatures Thot = 375 K (hot reservoir) and Tcold = 225 K (cold reservoir), respectively. Then the system was allowed to reach steady-state over a subsequent simulation period of 750 ps (see Figure   S2), during which the dynamics of all non-thermostated layers followed the NVE ensemble. For the larger model systems, the length of the NPT and Langevin stages was doubled (for the 32 and 48 layers systems) or tripled (for the 104 layers graphitic system) to ensure convergence of the obtained steady-state. Once steady-state was obtained, the last 500 ps were used to calculate the thermal conductivity of the twisted graphite and bulk h-BN. The statistical errors were estimated using ten different data sets, each calculated over a time interval of 50 ps.

Calculation of the interfacial thermal resistance
According to Fourier's law, the cross-plane thermal conductivity ( CP ) of a twisted graphitic interface of misfit angle θ can be calculated as: where ̇ is the heat flux, A is the cross-section area and Δ /Δ is the temperature gradient along the direction of heat flux (perpendicular to the basal plane in our case). Figure S1(b) shows a schematic temperature profile along the z-direction, where the vertical axis corresponds to the position of the various layers along the stack and the horizontal axis marks the temperature of the various layers. The actual temperature profiles extracted from the NEMD simulations for twisted graphite with different number of layers can be found in Figure S3 . For Bernal-stacked graphite (i.e., θ = 0°, red circles in Figure S3), only the linear region of the temperature profile was used to calculate CP and the points corresponding to the layers where the thermostats were applied were omitted (marked with green triangle in Figure S3). The CP of the system was calculated using Eq. (S1.1) by averaging over the two linear regions of the temperature profiles. For the twisted case (θ ≠ 0°), we found a sudden temperature decrease Δ at the position of the twisted interface (see black squares in Figure S3). CP , in this case, was calculated using the temperature gradient calculated for the same layer range as that for θ = 0°. To characterize the thermal properties of the twisted interface, the concept of interfacial thermal resistance (ITR), i.e., Kapitza resistance, 7,8 was introduced. Using the definition of the Kapitza resistance, 7 = Δ /̇, and noticing that Δ tot = ( /2 − 3)Δ AB + Δ and Δ = ( /2 − 3) AB + , where is the number of layers and AB , Δ AB and , Δ are the interlayer distance and temperature difference for adjacent AB-stacked and twisted graphene layers, respectively, Eq. (S1.1) can be rewritten as follows [see Figure S1(b)]: where AB and are the ITRs of adjacent AB-stacked and twisted graphene layers, respectively.
The first term on the left-hand side is just the sum of resistances of the various interfaces within the two optimally stacked slabs. The number of layers is divided by two to account only for one part of the system that is located between the two thermostats (see Figure S1) and we remove the three interfaces corresponding to the thermostats and the twisted interface. The second term on the lefthand side is the resistance of the twisted interface itself. On the right-hand side we have the overall S5 resistance (not to be confused with the resistivity) expressed as the inverse of the overall junction conductivity multiplied by the thickness. For the aligned contact ( = 0°, = ), the ITR can be simply calculated as AB = AB / AB . Once AB ( ) is known, ( ) can be calculated from Eq. (S1.2) using the value of CP ( ). We note that the sharp temperature drop at the twisted interface indicates that should be much larger than AB . Figure S3. Temperature profiles for graphitic stacks consisting of (a) 8 and (b) 16 layers. The red circles and black squares represent the temperature profiles for the aligned ( = 0°) and twisted ( = 30.16°) junctions, respectively. Green triangles represent data points that were omitted in the CP calculation.

Effect of NPT simulation time on the calculated thermal conductivity
According to our simulation protocol, the simulation box of the systems were relaxed using the NPT esemble for at least 250 ps (see Section 1.2) prior to the non-equilibrium simulation stage. To check the effect of the NPT equilibration time on the evaluated thermal conductivity, we performed additional convergence tests for the 8-layer graphite stack with θ = 0° by increasing the NPT simulation time to 1 ns. The results are summarized in Figure S4a,d. We find that the box size is already fully relaxed at 250 ps. More specifically, the averaged lattice constants for NPT simulation times of 250 ps and 1 ns were 2.4604 ± 0.0004 Å and 2.4605 ± 0.0003 Å, respectively. The relative difference between them is smaller than 0.01%. The corresponding calculated thermal conductivities were 0.196 ± 0.024 Wm -1 K -1 and 0.184 ± 0.023 Wm -1 K -1 , respectively (see Figure S5a). Furthermore, the residual in-plane stresses obtained at a simulation time of 250 ps in both x (-0.025 ± 0.138 GPa) and y (-0.012 ± 0.128 GPa) directions also indicate that the system is satisfactorily relaxed. Similar behavior was found for θ = 30.16° (see Figure S5b). Therefore, we conclude that the NPT simulation time used in our protocol (250 ps) is sufficiently long to obtain converged results.

Effect of residual in-plane stress on the calculated thermal conductivity
Our original simulation protocol starts from an initial configuration, where all graphene layeres have the equilibrium carbon-carbon covalent bond-length obtained from the REBO intralayer potential (1.420 Å ) and an interlayer distance of 3.4 Å . First, the system is equilibrated using an NPT ensemble simulation at 300 K and zero pressure with the REBO potential augmented by the ILP. This is followed by non-equilibrium NVE simulations with Langevin thermostats coupling the system to implicit heat baths, during which the simulation box dimensions are kept fixed and the thermal conductivity is evaluated. Since a previous study argued that the cross-plane thermal conductivity of graphite may be sensitive to in-plane stress, 10 we validated that any residual stress due to the difference in equilibration and non-equillibrium simulation protocols has minor effect on the calculated thermal conductivity. To this end, we performed additional simulations for the 8-layer graphitic stacks with θ = 0°, replacing the NPT ensemble equilibration step by an NVT equilibration step (see Figure S4c,f). Consequently, the unit-cell was not relaxed during the equilibration step and the residual in-plane stress was ~0.11 ± 0.06 GPa in both the x and y directions. The average lattice parameter obtained using the NPT and NVT equilibration protocols were 2.4605 ± 0.0003 Å, and 2.4602 Å, respectively, with a relative difference of ~0.01%. The corresponding calculated thermal S8 conductivities were 0.184 ± 0.023 Wm -1 K -1 and 0.186 ± 0.023 Wm -1 K -1 , respectively (see Figure   S5a). The small variations (within the error bar) in lattice parameters and calculated thermal conductivity values indicates that the residual stress due to the different protocols used during the equillibration and non-equillibrium dynamics simulation stages has minor effect on our predictions.

Effect of NPT equillibration temperature
Our original simulation protocol starts from an initial configuration, where all graphene layeres have the equilibrium carbon-carbon covalent bond-length obtained from the REBO intralayer potential (1.420 Å ) and an interlayer distance of 3.4 Å . First, the system is equilibrated using an NPT ensemble simulation at 300 K and zero pressure with the REBO potential augmented by the ILP. This is followed by non-equilibrium NVE simulations with Langevin thermostats coupling the system to implicit heat baths, during which the thermal conductivity is evaluated. To evaluate the effect of the equilibration step temperature, we repeated the calculations while setting the equilibration step temperature to 200 K and keeping the average temperature during the non-equilibrium calculation at 300 K. The results are presented in Figure S4b,e. We found that the time-averaged lattice parameter at steady-state during the thermal conductivity simulation stage was 2.4601 Å and the residual inplane stress was 0.21 ± 0.06 GPa (0.43 ± 0.06 GPa) along x (y) direction. The small differences in lattice parameters between the original (2.4605 ± 0.0003 Å) and control simulations (~0.016%) have a minor effect on the calculated thermal conductivity (see Figure S5a).

Effect of supercell cross-section area on the calculated thermal conductivity
To evaluate the convergence of our results with respect to the simulated supercell dimensions, we performed additional calculations for the 8-layer graphite stacks with θ = 0° and θ = 30.16° by increasing the supercell sizes (see Figure S5b). Supercells of cross-section areas of 26.83 nm 2 , 60.40 nm 2 , and 120.77 nm 2 were used for the aligned (θ = 0°) interface and 18.98 nm 2 and 37.94 nm 2 sized models were used for the twisted interface (θ = 30.16°). For the θ = 0° interface, following a 1 ns equilibration step under NPT conditions at 300 K and zero pressure, the corresponding calculated thermal conductivities were 0.183 ± 0.032 Wm -1 K -1 , 0.184 ± 0.023 Wm -1 K -1 , and 0.184 ± 0.010 Wm -1 K -1 , respectively. Similarly, for the θ = 30.16° twisted interface, we equilibrated the smaller and larger systems under NPT conditions at 300 K and zero pressure for 0.25 ns and 1 ns, respectively.
The corresponding thermal conductivities were 0.035 ± 0.002 Wm -1 K -1 and 0.034 ± 0.003 Wm -1 K -1 , respectively. These results clearly indicate that our calculated thermal conductivities are well converged with respect to the model interface cross-section area.

Comparision of the phonon spectrum and density of states calculated using the ILP and the Lennard-Jones potential
A proper description of the phonon dispersion is very important for studying the thermal transport properties. We already showed in the main text that the ILP gives more accurate values of the crossplane thermal conductivity for graphite than those predicted by the Lennard-Jones (LJ) potential. To further understand this finding, we repeated the phonon spectrum calculations for AB-stacked graphite obtained using the REBO interlayer potential and the ILP, 4 now using the AIREBO potential (with LJ parameters, ε = 2.84 meV and σ = 3.4 Å). The results are illustrated in Figure S6a. To identify the effect of interlayer potential on the phonon dispersions, we also calculated the phonon spectrum of monolayer graphene with the second generation REBO and AIREBO potential, respectively (see Figure S6b). Comparing Figure S6a and Figure S6b, we find that the interlayer potential mainly influences the phonon properties at the low energy regime. Therefore, the differences between the phonon dispersion curves calculated using the two potentials at the high energy regime mainly results from the intralayer potential terms. The comparison with the experimental phonon spectrum of graphite 11 shows that the phonon spectrum calculated by REBO+ILP agrees better with the experimental data in the low phonon energy regime, which is relevant for the cross-plane thermal conductivity calculation, than the spectrum obtained using AIREBO (see Figure S6c). The corresponding phonon density of states (DOS), as plotted in Figure S6d, also demonstrates the difference between the two force fields. This supports the reliability of the ILP for performing crossplane thermal conductivity calculations.

Transient MD protocol
For comparison purposes we also calculated the interfacial thermal conductivity (ITC) and resistance (ITR) of twisted bilayer graphene (tBLG) with the transient MD simulation approach 12-14 since the NEMD simulation protocol used in the main text becomes invalid in this case. In this protocol, the system was first equilibrated within the NPT ensemble at T = 200 K and zero pressure for 100 ps, which was followed by a 100 ps NVT ensemble equilibration stage and a 100 ps of NVE ensemble equilibration stage. After the system reached equilibrium, an ultrafast heat impulse was imposed on the top layer of the t-BLG for 50 fs to increase the temperature of the top layer from 200 K to 400 K, while that of bottom layer of tBLG remained unchanged. After the external heat source was removed, thermal energy flowed from the top layer to the bottom layer due to the temperature difference and the temperature of both layers approached 300 K when quasi-steady-state was reached. During the thermal relaxation time interval (500 ps), the temperature and energy of the system sections were recorded. The ITR could then be extracted using the following equation: [12][13][14] where t is the total energy of the top graphene layer, R is the ITR of the tBLG, A is the interfacial cross-section area, and bot and top are the instantaneous temperatures measured for the bottom and top layers of the tBLG, respectively. Note that in Eq. (S4.1) we assume a linear dependence of the heat flux on the temperature difference between the layers. The ITC of the tBLG is simply defined as ITC ≡ /ITR, where d is the average interlayer distance. Note that, for tBLG, the ITC is equivalent to the CP in the main text since there exists only one interface.
The ITC and ITR of tBLG as functions of misfit angle calculated with the transient MD simulation protocol are illustrated in Figure S7, demonstrating similar misfit-angle dependence as that for the NEMD protocol with Langevin thermostats exercised to obtain the results presented in the main text.
This further validates the reliability of the simulation protocol adopted in the main text, which is more suitable to treat thick slabs and allows to obtain a true steady-state.

Comparison between the transient MD and Langevin thermostat protocols
Comparing calculated for the multi-layered graphitic systems using the Langevin thermostat protocol ( Figure 2a of the main text) and the ITC (= CP in the bilayer case) calculated for bilayer graphene using the transient MD approach (Figure S7), a ~10 fold difference in magnitude is observed. We identify two main reasons for this discrepancy: (i) The dependence of CP on the number of layers (see Figure 2a in the main text); and (ii) the different simulation protocols implemented for the two systems due to their different thickness. Regarding the latter issue we would like to stress that in the bilayer system we cannot use the Langevin thermostat approach as coupling the interface layer directly to the implicit heat bath without any buffer layers may result in an unphysical behavior. To estimate the relative importance of the two contributions we extrapolated the CP data obtained using the Langevin thermostat approach (see Figure 4a in the main text) to the bilayer limit using a power-law fitting function. The extrapolated CP values obtained were 0.076 Wm -1 K -1 and 0.0061 Wm -1 K -1 for the 0 o and 30.16 o rotated interfaces, respectively (see Figure S8).
Comparing to the values obtained for the bilayer system using the transient MD approach (0.028 Wm -1 K -1 and 0.0033 Wm -1 K -1 , respectively), we find a difference of a factor of 2.7 and 1.8 for the 0 o and 30.16 o rotated interfaces, respectively. The remaining difference can be therefore attributed to the difference between the two simulation approaches. were obtained with the Langevin thermostat approach and transient MD approach, respectively.
To explain the stronger thickness dependence of CP observed for twisted interfaces, we note that there are two main factors (see Eq. (S1.2)) influencing the dependence of the cross-plane thermal conductivity on the thickness: (i) The ITR, , which measures the resistance of the twisted interface alone; and (ii) the resistance, AB , of the interface between each two optimally stacked layers within the slabs residing above and below the twisted interface. Thus, the equation for the thermal conductivity as a function of the number of layers, , may be written as follows (see Eq. (S1.2) for a full explanation): To investigate the origin of the higher slope of CP ( ) for the twisted interface, we varied the value of artificially. In Figure S9

Brillouin zone of supercell in tBLG
For tBLG, the lattice structure is rigorously periodic only at some specific misfit angles, θ, where the lattice vector 1 = 1 1 + 2 2 in the bottom layer equals the vector 2 = 1 1 + 2 2 in the top layer with certain integers m1, m2 and n1, n2. Here, 1 = (1,0) and 2 = �1/2, √3/2� are the primitive lattice vectors of the bottom layer and a is the lattice constant of monolayer graphene.
Thus, the exact superlattice period is then given by: 15 where θ is the angle between two lattice vectors 1 and 2 . In the simulations below, we always rotated the supercell such that its lattice vector is 1 = (1,0) and 2 = �1/2, √3/2�.

S16
In this case, the corresponding reciprocal lattice vector of the moiré superlattice satisfies the relation • = 2 , such that: Both the lattice vectors and the corresponding reciprocal lattice vectors of the superlattice of the tBLG are presented in Figure S10. Table S1 reports the parameters used to construct rhombus periodic supercells of different misfit angles that can be duplicated to construct a rectangular periodic supercell.

Special points for Brillouin zone integration
The calculation of the sum over wave vector q in eq 4 in the main text can be transformed to an integral using the relation Brillouin zone (BZ) and V is volume of the real-space unit-cell. The calculation of integral is usually inefficient since it requires calculating the value of the function over a large set of k points in the first BZ. To calculate such integrations more efficiently, simple k-point meshes can be replaced by a carefully selected set of special points in the BZ, , [16][17][18][19] over which the function is evaluated. The integral can then be estimated via: where is the weight of the i th data point, and = ∑ normalizes the weighting factors to unity.
The set of selected { } forms a grid in the irreducible Brillouin zone (IBZ), as is illustrated by the red points in Figure S10b. The coordinates of these points for a hexagonal lattice are presented in Eq. (S5.4).
This equation was used to calculate the transition rate presented in Figure 3 of the main text.

Temperature dependence of interfacial thermal conductivity
In the main text, the target temperatures of the Langevin thermostats for the bottom and middle layers of graphene and h-BN were set to 225 K and 375 K, respectively. After reaching the steady-state, the average temperature of the system was found to be ~300 K. To check the effect of average temperature on our results, we calculated the cross-plane thermal conductivity ( CP ) and the corresponding interfacial thermal resistance (ITR) at a different temperature gradient (325 K -475 K), resulting the average steady-state temperature of ~400 K. The protocol described in Section 1 above was used to perform these calculations, as well. Both ILP and Lennard-Jones (LJ) potential were tested for graphite whereas for the bulk h-BN simulations only the ILP was used. The results for graphite and bulk h-BN are illustrated in Figure S11 and Figure S12, respectively. For the ILP we find that the overall values of CP (ITR) decrease (increase) slightly with increasing average temperature. The LJ potential calculations, as well, exhibit very week dependence on average temperature within the range studied. Altogether, the thickness dependences of both quantities remain mostly insensitive to the average temperature.
where ( ) is its displacement from its equilibrium position.
The Lagrangian for this classical problem can be written as where the second term is the sum of interactions between all pairs of atoms.
Under the harmonic approximation, i.e., expanding the total potential energy around the equilibrium positions, The Lagrangian can be simplified as where , = 1,2,3 are the Cartesian coordinates of the displacement ( ) and Note that the first order term vanishes because we are expanding around the equilibrium positions.

Dynamical matrix
The equation of motion of the s th atom in the n th unit cell can be derived using the Euler-Lagrange equation as follows: If we displace all atoms equally, i.e. shifting ′ ′ ′ to ′ ′ ′ + , the total force on the s th atom in S20 the n th unit cell does not change. From the above equation we have We are looking for normal modes (because any general solution can be written as a linear combination of them); these are solutions where all atoms oscillate with the same frequency. Moreover, because of the lattice structure, we expect solutions to reflect this periodicity. So we guess solutions of the form where are real-space solutions that will be determined later and is the wave-vector in reciprocal space. Substituting Eq. (S7.7) into Eq. (S7.5), we can get is called dynamical matrix (dimension 3 × 3 ). Note that we have defined the relative cell distance It's clear that 2 ( ) and 2 (− ) obey the same equation, thus we have: The corresponding eigenvectors are orthonormal: The complex conjugate of Eq. (S7.8) gives

Second quantization
The general solution is a linear combination of all these normal modes, thus we have where the we define the normal coordinates as

Potential energy term
Similarity, the potential energy can be rewritten in terms of the normal coordinates as follows: i.e., where in above deriviation, the orthogonality of its eigenvectores and the symmetry of its eigenvalues are used. Thus the Lagrangian reads as: Now, we quantize H by asking the momenta and coordinates to be operators: here � and � obey the following commutation relations: Similar to the case of ordinary quantum harmonic oscillators, it is convenient to define ladder operators for each mode as follows: where � † and � are Bosonic creation and annihilation operators for phonons with momentum , branch index , and frequency , which obeys the Bosonic commutation relation: Substituting Eqs. (S7.27) into (S7.25) and using the properties Eq. (S7.26) and (S7.28), we have Therefore, we finally get the quantized representation of non-interacting phonons: The operator of atom displacements (Eq. 1.14) is expressed in terms of the phonon operators by: These equations will be used below.

Inter-phonon coupling within harmonic approximation 7.3.1 Hamiltonian with inter-phonon coupling
In this section, we consider systems that consist of two (or more) covalently bonded units that are weakly coupled between them. Unlike previous studies that considered phonon-phonon couplings resulting from anharmonicity effects, 20 here all phonon mode considered are Harmonic and the couplings arise from the division of the entire system into subunits. The Hamiltonian of the whole system can be written as: To derive the expressions of 1 , 2 and 12 in Eq. (S7.31), we consider the Hamiltonian of the whole system written as the function of the atomic displacements: Let's assume that subsystem I and II contain atoms with indeices ranging from = 1,2, ⋯ , /2 and = 2 + 1, 2 + 2, ⋯ , , respectively. Then the kinetic energy term in Eq. (S7.32) can be rewritten as: While the potential energy term is: Or equivalently, (S7.34)

Second quantization
We may now quantize this Hamiltonian and write it in the basis of the eigenstates of the coupled subunits. In second quantization, the atomic displacement operators of the two subunits are given in the following form: and the corresponding momenta operators as: Eq. (S7.35) reads: The second quantized kinetic energy operator is then written as: Correspondingly, the various potential energy terms are obtained by substituting Eq. (S7.38) in Eq. (S7.34): i.e., Going from the second line to the third line we used the fact that the sum over ′ runs between ±∞ and the summand depends only on the difference between and ′ , hence the sum is independent of the value of the index . Therefore, we can replace the sum over ′ by a sum over ≡ − ′ amd define ≡ − ′ . Following the same procedure, we can get the corresponding expressions for the second diagonal term: Similarly, for the off-diagonal terms we get: Where we have defined: It's easy to show that ′ ( ) has the following property: Following the same procedure, we have: we finally get the expressions of � 1 , � 2 and � 12 in Eq. (S7.31) as follows: Here h.c. means the Hermitian conjugate. Since the indices of � , � � ′ and � , � � ′ belong to the two system sections, we can define an abbreviated notation � and � using the index to identify which subsystem they belong to. In this case, the Hamiltonian of the coupled system can be S28 simplified as follows:

Dynamics of the ladder operators
In order to describe the dynamics of the ladder operators appearing in the Hamiltonian of Eq. (S7.47) we express them in the Heisenberg picture as follows: where we define the imaginary time ≡ and introduce the notation ≡ ( , ) and � ≡ (− , ).
The corresponding equation of motion for the ladder operators is given by: For the uncoupled system ( � = 0 � ) this gives: The commutator on the right-hand-side of Eq. (S7.51) can be evaluated as follows: i.e., Therefore, we have: The solution of Eq. (S7.53) is In summary, we have

Green's function
To describe thermal transport properties between different system sections, we use the formalism of thermal (or imaginary time) phonon Green's function. 21 To this end, we define the thermal Green's function for phonons as: where � is the time ordering operator and � is the statistical operator for the grand canonical ensemble (note that the chemical potential for phonons is zero): where the partition function is given by ≡ Tr� − � � , and To show this we shall first assume that > ′ such that Where we have used the commutativity of � and ± � ′ ℏ and the invariance of the trace operation towards cyclic permutations. Similarly, for < ′ we have: For simplicity, we will introduce the notation ′ ( , 0) ≡ ′ ( ). We further note that when using To show this, we shall again assume first that − ℏ < < 0 to write, that is Similarly, for ℏ > > 0 we have: Therefore, ′ ( ) can be expanded as a Fourier series in the domain [0, ℏ] as follows: (S7.59) where = 2 ℏ and the associated Fourier coefficient is given by Having proven the translational time invariance and the periodicity of the Green's functios we can now calculate it for the uncoupled system: such that: . Hence, we have: where � � is the Bose-Einstein distribution for phonons. Substituting Eq. (S7.66) in Eq. (S7.61) we obtain: Note that only those Green's functions of the form

Interaction picture
Next, we can proceed with calculating the Green's function of the coupled system. To this end, we define the coupling Hamiltonian operator term in the interaction picture as: The time evolution of � ( ) is then given by: Note that the Green's function defined above was given in the Heisenberg picture. To proceed, we need to transform it to the interaction picture: is the operator in interaction picture and the operator � is defined by: Note that while � is not unitary, it satisfies the following group property: and the boundary condition � ( 1 , 1 ) = 1 � . In addition, the derivative of � is simply: i.e., The exact thermal Green's function now may be rewritten in the interaction picture as: Where we used the fact that we are free to change the order of the operators within the time ordering operation (see pages 241-242 of Ref. 21 ).

Wick's theorem
Thus the Green's function can be expanded as 21 where 〈⋯ 〉 0 represents the ensemble average with respect to the non-interacting basis For consiceness we introduce the following notation: To simplify the calculation in Eq. (S7.77), we adopt Wick's theorem (see pages 237-241 of Ref. 21 ),which can be expressed as follows: Therefore, the second term in the numerator of Eq. (S7.78) can be calculated as: The first term in above equation is called disconnected part since the pairing is performed separately on � ( 1 ) and � ( ) � ′ † (0). All other terms have pairs that mix creation and annihilation operators of the Hamiltonian with � ( ) or � ′ † (0) and are said to have connected party. For simplicity we include all these terms in the notation Using this method, the third term in the numerator of Eq. (S7.78) can be calculated as: Here, the last term is denoted as Using Eq. (S7.82), the second and third terms on the right hand above equation can be simplified as follows: Where, we performed the following integration variables interchange 1 ⟷ 2 . Thus we have Similarity, we can calculate the fourth term of Eq. (S7.78) as: Higher order terms can be treated in the same manner (see pages 95-96 of Ref. 21 ). Then when substituting 2 , 3 , 4 and all higher order terms into Eq. (S7.78), we can simplify the numerator as follows: Noting that the expression (1 + 1 + 2 + ⋯ ) is exactly canceled with the denominator, the Green's function can be simplified as: where, (1) ( ) is calculated as: According to Wick's theorem 21 , the terms that contain − ′ and − † ′ † are equal to zero, thus the first term of above equation are calculated as:

S38
Simplification of above equation gives: Similarity, the second term of ′ (1) ( ) is calculated as: To derive the Green's function in frequency domain, we use the expression for Fourier series, then we have:

Second order approximation
The second order approximation of the Green's function in Eq. (S7.85), ( ), is calculated as:

Fermi's golden Rule
To get the Fermi's golden Rule, we need to calculate the retarded Green's function, which can be calculated from