Mechanical Controlled Thermal Switch and Hysteresis with Domain Boundary Engineered Phonon Transport

Heat flow control in phononics has received significant attention recently due to its widespread applications in energy transfer, conversion and utilization. Here, we demonstrate that by applying external stress or strain we can effectively tune the thermal conductivity through changing the density of twin boundaries, which in turn offers the intriguing mechanical-controlled thermal switch and hysteresis operations. Twin boundaries perpendicular to the transport direction strongly scatter phonons. As such, the heat flow is in inverse proportional to the density of twin boundaries and hence allows an excellent way to switch thermal conductivity mechanically and even leads to the interesting hysteresis behavior as a thermal memory. Our study provides a novel mechanism to couple thermal and mechanical properties of materials as a matter of"domain boundary engineering"and can have substantial implications in flexible thermal control and thermal energy harvesting.

Heat conduction and electric conduction are two fundamental energy transport phenomena in nature. The development of modern electronics to control charge transfer has impacted every aspect of our everyday life. However, heat conduction has never been treated on the equal footing and its flexible control is largely absent, in spite of the fact that almost 90 percent of the world's energy are utilized through heat processes. [1] To meet the demand, a new discipline, phononics emerges, [2] which is the science and technology of phonons that are main carriers of energies in non-metals. Phononics aim to manipulate phonons and to control heat flow and thermal energy as flexibly as in electronics. To achieve this ultimate goal, various functional thermal devices, like thermal diodes, thermal transistors, thermal logic gates have been proposed theoretically and partially been realized in experiments. [2] It was demonstrated that phononic logic gates can perform similar operations as their electronic counterparts which makes it possible that phonons can be even used as the information carriers. [3] In this sense, thermal memory materials [4] can be designed if the micro-structural conditions for functional materials are better understood. Nevertheless, candidate materials for potential applications have yet to be identified. Moreover, the controls of those phononic devices are mainly based on tuning temperatures. More flexible manipulating protocols are in demand. Considering that microstructure changes can tune the lattice vibrations, we can effectively manipulate the microstructures to control the phonon fluxes so as to generate diverse device applications based on the heterogeneous phonon properties of a material.
It is the purpose of this paper to show that by mechanically governing 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 and the other with low heat conduction, for thermal information storage and thermal switch. This idea is intuitive since it is known that phonon fluxes are strongly influenced by internal defects. Among all the types of defects, such as point defects, line defects, etc, interfacial boundaries are the ones that can best influence the thermal transport and are susceptible to external fields such as elastic stresses and strains in ferroelastics, [5] electric fields in ferroelectrics [6] and magnetic fields in ferromagnets. Mobile interfaces will, in each case, modify the phonon scattering and hence allows the modulation of the heat flux. For instance, recent experiment shows that the conductivity of Si-Ge alloy thin films is 3 ~ 5 times lower than the bulk values. [7] The reason lies in that long wavelength phonons are scattered strongly in film boundaries. Previous research also focused on the thermal properties of other 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   4 types of interfaces, [8][9][10] such as bi-crystal interfaces, tilt grain boundaries, and free surfaces.
None of the above mentioned interfaces can be tuned as we propose in this paper.
We investigated coherent twin boundaries because they are most easily manipulated mechanically. Our recent work has already shown that the morphologies and properties of twin patterns (TPs) could be changed greatly by mechanically deformation, [11][12][13] which leads us directly towards ways to tune thermal conductivity via optimizing twin patterns and further develop thermal memory devices. The important roles of domain boundaries on materials properties have already shown in many other fields, [14][15][16] which is known as "domain boundary engineering". [17] For example, weak doping generates superconducting twin walls below 3K in WO 3 while the matrix is non-superconducting. [14] In ferroelastic systems, certain alloys can take novel shape memory behaviors with the aid of domain boundary motion. [18] The interaction of domain wall and point defects largely determines the lifetime of ferroelectric memory devices, [19] while ferroelectric domain boundaries were found in paraelectric bulk materials. [20] We will show in this paper that within the same framework of "domain boundary engineering", [17] it is possible to manipulate the thermal conductivity and obtain two stable states for information storage by optimizing domain structures and domain patterns. We studied the effect of mechanically driven domain boundaries evolution on heat transfer by using non-equilibrium molecular dynamics (NEMD) technique. [21] We show that thermal conductivity is strongly related to the morphology of twin patterns. By applying different strain/stress state, we can effectively manipulate morphology of twin structure and eventually tune thermal conductivity by a factor 2 ~ 4.
The size of 2D simulation box is 160a × 100a (= 16 nm × 10 nm) in xy plane, where a (= 1 Å) is the lattice repetition unit. Periodic boundary condition is applied along the x direction and free boundary condition is used in the y direction. To study the structure evolution upon the external load, the top and bottom several layered atoms were fixed rigidly as the loading grip. We applied different kinds of strain state to examine the twinning evolution in mechanically driven system.
The strain tensor in 2D system can be described as [ε xx , ε yy , γ xy ] T . In our work, the deformation was performed at different ε yy /γ xy ratio with ε xx = 0. The dynamic loading was taken at T = 70 K (≈ 0.2 T m , where T m is melting point) by using Nosé-Hoover thermostat. [22,23] Under such low temperature, the diffusion process can be greatly suppressed. All the calculations are performed using the LAMMPS code. [24] For each loading step, we used the NEMD technique to calculate thermal conductivity in a given configuration. [21] 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 κ along x direction is calculated via Fourier law as κ = −J x /(∂T/∂x), where J x is the heat flux from heat source to heat sink and T is temperature. After applying heat transfer, the induced temperature gradient will range from 101 K (= 0.28 T m ) to 47 K (= 0.13 T m ), corresponding to the heat source and the heat sink (see Fig. S1 in Supplementary Information).
We first apply a simple shear deformation to the sample with strain tensor as [0, 0, γ xy ] T . Figure   2a shows the response of shear stress with shear strain (γ xy ) in a cycle. For the loading process (see black curve), the sample yields when γ xy reached 0.6%. The sample was unloaded to zerostrained state from γ xy = 1.6% (see red curve). Figure 2b shows the variation of κ during the loading/unloading loop. In elastic regime, the magnitude of κ does not change significantly ( ~ 140 to 150 W/m/K). When loading into the plastic regime, κ undergoes an abrupt drop to ~ 90 W/m/K, almost one half of the initial value. This magnitude is sustained under further stress.
When the system is unloaded, κ increases again and shows interesting hysteresis behavior during the loading/unloading cycle. This hysteresis of thermal transport could lead to the novel thermal memory device application [4,25] and is based on the underlying twin boundary dynamics.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   6 We examined the evolution of twin pattern for several typical strains (marked as (c), (d), (e), (f) in Figure 2b) upon loading/unloading, as shown in Figure 2c-2f. 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 boundaries (HTBs), 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 of twin pattern is quite stable and the twinning morphology will not undergo large changes even under the applied temperature gradient in NEMD calculations (see Fig. S2 in Supplementary Information). Upon unloading, horizontal twins will be preserved and produce permanent deformation of the sample. Vertical twins are unstable and will vanish gradually when removing the external load, which leads to an increase of κ.
These results indicate that it is very likely that the reduction of thermal conductivity results from the existence of vertical twins. To further probe the exact role that VTBs and HTBs play in thermal transfer, we examined the variation of both VTBs density (ρ VTB ) and HTBs density (ρ HTB ) in loading/unloading cycle. The units of ρ VTB and ρ HTB are a -1 . Figure 2g shows clearly that the production and annihilation of vertical twins are highly correlated with changes of the thermal conductivity, i.e., the formation of VTBs reduces thermal transport. The presence of horizontal twins does not seem to influence the heat transport and thermal conductivity is independent of the density of HTBs (Figure 2h). In addition, the dominant effects of VTBs, but not the HTBs, on the reduction of thermal conductivity can be seen in the phonon density of states. We found that the existence of VTBs can severely suppress the contribution of some medium-frequency phonons on the heat transport due to the twin boundary-phonon scattering (see Fig. S3 in Supplementary Information). To distinguish the two twin structures with different heat transfer properties, we term the systems only containing HTBs as twinning pattern 1 (TP1) and that containing VTBs (purely vertical twin boundaries or mixed twin structure) as twinning pattern 2 (TP2). For the present case, the magnitude of κ in TP2 can be lowered twice compared with TP1.
In kinetics theory of phonon transport, the thermal conductivity is proportional to the mean free path of phonons as κ = ΛCv/3, where Λ 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 kinds of phononscattering events, involving (a) phonon-phonon collision and (b) interactions of phonons and domain boundaries. The effective Λ can be given by P-P P-TB where Λ P-P and Λ P-TB are lengths of bulk phonon-phonon scattering and phonon-TBs scattering, respectively. Usually, the magnitude of Λ P-P approaches some μm [26,27] that is much larger than Λ P-TB which is restricted to some tenths nm (Λ P-P >> Λ P-TB ). Thus, Λ P-TB takes a dominant role in determination of the total mean free path Λ. For a constant Λ P-P , 1/κ should be proportional to 1/Λ P-TB (1/κ ∝ 1/Λ P-TB ). On the other hand, Λ P-TB here should be equal to the average twin boundary spacing (λ), which is inverse proportional to the density of VTBs as λ = 1/ρ VTB .
Although other studies have shown that the existence of interfaces can reduce heat transfer due to phonon scattering, [7] they did not show that this reduction can be used to flexibly tune the thermal conductivity, which is the key requirement for device applications. In this paper, we propose an alternative strategy, namely that thermal conductivity can be tuned by controlling twin patterns and twin boundary densities thought external strain field. The same conclusions will hold for ferroelectric materials where the polar 90 o twin boundaries can be changed by electric fields.

Calculation thermal conductivity by non-equilibrium molecular dynamics (NEMD) simulation
We used NEMD method to calculate thermal conductivity. [1] Specifically, the whole system are divided into 40 slabs along x direction (Fig. S1a). For calculating thermal conductivity in NEMD method, there will induce a temperature gradient on system. The range is dependent on the frequency of kinetic exchange between heat sources and sinks. In our simulation, the dynamics loading is performed at T = 70 K (~ 0.2 T m , where T m is melting point). After applying heat transfer, the induced temperature gradient will range from 101 K (= 0.28 T m ) to 47 (= 0.13 T m ). Such temperature gradient will not make thermal evolution of twinning morphologies, as two atomic configurations before and after applying heat flow shown in Fig. S2. Actually, the temperature range is dependent on the frequency of heat exchange. [1] 1 2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   15 3. Phonon spectrum analysis We further probe the role of VTBs on the reduction of thermal conductivity by analyzing phonon spectrum. Figure S3a shows the phonon frequency distribution of three absolutely different systems, i.e., the bulk one in absence of TBs (Fig. S3b), the one solely containing HTBs (Fig.   S3c) and the one primarily containing VTBs (Fig. S3d). The phonon spectrum was obtained by taking the Fourier transformation of auto-correlation function of velocity in x component. The vibrational frequency ranges from 0 THz to 40 THz. We find that for bulk system, the most probable frequencies are distributed around 20 THz, as two major peaks shown in Fig. S3b.
Thus, the phonons with the mediate-frequency and wavelength dominate the thermal transport.