The Impact of Tilt Grain Boundaries on the Thermal Transport in Perovskite SrTiO 3 Layered Nanostructures. A Computational Study.

Thermal management at solid interfaces presents a technological challenge for modern thermoelectric power generation. Here, we define a computational protocol to identify nanoscale structural features that can facilitate thermal transport in technologically important nanostructured materials. We consider the highly promising thermoelectric material, SrTiO3, where tilt grain boundaries lower thermal conductivity. The magnitude of the reduction is shown to depend on compositional and structural arrangements at the solid interface. Quantitative analysis indicates that layered nanostructures less than 10 nm will be required to significantly reduce the thermal conductivity below the bulk value, and it will be virtually independent of temperature for films less than 2 nm depending on the orientation with a reduction of thermal transport up to 75%. At the nanoscale, the vibrational response of nanostructures shows concerted vibrations between the grain boundary and inter-boundary regions. As the grain boundary acts markedly as a phonon quencher, we predict that any manipulation of nanostructures to further reduce thermal conductivity will be more beneficial if applied to the inter-boundary region. Our findings may be applied more widely to benefit other technological applications where efficient thermal transport is important.

. Lattice parameters of the simulation cells for the 10nm-GB, 2nm-GB, and lattice dynamics (LD-GB) layered nanostructures separated by tilt grain boundaries. The a, b and c cell dimensions correspond to the direction x, y and z respectively.
System a x b x c (nm) α = β = γ (degree) 10nm-GB ∑ 3{111}/ [1 ̅ 10] 20.3 x 2.2 x 1.9 90 10nm-GB ∑ 3{112}/ [1 ̅ 10] 19. 3  Section S1. Calculation of the heat flux and thermal conductivity We used the Green-Kubo method 1, 2 as implemented in LAMMPS. 3 First the heat-flux of the system is calculated using: where is the heat-flux, is the energy of atom , is the velocity vector for atom and is the stress tensor. This may also be expanded to give: and further: where represents the force and the separation between atoms and . We have implemented a simple program to compute the autocorrelation of the heat-flux and then to integrate it, yielding a value proportional to the thermal conductivity: where is the integer number of steps between heat-flux sampling, ∆ is the timestep of the simulation, is the volume of the system, is Boltzmann's constant and is the temperature of the system. To achieve a complete accounting of the heat flux, we have used the method of Sirk et al. 4 (where a cubic 40Å simulation box of 2180 water molecules was deemed sufficient to avoid dependence on box length) of including a k-space contribution to the energy and virial, which is developed as an extension of particle-particle particle-mesh technique (PPPM).
Convergence test. The Green-Kubo method has been used in these calculations as convergence with respect to system size is reached relatively quickly compared to NEMD 5 as finite size effects are much reduced by the absence of a heat source/sink. Additionally, data for independent dimensions of a system can be obtained simultaneously during a single simulation using the Green-Kubo method. A final consideration is that in calculations using the Green-Kubo method there is no heat source/sink and so its placement relative to the grain boundary and the difficulty of extrapolating to infinite length in NEMD calculations is avoided. 5 We have performed convergence tests on the 2nm-GB ∑ 3{111}/ [1 ̅ 10]. We have calculated the thermal conductivity and the heat-flux autocorrelation function (HFACF) spectra for a 1 x 2 x 2 expansion of the 2nm-GB ∑ 3{111}/ [1 ̅ 10] (Table S1 present the 1 x 1 x 1 unit cell). The thermal conductivity of a 1 x 1 x 1 and of the 1 x 2 x 2 unit cells are all within the fluctuation bounds. Figure S1 presents the comparison. Fluctuation bounds for calculations using the Green-Kubo method can be very large, greater than 20% even in well converged data as demonstrated by Schelling et al.. 5 Additionally, specifying the point at which the Green-Kubo integral has converged can be extremely challenging as stated by McGaughey and Larkin. 6 We have used the bottle-neck regime to calculate the fluctuation bounds as explained by McGaughey and Kaviany. 7 This is currently the only way available to do so using the Green-Kubo method to calculate fluctuation bounds of thermal conductivity values.
The size of the simulation cell imparts periodic constraint on how atoms may vibrate and that other atoms may display spurious correlations due to the periodic images, but this effect is minimised by two key factors: (1) the thermal conductivity of bulk SrTiO3 is converged with respect to simulation cell size, indicating spurious vibrational correlations are negligible in bulk SrTiO3 at this size, and (2) the scattering in the grain boundary systems is significantly increased compared to the bulk material due to the presence of grain boundaries and thus the phonon mean free path is reduced, allowing a smaller simulation size than is needed for a pure bulk material. 6 To support that our values for the thermal conductivity are also converged, one can think about the mean free path of phonons. Calculation of mean-free-path is performed by using a simple kinetic approximation as explained by Dove: 8 This expression may be rearranged to give the phonon mean-free-path, . The thermal conductivity, , is calculated via molecular dynamics calculations of bulk SrTiO3 using the LAMMPS code. 3 A thermal conductivity value of 8 W/(m.K) at 500 K was obtained. 9 The heat capacity per unit volume, , and the average speed of the phonons, 〈 〉, are calculated from lattice dynamics calculations as implemented in the GULP code. 10 All properties were calculated using the Teter potential model 11 and an energy minimised unit cell of stoichiometric bulk SrTiO3. The heat capacity per unit volume for bulk SrTiO3 is calculated to be 3,085,136 J/(K.m 3 ). Lattice dynamics calculations also yield the velocities of two types of phonons, S-waves ( , transverse, 4,822 m/s) and P-waves ( , longitudinal, 8,352 m/s) computed from the bulk ( ) and shear modulus ( ) of the material along with the density ( ): Using these values the faster P-waves indicate a phonon mean-free-path in bulk SrTiO3 of ~9.3 Å and the slower S-waves give a value of ~16.1 Å. Both of these values calculated for bulk SrTiO3 are shorter than the size of all our simulation cells (in all directions) (Table S1) used for calculation of thermal conductivity in our layered nanostructures.

Section S2. Identification of peaks in the heat-flux autocorrelation function (HFACF) spectra
Models. To analyse the contribution of each vibrational mode to the total heat-flux from the different regions (i.e. grain boundary and inter-boundary regions) of our simulated systems (i.e. layered nanostructures separated by tilt grain boundaries), it is first necessary to define the regions within the system that contribute to the "Grain Boundary" and to the "Inter-Boundary", referred to as GB and IB, respectively.  Table S1. There are two of each region in the layered nanostructures, due to the way we construct the configurations. 12 This is shown in Figure S2 for the layered nanostructure containing the ∑ 3{111}/[1 ̅ 10] grain boundary. We have imposed an additional constraint on the definition of the regions so that it must be possible to construct a single crystal SrTiO3 using only the inter-boundary regions IBA or IBB. Thus, by taking either region IBA or IBB alone, it is possible to reconstruct a bulk SrTiO3 system orientated in the [111] direction.

Methodology.
To analyse the contribution of the different regions, GB and IB, on the scattering effect of the different vibrational modes, first the Г-point vibrational modes are calculated along with their corresponding eigenvectors using the METADISE 13 and the Phonopy 14-16 codes. The eigenvectors for each vibrational mode given by Phonopy are then summed in each dimension (i.e. x, y and z). Modes with an eigenvector sum below a small cutoff of intensity 0.01 a.u. are not included in our analysis in order to avoid background noise. This enables us to generate the phonon density of state (PDOS) for the three x, y and z dimensions separately, expressing only the vibrational mode with non-zero sum of the eigenvectors. In our previous work, 9 we have demonstrated that this procedure allows us to compare the PDOS calculated using lattice dynamics with the spectrum of the heat-flux autocorrelation function (HFACF) derived from molecular dynamics calculations. The vibrational modes with an eigenvector sum of zero in all three dimensions will not appear in the spectrum of the heatflux autocorrelation function. This is because all vibrational motions are exactly cancelled by an anti-symmetric motion of a symmetry equivalent atom elsewhere in the system. Therefore in these vibrational modes the centre of mass of the system does not change (zero momentum), and additionally there cannot be net heat-flux during the vibration as the energy moving in a given direction is exactly cancelled by energy moving in the opposite direction. These vibrational modes are "Inactive Modes". When the sum of eigenvectors is non-zero, there is no exact cancellation of vibrational motions over the entire system. While the momentum of the simulated system is still zero, there is now a disparity in the number of species moving in a given direction at any given time. This disparity gives rise to a net system heat-flux oscillating in a specific direction over time. These vibrational modes are "Active Modes". A rigorous mathematical derivation of this principle is given in the appendix of Landry et al. 17 Figure S3 shows the phonon density of states at the Γ-point before and after removing the inactive modes for stoichiometric bulk SrTiO3 (i.e. equivalent to a stoichiometric infinite single crystal). Figure S3a shows the PDOS with four Longitudinal Optical (LO) Modes at 5.7 THz, 11.0 THz, 14.1 THz and 22.4 THz. One mode is removed as its eigenvector sum is equal to zero, and therefore the mode will not be associated with either a net change in heat-flux or in dipole moment (IR inactive mode). Figure   In order to provide some statistical analysis of the nature of the vibrational modes in the layered nanostructures, we have devised analyses that can divide the vibrational modes depending on whether they have a predominant character arising from the IB or GB regions, referred as "Per-Region Analysis", whether they have an in-plane (parallel to the grain boundary) or out-of-plane (perpendicular to the grain boundary) character, refereed as "Per-Direction Analysis", and whether they involved predominantly Sr, Ti or O species, refereed as "Per-Species Analysis". These analyses are based on Figure 5 in the main manuscript. All data is presented in Tables S1, S2 and S3 outline the frequency of each vibrational mode, the "Direction", the "Region" and the "Species" for each layered nanostructures containing the ∑ 3{111}/[1 ̅ 10], ∑ 3{112}/[1 ̅ 10] and ∑ 5{310}/[001] grain boundary respectively.
First of all, discussion from this point onwards is concerned only with active modes shown in the HFACF spectra. However, these spectra are of difficult interpretation and thus why we have evaluated the Phonon Density of States (PDOS) using lattice dynamics calculations. As we mentioned above, lattice dynamics calculations provide both the eigenvalues (the vibrational frequencies, ) and the associated eigenvectors to the vibrational frequency. These associated eigenvectors for a given vibrational frequency (eigenvalue, ), can be seen as a × 3 matrix, , where is the number of atoms in the system and 3 corresponds to the directions of the simulation cell , and . The eigenvectors correspond physically to the motion of the atoms in the system associated to a particular vibrational mode, and thus the matrix is the collective vibrational motion of the entire structure. By defining , the different components of the vibrational mode can be isolated. may be broken down in many ways but we have chosen the "Per-Region Analysis", "Per-Direction Analysis", and "Per-Species Analysis".
The matrix takes the form: We can express in terms of a list of vectors, each of length 3, corresponding to each atom in the system: If any component of is non-zero, then the corresponding frequency, , will appear in the spectrum of the HFACF. Using this approach a phonon density-of-states (PDOS) may be easily constructed from the lattice dynamics calculation, yielding direct correspondence to the spectrum of the HFACF. The vibrational mode is thus assigned a directional component of either , or depending on the maximum absolute component of . If | , | > | , |, | , | then the vibration has a predominant direction in x, if | , | > | , |, | , | then the vibration has a predominant direction in y, and if | , | > | , |, | , | then the vibration has a predominant direction in z. This predominant direction component is expressed as , and in the "direction" column of Tables S2-S4. Figure 5a in the main manuscript is simply constructed by counting all vibrational frequencies that have a given directional component and expressing them as a percentage. The direction is expressed as perpendicular (⊥), while the and directions are summed into the parallel component (∥). The parallel component (∥) has been normalised to account for the existence of two directions (the count is divided by two).
Per-Region Analysis. To evaluate the contributions of different regions, , of the system, we first classify each atom as belonging to either the Inter-Boundary (IB) or Grain-Boundary (GB) region. The IB region is the region between two boundaries, whereas the GB region is the region comprising the grain boundary. The number of atoms in a region is denoted by . For each vibrational frequency, a summation of the vectors , can be performed on the atoms belonging to IB and GB regions separately. , = [ , , , , , , , , ] = ∑ ,

=1
By calculating the magnitude of the vector , it is possible to define which region has the stronger contribution. , = ‖ , ‖ This will define whichever contribution is greater, , or , , which will determine whether the vibrational mode ( ) is classified as being dominant in the GB or IB region. Thus, for every vibrational frequency if , > , , the vibration has a predominant contribution from the grain boundary (GB), whereas if , < , , the vibration has a predominant contribution from the inter-boundary (IB) region. This classification can be found under "Region" in Tables S2-S4. As before the number of vibrational modes of each character (GB and IB) are counted and expressed as a percentage, this represents the location contribution in Figure 5b in the main manuscript.
It must be noted that although the different characters may arise from the IB region, the vibrational mode may only occur due to the presence of the GB region. This is because of the presence of extended defects (i.e. grain boundaries) induce complex vibrational motions in the layered nanostructures. Our analysis is able to differentiate between vibrational modes active mainly in the grain boundary (GB) or the inter-boundary (IB) regions under the constraints of using equal numbers of atoms in the two regions and self-periodicity of the inter-boundary regions, as stated in the "model" section. The sums of eigenvectors for the two regions (i.e. GB and IB) may also agree or disagree in their sign. We term this "directional agreement" and for each mode it is either "TRUE" (i.e. agree) or "FALSE" (i.e. disagree). We can also indicate the relative strength of the contribution of the different regions in terms of a simple ratio of the sums of the eigenvectors in each region (i.e. IB and GB) for the same vibrational mode, where the maximum value for the contribution is set to 1.0. For example, it can be seen in Table S1 that the two modes at 2.71 THz occur in the "direction" and , which are equivalent.
They are dominant in the inter-boundary region "IB" (set to +1.0), but the grain boundary region "GB" has a response nearly equal in magnitude in the opposite direction indicated by a negative sign (-0.975).
Per-Species Analysis. We can also extend the analysis of the vibrational modes by performing the summation of the eigenvectors , over the different species (i.e. Sr, Ti or O) in the system, rather than in regions. For each vibrational mode ( ) we can sum the eigenvectors of all Sr, or Ti or O species separately; we refer to the species with the index .

=1
It is then possible to assign a vibration type (i.e. the species involved the most in the vibration Sr, Ti or O) to each vibrational mode . This is not as simple as calculating which species has the largest magnitude (i.e. contributes the most to the vibration), as this is influenced by the stoichiometry of the system. Instead, the value , , is extracted from , , for each species, where is the direction of the vibrational modes (whether , or ). For each species, the sign of , , indicates whether the atoms of that species are moving in the + ordirection on average. By comparing the signs of , , for all species at the same vibrational frequency, it is possible to group together species which on average move together. An example is provided. For a vibrational mode : From these values it is clear that, on average, Ti and O move together while Sr moves in the opposite direction. Thus this mode will be assigned as Sr/Ti-O in "Species" in Tables S2-S4. In this example, the Sr atoms move in the opposite direction to the Ti and O atoms, and so this may be considered akin to a "Sr rattling" type mode. However, the greatest contribution to the heat-flux is actually form the Ti-O framework. This can be rationalised in terms of the total momentum which must be zero for any vibrational mode. Therefore any "rattling" type vibration of an atom must be counterbalanced by the opposite movement of the Ti-O framework, ensuring zero momentum. As the Ti-O network has a larger number of atoms and the time average kinetic energy of atoms is equal irrespective of species, then it follows that the biggest contribution to the heat-flux is from the Ti-O network. As before the number of vibrational modes of each character (primary Sr, Ti or O character) are counted and expressed as a percentage. This represents the "Species" contribution in Figure 5c in the main manuscript.   Section S3. Thermal conductivity of layered nanostructures -the directional components Figure S4 presents the in-plane (∥) and the out-of-plane (⊥) thermal conductivities for all the layered nanostructures studied. The in-plane (∥) thermal conductivity corresponds to the average thermal conductivity calculated in the Y (i.e. the direction of the grain boundary pipe) and Z directions of the simulation cells, thus parallel to the grain boundary plane. The out-of-plane (⊥) thermal conductivity corresponds to the contribution calculated across (i.e. perpendicular to) the grain boundary plane, which is the X direction in the simulation cells.