# Ab initio phonon transport based on nonequilibrium green's function formalism: a practical approach

Van-truong Tran, Roberto D'agosta, Marc Bescond, Sebastian Volz

► To cite this version:

Van-truong Tran, Roberto D'agosta, Marc Bescond, Sebastian Volz. Ab initio phonon transport based on nonequilibrium green's function formalism: a practical approach. physica status solidi (b), 2024, 2024, pp.2400353. 10.1002/pssb.202400353. hal-04889742

HAL Id: hal-04889742

https://hal.science/hal-04889742v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ab initio phonon transport based on Non-Equilibrium Green's Function formalism: A practical approach

Van-Truong Tran<fnref n="1" />*, Roberto D'Agosta<fnref n="2" />, Marc Bescond<fnref n="3" />,<fnref n="4" />, and Sebastian Volz<fnref n="5" />

Email: *vantruongtran.nanophys@gmail.com

## Abstract

We present a highly practical first-principles approach for studying phonon transport across various materials, ranging from low-dimensional systems to three-dimensional bulk structures, as well as nanostructures. This method integrates the Non-Equilibrium Green's Function (NEGF) formalism with ab initio calculations of dynamical matrices generated by Quantum Espresso. A detailed technical approach for extracting atomistic force constants and forming dynamical matrices that adhere to the tri-diagonal technique within the NEGF framework is provided. Additionally, an efficient method for constructing dynamical matrices of conventional cells from those of primary cells is introduced, remarkably reducing the computational costs associated with Density Functional Theory calculations. The robustness and applicability of this methodology are validated through several case studies.

## 1. Introduction

Understanding phonon transport is crucial in materials science due to its significant impact on thermal management and thermoelectric applications [1–5]. Efficient thermal management is vital for the performance and reliability of electronic devices [5–7], while enhancing thermoelectric materials' efficiency can lead to better energy conversion technologies [4,8–10]. However, studying phonon transport accurately and efficiently remains a significant challenge due to the limitations of existing computational methods. Density Functional Theory (DFT) calculations, while accurate, are computationally expensive, making it challenging to study large systems or complex materials. This high computational cost limits the ability to explore a wide range of materials and conditions [11–13]. Semi-empirical models, such as Force

<footnote n="1">¹Université Paris-Saclay, Centre de Nanosciences et de Nanotechnologies, 91120, Palaiseau, France</footnote>

<footnote n="2">²Nano-bio Spectroscopy group, Departamento de Polímeros y Materiales Avanzados: Física, Química y Tecnología, Universidad del País Vasco (UPV/EHU), Avenida de Tolosa 72, E-20018 San Sebastian, Spain and Ikerbasque, Basque Foundation for Science, Plaza de Euskadi 5, E-48009 Bilbao, Spain</footnote>

<footnote n="3">³Aix Marseille Université, CNRS, Université de Toulon, IM2NP UMR 7334, 13397, Marseille, France</footnote>

<footnote n="4">⁴Institute of Industrial Science, University of Tokyo, 4-6-1 Komaba, Meguro-ku Tokyo 153-8505, Japan</footnote>

<footnote n="5">⁵LIMMS, CNRS-IIS UMI 2820, The University of Tokyo, Tokyo 153-8505, Japan</footnote>
Constant models for phonons [14–16] or tight-binding models for electrons [17–19], require determining a large number of empirical parameters [14,17,20], which can be difficult and often result in low accuracy and lack of transferability when applied to other systems. The approximations involved may fail to capture intricate electron or phonon interactions and material-specific properties such as edge-deformation effect [21], or the presence of vacancies [22].

Regarding computations of transport properties, one often relies on methods such as the Boltzmann Transport Equation (BTE), Molecular dynamics (MD), or the Non-Equilibrium Green's Function technique (NEGF).

The BTE method is a powerful tool for modeling both electron and phonon transport, particularly in non-equilibrium systems, offering insights into microscopic details and capturing complex scattering events [23–26]. However, the BTE suffers from computational complexity, particularly for systems with strong anharmonicity or complex geometries [27–30], and often relies on approximations that may not capture all relevant physical phenomena [27–30]. In addition, although quantum generalizations within the BTE method are discussed, for example, the Wigner function provides a quantum analog to the classical phase-space distribution function used in the BTE [31,32], this method has limited ability to incorporate quantum effects [28,33].

MD simulations can also offer a powerful approach to studying thermal transport at the atomic level, providing detailed insights into the interplay between microscopic dynamics and macroscopic properties [34,35]. Advantages of MD include its ability to capture complex interactions between atoms and molecules, allowing for the investigation of non-equilibrium phenomena [36], the study of systems with varying geometries and boundary conditions [37,38], the exploration of transport mechanisms in diverse materials [39,40], and the accurate modeling of time-dependent phenomena, capturing the evolution of the system over time [41,42]. On the other hand, MD simulations are computationally demanding, often requiring significant computational resources and time [41,43]. Furthermore, they are limited by the accuracy of the interatomic potentials used to model interactions [44,45], making it crucial to carefully choose appropriate potentials for the system under study. Additionally, MD simulations can be challenging for systems with long-range interactions or large size differences between particles [41,46,47]. Clearly, quantum mechanical effects, for example, electron-phonon scattering, cannot be evaluated.

The NEGF formalism is another powerful tool for investigating transport properties, in particular in nanomaterials research [33,48–51]. Unlike the BTE and MD methods, which are based on classical approaches, the NEGF method operates at the quantum level, allowing for the inclusion of quantum interferences and detailed scattering processes [33,50,52–54]. Computations using NEGF can be challenging for large systems due to the computational demand posed by large matrices [48,50]. Fortunately, efficient algorithms, such as Sancho-Rubio iterative technique [55–57] and recursive technique [50], were developed to make it possible to reduce the size of the matrices used in NEGF, enhancing computational efficiency and feasibility. Importantly, NEGF works for both periodic and open-boundary systems and can easily incorporate defects into the Hamiltonian, enabling effective descriptions of quantum effects and phonon interferences [58–60].

The main objective of this work is to develop a numerical methodology to evaluate the phonon transport coefficients that integrate the NEGF formalism with DFT calculations. Here, we demonstrate this integration using the Quantum Espresso suite of codes (DFT-QE) [61]. This combination leverages the accuracy of first-principles calculations with the robustness of NEGF for phonon transport analysis at the quantum level. We focus on extracting atomistic force
The Green's function method for phonons is developed based on an analogy between the eigenvalue equation used to solve energy states derived from the Schrödinger equation and the equation used to solve angular frequency of vibration derived from Newton's second law of motion:

$$ \omega^2 U = DU, \qquad (2) $$

where $U$ is the column matrix containing the amplitude vectors of vibration of all the lattice sites, $\omega$ is the angular frequency, and $D$ is the dynamical matrix [64].

The retarded Green's function for phonon can be written as

$$ G = [\omega^2 - D_D - \Sigma_L^s - \Sigma_R^s]^{-1}, \qquad (3) $$

here, $\omega^2$ takes a role similar to energy $E$ in the Green's function for electrons [64]. $D_D$ is the dynamical matrix of the device (active region), which must be a tri-diagonal matrix block structure as in Eq. (1); and

$$ \begin{aligned} \Sigma_L^s &= D_{DL} G_L^0 D_{LD}, \\ \Sigma_R^s &= D_{DR} G_R^0 D_{RD}, \end{aligned} \qquad (4) $$

are the surface retarded self-energies contributed from the left and right leads ("s" stands for "surface"). $G_{L(R)}^0$ is the surface Green's function of the isolated left (right) contact. The self-energies were computed based on the Sancho-Rubio iterative technique [57,65]. It is also worth mentioning that this study primarily focused on the harmonic regime, consequently, many-body self-energies are absent from Eq. (3).

The transmission coefficient is defined as [33,50]

$$ T_p = \mathrm{Tr}\{\Gamma_L G \Gamma_R G^\dagger\}, \qquad (5) $$

where $Tr$ stands for the trace of the matrix, $\Gamma_{L(R)} = i(\Sigma_{L(R)} - \Sigma_{L(R)}^\dagger)$ is the injection rate at the left (right) contact. To compute efficiently phonon transmission in Eq. (5), a recursive technique [66] was employed to handle the size of the dynamical matrix of the device in the Green's function calculation. Phonon transmission was then computed as [64,67,68]:

$$ T_p = \mathrm{Tr}\{\Gamma_L^s [i(G_{11} - G_{11}^\dagger) - G_{11} \Gamma_L^s G_{11}^\dagger]\}, \qquad (6) $$

where $\Gamma_{L(R)}^s = i(\Sigma_{L(R)}^s - \Sigma_{L(R)}^\dagger)$ is the surface injection rate at the left (right) contact with the surface self-energies $\Sigma_{L(R)}^s$ are computed by Eq. (4). Here the subscript '11' in the Green functions refers to the same block form as in Eq. (1). It should also be noted that an equivalent equation to Eq. (6), involving the right surface injection rate and $G_{NN}$, can also be applied [68]. The heat flux can be computed as [69]

$$ I_{th} = \int_{-\infty}^{+\infty} \frac{d\omega}{2\pi} \mathrm{h}\omega T_p(\omega) [n_L(\omega, T_L) - n_R(\omega, T_R)], \qquad (7) $$

and the linear response phonon thermal conductance $K_p$ is obtained from [64]

$$ K_p = \lim_{\Delta T \to 0} \frac{I_{th}}{T_L - T_R} = \int_0^\infty \frac{d\omega}{2\pi} \mathrm{h}\omega T_p(\omega) \frac{\partial n}{\partial T}, \qquad (8) $$
with $n=1/(e^{h\omega/k_B T}-1)$ the Bose-Einstein's distribution function at temperature $T$, and $k_B$ the Boltzmann constant. By explicitly performing the derivative, we can recast Eq. (8) as

$$K_p = \frac{k_B}{2\pi} \int_0^\infty d\omega T_p(\omega) g^p(\omega, T), \qquad (9)$$

where $g^p(\omega, T) = \left(\frac{h\omega}{2k_B T}\right)^2 / \sinh^2\left(\frac{h\omega}{2k_B T}\right)$.

It is worth mentioning that the NEGF technique is applicable to not only 1D but also 2D and 3D structures, by classifying the transport direction and the other ones in which the structures are considered periodic. Specifically, the dynamical matrix $D_{n_1 n_2}$ of the $n_1$-th and $n_2$-th blocks along the transport direction can be added terms describing couplings between the $n_1$-th block and other blocks at the $n_2$-th slice perpendicular to the transport directions (along the periodic direction of the system). This description can be written as

$$D_{n_1 n_2} = d_{n_1 n_2} + \sum_m d_{n_1 m} e^{i q \mathbf{R}_{n_1 m}}, \qquad (10)$$

now $d_{n_1 n_2}$ describes only the couplings between the atoms of the $n_1$-th block and those of the $n_2$-th block. In this equation, **q** refers to the momentum transfer, while **R** is the vector defined from the $n_1$-th block to the $m$-th block. It is worth noting that to satisfy NEGF algorithms, only nearest-block interactions along the transport direction are required (resulting in a tri-diagonal block structure in the dynamical structure). Therefore, the condition of nearest-neighbor block interactions in other directions is not necessary.

If the system is 2D, for instance, in the case of graphene sheets, with transport along the x-axis and periodicity along the y-axis with lattice constant $L_y$, the Green's function includes $q_y$ as a consequence of Eq. (3) and Eq. (10).

$$G(E, q_y) = \left[ \omega^2 - D_D(q_y) - \sum_L^s (\omega, q_y) - \sum_R^s (\omega, q_y) \right]^{-1}, \qquad (11)$$

$D_D(q_y)$ is the dynamical matrix of the device region depending on the component $q_y$ of the wave vector.

The total transmission is now the sum of all transmissions over $q_y$ in the Brillouin zone:

$$T_p(\omega) = \int_{q_y \subset \text{BZ}} \frac{dq_y}{2\pi / L_y} T_p(\omega, q_y). \qquad (12)$$

Similarly, when the structure is 3D (bulk structure), the total transmission can be computed as,

$$T_p(\omega) = \int_{q_y, q_z \subset \text{BZ}} \frac{dq_y}{2\pi / L_y} \frac{dq_z}{2\pi / L_z} T_p(\omega, q_y, q_z), \qquad (13)$$

with $L_i$ the lattice constant along the $i$-th direction.

### 3. Method for phonon transport based on DFT-QE and NEGF

As seen from Eq. (12) and above, the evaluation of the transmission coefficient relies on the evaluation of the dynamical matrix $D_D$. Recent progress in computational materials science has
made it possible to obtain this matrix from state-of-the-art ab initio calculations, where the
interatomic forces are evaluated by solving the (Kohn-Sham) Schrödinger equation of the
many-body problem and the forces are obtained via the Hellmann-Feynman theorem. In the
following, we show how to use a nowadays mainstream DFT suite, Quantum Espresso, to
obtain the dynamical matrices and recast them in a form appropriate for the evaluation of
Green's functions. It is also worth mentioning that the proposed "ab initio" treatment for the
dynamical matrices is applicable to both the device region and the contacts.

**3.1. Extracting dynamical matrices from DFT-QE calculations**

This part is developed since DFT-QE generates data on force constants for each pair of atoms *i* and *j* atoms. However, to apply the NEGF technique, Hamiltonians must be cast in matrix form, representing the coupling between blocks or cells in real space. In the following, a block used in transport calculations can be a unit cell, supercell, or conventional cell containing multiple atoms.

**3.1.1. Reading data of force constants from ph.x calculations of QE**

3.1.1.1 Steps to do phonon calculations with QE

Phonon calculations within QE [61,70] are carried out by the module **PHonon**, which employs executables such as **ph.x**, **q2r.x**, **matdyn.x**. These calculations are preceded by relaxation using **pw.x** to obtain converged charge densities for the relaxed structures. The relaxation step reduces the residual forces between the atoms thus seeking for a stable atomic configuration. We assume, furthermore, that a convergence analysis has been performed, for example, in the wave-function and energy cut-off, to obtain reproducible results. The flowchart outlining the steps to obtain phonon bands is illustrated in Fig. 2.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Flowchart for conducting phonon calculations using Quantum Espresso.</figcaption></figure>
We can utilize either the force constants in reciprocal space obtained in step 2 or the interatomic force constants in real space obtained in step 3, derived from a Fourier transform implemented within **q2r.x**. We choose to exploit DFT data in step 2 (force constants in reciprocal space), which grants us greater flexibility in applying corrections (such as the method for Fourier transform and enforcing the acoustic sum rule) to the dynamical matrices in real space, constructing them in a format best suited for the tri-diagonal-blocks lines in the NEGF technique. However, it is noteworthy that a Matlab code has also been developed to read force constants data in real space from step 3 and is available for use.

### 3.1.1.2 Understanding the format of the output files *.dyn1, *.dyn2, ... generated by ph.x calculations

The **ph.x** calculations enable us to compute dynamical matrices in reciprocal space. However, the data written in files *.dyn1, *.dyn2, etc., do not include atomic masses, indicating that these data primarily represent force constants. The arrangement and content of the force constants data files are illustrated in Fig. 3.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Example of the structures of data files containing the force constants in reciprocal space generated in step 2 by the ph.x calculations.</figcaption></figure>

An example of the data for the force constants between two atoms *i* and *j* at a specific *q*-point in the BZ is presented in Fig. 3. The analysis of the data proceeds as follows:

At the *q*-point, for example, (0,0,0), the total coupling between atom 1 and other atoms of the same type is written as:

$$q = ( 0.000000000 0.000000000 0.000000000 )$$

1 1

0.98302408 0.00000000 0.00000000 0.00000000 -0.00000000 0.00000000

0.00000000 0.00000000 0.98302408 0.00000000 -0.00000000 0.00000000

-0.00000000 0.00000000 -0.00000000 0.00000000 0.35261976 0.00000000

In this structure, the total force constant between atom 1 and all neighbor atoms 1 is written in the form:

A B C

D E D
As expected, any pair of atoms has a force-constant matrix of size 3×3. However, due to each term containing real and imaginary parts, the resulting matrix is of size 3×6. Thus, we construct the 3×3 matrix for the force constants between a pair of atoms, (i,j) as follows,

$$K_{(i,j)}(1,1)=Row1(1)+i*Row1(2); K_{(i,j)}(1,2)=Row1(3)+i*Row1(4); K_{(i,j)}(1,3)=Row1(5)+i*Row1(6)$$

$$K_{(i,j)}(2,1)=Row2(1)+i*Row2(2); K_{(i,j)}(2,2)=Row2(3)+i*Row2(4); K_{(i,j)}(2,3)=Row2(5)+i*Row2(6)$$

$$K_{(i,j)}(3,1)=Row3(1)+i*Row3(2); K_{(i,j)}(3,2)=Row3(3)+i*Row3(4); K_{(i,j)}(3,3)=Row3(5)+i*Row3(6)$$

where Row1(k) is the k-th element of the first row associated with the pair (i,j). Similarly, one defines the other elements to construct a matrix $K_{ij}(\mathbf{q}_n)$ of dimensions 3×3 for each q-point.

### 3.1.1.3 Constructing dynamical matrix $D(\mathbf{q})$ of the unit cell

If the unit cell considered within DFT calculations has $N_a$ atoms, we will have $N_a \times N_a$ matrices $K_{ij}$ (3×3) force constant at each q-point. To finally form the matrices of blocks/cells used within the NEGF, we must construct the dynamical matrices of size $3N_a \times 3N_a$. To construct the dynamical matrices $D(\mathbf{q}_n)$ of size $3N_a \times 3N_a$ at any $\mathbf{q}_n$ we make use of the equation [71,72],

$$D(\mathbf{q}_n) = [D_{3 \times 3}^{ij}] = \left[ \frac{K_{ij}(\mathbf{q}_n)}{\sqrt{M_i M_j}} \right], \qquad (14)$$

where $M_i, M_j$ are the atomic masses.

### 3.1.1.4 Matlab codes

A Matlab function, **Read_DynMat_QE.m**, was developed to enable the reading of data from a *.dyn file and the construction of the dynamical matrix for the specified cell at each q-point.

Inputs and outputs of the function **Read_DynMat_QE.m** are:

**Inputs**: the name of the *.dyn file (single file), atomic masses of all atoms in the considered cell.

**Outputs**: List q-points read from the data file, and the $D(\mathbf{q}_n)$ $3N_a \times 3N_a$ matrix associated to each q-point.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Matlab code reads data from *.dyn files and constructs the full dynamical matrix.</figcaption></figure>

A Matlab function, **Read_DynMat_QE_in_individual_files_from_phari計算ation.m**, has also been developed. It utilizes the **Read_DynMat_QE.m** function internally to read data from multiple *.dyn files.

**Note on units**: in the input files *.dyn, the q points are in units $2\pi/\text{alat}$, where *alat* is the lattice constant used by DFT-QE, and the force constants are in units of Ryd/bohr$^2$. The output dynamical matrices $D(\mathbf{q}_n)$ $3N_a \times 3N_a$ are in the SI units, i.e., $N/(\text{kg.m})=1/\text{s}^2$. The q points retain the same units as in the DFT data files, which are in $2\pi/\text{alat}$.
### 3.1.1.5 Validation

One way to test the **Read_DynMat_QE.m** code is to use it for reading the dynamical matrix file generated by the **matdyn.x** calculations (phonon bands or phonon DOSs), using the “fldyn” setup in the input file, and then diagonalize the obtained dynamical matrix $D(q)$ $3N_a \times 3N_a$ at each given q-point to obtain the eigenvalues $\omega^2(q)$ of Eq. (2). Afterwards, the phonon bands can be plotted and compared with those obtained directly from DFT-QE. Figure 5 presents results for two case-studies, 2D graphene and MoS$_2$, demonstrating the reliability of the **Read_DynMat_QE.m** code.

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Validation of the **Read_DynMat_QE.m** code by comparison of DFT-QE phonon bands with phonon bands obtained from the output dynamical matrices $D(q)$ of size $3N_a \times 3N_a$ produced by **Read_DynMat_QE.m**.</figcaption></figure>

### 3.1.2. Fourier transform to obtain dynamical matrices in real space

To be utilized within the lattice Green's function (NEGF technique) for transport calculations, the dynamical matrices must be represented in real space to describe interactions between atoms. To achieve this, Fourier transforms $q \to r$ will be applied to the dynamical matrices $D(q)$.

#### 3.1.2.1 Fourier transforms $q \to r$

First, it is important to note that we perform a Fourier transform on the dynamical matrix of the selected cell to obtain dynamical matrices that describe interactions between the cell (chosen as the central cell or cell 0) and its neighboring cells. This Fourier transform is not for the coupling between two atoms as implemented in the DFT code. The dynamical matrix in real space, showing the interaction between cell 0 and cell $\alpha$ at the position defined by the vector $\mathbf{R}_{0\alpha}$ can be calculated by

$$D(\mathbf{R}_{0\alpha}) = \frac{1}{N_q} \sum_{i=1}^{N_q} D(q_i) e^{-i\mathbf{q}_i \cdot \mathbf{R}_{0\alpha}} \quad (15)$$
With {$D(\mathbf{q}_i), \mathbf{q}_i$} are obtained by **Read_DynMat_QE.m** and **Read_DynMat_QE_in Brady Files from ph calculation.m.** $R_{0\alpha}$ can be determined by

$$R_{0\alpha} = i_1 a_1 + i_2 a_2 + i_3 a_3 \qquad (16)$$

in which $i_1, i_2, i_3$ are integer numbers and $a_1, a_2, a_3$ are the lattice vectors whose information can be found in the setup of DFT calculations.

**Cut-off in real space for $R_{0\alpha}$:** To be consistent with the DFT-QE implementation for Fourier transforms, the cut-off of $R_{0\alpha}$ in the Fourier transform in Eq. (15) will be defined via the range of values of $i_1, i_2, i_3$ that must be consistent with the setup of the q-mesh (via the parameters $nq_1, nq_2, nq_3$) in the DFT-QE **ph.x** calculations, i.e.,
$$i_1 = \left[-\frac{nq_1}{2}, \frac{nq_1}{2}\right], i_2 = \left[-\frac{nq_2}{2}, \frac{nq_2}{2}\right], i_3 = \left[-\frac{nq_3}{2}, \frac{nq_3}{2}\right].$$

### 3.1.2.2 Checking $D(R_{0\alpha})$ obtained from Fourier transforms

To check if $D(R_{0\alpha})$ obtained from the Fourier transform are correct, we can use them to construct again $D(\mathbf{q}_i)$ (or it can say doing an inverse Fourier transform):

$$D(\mathbf{q}_i) = \sum_{\alpha} D(R_{0\alpha}) e^{i\mathbf{q}_i \cdot R_{0\alpha}} \qquad (17)$$

and then diagonalizing $D(\mathbf{q}_i)$ to obtain the phonon bands for comparison with the DFT-QE result. Results obtained by our method compared with the DFT-QE ones are shown in Fig. 6.

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Phonon bands were obtained from $D(R_{0\alpha})$ after Fourier transforms in the proposed approach (black lines) and the ones obtained by DFT-QE (color lines).</figcaption></figure>
### 3.1.2.3 Correcting diagonal block matrices: Acoustic sum rule

As depicted in Fig. 6, the acoustic waves at the Gamma point, q = (0,0,0), in our approach exhibit divergence, whereas the DFT results correctly show convergence of these waves. This discrepancy arises because our approach does not yet incorporate a correction for the acoustic sum rule (ASR). Specifically, after Fourier transforms, the diagonal 3×3 blocks in the matrix $D_{00} \equiv D(\mathbf{R}_{00})$ must be adjusted according to the equation [71,72],

$$D = [D_{3 \times 3}^{ij}] = \begin{bmatrix} -\frac{K_{ij}}{\sqrt{M_i M_j}} & \text{for } j \neq i \\ \sum_{j \neq i} \frac{K_{ij}}{M_j} & \text{for } j = i \end{bmatrix}, \qquad (18)$$

Thus, $D_{00}^{ii}(3 \times 3) = \sum_{j \neq i} \frac{K_{ij}}{M_i}$. It should be noticed that the sign “-” in Eq. (18) was already included in data $D(\mathbf{R}_{0\alpha})$ obtained by the Fourier transform, so a correction for the diagonal blocks 3×3 of $D(\mathbf{R}_{0\alpha})$ must be: $D_{00}^{ii}(3 \times 3) = -\sum_{j \neq i} \frac{K_{ij}}{M_i}$ or $K_{ii} = -\sum_{j \neq i} K_{ij}$ with $K_{ij} = \sqrt{M_i M_j} D_{3 \times 3}^{ij}$. This is identical to the acoustic sum rule applied within the DFT calculations.

To correct the diagonal-blocks 3×3 of $D(\mathbf{R}_{0\alpha})$, we developed a Matlab function `Acoustic_Sum_check_diaginal_blocks_of_DynMat_real_space.m`. When such a correction is applied, the new results of phonon bands are shown in Fig. 7.

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Phonon bands obtained using our approach, including the acoustic sum rule, compared with those obtained using DFT-QE for two systems: (a) a chain of B and N atoms with two atoms per cell, (b) 2D graphene sheets.</figcaption></figure>
3.1.2.4. Important note: using odd-q points in ph.x calculations

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Comparing the phonon bands of MoS<sub>2</sub> obtained by the two approaches.</figcaption></figure>

The phonon bands obtained by our approach show good overall agreement with the DFT data, as depicted in Fig. 7. However, discrepancies are still noticeable, particularly highlighted in Fig. 8 for the case of MoS<sub>2</sub>.

To understand the origin of these discrepancies and improve the accuracy of our approach, we compare the methods used for Fourier transform in the DFT-QE code and our approach. A summary is provided in Fig. 9:

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Comparison of the real-space zones where Fourier transforms are applied in DFT-QE and our approach.</figcaption></figure>

It is worth mentioning that in our code, $N_1$, $N_2$, $N_3$ represent the counts of the maximum nearest neighbor cells of the central cell 0 along each direction. These values are integers and should satisfy: $N_1 \le nq_1/2$, $N_2 \le nq_2/2$, $N_3 \le nq_3/2$.

A direct observation of Fig. 9 suggests that if we choose odd *nq* in DFT calculations, the results of our method should be closer to those of DFT-QE. Upon recalculating results for both approaches, a significant improvement can be observed in Fig. 10 compared to that shown in Fig. 8.
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Our approach reproduces very well the DFT-QE results if data in **ph.x** calculation were obtained with the setup of odd *nq*.</figcaption></figure>

In summary, it is crucial to set *nq* as odd numbers in **ph.x** calculations to maintain symmetry along each direction for the interaction between a central cell and its neighboring cells, thereby ensuring that the output from our approach closely agrees with the results from DFT-QE.

**3.2. Methods for phonon transport based on DFT-QE + NEGF**

In this section, we will focus on phonon transport using DFT-QE + NEGF. The dynamical matrices $D(\mathbf{R}_{0\alpha})$ obtained in the previous steps will be utilized to construct a dynamical matrix that adheres to the tridiagonal-block lines within the NEGF technique, where each *n*-th block interacts only with its nearest neighbors *n*-1 and *n*+1 (see Fig. 1 as a reminder).

**3.2.1. Two approaches to construct dynamical matrices adapted to tri-diagonal NEGF technique**

As seen in Fourier transform, $N_1, N_2, N_3$ are usually larger than 1, indicating that a primary cell interacts not only with its nearest neighbor cells but also with next-nearest or further neighbor cells. Therefore, in most cases, primary cells cannot be directly used within the NEGF framework, necessitating the construction of supercells (from primary cells) or conventional cells that are sufficiently large to ensure that each cell/block interacts only with its nearest neighbors along each direction.

It can be noted that a supercell lattice is a specific case of conventional cells. Henceforth, we will generally refer to them as conventional cells.

To obtain the dynamical matrices $D(\mathbf{R}_{0\alpha})$ of conventional cells, there are two methods as illustrated in Fig. 11.
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Two approaches to get $D(\mathbf{R}_{0\alpha})$ of cells satisfied the tri-diagonal block lines in the NEGF technique.</figcaption></figure>

In Method 1, further DFT calculations are required for supercells/conventional cells, which can be computationally expensive due to the larger cell size. In contrast, Method 2 does not necessitate additional DFT calculations but requires careful handling to accurately derive dynamical matrices for supercells/conventional cells from those of the primary cells.

In practice, the rapid increase in computation cost with the size of the cell considered within DFT methods is illustrated in the example shown in Fig. 12 for graphene. Due to limited computational resources, it is often impractical to complete DFT calculations when the size of the cell is relatively large, even with just a few tens of atoms. Therefore, to enable full ab initio phonon transport coefficient calculations based on DFT-QE + NEGF for large cell systems, it is essential to develop a general method to directly construct dynamical matrices $D(\mathbf{R}_{0\alpha})$ of conventional cells from those of the primary cells.

<figure><img src="image_12.png" /><figcaption>Method 1: Computing $D(r)$ of conventional cells by DFT-QE</figcaption></figure>

<caption>**Figure 12.** Computation cost for DFT-QE ph.x calculations with different cell sizes. Machine: Ganymede (C2N, University Paris Saclay), using simultaneously 30 CPUs.</caption>
### 3.2.2. Fast dynamical matrix construction approach for conventional cells

In this subsection, we present the development of a general method that allows for the construction of dynamical matrices for conventional cells directly from dynamical matrices of primary cells. It is noteworthy that this method can also be extended to electron systems, where Hamiltonians of primary cells are derived from DFT calculations or known in the Tight-Binding method.

In the primary-cell lattice, the vector from cell 0 to a neighbor cell can be defined as $\mathbf{R}_{0,m_1m_2m_3} = m_1\mathbf{a}_1 + m_2\mathbf{a}_2 + m_3\mathbf{a}_3$. We assume that dynamical matrices and force constants between atoms and cells are known in the primary cell lattice, obtained through Fourier transforms in the previous steps (section 3.1). Thus, we have the information of all the 3×3 matrices $K(m_1, m_2, m_3, n_a, n_b)$ presenting the force constant coupling between atom n_a-th in cell 0, and atom n_b-th in cell $(m_1, m_2, m_3)$ in the primary-cell lattice.

In the conventional-cell lattice, the lattice vectors are now denoted as $\mathbf{a}_1', \mathbf{a}_2', \mathbf{a}_3'$. The vector from the cell 0' to a neighbor cell is defined by $\mathbf{R}_{0,m_1'm_2'm_3'} = m_1'\mathbf{a}_1' + m_2'\mathbf{a}_2' + m_3'\mathbf{a}_3'$. To define the dynamical matrix describing the coupling of cell 0' with all its neighboring cells, we need to specify all force constant couplings $K(m_1', m_2', m_3', n_a', n_b')$ between atom n_a'-th in cell 0', and atom n_b'-th in cell $(m_1', m_2', m_3')$. The vector between these two atoms can be written as $\mathbf{r} = xyz(m_1', m_2', m_3', n_a', n_b') - xyz(m_1' = 0, m_2' = 0, m_3' = 0, n_a' = 0)$. If one can determine the cell indices $(m_1, m_2, m_3)$ in the primary-cell lattice of atom n_b'-th compared to atom n_a'-th, then $K(m_1', m_2', m_3', n_a', n_b')$ can be identified to $K(m_1, m_2, m_3, n_a, n_b)$ which is already known. The set of indices $(m_1, m_2, m_3)$ in the primary-cell lattice of the vector $\mathbf{r}$ can be defined by the relationship

$$ \mathbf{r} = m_1\mathbf{a}_1 + m_2\mathbf{a}_2 + m_3\mathbf{a}_3, \quad (19) $$

which thus can be written in the matrix form:

$$ \begin{pmatrix} r_x \\ r_y \\ r_z \end{pmatrix} = \begin{pmatrix} a_{1x} & a_{2x} & a_{3x} \\ a_{1y} & a_{2y} & a_{3y} \\ a_{1z} & a_{2z} & a_{3z} \end{pmatrix} \begin{pmatrix} m_1 \\ m_2 \\ m_3 \end{pmatrix} \Rightarrow \begin{pmatrix} m_1 \\ m_2 \\ m_3 \end{pmatrix} = A^{-1} \begin{pmatrix} r_x \\ r_y \\ r_z \end{pmatrix}, \quad (20) $$

$$ \text{where } A = \begin{pmatrix} a_{1x} & a_{2x} & a_{3x} \\ a_{1y} & a_{2y} & a_{3y} \\ a_{1z} & a_{2z} & a_{3z} \end{pmatrix}. $$

We summarize the steps to determine all $K(m_1', m_2', m_3', n_a', n_b')$ in the conventional-cell lattice as follows:

**Step 1: Construct information of cell 0 and cell 0' in the primary-cell and conventional-cell lattices**

We construct matrices to contain the information of the coordinates, the label of each atom by their order in the primary cell. For the conventional-cell lattice, the labels and the cell index of atoms are the labels in the primary cell, thus will help us to identify the atoms in the reference lattice (primary-cell lattice):
$$ \mathbf{X}_0 = \begin{pmatrix} x & y & z & \text{label} \end{pmatrix}, \quad \mathbf{X}'_0 = \begin{pmatrix} x & y & z & \text{label} \end{pmatrix} $$

Primary cell
2 atoms/cell

<figure><img src="image_13.png" /></figure>

a_1

a_2

a_3

conventional cell
8 atoms/cell

<figure><img src="image_14.png" /></figure>

a'_1

a'_2

a'_3

For example, in the primary cell of graphene, the atoms are labeled 1 and 2 within a unit cell. Correspondingly, in the conventional cell, these atoms are also labeled 1 and 2, as illustrated in the figure above.

**Step 2: Determine the cell indices ($m_1, m_2, m_3$) in the primary cell lattice that all the atoms in cell 0' belong to i-th atom in cell 0', for a specific order n is actually a translation of the atom n-th (the same order level) in the cell 0. Therefore, the cell indices ($m_1, m_2, m_3$) that the atom i-th of cell 0' belong to can be defined by the equation:**

$$ \mathbf{r} = \mathbf{X}'_0 (-i\mathbf{th}) (x, y, z, \text{order=n}) - \mathbf{X}'_{0-n\mathbf{th}} (x, y, z, \text{order=n}) $$

$$ (m_1, m_2, m_3)_{i-th} = A^{-1} \mathbf{r}^T $$

Thus, we can define the cell indices ($m_1, m_2, m_3$) of all the atoms in the cell 0. These parameters can be added as new columns in the matrix $\mathbf{X}'_0$:

$$ \mathbf{X}'_0 = \begin{pmatrix} x & y & x & \text{label} & m_1 & m_2 & m_3 \end{pmatrix} $$

This matrix has the dimensions $N_a \times 7$.

**Step 3: Do a translation to determine position of the atoms in any cell ($m_1', m_2', m_3'$) in the conventional-cell lattice**

The information of atoms in any cell ($m_1', m_2', m_3'$) can be determined by a translation

$$ \mathbf{R}'_{0,m_1', m_2', m_3'} = m_1' \mathbf{a}'_1 + m_2' \mathbf{a}'_2 + m_3' \mathbf{a}'_3 \quad \text{based on the information of the cell 0' written in } \mathbf{X}'_0 $$

$$ \text{xyz}(m_1', m_2', m_3', n_b', \text{label } n) = \text{xyz}(m_1' = 0, m_2' = 0, m_3' = 0, n_b', \text{label } n) + \mathbf{R}'_{0,m_1', m_2', m_3'} $$

**Step 4: Determine the relative vector between the atom $n_b'$ of the cell ($m_1', m_2', m_3'$) and the atom $n_a'$ of the cell ($m_1'=0, m_2'=0, m_3'=0$) in the conventional-cell lattice**

$$ \mathbf{r} = \text{xyz}(m_1', m_2', m_3', n_b', \text{label}=n_2) - \text{xyz}(m_1' = 0, m_2' = 0, m_3' = 0, n_a', \text{label}=n_1) $$

**Step 5: Determine the relative cell index ($m_1, m_2, m_3$) in the primary-cell lattice of the atom $n_b'$ of the cell ($m_1'=0, m_2'=0, m_3'=0$) compared to the atom $n_a'$ of the cell ($m_1'=0, m_2'=0, m_3'=0$)**
$$ \begin{pmatrix} m_1 \\ m_2 \\ m_3 \end{pmatrix} = A^{-1} \begin{pmatrix} r_x \\ r_y \\ r_z \end{pmatrix} $$

Step 6: Identify the force constant coupling between atom n<sub>b</sub>' of the cell (m<sub>1</sub>, m<sub>2</sub>, m<sub>3</sub}') and n<sub>a</sub>' of the cell 0' to the one of atom n<sub>2</sub> of cell (m<sub>1</sub>, m<sub>2</sub>, m<sub>3</sub>) and n<sub>1</sub> of cell 0

$$ K(m_1', m_2', m_3', n_a', n_b') \equiv K(m_1, m_2, m_3, n_1, n_2) $$

Once the force constants between atoms in cell 0' and atoms in other neighboring cells (including interactions between atoms within cell 0') are defined in the conventional cell lattice, it becomes possible to construct dynamical matrices for cell 0' and between this cell and its nearest neighbor cells by Eq. (18).

### 3.3. Notes on the codes and the DFT data

The DFT data for all structures used in this study are provided. As the primary focus of this work is on method development, detailed information on the setup of DFT calculations is not within the scope here. Specific details can be found in the input files used for each calculation. The Matlab codes developed for Fourier transforms and for constructing dynamical matrices of conventional cells from primary cells are tailored for each structure to facilitate practical implementation, with some input parameters imported manually. However, more generalized codes capable of handling various structures can be easily derived from these codes. The DFT data and Matlab codes can be accessed at the following link: Link_Codes_DFTdata.

Figure 13 summarizes the main steps of the proposed method for computing thermal transport.

<figure><img src="image_15.png" /><figcaption>Figure 13. summary of steps in the proposed approach for phonon transport based on DFT-QE + NEGF.</figcaption></figure>

### 4. Validation of the methodology for various systems

In this section, we validate the method presented in Section 3 by examining phonon transport in 1D graphene nanoribbons, 2D graphene sheets, and 3D bulk silicon. Additionally, data from the literature were extracted for comparative analysis.
## 4.1. One-dimensional systems: graphene ribbons

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Thermal conductance of 1D structures: graphene ribbons. Data from Tan [73] [73] were extracted for comparison.</figcaption></figure>

Initially, we applied the method to armchair graphene ribbons (AGNRs) with widths of $M = 5$ and 6 dimer lines, respectively, where each primary cell contains $2M$ atoms. Dynamical matrices for the primary cell and their couplings with neighboring cells were obtained, and conventional cells twice the size of the primary cells were constructed. The dynamical matrices of these conventional cells were then computed using data from the primary cells and subsequently utilized in NEGF calculations. The thermal conductance as a function of temperature is shown in Fig. 14. Data at 300 K from Tan et al. [73] were extracted for comparison, where the DFT Siesta method and code were utilized in their work. The thermal conductance results at 300 K from both methods exhibit good agreement, validating the proposed approach. It is also worth noting that in the ballistic regime, at higher temperatures, more phonon modes contribute to transport, leading to an increase in thermal conductance with rising temperature. Additionally, the number of phonon bands is proportional to the number of atoms in the unit cell, meaning a phonon dispersion contains $3N_{at}$ bands for a unit cell with $N_{at}$ atoms. Consequently, a wider ribbon exhibits higher conductance due to the increased number of phonon modes contributing to transport [74].

## 4.2. Two-dimensional systems: 2D graphene

We further validate the developed method using 2D structures, specifically graphene sheets. In Fig. 15(a), phonon transmission obtained by different methods is presented: solid black represents the Force Constant model (FC) and NEGF with force constant parameters from reference [14], solid red shows results from DFT calculations of dynamical matrices extracted from conventional cells containing 8 atoms per cell. The dashed blue line corresponds to the fast dynamical matrix construction approach outlined in section 3.2, using dynamical matrices of conventional cells with 8 atoms per cell constructed from primary cell data containing 2 atoms per cell. The results demonstrate good agreement among the different methods.

In Fig. 15(b), we computed the thermal conductance per unit surface perpendicular to the transport direction, $S = L_y \times d$ where $d = 0.34$ nm is the thickness of the graphene layer, $L_y = 0.426$ nm is the lattice constant along the y axis. Our approach results are compared with those from the work by Jiang et al. [75], which employed the valence force field model (VFFM). The results exhibit close agreement over a wide range of temperatures.
<figure><img src="image_17.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_18.png" /></figure>

**Figure 15.** (a) Total phonon transmission $T(\omega)$ calculated by 3 methods: force constant model (solid black), DFT with $D(r)$ obtained by Fourier transform for the conventional cells of 8 atoms that were used in DFT-QE calculations (solid red), and constructed from the DFT results of the primary cells of 2 atoms (dashed blue). (b) Thermal conductance per unit area of the crossing surface. The curve with symbols presents data from the work by Jiang et al. [75].

**4.3. Bulk (3D) systems: Silicon**

Finally, we validated our method for bulk structures, focusing on silicon. Phonon transport calculations involved constructing conventional cells containing 8 atoms, with dynamical matrices derived from DFT data of the primary cell (2 atoms/cell). In Fig. 16, the first row presents phonon transmission results from our approach (right panel) and compares them with those from Tian et al. [76], who utilized MD combined with the NEGF technique. It is important to note that due to unknown factors affecting the transmission function in ref. [70], comparisons between Figs. 16(a) and (b) should focus solely on the shape of the transmission profiles rather than exact numerical values. Nevertheless, both methods exhibit good qualitative agreement.

In Fig. 16(c), the thermal conductance per unit crossing area ($S = L_x \times L_y = 0.543^2 \text{ nm}^2$) is shown. It can be observed that at low temperatures ($T < 130$ K), the MD approach provides results close to our DFT + NEGF method. Above 130 K, however, the MD approach overestimates the thermal conductance of silicon. This difference can be attributed to the phonon bands obtained
by the MD method, which extend up to 600 cm⁻¹ (see Fig. 16(a)), while the frequency range obtained by DFT is shorter, not exceeding 510 cm⁻¹.

<figure><img src="image_19.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_21.png" /></figure>

**Figure 16.** Transmission and thermal conductance of purely bulk silicon. Data from Tian et al. [76] were extracted from comparison.

**5. Conclusion**

We have presented a detailed method for studying phonon transport using outputs from DFT-QE phonon calculations combined with the NEGF technique. Initially, we conducted a thorough analysis of the force constants obtained from phonon calculations and subsequently developed a Matlab code to extract and organize the dynamical matrix data for the entire unit cell at each q-point. A complete method then was developed to do Fourier transforms to obtain the dynamical matrices in real space, essential for describing interactions between cells required by the NEGF technique. The method included a correction of diagonal block matrices from the acoustic sum rule as well as a demonstration of the necessity of using odd-q meshes for accurate applications within NEGF.

A vital part was then developed to present an efficient approach to construct dynamical matrices of conventional cells directly from those of primary cells, ensuring compliance with the tridiagonal rule of the NEGF technique. This approach significantly reduces the computational burden associated with large cell sizes in DFT calculations.
We validated our methodology by applying it to phonon transport in various structures: from 1D structures such as armchair graphene ribbons, to 2D structures like graphene sheets, and 3D structures such as bulk silicon. Each case demonstrated the capability of our method to accurately predict phonon transport properties and its agreement with existing literature, affirming its reliability and effectiveness.

The proposed routine for computing phonon transport offers a robust framework for studying thermal transport in diverse materials with accuracy at the level of DFT calculations. Moreover, the fast dynamical matrix construction approach for conventional cells enables the construction of dynamical matrices for significantly larger cells. This capability facilitates the introduction of vacancies at reasonable concentrations, which would be unrealistic with smaller cells typically used in DFT calculations. An approximation, such as setting the mass of vacancies to infinity (effectively a very large number, numerically), allows for the study of the impact of defects on phonon transport properties.

Looking forward, the methodology presented in this work can also be extended to study time-resolved, nonequilibrium phonon transport phenomena. By integrating the DFT + NEGF framework with time-dependent perturbations, our approach could be adapted to investigate transient phonon dynamics in response to pulsed excitations or varying thermal gradients. This aligns with recent advances in the field, where time-resolved techniques have provided new insights into phonon dynamics in graphene ribbons, atomic wires, and molecular junctions [77-82]. Expanding into this area would open new avenues for exploring dynamic thermal transport in nanostructures and deepen our understanding of phonon behavior under non-equilibrium conditions, a topic of growing interest in both experimental and theoretical research.

## References

[1] D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, A. Majumdar, H. J. Maris, R. Merlin, and S. R. Phillpot, *Nanoscale Thermal Transport*, J Appl Phys **93**, 793 (2003).

[2] Y. Zhao, Y. Cai, L. Zhang, B. Li, G. Zhang, and J. T. L. Thong, *Thermal Transport in 2D Semiconductors—Considerations for Device Applications*, Adv Funct Mater **30**, (2020).

[3] Y.-J. Zeng, Y.-Y. Liu, W.-X. Zhou, and K.-Q. Chen, *Nanoscale Thermal Transport: Theoretical Method and Application*, Chinese Physics B **27**, 036304 (2018).

[4] L. Shi, *Thermal and Thermoelectric Transport in Nanostructures and Low-Dimensional Systems*, Nanoscale and Microscale Thermophysical Engineering **16**, 79 (2012).

[5] F. Duan, D. Wei, A. Chen, X. Zheng, H. Wang, and G. Qin, *Efficient Modulation of Thermal Transport in Two-Dimensional Materials for Thermal Management in Device Applications*, Nanoscale **15**, 1459 (2023).

[6] E. Pop, *Energy Dissipation and Transport in Nanoscale Devices*, Nano Res **3**, 147 (2010).
[7] K. M. Hoogeboom-Pot et al., *A New Regime of Nanoscale Thermal Transport: Collective Diffusion Increases Dissipation Efficiency*, Proc Natl Acad Sci U S A **112**, (2015).

[8] A. Shakouri, *Recent Developments in Semiconductor Thermoelectric Physics and Materials*, Annu Rev Mater Res **41**, (2011).

[9] A. I. Hochbaum, R. Chen, R. D. Delgado, W. Liang, E. C. Garnett, M. Najarian, A. Majumdar, and P. Yang, *Enhanced Thermoelectric Performance of Rough Silicon Nanowires*, Nature **451**, (2008).

[10] M. S. Dresselhaus, G. Chen, M. Y. Tang, R. Yang, H. Lee, D. Wang, Z. Ren, J. P. Fleurial, and P. Gogna, *New Directions for Low-Dimensional Thermoelectric Materials*, Advanced Materials **19**, 1043 (2007).

[11] N. Kazeev, A. R. Al-Maeeni, I. Romanov, M. Faleev, R. Lukin, A. Tormasov, A. H. Castro Neto, K. S. Novoselov, P. Huang, and A. Ustyuzhanin, *Sparse Representation for Machine Learning the Properties of Defects in 2D Materials*, NPJ Comput Mater **9**, (2023).

[12] W. Ye, C. Chen, Z. Wang, I. H. Chu, and S. P. Ong, *Deep Neural Networks for Accurate Predictions of Crystal Stability*, Nat Commun **9**, (2018).

[13] K. Burke, *Perspective on Density Functional Theory*, Journal of Chemical Physics **136**, (2012).

[14] L. Wirtz and A. Rubio, *The Phonon Dispersion of Graphite Revisited*, Solid State Commun **131**, 141 (2004).

[15] J. Zimmermann, P. Pavone, and G. Cuniberti, *Vibrational Modes and Low-Temperature Thermal Properties of Graphene and Carbon Nanotubes: Minimal Force-Constant Model*, Phys Rev B Condens Matter Mater Phys **78**, 1 (2008).

[16] A. P. P. Nicholson and D. J. Bacon, *A New Force-Constant Model for Graphite*, Journal of Physics C: Solid State Physics **10**, 2295 (1977).

[17] V.-T. Tran, J. Saint-Martin, P. Dollfus, and S. Volz, *Third Nearest Neighbor Parameterized Tight Binding Model for Graphene Nano-Ribbons*, AIP Adv **7**, 075212 (2017).
[18] K. Zhao, M. Zhao, Z. Wang, and Y. Fan, *Tight-Binding Model for the Electronic Structures of SiC and BN Nanoribbons*, Physica E Low Dimens Syst Nanostruct **43**, 440 (2010).

[19] Z. Liu, G. Yu, H. Yao, L. Liu, L. Jiang, and Y. Zheng, *A Simple Tight-Binding Model for Typical Graphyne Structures*, New J Phys **14**, (2012).

[20] P. Vogl, H. P. Hjalmarson, and J. D. Dow, *A Semi-Empirical Tight-Binding Theory of the Electronic Structure of Semiconductors†*, Journal of Physics and Chemistry of Solids **44**, 365 (1983).

[21] Y. W. Son, M. L. Cohen, and S. G. Louie, *Energy Gaps in Graphene Nanoribbons*, Phys Rev Lett **97**, 1 (2006).

[22] G. Bouzerar, S. Thébaud, S. Pecorario, and C. Adessi, *Drastic Effects of Vacancies on Phonon Lifetime and Thermal Conductivity in Graphene*, Journal of Physics Condensed Matter **32**, (2020).

[23] G. K. H. Madsen and D. J. Singh, *BoltzTraP. A Code for Calculating Band-Structure Dependent Quantities*, Comput Phys Commun **175**, 67 (2006).

[24] A. Cepellotti and N. Marzari, *On the Kinetic Theory of Thermal Transport in Crystals*, in *Handbook of Materials Modeling: Applications: Current and Emerging Materials*, Second Edition (2020).

[25] N. Neophytou, *Boltzmann Transport Method for Electronic Transport in Complex Bandstructure Materials*, in *SpringerBriefs in Physics*, Vol. Part F1003 (2020).

[26] J. I. Ramos, *The Boltzmann Equation and Its Applications*, Appl Math Model **14**, (1990).

[27] A. Nabovati, D. P. Sellan, and C. H. Amon, *On the Lattice Boltzmann Method for Phonon Transport*, J Comput Phys **230**, (2011).

[28] L. E. F. F. Torres, S. Roche, and J. C. Charlier, *Introduction to Graphene-Based Nanomaterials: From Electronic Structure to Quantum Transport* (Cambridge University Press, 2014).

[29] L. Lindsay, A. Katre, A. Cepellotti, and N. Mingo, *Perspective on Ab Initio Phonon Thermal Transport*, J Appl Phys **126**, (2019).
[30] M. Raya-Moreno, X. Cartoixà, and J. Carrete, *BTE-Barna: An Extension of AlmaBTE for Thermal Simulation of Devices Based on 2D Materials*, Comput Phys Commun **281**, (2022).

[31] *Quantum Kinetics in Transport and Optics of Semiconductors* (2008).

[32] M. Bonitz, *Quantum Kinetic Theory* (2016).

[33] S. Datta, *Quantum Transport: Atom to Transistor* (Cambridge University Press, Cambridge, 2005).

[34] Q. Mao, M. Feng, X. Z. Jiang, Y. Ren, K. H. Luo, and A. C. T. van Duin, *Classical and Reactive Molecular Dynamics: Principles and Applications in Combustion and Energy Systems*, Prog Energy Combust Sci **97**, 101084 (2023).

[35] J. Chen, *The Development and Comparison of Molecular Dynamics Simulation and Monte Carlo Simulation*, in *IOP Conference Series: Earth and Environmental Science*, Vol. 128 (2018).

[36] F. Djurabekova and K. Nordlund, *Molecular Dynamics Simulations of Non-Equilibrium Systems*, in *Handbook of Materials Modeling: Applications: Current and Emerging Materials*, Second Edition (2020).

[37] G. Mancini, M. Fusè, F. Lipparini, M. Nottoli, G. Scalmani, and V. Barone, *Molecular Dynamics Simulations Enforcing Nonperiodic Boundary Conditions: New Developments and Application to the Solvent Shifts of Nitroxide Magnetic Parameters*, J Chem Theory Comput **18**, (2022).

[38] R. Namakian, W. J. Meng, and D. Moldovan, *Role of Geometry and Coherent Twin Boundaries in Mechanical Response of Cu(110) Nanopillars under Tensile Loading: Insights from Molecular Dynamics Simulations*, J Mater Sci **59**, 7489 (2024).

[39] F. Khalili-Araghi, J. Gumbart, P.-C. Wen, M. Sotomayor, E. Tajkhorshid, and K. Schulten, *Molecular Dynamics Simulations of Membrane Channels and Transporters*, Curr Opin Struct Biol **19**, 128 (2009).

[40] P. Zhang, F. Jiao, L. Wu, Z. Kong, W. Hu, L. Liang, and Y. Zhang, *Molecular Dynamics Simulation of Transport Mechanism of Graphene Quantum Dots through Different Cell Membranes*, Membranes (Basel) **12**, (2022).
[41] M. S. Badar, S. Shamsi, J. Ahmed, and Md. A. Alam, *Molecular Dynamics Simulations: Concept, Methods, and Applications*, in (2022).

[42] R. Schneider, A. R. Sharma, and A. Rai, *Introduction to Molecular Dynamics*, Lecture Notes in Physics.

[43] J. Jung, C. Tan, and Y. Sugita, *GENESIS CGDYN: Large-Scale Coarse-Grained MD Simulation with Dynamic Load Balancing for Heterogeneous Biomolecular Systems*, Nat Commun **15**, 3370 (2024).

[44] T. Mueller, A. Hernandez, and C. Wang, *Machine Learning for Interatomic Potential Models*, J Chem Phys **152**, (2020).

[45] G. Kumar, R. R. Mishra, and A. Verma, *Introduction to Molecular Dynamics Simulations*, in *Lecture Notes in Applied and Computational Mechanics*, Vol. 99 (2022).

[46] E. Prašnikar, M. Ljubič, A. Perdih, and J. Borišek, *Machine Learning Heralding a New Development Phase in Molecular Dynamics Simulations*, Artif Intell Rev **57**, (2024).

[47] T. Cui, C. Tang, M. Su, S. Zhang, Y. Li, L. Bai, Y. Dong, X. Gong, and W. Ouyang, *Geometry-Enhanced Pretraining on Interatomic Potentials*, Nat Mach Intell **6**, 428 (2024).

[48] S. Datta, *Nanoscale Device Modeling: The Green's Function Method*, Superlattices Microstruct **28**, 253 (2000).

[49] X. Chen, Y. Liu, and W. Duan, *Thermal Engineering in Low-Dimensional Quantum Devices: A Tutorial Review of Nonequilibrium Green's Function Methods*, Small Methods **2**, 1700343 (2018).

[50] C. H. Lewenkopf and E. R. Mucciolo, *The Recursive Green's Function Method for Graphene*, J Comput Electron **12**, 203 (2013).

[51] M. Pourfath, *The Non-Equilibrium Green's Function Method for Nanoscale Device Simulation*, Vol. 2014 (Springer Vienna, Vienna, 2014).

[52] J. S. Wang, J. Wang, and J. T. Lü, *Quantum Thermal Transport in Nanostructures*, European Physical Journal B **62**, 381 (2008).
[53] Y. Guo, M. Bescond, Z. Zhang, M. Luisier, M. Nomura, and S. Volz, *Quantum Mechanical Modeling of Anharmonic Phonon-Phonon Scattering in Nanostructures*, Phys Rev B **102**, (2020).

[54] Y. Guo, M. Bescond, Z. Zhang, S. Xiong, K. Hirakawa, M. Nomura, and S. Volz, *Thermal Conductivity Minimum of Graded Superlattices Due to Phonon Localization*, APL Mater **9**, (2021).

[55] Y. Wu and P. a. Childs, *Conductance of Graphene Nanoribbon Junctions and the Tight Binding Model*, Nanoscale Res Lett **6**, 1 (2011).

[56] T. Ozaki, K. Nishio, and H. Kino, *Efficient Implementation of the Nonequilibrium Green Function Method for Electronic Transport Calculations*, Phys Rev B **81**, 035116 (2010).

[57] M. P. L. Sancho, J. M. L. Sancho, and J. Rubio, *Quick Iterative Scheme for the Calculation of Transfer Matrices: Application to Mo (100)*, Journal of Physics F: Metal Physics **14**, 1205 (2000).

[58] Y. Yoon, G. Fiori, S. Hong, G. Iannaccone, and J. Guo, *Performance Comparison of Graphene Nanoribbon FETs with Schottky Contacts and Doped Reservoirs*, IEEE Trans Electron Devices **55**, (2008).

[59] K. Rallis, P. Dimitrakis, I. G. Karafyllidis, A. Rubio, and G. C. Sirakoulis, *Electronic Properties of Graphene Nanoribbons with Defects*, IEEE Trans Nanotechnol **20**, (2021).

[60] Y.-M. Niquet, V.-H. Nguyen, F. Triozon, I. Duchemin, O. Nier, and D. Rideau, *Quantum Calculations of the Carrier Mobility: Methodology, Matthiessen's Rule, and Comparison with Semi-Classical Approaches*, J Appl Phys **115**, 054512 (2014).

[61] P. Giannozzi et al., *Quantum ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials*, J Phys Condens Matter **21**, 395502 (2009).

[62] G. Mil'nikov, J. ichi Iwata, N. Mori, and A. Oshiyama, *RSDFT-NEGF Transport Simulations in Realistic Nanoscale Transistors*, J Comput Electron **22**, (2023).

[63] P. S. Drouvelis, P. Schmelcher, and P. Bastian, *Parallel Implementation of the Recursive Green's Function Method*, J Comput Phys **215**, 741 (2006).
[64] V.-T. Tran, J. Saint-Martin, and P. Dollfus, *High Thermoelectric Performance in Graphene Nanoribbons by Graphene/BN Interface Engineering*, Nanotechnology **26**, 495202 (2015).

[65] Y. Wu and P. A. Childs, *Conductance of Graphene Nanoribbon Junctions and the Tight Binding Model*, Nanoscale Res Lett **6**, 1 (2011).

[66] C. H. Lewenkopf and E. R. Mucciolo, *The Recursive Green's Function Method for Graphene*, J Comput Electron **12**, 203 (2013).

[67] S. Datta, *Nanoscale Device Modeling: The Green's Function Method*, Superlattices Microstruct **28**, 253 (2000).

[68] V.-T. Tran, Electronic and Thermoelectric Properties of Graphene/Boron Nitride in-Plane Heterostructures, doctoral thesis, Université Paris-Saclay, 2015.

[69] T. Yamamoto and K. Watanabe, *Nonequilibrium Green's Function Approach to Phonon Transport in Defective Carbon Nanotubes*, Phys Rev Lett **96**, 255503 (2006).

[70] P. Giannozzi et al., *Advanced Capabilities for Materials Modelling with Quantum ESPRESSO*, Journal of Physics: Condensed Matter **29**, 465901 (2017).

[71] V.-T. Tran, J. Saint-Martin, and P. Dollfus, *High Thermoelectric Performance in Graphene Nanoribbons by Graphene/BN Interface Engineering*, Nanotechnology **26**, 495202 (2015).

[72] V.-T. Tran and A. Cresti, *Thermoelectric Properties of In-Plane 90°-Bent Graphene Nanoribbons with Nanopores*, Nanotechnology **32**, 395401 (2021).

[73] Z. W. Tan, J.-S. Wang, and C. K. Gan, *First-Principles Study of Heat Transport Properties of Graphene Nanoribbons*, Nano Lett **11**, 214 (2011).

[74] V.-T. Tran, J. Saint-Martin, P. Dollfus, and S. Volz, *Optimizing the Thermoelectric Performance of Graphene Nano-Ribbons without Degrading the Electronic Properties*, Sci Rep **7**, 2313 (2017).

[75] J. W. Jiang, J. S. Wang, and B. Li, *Thermal Conductance of Graphene and Dimerite*, Phys Rev B Condens Matter **79**, 055102 (2009).
[76] C. Li and Z. Tian, *Phonon Transmission across Silicon Grain Boundaries by Atomistic Green's Function Method*, Front Phys **7**, (2019).

[77] J. S. Wang, B. K. Agarwalla, and H. Li, *Transient Behavior of Full Counting Statistics in Thermal Transport*, Phys Rev B Condens Matter Mater Phys **84**, (2011).

[78] R. Tuovinen, N. Säkkinen, D. Karlsson, G. Stefanucci, and R. Van Leeuwen, *Phononic Heat Transport in the Transient Regime: An Analytic Solution*, Phys Rev B **93**, (2016).

[79] X. Chen, J. Yuan, and J. Wang, *Towards Ultrafast Cooling through Transient Phonon Currents: A Closed-Form Solution*, Phys Rev B **103**, (2021).

[80] L. Medrano Sandonas, A. Croy, R. Gutierrez, and G. Cuniberti, *Atomistic Framework for Time-Dependent Thermal Transport*, Journal of Physical Chemistry C **122**, (2018).

[81] Y.-J. Zeng, Z.-K. Ding, H. Pan, Y.-X. Feng, and K.-Q. Chen, *Nonequilibrium Green's Function Method for Phonon Heat Transport in Quantum System*, Journal of Physics: Condensed Matter **34**, 223001 (2022).

[82] M. Ridley, N. W. Talarico, D. Karlsson, N. Lo Gullo, and R. Tuovinen, *A Many-Body Approach to Transport in Quantum Systems: From the Transient Regime to the Stationary State*, J Phys A Math Theor **55**, 273001 (2022).