⇦Chinese Physics Letters, 2023, Vol. 40, No. 9, Article code 090502 Inverse Design of Phononic Crystal with Desired Transmission via a Gradient-Descent Approach Yuhang Wei (魏宇航) and Dahai He (贺达海)* Affiliations Department of Physics and Jiujiang Research Institute, Xiamen University, Xiamen 361005, China Received 22 July 2023; accepted manuscript online 24 August 2023; published online 10 September 2023 *Corresponding author. Email: dhe@xmu.edu.cn Citation Text: Wei Y H and He D H 2023 Chin. Phys. Lett. 40 090502    Abstract We propose a general approach based on the gradient descent method to study the inverse problem, making it possible to reversely engineer the microscopic configurations of materials that exhibit desired macroscopic properties. Particularly, we demonstrate its application by identifying the microscopic configurations within any given frequency range to achieve transparent phonon transport through one-dimensional harmonic lattices. Furthermore, we obtain the phonon transmission in terms of normal modes and find that the key to achieving phonon transparency or phonon blocking state lies in the ratio of the mode amplitudes at ends.
cpl-40-9-090502-fig1.png
cpl-40-9-090502-fig2.png
cpl-40-9-090502-fig3.png
 DOI:10.1088/0256-307X/40/9/090502 © 2023 Chinese Physics Society Article Text The advancement of materials science has exerted a significant influence on human society. With the rise of first-principle theory and recent strides in micro-nanotechnology, computational prediction of material properties for specific structures has become feasible, reducing reliance solely on experiments. Currently, various machine learning algorithms can be employed in materials research to replace traditional methods, which represents a new paradigm to predict and design new metamaterials.[1-3] A more efficient approach to developing new materials is based on inverse design using known material properties.[4] This inverse approach not only enables the design of materials with desired performance but also provides valuable insights into how microstructures can cooperate to achieve macroscopic properties. Inverse design of materials generally involves two approaches: optimization in parameter space[5] and obtaining a generative model through data-driven methods.[6-9] As a scenario of its application, thermal metamaterials have experienced rapid development in recent years, whose main principle is to realize macroscopic thermal management by adjusting the thermal transport of microstructures.[10-14] Thermal devices made of thermal metamaterial such as thermal diode,[15] thermal transistor,[16] thermal cloak[17,18] and phonon focusing[19] have been theoretically proposed. With the advancement of various machine learning and intelligent algorithms, researchers can more conveniently use computer-aided methods to predict the thermal transport properties of materials[20] or develop new thermal metamaterials.[21,22] For example, multi-objective optimization,[23] active machine learning,[24] genetic algorithm optimization,[25] and particle swarm optimization[26] have been applied in this field. However, machine learning often necessitates a substantial training dataset, resulting in time-consuming learning processes. The development of efficient approaches for the inverse design of thermal metamaterials and the clarification of underlying physical mechanisms remains an intriguing challenge. In this Letter, we study heat conduction in one-dimensional harmonic lattices connected to Langevin heat baths at both ends. The phonon transmission can be calculated using Langevin Green's function theory.[27] However, the problem of how to inverse design the corresponding configuration based on given transmission remains to be explored. A remarkable approach to designing materials lies in exploring the properties of unit cells and then combining them to construct the entire system.[28] Conversely, we propose an inverse design approach based on gradient descent, which allows for overall optimization starting from random configurations, leading to rapid convergence to desired configurations. When the transmission of a phonon mode reaches $100\%$ or $0$, it is called a phonon transparent or blocked state. We successfully apply the gradient descent method to find configurations for phonon transparent states in one-dimensional harmonic lattices. By deriving the normal mode form of the transmission, we discover that in the case of weak coupling between the system and heat baths, each normal mode corresponds to a phonon transport channel, and the opening and closing of the channels are determined by the amplitudes of the particles connected to the heat bath in that mode. We explain the origins of phonon transparent state and blocking state, provide theoretical guidance for designing super thermal conductivity metamaterials[28-30] and acoustic frequency combs, and depict new insights into the physical mechanisms of phonon transport. Model and Method. Consider a general one-dimensional harmonic system with the Hamiltonian \begin{align} H=\frac{1}{2}P^{\scriptscriptstyle{\rm T}}M^{-1}P+\frac{1}{2}X^{\scriptscriptstyle{\rm T}}KX, \tag {1} \end{align} where the column vectors $P$ and $X$ denote the particle momentum and displacement relative to the equilibrium position, respectively. The diagonal matrix $M$ represents the particle mass, and the real symmetric matrix $K$ is the force matrix satisfying $M\ddot{X}+KX=0$. The $1$st and $N$th particles are connected to the Langevin heat baths with the temperature $T_{\scriptscriptstyle{\rm L}}$ and $T_{\scriptscriptstyle{\rm R}}$, and the heat flow between the two baths can be given by the Langevin Green's function theory, \begin{align} J=\frac{k_{\scriptscriptstyle{\rm B}}(T_{\scriptscriptstyle{\rm L}}- T_{\scriptscriptstyle{\rm R}})}{2\pi}\int_{0}^{+\infty}\mathcal{T}(\omega){\rm d}\omega, \tag {2} \end{align} where $\mathcal{T}(\omega)=4\gamma_{\scriptscriptstyle{\rm L}}\gamma_{\scriptscriptstyle{\rm R}}\omega^{2}|G_{\scriptscriptstyle{1N}}|^{2}$ is the transmission of heat current with Green's function $G=Z^{-1}=(K-\omega^{2}M-\varSigma)^{-1}$. The only nonzero elements of $\varSigma$ are $\varSigma_{11}=i\omega \gamma_{\scriptscriptstyle{\rm L}}$ and $\varSigma_{\scriptscriptstyle{NN}}=i\omega \gamma_{\scriptscriptstyle{\rm R}}$, where $\gamma_{\scriptscriptstyle{\rm L,R}}$ represent the dissipation coefficients. The inverse problem we need to solve is to determine the suitable force matrix $K$ and mass matrix $M$ that satisfy the desired phonon transmission. However, directly setting the elements of the force matrix $K$ as adjustable parameters may result in the eigenfrequencies of the system being negative or even complex when solved reversely. Moreover, if the long-range interaction is concerned, the force matrix $K$ cannot well reflect the distances of interactions. To avoid the aforementioned problems, one may set the elastic coefficients of the system as adjustable parameters. The potential energy between the $i$th and $j$th particles is given by $V_{ij}=k_{ij}(x_i-x_j)^2/2$. The on-site potential energy of the $i$th particle is given by $V_{ii}=k_{ii}x_i^2/2$. Furthermore, we can define the elastic coefficient matrix $k$, where the diagonal elements $k_{ii}$ represent the on-site elastic coefficients of the $i$th particle, and the off-diagonal elements $k_{ij}$ represent the elastic coefficients between the $i$th and $j$th particles. The relationship between the elastic coefficient matrix $k$ and the force matrix $K$ is \begin{align} K_{ij}=-k_{ij}~(i\ne j),~~~K_{ii}=\sum_{j=1}^{N}k_{ij}, \tag {3} \end{align} where $k_{ij}\ge 0(i,\,j=1,\,2,\,\dots ,\,N)$. At this point, the inverse problem is transformed into finding the configuration parameters $\{M_0,k_0\}$ that satisfy the given transmission $\mathcal{T}_0(\omega)$. Since the transmission of the system will approach $0$ at sufficiently high frequencies, it is necessary to choose a cutoff frequency of $\omega_{\max}$ to set the target transmission to $0$ for frequencies greater than $\omega_{\max}$. We randomly select an initial set of parameters $\{M_1,k_1\}$, which are denoted as $\{M_j,k_j\}$ in the $j$th iteration. The loss function is defined as \begin{align} \delta (M,k;\omega)=\int_{0}^{\omega _{{\max}}}[\mathcal{T}(M,k;\omega)-\mathcal{T}_0(\omega) ]^2 {\rm d}\omega. \tag {4} \end{align} Then the problem to find configuration parameters that satisfy the target transmission is transformed into finding the zeros of the loss function. We employ the gradient descent method, as illustrated in Fig. 1, to optimize the parameters. During the iteration process, the loss function gradually decreases until it falls below the threshold value $\delta_{{\min}}$ at the $n$th step. At this point, we obtain the optimized results $\{M_n,k_n\}$, which enable the transmission $\mathcal{T}(M_n,k_n;\omega)$ to approximate the target transmission $\mathcal{T}_0(\omega)$.
cpl-40-9-090502-fig1.png
Fig. 1. Flowchart of the gradient descent method. Firstly, all $s$ adjustable parameters in the physical quantity $\{M,k\}$ are uniformly recorded as vector ${\boldsymbol a}$. Then, we continuously move the vector ${\boldsymbol a}$ in the direction of the gradient of the loss function $\delta({\boldsymbol a})$, with the step size controlled by $h$. Finally, after $n$ iterations, when the loss function becomes smaller than the threshold value $\delta_{{\min}}$, we obtain the optimized physical quantities $\{M,k\}$ from the vector ${\boldsymbol a}$.
The most time-consuming part of this algorithm is the calculation of $G_{\scriptscriptstyle{1N}}$ with given lattice configuration parameters $\{M,k\}$. This is because we need to compute the inverse of the $N$-order matrix $Z=K-\omega^2 M-\varSigma$. Traditionally, the transfer matrix method is used to calculate $G_{\scriptscriptstyle{1N}}$. When considering only nearest neighbor particles with harmonic interactions, $Z$ becomes a tridiagonal matrix. Thus, we have \begin{align} G_{\scriptscriptstyle{1N}}=\frac{Z^{*}_{\scriptscriptstyle{1N}}}{{\rm Det}(Z)}=\frac{\prod_{i=1}^{N-1}K_{i,i+1} }{{\rm Det}(Z)}. \tag {5} \end{align} The determinant of the tridiagonal matrix $Z$ can be expressed as \begin{align} {\rm Det}(Z)=(1-\varSigma _{11}) \bigg[\prod_{i=1}^{N}T_i \bigg] \begin{pmatrix} 1\\ \varSigma _{\scriptscriptstyle{NN}} \end{pmatrix}, \tag {6} \end{align} where \begin{align} T_i=\begin{pmatrix} K_{ii}-m_i\omega ^2&-K_{i,i+1} \\ K_{i-1,i}&0 \end{pmatrix}, \tag {7} \end{align} with $K_{0,1}=K_{N,N+1}=1$. When the lattice configuration exhibits periodicity, the conventional transfer matrix method can reduce the time complexity from $O(N)$ to $O(\log N)$ using the fast power method. However, when the lattice configuration lacks periodicity, the disadvantage of the traditional method lies in its difficulty to perform parallel computations for different independent variable frequencies. Here, we propose an iterative method that vectorizes and parallelizes the independent variable frequencies for computation. Since only $G_{\scriptscriptstyle{1N}}$ is desirable, solving the inverse of matrix $Z$ is unnecessary. Instead, one can only focus on solving the column where $G_{\scriptscriptstyle{1N}}$ resides. Following this idea, we can simplify $ZG=I$ as \begin{align} \begin{pmatrix} Z_{11} & Z_{12} & & &\\ Z_{21} & Z_{22} & Z_{23} & &\\ & Z_{32} & Z_{33} & &\\ & & & \ddots &Z_{\scriptscriptstyle{N-1,N}}\\ & & & Z_{N,N-1}& Z_{\scriptscriptstyle{NN}} \end{pmatrix} \begin{pmatrix} G_{\scriptscriptstyle{1N}} \\G_{2\,N} \\G_{3\,N} \\\vdots \\G_{\scriptscriptstyle{NN}} \end{pmatrix} =\begin{pmatrix} 0 \\0 \\0 \\0 \\1 \end{pmatrix}. \tag {8} \end{align} We can solve $G_{\scriptscriptstyle{1N}}$ iteratively using the set of $N$ linear equations in Eq. (8). By defining $G_{\rm iN}=f_i\cdot G_{\scriptscriptstyle{1N}}$, we have $f_1 =1$, $f_2=-Z_{11}/Z_{12}$ and \begin{align} f_i=-\frac{Z_{i-2,i-2}\cdot f_{i-2}+Z_{i-1,i-1}\cdot f_{i-1}}{Z_{i-1,i}}, \tag {9} \end{align} with $i=3,\,4,\,\dots ,\,N$. $Z_ {i-1,i}=-k_{i-1,i}\ne 0$ ensures the rationality of Eq. (9). After finding $f_{\scriptscriptstyle{N-1}}$ and $f_{\scriptscriptstyle{N}}$ through iterative algorithm, $G_{\scriptscriptstyle{1N}}$ can be represented as \begin{align} G_{\scriptscriptstyle{1N}}=\frac{1}{Z_{\scriptscriptstyle{N-1,N}}\cdot f_{\scriptscriptstyle{N-1}}+Z_{\scriptscriptstyle{NN}}\cdot f_{\scriptscriptstyle{N}}}. \tag {10} \end{align} Compared to the traditional transfer matrix method, the advantage of the linear equation method is using a more efficient vectorized parallel algorithm. Additionally, the linear equation method can be extended to disordered lattices with next-nearest-neighbor interactions, while maintaining a time complexity of $O(N)$, which is the transfer matrix method we cannot achieve. In contrast, the transfer matrix method cannot be applied to non-nearest neighbor coupled lattices. Result and Discussion. Using the gradient descent method, one can easily design inversely a set of configuration parameters $\{M,k\}$ that satisfy a given target transmission $\mathcal{T}_0(\omega)$. Figure 2 presents an example of finding a configuration with the property of transparent phonon transmission. Without losing generality, we set the target transmission as \begin{align} \mathcal{T}_0(\omega)=\left\{\begin{matrix} 1,~~~~0\le \omega \le \omega_{\rm b},\\ 0,~~~~~~~~~\omega > \omega_{\rm b}, \end{matrix}\right. \tag {11} \end{align} where the upbound of the frequency range $\omega_{{\rm b}}$ is set to $1.5$. In this case we have $\langle m\rangle, \langle k \rangle \approx 1$, which is an unnecessary condition and only for convenience. We set the adjustable parameters as the masses of all $N$ particles and the elastic coefficients of the $N-1$ springs in a harmonic chain consisting of $N$ particles, considering only nearest-neighbor interactions and using free boundary conditions.
cpl-40-9-090502-fig2.png
Fig. 2. (a) Comparison of the target transmission with the transmission corresponding to the optimized configuration shown in (b). (b) One of the optimized configuration parameters. The blue line represents the masses of particles and the red line represents the elastic coefficients of springs. Here we consider the case of nearest-neighbor interactions and set $N=100$, $\gamma_{\scriptscriptstyle{\rm L,R}}=1$.
In fact, the exact solution $\{M_0,k_0\}$ of the inverse problem is continuously distributed on a hypersurface in the parameter space. The more adjustable parameters there are, or the smaller the $\delta_{{\min}}$ is, the more accurate the numerical solution will be. Figure 2(b) shows one of the approximate solutions for $\{M_0,k_0\}$. It is worth noting that the regular transmission shown in Fig. 2(a) corresponds to the irregular (non-periodic) configuration shown in Fig. 2(b). The reason why such configurations can achieve quasi-transparent phonon transport within a certain frequency range lies in the synergistic effect among the particles. In order to explain the reason behind the quasi-transparent phonon transport, we further derive the expression of the transmission in terms of normal modes as (see the Supplemental Material for details) \begin{align} \mathcal{T}(\omega)=\frac{4\omega ^2 S_1^2}{ (1-\omega ^2 S_2)^2 + \omega ^2 S_3^2}, \tag {12} \end{align} where \begin{align} &S_1=\sum_{i=1}^{N} \frac{E_1^i E_{\scriptscriptstyle{N}}^i}{\omega _i^2-\omega ^2},\notag\\ &S_2=\sum_{1=i < j}^{N}\frac{(E_1^i E_{\scriptscriptstyle{N}}^j - E_1^j E_{\scriptscriptstyle{N}}^i)^2}{(\omega _i^2-\omega ^2)(\omega _j^2-\omega ^2)},\notag\\ &S_3=\sum_{i=1}^{N}\frac{ (E_1^i)^2+(E_{\scriptscriptstyle{N}}^i)^2 }{\omega _i^2-\omega ^2}. \tag {13} \end{align} The reduced amplitude $E_{1,N}^i=e_{1,N}^i\sqrt{ \gamma_{\scriptscriptstyle{\rm L,R}}/m_{{\scriptstyle 1,N}}}$ represents the amplitude $e_{{\scriptstyle 1,N}}^i$ of $i$th mode incorporating the effect of $\gamma_{\scriptscriptstyle{\rm L,R}}$ and $m_{{1,N}}$. Under the condition of weak system-bath coupling and absence of eigenfrequency degeneracy, we have \begin{align} \lim_{\gamma_{\scriptscriptstyle{\rm L}}\gamma_{\scriptscriptstyle{\rm R}}\to 0}\mathcal{T}(\omega)=\sum_{i=1}^{N}\mathcal{T}_i(\omega), \tag {14} \end{align} where \begin{align} \mathcal{T}_i(\omega)=\frac{4(E_1^iE_{\scriptscriptstyle{N}}^i)^2}{d\cdot(\omega_i-\omega)^2+[(E_1^i)^2+(E_{\scriptscriptstyle{N}}^i)^2]^2}, \tag {15} \end{align} with $\omega_i=0$, $d=1$ or $\omega_i\ne 0$, $d=4$. From Eq. (15), one can observe that the contributions of each normal mode to the transmission exhibit independent Lorentzian line shapes. Then the height $\mathcal{T}_i(\omega_i)$ and half-height width $\Delta\omega_i$ of the Lorentzian line shape shown in Eq. (15) can be expressed as \begin{align} &\mathcal{T}_i(\omega_i)= \frac{4}{(E_1^i/E_{\rm N}^i)^2+(E_{\rm N}^i/E_1^i)^2+2},\notag\\ &\Delta \omega_i= (E_1^i)^2+(E_{\rm N}^i)^2. \tag {16} \end{align} We stress that $\Delta \omega_i$ in Eq. (16) comes from the dissipation, which is different from the mean free path coming from the anharmonicity.[31] It is worth noting that Eq. (16) is irrelevant with $d$ in Eq. (15). It comes from the fact that the half-height width of the Lorentzian line shape corresponding to $\omega_i=0$ is twice the value of the other modes. We only take half of the width since the frequency ranges from $0$ to positive infinity. Interestingly, if defining the left-end and right-end thermal conductances of the $i$th mode as $g_{_{\scriptstyle \rm L,R}}^i=k_{\scriptscriptstyle{\rm B}}(E_{1,N}^{i})^2$, one can observe that thermal conductance and electrical conductance exhibit similar series connection properties shown by \begin{align} \frac{1}{g^i}=\frac{1}{g_{\scriptscriptstyle{\rm L}}^i}+\frac{1}{g_{\scriptscriptstyle{\rm R}}^i}, \tag {17} \end{align} where the total thermal conductance of the $i$th mode in the harmonic lattice is represented by \begin{align} g^i=\frac{J_i}{T_{\scriptscriptstyle{\rm L}}-T_{\scriptscriptstyle{\rm R}}}=k_{\scriptscriptstyle{\rm B}}\int_{0}^{+\infty} \mathcal{T}_i(\omega){\rm d}\omega=\frac{k_{\scriptscriptstyle{\rm B}}}{2}\mathcal{T}_i(\omega_i)\Delta\omega_i. \tag {18} \end{align} It can be observed that the contribution of the $i$th mode to the heat flux $J_i$ in Eq. (18) agrees with previous studies.[32] Equation (15) gives the contribution of each mode to the transmission in the case of weak coupling. According to Eq. (15), when $|E_1^i/E_{\rm N}^i| =1$, the height of the Lorentzian line shape of the $i$th mode becomes $1$, representing transparent phonon transport. On the other hand, when there is a significant difference between $E_1^i$ and $E_{\scriptscriptstyle{N}}^i$, the height of the Lorentzian line shape vanishes, indicating that the phonon mode is blocked. The insulation of phonon transport is peculiarly different from the phonon Anderson localization occurring in disordered lattices.[33,34] For a localized mode, as long as the reduced amplitude ratio at the ends is close to $1$, transparent phonon transport can still be achieved, corresponding to a Lorentzian line shape with a height close to $1$ and a vanishing half-height width.
cpl-40-9-090502-fig3.png
Fig. 3. (a) Comparison of the target transmission (black dashed line) with the transmission corresponding to the optimized configuration shown in (b) (red solid line). (b) One of the optimized configuration parameters. (c) The ratio of the reduced mode amplitudes at the ends $|E_1^i/E_{\scriptscriptstyle{N}}^i|$ of the $i$th mode as a function of frequency $\omega$. One can observe that $|E_1^i/E_{\scriptscriptstyle{N}}^i|=1$ corresponds to the transparent phonon transport in (a). (d) The eigenvectors of the $20$th and $80$th modes in the blocked state are shown. The $20$th mode with a lower eigenfrequency does not exhibit Anderson localization, while the $80$th mode with a higher eigenfrequency undergoes Anderson localization. Here we consider the case of nearest-neighbor interactions and set $N=100$, $\gamma_{\scriptscriptstyle{\rm L,R}}=1$.
To verify the above conclusions, we design a new target transmission without loss of generality as \begin{align} \mathcal{T}_0(\omega)=\left\{\begin{matrix} \!1,~~~~\omega\in [0,\omega_{{\rm b}_1}]\cup [\omega_{{\rm b}_2},\omega_{{\rm b}_3}],\\ ~0,~~~~{\rm otherwise},~~~~~~~~~~~~~~~~~~~~~ \end{matrix}\right. \tag {19} \end{align} where $\omega_{{\rm b}_1}=0.5$, $\omega_{{\rm b}_2}=1$, and $\omega_{{\rm b}_3}=1.5$. The adjustable parameters are set as the masses of all $N$ particles in a harmonic chain, with the elastic coefficients of the $N-1$ springs all set to $1$, considering only nearest-neighbor interactions under free boundary conditions. As shown in Figs. 3(a) and 3(b), we have found a set of configuration parameters that can achieve the target transmission. It is worth noting that our findings challenge the traditional notion that transparent phonon transport requires a coordinated combination of mass and elastic coefficients.[28] Comparing Figs. 3(a) and 3(c), we can find that, when the reduced amplitude ratio $|E_1^i/E_{\scriptscriptstyle{N}}^i|$ is close to 1 at both ends, the normal mode corresponds to the transparent phonon state, while its deviation from $1$ corresponds to the phonon blocking state. Figure 3(d) illustrates the eigenvectors of two different phonon blocking states at low and high frequencies, indicating that the cause of phonon blocking state is not necessarily Anderson localization. It is worth noting that Fig. 3 is given far from the limit $\gamma_{\scriptscriptstyle{\rm L}}\gamma_{\scriptscriptstyle{\rm R}}\to 0$ in Eq. (14). However, the value of $|E_1^i/E_{\scriptscriptstyle{N}}^i|$ can still provide insights into whether the transmission around the eigenfrequency of the $i$th mode approaches $1$ or $0$. In summary, we have proposed a approach based on the gradient descent method to solve the inverse problem of designing the configuration of harmonic lattices that matches the target transmission. Compared to traditional methods, our approach enables quick exploration of various solutions via global optimization in addition to improving the calculation steps for transmission. It should be noted that the efficiency of optimization is affected by the setting of random initial values. Based on the characteristics of the target function, we suggest setting the initial value uniformly and adding a certain level of noise to it. Using this method, we have successfully constructed the harmonic lattice with transparent phonon transport. Furthermore, we obtain a form of transmission with respect to normal modes and find that the contribution of each mode to the transmission exhibits a Lorentzian line shape, with its height determined by the ratio of reduced mode amplitudes at the ends. Our findings provide a novel approach to blocking phonon transport which is distinct from the Anderson localization, and explain the mechanism behind the phonon transparency. Overall, our study offers new insights into phonon transport and designing novel materials for thermal and acoustic applications. Acknowledgments. This work was supported by the National Natural Science Foundation of China (Grant No. 12075199), the Natural Science Foundation of Fujian Province (Grant No. 2021J01006) and Jiangxi Province (Grant No. 20212BAB201024). References Machine learning in materials informatics: recent applications and prospectsMachine learning for molecular and materials scienceDeep Generative Models for Materials Discovery and Machine Learning-Accelerated InnovationMachine-enabled inverse design of inorganic solid materials: promises and challengesSearch for Catalysts by Inverse Design: Artificial Intelligence, Mountain Climbers, and AlchemistsInverse design with deep generative models: next step in materials discoveryGenerative Deep Neural Networks for Inverse Materials Design Using Backpropagation and Active LearningInverse Design of Solid-State Materials via a Continuous RepresentationMachine‐Learning Microstructure for Inverse Material DesignFull Control and Manipulation of Heat Signatures: Cloaking, Camouflage and Thermal MetamaterialsThermal Metamaterial: Fundamental, Application, and OutlookControlling macroscopic heat transfer with thermal metamaterials: Theory, experiment and applicationTransforming heat transfer with thermal metamaterials and devicesEmerging theory and phenomena in thermal conduction: A selective reviewThermal Diode: Rectification of Heat FluxNegative differential thermal resistance and thermal transistorExperimental Demonstration of a Bilayer Thermal CloakMetamaterials for Manipulating Thermal Radiation: Transparency, Cloak, and ExpanderPhonon Focusing Effect in an Atomic Level Triangular StructureMachine learning approach for the prediction and optimization of thermal transport propertiesMachine-Learning-Optimized Aperiodic Superlattice Minimizes Coherent Phonon Heat ConductionMachine learning-optimized Tamm emitter for high-performance thermophotovoltaic system with detailed balance analysisDirectional Design of Materials Based on Multi-Objective Optimization: A Case Study of Two-Dimensional Thermoelectric SnSeUnexpected thermal conductivity enhancement in aperiodic superlattices discovered using active machine learningDesigning thermal demultiplexer: Splitting phonons by negative mass and genetic algorithm optimization*Designing radiative cooling metamaterials for passive thermal management by particle swarm optimizationHeat transport in low-dimensional systemsAnomalous transparency induced by cooperative disorders in phonon transportTotal-transmission and total-reflection of individual phonons in phononic crystal nanostructuresEnhancing thermal transport in multilayer structures: A molecular dynamics study on Lennard-Jones solidsApproach to Phonon Relaxation Time and Mean Free Path in Nonlinear LatticesThermal conduction in classical low-dimensional latticesAbsence of Diffusion in Certain Random LatticesLocalization of Eigenstates and Transport Phenomena in the One-Dimensional Disordered System
[1] Ramprasad R, Batra R, Pilania G, Mannodi-Kanakkithodi A, and Kim C 2017 npj Comput. Mater. 3 54
[2] Butler K T, Davies D W, Cartwright H, Isayev O, and Walsh A 2018 Nature 559 547
[3] Fuhr A S and Sumpter B G 2022 Front. Mater. 9 865270
[4] Noh J, Gu G H, Kim S, and Jung Y 2020 Chem. Sci. 11 4871
[5] Freeze J G, Kelly H R, and Batista V S 2019 Chem. Rev. 119 6595
[6] Lu S H, Zhou Q H, Chen X Y, Song Z L, and Wang J L 2022 Natl. Sci. Rev. 9 nwac111
[7] Chen C T and Gu G X 2020 Adv. Sci. 7 1902607
[8] Noh J, Kim J, Stein H S, Sanchez-Lengeling B, Gregoire J M, Aspuru-Guzik A, and Jung Y 2019 Matter 1 1370
[9] Pei Z R, Rozman K A, Doğan O N, Wen Y H, Gao N, Holm E A, Hawk J A, Alman D E, and Gao M C 2021 Adv. Sci. 8 2101207
[10] Han T C, Bai X, Thong J T L, Li B W, and Qiu C W 2014 Adv. Mater. 26 1731
[11] Wang J, Dai G L, and Huang J P 2020 iScience 23 101637
[12] Yang S, Wang J, Dai G L, Yang F B, and Huang J P 2021 Phys. Rep. 908 1
[13] Li Y, Li W, Han T C, Zheng X, Li J X, Li B W, Fan S H, and Qiu C W 2021 Nat. Rev. Mater. 6 488
[14] Chen J, He J, Pan D K, Wang X T, Yang N, Zhu J J, Yang S A, and Zhang G 2022 Sci. Chin. Phys. Mech. & Astron. 65 117002
[15] Li B W, Wang L, and Casati G 2004 Phys. Rev. Lett. 93 184301
[16] Li B W, Wang L, and Casati G 2006 Appl. Phys. Lett. 88 143501
[17] Han T C, Bai X, Gao D L, Thong J T L, Li B W, and Qiu C W 2014 Phys. Rev. Lett. 112 054302
[18] Xu L J and Huang J P 2019 Phys. Rev. Appl. 12 044048
[19] Jiang J H, Lu S, and Chen J 2023 Chin. Phys. Lett. 40 096301
[20] Ouyang Y L, Yu C Q, Yan G, and Chen J 2021 Front. Phys. 16 43200
[21] Hu R, Iwamoto S, Feng L, Ju S, Hu S, Ohnishi M, Nagai N, Hirakawa K, and Shiomi J 2020 Phys. Rev. X 10 021050
[22] Hu R, Song J, Liu Y, Xi W, Zhao Y, Yu X, Cheng Q, Tao G, and Luo X 2020 Nano Energy 72 104687
[23] Yan S S, Wang Y, Gao Z B, Long Y, and Ren J 2021 Chin. Phys. Lett. 38 027301
[24] Roy Chowdhury P and Ruan X L 2022 npj Comput. Mater. 8 12
[25] Tan Y, Wang L, Wang Z, Peng J, and Ren J 2021 Chin. Phys. B 30 036301
[26] Yan S S, Liu Y, Wang Z, Lan X H, Wang Y, and Ren J 2023 Chin. Phys. B 32 057802
[27] Dhar A 2008 Adv. Phys. 57 457
[28] Zhai J X, Zhang Q Y, Cheng Z H, Ren J, Ke Y Q, and Li B W 2019 Phys. Rev. B 99 195429
[29] Jiang P F, Ouyang Y L, Ren W J, Yu C Q, He J, and Chen J 2021 APL Mater. 9 040703
[30] Yu C Q, Ouyang Y L, and Chen J 2022 Front. Phys. 17 53507
[31] Liu Y and He D H 2021 Chin. Phys. Lett. 38 044401
[32] Lepri S, Livi R, and Politi A 2003 Phys. Rep. 377 1
[33] Anderson P W 1958 Phys. Rev. 109 1492
[34] Ishii K 1973 Prog. Theor. Phys. Suppl. 53 77