Thermalization in Asymmetric Harmonic Chains

TL;DR

This paper studies how asymmetry in particle interactions affects thermalization in solids, showing that asymmetry enhances relaxation dynamics differently than nonlinearity.

Contribution

The study introduces a one-dimensional asymmetric harmonic model to isolate the effect of asymmetry on thermalization.

Findings

Thermalization time follows a power-law with an exponent larger than the inverse-square law in the thermodynamic limit.

Matthiessen’s rule accurately estimates thermalization time when asymmetry and nonlinearity are combined.

Asymmetry enhances higher-order effects and governs relaxation dynamics distinct from nonlinearity.

Abstract

The symmetry of the interparticle interaction potential (IIP) plays a critical role in determining the thermodynamic and transport properties of solids. This study investigates the isolated effect of IIP asymmetry on thermalization. Asymmetry and nonlinearity are typically intertwined. To isolate the effect of asymmetry, we introduce a one-dimensional asymmetric harmonic (AH) model whose IIP possesses asymmetry but no nonlinearity, evidenced by energy-independent vibrational frequencies. Extensive numerical simulations confirm a power-law relationship between thermalization time (Teq) and perturbation strength for the AH chain, revealing an exponent larger than the previously observed inverse-square law in the thermodynamic limit. Upon adding symmetric quartic nonlinearity into the AH model, we systematically study thermalization under combined asymmetry and nonlinearity. Matthiessen’s rule provides a good estimate of Teq in this case. Our results demonstrate that asymmetry plays a distinct role in enhancing higher-order effects and governing relaxation dynamics.

Figures

Full text

1. Introduction

In the 1950s, Fermi, in collaboration with Pasta, Ulam, and Tsingou (FPUT), conducted the first numerical experiments to test the ergodic hypothesis using a simple mechanical system of springs and masses [1]. Their pioneering work unexpectedly revealed that the system, far from equilibrium, did not evolve into the anticipated thermalized state but instead returned to a state of near non-equilibrium, a phenomenon now known as the FPUT recurrences. This finding sparked extensive research aimed at explaining and understanding this phenomenon through longer simulations and larger system sizes in various nonlinear chains [2,3,4,5,6,7,8,9,10,11].

In recent years, the application of wave turbulence theory has significantly advanced the understanding of this problem [12,13,14,15,16,17]. Extensive numerical simulations of various models have revealed that, in the thermodynamic limit, there exists a universal scaling law for the thermalization behavior of near-integrable systems: thermalization time (Teq) is inversely proportional to the square of perturbation strength [18,19,20,21,22,23,24,25]. Recent studies indicate that this inverse-square law holds true even in high-dimensional lattice systems, unaffected by lattice structures, interaction potentials, or whether the lattice is being ordered or not [26]. To observe this universal law, a suitable reference integrable system must be chosen to define the perturbation strength, ensuring an accurate description of the system’s ability to thermalize. For example, a general nonlinear monatomic chain can be considered as a perturbation of the Toda model, with the perturbation strength defined as the distance from the Toda system [18]. When the Toda integrability is broken, such as in diatomic chains [21,22], mass-disordered chains [23,24], chains with on-site potentials [14], or high-dimensional systems [26], the system is considered a perturbation of the harmonic one, with perturbation strength defined relative to the harmonic reference point.

However, deviations from this universal law, such as steeper slopes, have been observed in chains with cubic or quintic nonlinearity [19,20]. Two primary explanations for this deviation exist: one suggests a finite-size effect (discreteness) [20], while the other points to higher-order effects [19]. Notably, both third- and fifth-order nonlinearities involve asymmetric interaction potentials, meaning that the amplitude of forces corresponding to the same displacement is not equal in the tension and compression states.

Asymmetric interparticle interaction potentials (IIPs) play a critical role in lattice models, influencing thermal expansion effects, which symmetric potentials cannot produce [27]. Furthermore, IIP asymmetry significantly affects transport properties [28]. For instance, in a 1D momentum-conserving system, the bulk viscosity of a system with symmetric IIP remains finite in the thermodynamic limit, while it diverges for an asymmetric IIP [29,30,31,32]. When the IIP is symmetric, the heat conductivity ( $[eqn]$ ) scales with system size as $[eqn]$ [32,33], although mode coupling theory [34,35], Peierls–Boltzmann kinetic theory [36], and wave turbulence theory [37] predict $[eqn]$ , which is supported by recent ultra-large-scale nonequilibrium simulations [38]. In contrast, for asymmetric IIP, $[eqn]$ [29,32,33,39,40,41]. More surprisingly, normal heat conduction (i.e., $[eqn]$ independent of system size in the thermodynamic limit) has been observed in chains with asymmetric IIP in the near-integrable regime [42,43,44,45], challenging conventional views. While some researchers attribute this to finite-size effects [46,47], it has been shown that systems with asymmetric IIP exhibit a larger kinetic region [44,45,48,49], suggesting that asymmetric IIP may lead to diffusive kinetic behaviors, whereas symmetric IIP does not [50].

This raises the question of whether or not asymmetric IIP plays a distinct role in thermalization. In the models discussed above, such as those with odd-order nonlinearity, the effects of asymmetry and nonlinearity are intertwined, making it challenging to separate their contributions. To address this, we introduce an asymmetric harmonic (AH) model, which is purely asymmetric and not nonlinear (since nonlinearity typically involves the frequency of the motion of particles being dependent on the amplitude or, equivalently, the input energy [51]), to investigate the influence of pure asymmetry on thermalization. Subsequently, we introduce quartic nonlinearity into the AH model to explore the thermalization behavior when asymmetry and nonlinearity are interwoven. In the following sections, we first present the models and methods in Section 2, followed by the numerical results in Section 3, and conclude with a summary and discussion in Section 4.

2. Models and Method

We consider a homogeneous lattice consisting of N particles of unit mass, with the Hamiltonian given by

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the momentum and displacement from the equilibrium position of the j-th particle, respectively, and V is the nearest-neighbor interaction potential. To investigate the effects of asymmetry and nonlinearity on thermalization, we consider two types of interaction potentials. The first is the AH potential [43], defined as

$$[eqn]$$

where r is a free parameter that controls the degree of asymmetry. Actually, Equation (2) can be rewritten as

$$[eqn]$$

where $[eqn]$ is the sign function. It is clearly shown that r also is the perturbation strength relative to the harmonic (reference integrable) system. The dynamics of AH model is independent of the system’s energy (temperature) and depend solely on the parameters r and N [43]. That is, the AH model does not have nonlinearity in the sense of frequency-energy dependence [51]. The nonintegrability of the AH model originates from the nonanalyticity of the potential $[eqn]$ at $[eqn]$ . To be precise, from the perspective of the force-displacement relationship, the AH model is a piecewise linear model, which exhibits a nonlinear relationship overall. The AH potential is similar to the broken linear potential studied in the original FPUT work [1], which corresponds to a piecewise-smooth dynamical system [52]. In practice, non-smooth dynamics arising from switches, impacts, sliding, and abrupt changes are common in various fields of physics, biology, and engineering [53].

The AH model has the advantage that its dynamics are independent of the system’s energy (temperature) and depend solely on the parameters r and N [43]. We next introduce a symmetric fourth-order nonlinearity into the AH potential, yielding the AH- $[eqn]$ potential

$$[eqn]$$

where $[eqn]$ is a positive parameter controlling the strength of nonlinearity. Note that when $[eqn]$ , the AH- $[eqn]$ potential reduces to the FPUT- $[eqn]$ potential. The dimensionless parameter $[eqn]$ governs the strength of nonlinearity for the AH- $[eqn]$ model, where $[eqn]$ is the energy density, with E being the total energy of the system. For simplicity, we omit the tilde in subsequent expressions.

In this work, we consider fixed boundary conditions, i.e., $[eqn]$ . The normal modes of the chain are defined by

$$[eqn]$$

where $[eqn]$ are the mode indices. The frequency $[eqn]$ and energy $[eqn]$ of the k-th mode are given by

$$[eqn]$$

To each mode k, we associate a phase $[eqn]$ defined by

$$[eqn]$$

In the thermalized state, energy equipartition is achieved, meaning

$$[eqn]$$

where $[eqn]$ represents the time average of $[eqn]$ over the time interval $[eqn]$ , with $[eqn]$ controlling the time average window, i.e.,

$$[eqn]$$

For numerical simulations, we use $[eqn]$ , which accelerates the calculations and reduces memory effects from the initial state, as noted in Ref. [54].

To measure how close the system is to thermal equilibrium, we introduce a modified version of the normalized effective relative number of degrees of freedom [55], denoted as $[eqn]$ , which is sensitive to the early energy growth of high-frequency modes ( $[eqn]$ ), as

$$[eqn]$$

where $[eqn]$ is the spectral entropy, defined by

$$[eqn]$$

with

$$[eqn]$$

As the system approaches thermal equilibrium, $[eqn]$ saturates at the value 1.

In our numerical simulations, the equations of motion are integrated using the eighth-order Yoshida algorithm [56], with a typical time step of $[eqn]$ . The relative error in energy conservation is less than 10^−5^, and further reducing the time step by an order of magnitude does not yield significant differences. To suppress fluctuations, we average over 120 random choices of initial phases uniformly distributed in $[eqn]$ . Energy is initially distributed among 10% of the lowest-frequency modes ( $[eqn]$ ) in all simulations, and we have verified that varying the percentage of excited modes does not produce qualitative differences. In this study, the energy density is kept constant at 10^−3^, i.e., $[eqn]$ .

3. Numerical Results

Figure 1 presents the numerical results of $[eqn]$ as a function of $[eqn]$ at various selected times for the AH chain with different perturbation strengths r, and with fixed system size $[eqn]$ . The results show that the energy initially concentrated in the excited modes gradually redistributes to the other modes over time, with the energy of the remaining modes increasing continuously. This behavior contrasts with that observed in the FPUT model [7,54,57] and the perturbed Toda model [18], where $[eqn]$ maintains an exponentially distributed profile over a large initial time scale, corresponding to the so-called metastable state. Due to the non-smoothness of the AH model at $[eqn]$ , it lacks Toda integrability, and thus no metastable state is observed. As shown in Figure 1a–c, the system reaches equipartition more rapidly as r increases.

To observe the overall thermalization dynamics and determine the thermalization time, we examine the evolution of $[eqn]$ , as defined in Equation (10). Figure 2 presents the numerical results for the AH chain with different values of r. It is evident that, over a sufficiently large time scale, all $[eqn]$ values increase from 0 to 1, following very similar sigmoidal profiles, indicating that energy equipartition is eventually achieved. Additionally, the time required to reach the equipartition state increases as r decreases.

Physically, the thermalization time Teq is defined as the time at which $[eqn]$ reaches the threshold value of 1. However, for practical purposes, we are generally concerned with the scaling behavior of Teq, rather than its exact value. To reduce computational cost, we define Teq as the time when $[eqn]$ reaches the threshold value of 0.5. While this is an arbitrary choice, it does not affect the scaling behavior of Teq [3]. As shown in Figure 2b, the sigmoidal profiles in Figure 2a can be made to overlap completely upon suitable shifts, confirming that the specific threshold value does not influence the scaling law of Teq. In the following, we will explore how Teq depends on r and N for the AH chain, and on r and $[eqn]$ for the AH- $[eqn]$ chain.

In Figure 3a, we show the dependence of Teq on r for the AH chain with various values of N, on a log–log scale. For the range of r explored, Teq becomes nearly independent of N as N increases further. The numerical results suggest that $[eqn]$ is sufficiently large for the thermodynamic limit to be effectively reached. It is observed that Teq versus r follows a power-law behavior,

$$[eqn]$$

Figure 3b,c show, respectively, the dependence of the slope $[eqn]$ and the intercept of the linear fitting for the data in Figure 3a on N. It can be seen that $[eqn]$ quickly saturates at −2.65, while the intercept stabilizes at 1.35. From this, we can estimate the thermalization time for the AH model in the thermodynamic limit as

$$[eqn]$$

which deviates from the previously observed inverse-square law [18,19,20,21,22,23,24]. This phenomenon of a steeper slope is also observed in models with odd-order nonlinearity (e.g., the FPUT- $[eqn]$ chain) [19,20]. The mechanism behind this deviation remains unclear. Ref. [20] suggests that the deviation is due to finite-size effects (discreteness), while Ref. [19] shows that the results for different system sizes nearly coincide within the parameter range studied, yet the steeper slope persists. Additionally, it has been pointed out that asymmetric IIP leads to asymmetric spectral peaks of modes (i.e., a higher-order effect), which is considered the root cause of the deviation. Despite the AH model being purely asymmetric and devoid of nonlinearity, it still exhibits this deviation, further confirming that asymmetry results in a steeper slope. In other words, the asymmetry of the IIP enhances the contribution of higher-order effects. Next, we will investigate the role of nonlinearity, exemplified by the AH- $[eqn]$ model [see again Equation (4)].

Figure 4 presents the dependence of Teq on r in log–log scale for the AH- $[eqn]$ chain at fixed system size $[eqn]$ , under various values of $[eqn]$ . As r decreases, Teq approaches a saturation value, and this value declines with increasing $[eqn]$ . This behavior indicates that, in the small-r regime, thermalization is predominantly governed by the quartic nonlinearity, whereas for large r, the asymmetry becomes the dominant factor. These observations suggest a competition between the asymmetric potential and the quartic nonlinearity in driving thermalization. Assuming their contributions are independent, as proposed in Ref. [19], Matthiessen’s rule (MR) [58] can be applied to estimate the total thermalization time, which is,

$$[eqn]$$

where $[eqn]$ is given by Equation (14) and $[eqn]$ in the thermodynamic limit [19,20]. Based on Ref. [19], we adopt the empirical estimate

$$[eqn]$$

Combining Equations (14)–(16) yields an explicit expression for Teq in the AH- $[eqn]$ chain

$$[eqn]$$

The solid lines in Figure 4 correspond to the prediction of Equation (17), showing good agreement with numerical data. However, the deviation increases as $[eqn]$ decreases, which is attributed to finite-size effects. Prior works [12,13] have shown that exact four-wave resonances, responsible for $[eqn]$ , require large system sizes ( $[eqn]$ 163,264). Nonetheless, due to nonlinearity-induced spectral broadening, quasi-resonant four-wave interactions can still occur in finite systems [59,60], especially at stronger nonlinearities. Thus, the $[eqn]$ scaling emerges only for either large systems under weak nonlinearity or moderate systems under stronger nonlinearity [19], explaining the discrepancy in the small- $[eqn]$ regime and its disappearance as $[eqn]$ increases. Furthermore, Equation (17) implies that Teq decreases monotonically with both r and $[eqn]$ since $[eqn]$ and $[eqn]$ for fixed $[eqn]$ and r, respectively. However, the actual system behavior may exhibit deviations from this monotonic trend due to intertwined interactions.

Figure 5 shows the thermalization time Teq as a function of $[eqn]$ on a log–log scale for the AH- $[eqn]$ chain at fixed system size $[eqn]$ , with $[eqn]$ (red circles) and $[eqn]$ (blue triangles). The dashed-dotted lines represent theoretical predictions from Equation (17). While the general trend of the theoretical curves aligns with the simulation data, notable deviations emerge in the intermediate $[eqn]$ regime. Specifically, Teq exhibits a non-monotonic dependence on $[eqn]$ , i.e., it initially increases with increasing nonlinearity, then decreases, contrary to the common expectation that stronger nonlinearity facilitates thermalization. This discrepancy raises questions about the applicability of MR in this context. In condensed matter theory, MR assumes that scattering mechanisms are independent. However, the fourth-order nonlinearity appears to suppress the effective asymmetry of the IIP, thereby invalidating this assumption. As discussed in Ref. [61], in the limit $[eqn]$ , the potential becomes effectively symmetric and the effective asymmetry degree $[eqn]$ approaches zero.

To quantify this effect, we adopt the method in Ref. [61] to compute $[eqn]$ for the AH- $[eqn]$ model. As shown in the inset of Figure 5, $[eqn]$ decreases approximately exponentially with increasing $[eqn]$ . Based on this observation, we introduce a correction $[eqn]$ , where $[eqn]$ is a fitting parameter, and modify Equation (17) accordingly

$$[eqn]$$

The solid lines in Figure 5 represent fits using Equation (18). These curves closely follow the numerical data, validating the correction to a certain extent. Notably, the fitted values of $[eqn]$ (e.g., $[eqn]$ for $[eqn]$ ; $[eqn]$ for $[eqn]$ ) differ from those derived directly from $[eqn]$ (3.35) and show dependence on r, while the original $[eqn]$ does not (see the inset in Figure 5). This discrepancy may arise from two approximations: the exponential form assumed for $[eqn]$ and the application of MR under the implicit assumption that $[eqn]$ and $[eqn]$ are independent. The precise origin of the difference remains an open question that requires further investigation.

Nevertheless, the improved agreement between theory and simulation suggests that accounting for the reduction of effective asymmetry due to nonlinearity is crucial. For fixed r, increasing $[eqn]$ suppresses the asymmetry of the IIP, leading to the observed non-monotonic behavior of Teq. It is also worth noting that, for fixed $[eqn]$ , variations in r could influence the effective strength of nonlinearity through changes in the distribution of relative displacements and energy partition. However, Figure 4 shows that this effect is minimal, as the numerical data remain consistent with the uncorrected MR prediction.

4. Conclusions

In summary, we have investigated the thermalization behavior of a 1D lattice system, focusing on the role of interaction asymmetry and nonlinearity. We first introduced the AH model, where particles interact through a purely asymmetric, non-smooth potential without explicit nonlinearity, to isolate the effects of IIP asymmetry. Numerical results show that, unlike the perturbed Toda chain, the AH model exhibits no metastable state (see Figure 1), as the non-smoothness at $[eqn]$ breaks Toda integrability. Thus, the AH model is best regarded as a perturbation of the harmonic chain, with the degree of asymmetry serving as the perturbation strength.

In the thermodynamic limit, we find that the thermalization time Teq follows a power-law dependence on the perturbation strength, though with an exponent larger than the universal value reported in earlier studies (see Figure 3). This result is qualitatively consistent with findings for the FPUT- $[eqn]$ (cubic) and quintic nonlinear chains [19,20]. Our results—as well as those in Ref. [19]—indicate that finite-size effects are negligible in models with asymmetric IIPs, with rapid convergence observed as system size increases. Therefore, the enhanced slope is more likely due to higher-order contributions arising from asymmetry. Additionally, the AH model avoids the numerical blowup encountered in simulations of odd-order nonlinear systems [20,62], making it a suitable framework for studying higher-order effects. Nevertheless, the mechanisms by which asymmetry amplifies higher-order contributions and leads to a steeper scaling exponent remain open questions.

We further examined the AH- $[eqn]$ model, in which asymmetry and quartic nonlinearity are simultaneously present. The thermalization time in this case is well described by MR, and in the strong nonlinearity limit the inverse-square law for Teq is recovered. However, a clear deviation between numerical data and the uncorrected MR prediction is observed (see Figure 5), particularly in the intermediate regime. This deviation arises from the interdependence of asymmetry and nonlinearity, which violates the independence assumption underpinning MR. Such interwoven effects often give rise to complex behaviors, highlighting the importance of isolating and characterizing individual mechanisms. The AH and AH- $[eqn]$ models together offer a clear and tractable framework for this purpose.