Introduction to 2D Materials and Thermal Conductivity

Two-dimensional (2D) materials, such as graphene, hexagonal boron nitride (h-BN), and transition metal dichalcogenides (TMDs) like molybdenum disulfide (MoS₂), have attracted intense research interest due to their exceptional electronic, optical, and mechanical properties. Their atomic-scale thickness — often just one or a few atoms — leads to quantum confinement effects that give rise to novel physical phenomena. Among these, thermal conductivity plays a critical role in determining the performance and reliability of devices built from these materials. For example, graphene exhibits an extremely high intrinsic thermal conductivity (around 3000–5000 W/m·K near room temperature), making it a promising candidate for heat spreading in microelectronics. Conversely, materials like MoS₂ have much lower thermal conductivity, which can be advantageous for thermoelectric energy conversion applications where a low thermal conductivity is desirable. Understanding and predicting how heat propagates through these atomically thin layers is essential for designing next‐generation electronic, photonic, and energy devices. However, direct experimental measurement of thermal transport in 2D materials is challenging due to their small dimensions and the influence of substrate, interfaces, and sample quality. Computational simulation, particularly using Density Functional Theory (DFT), provides a powerful and complementary approach to explore the underlying phonon‐mediated heat transport mechanisms with atomic resolution and without the limitations of sample preparation.

The Role of Density Functional Theory

Foundations of DFT

Density Functional Theory is a quantum mechanical modeling method that maps the many‐electron problem onto a system of non‐interacting electrons moving in an effective potential. The core of DFT is the Kohn–Sham approach, which expresses the total energy of a system as a functional of the electron density. By solving the Kohn–Sham equations self‐consistently, researchers obtain the electronic ground state, including the electron density, total energy, and forces on atoms. Modern DFT implementations employ exchange–correlation functionals such as the local density approximation (LDA) or the generalized gradient approximation (GGA, e.g., PBE or PBEsol). For 2D materials, van der Waals corrections (e.g., DFT‑D3 or vdW‑DF) are often necessary to accurately describe interlayer interactions when modeling few‐layer systems. DFT codes widely used for this purpose include Quantum ESPRESSO, VASP, and CP2K. While DFT itself does not directly compute thermal conductivity, it provides the essential input — namely, the interatomic force constants (IFCs) — that determine lattice dynamics and phonon transport.

Phonon Calculations with DFT

The temperature‐dependent thermal conductivity of crystalline materials is dominated by lattice vibrations, or phonons. To compute phonon properties from first principles, one must calculate the force constants of the crystal. Two main methods are used: the finite‐displacement (frozen phonon) method and density functional perturbation theory (DFPT). In the finite‐displacement approach, small atomic displacements are introduced in a supercell, and the resulting forces are computed from DFT to extract harmonic and anharmonic IFCs. DFPT, on the other hand, directly computes the dynamical matrix by linear response theory, avoiding the need for large supercells for long‐wavelength phonons. Both methods require careful convergence with respect to k‑point sampling, plane‐wave cutoff, and supercell size. The phonon dispersion relations are obtained by diagonalizing the dynamical matrix, and the group velocities (v_g) and phonon lifetimes (τ) can be derived from the harmonic and anharmonic force constants. For thermal conductivity calculations, the third‐order (cubic) anharmonic IFCs are particularly important because they determine the phonon–phonon scattering rates (Umklapp and normal processes). Tools like Phonopy (for harmonic calculations) and Phono3py (for anharmonic calculations) integrate seamlessly with DFT outputs and automate much of the workflow.

Methodology for Thermal Conductivity Simulation

Structure Optimization

The first step in any DFT‐based thermal conductivity study is to obtain the fully relaxed atomic structure of the 2D material. For monolayers, a vacuum region of at least 15 Å is added perpendicular to the plane to eliminate spurious interactions between periodic images. The in‐plane lattice parameters and atomic positions are optimized until residual forces are below a threshold (typically 10⁻⁴ eV/Å) and stresses are negligible. The choice of exchange–correlation functional and van der Waals correction can significantly affect the relaxed lattice constants and, consequently, the phonon frequencies. For example, the PBE functional tends to overestimate lattice constants of graphene slightly, while PBEsol yields better agreement with experiment for some TMDs. A careful convergence test on the energy cutoff and k‑point grid (e.g., 24×24×1 for a unit cell of graphene) is essential to ensure structural stability and accurate forces.

Phonon Dispersion and Force Constants

Once the equilibrium structure is determined, harmonic phonon calculations are performed. Using DFPT, the dynamical matrix is computed on a coarse q‑point grid (e.g., 6×6×1) and then Fourier interpolated onto a finer mesh for phonon dispersion plots. Acoustic sum rules must be enforced to eliminate imaginary modes (negative frequencies) that would indicate structural instability. For 2D materials, the ZA (out‐of‐plane acoustic) branch often has a quadratic dispersion near the Γ point, a signature of the reduced dimensionality. For anharmonic calculations, a supercell (e.g., 4×4×1 for graphene) is constructed, and a set of displaced configurations is generated to compute third‐order IFCs via finite differences. The cutoff radius for anharmonic interactions is typically chosen to include all neighbors within a few angstroms; an insufficient cutoff can lead to inaccurate scattering rates, especially for materials with long‐range interactions like graphene.

Solving the Boltzmann Transport Equation

The lattice thermal conductivity κ_lat is obtained by solving the Boltzmann Transport Equation (BTE) for phonons. In the relaxation time approximation (RTA), κ = (1/N_q V) Σ_{q,ν} C_{q,ν} v_{q,ν}² τ_{q,ν}, where C is the mode‐specific heat capacity, v the group velocity, and τ the phonon lifetime from all scattering mechanisms (anharmonic, isotope, boundary). However, the RTA neglects the correction from the phonon distribution function’s deviation from equilibrium, which can be significant for materials with strong normal scattering (e.g., graphene). An iterative solution of the BTE, as implemented in codes like ShengBTE (which uses the full phonon dispersion and third‐order IFCs as input), provides more accurate results. For graphene at room temperature, iterative BTE calculations yield thermal conductivity values in the range 3000–4000 W/m·K, with the dominant contribution from long‐mean‐free‐path phonons in the ZA branch. The convergence of κ with respect to the q‑point grid (e.g., 100×100×1 or finer) is critical because the thermal conductivity integrals involve an integration over the entire Brillouin zone.

Practical Considerations

Several numerical parameters affect the accuracy of computed thermal conductivity:

  • Supercell size: For harmonic calculations, a sufficiently large supercell (or a dense DFPT q‑mesh) is needed to avoid wrapping errors. For anharmonic force constants, the supercell must be larger than the spatial extent of the interatomic interactions considered.
  • k‑point sampling: A dense k‑point mesh (e.g., 12×12×1 for a 2D supercell) ensures converged forces and total energies.
  • Smearing scheme: Gaussian or Methfessel–Paxton smearing with a width of 0.01–0.02 Ry is often used to handle metallic or semimetallic systems like graphene. Insulators are less sensitive.
  • Isotope scattering: Natural isotopic mixtures (e.g., 98.9% ¹²C and 1.1% ¹³C in carbon) reduce thermal conductivity by 10–20% in graphene. This contribution is routinely included via the isotope scattering rate formula (e.g., Tamura’s model).
  • Temperature dependence: The phonon lifetimes and heat capacities are temperature‐dependent, so κ_lat is computed at multiple temperatures (e.g., 100 K, 300 K, 500 K) to obtain the full thermal conductivity curve.

Challenges in DFT-Based Thermal Conductivity Simulations

Despite their success, DFT‐based simulations of thermal transport in 2D materials face several challenges. Computational cost grows rapidly with the number of atoms in the supercell and with the density of the q‑point grid required for converged BTE solutions. For many 2D materials with large unit cells (e.g., MoS₂ with 6 atoms per monolayer), full anharmonic calculations become expensive. Additionally, the accuracy of DFT depends on the choice of exchange–correlation functional; thermal conductivity values can differ by 20–30% between LDA and GGA due to changes in phonon frequencies and anharmonicity. Size effects are another concern: in experimental samples, sample length (typically micrometers) is much smaller than the mean free path of long‐wavelength phonons in high‐conductivity materials, leading to ballistic or quasi‐ballistic transport that the BTE with periodic boundary conditions cannot capture without special treatments (e.g., inclusion of boundary scattering). Furthermore, higher‐order phonon scattering (fourth‐order and beyond) is usually neglected, but it may become important at high temperatures or for materials with strong anharmonicity. The role of phonon–electron scattering in metallic or doped 2D materials is also often omitted, though it can reduce thermal conductivity in systems like doped graphene. Finally, the presence of defects, grain boundaries, and strain in realistic samples is challenging to model without introducing significant approximations. DFT supercell calculations with a single vacancy or a substitutional defect can provide insights, but the spatial and statistical complexity of real disorder remains a difficulty.

Applications and Future Directions

Thermal Management in Electronics

The extremely high thermal conductivity of graphene makes it a natural candidate for heat spreaders in high‐power electronics. DFT simulations help optimize the use of graphene by predicting how substrate coupling, number of layers, and grain size affect heat dissipation. For instance, simulations show that suspending graphene eliminates the substrate‐induced phonon scattering, restoring its intrinsic conductivity, whereas placing it on SiO₂ reduces the conductivity to 500–1000 W/m·K due to interface scattering. Similar studies on hexagonal boron nitride and transition metal dichalcogenides guide the design of all‐2D heterostructures where thermal management is critical.

Thermoelectric Materials

Two‐dimensional materials with low thermal conductivity are attractive for thermoelectric energy conversion, which requires a high Seebeck coefficient, high electrical conductivity, and low thermal conductivity. DFT‐based thermal conductivity screening has identified promising 2D thermoelectrics such as SnSe (which has an extremely low κ_lat of ~0.4 W/m·K in the out‐of‐plane direction), Bi₂Te₃, and few‐layer WSe₂. By computing the phonon dispersion and scattering rates, researchers can pinpoint the mechanisms leading to low κ (e.g., soft acoustic modes, rattling atoms, or complex unit cells) and propose chemical substitutions to further reduce thermal conductivity without degrading electronic transport. The predictions guide experimental synthesis, accelerating the search for high‐efficiency thermoelectric materials.

Machine Learning and High-Throughput Screening

The high computational cost of DFT‐based BTE calculations has motivated the development of machine‐learning (ML) interatomic potentials that can faithfully reproduce DFT forces while being orders of magnitude faster. Models such as the Behler–Parrinello neural network, the Gaussian approximation potential (GAP), and the deep potential (DP) have been trained on DFT data for 2D materials and used to compute phonon transport properties with near‐DFT accuracy. These ML potentials enable thermal conductivity calculations for larger supercells (e.g., 10 000 atoms) and longer time scales, allowing the study of grain boundaries, nanotubes, and disordered alloys. Coupled with high‐throughput frameworks like the Materials Project and AFLOW, researchers are screening thousands of hypothetical 2D compounds for extreme thermal transport properties, identifying candidates such as borophene, phosphorene, and silicene. The combination of DFT, ML, and automation promises to accelerate the discovery of 2D materials with tailored thermal properties for a wide range of applications.

Conclusion

Simulating thermal conductivity in 2D materials with Density Functional Theory provides a detailed, atomistic understanding of heat transport mechanisms that underpin the performance of nanoscale devices. By combining DFT‐derived force constants with the Boltzmann Transport Equation, researchers can predict intrinsic thermal conductivity values, explore the influence of defects and strain, and guide the engineering of materials for thermal management or thermoelectric conversion. Despite computational challenges, advances in algorithms, machine learning, and high‐performance computing are steadily improving the accuracy and accessibility of these simulations. As the library of known 2D materials continues to expand, DFT‐based thermal conductivity simulations will remain an indispensable tool for both scientific discovery and technological innovation. The insights gained from these simulations are already informing the design of next‐generation electronics, energy converters, and flexible sensors, turning atomic‐scale understanding into real‐world applications.