The second-order IFCs describe pairwise harmonic interactions and thus define phonons. These can be easily obtained using the now widely known methods such as density-functional perturbation theory (DFPT) and the finite-displacement method at minimal to moderate computational expense. Macroscopic thermal properties such as the vibrational free energy, entropy, and heat capacity at the harmonic level can then be derived directly from the theory of phonons. Anharmonic IFCs of third-and-higher order serve as perturbation to the harmonic phonons. They are directly related to other critical thermal properties such as thermal expansion, thermal conductivity, and dynamical phase stabilization and provide corrections to the harmonic free energy. In contrast to the second-order IFCs, the direct calculation of anharmonic IFCs has long been formidable. Usually, DFPT is applied only for second order and could be done for third order , but is essentially infeasible beyond this due to the 2n+1 theorem of the linear response theory that requires beyond second-order derivatives of electronic wave functions . As for finite-displacement, the combinatorial explosion of tensor elements in Φ d (d ≥ 3) with increasing order of IFCs and growing numerical errors severely hampers its application to the anharmonic terms. The finite-displacement approach has been taken for thirdorder IFCs , but this is inadequate for fourth-order IFC. Utilizing polynomial machine learning potentials can reduce the computational costs by one order of magnitude compared with the conventional finite-displacement approach for lattice thermal conductivity calculations, as demonstrated recently by Togo et al. A high efficiency even at the third-order level is the preferred guarantee for a truly wide-scale deployment to small and large systems.