{'fraction_non_alphanumeric': 0.04674539575425277, 'fraction_numerical': 0.018206101504287923, 'mean_word_length': 4.16031918752267, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Validity of Fourier's law in one-dimensional momentum-conserving lattices with asymmetric interparticle interactions

We have numerically studied heat conduction in a few one-dimensional momentum-conserving lattices with asymmetric interparticle interactions by the nonequilibrium heat bath method, the equilibrium Green-Kubo method, and the heat current power spectra analysis. Very strong finite-size effects are clearly observed. Such effects make the heat conduction obey a Fourier-like law in a wide range of lattice lengths. However, in yet longer lattice lengths, the heat conductivity regains its power-law divergence. Therefore the power-law divergence of the heat conductivity in the thermodynamic limit is verified, as is expected by many existing theories.


I. INTRODUCTION
Heat conduction induced by a small temperature gradient in the stationary state is expected to satisfy Fourier's law: j = −κ∇T . A complete understanding of the microscopic mechanism that determines this law remains a challenging problem of nonequilibrium statistical mechanics [1]. Fourier's law implies that a fixed temperature difference ∆T that is applied to a homogeneous material with length L induces a stationary state heat current j that should be inversely proportional to L: j = −κ∆T /L. On the other hand, j is found, by numerical simulation, to decay as L −1+α with a positive α in many one-dimensional (1D) lattice models [2]. This implies an infinite κ which diverges with L as L α in the thermodynamic limit. When an on-site potential is absent, i.e., the particles in the lattice only interact with other particles, the total momentum is conserved. It has been generally accepted for many years, with only very few exceptions such as the coupled rotator lattice [3], that heat conduction in such a lattice belongs to this class [1,4]. Many theories, e.g., renormalization group [5] and mode coupling [6,7], support such a conclusion. As for the detailed value of the divergency exponent α, however, different theories suggest differently. Mainstream expectations include 1/3 [5,6], 2/5 [7], and 1/2 [8]. More details can also be found in the review articles Ref. [1].
It is worth mentioning that such a topic is not merely a purely academic issue. The rapid progress nowadays in nanotechnology has already enabled us to experimentally measure the size dependence of the heat conductivity in many one-dimensional [9] and two-dimensional [10] microscopic materials. The topic thus also has extensive application values [11].

II. MODEL AND SIMULATION
Among various models, the lattice model, due to its simplicity, has been studied the most. The Hamiltonian of a 1D lattice takes the form where x i is the displacement of the ith particle relative to its equilibrium position. Without an on-site potential, the lattice is total momentum conserving. In this paper we focus on the lattices with asymmetric interparticle interactions, namely V (x) = V (−x). The first one we study is the Fermi-Pasta-Ulam (FPU)αβ lattice, Throughout the paper we set k 2 = 1 ,k 3 = 2, and k 4 = 1. This is the simplest model with asymmetric interparticle interactions. We study its heat conduction by the nonequilibrium heat bath method first. Fixed boundary conditions are applied to a lattice with length L. The left-and rightmost ends of the lattice are coupled to Langevin heat baths with temperatures T L = 1.5 and T R = 0.5, respectively. In the nonequilibrium stationary state, the heat conductivity is defined as: where J is the heat current and ∇T is the temperature gradient in the lattice. A number of independent runs starting from different randomly chosen initial states are performed. The simulation time depends on the lattice size. For the longest lattice (L = 65536) the average is performed over 4 × 10 8 dimensionless time units. A very flat heat conductivity is observed in the lattice length L ranging from several hundred to ten thousand. However, in yet longer lattices, say, L > 1 × 10 4 , the running slope of the heat conductivity κ grows again. In fact such a curving up of κ was reported in our previous work [12]. But in this study the asymmetric term k 3 = 2 is much greater than that in the previous work. By comparing the two results, we see that even the asymmetry is increased to so high a value (k 3 = 2 is the maximum value that keeps the potential single well, at the given k 2 = k 4 = 1.), the curving up of κ can only be postponed, but its asymptotic behavior is not affected. We next use the equilibrium Green-Kubo method [13], which provides an alternative way to determining α with much higher cost efficiency [12,14,15]. The rescaled heat current autocorrelation function is defined as where J(t) ≡ i j i (t) is the instantaneous global heat current and N is the total number of particles. The Boltzmann constant k B is set to unity. The simulations are carried out in lattices with periodic boundary conditions since they provide the best convergence to the thermodynamic limit. Microcanonical simulations are applied with zero total momentum and identical energy density ǫ = 0.846 that corresponds to the temperature T = 1. Fig. 2 depicts c(τ ) as a function of the time lag τ for various particle numbers N . We see that in the shortτ regime, c(τ ) follows a very fast exponential-like decay, just like what was reported in Ref. [16]. However, when τ is greater than a threshold value τ c = 700, c(τ ) recovers a slow power-law decay. A reasonable estimate of the asymptotic exponent for this certain model, which is also supported by a theoretical expectation [5], is −2/3. In Ref. [16] because c(τ ) is studied in a short range of values of τ (τ < 400), the real asymptotic decay is not observed. Here we emphasize that such a crossover at τ c is not induced by an insufficient particle number of the lattice in the calculation of the heat current correlation c(τ ). First, the curves with different N well overlap with each other around the crossover regime. Second, and even more important, it is clearly seen by comparing those curves that the smaller the N the faster is the c(τ ) decay. Thus the slowing down of the decay is not a consequence of an insufficient particle number N .
It is commonly accepted that the length dependence of the heat conductivity κ NE (L) follows the κ GK (L) defined below, where t = L/v s with at most a constant-factor difference as L → ∞ [1]. If c(τ ) decays asymptotically with τ as τ γ (γ > −1), then the heat conductivity should diverge with L as L α (α = 1 + γ). Otherwise, if c(τ ) decays faster than τ −1 , then the heat conductivity converges.
The speed of sound, v s , can be obtained by simulating the heat diffusion process [17]. In the FPU-αβ lattices v s ≈ 1.5. At first glance, since c(τ ) displays a slow power-law decay τ −2/3 from τ = τ c = 700 upwards, it seems reasonable to predict that κ follows L 1/3 from L = L c ≡ τ c v s upwards. However, the finding in Fig. 1 is evidently inconsistent with this expectation. To explain this inconsistency, we suppose, as is observed in Fig. 2, when τ is greater than a threshold value τ c , c(τ ) follows Bτ γ .
Then for L > v s τ c , the length dependence of the heat conductivity is where A ≡ τc 0 c(τ )dτ and D ≡ A(1+γ) B − τ c 1+γ are constants independent of L. The running slope of κ GK (L) in double logarithmic scale equals When L → v s τ c +, α(L) → B A τ 1+γ c . In the thermodynamic limit, α → 1 + γ.
Therefore, in order to observe a value of α ∈ ( B A τ 1+γ c , 1 + γ) by the nonequilibrium heat bath method, we need to simulate a lattice with a length of at least For the FPU-αβ lattice, v s = 1.5. According to Fig. 2, τ c = 7 × 10 2 , A = 59, γ = −2/3, and B = 0.7. Thus D = 19.2 and α(L) ∈ (0.1, 1/3). To observe α close to its asymptotic value 1/3, say, α = 0.3, a simulation with length L(0.3) = 7.7 × 10 6 is necessary. In contrast, the longest lattice that has already been simulated is only two orders of magnitude shorter. This clearly explains why a flat κ has been observed in Fig. 1.
The third method we use to give support to such a crossover is the power spectrum S(ω) of the global heat current. In momentum-conserving lattices, heat conduction in the thermodynamic limit is dominated by low-frequency phonons. According to the Wiener-Khinchin theorem, in the low-frequency limit, the heat current autocorrelation function c(τ ) ∼ τ γ corresponds to S(ω) ∼ ω −γ−1 if γ > −1, and a flat S(ω) otherwise [7]. The power spectrum S(ω) of the FPU-αβ lattices with different particle numbers are depicted in Fig. 3. The conditions of the calculations are the same as those we use to calculate c(τ ). A flat S(ω) regime, ω ∈ (1 × 10 −4 ∼ 5 × 10 −3 ) is observed [18]. This regime clearly corresponds to the flat κ regime shown in Fig. 1, as well as the fast exponential-like decay regime of c(τ ) shown in Fig. 2. However, in yet a lower ω regime, S(ω) regains its slope, approaching ω −1/3 . This corresponds to the curving up of κ shown in Fig. 1, as well as the slow power-law decay of c(τ ) shown in Fig. 2. By comparing curves for different N , we conclude that such a crossover is definitely not induced by a finite N effect.
The conclusions by the three different methods are well consistent with each other. All of them confirm the theoretical expectation that heat conduction generally displays a power-law divergence in 1D nonlinear lattices with total momentum conservation. Asymmetric interparticle interaction can only induce a transient finite-size flat κ whereas its asymptotic divergence remains unchanged. The temperature is quite important in this model because when it is higher the interparticle interaction is dominated more by the symmetric quartic term. To check the role of the temperature, we have also calculated c(τ ) in the same lattice but with a much lower energy density ǫ = 0.0977, which corresponds to T = 0.1. Similarly, a fast exponential-like decay is observed in a range of τ until τ reaches a critical value τ c ≈ 2 × 10 3 , which is much longer (see the inset in Fig.2). However, the asymptotic decay remains unchanged again.
Finally we study lattices with asymmetric interparticle interactions, called the LWAII model in Ref. [19], Such an asymmetry induces an even stronger finite-size effect. We have also calculated its rescaled heat current autocorrelation function c(τ ) at temperature T = 2.5, which corresponds to ǫ = 2.403 (see Fig. 4). According to Ref. [19], heat conduction at this temperature obeys the Fourier law. Once more we observe an exponential-like decay in an even wider range of τ . The heat current correlation c(τ ) also ceases this fast decay and changes to a slow power-law decay as τ > τ c ≈ 7×10 3 . Calculation for each lattice was performed on one or more Nvidia Tesla-2075 graphics processing units (GPUs), each of which has 448 Cuda processing cores on board. This part of the calculation alone cost several months of wall time in our eight-GPU graphics workstation. Due to rather big fluctuations, we are still not able to obtain the asymptotic power exponent γ with satisfactory accuracy. However, γ < −1 is highly unlikely. That c(τ ) decays asymptotically more slowly than τ −1 is highly unlikely.
The theoretical expected value γ = −2/3 is again a reasonable estimate. For the LWAII model v s = 2.0. According to Fig. 4, τ c = 7 × 10 3 , A = 640, and B = 0.5. Thus D = 408 and α(L) ∈ (0.015, 1/3). In Fig. 5, we plot the running slope α as a function of the required length L. The longest lattice that is studied in Ref. [19] is L = 5 × 10 4 (indicated by the black vertical dashed line in the figure) where α(L) = 0.022. The value is very close to zero. The Fourier-like law that has been observed is thus understandable. To observe α approaching 0.3, simulating a lattice with length L(0.3) = 1 × 10 11 , which is more than six orders of magnitude longer, is required. We have to say that it is unlikely to achieve such a goal in the near future. Such a length is in the macroscopic scale already (10 m if the lattice constant is 1Å). The result for the FPU-αβ model is also plotted for reference.

III. DISCUSSION
In summary, we studied heat conduction in a few 1D lattices with asymmetric interparticle interactions. Strong finite-size effects are observed in all cases, which make the heat conductivity κ obey a Fourier-like law in a wide regime of the lattice length L. The heat current autocorrelation function c(τ ) accordingly displays a fast exponential-like decay in a wide range of values of τ . These agree with the observations made in existing studies [16,19] very well. This phenomenon can be attributed to the thermal pressure, which only exists in lattices with asymmetric interparticle interactions. The thermal pres-sure makes the particle densities vary with instantaneous local temperatures. This effect induces a kind of disorder, which provides a mechanism to scatter phonons that process heat transport. However, in momentum-conserving lattices the asymptotic behavior of heat conduction is The running slope α of the heat conductivity as a function of the required length of lattice, L, in nonequilibrium heat bath method for the FPU-αβ lattice and LWAII. The longest lattice that has been simulated for LWAII (L = 5 × 10 4 ) is indicated by the black vertical dashed line. To observe α approaching 0.3, simulating a lattice which is greater than six orders of magnitude longer is necessary. dominated by long-wavelength (low-frequency) phonons. The particle density fluctuation is negligible in large spatial scales. As a consequence, in a yet longer-τ regime, c(τ ) recovers a slow power-law decay with an exponent less negative than −1. Accordingly the heat conductivity κ should asymptotically recover the power-law divergence as well. This conclusion is also consistent with some very recent studies [20]. However, it is very hard to directly observe this asymptotic behavior by the nonequilibrium heat bath method because the required lengths of lattice are very long. By analyzing the data of c(τ ), we are able to predict quantitatively those necessary lengths. Their values are so long that the computational demand is beyond the present and near-future capabilities.