Title: Hierarchy of Exchange-Correlation Functionals in Computing Lattice Thermal Conductivities of Rocksalt and Zincblende Semiconductors

Abstract: Lattice thermal conductivity ($\kappa_{\rm L}$) is a crucial characteristic of crystalline solids with significant implications for practical applications. While the higher order of anharmonicity of phonon gas model is commonly used for explaining extraordinary heat transfer behaviors in crystals, the impact of exchange-correlation (XC) functionals in DFT on describing anharmonicity has been largely overlooked. Most XC functionals in solids focus on ground state properties that mainly involve the harmonic approximation, neglecting temperature effects, and their reliability in studying anharmonic properties remains insufficiently explored. In this study, we systematically investigate the room-temperature $\kappa_{\rm L}$ of 16 binary compounds with rocksalt and zincblende structures using 8 XC functionals such as LDA, PBE, PBEsol, optB86b, revTPSS, SCAN, rSCAN, r$^2$SCAN in combination with three perturbation orders, including phonon within harmonic approximation (HA) plus three-phonon scattering (HA+3ph), phonon calculated using self-consistent phonon theory (SCPH) plus three-phonon scattering (SCPH+3ph), and SCPH phonon plus three- and four-phonon scattering (SCPH+3,4ph). Our results show that the XC functional exhibits strong entanglement with perturbation order and the mean relative absolute error (MRAE) of the computed $\kappa_{\rm L}$ is strongly influenced by both the XC functional and perturbation order, leading to error cancellation or amplification. The minimal (maximal) MRAE is achieved with revTPSS (rSCAN) at the HA+3ph level, SCAN (r$^2$SCAN) at the SCPH+3ph level, and PBEsol (rSCAN) at the SCPH+3,4ph level. Among these functionals, PBEsol exhibits the highest accuracy at the highest perturbation order. The SCAN-related functionals demonstrate moderate accuracy but are suffer from numerical instability and high computational costs.