OPEN

SUBJECT AREAS:

ATOMISTIC MODELS

SURFACES, INTERFACESAND THIN FILMS

Received 21 February 2014

Accepted 20 August 2014

Published 16 September 2014

Correspondence and requests for materials should be addressed to X.D. (dingxd@mail. xjtu.edu.cn); J.R. (jieustc@gmail.com) or E.K.H.S. (ekhard@esc.

cam.ac.uk)

Strain-controlled thermal conductivity in ferroic twinned films

Suzhi Li1, Xiangdong Ding1, Jie Ren2,6, Xavier Moya3, Ju Li1,4, Jun Sun1 & Ekhard K. H. Salje1,5

1State Key Laboratory for Mechanical Behavior of Materials, Xi'an Jiaotong University, Xi'an 710049, China, 2Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA, 3Department of Materials Science, University of Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, UK, 4Department of Nuclear Science and Engineering and Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA, 5Department of Earth Sciences, University of Cambridge, Cambridge CB2 3EQ, UK, 6Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.

Large reversible changes of thermal conductivity are induced by mechanical stress, and the corresponding device is a key element for phononics applications. We show that the thermal conductivity k of ferroic twinned thin films can be reversibly controlled by strain. Nonequilibrium molecular dynamics simulations reveal that thermal conductivity decreases linearly with the number of twin boundaries perpendicular to the direction of heat flow. Our demonstration of large and reversible changes in thermal conductivity driven by strain may inspire the design of controllable thermal switches for thermal logic gates and all-solid-state cooling devices.

Phonons are responsible for the transmission of sound and heat in the solid state. The ability to control phonon flow promises major technological developments in acoustics and thermal management1,2, cf. transistors that control charge flow in electronics. Solid-state thermal switches and diodes1 are particularly desirable for room-temperature cooling systems based on caloric materials3, as losses associated with the fluids that are currently employed to exchange heat between the working body and the heat sink or the heat load, respectively, reduce the performance of existing prototypes, and therefore cooling power4. Also, tunable thermal conductivities may be exploited for the design of advanced thermoelectric materials5 that display high values of the thermoelectric figure of merit ZT = SFoT/k over broad temperature ranges, where S, s, and k are the Seebeck coefficient, electrical conductivity, and thermal conductivity of the material, respectively, and T is the absolute temperature.

In recent years, a number of strategies have been exploited to demonstrate solid-state passive thermal diodes at room temperature6-10, but the rectification coefficients khigh/klow are small. Inhomogeneously mass-loaded carbon and boron nitride nanotubes6 display values of khigh/klow ~ 1 (1.02 for carbon nanotubes, 1.07 for boron nitride nanotubes), whereas few-layer graphite with asymmetric shapes display values khigh/klow < 2 (1.28 for triangular-shaped specimens7,1.6 for Y-shaped specimens8). Passive thermal rectification has been also demonstrated in bulk materials made of two oxides with different thermal conductivities9, but khigh/klow = 1.46 at best10.

Tunable thermal diodes offer active control of heat flow and therefore are highly attractive for the implementation of thermal circuits and solid-state coolers1. Besides various thermal diodes, control of thermal conductivity has been recently demonstrated in telescopically extended multiwall carbon nanotubes11 (khigh/klow = 6.7), and electrostatically actuated etch-released silicon nitride microfilaments12 (khigh/klow = 1.9). However, these intricate devices provide small areas for heat exchange, and are likely prone to fatigue during prolonged operation. Tunable thermal conductivity has been also demonstrated in VO2 nanorods13 (khigh/klow = 1.6), but substantial Joule heating associated with the electric currents that drive the metal-insulator transition in the nanorods compromise rectification performance.

Multiferroic materials constitute fertile ground to achieve sought-after large and reversible changes in thermal conductivity, as applied fields (electric field, magnetic field or stress) can control the number and orientation of phonon barriers associated with ferroic twin boundaries (TBs), whose density depends primarily on the magnitude and direction of the applied field14. Here we show that the thermal conductivity of ferroelastic twinned films can be reversibly controlled by strain to yield large thermal switching rate, due to the introduction and removal of twin boundaries that lay perpendicular to the direction of heat flow. Our strategy can be readily extended to ferroelectric and ferromagnetic twinned films whose twin boundary configuration can be controlled by electric and magnetic fields. Switchable domain boundaries can therefore replace previously exploited topo-logical defects (e.g. grain boundaries and kinks)15,16 that cannot be controlled by external fields, and lead to a new

www.nature.com/scientificreports \ > _

paradigm for the design of thermal devices. Also, unlike competing mechanisms to control heat flow, we avoid the use of power-dissipating currents (cf. Peltier elements), and moving parts.

We investigate coherent twin boundaries because they are most easily manipulated mechanically, electrically and magnetically. Our recent work has already shown that the morphologies and properties of twin patterns could be changed greatly by mechanically deformation1719. The important role of domain boundaries on materials properties20-23 has been previously evidenced beyond the field of phononics1. For example, superconducting twin walls exist in WO3 while the matrix is non-superconducting20; and ferroelastic systems such as Pb3(PO4)2 (ref. 24) and BaTiO3 (ref. 25), and martensitic alloys26 such as Ni-Ti and CuZnAl can take novel shape memory functionalities with the aid of domain boundary motion. Also, the interaction of domain wall and point defects largely determines the lifetime of ferroelectric memory devices27, while some ferroelectric twin boundaries are found in paraelectric bulk materials28, and may be tuned electrically.

We simulate the effect of mechanically driven domain boundaries evolution on heat transfer by using non-equilibrium molecular dynamics (NEMD) technique29. The simulation details are shown in Method. The initial condition is a sandwich twinned structure containing two horizontal twin boundaries, as shown in Fig. 1. The surface ratio of intermediate layer to the whole sample is fixed to be 0.7. The size of 2D simulation box is 160a X 100a in xy plane, where a is the lattice repetition unit. Periodic boundary condition is applied along the x direction and free boundary condition is used in the y direction. We apply different kinds of strain state to examine the twinning evolution in the mechanically driven system. The strain tensor is composed of [exx, Syy, yXy]T, where exx is tensile strain along x direction, Syy is tensile strain along y direction, and yxy is shear strain. In our work, the deformation is performed at different Syy/ yxy ratio with exx = 0. The simulations details are shown in Methods and Supplementary Section 1.

We first apply a simple shear deformation to the sample with strain tensor as [0, 0, yXy]T. Figure 2a shows the response of shear stress with shear strain (yxy) in a cycle. Under loading (black curve) the sample yields when yxy reached ~0.6%. The sample was unloaded to zero-strained state from yxy = 1.6% (red curve). Figure 2b shows

£yy A

Figure 1 | Thin-film model. A two-dimensional sandwich model with two fixed horizontal twin boundaries (HTBs). The lattice unit is in the shape of a parallelogram with tilt angle of 4 degrees. The middle layer has an area ratio of 0.7 to the whole sample. The strain tensor in 2D system can be described as [e„, e^, yx/]T. In our work, the deformation was performed at different ratio with exx = 0. Heat flow will be applied along x

direction. Atoms are colored according to values of their centro-symmetry parameters in range of 0 ~ 1 (ref. 41).

the variation of k during the loading/unloading loop. Here k is given in relative unit, r.u., related to the model potential we used here. We have not calibrated the absolute values to those which are obtained experimentally. In elastic regime, the magnitude of k does not change significantly (~140 to 150 r.u.) under stress. When loading into the plastic regime, k undergoes an abrupt drop to ~90 r.u., almost one half of the initial value. This magnitude is sustained under further stress. When the system is unloaded, k increases again and shows hysteretic behavior during the loading/unloading cycle. Our work here shows that by mechanically controlling the microstructures rather than the intrinsic properties of the single domain state, it is possible to generate two logic states, one with high heat conduction (such as points c, f in Fig. 2b) and the other with low heat conduction (such as points d, e in Fig. 2b), which can be used for thermal information storage and thermal switching13,30.

We examined the evolution of twin pattern for several typical strains (marked as (c), (d), (e), (f) in Fig. 2b) upon loading/unloading, as shown in Figs. 2c-f. When the system deforms plastically, new horizontal twin layers are induced by the applied shear strain, accompanied with the formation of a certain amount of vertical twin boundaries (VTBs). VTBs nucleate from one horizontal twin boundary (HTB), propagate and terminate in another horizontal twin boundary. The new-formed VTBs and HTBs superimpose and finally evolve into a much more complicated twin pattern. This pattern then shows much reduced thermal conductivity. Moreover, the structure oftwin pattern is quite stable and the twinning morphology will not undergo large changes even under the applied temperature gradient in NEMD calculations (see Supplementary Section 2). Upon unloading, horizontal twins are preserved and produce a permanent deformation of the sample. Vertical twins, on the other hand, are unstable and vanish gradually when removing the external load. This leads to an increase of k.

These results show that the reduction of thermal conductivity results from the existence of vertical twins with an orientation perpendicular to the heat flow. To further probe the exact role that VTBs and HTBs play in thermal transfer, we examined the variation of both VTBs density (rVTB) and HTBs density (rHTB) in loading/unloading cycle. Twin densities in our 2D system are calculated as the ratio of twin boundary area to total volume of system, which is in unit of a21. Figure 2g shows that the production and annihilation of vertical twins are correlated with changes of the thermal conductivity, i.e. the formation of VTBs reduces thermal transport. As the stress reaches the upper yield point (yxy = 0.8%), VTBs begin to nucleate, which corresponds to an abrupt drop of thermal conductivity. The magnitude of rVTB reaches to the maximum value at lower yield point (yxy = 1.0%) and is kept around that level under further shearing. Correspondingly, the thermal conductivity is also almost unchanged in this strain range of 1.0%-2.0%. In contrast, the presence of horizontal twins does not influence the heat transport, and thermal conductivity is independent of the density of HTBs (Fig. 2h). To distinguish the two twin structures with different heat transfer properties, we term the system containing only HTBs as twinning pattern 1 (TP1) and the system containing VTBs (purely vertical twin boundaries or mixed twin structure) as twinning pattern 2 (TP2). For the present case, the magnitude of k in TP2 is halved compared with TP1. The above calculations agree well with the experimental observation that the existence of twins in ferroelastic Gd2(MoO4)3 reduces thermal conductivity largely in comparison with a monodomain

sample31.

We have shown that VTBs density is the dominant factor for heat transportation properties. Since the twin patterns are sensitive to external strain state, it is possible to manipulate VTBs densities by applying different strain fields to tune thermal conductivity. When replacing the simple shear to a mixture of loading state [0, Syy, yXy]T, i.e. applying certain amount of tensile strain (Syy) along y direction besides shear strain (yxy), a much higher magnitude of resolved shear

Figure 2 | Strain-control of thermal conductivity in twinned films. (a) Stress-strain curve in a simple shear loop. (b) Variation of thermal conductivity with shear strain. (c)-(e) Representative atomic images corresponding to the different strain-scenarios marked in (b). Twin pattern 1 (TP1) in absence of vertical twin boundaries (VTBs) has a lower thermal conductivity than twin pattern 2 (TP2) with VTBs. Density of (g) VTBsand (h) HTBsas a function of shear strain, in units of a"1. Atoms are colored according to values of their centro-symmetry parameters41.

Figure 3 | Enhancing thermal-conductivity tunability using multiple strains. (a) Thermal conductivity variation with shear strain under different £yy/yxy ratio. (b) Variation ofVTBs density with shear strain. (c)-(e) Representative atomic images corresponding to the different strain scenarios marked in (b). Atoms are colored according to values of their centro-symmetry parameters41.

stress will be imposed on the VTBs planes and the density of VTBs can be enhanced greatly. Figure 3a and 3b show the variation of thermal conductivity and VTBs density with applied shear strain under different £yylyxy ratio, respectively. The VTBs density increases continuously when a tensile strain is applied. Correspondingly, the thermal conductivity is reduced dramatically. Figures 3c-e show the atomic configurations at yxy = 1.8% under different £yylyxy ratios. We find that the application of tensile strain effectively suppresses the formation of HTBs, and enhances the formation of VTBs. Furthermore, thermal conductivity can be reduced to 46 r.u. (see orange point in Fig. 3a), which is more than 30 smaller than the bulk values 1488 r.u. (see Supplementary Section 3).

The atomic mechanism of the reduction of thermal conductivity by VTBs is seen from the analysis of phonon spectra and vibrational eigenmodes32. We choose a sample with the system size of 80a X 50a to calculate the phonon density of states (DOS), keeping all other simulation conditions the same. Figure 4a shows the DOS for different values of pVTB extending to 27 THz. With increasing pVTB the spectral weight at high frequencies (15-27 THz) decreases. This

implies the twin boundary-phonon scattering by VTBs are at high frequencies while low frequency phonons are much less affected by VTBs (increasing is due to the normalization of the phonon spectra). To quantitatively describe the impact ofVTBs on the phonon transport, we calculate the participation ratio pi of each mode, which is

defined as33,34

pi {1=NE (£ e^i)2,

where N is total number of atoms and eia,i is the ath eigenvector component of each eigenmode i. The participation ratio describes the fraction of atoms participating in a mode and provides the information on localization effect of each mode. The value of pi varies from the O(1lN) for localization to O(1) for delocalization. We find that more phonons become localized when rVTB increases and correspondingly thermal conductivity is reduced (Fig. 4b). In particular, significant phonon localization occurs at high frequencies (around 18-27 THz). This is related to the highly distorted atomic structure at VTBs. During the twin patterns evolution, atoms inside

kinds of phonon scattering events, involving (a) phonon-phonon collision and (b) interactions of phonons and twin boundaries. The effective A is

1 _ 1 1 A Ap.p Ap -tb

Figure 4 | Phonon analysis of the influence of VTBs on thermal conductivity. Variation of (a) phonon spectrum and (b) participation ratio of systems under different pVTB during shear plasticity. (c) Bulk phonon dispersion relation. Red and black lines refer to the longitudinal and transverse branches, respectively.

twin boundaries displace and lead to phonon scattering at the smallest length scale of k = 0. 5 — 1p/a and corresponding frequencies around 18-27 THz (see phonon dispersion relation in Fig. 4c). In addition, Fig. 4c clearly shows that phonons of 18-27 THz, beyond the transverse mode, are specific to the longitudinal mode. Therefore, we see clearly that the emergence of VTBs strongly scatters the transport of the longitudinal phonon transport that in turn severely suppresses the longitudinal thermal transport.

In kinetics theory of phonon transport, the thermal conductivity is proportional to the mean free path of phonons as k = ACv/3, where A is the mean free path, C is the heat capacity and v is the acoustic velocity. In our system, the mean free path is determined by two

where AP-P and AP-TB are the characteristic lengths of the bulk pho-non-phonon scattering and the phonon-TBs scattering. The magnitude of AP-P in ferroelastic materials is usually on the order of 0.1 micrometers at low temperature31, which is much larger than AP-TB here. Thus, AP-TB plays the dominant role in determining the total mean free path A. For a constant AP-P, 1/k should be proportional to 1/AP-TB (1/k / 1/AP-TB). On the other hand, here AP-TB should be equal to the average twin boundary spacing (d), which is inversely proportional to the density of VTBs d = 1/pVTB. Therefore, we can finally obtain an approximate relationship of 1/k and pVTB as

K !Pvtb, (3)

To test this proportionality, we extract data within the strain range of 1.2-1.8% for different strain states in Fig. 3a and plot them in Fig. 5a. The data fit the linear relationship between 1/k and pVTB in Eq. 3 very well. It also manifests that the correspondence of 1/k and pVTB is not affected by the strain state we applied. However, it should be noticed that the intersection point of linear line with y axis (1/k0) does not refer to the bulk value (1/Kbulk) because systems with shorter box lengths (Lx) will underestimate the value of Kbulk due to the effect of periodic boundary condition on mean free path of phonons (Supplementary Fig. S3). The inset in Fig. 5b shows the dependence of simulation length on 1/k0 and the value converges to 1/Kbulk( = 6.7 X 10"4 r.u.). Moreover, the linear relationship described in Eq. 3 can be further verified in the twinned system with larger lengths. The slope 1/k to rVTB increases correspondingly and reaches 0.14 (orange line) in our calculations. Referring the slope to the bulk system as reference, this slope is predicted to become even steeper when removing the limitation of small length scales in our present MD simulations. Our results are compiled in Fig. 5b where the reference point of the untwinned materials is 1/Kbulk and all other values relate to nanostructured materials with a defined number of vertical twins during cold shearing and different sample sizes. Twin densities of 0.1a"1 (i.e. twin walls separated by 10 unit cells) are plausible in view of wall separations of 15 nm (some 40 unit cells) in SrTiO3 when the walls are thermally induced35. Cold shearing leads to much higher wall densities so that device applications become very likely in the light blue triangle in Fig. 5b. At a likely density of 0.1a"1 in light blue triangle, 1/k = 0.015 r.u. which is larger by a factor of 22 than the equivalent bulk value. An upper limit of the effect can be estimated when the twin density increases to 0.16a"1. The simulate contrast of the thermal conductivity reaches a factor Khigh/Kbulk = 30 in this case. However, it should be emphasized that the thermal switching rate is fixed for a given sized system. For example, Khigh/Kjow < 2.9 when Lx = 160a (see the black open symbols in Fig. 5b). Based on these results, we are able to approach the goal oftuning thermal properties in the framework of "domain boundary engineering" and estimate a maximum possible reduction of thermal conductivity by applied stress where the density of VTBs can be manipulated very precisely by external strain (or stress). This task becomes even easier with recently developed nanotechnology which fabricates materials with nanotwinned structure36. Since VTBs reduce thermal transport, we can hence tune thermal conductivity via the applied strain to the sample, which, in turn, can be manipulated by electric or magnetic fields in multiferroic materials.

Overall, our work shows how to achieve large and reversible changes in the thermal conductivity of multiferroic twinned films using applied fields. By controlling magnitude and direction of the applied field, we manipulate the density of twin boundaries that lay

www.nature.com/scientificreports \ > _

Figure 5 | The variation of thermal conductivity with twin boundary density. (a) Linear response of 1/k with VTBs density (pvtb) for different strain states at given box length (Lx = 160 a). (b) Variation of 1/k with Pvtb at different Lx.

perpendicular to the direction of heat flow and act as phonon barriers. Our strategy may be exploited for the implementation of fast controllable thermal devices, where heat-exchange surfaces can be largely increased by using multilayer geometries.

Methods

Our potential is constructed to be generic for all ferroelastic phase transitions by choosing the shear angle as "order parameter'' and can be used to study twinning and mobility of twin boundaries14 in ferroelastic materials. Comparing this configuration with the ferroelastic transition from cubic to tetragonal in SrTiO3 at 105 K, we focus on the modeling of a plane perpendicular to the twin plane (i.e. our axes are [101] and [1 01] in the Pm 3 m setting). Equivalent planes can be found in all ferroelastic systems. The Landau potential of such materials is implemented in the model to represent the interactions in a mono-atomic, two-dimensional lattice. The potential energy U(r) contains three parts: first-nearest atomic interactions of 20(r — 1)2, second-nearest interactions — 10(r — \/2)2 1 2000(r — \/2)4, and third-nearest interactions — (r — 2)4, where r is the atomic distance. Details of this potential are described in our previous work17-19,21. The equilibrium unit cell is in the shape of a parallelogram with shear angle of 4 degrees from square.

To study the structure evolution upon the external load in the sandwich sample, the top and bottom three atomic layers were fixed rigidly as the loading grip. We apply different kinds of strain state to examine the twinning evolution in the mechanically driven system. The strain tensor is composed of [£xx, £yy, CxyV, where £xx is tensile strain along x direction, Syy is tensile strain alongy direction, and yxy is shear strain. In our work, the deformation is performed at different eyy/yxy ratio with exx = 0. The dynamic loading is taken at T < 0.2Tm by using Nose-Hoover thermostat37,38, where Tm is melting point. Under such low temperature, any atomic diffusion process can be greatly suppressed. All the calculations are performed using the LAMMPS code39. The atomic configurations are displayed by AtomEye40 and the colors are shown according to the centro-symmetry parameters41.

For each loading step, we used the NEMD technique to calculate thermal conductivity in a given configuration29. The idea is to apply a heat flux in the system along the x direction, which will result in a temperature gradient. When the heat transfer becomes a steady flow, the thermal conductivity k along x direction is calculated via Fourier law as k = —Jxl(hT/hx), where Jx is the heat flux from heat source to heat sink and T is temperature. The induced temperature gradient at steady-flow ranges from 0.28Tm to 0.13 Tm, corresponding to the heat source and the heat sink (Supplementary Fig. S1). Note that real films, with larger unit cells and multiple atomic species, will exhibit much higher melting points and therefore permit room temperature operation.

Finally, we note that the present simulations were performed in 2D systems. Although our recent work shown that 2D simulations can reproduce the same thermodynamics properties in ferroelastic systems as in 3D simulations42, a study of three-dimensional materials would be highly desirable but also extremely costly in computer time.

1. Li, N. et al. Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond. Rev. Mod. Phys. 84, 1045-1066 (2012).

2. Maldovan, M. Sound and heat revolutions in phononics. Nature 503, 209-217 (2013).

3. Moya, X., Kar-Narayan, S. & Mathur, N. D. Nature Mater. 13, 439-450 (2014).

4. Yu, B., Liu, M., Egolf, P. W. & Kitanovski, A. A review of magnetic refrigerator and heat pump prototypes built before the year 2010. Int. J. Refrig. 33, 1029-1060 (2010).

5. Chang, H.-C., Chen, C.-H. & Kuo, Y.-K. Great enhancements in the thermoelectric power factor of BiSbTe nanostructured films with well-ordered interfaces. Nanoscale 5, 7017-7025 (2013).

6. Chang, C. W., Okawa, D., Majumdar, A. & Zettl, A. Solid-state thermal rectifier. Science 314, 1121-1124 (2006).

7. Tian, H. et al. A novel solid-state thermal rectifier based on reduced graphene oxide. Scri. Rep. 2, 523 (2012).

8. Zhang, G. & Zhang, H. Thermal conduction and rectification in few-layer graphene Y Junctions. Nanoscale 3, 4604-4607 (2011).

9. Kobayashi, W., Teraoka, Y. & Terasaki, I. An oxide thermal rectifier. Appl. Phys. Lett. 95, 171905 (2009).

10. Sawaki, D., Kobayashi, W., Moritomo, Y. & Terasaki, I. Thermal rectification in bulk materials with asymmetric shape. Appl. Phys. Lett. 98, 081915 (2011).

11. Chang, C. W. etal. Tunable thermal links. Appl. Phys. Lett. 90, 193114 (2007).

12. Supino, R. N. & Talghader, J. J. Electrostatic control of microstructure thermal conductivity. Appl. Phys. Lett. 78, 1778-1780 (2001).

13. Xie, R. et al. An electrically tuned solid-state thermal memory based on metal-insulator transition of single-crystalline VO2 nanobeams. Adv. Funct. Mater. 21, 1602-1607 (2011).

14. Salje, E. K. H. Ferroelastic materials. Annu. Rev. Mater. Res. 42, 265-283 (2012).

15. Jiang, J.-W., Yang, N., Wang, B.-S. & Rabczuk, T. Modulation of thermal conductivity in kinked silicon nanowires: Phonon interchanging and pinching effects. NanoLett. 13, 1670-1674 (2013).

16. Bagri, A., Kim, S.-P., Ruoff, R. S. & Shenoy, V. B. Thermal transport across twin grain boundaries in polycrystalline graphene from nonequilibrium molecular dynamics simulations. Nano Lett. 11, 3917-3921 (2011).

17. Salje, E. K. H., Ding, X., Zhao, Z., Lookman, T. & Saxena, A. Thermally activated avalanches: Jamming and the progression of needle domains. Phys. Rev. B 83, 104109 (2011).

18. Ding, X., Zhao, Z., Lookman, T., Saxena, A. & Salje, E. K. H. High junction and twin boundary densities in driven dynamical systems. Adv. Mater. 24, 5385-5389 (2012).

19. Zhao, Z., Ding, X., Lookman, T., Sun, J. & Salje, E. K. H. Mechanical loss in multiferroic materials at high frequencies: Friction and the evolution of ferroelastic microstructures. Adv. Mater. 25, 3244-3248 (2013).

20. Aird, A. & Salje, E. K. H. Sheet superconductivity in twin walls: Experimental evidence of WO3-x. J. Phys-condens. Mater. 10, L377-L380 (1998).

21. Ding, X. et al. Dynamically strained ferroelastics: Statistical behavior in elastic and plastic regimes. Phys. Rev. B 87, 094109 (2013).

22. Catalan, G., Seidel, J., Ramesh, R. & Scott, J. F. Domain wall nanoelectronics. Rev. Mod. Phys. 84, 119-156 (2012).

23. Salje, E. & Zhang, H. Domain boundary engineering. Phase Transit. 82, 452-469 (2009).

24. Aktas, O., Salje, E. & Carpenter, M. Resonant ultrasonic spectroscopy and resonant piezoelectric spectroscopy in ferroelastic lead phosphate, Pb3(PO4)2. J. Phys-condens. Mater. 25, 465401 (2013).

25. Aktas, O., Carpenter, M. A. & Salje, E. K. H. Polar precursor ordering in BaTiO3 detected by resonant piezoelectric spectroscopy. Appl. Phys. Lett. 103, 142902 (2013).

26. Otsuka, K. & Ren, X. Physical metallurgy of Ti-Ni-based shape memory alloys. Prog. Mater. Sci. 50, 511-678 (2005).

27. Scott, J. F. Fatigue as a phase transition. Integr. Ferroelectr. 38, 769-777 (2001).

28. Van Aert, S. et al. Direct observation of ferrielectricity at ferroelastic domain boundaries in CaTiO3 by electron microscopy. Adv. Mater. 24, 523-527 (2012).

29. MullerPlathe, F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys. 106, 6082-6085 (1997).

30. Wang, L. & Li, B. Thermal memory: A storage ofphononic information. Phys. Rev. Lett. 101, 267203 (2008).

31. Mielcarek, S. et al. Low-temperature thermal conductivity of ferroelastic Gd2(MoO4)3. Physica B 299, 83-87 (2001).

32. Kong, L. T. Phonon dispersion measured directly from molecular dynamics simulations. Comput. Phys. Commun. 182, 2201-2207 (2011).

33. Allen, P. B., Feldman, J. L., Fabian, J. & Wooten, F. Diffusons, locons and propagons: character of atomic vibrations in amorphous Si. Philos. Mag. B 79, 1715-1731 (1999).

34. Bodapati, A., Schelling, P. K., Phillpot, S. R. & Keblinski, P. Vibrations and thermal transport in nanocrystalline silicon. Phys. Rev. B 74, 245207 (2006).

35. Chrosch, J. & Salje, E. K. Near-surface domain structures in uniaxially stressed. J. Phys-condens. Mater. 10, 2817 (1998).

36. Lu, L., Chen, X., Huang, X. & Lu, K. Revealing the maximum strength in nanotwinned copper. Science 323, 607-610 (2009).

37. Nose, S. A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 81, 511-519 (1984).

38. Hoover, W. G. Canonical dynamics - equilibrium phase-space distributions. Phys. Rev. A 31, 1695-1697 (1985).

39. Plimpton, S. Fast parallel algorithms for short-range molecular-dynamics. J. Comput. Phys. 117, 1-19 (1995).

40. Li, J. AtomEye: an efficient atomistic configuration viewer. Modell. Simul. Mater. Sci. Eng. 11,173 (2003).

41. Kelchner, C. L., Plimpton, S. J. & Hamilton, J. C. Dislocation nucleation and defect structure during surface indentation. Phys. Rev. B 58, 11085-11088 (1998).

42. Zhao, Z., Ding, X., Sun, J. & Salje, E. K. H. Thermal and athermal crackling noise in ferroelastic nanostructures. J. Phys-condens. Mater. 26, 142201 (2014).

Acknowledgments

X.D and J.S. appreciate the support of NSFC (51171140, 51231008, 51321003, 51320105014), the 973 Program ofChina (2010CB631003,2012CB619402) and 111 project (B06025). J.R. acknowledges the support from National Nuclear Security Administration of the U.S. DOE at LANL under Contract No. DE-AC52-06NA25396 through the LDRD Program. X.M. is grateful for support from the Royal Society. J.L. acknowledges the support by NSF DMR-1240933 and DMR-1120901. E.K.H.S. is grateful for support by the Leverhulme fund (RG66640) and EPSRC (EP/K009702/1). Zhao Wang is thanked for helpful discussion.

Author contributions

X.D., J.R., J.S. and E.K.H.S. conceived and designed the project. S.L., X.D. and J.L. performed molecular dynamics simulations. J.S. provided the computational source. S.L., X.D., J.R., X.M. and E.K.H.S. wrote the paper. All authors contributed to discussions of the results.

Additional information

Supplementary information accompanies this paper at http://www.nature.com/

scientificreports

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Li, S. et al. Strain-controlled thermal conductivity in ferroic twinned films. Sci. Rep. 4, 6375; DOI:10.1038/srep06375 (2014).

ir m 11

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder in order to reproduce the material. To view a copy of this license, visit http:// creativecommons.org/licenses/by-nc-nd/4.0/