Temperature oscillation in one-dimensional superlattice induced by phonon localization

Thermal transport properties and thermodynamic quantities often present anomalous behaviors in low-dimensional systems. In this paper, it is found that temperature oscillates spatially in one-dimensional harmonic and weakly anharmonic superlattice. With the increase of anharmonicity, the temperature oscillation gradually disappears and a normal temperature gradient forms. Further analysis reveals that the formation of temperature oscillation is due to the localization of high frequency phonons which cannot be thermalized. Moreover, the localized modes interact weakly with heat reservoirs, thus, their contributions to local temperature remain negligible while varying the temperatures of heat reservoirs. The oscillated temperature profile is in a good agreement with Visscher’s formula. The temperature oscillation discovered here has great potential in applications of phononic devices for heat manipulation.


Introduction
The last few years have witnessed a growing interest in regulating thermal transport by nanostructures [1,2]. Due to the subtle sizes and structures, artificial nanoscale systems manifest unique thermal properties [3,4], which have great potentials in applications of thermal dissipation in microelectronics and thermoelectric materials for solid-state refrigeration [5]. Many intriguing phenomena like anomalous thermal transport [6], thermal rectification effect [7,8], negative differential thermal resistance [9] have been discovered. These phenomena are further discussed in order to control thermal dissipation at interfaces [10,11] and design new thermal devices [12][13][14]. Recently, the wave nature of phonons has been utilized to control thermal transport [15][16][17][18][19]. A series of new structures such as phononic crystals [20], nanophononic metamaterial [21,22] have been brought up and their thermal properties are investigated by simulations and experiments. However, most of previous studies were focused on thermal conductivity and thermal current, neglecting the characteristics of temperature gradient.
The temperature profile is anomalous in some low dimensional systems. Recently, a counterintuitive phenomenon referred as negative temperature jump was found in one-dimensional (1D) interface model [23]. What is more, an overcooling effect, meaning the local temperature being lower than both heat reservoirs, were found in a multi-segment system [24]. Later, thermal siphon effect [25] was discovered in a complex network based on the FPUT-β model. The concept of negative thermal transport [26] was brought up in analogy with negative refractive index [27]. From the perspective of statistic theory, so-called non-Clausius heat transfer was discussed in ergodic and non-ergodic systems under non-equilibrium states [28]. In this case, it is important to investigate and reveal the intrinsic physics of anomalous temperature gradient.
Superlattice, an excellent candidate for thermoelectric application [29], has attracted great attention during the past few decades. Many novel phenomena have been discovered in superlattice [30][31][32], which makes it the perfect test bed to explore physical mechanism of thermal transport in novel nanostructures. Especially, due to its solvability, harmonic systems with impurities or artificial structures were investigated intensively [33][34][35][36]. Analogy to Anderson localization, the effect of phonon localization on thermal transport properties was studied [32,37,38]. We find that the negative temperature jump [24] which seems to violate the second law of thermodynamics not only exists in interfacial structure but also in superlattice. Although the counterintuitive phenomenon has been reported, the underlying mechanism still needs to be further explored.
In this paper, the anomalous temperature profile of a 1D binary superlattice model is investigated. Firstly, temperature oscillation is observed in a harmonic superlattice through non-equilibrium molecular dynamics (NEMD) simulations. The effect of anharmonicity and on-site interaction on negative temperature jump are discussed. Then, through temperature spectrum decomposition, the contributions of phonons from different frequencies to local temperature are revealed. The mechanism of phonon localization is demonstrated by conducting normal mode analysis and inverse participation ratio (IPR) calculation. The local temperature is theoretically calculated in a perturbation approximation by employing Visscher's formula [39], and compared with that by NEMD simulation. Finally, the dependence of temperature profile on temperatures of heat reservoirs is studied.

Model and method
For a general 1D lattice model, the Hamiltonian is given by where m i , p i and u i are the mass, momentum and displacement of the ith atom, respectively. V (u i−1 , u i ) is the interaction potential between two adjacent atoms, U (u i ) is the on-site potential. In this work, we consider a binary-mass lattice, where the mass of an atom can be either m a = m + Δm or m b = m − Δm.
The interaction potential takes the form of widely used FPUT-β model: , and on-site potential takes the form of widely used φ 4 model: U (u i ) = ηu 4 i . The system is coupled to Langevin heat reservoirs with temperature T L and T R at two ends, respectively. In our NEMD simulations, we set m = 1.0, Δm = 0.2, k B = 1.0, T L = 0.6, T R = 0.4, coupling strength λ = 1.0. The local temperature is defined by average kinetic energy as T i = m i u 2 i . The equations of motion are integrated by Verlet-velocity method with a time step of 0.005. T i is computed for an average of 10 10 timesteps after the system reaches steady state.
In order to get the information of temperature in frequency space, the temperature spectrum decomposition is conducted [24]. The power spectrum density (PSD) is defined as: (2) Therefore, we have the relation between PSD and local temperature as: On the other hand, if considering an isolated system with no interactions to heat reservoirs, the normal mode analysis can be used: where K is the force matrix and M is the mass matrix. ω p and e p are the eigenfrequency and eigenvector of the pth mode. For a given mode, the degree of localization can be characterized by the IPR [36]: where a n is component of the eigenvector of a given mode. N is the total number of atoms. p −1 takes the value from 1/N to 1, corresponding from fully propagating to fully localized mode.

Results
Firstly, we consider a one-dimensional harmonic superlattice. Each period contains the same amount of atom A and atom B. The temperature profiles of systems with different period lengths are shown in figure 1. It is found that the local temperature oscillates and the heavy atoms are always hotter. As the period length increases, the amplitude of temperature oscillation increases accordingly. Especially when the period length reaches 16 atoms as a quarter of system length (in figure 1(d)), the temperature of light atoms is even lower than that of right heat reservoir (heat sink). This phenomenon is referred as overcooling effect. It is also noticed that the last period of light atoms connecting directly to the right reservoir does not overcool. Particularly when the period length equals to the system length, there is just a temperature jump instead of temperature oscillation. That is because in our numerical simulations, the atoms at both ends are strongly coupled to heat reservoirs. We also confirm that the overcooling effect appears regardless of the system length. It helps us to rule out that the temperature oscillation is due to fluctuations in small systems. As a matter of fact, the temperature oscillation that increases with period length is a result of phonon localization which will be further discussed later [40].
In the above harmonic model, phonons transport ballistically and do not interact with others. After introducing anharmonicity to the system and phonon-phonon interactions are induced, the system will undergo different thermalization processes. In figure 2(a), overcooling effect gradually disappears as the anharmonicity increases and the temperature oscillation still persists with weak anharmonicity. We believe that the delocalization process is positively related to the anharmonic strength. A system with anharmonicity on the order of 0.1-1 is thermalized, and a normal temperature gradient forms. Except nonlinear interactions, the effect of on-site potential on temperature oscillation is also studied. Because nonlinear interaction only brings chaos, while the on-site potential breaks the momentum conservation of the system. As shown in figure 2(b), temperature oscillation and overcooling effect weakens quickly when η = 0.001, and diminish when η only equals to 0.01. This implies that temperature oscillation and overcooling effect are sensitive to the breaking of momentum conservation.
As thermal transport is contributed by phonons from different frequencies. To clarify the underlying mechanism of temperature oscillation and overcooling effect, it is necessary to quantify the contributions of phonons from different frequencies to the temperature profile. Considering this, we plot the power density spectrum (PSD) patterns of harmonic superlattice and FPUT-β superlattice with β = 1. The latter system is proven to show normal temperature profile (shown in figure 2(a)). In figure 3(a), the PSD patterns at light-atom sites are much smaller than those at heavy-atom sites, especially for middle and high frequency phonons. This implies that such phonons do not contribute to temperature at light-atom sites, which results in the lower local temperature at light-atom sites and leads to temperature oscillation. In figure 3(b), the cut-off frequency is up to 0.5 owing to the existence of high-order interactions. More interestingly, the PSD patterns of FPUT-β superlattice are similar to those of harmonic superlattice shown in figure 3(a) at low and middle frequency regions. However, it is obvious that more modes of FPUT-β superlattice at high frequencies above 0.4 are excited and make contributions to the local temperature due to phonon-phonon interactions. In other words, the energy of low frequency phonons can be transported to high frequency phonons, and then contributes to temperature.
To reveal the reason why some middle and high frequency phonons do not contribute to temperature at light-atom sites. We take a step further and carry out analytical approaches to get the localization information of each eigenmode. Since phonon localization may occur due to destructive interference at such system, which is derived from the wave nature of phonons [30][31][32]. We calculated the IPR of harmonic system using equation (5). As shown in figure 4(a), the IPR of last several modes are relatively large which means they are localized modes. In figure 4(b), the real space distribution of the last four modes is presented. ω 1 which has the largest IPR has only one peak at the light atom sites, ω 3 which has smaller IPR has two peaks. ω 2 and ω 4 are degenerated and have one peak whose magnitude equals to ω 3 and two small peaks at both sides. Thus, it can be concluded that the localized modes are located at light atom sites.
It also needs to be noticed that, in figure 1, temperature oscillation increases when periodic length increases and overcooling effect occurs when the periodic length equals to 16. As temperature oscillation is directly induced by phonon localization, the amplitude of temperature oscillation can be used to measure the intensity of localization. For cases shown in figure 1(a) and (b), temperature oscillation is small, which means localization is weak. This is because phonons, with wave length much large than periodic length, can not be localized and will transport through the system. This point has been verified by previous studies [41]. So, when the periodic length increases, more phonons with large wave length can be localized. As a result, localization is strengthened, temperature oscillation increases and overcooling effect occurs.
So far, we have discussed the results from numerical simulations and lattice dynamics, revealed the underlying mechanism of phonon localization. To further show the physical mechanism from theoretical point, the temperature profile of harmonic system is calculated through the formalism developed by Matsuda and Ishii [42] and improved by Visscher [39]. The formalism was first developed in order to investigate the effect of disorder on energy transport in harmonic chain, which can also be applied to our system. In nonequilibrium state, the coupling strength λ between heat reservoirs and atoms is treated as a  perturbation. The equations of motion are: and the local temperature is defined from average square velocity: where C is a dimensionless constant anḋ To solve equation (6), we shall expand u: where γ is chosen to satisfy equation (8). Then, we have For unperturbed case u α i satisfy: and we set: Then the equation (7) can be written as: Keeping terms in equation (12) up to first order in λ, it can be found that: so the local temperature of an atom in the bulk can be written as: where the value of ε has been fixed in order to keep T n equals to T L and T R under equilibrium states, u α,β i is the eigenvector of each atom, ω is the eigenfrequency, λ is the coupling strength. We consider the heat reservoir as a perturbation in equation (14) and calculate the local temperature of each atom. Equation (14) shows that the contributions from each eigenmode to the local temperature are proportional to their square amplitudes on site and at both ends. From figure 4(b), the square amplitudes of localized modes at both ends are close to zero, which means they are barely thermalized and make little contributions to the local temperature.
As shown in figure 5, the analytical temperature profile oscillates with mass alternation, which matches well with numerical simulations. Both of them are in correspondence to space distribution of localized modes (red dot line). This further verifies the underlying mechanism of phonon localization. It also should be noticed that temperatures of numerical and analytical results show some difference for atom sites adjacent to right reservoir. This is because, in numerical simulations, the coupling with heat reservoirs is strong (λ = 1) and atoms adjacent to heat reservoirs are driven to be thermalized at the temperature of heat reservoirs. While in our perturbation approximation, the local temperatures are calculated under the circumstance that we treat all atoms in the bulk. Therefore, while the summation of square amplitudes of localized modes on the right end is not zero, it has a sharp tendency toward zero. Now that anomalous temperature oscillation and underlying phonon localization mechanism are verified through different approaches. Finally, we conduct numerical simulations to investigate how temperatures of right heat reservoir affects this phenomenon. Here, the temperature of left heat reservoir is set as 0.6 (T L = 0.6). As shown in figure 6, temperature oscillation and overcooling effect occur in both equilibrium (T R = 0.6) and non-equilibrium (T R = 0.4) cases. However, when temperature of right heat reservoir is set as zero (T R = 0.0), temperature oscillates, but overcooling effect does not occur. This is because there are propagating phonons in the system, which can contribute to temperature. All in all, phonons can not interact with others in harmonic system. That is to say, propagating phonons can not 'see' temperatures contributed by localized phonons. So, propagating phonons, as well as thermal energy carried by them, are still transported from 'high temperature' to 'low temperature'. This counterintuitive phenomenon does not break the second law of thermodynamics.

Conclusion and discussion
In this paper, we investigated the temperature profile of an one-dimensional superlattice model. It is found that temperature oscillates in harmonic and weakly anharmonic superlattice. In some cases, the temperature profile exhibits overcooling effect in light atom sites. As we increase the anharmonicity of the system, the temperature oscillation gradually vanishes due to phonon-phonon interactions. The underlying phonon localization mechanism is verified by lattice dynamics analysis and a theoretical formula. It is revealed that localized modes located on light atom sites almost do not thermalize and make little contributions to local temperature. Thus, temperature oscillation is an intrinsic phenomenon which is caused by phonon localization. Our study gains a deep insight into anomalous temperature profile in nanostructures and provides guidance for designing phononic devices for heat manipulation.
As a finial remark, the phenomena of temperature oscillation and overcooling studied here may also be extended to superlattices of 2D or 3D systems [32]. Instead of studying on conceptual system, researches on real materials are expected which will provide more information on experiment observation of the phenomena reported here. To achieve this goal, materials with long phonon mean free path are needed to realize ballistic or quasi-ballistic transport, for example graphene. Besides that, low temperature is another choice. Because only low frequency phonons can be excited at low temperature due to Bose-Einstein distribution, which generally have long mean free path. Figure A1. The evolution of temperature distribution for harmonic superlattice with β equals to 10 −5 . Here, the periodic length is 16. Figure A2. The temperature distribution of steady state (T s ) for harmonic superlattice with different initial temperatures (T i ). Here, the periodic length is 16.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

Appendix A. The evolution of temperature distribution for a long simulation time
The temperature in our paper is averaged over 10 10 timesteps after the system reaches steady state. To check whether the temperature will change slowly in the weakly anharmonic superlattice, the simulation time is increased to 5 ×10 10 timesteps after the system reaches steady state. The temperature distributions are averaged over every 10 10 timesteps, which are shown in figure A1. It is found that the amplitude of temperature oscillation does not decreases. This implies that the temperature oscillation is robust. However, it needs to point out that the system will be thermalized by arbitrarily weak nonlinear interactions when its size is large enough, which has been clarified in reference [40].

Appendix B. The temperature distribution with different initial temperatures
To check whether the energy trapped by localized phonons depends on the initial condition, the temperature distributions of the same system are calculated with different initial temperatures. As shown in figure A2, we consider two cases with initial average temperatures equal to 1 and 0.25, respectively. After the systems reach steady state, the temperature distributions are similar to each other. So, the finial temperature distribution does not depend on initial temperature. This also implies that the initial temperature is contributed mainly by propagated modes, and the localized modes are not excited. If the localized modes can be excited, their energy will not be able to decay, which has been clarified in reference [24].