Effect of anisotropic interactions on the heat conduction of one-dimensional chains

One-dimensional (1D) chain models are frequently employed to investigate heat conduction in various systems. Despite their widespread use, there has been limited exploration of 1D chain models with anisotropic interactions. In this study, we examine heat conduction in a 1D chain model with orientation–position coupled interaction, namely the compressible XY model, and compare it with isotropic Fermi-Pasta-Ulam-Tsingou β (FPUT-β) systems. At low temperatures, the local temperatures in the translational and rotational degrees of freedom differ due to the difference in the contact thermal resistance in the two degrees of freedom. The system maintains orientational order, and the orientation–position-coupling effect lowers thermal conductivity in translational degrees of freedom. As the temperature rises above a transition point, the rotation of particles switches from oscillation to diffusion, leading to a solid-like to fluid-like transition in the rotational degree of freedom. The anisotropic interactions become negligible under time averaging, making the system’s properties similar to those of isotropic FPUT-β systems. Additionally, we investigate the system’s orientational structure to elucidate this transition. These findings will enhance our understanding of the dynamics of nanoscale anisotropic systems, such as the heat conduction of spin chains.

1D chain models are widely employed to investigate heat conduction in low-dimensional systems with various types of interactions [17][18][19][20].However, few studies have examined heat conduction in 1D lattices with anisotropic interactions [21,22], which are ubiquitous in nature.Anisotropy affects systems' physical properties and dynamical behaviors at different scales, including crystal growth [23], material microstructure [24], self-assembly [25], and materials performance [26].To accurately model anisotropically interacting particles, it is necessary to incorporate rotational degrees of freedom into the system.This alters the system's symmetry and enriches its phase behaviors and dynamics [27][28][29].In dynamical models of anisotropic systems, particles' rotational and translational motions are typically coupled.This orientation-position coupling is closely related to numerous research topics, such as spin-lattice coupling [30][31][32], spin-phonon coupling [33][34][35], and magnetoelastic interactions [36][37][38].Among the various spin models, the Ising model is a classic statistical physics model that describes interactions between neighboring spins [39,40].Landau and others introduced the compression effect to the Ising model to provide a more realistic description of spin systems [41][42][43].In their modification, orientation-position coupling was introduced to the Ising model.
The XY model is another widely studied spin model that describes magnetic materials' dynamics and phase transitions [44][45][46][47].In the classic XY model, a spin is represented by a unit vector that rotates continuously in a plane, allowing researchers to apply molecular dynamics simulations to study the spin systems.This article comprehensively investigates heat conduction in a 1D chain of particles with anisotropic interactions.Drawing on the form of the compressible Ising model [43,[48][49][50], we extend the FPUT-β model by incorporating a factor in the form of the XY model to account for orientation-position coupling, namely the compressible XY model.We derive analytical expressions for heat flux and other quantities related to heat conduction in the anisotropic system and use molecular dynamics simulations to compute thermal conductivity, temperature gradient, heat flux, and other properties in the heat conduction of the 1D chain in both the compressible XY model and FPUT-β model.Our findings indicate that translational thermal conductivity exhibits different behaviors depending on the anisotropy.We explain the physical mechanisms underlying our results and compare them with previous studies on heat conduction in 1D chain models.Additionally, we provide a detailed characterization of the system's orientational structure at different temperatures to support our explanation.
This article is structured as follows.Section 2 describes the new compressible XY model and the detailed parameter settings for our simulations.In section 3, we sequentially discuss the temperature profile, heat flux, and thermal conductivity of the new model 1D chains.The system's orientational structure is then carefully analyzed, proving that there is a transition from order to disorder in the rotational degree of freedom.Finally, we summarize and discuss our findings and provide an outlook for future research.

Model and simulation methods
Drawing on the compressible Ising model [41], we extend the FPUT-β model to consider the orientation-position coupling.The pair interaction Φ i,i−1 between two neighboring particles i and i − 1 is In equation ( 1), the position-dependent factor U takes the form of FPUT-β potential [51], where x i is the displacement of the ith particle from its lattice position.The orientation-dependent factor f is in the form of the classic XY model [52], and θ i is the orientation of (the spin of) the ith particle.Analogous to the compressible Ising model, the model defined by equation ( 1) can be called a compressible XY model.Note that in the equation for the rotational degrees of freedom, equation ( 2), the constant α can enhance the strength of interactions and alter the temperature scale.To explain the role of α, we investigate the interaction potential Φ i,i−1 between particles within the system, expand upon it and obtain From the expansion, it can be seen that the constant α contributes as a factor to the interaction potential of the system, which strengthens the interaction strength and changes the temperature scale in the rotational degree of freedom.In real physical systems, the interaction strengths in these two degrees of freedom can be different.For example, when a ferromagnetic material is heated above the Curie temperature, it loses its ferromagnetic property (spontaneous demagnetization phase transition) without disrupting its lattice structure.It does not qualitatively change the nature of the system.In order to observe the effects of anisotropic interactions better, and reduce the numerical error, we set α = 1/2.Note that introducing a constant α greater than 0 results in an enlargement of the temperature scale within the rotational degrees of freedom.The Hamiltonian of the 1D chain of anisotropic interacting particles reads The first two terms under summation on the RHS of equation ( 4) are the translational and rotational kinetic energies, where v i and w i are the velocity and angular velocity of the ith particle, respectively.The first and last particles (i = 1, N) in the chain are connected to the heat baths of high and low temperatures [53], respectively.We use the fixed boundary conditions, i.e. x 0 = x N+1 = 0 and θ 0 = θ N+1 = 0.The equations of motion of the system are where i = 1, 2, . . ., N and j = 1, N. ξ j and η j are Gaussian white noise satisfying where ⟨• • • ⟩ represents the ensemble average, and λ 1 and λ 2 are the coupling strengths of the translational and rotational degrees of freedom of the system to the heat bath, respectively.As known in the context of molecular dynamics simulations [54,55], the precise values of λ 1 and λ 2 do not alter the steady-state properties of the system.We set λ 1 = λ 2 = 1 in the Langevin thermostat.In this work, we assume a dimensionless unit mass m = 1 and rotational inertia I = 1 for simplicity.k B is the Boltzmann constant, and in the simulation, we set k B = 1.Thus the temperature takes the dimension of energy.T j represents the temperature of the heat bath, with T 1 > T N .
For finite sizes or under low temperatures, the temperature profile in the system may have significant boundary jumps caused by the conduct thermal resistance, which plays a crucial role in determining thermal conductivity [56].To avoid boundary jumps caused by the finite-size effects, a system size of N ⩾ 1000 is typically chosen [57].In our simulation, after testing different N, we find N = 1000 is large enough to avoid a finite-size effect.So, we use molecular dynamics methods [58,59] for simulation with a system size of N = 1000 and solve the stochastic equations of motion using the second-order BAOAB algorithm [60] with a time step of ∆t = 0.005.We set the following parameters for the system: the average temperature was denoted as T, with a temperature difference of 2∆T = 0.6T between the two heat baths at the left and right ends of the system, denoted as T + = T + ∆T and T − = T − ∆T, respectively.In the system's initial state, we set x i = 0 and θ i = 0 for i = 1, 2 . . .N.Then, the system evolves with the heat baths to reach a steady state.Each simulation runs for at least 5 × 10 7 time units to ensure the steady state is reached.After the system reached a steady state, we continued the evolution for an additional 5 × 10 7 time units to collect data for analyzing relevant physical quantities during this interval.

Results and analysis
Abnormal heat conduction is found in the nanoscale low-dimensional lattice [15,56,61].When the rotational degree of freedom is introduced, the system's dynamic behavior becomes more complicated due to the orientation-position coupling.The heat conduction properties require reexamination.To compare heat conduction with isotropic models, we need to extend the form of relevant physical quantities, such as temperature, heat flux, and thermal conductivity, to include the new rotational degrees of freedom.

Temperature profile
As in the FPUT-β model, we can define local temperature by the time average of kinetic energy for both translational and rotational degrees of freedom.The local temperatures in translational and rotational degrees of freedom of the ith particle are which can be calculated directly from the simulation data.Figures 1(a)-(c) show the temperature profiles of the system at different temperatures.T β represents the temperature in the FPUT-β model, while T t and T r are the temperatures on translational and rotational degrees of freedom in the compressible XY model, respectively.For this finite-size non-equilibrium system, the boundary temperature jump caused by contact thermal resistance exists in all temperatures.At low temperatures, the ratio of boundary temperature jump to bath temperature difference is huge (figures 1(a) and (b)); at high temperatures, this ratio is small (figure 1(c)).For a given temperature, the temperature jumps in different degrees of freedom are not equal (figures 1(b) and (c)).This leads to misalignment of the temperature profiles of different degrees of freedom in the bulk region; that is, the local translational temperature T t is not equal to the rotational temperature T r .This result is inconsistent with the equipartition theorem, which holds strictly in the equilibrium systems [62], while our system is in the non-equilibrium steady state.The temperature gradients ∇T t and ∇T r can be well established in the central linear region of [N/5, 4N/5].In the following, this linear central region is used to compute the temperature gradient and obtain an effective thermal conductivity [63], as in [64].At low temperatures, the ratio ∇T t /∇T r deviates significantly from 1, and as the temperature rises, it tends to 1.

Heat flux
The instantaneous local heat flux J i is derived from the discrete version of the continuity equation where the local energy density h i [15,22] is Substituting equation ( 9) into equation ( 8), the local heat flux J i can be obtained: where ) in equations ( 2) and (3), and is the force (torque) exerted by the (i − 1)th particle on the ith particle.The local average heat flux is:

Ji = ⟨v
where ⟨• • • ⟩ denotes the time averaging, and are the heat flux of the ith particle in the translational and rotational degrees of freedom, respectively.Hereafter, we denote Ji as J i because only time-averaged quantities are considered.After the thermalization of the system, the local average total heat flux is equal everywhere, i.e.
Like the local temperature, the heat flux can be calculated directly from the simulation data.The heat flux's dependence on temperature is shown in figure 2. The total heat flux J tot coincides with J β at both low and high-temperature limits (figure 2(a)), where J β is the heat flux in the corresponding FPUT-β model, i.e.
and that J t i ⩽ J β i always holds.At low temperature, θ i ∼ 0 and f ∼ 1  2 , so J t < J β ; at high temperature, θ i is distributed randomly and uniformly in [0, 2π] and f ∼ 1, so J t ≈ J β .In figure 2(b) J r /J t ∼ 1 at low temperatures indicates that the heat flux on both degrees of freedom are of the same order; near T c = 10 −1.4 (purple squares in figure 2), J r /J t begin to decrease against temperature; at high temperature J r /J t → 0 means there is almost no heat flux in the rotational degree of freedom.The similar behavior in translational/rotational degree of freedom can be understood theoretically.Expanding the orientation factor f i,i−1 into a Taylor series, and since the fourth-order term becomes a higher-order negligible quantity at low temperatures, only the quadratic term is retained.Thus, we obtain in the low-temperature FPUT-β model, the quadratic terms play a major role [9]).In this case, f i,i−1 and U i,i−1 have the same mathematical structure, i.e. at low temperatures, both translational and rotational degrees of freedom follow the same laws of physical motion.Therefore, at low temperatures, J r /J t ∼ 1.

Thermal conductivity
Under small temperature differences (here ∆T = 0.3T), Fourier's law [65] always holds for both translational and rotational degrees of freedom where J l , K l and ∇T l are heat flux, thermal conductivity, and temperature gradient, respectively.The subscript l = β, t, r, tot means that the quantity is in the FPUT-β model, in translational/rotational degrees of freedom, or for the total effect of the system (translational + orientational), respectively.Thermal conductivity can be calculated from the heat flux and the temperature gradient.Figure 3(a) shows the thermal conductivity in two degrees of freedom of the compressible XY model against temperature, and compared to the FPUT-β model.In our simulations of the FPUT-β model, the behavior of the thermal conductivity at the high and low temperatures replicates the results in [56]: K β at high and low temperatures is consistent with previous results: the low-temperature limit satisfies K β ∼ T −1 , while the high-temperature limit is K β ∼ T 1 4 .For the compressible XY model, the trend of translational thermal conductivity K t at the high and low-temperature limits is consistent with the FPUT-β model, but K t is always smaller than K β .This is due to the coupling between the translational and rotational degrees of freedom.In equations ( 1) and ( 3), the f factor in the rotational degree of freedom is equivalent to changing the energy scale of the translational degree of freedom, thereby affecting the translational thermal conductivity.
The trend of K r is even more noteworthy.Figure 3(b) shows the dependence of K r /K t on T. K r and K t are of the same order of magnitude at low temperatures.As the temperature increases, K r significantly decreases and becomes much smaller than K t (figure 3(a)).As T increases, this ratio rapidly decreases to the order of 10 −3 .The purple squares in figure 3(a) at T c = 10 −1. 4 show a sudden drop of K r .The purple squares in figures 3(b) and (c) corresponding to the plunge of K r /K t are also at T c = 10 −1. 4 .Figure 2(a) shows that the trends of J t , J tot , and J β are similar, but J r does not always increase with temperature and has a decreasing process.The decrease in J r directly is because of the decrease in K r .In figure 2 purple squares are also located at T c = 10 −1. 4 and near the inflection point of J r .Although K is calculated from ∇T and J, it actually reflects the system's structural properties and determines J's behavior.It is reasonable to speculate that a transition in the structure and the system's dynamics in the rotational degree of freedom may have occurred at near T c = 10 −1.4 .

Orientational structure
Note that in the simulation in this study, the system is in a non-equilibrium steady state and has a finite size.Though it is subtle and difficult to talk about 'phase transitions' strictly for the non-equilibrium system in finite size, it is still possible to examine the order of the system and its dynamic transition by observing the configuration of the system directly.In the 1D FPUT-β chain, the positions of the particles are restricted within the lattice structure by the fixed boundary conditions.Therefore, changes in the system's structure are more likely attributable to the changes in the structure of the particles' orientation.
The inflection point of heat flux and thermal conductivity in the rotational degree of freedom in the low-temperature range implies that the system undergoes a structural transition in the rotational degree of freedom at near T c = 10 −1.4 .To verify this, we first observe the rotation trajectories of the intermediate particle (i = 500) in the chain at three typical temperatures, which are extremely low temperature T = 0.001, low temperature T and a temperature above the transition T 0.1 > T c (figure 4).At T the spin oscillates around the equilibrium position; as the temperature increases, at T = 0.01, the oscillatory motion of the higher fluctuations but the orientational order remains; when the temperature continues to rise over the transition point to T = 0.1 > T c , the spin no longer oscillates around the equilibrium position but starts to rotate randomly, which is similar to diffusive motion.Therefore, system's the rotational motion at low temperatures is similar to the lattice vibration in the solid and turns into the diffusion behavior in the fluid when T > T c .This solid-like to fluid-like transition in the rotational degree of freedom is similar to the melting transition.
To further investigate the dynamic process of particles' orientation, we calculate the mean square displacement (MSD) [66] in the rotational degree of freedom for the intermediate particle: where ⟨• • • ⟩ denotes taking the average over t 0 .The corresponding rotational diffusion coefficient is Figures 5(a)-(c) represent the rotational MSD of particles at T = 0.001, 0.01, and 0.1, respectively.In figures 5(a) and (b), the MSD exhibits a periodic oscillatory behavior.This indicates that the system is crystalline-like in the rotational degrees of freedom.In figure 5(c), the particles' orientation undergoes rotational diffusion at T = 0.1, where the rotational melting happens.Figures 5(d) and (e) exhibit the relationship between the rotational diffusion coefficient D and the temperature T. For T < T c = 10 −1. 4  (purple square in figure 5(e)), D ∼ 0, the system's orientational structure is solid-like.For T > T c , D increases sharply against T, and the orientational structure is fluid-like.It is because of this transition that J r and K r decrease dramatically, as shown in figures 2 and 3.
The solid-like to fluid-like transition in the orientational structure is an order/disorder transition.The spatial correlation C i of the particles' orientation helps study the system's orientational order, i.e. it quantifies the degree of orientation consistency within the system.It can be calculated as: where θ is the unit vector of the particle's orientation, and ⟨• • • ⟩ represents averaging over time and j. 0 ⩽ C i ⩽ 1. C i = 0 indicates that the system is completely disordered in the orientational space, while C i = 1 corresponds to a completely orientational order.Figure 6 displays the snapshots and correlations of particles' orientation within the system averaging over a short time.The snapshots in figures 6(a) and (c) are colored by the particles' orientation in the HUE cyclic color scheme.At T = 0.001, the color changes continuously and smoothly from left to right (figure 6(a)).The colors are mainly pink and red, indicating θ i ∼ 0, which is determined by the fixed boundary conditions.The particle orientation in the system is highly ordered.This can also be seen from the spatial correlation of particles' orientation.The black C i curve in figure 6(d) decreases slightly over an extensive spatial range and stabilizes near i = 300.Particles' orientation within the system exhibits a high spatial correlation, indicating that most particles' orientation vectors align with high direction synchronization.As the temperature increases, at T = 0.01, the color changes continuously but less smoothly than at T = 0.001 (figure 6(b)).In  addition to pink and red, green and blue also appear somewhere.This indicates that the system's orientation can only be maintained within a limited spatial range.The system's orientational order is decreasing, which is supported by the blue C i curve decreasing quickly to a small but nonzero value.When T = 0.1, the colors alternate rapidly and randomly (figure 6(c)).The red C i curve drops to zero dramatically.No orientational order is exhibited at that temperature, and the system's orientation becomes chaotic and almost completely disordered.
Here it is instructive to look at the probability distribution of particle orientations in the system at different temperatures (figure 7).It helps to understand the dynamics and the orientational structure analyzed above.At T = 0.001, most particles' orientation is around 0. This is a characteristic of crystals.All particles' orientations slightly oscillate around 0. As the temperature increases to T = 0.01, the curve of the distribution of particles' orientation becomes broader and lower.When T = 0.1, the curve becomes flat.These probability distributions also explain whether the system's orientation is ordered.This conclusion is consistent with that shown in figure 6.It is worth pointing out that, a uniform distribution at T = 0.1 means that the long-term average of any physical quantity will not contain a particular orientation dependence.The system will behave as if it were isotropic.

Summary and discussion
The heat conduction of the 1D compressible XY model is studied in this article.This model has different dynamics in the translational and rotational degrees of freedom, so the contact thermal resistance differs in these degrees.This directly leads to the unequal local temperatures on the two degrees of freedom, and the equipartition theorem of energy is no longer valid.The heat conduction property in the translational degree of freedom is similar to that of a 1D FPUT-β chain.However, due to the orientation-position coupling, the orientation factor f in the interaction equation ( 1) is equivalent to changing the energy scale of U with temperature.Then in the low-temperature range, the heat flux in the translational degree of freedom of the compressible XY model (J t ) is smaller than that in the FPUT-β (J β ).In the high-temperature range, since the particle orientation is completely disordered, f → 1, the two heat fluxes are equal.
Even more interesting is that the heat conduction properties in the orientation degrees of freedom differ significantly from those in the translational degrees.We find that at near T c = 10 −1. 4 there is a transition in the dynamics of the rotational degree of freedom.When T < T c , the motion in the rotational degree of freedom is the vibration of the spin around the equilibrium position, the system's orientation structure is solid-like, and the particles' orientation is highly ordered.When T > T c , the orientation structure is completely disordered, and the particles do rotational diffusion in rotational degrees of freedom.For non-equilibrium systems, we avoid using the theory of equilibrium state statistical physics and confirm this transition with various methods, such as the system's orientation configurations, the rotational trajectories, and the orientation correlation function.
It is important that our model is simplified and idealized and does not directly reflect the properties of real materials.However, in the real world, particles' interaction is usually anisotropic.Introducing rotational degrees of freedom to the system's interactions is a more realistic approach.For example, spintronic systems [67] in magnetic materials have interactions between spins and lattices vibrations, leading to coupling between spins and lattices.Some superconductors [68], superconducting electron pairs, interact with lattice vibrations, leading to coupling between spins and phonons.To our knowledge, there are few studies on the effect of the anisotropic interactions in heat conduction in 1D chains [21].The present work preliminarily explores the effect of position-orientation coupling on the heat conduction of 1D chains.It leads to deeper insights into the dynamics and the phase of the anisotropic interacting systems.It sheds light on further studies of the potential applications in the low-dimensional heat conduction of spin chain systems.

Figure 3 .
Figure 3. (a) The dependence of the thermal conductivities K β (black), Kt (blue) and Kr (red) on T. (b) The dependence of Kr/Kt on T. The inset (c) is the enlarged image of the low-temperature range 10 −1.7 ⩽ T ⩽ 10 −1.3 in (b).The purple squares in (a)-(c) are all located at Tc = 10 −1.4 .

Figure 5 .
Figure 5. Rotational MSD at (a) T = 0.001, (b) T = 0.01 and (c) T = 0.1.(d) The dependence of the rotational diffusion coefficient D on T. The inset (e) is the enlarged image of the low-temperature range 10 −1.7 ⩽ T ⩽ 10 −1.3 in (d).The purple square is located at Tc = 10 −1.4 .

Figure 6 .
Figure 6.A snapshot of the orientation of particles after the system reaching a steady state at (a) T = 0.001, (b) T = 0.01 and (c) T = 0.1.(d) Correlation of particles' orientation at corresponding temperatures in (a)-(c).

Figure 7 .
Figure 7.The probability distribution of particles' orientation after the system reaches a steady state at T = 0.001, 0.01 and 0.1.