{'fraction_non_alphanumeric': 0.03698413311199475, 'fraction_numerical': 0.017179477280623866, 'mean_word_length': 4.456920160101117, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous Heat Conduction in Three-Dimensional Nonlinear Lattices

Heat conduction in three-dimenisional nonlinear lattice models is studied using nonequilibrium molecular dynamics simulations. We employ the FPU model, in which there exists a nonlinearity in the interaction of biquadratic form. It is confirmed that the thermal conductivity, the ratio of the energy flux to the temperature gradient, diverges in systems up to 128x128x256 lattice sites. This size corresponds to nanoscopic to mesoscopic scales of several tens of nanometers. From these results, we conjecture that the energy transport in insulators with perfect crystalline order exhibits anomalous behavior. The effects of lattice structure, random impurities, and natural length in interactions are also examined. We find that face-centered cubic (fcc) lattices display stronger divergence than simple cubic lattices. When impurity sites of infinitely large mass, which are hence fixed, are randomly distributed, such divergence vanishes.


I. INTRODUCTION
It is a long-standing problem to reproduce irreversible heat conduction phenomena, which are described by the Fourier law, J = −κ∇T , where κ denotes the heat conductivity, on the basis of time-reversible microscopic dynamics. Recent studies have succeeded in reproducing the Fourier law [1,2,3,4], and other linear nonequilibrium transport phenomena [5,6]. However, there remains an unsolved problem. The studies cited above mainly used particle systems, such as hard spheres and Lennard-Jones particles, for example, but historically, nonlinear lattice systems have also been studied [7,8,9,10,11,12]. Shimada et al. showed that the Fourier law is reproduced in a three-dimensional (3D) polymer-like lattice with an Fermi-Pasta-Ulam-β (FPU-β) type interaction [1], but Shiba et al. observed the divergence of thermal conductivity for 3 and 4D FPU-β lattices [13]. The purpose of this article is to study the energy transport behavior in nonlinear lattice systems.
Let us start with a brief review of the study of heat conduction in nonlinear lattices. It is well-known that no temperature gradients are formed if we use integrable systems, such as harmonic chains [14]. To realize linear temperature gradients, we need features that give rise to thermalization processes. One possibility for generating such a feature is represented by nonintegrable nonlinear interactions [15]. Nonintegrability results in macroscopically irreversible processes, even in conserved systems, even though the basic equation is reversible.
In recent years, heat conduction in 1D nonintegrable chains has been widely investigated [7]. In these studies, it was found that the heat conductivity of nonlinear lattices depends on the system size as κ ∼ N α , where N denotes the length of the system, when the total momentum of the system is conserved. Although the value of the exponent α seems to depend on the model, it is approximately 0.4 for several models, for example, FPU-β lattices [8], FPU-α lattices [9], diatomic Toda lattices [10], and 1D binary hard-sphere gases [16,17] exhibit similar power-law behavior. This anomalous behavior originates in the power-law decay of the energy-flux autocorrelation function, C(t) = J(t) · J(0) , which is the integrand of the Green-Kubo formula [18,19,20], where J(t) is the total heat flux and V is the volume of the system. In 1D nonlinear lattices, V = N is the chain length. The autocorrelation functions of systems with α ∼ 0.4 decay approximately as C(t) ∼ t −β with β ∼ 0.6. Phenomenological explanations for this decay exponent have been proposed using mode-coupling theory [21], the renormalization group analysis of 1D fluctuating hydrodynamics [22,23], and kinetic theory [24]. The same types of divergence of κ have been observed in 1D and 2D fluid systems. For simple classical fluids, it has been well established that the autocorrelation function exhibits power-law decay called "longtime tail" [25,26,27,28,29,30], which takes the form C(t) ∼ t −d/2 for d ≥ 2. This was discovered by Alder and Wainwright [31,32] using computer simulations, and the study of long-time tail is still ongoing [33,34]. In 3D systems, the conductivity converges in the limit N → ∞, as demonstrated by the Green-Kubo formula, eq. (1). In 2D systems, it diverges as ln N . Such dimensionality dependence has been confirmed by computer simulations with hard-particle systems [1,2], Lennard-Jones fluids [3], and other systems. The system dimensionality plays a key role in the behavior of transport coefficents in fluid systems.
How does the conductivity of nonlinear lattice systems depend on their dimensionality? One may naively expect a dependence that is similar to that in fluids, assuming that a phonon gas can be described by hydrodynamic equations of motion. But this appears to be a too naive assumption. Recently, it was shown that the conductivity of a 3D extension of FPU-β lattices displays divergence and that this divergence is consistent with the power-law decay of the autocorrelation function [13].
Although the derivation of the long-time tail from the linearized hydrodynamic equations using the modecoupling hypothesis [22,23,27,28,29,30] implies its validity for mesoscopic to macroscopic scales, but not for microscopic scales, particle systems also possess longtime tails on microscopic scales. In simple particle systems, distances on the order of ten times the meanfree path are sufficient to observe the characteristic size dependence of the conductivity [1,2,3]. In nonlinear lattices [13], however, it has been found that anomalous divergence continues up to systems of at least 64 × 64 × 128 lattice sites. This size corresponds to mesoscopic systems. This means that the energy transport in nanoscopic to mesoscopic crystals is more complicated than that in fluid systems. From these findings, it is clear that further study of this problem is interesting both theoretically and with regard to nano technology. Such anomalous features should contribute to explanation of the anomalous features in heat transport observed in experiments and numerical simulations of nano scale systems such as carbon nanotubes [35,36,37].
Our model and simulation method are presented in the next section. The robust divergence of the thermal conductivity in 3D nonlinear lattices is demonstrated in the following two sections: In §3, a simple cubic lattice is treated, and in §4, an fcc lattice is treated. Such divergence is also shown to exist in a system with a nonlinear interaction possessing a natural length scale in §5. In §6, the effect of the randomness of the mass is studied. In §7, the results for a 2D nonlinear lattice are given. The last section contains a summary and conclusion.

II. MODEL AND SIMULATION METHOD
The first model we study is a simple 3D extension of FPU-β lattices [13]. The model Hamiltonian is where the 3D vectors p i and r i denote the momentum and displacement of a particle at lattice point i, respectively, and the mass m is taken to be unity. The summation over i, j is a sum over nearest neighbors. k and g are parameters indicating the strength of interactions between these nearest-neighbor particles, and they are fixed as k = 1.0 and g = 0.1 in this paper. We point out here that the spatial degrees of freedom of these dynamical models consist of the displacement vector r i = q i − q 0 i , where q i represents the real spatial coordinates of the particles, and q 0 i represents the equilibrium positions of the particles, which are at simple cubic lattice points. In other words, this model system possesses a crystal structure. In this model the longitudinal and transverse modes are treated identically, and thus have the same dispersion relation in the harmonic limit.
The system size is denoted by N x × N y × N z . We express the system size in terms of the number of particles (N x , N y , N z ), not the lengths (L x , L y , L z ). Thus, the system size is dimensionless N c = L c /a (c = x, y, z). Here, a is the lattice spacing constant, which does not appear in the Hamiltonian given by eq. (2). The only characteristic length scale that appears in eq. (2) is k/g, which is a typical length of phonon-phonon interactions.
We use periodic boundary conditions in the x-and ydirections. The particles on both ends in the z-direction are attached to rigid walls, which are separated by one lattice space from the walls. These particles interact with the wall through a potential that is identical to that for particle-particle interactions, and their local temperature is controlled by the Nosé-Hoover method [38]. Thus, the equations of motion are modified for these particles aṡ Here, ζ i denotes the Nosé-Hoover thermostat variables, which obeyζ where T i denotes the temperatures of the heat baths, which are set to T L on the left and T R on the right. (Here, i z increases to the right.) The difference between T L and T R drives the energy transport in the system. The quantity Q is the relaxation time of the heat baths, which is set to unity. For particles in the bulk, the equations of motion arė Particle dynamics simulations were used to study this system. The initial displacements were zero for all particles (i.e., r i = 0) and the initial momenta p i were randomly assigned for each particle. Starting from this initial state, initial relaxation steps were discarded. After the system reached a steady state, the temperature distribution, energy flux, and thermal conductivity were computed. a. Temperature Using the Virial theorem, we define the local temperature T (i) = T (i x , i y , i z ) of a particle on lattice point i as the long-time average of the kinetic energy. (Here, i x , i y . and i z are the labels of the lattice sites, To obtain better accuracy, we averaged the temperature over the N x × N y particles in the same cross section, and thus define T (i z ) as b. Heat flux In this paper, we refer to energy flux as the heat flux, because the system itself is a thermodynamic system that relaxes to equilibrium due to its nonintegrability. Because the heat flux is conserved in the bulk region, there are two ways to measure it. One way is to sum up the local contributions to the energy flow due to the interactions. The microscopic energy transfer from site i to site k is given by where V ik is the interaction potential between sites i and k. By performing this summation, we could obtain the total heat flow density as j z = N −1 j,k j i→k . However, in this study, we used an alternative method; we calculated the work done by the heat baths as Here, the summation is over all the particles to which the left (or right) heat baths are attached. We have confirmed that these two methods yield the same value of the heat flow to at least two significant figures. c. Thermal conductivity In the simulations carried out for this study, the temperature distribution was found to be linear in the system and no temperature gap was observed. Therefore, the thermal conductivity κ(N z ) can be estimated by For the numerical integration of the equations of motion, we used a Störmer-Verlet difference scheme, given by Here, n represents the number of time steps. This scheme exhibits second-order convergence with ∆t. For the bulk particles, without temperature control, this scheme reduces to the leap-frog integrator, which is symplectic. Therefore, the integration of the bulk area is numerically stable for all ∆t satisfying ∆t ≤ 0.02. However, this scheme is not symplectic for the thermostatted particles. The time step ∆t is taken to be small enough for the integration of these particles connected to the heat bath to be stable. Throughout this paper, ∆t is set to 0.02.

III. SIMPLE CUBIC LATTICE
The results for the simple cubic FPU-β lattice are given in this section up to N z = 256. (A previous study gave results up to 128 [13].) For each data point, the results are averaged over simulations starting from five independent initial conditions. First, we present the results for the temperatures (T L , T R ) = (20.0, 10.0). For this case, the spatial temperature profile is plotted in Fig. 1 for systems with N z = 48, 96, and 192 with the positions of the particles rescaled by the system size N z . No boundary gap is observed in this temperature profile, and the results T (i z /N z ) for all values of N z nearly coincide. We thus conclude that we can estimate the thermal conductivity using eq. (10). The system size dependence of the estimated values of κ(N z ) are plotted in Fig. 2 on semi-log and log-log scales, up to the system size 128 × 128 × 256. These results for κ(N z ) suggest logarithmic divergence, although power-law divergence is not excluded. When a power law is fitted to the data for 96 ≤ N z ≤ 256, we obtain κ ∼ N 0.221(4) z . Such a divergence was previously observed up to N z = 128 [13], and it is now confirmed up to N z = 256. Note that the results show no sign of convergence.
The behavior of κ(N z ) when the heat bath temperatures are changed is shown in Fig. 3. For the same N z , thermal conductivity becomes lower for higher (T L , T R ). This is expected because thermalization is enhanced when the temperature is high because of the stronger nonlinear interactions. However, we observe that the divergence of κ(N z ) with increasing system size is observed even at higher temperatures such as (T L , T R ) = (40.0, 20.0).
After beginning the simulation with the initial conditions described in §2, we waited for the time period t w ∼ 10 5 . By this time, the system had relaxed to a nonequilibrium steady state, where the amount of energy flow per time was stationary. For the 128 × 128 × 256 system, which is the largest system considered in the present study, and includes about 4 × 10 6 particles, the total number of simulation steps multiplied by the number of particles is of the order of 10 13 . For the 128 × 128 × 256 system, approximately 20h of CPU time was required to numerically calculate one sample, using an SX8 cluster with a peak performance of 128 GFlops.

IV. FCC LATTICE
We also studied the thermal conductivity of an fcc system in order to determine its dependence on the lattice structure. The Hamiltonian of this model is also that given in eq. (2), but in this case, each particle interacts with its twelve nearest neighbors. Lattices of sizes N × N × (N + 1) are used, and each lattice point is again denoted by (i x , i y , i z ). In these simulations, particles only occupy sites where (i x + i y + i z ) is even. Therefore, the total number of particles is N 2 (N + 1)/2. Nosé-Hoover heat baths were attached to the 1st and (N + 1)th layers in the z-direction. In these layers, the particles are connected to the walls with the same interactions as those in  (1) . This value of the exponent is larger than that for the simple cubic lattice 0.221(4). This suggests that the power-law exponent of the divergence is not determined only by the lattice dimensionality. It is interesting that this divergence seems to be stronger than that in 1D models, for which the exponent was found to be approximately 0.4 [7].

V. NATURAL LENGTH
Using the Hamiltonian in eq. (2), we cannot realize a system in which longitudinal and transverse modes are treated differently. However, in real solids, the dispersion relations for these modes generally differ, and they can have very different structures. To include such effects, we need to consider the interactions among the particles more carefully.
Here, we consider the inclusion of the natural length of the interaction as a first step, and we consider the modified Hamiltonian where dq ij = |q i − q j | − l 0 , and l 0 is the natural length of the interaction. A system comprising a simple cubic lat-tice was investigated using this Hamiltonian by Shimada et al., and they observed normal heat conduction in three dimensions [1]. However, this model has the flaw that its transverse modes are softened, and as a result, the crystal structure was not maintained in their simulations.
In order to realize a system that maintains a crystal structure in thermal equilibrium employing such a class of Hamiltonians, we need to employ a stable structure. One possibility is the fcc lattice, which can maintain its structure with only nearest-neighbor interactions. In this section, we present simulation results for systems comprising fcc lattices, employing the Hamiltonian in eq. (11). Here, we set the parameters as  The temperature profiles are plotted in Fig. 5. We again observe that the temperature gaps near the walls are small and that the results T (i z /N z ) are the same for different system sizes N . We can therefore obtain the system size dependence of thermal conductivity using eq. (10). As in the previous case, we consider N ×N ×(N +1) lattices with N 2 (N + 1)/2 particles.
The system size dependence of thermal conductivity is plotted in Fig. 6 on both semi-log and log-log scales. It is clearly seen that the divergence is not logarithmic but of a power-law type. Fitting the data for the region N ≥ 48, we obtain a tentative value for the power-law exponent of α ∼ 0.295 (4).
This result suggests that the anomalous behavior we have observed is not peculiar to a Hamiltonian of the form of eq. (2), and we conjecture that it is a general property of insulating solids. So far, all the results we have obtained for nonlinear lattice models with perfect crystalline order are consistent with anomalous thermal conductivity. It is a future problem to investigate the extent to which this anomalous behavior can be observed.

A. Random mass
The purpose of this paper is to show that the divergence of thermal conductivity is a robust property of nonlinear lattices, even in 3D cases. It has been a longstanding belief that the following play important roles in thermal and transport phenomena: (i) nonlinearity in the interactions, (ii) the dimension of the system, and (iii) impurities or disorders. In the previous sections, we showed that anomalous behavior exists if we include only (i) and (ii). In this section, we investigate some systems with disorder.
In 1D cases, disordered harmonic chains have been investigated for many years [39,40,41,42]. For such systems driven by a Langevin thermostat, it has been proved that there is a unique steady state with κ(N ) ∼ N 1/2 [40,41,42].
Li et al. investigated heat transport in disordered FPU chains using Nosé-Hoover thermostats [43]. They investigated heat transport at various temperatures and found that for sufficiently high temperatures, the thermal conductivity diverges as κ(N ) ∼ N 0.43 . Although in the case of Nosé-Hoover thermostats, at low temperatures a unique nonequilibrium steady state might not exist, when the temperature T is much higher than 0.5, there should be no such problem. To eliminate this possible problem, therefore, we only investigated 3D systems at high temperatures. The model we investigated has the following Hamiltonian: The only difference between this Hamiltonian and that appearing in eq. (2) is that, in the above form, the masses of the particles m i vary among lattice sites. This allows us to include disorder in the form of random masses. We chose each m i as a random number given by where λ is a parameter adjusting the amplitude of randomness, and R i is a random number distributed uniformly on the interval of [0, 1). We set the average mass m 0 to unity. Using the Hamiltonian in eq. (13), we carried out simulations for a system with a simple cubic lattice. The setup for the simulation was similar to that for the simulations discussed in §3, with N x : N y : N z = 1 : 1 : 2. The averages of the thermal conductivities obtained in 5 mass configurations were computed for various system sizes and two values of λ. The system size dependence of the thermal conductivity is plotted in Fig. 7. It is observed that the thermal conductivity is lower when the effect of the disorder is stronger and that, here too, thermal conductivity diverges with increasing system size. We conclude that disorder does not destroy the anomalous thermal conductivity.

B. Randomly fixed sites
Because we could find no indication of normal heat conduction in the simulations considered so far, we also studied one extreme situation when the disorder effect should be very strong, that in which the system possesses impurities of infinite mass. Although the form of the Hamiltonian in this case is also given by eq. (13), the mass was chosen as with a certain fraction of the masses set to infinity. In other words, we fixed a certain fraction of the particles at their mechanical equilibrium position with r i = 0. We selected the fixed particles randomly from all the particles that are not connected to the heat baths with proportions of 10% and 20%. These values are both below the critical percentage above which the paths of the heat flow are completely blocked.
The system size dependence of the thermal conductivity is plotted in Fig. 8. It is seen that although the model with 10% fixed sites exhibits some divergence, the thermal conductivity of the model with 20% fixed sites clearly converges to a finite value as N z is increased. These results suggest that defects in crystalline solids play an important role in thermal conductivity behavior. Moreover, the conductivity clearly depends on the density of the fixed sites.
We now remark on past works investigating the 2D harmonic model with missing bond defects [44] (not missing particles), in which each particle has a scalar (not vector) dynamical variable, (p i , q i ). This model exhibits convergent thermal conductivity when there are sufficiently many missing bond defects, similarly to the system studied here. More accurate, quantitative, and detailed investigations of the role of the various type of defects will be studied in the future.

VII. TWO-DIMENSIONAL LATTICE
Here, we comment on some investigations of 2D systems. We have already presented evidence that a 3D nonlinear lattice system generally exhibits power-law dependence on its system size. Because it is unlikely that the thermal conductivity of a 2D system has weaker divergence than that of a three-dimensional system, we conjecture that 2D nonlinear lattice systems also exhibit power-law divergence of the thermal conductivity.
There have not been many studies on 2D systems. Although it is not yet proven, it is widely believed that the thermal conductivity for such systems diverges logarithmically as the size increases [7,45,46]. We simulated a system consisting of a simple square 2D  FPU-β lattice whose Hamiltonian is of the same form as eq. (2). We fixed the aspect ratio to N x : N y = 1 : 2, where the y-axis is taken to be the direction of heat flow. For the size 384 × 768, we waited for approximately t w ∼ 5 × 10 5 for the system to reach the nonequilibrium steady state.
The temperature profile of the steady state is plotted in Fig. 9. It is seen that there are no temperature gaps near the walls. Thus, we can define the thermal conductivity by eq. (10).
The system size dependence of the thermal conductivity is plotted in Fig. 10. There, we see that the conductivity exhibits power-law divergence. Because the values of N y used here are much larger than the values of N z used in the three-dimensional systems that we studied, we can clearly distinguish the behavior here from logarithmic divergence. The data presented here are sufficient to conclude that the divergence is stronger than conventional logarithmic divergence. Our estimation of the power-law exponent is α ∼ 0.268 (3), obtained by fitting in the region N y ≥ 128.

VIII. SUMMARY AND CONCLUSION
Heat conduction in FPU-β lattices was studied using nonequilibrium molecular dynamics simulations. The divergence of the thermal conductivity was observed in simple cubic lattices up to the size 128 × 128 × 256. Such divergence was also observed in fcc lattice systems up to 192 × 192 × 193, in fcc lattice systems with a natural length included in the nonlinear interaction, and in simple cubic lattice systems with randomly distributed particle masses with a variance of 20% (λ = 0.4). Similar divergence for square lattice systems was also found for sizes up to 384×768. These divergences are characterized by power-law dependences on the system size, and the power exponents do not appear to be unique, ranging from 0 (logarithmic) to 0.22 for simple cubic lattice systems, and being approximately 0.5 for fcc lattice systems. Convergence was observed in simple cubic lattice systems with randomly fixed particles.
The present results imply that the long-time-tail decay exponent −d/2 of the energy flux autocorrelation function is not observed in the present systems with diverging thermal conductivity, although it has been observed in particle systems with an order-of-magnitude smaller number of degrees of freedom [1,2,3]. Of course, the size limitation is always an issue to consider in finite-system analysis. Decay with an exponent of −3/2 may appear in larger 3D nonlinear lattices, and thermal conductivity may converge. For this reason, it would be worthwhile carrying out simulations on larger lattices using more powerful computers. However, it should be remarked here that the 128×128×256 lattice is already mesoscopic to macroscopic in the present technological sense: this lattice corresponds to insulators of size 50 × 50 × 100nm 3 when the lattice constant is assumed to be of the scale of diamonds, for example. Experimental observations of heat transport in crystals of this scale should be carried out.
There remain several open problems in addition to the system size problem. One is to elucidate the direction dependence of the heat transport; in the present study, the heat flow was always in the [001] direction. Another is to make the interaction potential function more realistic.