$F^1_{\omega}(\widehat{x}_1)$ and $F^2(\widehat{x}_1)$ are also expanded in term of cosines. Following the steps in the above reference, we can express the temperature as $$\label{eq:FilmTemp}
\Delta \widetilde{T}_{1}(\widehat{x}_1) = [\underline{\underline{A}}^{-1}(\underline{\underline{f_1}}\underline{B}+\underline{f_2})]^{T}\underline{\phi (x)}$$ where the matrix $\underline{\underline{A}}$ contains elements $A_{00}=1-\frac{k_{00}}{4}$, $A_{0n}=-\frac{1}{2}k_{0n}$, $A_{n0}=-\frac{k_{n0}}{4}$, $A_{nn}=1-\frac{1}{2}k_{nn}$ and $A_{nm}=-\frac{1}{2}k_{nm}$ ($m\neq n \neq 0$) and $\underline{B}$ is a $N_{\omega}$ column vector of the unknown coefficients $B_{\omega}$, where $N_{\omega}$ is the number of discretization in phonon frequency. $\underline{\underline{f_1}}$ is a $N \times N_{\omega}$ matrix, consisting of the Fourier coefficients of $F^1_{\omega_i}(\widehat{x}_1)$ evaluated at each phonon frequency $\omega_i$ and $\underline{f_2}$ is a N column vector, consisting of the Fourier coefficients of $F^2(\widehat{x}_1)$.

Then, $\widetilde{g}^+_{1\omega}(x_1,\mu)$ and $\widetilde{g}^-_{1\omega}(x_1,\mu)$ can be expressed in terms of the unknown coefficients $B_{\omega}$ by plugging Eq. (\[eq:FilmTemp\]) into Eqs. (\[eq:ODEsolution\_gPlus\]) and (\[eq:ODEsolution\_gMinus\]).

Substrate {#subsec:substrate}
---------

The substrate can be treated as a semi-infinite region subject to a surface heat flux. Therefore, the BTE for the substrate becomes $$\label{eq:SubstrateBTE}
i\eta\widetilde{g}_{2\omega}+ \mu v_{2\omega} \frac{\partial \widetilde{g}_{2\omega}}{\partial x_2} = -\frac{\widetilde{g}_{2\omega}}{\tau_{2\omega}}+\frac{C_{2\omega}}{4\pi \tau_{2\omega}}\Delta \widetilde{T}(x_2) + \frac{1}{2}P_{\omega}v_{2\omega}|\mu | \delta (x_2),$$ where the unknown coefficients $P_{\omega}$’s are determined through the interface conditions.

We then apply the Green’s function method given in Ref. 61. The unknown distribution function in spatial frequency domain is then written as $$\label{eq:SubstrateDistribution}
