Phononic Structure Engineering: the Realization of Einstein Rattling in Calcium Cobaltate for the Suppression of Thermal Conductivity

Phonons in condensed matter materials transmit energy through atomic lattices as coherent vibrational waves. Like electronic and photonic properties, an improved understanding of phononic properties is essential for the development of functional materials, including thermoelectric materials. Recently, an Einstein rattling mode was found in thermoelectric material Na0.8CoO2, due to the large displacement of Na between the [CoO2] layers. In this work, we have realized a different type of rattler in another thermoelectric material Ca3Co4O9 by chemical doping, which possesses the same [CoO2] layer as Na0.8CoO2. It remarkably suppressed the thermal conductivity while enhancing its electrical conductivity. This new type of rattler was investigated by inelastic neutron scattering experiments in conjunction with ab-initio molecular dynamics simulations. We found that the large mass of dopant rather than the large displacement is responsible for such rattling in present study, which is fundamentally different from skutterudites, clathrates as well as Na analogue. We have also tentatively studied the phonon band structure of this material by DFT lattice dynamics simulation, showing the relative contribution to phonons in the distinct layers of Ca3Co4O9.

In addition to materials engineering with phonon rattlers and nanostructures, phononic crystals (the artificial periodic structures made of two elastic materials) have been used to effectively suppress the thermal conductivity 9 . Periodic changes in density and/or elastic constants in phononic crystals can be manipulated in order to engineer the phononic bandgaps (ranges of the wavelength or frequency within which waves cannot propagate through the structure), filtering out phonons with forbidden frequencies 10 . The ZT value of 2.4 achieved in the Bi 2 Te 3 /Sb 2 Te 3 superlattice (a typical one-dimensional phononic crystal) was attributed to the significant suppression of thermal conductivity 11 . The assembly of PbTe and Ag 2 Te nanocrystals into a binary 3D superlattice not only reduces the thermal conductivity drastically, but also enhances electrical conductivity by about 2 orders of magnitude over the sum of individual conductivity of single-component PbTe and Ag 2 Te films 12 . Recently, a high-performance centimeter-scale 3D superlattice was intentionally structured with La doped SrTiO 3 nanocubes coated by 2D electron gas surface made of ultrathin Nb doped SrTiO 3 13 . Designing such 3D phononic crystals appears to be a promising approach to control the phonon transport while enhancing the electrical conductivity simultaneously. This raises the question of whether a phononic bandgap can be engineered locally by chemical doping, rather than by constructing phononic crystals.
To address this question, we use a layered thermoelectric material [Ca 2 CoO 3 ] RS [CoO 2 ] H δ (δ ≈ 1.61, commonly referred to by its approximate composition Ca 3 Co 4 O 9 ), as an example to explore the possibility of engineering the phononic structure locally. This compound has a layered structure consisting of disordered rock-salt type [Ca 2 CoO 3 ] layers alternating with electrically conductive hexagonal CdI 2 -type [CoO 2 ] layers [ Fig. 1(a)]. We were prompted to study this material following the observation of Einstein rattling in Na 0.8 CoO 2 , which has similar layered structure. The rattling of Na was found between [CoO 2 ] layers 14  ). All the as-prepared materials have similar grain size and relative densities (above 95% of their bulk theoretical density), minimizing the external effects such as porosity and grain size on their intrinsic phononic transport properties. Characterization by laboratory-based X-ray powder diffraction with Cu Kα radiation indicated that all these samples were single-phase with the Ca 3 Co 4 O 9 structure type (Supplementary Figure s1). In order to further study the dopant position in this compound, we carried out synchrotron XRD powder diffraction (with an X-ray wavelength of 0.82565 Å) on the undoped sample and co-doped sample. Rietveld refinement was carried out using Jana 2006 program. The initial model and refinement procedures have been discussed in our previous structural studies 15,16 on this material. The results indicated that Bi was preferentially doped onto Ca sites and Ga onto Co sites in the rock-salt layer [Ca 2 CoO 3 ] (supplementary information).
To study the influence of Bi and Ga substitutions on phonon transports, we have first measured the overall lattice dynamics using inelastic neutron scattering (INS) experiment on CCO and BG-CCO powder samples. Surprisingly, the measured GDOS spectra shown in Fig. 2(a) demonstrate very small modification of GDOS by Bi or Ga doping.
In order to reveal the underlying mechanism, in particular the role of Bi and Ga on the lattice dynamics, we also calculated the phonon density of states using ab-initio molecular dynamics (MD) simulation. Distinct from Na 0.8 CoO 2, Ca 3 Co 4 O 9 has more complicated structure, in which a misfit between b axis of the two layer types gives rise to the incommensurability δ = b RS /b H 17 . For the purposes of simulation, we have approximated the misfit structure to a commensurate model with δ = 8/5 = 1.6 and an overall composition of Ca 3 Co 3.9 O 9.3 . This structure is shown in Fig. 1(a). The starting configuration was then obtained by geometry optimization of the commensurate approximate structure shown in Fig. 1(a) and an equivalent model of BG-CCO with both dopants confined in the [Ca 2 CoO 3 ] block (two Bi atoms on Ca sites and one Ga on a Co site), shown in Fig. 1(b). To avoid the uncertainties of empirical force-field methods we have used density function theory (DFT) to determine the forces in MD simulations. The GDOS spectra calculated from the MD atomic trajectories are shown in Fig. 2(b). Although the incommensurate structure of Ca 3 Co 4 O 9 has been approximated to a commensurate model, the calculated GDOS spectra from MD simulation match the experimental GDOS spectra [ Fig. 2(a)] sufficiently well. Similar calculated phonon DOS spectra have been reported by Baran 18 and Rébola 19 , who carried out the DFT lattice dynamics simulations using a commensurate structure of 5/3 or even 3/2. The general agreement between experiment and simulation suggests that the influence of incommensurability on the lattice dynamics is most likely to be at rather small momentum transfer, which is beyond the current region of study. It also indicates that our MD trajectories can be used to deduce the influence of cation doping on the molecular dynamics of the system.
The partial GDOS spectra of the individual elements in both pristine CCO and co-doped BG-CCO are shown in Fig. 3. As a reference, the spectra of the individual elements in pristine CCO are analyzed. Figure 3(a) demonstrates that Ca in CCO mainly contributes to the energy range 10-40 meV. In contrast, the spectra of Co shown in Fig. 3(b) can be divided into two segments: 10-40 meV with a peak maximum at 28 meV, and a broad feature between 50-80 meV with a maximum at 65 meV, although its phonon density of states in the range 50-80 meV is relatively weak. Similarly, the partial GDOS spectra of O in pristine CCO [ Fig. 3(c)] also consists of two parts, 10-50 meV and 50-80 meV with maxima at 35 meV and 65 meV, respectively. The strong GDOS in the energy range 50-80 meV dominates the dynamic behavior of this particular element.
It may also worth looking at the dynamics from individual layers in this material. Since the overall lattice dynamics was almost not changed by doping (Fig. 2), we then investigated the dynamics from two distinct layers by simply considering the undoped material. We have carried out the DFT lattice dynamics simulation with harmonic approximation. The calculated phonon dispersion curves and phonon DOS of individual layers are present in Fig. 4. It is found that the obtained phonon DOS spectra from DFT lattice dynamics calculation [ Fig. 4(b,c) is quite similar with the spectra from MD simulation (Fig. 3a-c). It also indicates that there is good consistency between experiment, MD simulation and lattice dynamics calculation. The dispersion curves (Fig. 4a) show a very large number of branches, and a full analysis of any individual branches is beyond the scope of present work. However, it is encouraging to note that despite the use of approximate commensurate structure, there are only slight negative eigenvalues, which are at the gamma point (cut-off at zero frequency). The dispersion curve shows a gap at around 50 meV, which is also clear in the partial DOS plots illustrated in Fig. 4b,c. Inspection of these plots for Ca and O in the [Ca 2 CoO 3 ] layer (Fig. 4b) show rather more intensity below the gap, whilst for O in the [CoO 2 ] layer (Fig. 4c), the vast majority is above the gap. To show the relative contribution of Co in each layer, we have normalised the spectra of these two types of Co, subtracted Co_H from Co_RS and plotted the intensity difference in Fig. 4d. It shows that Co in the [Ca 2 CoO 3 ] layer has more intensity below 25 meV and less intensity in the range of 25-70 meV.  Returning to the GDOS spectra, a single peak at 10 meV dominates the GDOS of Bi [ Fig. 3(d)], which is fairly isolated from the vibrations of other elements in the compound, this being the hall-mark of an Bi Einstein-oscillator. It is also interesting to note that both experimental and calculated spectra show that the GDOS of BG-CCO is higher than that of pristine CCO in the energy range below 10 meV (Fig. 2). The concentration of Bi is low, but nevertheless, the higher intensity in this energy range for BG-CCO is consistent with the Einstein oscillation of Bi. Beyond this energy range, not only Bi but also Ga, Ca, and Co can contribute to the total phonon density of states.
In general, the vibrational kinetic energy E can be expressed as where m is the atomic mass, f is the frequency, and ADP is the atomic displacement parameters. Here, we assign U as ADP 2 to simplify the equation. Therefore, the vibrational kinetic energy E is proportional to m and U at a particular frequency. This makes m × U a key parameter to study the individual contributions of each atom to the total energy E. Figure 5 plots the m × U values against the anisotropic atomic displacement parameters U of each atom in the BG-CCO. It demonstrates that Bi has the largest m × U value among all the elements in BG-CCO, indicating that the Bi rattler absorbs a large portion of kinetic energy in the system. It is also worth noting that the parameters of Bi are comparable to those of the other atoms, implying that the large atomic mass is responsible for the rattling effects, not the large displacements. It is reasonable to consider that, since the mass of Bi is so different from other atoms in the lattice, it oscillates in the time-average mean-field of any other atoms, rather than responding to their dynamics. There could be some weak couplings between Bi and lattice via Ca site at frequencies below 10 meV (the local-mode energy of Bi rattler). But since the interaction is very limited, we only observed very small modification in both experimental INS spectra (Fig. 2a) and MD spectra (Fig. 2b) before and after doping, and the only significant change in the spectra is the additional shoulder around 10 meV due to Bi rattler. This new-type rattling is fundamentally different from skutterudites, clathrates and Na 0.80 CoO 2 etc., in which the displacement parameters of Na rattlers can be an order of magnitude higher than other atoms in the framework, but the mass of Na is rather similar to the other atoms in the lattice 20 .
Although mass-mismatch strategy has been used in alloy-based thermoelectric materials to hinder phonon propagation 21,22 , our system is too complex to be treated in the theories. In order to further study such mass-mismatch effect in our system, we carried out the ab-initio MD simulation with a substitute atom having the same electronic structure, atomic size and chemical environment of Bi but an atomic mass of 122 (equivalent to Sb). The calculated GDOS for Bi and this substitute atom are shown in Fig. 6. It shows that the characteristic peak of GDOS for the nominal atom, which is near half of the atomic mass of Bi, becomes much broader and the corresponding frequency is increased by ~2. This agrees well with the expectation from Equation (1) that the vibration frequency f should increase by 2 when the dopant losses its half mass, indicating that the change of vibrational kinetic energy caused by atomic displacements is very limited. It also suggests that the oscillation mass could be simplified to the dopant mass in the current system. In addition, the light substitute-atom shifts its vibrational frequencies to the range in which other elements vibrate, resulting in a more effective vibrational mixing. Most importantly, it further demonstrates that it is the atomic mass that plays the key role in determining the phonon transport in this particular layered compound.
To elucidate the influence of temperature on the Einstein oscillation, we further experimentally investigated the dynamical behavior in the energy range 0-20 meV at three different temperatures, 400 K, 800 K and 1000 K [ Fig. 7(a-c)] at PSI. The result demonstrates that the GDOS of BG-CCO at E < 10 meV is higher than that of its pristine counterpart at all the measurement temperatures, which is in good agreement with the simulated results. In order to provide a direct comparison, we also plotted the ratio of BG-CCO experimental GDOS spectra to the pristine CCO experimental GDOS spectra over the full range of energies at the measured temperatures in Fig. 8. The difference between the doping engineered BG-CCO and the pristine CCO spectra is distinct at energies below 10 meV, which is mainly contributed by the Einstein-like rattling of Bi. It is also noteworthy that such a local vibration mode persists at all temperatures, suggesting that the Bi is a good rattler, which is fairly independent of host lattice.
To further understand the consequence of the observed phononic behaviors in the as-prepared materials, the thermal conductivities as a function of temperature for the samples of CCO, Ga-CCO, Bi-CCO and BG-CCO were measured and are plotted in Fig. 9(a). This shows that both Bi and Ga doping can reduce the thermal conductivity over the measurement temperature range. The reduction was found to be more significant with the presence of Bi. In particular, the BG-CCO presents the lowest thermal conductivity, which is ~40% lower than that of the pristine CCO. We used the Wiedemann-Franz law to separate the individual thermal conductivities contributed to by hole carriers (κ h ) and phonons (κ ph ). The results are plotted in Fig. 9(b), in which it can be clearly seen that the reduction of total thermal conductivity is mainly due to the suppression of phonon transport through Bi rattler. On the other hand, the mass of Ga is similar to both Ca and Co, and in this case we do not observe a single sharp Einstein mode in the partial phonon DOS of Ga (Fig. 3d). Thus the reduction of thermal conductivity from Ga doping is less significant. In addition, we have conducted the measurement of electrical resistivity on both undoped and doped materials (Supplementary Table s1), which shows that the electrical resistivity was actually decreased by doping over the measured temperature range. This suggests that not only thermal conductivity was reduced by doping, but there is also a slight enhancement in the electrical conductivity, both being favorable for thermoelectric applications. A full understanding of electronic structure by chemical doping should pave the way to considerable improvement in the thermoelectric performance.
In summary, we have realized a new type of rattler in thermoelectric material Ca 3 Co 4 O 9 by substituting Bi for Ca site in the central block between [CoO 2 ] layers. Inelastic neutron scattering experiment shows that doping has very limited influence on the overall lattice dynamics. However, MD simulations identify the Bi rattler at 10 meV, which is local and fairly independent of host lattice. It remarkably reduced the thermal conductivity while slightly enhancing the electrical conductivity. More importantly, the rattling is found to be resulted from the big contrast of atomic mass between dopant and other elements in the lattice rather than the large displacement of dopant. This is essentially and fundamentally different from skutterudites, clathrates and Na x CoO 2 , which may bring new insights into the design of better thermoelectric materials.  with ethanol for 12 h. The mixed powders were then dried and calcined twice at 1173 K for 20 h in air with intermediate regrinding. The products were ball-milled again and placed into a 20 mm graphite mould. Spark plasma sintering was carried out in a Dr Sinter SPS-825 (Syntex, Inc., Japan) system under vacuum. Prior to sintering, a moderate pressure of ~10 MPa was applied to the mould to ensure a closed electrical loop for the current to pass through. The sample was then heated to 1073 K under a uniaxial pressure of 50 MPa and held for 5 min. The sintered pellets were annealed at 1173 K in air for 20 h and cooled to room temperature at a rate of 5 K/min to remove the graphite foil on the surface and re-oxidise them to the same state.

Methods
Measurement of thermal properties. Thermal diffusivity, heat capacity and thermal expansion were measured by a laser flash system (NETZSCH LFA-427), differential scanning calorimeter (NETZSCH DSC-404C) and dilatometer (NETZSCH DIL-402C), respectively. Thermal conductivity was evaluated from κ = DρC p where D is thermal diffusivity, ρ is density and C p is heat capacity.
Inelastic neutron scattering experiment. The inelastic neutron scattering (INS) experiments were performed using the time-of-flight neutron spectrometer FOCUS with an incident neutron wavelength of λ = 4 Å at the Paul Scherrer Institute in Switzerland. Energy transfer was measured up to 100 meV. Data for both pristine and co-doped samples were collected at 400, 800 and 1000 K. All the measured spectra were normalized to a standard vanadium sample and corrected for the scattering from the empty scan. The obtained data were Molecular dynamics simulation. The plane wave density functional theory (DFT) code (VASP) 25 was used for all MD simulations. The calculation consists of two parts: geometry optimization to determine the minimum-energy starting configuration, and then MD simulations to determine the dynamics at a given temperature. Production simulations in the microcanonical ensemble at 298 K, were developed from the minimum-energy configuration in the following way: (i) 4 ps of equilibration in the isokinetic ensemble (i.e., velocity scaling) with the temperature fixed to T; (ii) 4 ps of equilibration in the microcanonical ensemble; and (iii) at least 8 ps of production in the microcanonical ensemble with the thermodynamic temperature  fluctuating about T. In all VASP calculations, we adopted the projector augmented wave (PAW) potential 26 and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation function 27 . For all MD simulations, the energy cutoff was reduced to 300 eV and a time step of 1 fs was used.
The dynamical structure factor, S(q, w) was calculated from the atomic trajectories of the MD simulations using the code NMoldyn 28 . Firstly, the incoherent intermediate scattering function, F(q, t) was calculated on a regular grid in both momentum transfer q and time t. Then F(q, t) was Fourier-transformed to obtain the scattering function, S(q, w). In-house software was used to take the appropriate cut was taken through this surface to coincide with the experimental INS data, and this was then convoluted with the instrumental resolution function to obtain the calculated INS spectrum, which is compared with the observed spectrum in Fig. 2. Lattice dynamics simulations. The DFT + U were calculated using the projected augmented wave (PAW) 29 method with the Perdew-Burke and-Ernzerhof (PBE) 27 functional and the Dudarev's scheme 30 within the VASP code. A value of U = 4 eV was included in the Co's 3d to account for the strongly corrected effect in Ca 3 Co 4 O 9 as described in the previous study 31 . To calculate the phonon dispersion and the contribution of the individual ions in the phonon density of states (PDOS), a supercell consisting of 108 atoms with a mismatch ratio of ~8/5 between the CoO 2 layer and the rocksalt layer, with the experimental lattice constants obtained from the synchrotron diffraction study. The internal coordinates of the individual ions were first relaxed till the forces less than 0.004 eV/atom and an energy convergence of 10 −7 . A cut-off energy of 550 eV and a dense Monkhorst-Pack 32 k-point set of 4 × 2 × 2. The set of forces was calculated using the Finite Displacement Method on the phonon band structure (Fig. 4a) and partial DOS (Fig. 4b,c) were then extracted through the PHONOPY 33 software from the first principle calculation. The partial DOS was calculated using a Γ -centred k-point mesh with 960 kpoints.