First-Principle Validation of Fourier’s Law: One-Dimensional Classical Inertial Heisenberg Model

The thermal conductance of a one-dimensional classical inertial Heisenberg model of linear size L is computed, considering the first and last particles in thermal contact with heat baths at higher and lower temperatures, Th and Tl ( Th>Tl ), respectively. These particles at the extremities of the chain are subjected to standard Langevin dynamics, whereas all remaining rotators ( i=2,⋯,L−1 ) interact by means of nearest-neighbor ferromagnetic couplings and evolve in time following their own equations of motion, being investigated numerically through molecular-dynamics numerical simulations. Fourier’s law for the heat flux is verified numerically, with the thermal conductivity becoming independent of the lattice size in the limit L→∞ , scaling with the temperature, as κ(T)∼T−2.25 , where T=(Th+Tl)/2 . Moreover, the thermal conductance, σ(L,T)≡κ(T)/L , is well-fitted by a function, which is typical of nonextensive statistical mechanics, according to σ(L,T)=Aexpq(−Bxη) , where A and B are constants, x=L0.475T , q=2.28±0.04 , and η=2.88±0.04 .


Introduction
Two centuries ago, Fourier proposed the law for heat conduction in a given macroscopic system, where the heat flux varies linearly with the gradient of temperature, J ∝ −∇T [1].For a simple one-dimensional system (e.g., a metallic bar along the x axis, J = Jx), the heat flux J (rate of heat per unit area) is given by where κ is known as thermal conductivity.In principle, κ may depend on the temperature, although most measurements are carried at room temperature, leading to values of κ for many materials (see, e.g., Ref. [2]).Usually, metals (like silver, copper, and gold) present large values of κ, being considered as good heat conductors, whereas poor heat conductors (like air and glass fibber) are characterized by small thermal conductivities; typically, the ratio between thermal conductivities of these two limiting cases may differ by a 10 4 factor.In most cases, good thermal conductors are also good electrical conductors, and obey the Wiedemann-Franz law, which states that the ratio of their thermal and electrical conductivities follow a simple formula, being directly proportional to the temperature [3].
In the latest years, many works were pursued for validating Fourier's law in a wide variety of physical systems, both experimentally and theoretically.Particularly, investigations for which microscopic ingredients may be responsible for the property of heat conduction were carried, and it has been verified that thermal conductivity may be generated by different types of particles (or quasi-particles).In the case of good electrical conductors the most significant contribution to the thermal conductivity comes from free electrons, whereas in electrical insulators such contributions may arise from quasi-particles, like phonons and magnons, or even from defects.As examples, for antiferromagnetic electrical insulators such as Sr 2 CuO 3 and SrCuO 2 , magnons yield the most relevant contribution for the thermal conductivity, which can be fitted by a 1/T 2 law, at high temperatures [4].Several experimental investigations have verified Fourier's law in a large diversity of systems [4,5,6,7,8], including coal and rocks from coalfields [6], as well as two-dimensional materials [7,8].On the other hand, some authors claim to have found anomalies [9], or even violations of this law, for silicon nanowires [10], carbon nanotubes [11], and low-dimensional nanoscale systems [12].Furthermore, a curious crossover, induced by disorder, was observed in quantum wires, where by gradually increasing disorder one goes from a low-disorder regime, where the law is apparently not valid, to another regime characterized by a uniform temperature gradient inside the wire, in agreement with Fourier's law [13,14].
Since the introduction of the entropy S q in Eq. ( 4), a large amount of works appeared in the literature defining generalized functions and distributions (see, e.g., Ref. [38]).Particularly, a recent study based on superstatistics, has found a stretched q-exponential probability distribution [40], as well as its associated entropic form.
As already mentioned, the latest advances in experimental techniques made it possible to investigate thermal and transport properties and consequently, Fourier's law, in low-dimensional (or even finite-size) systems, like two-dimensional materials [7,8], silicon nanowires [10], carbon nanotubes [11], and low-dimensional nanoscale systems [12].These measurements motivate computational studies in finite-size systems of particles that present their own equations of motion, e.g., systems of interacting classical rotators, whose dynamics may be followed through a direct integration of their equations of motion.In this way, one may validate (or not) Fourier's law, by computing the temperature and size dependence of the thermal conductance.A recent analysis of a system of coupled nearest-neighbor-interacting classical XY rotators [36], on d-dimensional lattices (d = 1, 2, 3) of linear size L, has shown that, for a wider range of temperatures, the temperature dependence of the thermal conductance was better fitted by a more general Ansatz than the q-Gaussian distribution of Eq. ( 2).In fact, Fourier's law was validated in Ref. [36] by fitting the thermal conductance in terms of the functional form of Eq. ( 8), with values of η(d) > 2.
In the present work we analyze the thermal conductance of a one-dimensional classical inertial Heisenberg model of linear size L, considering the first and last particles in thermal contact with heat baths at temperatures T h and T l (T h > T l ), respectively.All remaining rotators (i = 2, • • • , L − 1) interact by means of nearest-neighbor ferromagnetic couplings and evolve in time through molecular-dynamics numerical simulations.Our numerical data validate Fourier's law, and similarly to those of Ref. [36], the thermal conductance is also well-fitted by the functional form of Eq. ( 8).The present results suggest that this form should apply in general for the thermal conductance of nearest-neighbor-interacting systems of classical rotators.In the next section we define the model and the numerical procedure; in Section III we present and discuss our results; in Section IV we pose our conclusions.

Model and Numerical Procedure
The one-dimensional classical inertial Heisenberg model, for a system of L interacting rotators, is defined by the Hamiltonian, where ℓ i ≡ (ℓ ix , ℓ iy , ℓ iz ) and S i ≡ (S ix , S iy , S iz ) represent, respectively, continuously varying angular momenta and spin variables at each site of the linear chain, whereas ⟨ij⟩ denote summations over pairs of nearest-neighbor spins; herein we set, without loss of generality, k B , moments of inertia, and ferromagnetic couplings, all equal to unit.Moreover, spins present unit norm, S 2 i = 1, and at each site angular momentum ℓ i must be perpendicular to S i , yielding ℓ i • S i = 0; these two constraints are imposed at the initial state and should be preserved throughout the whole time evolution.
One should notice that, in contrast with a system of coupled classical XY rotators, where canonical conjugate polar coordinates are commonly used [36], in the Heisenberg case one often chooses Cartesian coordinates [41,42,43].The reason for this is essentially technical, since in terms of spherical coordinates (more precisely, θ, ϕ and their canonical conjugates ℓ θ , ℓ ϕ ), a troublesome term (1/ sin 2 θ) appears in the corresponding equations of motion, leading to numerical difficulties [44,45].However, some of the analytical results to be derived next recover those of the classical inertial XY model for S i = (sin θ i , cos θ i , 0) and ℓ i = ℓ i ẑ.
It is important to mention that previous researches on the thermal conductivity have been carried either for a classical one-dimensional Heisenberg spin model, by using Monte Carlo and Langevin numerical simulations [46], as well as for a classical one-dimensional spin-phonon system, through linear-response theory and the Green-Kubo formula [47].These investigations did not take into account the kinetic contribution in Eq. ( 9), so that in order to obtain the thermal conductivity they assumed the validity of Fourier's law.The main advantage of the introduction of the kinetic term in Eq. ( 9) concerns the possibility of deriving equations of motion, making it feasible to follow the time evolution of the system, through molecular-dynamics simulations, by a numerical integration of such equations.This technique allows one to validate Fourier's law, as well as to obtain its thermal conductivity directly.
In order to carry on this procedure we consider an open chain of rotators with the first and last particles in thermal contact with heat baths at higher and lower temperatures, T h and T l (T h > T l ), respectively (cf.Fig. 1), whereas all remaining rotators (i = 2, • • • , L − 1) follow their usual equations of motion (see, e.g., Refs.[41,42,43]).In this way, one has for sites i = 2, . . ., L − 1, whereas the rotators at extremities follow standard Langevin dynamics, Above, γ h and γ l represent friction coefficients, whereas η h and η l denote independent three-dimensional vectors, η h ≡ (η hx , η hy , η hz ), η l ≡ (η lx , η ly , η lz ), whose each Cartesian component stand for a Gaussian white noise with zero mean and correlated in time, with the indexes µ and ν denoting Cartesian components; from now on, we will set the friction coefficients γ h and γ l equal to unit.
The condition of a constant norm for the spin variables yields which should be used together with ℓ i • S i = 0 in order to eliminate li and calculate Si from Eqs. ( 10) and (11).One has for rotators at sites whereas for those at extremities, Rh Rl For the system illustrated in Fig. 1 we will consider the temperatures of the heat baths differing by 2ε, with ε representing a positive dimensionless parameter; moreover, the temperature parameter T = (T h + T l )/2 will be varied in a certain range of positive values.The set of Eqs. ( 14) and ( 15) are transformed into first-order differential equations (e.g., by defining a new variable V i ≡ Ṡi ) to be solved numerically through Verlet's method [48,49], with a time step dt = 0.005, for different lattice sizes L. The rotators at the bulk (i= 2, • • • , L − 1) follow a continuity equation, where so the stationary state is attained for (dE i /dt) = 0, i.e., J i = J i−1 .
Data are obtained at stationary states, which, as usual, take longer times to be reached for increasing lattice sizes.For numerical reasons, to decrease fluctuations in the bulk due to the noise, we compute an average heat flux by discarding a certain number of particles p near the extremities (typically p ≃ 0.15L).In this way, we define an average heat flux as with whereas ⟨..⟩ denotes time and sample averages, to be described next.Let us emphasize that for S i = (sin θ i , cos θ i , 0) and ℓ i = ℓ i ẑ, one recovers the expression for the heat flux of the classical inertial XY model, i.e., J i = 1 2 (ℓ i + ℓ i+1 )sin (θ i − θ i+1 ) [36,50], showing the appropriateness of the Cartesian-coordinate approach used herein for the classical inertial Heisenberg model.
Let us now describe the time evolution procedure; for a time step dt = 0.005, each unit of time corresponds to 200 integrations of the equations of motion.We have considered a transient of 5 × 10 7 time units to start computing the averages ⟨J i ⟩ in Eq. ( 18), checking that this transient time was sufficient to fulfill the condition J i = J i−1 (within a three-decimal digits accuracy at least), for all values of L analyzed.After that, simulations were carried for an additional interval of 2 × 10 8 time units (leading to a total time of 2.5 × 10 8 for each simulation).The interval 2 × 10 8 was divided into 80 equally-spaced windows of 2.5 × 10 6 time units, so that time averages were taken inside each window; then an additional sample average was taken over these 80 time windows, leading to the averages ⟨J i ⟩.
Using the results of Eq. ( 18) one may calculate the thermal conductivity of Eq. ( 1), and consequently, the thermal conductance, In the next section we present results for both quantities, obtained from the numerical procedure described above.

Results
We simulated the system of Fig. 1 for different lattice sizes, namely, L = 50, 70, 100, 140, considering the heatbath temperatures differing by 2ε, with ε = 0.125.The temperature parameter T = (T h + T l )/2 was varied in the interval 0 < T ≤ 3.5, such as to capture both low-and high-temperature regimes.The values of L (L ≥ 50) were chosen adequately to guarantee that the thermal conductivity κ did not present any dependence on the size L in the high-temperature regime, as expected.
In Fig. 2 we present numerical data for the thermal conductivity [panel (a)] and thermal conductance [panel (b)] versus temperature (log-log representations) and different sizes L. In Fig. 2(a) we exhibit κ(L, T ) (the dependence of the thermal conductivity on the size L, used herein, will become clear below) showing a crossover between two distinct regimes (for T ≃ 0.3), as described next.(i) A low-temperature regime, where κ depends on the size L, decreasing smoothly for increasing temperatures (L fixed).The plots of Fig. 2(a) show that, in the limit T → 0, an extrapolated value, κ(L, 0) ≡ lim T →0 κ(L, T ), increases with L. This anomaly is attributed to the classical approach used herein, indicating that for low temperatures a quantum-mechanical procedure should be applied.(ii) A high-temperature regime, where κ essentially does not depend on L (in the limit L → ∞), as expected from Fourier's law.Moreover, in this regime one notices that κ decreases with the temperature as generally occurs with liquids and solids.For increasing temperature, the thermal conductivity of most liquids usually decreases as the liquid expands and the molecules move apart; in the case of solids, due to lattice distortions, higher temperatures make it more difficult for electrons to flow, leading to a reduction in their thermal conductivity.The results of Fig. 2(a) indicate that the thermal conductivity becomes independent of the lattice size in the limit L → ∞, scaling with the temperature as κ(T ) ∼ T −2.25 at high temperatures.In spite of the simplicity of the one-dimensional classical inertial Heisenberg model of Fig. 1, the present results are very close to experimental verifications in some antiferromagnetic electrical insulators such as the Heisenberg chain cuprates Sr 2 CuO 3 and SrCuO 2 , for which the thermal conductivity is well-fitted by a 1/T 2 law at high temperatures [4].
The same data of Fig. 2(a) is exhibited in Fig. 2(b) where we plot the thermal conductance σ(L, T ) = κ(L, T )/L versus temperature, characterized by the two distinct temperature regimes described above.The low-temperature regime shows that the zero-temperature extrapolated value κ(L, 0) scales as κ(L, 0) ∼ L, leading to σ(L, 0) ≡ lim T →0 κ(L, T )/L ≃ 0.5.On the other hand, in the high-temperature regime the thermal conductance presents a dependence on L, as expected.
In Fig. 3 we exhibit the thermal-conductance data of Fig. 2(b) in conveniently chosen variables, yielding a data collapse for all values of L considered.The full line represents essentially the form of Eq. ( 8), so that one writes where x = L 0.475 T , q = 2.28 ± 0.04, η = 2.88 ± 0.04, A = 0.492 ± 0.002, and B = 0.33 ± 0.04.Notice that this value of η lies outside the range of what is commonly known as "stretched" [cf.Eq. ( 8)], so that the form above should be considered rather as a "shrinked" q-exponential.It should be mentioned that in the case of coupled nearestneighbor-interacting classical XY rotators on d-dimensional lattices (d = 1, 2, 3) [36], the thermal conductance was also fitted by the form of Eq. ( 21), with values of η(d) > 2. Particularly, in the one-dimensional case, such a fitting was attained for x = L 0.3 T , q = 1.7, and η = 2.335, showing that these numbers present a dependence on the number of spin components (n = 2, for XY spins and n = 3, for Heisenberg spins), as well as on the lattice dimension d.By defining the abscissa variable of Fig. 3 in the general form x = L γ(n,d) T , and using the q-exponential definition of Eq. ( 3), one obtains that the slope of high-temperature part of the thermal-conductance data scales with L as where we have introduced the dependence (n, d) on all indices.Since the thermal conductivity (κ = Lσ) should not depend on the size L (in the limit L → ∞), Fourier's law becomes valid for 0.4921e 2.28 The data of Fig. 3 lead to [η(3, 1)γ(3, 1)]/[q(3, 1) − 1] = 1.069 ± 0.083, whereas those for XY rotators on ddimensional lattices yield 1.0007, 0.95, and 0.93, for d = 1, 2, and 3, respectively [36], indicating the validation of Fourier's law for systems of coupled nearest-neighbor-interacting classical n-vector rotators, through the thermal conductance form of Eq. ( 21).

Conclusions
We have studied the heat flow along a one-dimensional classical inertial Heisenberg model of linear size L, by considering the first and last particles in thermal contact with heat baths at different temperatures, T h and T l (T h > T l ), respectively.These particles at extremities of the chain were subjected to standard Langevin dynamics, whereas all remaining rotators (i = 2, • • • , L − 1) interacted by means of nearest-neighbor ferromagnetic couplings and evolved in time following their own classical equations of motion, being investigated numerically through molecular-dynamics numerical simulations.
Fourier's law for the heat flux was verified numerically and both thermal conductivity κ(T ) and thermal conductance σ(L, T ) = κ(T )/L were computed, by defining T = (T h + T l )/2.We have found that the slope of high-temperature part of the thermal-conductance data scales with the system size as σ ∼ L −1.069 , indicating that in the limit L → ∞, one should get a thermal conductivity independent of L. Indeed, in this limit, we have found κ(T ) ∼ T −2. 25 , for high temperatures.The whole thermal-conductance data was well-fitted by the function σ(L, T ) = A exp q (−Bx η ), typical of nonextensive statistical mechanics, where A and B are constants, x = L 0.475 T , q = 2.28 ± 0.04, and η = 2.88 ± 0.04.Since the value of η found herein lies outside the range of what is commonly known as "stretched" (0 < η ≤ 1), herein we called this fitting function of a "shrinked" q-exponential.The present results reinforce those obtained recently for XY rotators on d-dimensional lattices [36], indicating that Fourier's law should be generally valid for systems of coupled nearest-neighbor-interacting classical n-vector rotators, through the "shrinked" q-exponential function for the thermal conductance, with the indices q(n, d) and η(n, d) presenting a dependence on both number of spin components and lattice dimension.
In spite of the simplicity of the model considered herein, the results for the thermal thermal conductivity at high temperatures (κ(T ) ∼ T −2.25 ) are very close to experimental verifications in some antiferromagnetic electrical insulators such as the Heisenberg chain cuprates Sr 2 CuO 3 and SrCuO 2 , for which the thermal conductivity is well-fitted by a 1/T 2 law at high temperatures [4].Since nonextensive statistical mechanics has been used in the description of a wide variety of complex systems, one expects that the present results should be applicable to many of these systems in diverse non-equilibrium regimes.

Figure 1 :
Figure 1: Illustration of the system defined in Eq. (9), where the rotators at extremities of the chain are subjected to heat baths at different temperatures.The hot (R h ) and cold (R l ) reservoirs are at temperatures T h = T (1 + ε) and T l = T (1 − ε), respectively, leading to an average heat flux J = Jx throughout the bulk (see text).The rotators at sites i = 2, . . ., L − 1 interact with their respective nearest neighbors.

Figure 2 : 3 .
Figure 2: (Color online) Numerical data for the thermal conductivity [panel (a)] and thermal conductance [panel (b)] are represented versus temperature (log-log plots) for different sizes (L = 50, 70, 100, 140) of the one-dimensional classical inertial Heisenberg model.One notices a crossover between the low-and high-temperature regimes for T ≃ 0.3.As expected, higher temperatures amplify the effects of the Gaussian white noise, leading to larger fluctuations on numerical data, as shown clearly on panel (a).All quantities shown are dimensionless.

Figure 3 :
Figure 3: The plots for the thermal conductance of Fig. 2(b) are shown in a log-log representation, for a conveniently chosen abscissa (x = L 0.475 T ), leading to a collapse of data for all values of L considered.The fitting (full line) is given by the function of Eq. (21).