Controlling thermal conductivity in nanocrystalline materials is of great interest in various fields such as thermoelectrics. However, its reduction mechanism has not been fully given due to the difficulty to assess local thermal conduction at grain boundaries (GBs) and grain interiors. Here, we calculated spatially decomposed thermal conductivities across and along MgO symmetric GBs using perturbed molecular dynamics, varying the GB separation from 2.1 to 20.0 nm. This reveals the different length scale of GB scattering for two directions: over hundreds of nanometers across GBs while within a few nanometers along GBs. Numerical analyses based on the spatially decomposed thermal conductivities demonstrate that the former is dominant upon suppressing thermal conductivity in polycrystalline materials, whereas the latter has a non-negligible impact in nanocrystalline materials because of a large reduction of intragrain thermal conductivity along GBs. These insights provide the exact mechanisms of heat transport in nanocrystalline materials toward more precise control of thermal conductivity.

Understanding and engineering of lattice thermal conductivity are important for a wide spectrum of applications, including thermal management in electronics,^{1,2} thermal barrier coatings,^{3} and thermoelectrics.^{4,5} Grain boundaries (GBs) are known to reduce the lattice thermal conductivity of polycrystalline materials because of their locally distinct atomic structures that scatter phonons incident from grain interiors.^{6–10} This fact has motivated many studies to synthesize nanocrystalline materials with grain sizes in a nano-meter order to achieve very low thermal conductivity, thereby improving the performance of thermal barrier coatings and thermoelectric materials.^{11–16} For example, yttria-stabilized zirconia with an average grain size $Dave$ = 10 nm exhibited a very low thermal conductivity of 1.2 W m^{−1 }K^{−1} at 480 K.^{11} In addition, the thermal conductivity of nanocrystalline silicon with $Dave$ = 9.7 nm is about 5% of that of single crystal at 300 K.^{16} These demonstrate that introducing a high number density of GBs is an effective strategy to reduce the thermal conductivity.

The lattice thermal conductivity of a solid can be interpreted by the following equation according to a kinetic theory,^{12}

where *ω* is the phonon frequency, *T* is the absolute temperature, and $C\omega ,T$, $v\omega ,T$, and $\Lambda eff\omega ,T$ are the spectral heat capacity, group velocity, and effective phonon mean free path (MFP), respectively. $\Lambda eff$ includes several phonon scattering mechanisms such as Umklapp and GB scatterings, as often combined by Matthiessen's rule: $\Lambda eff\u22121\omega ,T=\Lambda Umklapp\u22121\omega ,T+\Lambda GB\u22121\omega ,T$. Usually, thermal conductivity reduction by GBs is attributed to shortening of $\Lambda GB$, assuming that $\Lambda GB$ is equal to averaged grain size.^{12} GB scattering is also interpreted on the basis of the interfacial thermal resistance of GBs, $RGB$, as follows:

where $\kappa SC$ is the intragrain thermal conductivity, which is presumably equal to single-crystalline thermal conductivity.^{11,17} An assumption in this formula is that the thermal energy transmits through grain interiors as in single crystal until it reaches GBs.

Several experimental and theoretical studies have demonstrated, however, that thermal transport in nanocrystalline materials substantially deviates from the macroscopic description given above.^{12,16,18} Wang *et al.* and Jugdersuren *et al.* measured the thermal conductivities of nanocrystalline silicon with different grain sizes of 550, 144, 76, and 10 nm and analyzed them based on Eq. (1).^{12,16} This revealed that $\Lambda GB$ is always smaller than $Dave$ and decreased more rapidly than linear with decreasing $Dave$. Jugdersuren *et al.* suggested that this may be attributed to the reduction of intragrain thermal conductivity with decreasing grain sizes, which was not considered in their models. Indeed, Dong *et al.* modified Eq. (2) by substituting $\kappa SC$ by the grain-size dependent intragrain thermal conductivity, resulting in a better fit to $\kappa lat$ of nanocrystalline diamond and silicon.^{18} Spiteri *et al.* suggested that GBs along heat flow, which is not considered in Eq. (2), may reduce the intragrain thermal conductivity, based on the theoretical results for nanocrystalline diamond.^{19} These results indicate that more complex GB scattering occurs in nanocrystalline materials compared with ordinary polycrystalline materials. To elucidate the exact mechanisms behind thermal conduction in the nanocrystalline materials, it would be helpful to investigate nanoscale thermal transport at grain interiors and boundaries for both directions across and along GBs.

In this study, we examine thermal conduction in nanocrystalline MgO using the perturbed molecular dynamics (MD) method,^{20} which allows us to decompose the thermal conductivity to atomic (spatial) contributions. We choose MgO as a model material because of common sublattices to many oxides, high thermal conductivity (for ease of identifying the GB effect), and many experimental and theoretical works for comparison.^{9,10,21–26} Thermal conduction in MgO with grain sizes from 2.1 to 20.0 nm has been calculated. This reveals that, for the direction across GBs, thermal conductions at grain interiors and boundaries are strongly correlated, because GBs act as bottlenecks of heat transport and modify phonons spanning over dozens or hundreds of nanometers from the planes. Thermal conductivity also significantly decreases along GBs especially when the distance from GB is less than 10 nm, because transverse phonons that interact with GBs are well scattered along GBs. We also report how these nanoscale GB scattering mechanisms affect overall (averaged) thermal conduction in nanocrystalline materials toward further control of thermal conductivity.

The GB models used in this study include two identical GB planes separated regularly [Fig. 1(a)]. They are not nanocrystals but GB superlattices, because grains are infinite along GB planes under three-dimensional periodic boundary conditions applied in our simulations. However, we intentionally used these models in order to investigate GB scattering across and along GB planes separately, as previously done for diamond.^{19} After preliminary surveys over a wide range of GBs,^{9} we focus on the high-angle symmetric tilt Σ5(310)/[001] GB (tilt angle 2*θ* = 36.78°), whose structure was constructed by a simulated annealing method.^{9,26} The open GB core structure shown in Fig. 1(a) is in agreement with that obtained using first-principles calculations.^{27} To examine the GB separation dependence of the thermal conductivity, additional MgO slabs were added to the grain so that the separation between GB planes*, d*, is between 2.1 and 20.0 nm.

Effective thermal conductivities, which means overall thermal conductivities of the system including grain interiors and boundaries, were calculated using the perturb MD method at 300 K.^{20,28,29} To do this, a set of custom-written codes was added to the large-scale atomic/molecular massively parallel simulator (LAMMPS).^{30} In the perturbed MD, the lattice thermal conductivity along the *x* direction is calculated as

where *F*_{ext} is the magnitude of the perturbation and *J _{x}* is the heat flux along the

*x*direction. The microscopic heat flux

**J**is defined by Irving and Kirkwood as

^{31}

where **J**_{i} is the atomic contribution of atom *i* to the heat flux, *V* is the volume of the GB model, *m _{i}* and

**v**

_{i}are the mass and velocity of atom

*i*, respectively,

*Φ*is the interatomic potential energy between atoms

_{ij}*i*and

*j*,

**I**is a unit tensor of the second rank,

**F**

_{ij}is the force exerted by atom

*j*on atom

*i*, and

**r**

_{ij}is the position vector to atom

*j*from atom

*i*. The perturbed MD method is one kind of non-equilibrium MD (NEMD) as it exerts external forces to atoms for slightly enhancing the heat flux along a specific direction. The additional thermal energy is immediately removed from the system via thermostat, avoiding the temperature increase and gradient. We employed a set of Buckingham potentials for MgO reported by Landuzzi

*et al.*;

^{32}the thermal conductivity of MgO is well reproduced using this potential as discussed previously.

^{9}The detailed procedure for calculating the thermal conductivity is reported elsewhere.

^{9,10}

As seen in Eq. (3), the thermal conductivity is calculated from the time average of the heat flux in the perturbed MD. This means that the decomposing heat flux enables the estimation of partial contributions to $\kappa lat$, e.g., atoms and waves.^{33,34} Effective thermal conductivity can be considered as the summation of atomic thermal conductivities $\kappa i$ such that

where *J _{i,x}* is the contribution of atom

*i*to the heat flux in the

*x*direction. Atomic thermal conductivities must be normalized by multiplying the supercell volume

*V*when compared with different-sized GB models, because they are proportional to the inverse of

*V*as seen in Eq. (4). In this study, atomic thermal conductivities were projected onto the two-dimensional map and one-dimensional line as shown in Fig. 1 to investigate the spatial dependence of thermal conductivities. These projected thermal conductivities depend on the number of atoms in the depth direction ([001] here) or on the plane normal to GB, and thus were corrected by dividing the cell depth or the GB cross-sectional area.

In addition, the heat flux in the *x* direction $Jx$ can be divided into the contributions of a longitudinal wave $Jxx$ and transverse waves $Jxy$ and $Jxz$ according to

where *F _{ij,α}* and

*v*are the

_{i,α}*α*component of

**F**

_{ij}and

**v**

_{i}. Wave-decomposed thermal conductivity can be also calculated using Eqs. (6)–(8).

First, we examined the spatial and directional dependence of thermal conduction by decomposing effective thermal conductivity into atomic thermal conductivities and projecting them onto the plane perpendicular to the [001] direction. Figures 1(b) and 1(c) show the results of thermal conduction across and along GB planes, respectively, with *d* = 8.1 nm. For the direction across GB planes [Fig. 1(b)], atomic thermal conductivities were uniformly low in grain interiors and relatively high for atoms on the GB plane. The reason is that the single pairs of atoms at GBs are only the sites bonded across the plane, and heat is mainly transferred through these sites.^{10} Phonon thermal transport in the grain interiors is confined by the thermal conduction “bottlenecks” at GBs, resulting in uniformly low atomic thermal conductivities. This phenomenon is also seen in other high-angle tilt MgO GBs, where GBs consist of open structures with low-thermal-conductivity voids like in the Σ5(310)/[001] GB.^{9} On the other hand, along GB planes [Fig. 1(c)], atomic thermal conductivities significantly depend on the distance from the GB plane. The values were very low on the planes but increased with distance from the plane, demonstrating that GB scattering is effective within only a few nm along GBs.

To investigate how GB separation modifies nanoscale thermal transport near GBs and grain interiors, we projected atomic thermal conductivities onto the direction perpendicular to GBs, which we refer to one-dimensional (1D) thermal conductivity. The results for different *d* values were plotted in Figs. 1(d) and 1(e) as a function of the normalized *x* position, where GBs lie at *x* = 0.25 and 0.75. Across GBs [Fig. 1(d)], thermal conduction at grain interiors is significantly suppressed especially with small *d* values. The 1D thermal conductivity was found to be almost constant at grain interiors and boundaries when *d* is short, while that with relatively large *d* (8.1, 14.8, and 20.0 nm) increased with distance from the plane. The gradient is, however, very small even for *d* = 20.0 nm, demonstrating the long-range GB effect on the intragrain thermal conduction or phonons therein. This effect probably extends to hundreds of nanometers away from the boundary as roughly estimated from Eq. (2); the effective thermal conductivity of *d* = 100 nm is still as low as 65% of the single-crystalline thermal conductivity (Fig. S1). Interestingly, it was found that the thermal conductivity near GBs is not independent from intragrain thermal conductivity but increased simultaneously. This can be attributed to the delocalized nature of phonons. When *d* is sufficiently large, the thermal conductivity in grains should reach the single-crystalline thermal conductivity. In this case, the difference in the thermal conductivity between GBs and grains becomes constant regardless of GB separation. Then, GB scattering can be treated as the constant interfacial resistance as in Eq. (2).

In contrast to the thermal conduction across GB planes, that along GB planes was found to exhibit the significant spatial variation [Fig. 1(e)]. 1D thermal conductivity near GBs is very low and almost constant regardless of *d* values, because the open GB core structures that are repeatedly aligned determine the local thermal conduction along the plane. The 1D thermal conductivity sharply increased away from the GB plane, demonstrating the relatively short-range GB effect along the planes presumably less than dozens of nanometers. When *d* is very small, however, the 1D thermal conductivity at grain interiors is uniformly low even along the GBs. To be more quantitative, we estimated GB- and intragrain thermal conductivity along the GB planes by summing atomic thermal conductivities over the GBs and grain interiors separately and correcting them by the ratio of volumes between GBs and grains (Fig. 2).^{35} The GB region is defined as within 0.5 nm from the plane. The GB thermal conductivity decreased slightly with decreasing *d*, whereas the intragrain thermal conductivity exhibited a drastic reduction from 72.6 to 26.5 W m^{−1 }K^{−1}. This result indicates that GB scattering along the boundary is one of the origins of the reduction of the intragrain thermal conductivity in nanocrystalline materials.

To investigate more in-depth mechanisms behind GB scattering, we estimated the normalized wave contributions to effective thermal conductivities as shown in Fig. 3 together with those in single crystalline MgO. Across GB planes, thermal conductivities were decomposed into three kind of waves, namely, longitudinal waves and two kinds of transverse waves, whose amplitudes are for [001] and $[13\xaf0]$ directions ([001] and $[13\xaf0]$ transverse waves). The ratio of longitudinal waves in GB superlattices is larger than that of the single crystal, indicating that the scattering for longitudinal waves is less significant than those for transverse waves. This may be because longitudinal waves propagate as one-dimensional vibrations and transmit via the single pairs of atoms bonded across the GB plane. On the other hand, along GB planes, the ratio of the [310] transverse waves gradually decreased with decreasing *d*. This is because this kind of transverse waves propagates via atomic vibrations perpendicular to GB planes and is more likely to interact with the GBs, compared with the [001] transverse waves, in which the atom vibrates parallel to the planes. In addition, the wave contributions of *d* = 20.0 nm and single crystal are much different for the direction across GB planes while they are almost the same for the direction along GB planes. This revealed again that the spatial range of GB scattering is highly dependent on the direction. We note that the selective suppression of thermal transport seen in Fig. 3 must be related to the reduction of the group velocity and/or relaxation time of phonons [Eq. (1)]. See also Fig. S2 for how group velocities are decreased by the presence of GBs.

So far, the effect of GBs on thermal conduction has been analyzed separately for the directions across and along the planes using GB superlattice models. Finally, we estimated their combined effects by performing numerical simulations of thermal conduction over nanocrystalline MgO based on Fourier's second law and the obtained 1D thermal conductivities (see the supplementary material for details). We consider cubic and rectangular (plate-like and columnar) grains surrounded by six GBs as seen in Fig. 4(a) and the inset of Fig. 4(b). The separations between GBs are denoted as *d _{x}* for the direction parallel to the heat flux and

*d*for the direction perpendicular to the heat flux. The essence of this simulation is that direction-dependent GB effects on the local thermal conductivity was approximated by Matthiessen's rule as

_{yz}where *x*, *y*, and *z* are the positions inside the grain and correspond to the distance from the GB plane for each direction. $\kappa x1D$ is the 1D thermal conductivity for the direction across GB planes [Fig. 1(d)], and $\kappa y1D(y)$ and $\kappa z1D(z)$ are the 1D thermal conductivities for the direction along GB planes [Fig. 1(e)]. Using Eq. (9), heat transfer from the heat source to heat sink was calculated using Fourier's second law [Fig. 4(a)], and the effective thermal conductivity of nanocrystalline MgO was obtained from the heat flux and temperature gradient at the steady state.

Figure 4(b) shows effective thermal conductivities of nanocrystalline MgO depending on the grain size and grain shape. When *d _{x}* =

*d*= 20.0 nm (cubic), an effective thermal conductivity is as high as 36.4 W m

_{yz}^{−1 }K

^{−1}. With decreasing only

*d*, the shape of grains approaches a plate-like shape and thermal conductivity decreases to 4.6 W m

_{x}^{−1 }K

^{−1}(by 84%). This large decrease originates from the large GB impact on thermal conduction across the boundaries. On the other hand, with decreasing only

*d*from 20.0 to 2.1 nm, the shape of grains approaches a columnar shape and the thermal conductivity decreases to 24.9 W m

_{yz}^{−1}K

^{−1}(by 32%). This demonstrates that GB scattering along the planes is not negligible for thermal conductivities of nanocrystalline MgO. Such columnar-like grains are often seen in the microstructure formed using electron-beam physical-vapor deposition (EB-PVD).

^{36}Even for the case of

*d*= 2.1 nm, the thermal conductivity decreased by 23% with decreasing

_{x}*d*from 20.0 to 2.1 nm (from plate-like to cubic). This shows that the effect of GB scattering along the planes is large even for very small grains and probably causes a rapid reduction of $\Lambda GB$ and thermal conductivity in nanocrystalline materials.

_{yz}In summary, we have investigated nanoscale thermal transport in MgO GBs with different GB separations from 2.1 to 20.0 nm using the perturbed MD. Spatially decomposed thermal conductivities revealed that, for the direction across GBs, the thermal conduction at GBs and grain interiors are not independent from each other but have a strong correlation, because of the long-range GB effect that limits intragranular phonon transport over dozens or hundreds of nanometers. Although GB scattering along the planes was found to be short and weak compared with that across the planes, it caused the significant reduction of the intragrain thermal conductivity especially within a few nm from the planes. Our numerical analyses based on Fourier's second law revealed that not only grain size but also grain shape and morphology are important for determining the thermal conductivity of nanocrystalline materials because of non-negligible impact of GBs parallel to the heat flow. We would like to emphasize that the magnitude of GB scattering is different depending on materials and GB structures,^{10,19} and more extensive investigations using our approach will lead to a general understanding of thermal conduction in nanocrystalline materials and more accurate control of thermal conductivities.

See the supplementary material for additional discussions on thermal conduction in GB superlattice models and details of numerical simulations of thermal conduction in nanocrystalline MgO.

This work was supported by KAKENHI from the Japan Society for the Promotion of Science (JSPS) (Grant No. JP19H05786). S.F. was also supported by KAKENHI (Grant Nos. JP20H05195 and JP20K15034). T.Y. was supported by JST-CREST (Grant No. JPMJCR17J1). M.Y. was supported by KAKENHI (Grant No. JP20K05062).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.