{'fraction_non_alphanumeric': 0.041584865250758524, 'fraction_numerical': 0.018264025224582068, 'mean_word_length': 4.497056899934598, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solitary waves in nonequilibrium molecular dynamics simulations of heat flow in one-dimensional lattices

We study the use of the Evans Nonequilibrium Molecular Dynamics (NEMD) heat flow algorithm for the computation of the heat conductivity in one-dimensional lattices. For the well-known Fermi-Pasta-Ulam (FPU) model, it is shown that when the heat field strength is greater than a certain critical value (which depends on the system size) solitons can be generated in molecular dynamics simulations starting from random initial conditions. Such solitons are stable and travel with supersonic speeds. For smaller heat fields, no solitons are generated in the molecular dynamics simulations; the heat conductivity obtained via the NEMD algorithm increases monotonically with the size of the system.


I. INTRODUCTION
Heat conduction in one-dimensional (1D) lattices has attracted much research interest.
Surprisingly, it has been found that many 1D lattices do not obey the Fourier's law [1]- [17]: the thermal conductivity is divergent in the thermodynamic limit. For harmonically coupled oscillators, it was rigorously shown that the thermal conductivity λ is proportional to the number of oscillators N [1]. Such a divergence is founded in the existence of extended waves (phonons) freely traveling (and carrying thermal energy) along the lattice without attenuation. In later studies [2,3], impurities or defects in the chain were taken into account, since it was anticipated that phonon waves could be damped by the scattering processes due to such defects, thus possibly removing the N divergence of λ. However, it was demonstrated for isotopically disordered harmonic chains that the heat conductivity still diverged at a somewhat slower rate (λ ≈ N 1/2 ) [2,3]. Another way of trying to achieve normal heat conduction in one-dimensional lattices is by invoking anharmonicity [4]: here nonlinearity makes it possible for phonons to interact among themselves thus impeding their free propagation. However, Lepri et al [12] have found that even strong nonlinearity and chaotic behavior is insufficient to ensure the existence of normal heat conduction. In the well-known Fermi-Pasta-Ulam (FPU) β-model they found a power-law divergence of the thermal conductivity λ ∝ N γ for γ ≈ 0.4. This power-law divergence was qualitatively attributed to the long-time tail of the heat flux autocorrelation function, whose time integral gives the thermal conductivity of the system [13].
Previous studies of heat conductivity have used a straightforward simulation method [6]- [17]: In molecular dynamics (MD) simulations two heat reservoirs with high and low temperatures T H and T L respectively are located on each side of the lattice. The average heat flux and the internal temperature gradient are measured, with the thermal conductivity being the ratio of these two quantities. However, there are a number of disadvantages with this approach. In particular, the system is spatially inhomogeneous and one cannot define an intrinsic temperature T for the system due to the large temperature gradient. Consequently it is impossible to obtain the T dependence of the heat conductivity. In addition, problems associated with the use of the Nosé-Hoover thermostat for boundary particles have been discussed in Ref. [15].
Recently Maeda and Munakata [10] proposed a homogeneous nonequilibrium moleculardynamics (NEMD) method based on the Evans heat flow algorithm [18][19][20]. However,the system size used in Ref. [10] was too small (32 particles) to allow one to study the behavior of the NEMD algorithm in the thermodynamic limit.
In this paper we present a detailed study of the Evans NEMD heat flow algorithm for one-dimensional (1D) lattices. We demonstrate that when the heat field is sufficiently large, well-defined solitary waves (solitons) can be generated in simulations with random initial conditions. These soliton objects travel at a supersonic speeds, and they also appear in the corresponding Hamiltonian systems. When a soliton is generated, the normal process of homogeneous heat conduction is destroyed, as heat is transferred through the chain via the highly localized solitary wave. This results in a sharp increase in the effective thermal conductivity of the system. Due to this instability, progressively smaller fields are required, as the system size increases, to observe the linear regime of the thermal conductivity and thereby carry out the extrapolation of the thermal conductivity to zero field strength.
The paper is organized as follows. In section 2 we describe the NEMD equations of motion for one-dimensional lattices with nearest neighbor interactions, and we point out the possibility of solitary wave solutions in the system. In section 3 we carry out the nonequilibrium heat flow simulations on the well-known Fermi-Pasta-Ulam (FPU) model. We show that for the heat field strengths greater than a certain critical value, which depends on the system size and temperature, solitons can be generated from random initial conditions. For smaller heat fields, no solitons can be generated in molecular dynamics simulations with random initial conditions; in this case, the heat conductivity can be obtained via the NEMD algorithm. Some concluding remarks are presented in Section 4.

II. NEMD EQUATIONS AND SOLITARY WAVE SOLUTIONS
We consider a 1D system of N particles located along the x axis with lattice constant a = 1. Each particle is allowed to move in the y direction, perpendicular to the x-axis, and we denote by q i the displacement of the i-th particle and by p i the corresponding momentum.
The Hamiltonian of the system is expressed as where m i is the mass of the i-th particle, and U represents the interparticle interaction potential. For example, the FPU β−model is exemplified by U(r) = 1 2 r 2 + β 4 r 4 , r ≡ q i+1 − q i . The Evans thermal conductivity algorithm for general N-particle systems of fluid has been described in detail [18,19]. The basis of this algorithm is the Green-Kubo relation for the thermal conductivity. For a 1D system defined by the Hamiltonian model (1), the thermal conductivity coefficient is given by where k B is the Boltzmann's constant, T is the absolute temperature of the system, and L = Na is the system size (The lattice spatial constant a is taken to be 1 in this paper).
The notation · · · eq denotes an equilibrium ensemble average. The heat flux vector, J x , is given by In the Evans NEMD algorithm the N particle system is coupled to a "heat field" F e . The coupling is defined in such a way that the energy dissipation is proportional to J x F e , i.e., dH/dt=LJ x F e , and that the adiabatic incompressibility of phase space (AIΓ) condition is satisfied [19]. The thermal conductivity coefficient can then be found from the ratio of the heat flux to the applied heat field: replaced by a time average of J x (t) if the nonequilibrium steady state is unique.
The equations of motion which satisfy the above conditions arė where is the total force on particle i due to the nearest neighbor interaction, and α is the thermostat multiplier. For a Gaussian thermostat, α is This Gaussian thermostat ensures that the system's kinetic energy is fixed at the constant value K 0 . We note that if F e = 0, the total force exerted on the system by the heat field is exactly zero, i.e., F e D i = 0 from Eq. (7). Thus the total momentum, p i (t), will remain zero for t > 0 provided that its initial value vanishes.
We have simulated the above described NEMD equations (5)-(8) for a variety of interparticle interaction potentials, and we observe well-defined solitary wave excitations. Before describing in detail our numerical results we first analyze why such solitary waves can occur.
For simplicity, we set m i = 1 for all i. From Eqs. (5) and (6) we obtain This equation, together with Eq.(8) for the definition of α, forms a set of lattice equations for the variable Q i . Note that if F e and α are set equal to zero, then Eq.(10) becomes a lattice system which can support stable supersonic solitary waves for many types of nearest neighbor interaction potentials U (see, e.g., Refs. [21][22][23][24][25][26]). But exact analytical solutions for the solitary waves are not available except for some rather special potential functions U.
To analyze Eq.(10) we can use the so-called quasi-continuum approximation techniques of Refs. [21][22][23][24][25][26] to reduce Eq.(10) to a partial differential equation: which can be considered as a generalized Boussinesq equation [24]. Here due to the periodic boundary conditions used for the system (5,6), i.e., q N +1 = q 1 , The solution of Eq.(11) must We assume that Q(x, t) = Q s (x − V t) + Q 0 , where Q 0 is a constant background and Q s is a localized soliton profile which vanishes as z = x − V t → ±∞. Then from Eq.(12) we have Substituting Q(x, t) = Q s (x−V t) + Q 0 into Eq.(11) we find that the solitary wave profile But it does not seem possible to find an analytical solution for this partial differential equation for the general case of Q 0 = 0, F e = 0 and α = 0.
However, we note that if the system is large (N ≫ 1) then Q 0 must be very small according to Eq. (13). Moreover, assuming that F e and α are small parameters, we can set Q 0 = 0, F e = 0, and α = 0 in Eq. (14). After such a drastic simplification, it is possible to find approximate solitary wave solutions analytically for various types of interparticle potentials U(Q) [21][22][23][24][25][26]. Here we consider a particular example, the FPU potential U(Q) = 1 2 Q 2 + 1 4 βQ 4 . According to Refs. [21][22][23][24][25][26], the solitary wave profile Q(z) can be approximated by Once such an approximate soliton solution is obtained then we can go back to the original lattice variables (q i , p i ). Notice that Q i = q i − q i−1 , thus Consequently, the contribution to the kinetic energy from the soliton is In numerical simulation we observe that when a soliton is generated, small amplitude phonon waves give a negligible contribution to the system's kinetic energy. In such a case, we have where T is the absolute temperature of the system which is fixed by the Gaussian thermostat In the next section, we present numerical simulation results to confirm these analytical predictions.

III. SIMULATION RESULTS
Following Ref. [10] we apply the NEMD heat flow algorithm [Eqs. (5,6)] to the FPU βmodel for which U F P U (r) = r 2 /2 + βr 4 /4, and the parameter β is taken to be 1 without loss of generality. Periodic boundary conditions are always used. Unless indicated otherwise, the initial conditions for q i and p i are obtained by a random number generator. The equations of motion are integrated using a fourth-order operator-splitting integrator which conserves the system's kinetic energy [27]. The time stepsize is δt = 0.005 and the total simulation time is between 10 4 to 10 6 units for each trajectory.
The first feature of note is that for a given temperature T and particle number N, stable solitary waves can be generated during simulations (which start from random initial conditions) if the heat field strength F e is greater than a certain critical value F cr . The solitary wave travels in the direction of heat flow (to the right in Fig.1) with a supersonic speed (V s > 1). When the soliton is generated the normal process of homogeneous heat conduction is destroyed and the heat flux increases drastically (Fig. 2). In such a case, heat is transported in the form of a highly localized energy pulse carried by the soliton, and the average value of heat flux is nearly independent of F e .
We find that the soliton's velocity increases with temperature and system size (Fig.3), but it is nearly independent of the heat field strength F e (≤ 1.0). For example, for N = 100, T = 10, the soliton's velocities are found to be about 4.7, 4.7 and 4.8 for F e = 0.01, 0.1, and 1.0, respectively. This is in good qualitative agreement with that predicted by equation (18). However, it is found that the analytical result (18) overestimates the soliton's velocity.
For instance, for N = 200, T = 1, Eq.(18) gives V = 5.586; in the numerical simulation we found that the soliton's velocity is about 3.2. The discrepancy is mainly due to the fact that Eq. (18) is obtained by using a crude approximation of the soliton solution for the lattice system (10). Similar to the case of Hamiltonian lattices [21][22][23][24][25][26], here we also find that the soliton's amplitude increases with its velocity.
Although the spontaneous generation of solitary waves (from random initial conditions) is observed in the NEMD simulations only when the heat field is strong enough, such waves, once generated, continue to exist in the system after the heat field and thermostat are switched off (Fig.4). The reason for such a behavior is that the soliton is an inherent excitation in the FPU β-model, as explained in the previous section.
When the heat field strength F e is smaller than the critical value no soliton can be generated from random initial conditions (See Fig.5). In this case, the time-averaged heat flux J x (t) can be measured (See Fig.6), and the conductivity, lim t→∞ J x (t) /(T F e ), can be calculated. In Fig.7 the heat conductivity obtained through the NEMD simulations is plotted (Error bars are estimated to be within 10% at most). The NEMD heat conductivity increases with the system size. When F e → 0 the conductivity converges to a finite value, which, according to the generalized Green-Kubo realtion [19,28], should equal the conductivity obtained through Eq.(2). We have tested this convergence for the FPU system of N = 100 particles. In Fig.8 we plot the function where the ensemble average is obtained by using ten independent trajectories of the length 10 6 units in time. It is clear that when t → ∞, λ approaches to a value around 93, which is in good agreement with the heat conductivity obtained through NEMD algorithm (See Fig.7).
Finally, we investigate how the critical field strength F cr for generating solitons depends on the particle number N and the system temperature T . We found that F cr increases with T . For example, at N = 100 the critical value F cr is around 0.0054 and 0.0085 for T = 1 and T = 10, respectively. On the other hand, we found that F cr decreases monotonically with N as shown in Fig.9. For a system of 10000 particles, the critical field is as small as 0.0005 (accuracy is within ±10 −4 ). This means that an extremely small field has to be used in order to observe the linear regime (with no solitons) of heat conduction. However, when the heat field is too small the noise-signal ratio [when using Eq.(4)] would become too large and thus the efficiency of the NEMD algorithm will be drastically reduced.

IV. CONCLUDING REMARKS
In conclusion, we have shown that the Evans NEMD heat flow algorithm, which was designed originally for computing thermal conductivity in liquids, can generate solitons when it is applied to 1D lattices. In the well-known FPU model, we have shown that when the heat field strength is greater than a certain critical value a soliton can be generated from random initial conditions. Such a soliton is stable and it travels with a supersonic speed which is determined by the system size and temperature. Because of this instability, progressively smaller fields have to be used (as the system size increases) to observe the linear regime of the thermal conductivity and thereby carry out the extrapolation to zero field.
This greatly reduces the efficiency with which the algorithm can be used to compute the thermal conductivity of large 1D lattices. Nevertheless, for small systems, we have found that the heat conductivity increases with the size of the system, which is in qualitative agreement with the previous finding [12][13][14].
The present study has been primarily focused on FPU β-model, but we have checked numerically that similar phenomena also exist for other types of 1D lattices with distinct interparticle interaction potentials (e.g., the Toda potential and Morse potential) and even for diatomic lattices. In particular, the spontaneous generation of solitons can be observed not only for the nonequilibrium heat flow systems with Gaussian thermostat but also with the Nosé-Hoover thermostat and an isoenergetic thermostat. The mechanism underlying the observed chaos-soliton transition above the critical field strength still remains to be identified, particularly in terms of Lyapunov spectra-shift and phase space contraction [19] induced by the nonequilibrium heat flow algorithms.

ACKNOWLEDGMENT
We thank Debra J. Searles and Professor E G.D. Cohen for useful discussions and comments about this work. This work is supported by ARC large grant A69800064.