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

Abstract

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, $T_h$ and $T_l$ ($T_h>T_l$), respectively. These particles at the extremities of the chain are subjected to standard Langevin dynamics, whereas all remaining rotators ($i=2,\cdots ,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\to \infty$, scaling with the temperature, as $\kappa \left(T\right)\sim T^{-2.25}$, where $T=(T_h+T_l)/2$. Moreover, the thermal conductance, $\sigma (L,T)\equiv \kappa \left(T\right)/L$, is well-fitted by a function, which is typical of nonextensive statistical mechanics, according to $\sigma (L,T)=A{exp}_q(-B{x}^{\eta })$, where A and B are constants, $x=L^{0.475}T$, $q=2.28\pm 0.04$, and $\eta =2.88\pm 0.04$.

1. Introduction

Two centuries ago, Fourier proposed the law for heat conduction in a given macroscopic system, where the heat flux varies linearly with the temperature gradient, $\mathbf{J}\propto -\nabla T$ [1]. For a simple one-dimensional system (e.g., a metallic bar along the $\widehat{x}$ axis, $\mathbf{J}=J\widehat{x}$), the heat flux J (rate of heat per unit area) is given by

$$J=-\kappa \frac{dT}{dx},$$

(1)

where $\kappa$ is known as thermal conductivity. In principle, $\kappa$ may depend on the temperature, although most measurements are carried at room temperature, leading to values of $\kappa$ for many materials (see, e.g., Ref. [2]). Usually, metals (like silver, copper, and gold) present large values of $\kappa$, and are considered good heat conductors, whereas poor heat conductors (such as air and glass fiber) are characterized by small thermal conductivities; typically, the ratio between the 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 follows a simple formula, being directly proportional to the temperature [3].

In recent years, numerous studies have been conducted to validate 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 out, 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 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. For instance, for antiferromagnetic electrical insulators, such as Sr₂CuO₃ and SrCuO₂, which surprisingly behave as $S=1/2$.

Heisenberg chains, magnons yield the most relevant contribution for the thermal conductivity, which can be fitted by a $1/T^2$ law, at high temperatures [4]. In these materials, the low-temperature regime presents ballistic-like heat conduction, increasing as the size of the system increases, while the high-temperature regime presents normal heat conduction [5,6,7].

Being a classical result, there is, in principle, no reason why Fourier’s law should generally apply to physical systems. This aspect has generated controversies in the literature, both in experimental and theoretical studies (for a comprehensive theoretical discussion, see, e.g., Ref. [8]). As a typical anomaly, the thermal conductivity $\kappa$ (which should be an intensive quantity) appears, in many cases, to depend on the size of the system, e.g., on the total number of constituents, as it happens for chains of nonlinear oscillators, where $\kappa$ increases with the total number of elements [9]. This anomaly is usually considered a failure of Fourier’s law. In addition, non-Fourier heat conduction can also emerge from the Maxwell–Cattaneo–Vernotte hyperbolic heat equation, which represents the relativistic version of the heat equation [10]. Recent advances in non-Fourier heat conduction can be found in the work by Benenti et al. [11].

Several experimental investigations have verified Fourier’s law in a diverse range of systems [4,12,13,14,15], including coal and rocks from coalfields [13], as well as two-dimensional materials [14,15]. On the other hand, some authors claim to have found anomalies [16], or even violations of this law for silicon nanowires [17], carbon nanotubes [18], and low-dimensional nanoscale systems [19]. 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 [20,21].

From the theoretical point of view, many authors have investigated Fourier’s law in a wide diversity of models [9,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46], like a Lorentz gas [23], biological [30] and small quantum systems [29], chains of coupled harmonic [31] or anharmonic [9,28,34] oscillators, models characterized by long-range [38,46] or disordered [41] interactions, as well as systems of coupled classical rotators [42,43,44,45]. In the case of a coupled XY nearest-neighbor-interacting rotator chain [44], the temperature dependence of the thermal conductance was well-fitted by a q-Gaussian distribution,

$$P_q\left(u\right)=P_0{exp}_q(-\beta u^2),$$

(2)

defined in terms of the q-exponential function,

$${exp}_q\left(u\right)={[1+(1-q)u]}_{+}^{1/(1-q)};\left({exp}_1\left(u\right)=exp\left(u\right)\right),$$

(3)

where $P_0\equiv P_q\left(0\right)$ and ${\left[y\right]}_{+}=y$, for $y>0$ (zero otherwise). The distribution in Equation (2) is very common in the context of nonextensive statistical mechanics [47], since it appears from the extremization of the generalized entropy, known as $S_q$, characterized by a real index q [48],

$$S_q=k\sum _{i=1}^W{p}_i\left({ln}_q\frac{1}{p_i}\right),$$

(4)

where we introduced the q-logarithm definition,

$${ln}_{q}u=\frac{u^{1-q}-1}{1-q};\left({ln}_{1}u=lnu\right).$$

(5)

Therefore, one recovers Boltzmann–Gibbs (BG) entropy,

$$S_{\mathrm{BG}}=-k\sum _{i=1}^W{p}_{i}ln{p}_i,$$

(6)

as ${lim}_{q\to 1}{S}_q=S_{\mathrm{BG}}$, whereas in the microcanonical ensemble, where all microstates present equal probability, $p_i=1/W$, Equation (4) becomes,

$$S_q=k{ln}_{q}W.$$

(7)

Above, the q-exponential function in Equation (3) appears precisely as the inverse function of the q-logarithm of Equation (5), i.e., ${exp}_q\left({ln}_{q}u\right)={ln}_q\left({exp}_q\left(u\right)\right)=u$.

Since the introduction of the entropy $S_q$ in Equation (4), a large amount of works appeared in the literature, defining generalized functions and distributions (see, e.g., Ref. [47]). In particular, a recent study based on superstatistics has found a stretched q-exponential probability distribution [49],

$$P_q\left(u\right)=P_0{exp}_q{(-\beta |u|}^{\eta })(0<\eta \le 1),$$

(8)

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, hence, Fourier’s law, in low-dimensional (or even finite-size) systems, like two-dimensional materials [14,15], silicon nanowires [17], carbon nanotubes [18], and low-dimensional nanoscale systems [19]. 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 the 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 [45], 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 Equation (2). In fact, Fourier’s law was validated in Ref. [45] by fitting the thermal conductance in terms of the functional form of Equation (8), with values of $\eta \left(d\right)>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,\cdots ,L-1$) interact by means of nearest-neighbor ferromagnetic couplings and evolve in time through molecular-dynamics numerical simulations. For this classical model, we specifically concentrate on the high-temperature limit, where there is no need for a spin wave approach, such as the Holstein–Primakoff quantum transformations. Our numerical data validate Fourier’s law, and similar to those of Ref. [45], the thermal conductance is well-fitted by the functional form of Equation (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 Section 2, we define the model and the numerical procedure; in Section 3, we present and discuss our results; in Section 4, we present our conclusions.

2. Materials and Methods

The one-dimensional classical inertial Heisenberg model, for a system of L-interacting rotators, is defined by the Hamiltonian,

$$\begin{array}{c}\hfill {\displaystyle \mathcal{H}=\frac{1}{2}\sum _{i=1}^L{\mathsf{\ell }}_i^2+\frac{1}{2}\sum _{\langle ij\rangle }\left(1-{\mathbf{S}}_i·{\mathbf{S}}_j\right),}\end{array}$$

(9)

where ${\mathsf{\ell }}_i\equiv ({\mathsf{\ell }}_{ix},{\mathsf{\ell }}_{iy},{\mathsf{\ell }}_{iz})$ and ${\mathbf{S}}_i\equiv (S_{ix},S_{iy},S_{iz})$ represent, respectively, continuously varying angular momenta and spin variables at each site of the linear chain, whereas ${\sum }_{\langle ij\rangle }$ denote summations over pairs of nearest-neighbor spins; herein, we set, without loss of generality, $k_{\mathrm{B}}$, moments of inertia, and ferromagnetic couplings, all equal to the unit. Moreover, spins present the unit norm, ${\mathbf{S}}_i^2=1$, and at each site, angular momentum ${\mathsf{\ell }}_i$ must be perpendicular to ${\mathbf{S}}_i$, yielding ${\mathsf{\ell }}_i·{\mathbf{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 to a system of coupled classical XY rotators, where canonical conjugate polar coordinates are commonly used [45], in the Heisenberg case, one often chooses Cartesian coordinates [50,51,52]. The reason for this is essentially technical, since in terms of spherical coordinates (more precisely, $\theta ,\varphi$, and their canonical conjugates ${\mathsf{\ell }}_{\theta },{\mathsf{\ell }}_{\varphi }$), a troublesome term $(1/{sin}^2\theta )$ appears in the corresponding equations of motion, leading to numerical difficulties [53,54]. However, some of the analytical results to be derived next recover those of the classical inertial XY model for ${\mathbf{S}}_i=(sin{\theta }_i,cos{\theta }_i,0)$ and ${\mathsf{\ell }}_i={\mathsf{\ell }}_i\widehat{\mathbf{z}}$.

It is important to mention that previous research on the thermal conductivity has been carried out for a classical one-dimensional Heisenberg spin model, by using Monte Carlo and Langevin numerical simulations [55], as well as for a classical one-dimensional spin-phonon system, through linear-response theory and the Green–Kubo formula [56]. These investigations did not take into account the kinetic contribution in Equation (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 Equation (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 obtain its thermal conductivity directly.

In order to carry out 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. Figure 1), whereas all remaining rotators ($i=2,\cdots ,L-1$) follow their usual equations of motion (see, e.g., Refs. [50,51,52]). In this way, one has for sites $i=2,\cdots ,L-1$,

$$\begin{array}{c}\begin{array}{cc}\hfill & {\dot{\mathbf{S}}}_i={\mathsf{\ell }}_i\times {\mathbf{S}}_i,\hfill \\ \hfill & {\dot{\mathsf{\ell }}}_i={\mathbf{S}}_i\times ({\mathbf{S}}_{i+1}+{\mathbf{S}}_{i-1}),\hfill \end{array}\hfill \end{array}$$

(10)

whereas the rotators at extremities follow standard Langevin dynamics,

$$\begin{array}{c}\begin{array}{cc}\hfill & {\dot{\mathsf{\ell }}}_1=-{\gamma }_h{\mathsf{\ell }}_1+{\mathbf{S}}_1\times {\mathbf{S}}_2+{\mathsf{\eta }}_h,\hfill \\ \hfill & {\dot{\mathsf{\ell }}}_L=-{\gamma }_l{\mathsf{\ell }}_L+{\mathbf{S}}_L\times {\mathbf{S}}_{L-1}+{\mathsf{\eta }}_l.\hfill \end{array}\hfill \end{array}$$

(11)

Above, ${\gamma }_h$ and ${\gamma }_l$ represent friction coefficients, whereas ${\mathsf{\eta }}_h$ and ${\mathsf{\eta }}_l$ denote independent three-dimensional vectors, ${\mathsf{\eta }}_h\equiv ({\eta }_{hx},{\eta }_{hy},{\eta }_{hz})$, ${\mathsf{\eta }}_l\equiv ({\eta }_{lx},{\eta }_{ly},{\eta }_{lz})$, where each Cartesian component stands for a Gaussian white noise with zero mean and correlated in time,

$$\begin{array}{c}\begin{array}{cc}\hfill & \langle {\eta }_{h\mu }\left(t\right)\rangle =\langle {\eta }_{l\mu }\left(t\right)\rangle =0,\hfill \\ \hfill & \hfill \langle {\eta }_{h\mu }\left(t\right){\eta }_{l\nu }\left(t^{\prime }\right)\rangle =\langle {\eta }_{h\mu }\left(t^{\prime }\right){\eta }_{l\nu }\left(t\right)\rangle =0,\hfill \\ \hfill & \langle {\eta }_{h\mu }\left(t\right){\eta }_{h\nu }\left(t^{\prime }\right)\rangle =2{\delta }_{\mu \nu }{\gamma }_h{T}_h\delta (t-t^{\prime }),\hfill \\ \hfill & \langle {\eta }_{l\mu }\left(t\right){\eta }_{l\nu }\left(t^{\prime }\right)\rangle =2{\delta }_{\mu \nu }{\gamma }_l{T}_l\delta (t-t^{\prime }),\hfill \end{array}\hfill \end{array}$$

(12)

with the indexes $\mu$ and $\nu$ denoting Cartesian components; from now on, we will set the friction coefficients ${\gamma }_h$ and ${\gamma }_l$ equal to the unit. One should mention that different types of thermostats have been used to investigate transport properties in systems out of equilibrium (see, e.g., Ref. [42] for an application of Nosé–Hoover thermostats to a system of interacting planar rotators); however, for the present Heisenberg chain, we found it more convenient to use standard Langevin thermostats, as defined above.

The condition of a constant norm for the spin variables yields

$$\begin{array}{c}\hfill \frac{d{S}_i}{dt}=\frac{d{\left({\mathbf{S}}_i·{\mathbf{S}}_i\right)}^{1/2}}{dt}=0\Rightarrow {\mathbf{S}}_i·{\dot{\mathbf{S}}}_i=0,\end{array}$$

(13)

which should be used together with ${\mathsf{\ell }}_i·{\mathbf{S}}_i=0$ in order to eliminate ${\ddot{\mathsf{\ell }}}_i$ and calculate ${\ddot{\mathbf{S}}}_i$ from Equations (10) and (11). For rotators at sites $i=2,\cdots ,L-1$, one has

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \hfill {\ddot{\mathbf{S}}}_i=({\mathbf{S}}_{i+1}+{\mathbf{S}}_{i-1})-\left[{\mathbf{S}}_i·({\mathbf{S}}_{i+1}+{\mathbf{S}}_{i-1})+{\dot{\mathbf{S}}}_i^2\right]{\mathbf{S}}_i,\hfill \end{array}\end{array}$$

(14)

whereas, for those at extremities,

$$\begin{array}{c}\begin{array}{cc}\hfill & {\ddot{\mathbf{S}}}_1=-{\dot{\mathbf{S}}}_1+{\mathbf{S}}_2-\left[{\mathbf{S}}_1·{\mathbf{S}}_2+{\dot{\mathbf{S}}}_1^2\right]{\mathbf{S}}_1+{\mathbf{S}}_1\times {\mathsf{\eta }}_h,\hfill \\ \hfill & {\ddot{\mathbf{S}}}_L=-{\dot{\mathbf{S}}}_L+{\mathbf{S}}_{L-1}-\left[{\mathbf{S}}_L·{\mathbf{S}}_{L-1}+{\dot{\mathbf{S}}}_L^2\right]{\mathbf{S}}_L+{\mathbf{S}}_L\times {\mathsf{\eta }}_l.\hfill \end{array}\hfill \end{array}$$

(15)

For the system illustrated in Figure 1, we will consider the temperatures of the heat baths differing by $2\epsilon$, with $\epsilon$ representing a positive dimensionless parameter; moreover, the temperature parameter $T=(T_h+T_l)/2$ will vary in a certain range of positive values. Equations (14) and (15) are transformed into first-order differential equations (e.g., by defining a new variable ${\mathbf{V}}_i\equiv {\dot{\mathbf{S}}}_i$) to be solved numerically through the velocity Verlet method [57,58], with a time step $dt=0.005$, for different lattice sizes L (please, see the Appendix A). The rotators at the bulk (i $=2,\cdots ,L-1$) follow a continuity equation,

$$\begin{array}{c}\hfill \frac{d{E}_i}{dt}=-(J_i-J_{i-1}),\end{array}$$

(16)

where

$$\begin{array}{c}\hfill {\displaystyle E_i=\frac{1}{2}{\mathsf{\ell }}_i^2+\frac{1}{2}\sum _{j=i\pm 1}\left(1-{\mathbf{S}}_i·{\mathbf{S}}_j\right),}\end{array}$$

(17)

so the stationary state is attained for $(d{E}_i/dt)=0$, i.e., $J_i=J_{i-1}$. The derivation is simple, since from Equation (13) and ${\mathsf{\ell }}_i·{\mathbf{S}}_i=0$, we have ${\dot{\mathbf{S}}}_i^2={\mathsf{\ell }}_i^2$, hence,

$$\frac{d}{dt}{E}_i={\dot{\mathbf{S}}}_i·{\ddot{\mathbf{S}}}_i-\frac{1}{2}\left[{\dot{\mathbf{S}}}_i·\left({\mathbf{S}}_{i+1}+{\mathbf{S}}_{i-1}\right)+{\mathbf{S}}_i·\left({\dot{\mathbf{S}}}_{i+1}+{\dot{\mathbf{S}}}_{i-1}\right)\right].$$

(18)

This equation, together with Equation (14), yields

$$\frac{d}{dt}{E}_i=\frac{1}{2}\left[{\dot{\mathbf{S}}}_i·\left({\mathbf{S}}_{i+1}+{\mathbf{S}}_{i-1}\right)-{\mathbf{S}}_i·\left({\dot{\mathbf{S}}}_{i+1}+{\dot{\mathbf{S}}}_{i-1}\right)\right]=0$$

(19)

at the stationary state. Data are obtained at stationary states, which, as usual, take longer to reach 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\simeq 0.15L$). In this way, we define an average heat flux as

$$\begin{array}{c}\hfill J\equiv \frac{1}{L-2p}\sum _{i=p+1}^{L-p}\langle J_i\rangle ,\end{array}$$

(20)

$$\begin{array}{c}\hfill J_i=\frac{1}{2}\left({\mathbf{S}}_i·{\dot{\mathbf{S}}}_{i+1}-{\mathbf{S}}_{i+1}·{\dot{\mathbf{S}}}_i\right),\end{array}$$

(21)

whereas $\langle ..\rangle$ denotes time and sample averages, which will be described next.

Let us emphasize that for ${\mathbf{S}}_i=(sin{\theta }_i,cos{\theta }_i,0)$ and ${\mathsf{\ell }}_i={\mathsf{\ell }}_i\widehat{\mathbf{z}}$, one recovers the expression for the heat flux of the classical inertial XY model, i.e., $J_i=\frac{1}{2}({\mathsf{\ell }}_i+{\mathsf{\ell }}_{i+1})sin({\theta }_i-{\theta }_{i+1})$ [45,59], 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 considered a transient of $5\times {10}^7$ time units to compute the averages $\langle J_i\rangle$ in Equation (20), and checked that this transient time was sufficient to fulfill the condition $J_i=J_{i-1}$ (within, at least, a three-decimal digits accuracy), for all values of L analyzed. After that, simulations were carried out for an additional interval of $2\times {10}^8$ time units (leading to a total time of $2.5\times {10}^8$ for each simulation). The interval $2\times {10}^8$ was divided into 80 equally spaced windows of $2.5\times {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 $\langle J_i\rangle$.

Using the results of Equation (20), one may calculate the thermal conductivity of Equation (1), and consequently, the thermal conductance,

$$\begin{array}{c}\hfill \sigma =\frac{J}{T_h-T_l}=\frac{J}{2T\epsilon }\equiv \frac{\kappa }{L}.\end{array}$$

(22)

In the next section, we present the results of both quantities, obtained from the numerical procedure described above.

3. Results

We simulate the system of Figure 1 for different lattice sizes, namely, $L=50,70,100,$ 140, considering the heat-bath temperatures differing by $2\epsilon$, with $\epsilon =0.125$. The temperature parameter $T=(T_h+T_l)/2$ varied in the interval $0<T\le 3.5$, capturing both low- and high-temperature regimes. The values of L ($L\ge 50$) were chosen adequately to guarantee that the thermal conductivity $\kappa$ did not present any dependence on the size L in the high-temperature regime, as expected.

In Figure 2, we present numerical data for the thermal conductivity Figure 2a and thermal conductance Figure 2b versus temperature (log–log representations) and different sizes L. The similar qualitative behaviors of the data displayed in both properties of Figure 2, for different values of L, evidence that the sizes considered in the present analysis ($L\ge 50$) are sufficiently large, in the sense that finite-size effects do not play a relevant role. In Figure 2a, we exhibit $\kappa (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\simeq 0.3$), as described next. (i) A low-temperature regime, where $\kappa$ depends on the size L, decreasing smoothly for increasing temperatures (L fixed). The plots of Figure 2a show that, in the limit $T\to 0$, an extrapolated value, $\kappa (L,0)\equiv {lim}_{T\to 0}\kappa (L,T)$, increases with L. Such a low-temperature increase with L has been observed in other one-dimensional models (see, e.g., Refs. [42,43,44,45]) and is reminiscent of the behavior expected for a chain of coupled classical harmonic oscillators. 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 $\kappa$ essentially does not depend on L (in the limit $L\to \infty$), as expected from Fourier’s law. Moreover, in this regime, one notices that $\kappa$ decreases with the temperature as it 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 Figure 2a indicate that the thermal conductivity becomes independent of the lattice size in the limit $L\to \infty$, scaling with the temperature as $\kappa \left(T\right)\sim T^{-2.25}$ at high temperatures. Therefore, the system becomes a thermal insulator at high temperatures, approaching this state according to $\kappa \left(T\right)\sim T^{-2.25}$. Despite the simplicity of the one-dimensional classical inertial Heisenberg model of Figure 1, the present results are very close to experimental verifications in some antiferromagnetic electrical insulators, such as the Heisenberg chain cuprates ${\mathrm{Sr}}_2{\mathrm{CuO}}_3$ and ${\mathrm{SrCuO}}_2$, for which the thermal conductivity is well-fitted by a $1/T^2$ law at high temperatures [4]. We should note that the one-dimensional Heisenberg model with nearest-neighbor ferromagnetic interactions, defined by the Hamiltonian of Equation (9), does not present an equilibrium phase transition, being characterized by a paramagnetic state for all temperatures $T>0$. In this case, one may perform the following transformations in the Hamiltonian of Equation (9), leaving it unaltered: $1/2\to -1/2$ (which incorporates the coupling constant), as well as ${\mathbf{S}}_j\to -{\mathbf{S}}_j$, keeping ${\mathbf{S}}_i$ unchanged. Consequently, the Hamiltonian of Equation (9) applies to antiferromagnetic systems at high temperatures, as well. 0

The same data of Figure 2a are exhibited in Figure 2b where we plot the thermal conductance $\sigma (L,T)=\kappa (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 $\kappa (L,0)$ scales as $\kappa (L,0)\sim L$, leading to $\sigma (L,0)\equiv {lim}_{T\to 0}\kappa (L,T)/L\simeq 0.5$. Such low-temperature results are in full agreement with those obtained in previous simulations of coupled classical XY rotators [42,43,44,45]. On the other hand, in the high-temperature regime, the thermal conductance presents a dependence on L, as expected.

In Figure 3, we exhibit the thermal-conductance data of Figure 2b in conveniently chosen variables, yielding a data collapse for all values of L considered. The full line essentially represents the form of Equation (8), so that one writes

$$\sigma (L,T)=A{exp}_q(-B{x}^{\eta }),$$

(23)

where $x=L^{0.475}T$, $q=2.28\pm 0.04$, $\eta =2.88\pm 0.04$, $A=0.492\pm 0.002$, and $B=0.33\pm 0.04$. Notice that this value of $\eta$ lies outside the range of what is commonly known as “stretched” [cf. Equation (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 nearest-neighbor-interacting classical XY rotators on d-dimensional lattices ($d=1,2,3$) [45], the thermal conductance was also fitted by the form of Equation (23), with values of $\eta \left(d\right)>2$. In particular, in the one-dimensional case, such a fitting was attained for $x=L^{0.3}T$, $q=1.7$, and $\eta =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. It is important to mention that the generalized forms in Equations (8) and (23) have been used in the literature for an appropriate description of a wide variety of physical phenomena, like velocity measurements in a turbulent Couette–Taylor flow [60], relaxation curves of RKKY spin glasses, such as CuMn and AuFe [61], cumulative distribution for the magnitude of earthquakes [62], and more recently, for the thermal conductance of a system of interacting XY rotators [45]. Moreover, its associated entropic form has been studied in detail in Ref. [49].

By defining the abscissa variable of Figure 3 in the general form $x=L^{\gamma (n,d)}T$, and using the q-exponential definition of Equation (3), the slope of the high-temperature part of the thermal-conductance data scales with L, as

$$\sigma \sim L^{-\left[\eta \right(n,d\left)\gamma \right(n,d\left)\right]/\left[q\right(n,d)-1]},$$

(24)

where we introduce the dependence $(n,d)$ on all indices. Since the thermal conductivity ($\kappa =L\sigma$) should not depend on the size L (in the limit $L\to \infty$), Fourier’s law becomes valid for

$$\frac{\eta (n,d)\gamma (n,d)}{q(n,d)-1}=1.$$

(25)

The data of Figure 3 lead to $\left[\eta \right(3,1\left)\gamma \right(3,1\left)\right]/\left[q\right(3,1)-1]=1.069\pm 0.083$, whereas those for XY rotators on d-dimensional lattices yield $1.0007,0.95$, and $0.93$, for $d=1,2$, and 3, respectively [45], indicating the validation of Fourier’s law for systems of coupled nearest-neighbor-interacting classical n-vector rotators, through the thermal conductance form of Equation (23).

Recently, similar analyses were carried out for an XY Hamiltonian with anisotropies, in such a way to approach the Ising model in particular limits [63]. All the results for the quantity in Equation (25), computed up to the moment, are summarized in

Table 1. Values of the ratio $\eta \gamma /(q-1)$ (highlighted in blue color) analyzed up to the moment: $n=1$ ($d=1$) [63], $n=2$ (dimensions $d=1,2,3$) [45], together with the present results for $n=3$ ($d=1$). In all cases studied, the limit of Equation (25) is numerically approached.

$\frac{\eta \gamma }{q-1}$	$d=1$ (linear chain)	$d=2$ (square lattice)	$d=3$ (simple cubic lattice)
$n=1$ (Ising ferromagnet)	1.0063 $q=1.65,\eta =1.94,\gamma =0.336$	-	-
$n=2$ ($XY$ ferromagnet)	1.0007 $q=1.7,\eta =2.335,\gamma =0.3$	0.95 $q=3.2,\eta =5.23,\gamma =0.4$	0.93 $q=3.5,\eta =5.42,\gamma =0.43$
$n=3$ (Heisenberg ferromagnet)	1.069 $q=2.28,\eta =2.88,\gamma =0.475$	-	-

, where one notices that finite-size effects play an important role in increasing dimensions, as expected.

4. Conclusions

We 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 the extremities of the chain were subjected to standard Langevin dynamics, whereas all remaining rotators ($i=2,\cdots ,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 $\kappa \left(T\right)$ and thermal conductance $\sigma (L,T)=\kappa \left(T\right)/L$ were computed, by defining $T=(T_h+T_l)/2$. The slope of the high-temperature part of thermal-conductance data scales with the system size was $\sigma \sim L^{-1.069}$, indicating that in the limit $L\to \infty$, one should obtain a thermal conductivity independent of L. Indeed, in this limit, we found $\kappa \left(T\right)\sim T^{-2.25}$ for high temperatures. The thermal-conductance data were well-fitted by the function $\sigma (L,T)=A{exp}_q(-B{x}^{\eta })$, typical of nonextensive statistical mechanics, where A and B are constants, $x=L^{0.475}T$, $q=2.28\pm 0.04$, and $\eta =2.88\pm 0.04$. This fitting augments the applicability of such a function, which has been used for describing several physical phenomena in the literature, like velocity measurements in a turbulent Couette–Taylor flow [60], relaxation curves of RKKY spin glasses [61], cumulative distribution for the magnitude of earthquakes [62], and thermal conductance of a system of interacting XY rotators [45]. Since the value of $\eta$ found herein lies outside the range of what is commonly known as “stretched” $(0<\eta \le 1)$, herein, we refer to this fitting function of a “shrinked” q-exponential. The present results reinforce those obtained recently for XY rotators on d-dimensional lattices [45], 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 $\eta (n,d)$ presenting a dependence on the number of spin components and lattice dimension.

Despite the simplicity of the model considered herein, the results for the thermal thermal conductivity at high temperatures ($\kappa \left(T\right)\sim T^{-2.25}$) are very close to experimental verifications in some antiferromagnetic electrical insulators, such as the Heisenberg chain cuprates ${\mathrm{Sr}}_2{\mathrm{CuO}}_3$ and ${\mathrm{SrCuO}}_2$, for which thermal conductivity is well-fitted by a $1/T^2$ law at high temperatures [4]. At equilibrium, the present model exhibits a paramagnetic state for all temperatures, so that its Hamiltonian may be shown to cover both ferromagnetic and antiferromagnetic systems. The present results show that even for models exhibiting simple equilibrium properties, one may have out-of-equilibrium regimes characterized by transport properties typical of nonextensive statistical mechanics, like the ones found herein. Since nonextensive statistical mechanics have 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.

In summary, we demonstrated that (i) for the classical one-dimensional inertial ferromagnetic Heisenberg model, the (macroscopic) Fourier-law is validated from (microscopic) first principles, i.e., the temperature-dependent thermal conductivity is, in the high-temperature regime, finite and independent of the system size (the low-temperature regime is to be handled within a quantum grounding, which is out of the goal of the present paper); (ii) For all temperatures and sizes, the thermal conductivity appears to be consistent with q-statistics since it can be neatly collapsed within a shrunken q-exponential form; (iii) within this shrunken q-exponential form, a single universal condition, namely $\frac{\eta (n,d)\gamma (n,d)}{q(n,d)-1}=1$, validates the Fourier law for the n-vector models for $n=1,2,3$, which constitutes a numerical indication that this centennial macroscopic law is possibly valid for all values of $(n,d)$, where $n\to \infty$ and $n\to 0$ correspond to the spherical model and ‘self-avoiding walk’, respectively. It is not our present aim to review the rich existing literature on the validity of the Fourier law within diverse classical and quantum approaches, but we rather restrict our focus to analytical and numerical first-principle approaches of classical systems that are similar to the present one.