Thermal transport in the Fermi-Pasta-Ulam model with long-range interactions

TL;DR

This study investigates how long-range interactions affect thermal transport in a one-dimensional Fermi-Pasta-Ulam model, revealing a non-monotonic conductivity dependence on interaction decay and nearly ballistic transport at a specific decay exponent.

Contribution

It introduces a detailed numerical analysis of thermal transport in a long-range interacting FPU model, highlighting the non-monotonic behavior of conductivity and its relation to system integrability.

Findings

Conductivity peaks at decay exponent δ=2.0.

At δ=2.0, conductivity diverges linearly with system size.

Nearly ballistic transport observed at δ=2.0.

Abstract

We study the thermal transport properties of the one dimensional Fermi-Pasta-Ulam model ($\beta$-type) with long-range interactions. The strength of the long-range interaction decreases with the (shortest) distance between the lattice sites as ${distance}^{-\delta}$, where $\delta \ge 0$.Two Langevin heat baths at unequal temperatures are connected to the ends of the one dimensional lattice via short-range harmonic interactions that drive the system away from thermal equilibrium. In the nonequilibrium steady state the heat current, thermal conductivity and temperature profiles are computed by solving the equations of motion numerically. It is found that the conductivity $\kappa$ has an interesting non-monotonic dependence with $\delta$ with a maximum at $\delta = 2.0$ for this model. Moreover, at $\delta = 2.0$, $\kappa$ diverges almost linearly with system size $N$ and the temperature profile has a negligible slope, as one expects in ballistic transport for an integrable system. We demonstrate that the non-monotonic behavior of the conductivity and the nearly ballistic thermal transport at $\delta = 2.0$ obtained under nonequilibrium conditions can be explained consistently by studying the variation of largest Lyapunov exponent $\lambda_{max}$ with $\delta$, and excess energy diffusion in the equilibrium microcanonical system.

Equations

$$\mathcal{H}=\sum _i\frac{p_i^2}{2m}+\frac{1}{2}\sum _i(x_{i+1}-x_i{)}^2+\frac{1}{4\tilde{N}}\sum _{i,j\ne i}\frac{(x_j-x_i{)}^4}{{|i-j|}^{\delta }}.$$

$$\tilde{N}=\frac{1}{N}\sum _i\sum _{j\ne i}{|i-j|}^{-\delta }.$$

$$J_i=-⟨\frac{1}{2}({\dot{x}}_{i+i}+{\dot{x}}_i)(x_{i+1}-x_i)⟩,$$

$$T_i=⟨m{\dot{x}}_i^2⟩$$

$$\kappa =\frac{JN}{T_L-T_R}.$$

$${\rho }_E(r,t)=\frac{⟨\Delta E_j(t)\Delta E_i(0)⟩}{⟨\Delta E_i(0)\Delta E_i(0)⟩}+\frac{1}{N_b-1},$$

$$⟨\Delta r^2(t){⟩}_{{}_E}=\sum _{r=-N/2}^{N/2}{r}^2{\rho }_{{}_E}(r,t).$$

Full text

I Introduction

Long-range interactions are ubiquitous in nature, the most common examples being gravitation and Coulomb interaction. Over the last several decades long-ranged systems have been extensively studied and it is now known that such systems have a very rich, and often intriguing, thermo-statistical behavior (for reviews on long-range interacting systems see LR_Levin ; LR_Mukamel ; LR_Ruffo ). These systems exhibit breakdown of ergodicity, suppression of chaos, inequivalence of ensembles, long lived non-Gaussian quasi-stationary states, phase transitions in one dimension, negative specific heat, to name a few, and therefore often possess properties that deviate fantastically from “well behaved” short-ranged equilibrium systems. Needless to say, even after a lot of effort an in-depth understanding of long-range interacting systems seems to be lacking.

In non-equilibrium statistical physics, a profusely investigated topic is thermal transport, particularly in low dimensional microscopic models HT_Dhar ; HT_LLP . A flurry of research started after it was discovered that many simple one-dimensional (1D) models violate the celebrated Fourier’s law of heat conduction $J=-\kappa\nabla T$, where $\kappa\sim N^{0}$ is the constant of thermal conductivity ($N$ being the system size), and the heat current $J$ is proportional to the temperature gradient $\nabla T$ which is assumed to be small. Subsequently a large number of works demonstrated that $\kappa$ is actually not a constant independent of $N$ but often diverges as $\kappa\sim N^{\alpha}$, $0<\alpha<1$ HT_Dhar . These models are frequently momentum-conserving systems i.e., without external pinning potentials (a notable exception is the coupled rotor model Rotor ). One such model is the celebrated Fermi-Pasta-Ulam (FPU) model that has been studied extensively in the last 50 years in statistical mechanics, chaos theory, nonlinear dynamics and several other contexts FPUoriginal ; FPU50yr ; FPUparadox ; FPUproblem .

In this paper we wish to study thermal transport properties of the FPU model in 1D which has been appropriately modified to include long-range interactions. Transport studies of long-ranged systems are few (see LRDiode ; LRXYTrans ; LRRect ; LRTrans for some recent works) but are extremely important. For example, studies of thermal rectification ThDiode in short-ranged systems show that rectification generally ceases to exist as one approaches the thermodynamic limit Size . This raises a question on their applicability in the fabrication of real thermal diodes which are thought to be one of the crucial components of phononics Phononics . In LRRect ; LRDiode , the authors studied thermal rectification properties of a long-ranged nonlinear lattice model and made an important discovery that rectification is enhanced and also possibly survives even in the thermodynamic limit in the presence of long-range interactions. Thus, besides the purely theoretical interest of understanding the physics of long-ranged interacting systems and low-dimensional thermal transport, these studies also are technologically important. This is more so since presently it is possible to fabricate actual materials with long-range interactions, such as the Coulomb crystal crystal , the Ising pyrochlore magnets ${\rm Dy_{2}Ti_{2}O_{7}}$ and ${\rm Ho_{2}Ti_{2}O_{7}}$ Mat1 ; Mat2 and Permalloy nanomagnets Mat3 with long-range magnetic interactions.

The remainder of the paper is organized as follows. In Sec. II we describe the Fermi-Pasta-Ulam model with long-range interactions and the numerical scheme that has been employed to study its transport properties. We describe the numerical computation of our quantities of interest such as the heat current, conductivity, and temperature profiles. In Sec. III, we present the result of our non-equilibrium molecular dynamics simulation. Next, in Sec. IV, the equilibrium (microcanonical) version of the model is studied and by computing the largest Lyapunov exponent and spatiotemporal excess energy correlation function we explain the thermal transport features of the model under non-equilibrium conditions. We finally summarize our main results in Sec. V and conclude with a discussion.

II The long-ranged FPU model with heat baths

We consider a one-dimensional lattice of $N$ anharmonic oscillators, each with mass $m$, displacement $x_{i}$ and momentum $p_{i}$ ($1\leq i\leq N$) which interact with each other by long-range interaction. The Hamiltonian of our long-ranged Fermi-Pasta-Ulam (LR-FPU) model is FPU1D ; BagchiTsallis

[TABLE]

Thus the interaction potential has two parts - the second term on the right in Eq. (1) is the nearest neighbor (harmonic) interaction whereas, the third term corresponds to the long-ranged anharmonic interactions. The strength of the long-range interaction decays as a power-law $|i-j|^{-\delta}$ ($\delta\geq 0$) with the (shortest) distance between the lattice sites $i$ and $j$. Thus $\delta$ sets the range of the interaction; for $\delta=0$ we have a mean-field scenario where every oscillator interacts with all others with equal strength, whereas for $\delta\to\infty$ the third term reduces to nearest neighbor anharmonic interactions. Note that the prefactor $1/\tilde{N}$ makes the Hamiltonian Eq. (1) extensive for all values of $\delta$ and is given as

[TABLE]

For the mean-field ($\delta=0$) scenario $\tilde{N}=N$ which is the so called Kac scaling factor. Thus, inside the shaded (green) box in Fig. 1, every oscillator interacts with all others obeying Eq. (1) with interaction strengths depending on the value of $\delta$.

In order to study thermal transport in this kind of model we do the following alterations, as in Ref. LRDiode . We attach $n$ additional oscillators (we set $n=2$, as shown in Fig. 1) to both sides of the one dimensional chain and connect Langevin heat baths HT_Dhar at the two ends. We refer to these oscillator sites, having displacements denoted by $x_{-1},x_{0}$ at the left and $x_{N+1},x_{N+2}$ at the right in Fig. 1, as the lead sites (or simply leads) that connect the main long-ranged system under investigation to the heat baths. The leads interact with each other and the LR-FPU system harmonically i.e., $V(x)=x^{2}/2$ (spring constant set to unity). The two heat baths are set to unequal temperatures $T_{L}$ and $T_{R}$ (we set $T_{L}>T_{R}$ always) which will drive a heat current though the long-ranged system.

The advantage of attaching the lead sites at the two ends is the ease with which one can measure the heat current $J$. Since the interactions $\mathcal{H}$ in the main long-ranged system is quite complicated, computing the current in the bulk of the system can be cumbersome. However exploiting the fact that when the system (including the oscillators at the lead sites) attains a nonequilibrium steady state (NESS), the current $J_{i}$ at the $i$th site on the lattice is constant and independent of the lattice site index $i$ i.e., $J_{i}=J=constant$. Thus, one can measure the steady state current $J$ at the lead sites to compute the heat current propagating through the main system. Since the lead site oscillators interact with a simple potential $V(x)=x^{2}/2$, the current at these sites has a simple expression

[TABLE]

where $i$ is the site index of one of the leads. The steady state local temperature profile can be computed using the equipartition theorem (with Boltzmann constant $k_{B}=1$)

[TABLE]

and assuming that the long-ranged system achieves local thermal equilibrium at temperature $T_{i}$. The finite size conductivity $\kappa$ for a system of size $N$ is computed using the Fourier’s law

[TABLE]

To numerically integrate the equations of motion, we employ the second order velocity-Verlet algorithm with small time-step $\sim 10^{-2}$. For long-ranged systems, a naive force calculation is extremely expensive computationally because of the $\mathcal{O}(N^{2})$ operations that have to be performed at each time step. This part can be accelerated by a technique which exploits the convolution theorem and using efficient fast Fourier transformations (see SpinIce for implementation algorithms in periodic and open spin systems with dipolar interactions). This drastically reduces the computation time with only $\mathcal{O}(NlogN)$ operations to perform at each iteration.

Starting for random initial conditions ($x_{i}$’s are chosen from a uniform distribution and $p_{i}$’s from a Gaussian distribution, both centered at zero) we evolve the system by integrating the equations of motion until a NESS is attained. Thereafter, we compute the heat current $J$, local temperature $T_{i}$ profiles, and conductivity $\kappa$ using Eqs. (3), (4), and (5), respectively. We present our results obtained from the nonequilibrium simulation in the next section.

III Results from nonequilibrium molecular dynamics

We choose the Langevin heat baths to have temperatures $T_{L}=T_{0}+{\Delta T}/2$ and $T_{R}=T_{0}-{\Delta T}/2$ with small temperature difference $\Delta T$ between the two ends; thus the average temperature of the system is $\frac{1}{2}{(T_{L}+T_{R})}=T_{0}$. Unless mentioned otherwise, we set $T_{0}=1.0$ and $\Delta T=0.2$ ($T_{L}=1.1$ and $T_{R}=0.9$) for our simulations. We employ free boundary conditions for the oscillators at the ends, although we expect our results to remain practically unchanged even with fixed boundary conditions.

We evolve the system for a large number of iterations ($\sim 10^{6}$) and monitor the heat currents entering the system through the left $J_{L}$ and leaving the system from the right $J_{R}$ (as in Fig. 1). When a NESS is achieved one should obtain $J_{L}=J_{R}$ in magnitude; $J_{L},J_{R}$ are obtained using Eq. (3). As an additional check for stationarity, we also computed the velocity distribution of the individual oscillators and obtained a distribution that remains practically unaltered at large times and, as anticipated, fits to a Gaussian distribution. We thereafter compute the heat current as $J=\frac{1}{2}{(J_{L}+J_{R})}$, temperature profile $T_{i}$ using Eq. 4 (with $m=1$), and conductivity $\kappa$ using Eq. (5), typically for $\sim 10^{8}-10^{9}$ time iterations depending on the value of $N$ chosen.

In Fig. 2, we show the thermal conductivity $\kappa$ ($\Delta T=0.2$ is a constant) of the LR-FPU model as a function of $\delta$ for system sizes $N=1024,2048,4096,8192$ It is found that the conductivity $\kappa$ exhibits an interesting non-monotonic behavior with $\delta$. For small values of $\delta<2$, $\kappa$ increases monotonically, attains a maximum value at $\delta=2.0$ and thereafter decreases monotonically. Thus very long-ranged ($\delta\to 0$) as well as nearest neighbor ($\delta\to\infty$) systems have a low thermal conductivity and maximum conductivity is obtained for an optimal $\delta=2.0$. Thus, by tuning the long-range parameter $\delta$, one can manipulate the conductivity of the system over a wide range (roughly two orders of magnitude for $N=8192$ in Fig. 2). In the inset of Fig. 2 the heat current $J=\kappa\Delta T/N$ is shown as a function of $\delta$. This is, in some sense, counter-intuitive since the range (or strength) of the long-range interaction decreases smoothly with increasing $\delta$ and thus naively one would expect the heat current to decrease monotonically with $\delta$. A non-monotonicity at $\delta=1$ would also have been easy to understand, since for $\delta$ equal to the embedding dimension, these models go from being a true long-ranged system to a short (finite) ranged system LR_Ruffo ; LR_Levin . This in many cases, such as in the LR-FPU FPU1D 111The numerical result of FPU1D in the “weak chaos” regime at large $N$ exhibits an unexplained saturation effect, (see Fig. 1a in FPU1D , which was also later pointed out in Fig. 5(b) of BagchiTsallis ). This makes the result in FPU1D less reliable for large $N$ and small $\delta$. Also, here we perform our simulations for very different parameter values than what are used in FPU1D ., quartic LR-FPU BagchiTsallis , and LR-XY (coupled planar rotors) XYChaos ; LLExpXY1 models, is marked by a ‘transition’ from weak chaos to strong chaos at $\delta=1$ (for 1D microcanonical systems). Instead, in the present case, it is quite intriguing that $\delta=2.0$ appears to be a special value having maximum conductivity $\kappa$.

Next, in Fig. 3(a) we present the conductivity $\kappa$ as a function of the system size $N$ for different values of $\delta$. We find that for all values of $\delta$ (computed up to very large values of $N=16384$) the conductivity $\kappa$ shows a power-law divergence $\kappa\sim N^{\alpha}$ at large $N$, as it should for momentum conserving models. However $\alpha$ changes quite drastically (and nonmonotonically) with $\delta$, as can be clearly seen in Fig. 3(a). Surprisingly enough, corresponding to $\delta=2.0$, the conductivity diverges almost linearly ($\alpha\approx 0.98$) with the system size $N$ akin to ballistic heat transport in integrable models e.g., the coupled harmonic oscillators Harmonic .

For $\delta=2.0$, we have checked the $N$-divergence of $\kappa$ with two other average temperatures $T_{0}=0.5,2.0$ (and same $\Delta T=0.2$) and found that the exponent remains practically unaltered as shown in Fig. 3(b). We should mention that the size dependence of $\kappa$ for the LR-FPU model seems to be very different from the results for the LR-XY model studied recently LRXYTrans . As an example, they obtain a thermal insulator i.e., $\kappa(N)\to 0$ as $N\to\infty$ for $\delta=0$ (Fig. 2 in LRXYTrans ), unlike what we obtain here.

In Fig. 4, $\alpha$ is shown as a function of $\delta$ and an interesting picture emerges. We find that the $\alpha-\delta$ plane in Fig. 4 can be split into the following regions:

(i) A region of strong long-range interaction for $0<\delta<1$. This region probably can be further split into two sub-regimes, namely $0\leq\delta<1/2$ where $\alpha\approx 0.66$ and $1/2<\delta<1$ for which $\alpha$ shows a decreasing behavior. A similar suggestion has been put forward in some works LLExpXY1 ; dby2 since these two sub-regimes often exhibit different properties.

(ii) An intermediate region of weak long-range interaction for $1<\delta<2$ for which the $\alpha$ varies quite drastically from $0.3\lesssim\alpha\lesssim 1$.

(iii) A region of short-range (or finite range but not nearest neighbor) interactions for $\delta>2$ where $\alpha$ monotonically decreases from $\alpha\approx 1$ to $\alpha\approx 0.4$; the value $\alpha\approx 0.4$ for large $\delta\approx 5.0$ lies well in the range $0.3<\alpha<0.5$ and has been obtained for the nearest neighbor FPU model TwoFifth in the past.

Actually, such a picture has been proposed for both classical XYCelia3Regions and quantum ThreeRegimesQnt long-ranged interacting models (also see LR_Mukamel ). However, one needs more detailed simulation in each region to ascertain this accurately.

Next, the local temperature profile $T_{i}$ is shown for different values of $\delta$ and $N=2048$ in Fig. 5. For $\delta=0.0$, the temperature profile is highly nonlinear and becomes less nonlinear as $\delta$ is increased, except for $\delta=2.0$, for which the temperature profile is nearly flat (in fact, it has a small slope in the $wrong$ sense for reasons unclear to us). The curve marked ‘equilibrium’ in Fig. 5 corresponds to $\Delta T=0$ and $\delta=2.0$. This is to demonstrate that the long-range interacting system equilibrates properly ($T_{i}=T_{0}=1.0$) in absence of external thermal drive at $\delta=2.0$. Evidently, the local slope (say, around $i=N/2$) of the temperature profile also changes non-monotonically as $\delta$ is varied with a nearly flat profile for $\delta=2.0$, again similar to the coupled harmonic oscillators.

To summarize the results presented in this section, we find that $\delta=2.0$ has a very special behavior which leads to a non-monotonic variation of various quantities such as the conductivity, the slope of the temperature profiles in the bulk of the system, and the divergence exponent $\alpha$. For $\delta=2.0$, our long-range interacting nonlinear system exhibits transport properties similar to an integrable model. In the next section, we provide additional numerical evidence to demonstrate that the results obtained here are physical and not an artifact of the nonequilibrium setup e.g., the model of heat baths, boundary conditions, interaction potential of the leads, etc.

IV Results from equilibrium molecular dynamics

In order to understand the results obtained in the previous section, we next study the equilibrium (microcanonical) version of the LR-FPU model using equilibrium molecular dynamics simulation. We remove the baths (along with the leads) and close the 1D chain to obtain a periodic lattice where each oscillator interacts with all the others according to Eq. (1). The total energy of the system is a conserved quantity now, which is achieved by rescaling the momenta $p_{i}$ appropriately after they are assigned randomly at $t=0$. We choose the energy density $u\equiv\langle\mathcal{H}\rangle/N=1.0$ which is roughly the same as the average energy of each oscillator at $T_{0}=1.0$ (the average temperature for which we performed the nonequilibrium simulations).

First, we compute the largest Lyapunov exponent $\lambda_{max}$ as a function of the parameter $\delta$. The numerical computation of $\lambda_{max}$ here has been performed using the algorithm by Benettin et al. Benettin . Note that, the larger the value of $\lambda_{max}$ the greater is the chaoticity of the model and vice-versa. For integrable systems, such as the harmonic oscillators, $\lambda_{max}=0$.

We show the largest Lyapunov exponent $\lambda_{max}$ as a function of the parameter $\delta$ in Fig. 6. As can be clearly seen there is a nonmonotonicity in the behavior of $\lambda_{max}$ exactly at $\delta=2.0$ for which we have obtained nonmonotonic behavior for $\kappa$ from nonequilibrium simulations (Fig. 2). For fixed $N$, $\lambda_{max}$ has a minimum at $\delta=2.0$, implying a suppression of chaos. We speculate that this lack of chaoticity pushes the nonlinear long-ranged system towards an integrable limit which, in turn, is responsible for the nearly ballistic transport features ($\kappa\sim N$ and negligible temperature gradient) at $\delta=2.0$ under nonequilibrium conditions. Since in ballistic transport the heat carriers propagate with minimum scattering (least resistance), the thermal conductivity is a maximum at $\delta=2.0$, as was obtained in the previous section.

Let us extend the connection between $\lambda_{max}$ and $\kappa$ further. In Fig. 7, we have shown $\lambda_{max}^{-1}$ (obtained from Fig. 6) and $\kappa$ (from Fig. 2) by shifting and rescaling these two quantities as $(\lambda^{-1}_{max}-\theta_{1})$ and $\theta_{2}(\kappa\Delta T-\theta_{3})$; $\theta$’s are real numbers. It can be seen that in both the regions $0\leq\delta<2$ and $\delta>2$, the two quantities follow each other approximately. Thus using the equilibrium results ($\lambda_{max}~{}vs.~{}\delta$) one can understand the nonequilibrium results ($\kappa~{}vs.~{}\delta$) qualitatively and, to an appreciable extent, quantitatively.

Another set of numerical experiments can be performed by exploiting the connection between thermal transport and energy diffusion since both essentially describe the same microscopic process. It is well established that for normal and ballistic heat transport the corresponding energy diffusion should also be normal and ballistic respectively. In recent times a lot of effort has been dedicated to extend such a connection to the regime of anomalous heat transport (see NJPSpread ; Espread , and references therein). In the following, we study energy diffusion in our equilibrium model by computing the spatiotemporal excess energy correlation function using the relation (for a microcanonical system) FC ; Methods

[TABLE]

where we coarse grain the lattice as $N_{b}=N/b$, $r=(i-j)b$ and $b$ is the number of sites (oscillators) in each bin and $\Delta E_{k}=E_{k}-\langle E_{k}\rangle$ is the excess energy of the $k-th$ bin ($1\leq k\leq N_{b}$) Methods .

For an equilibrated system, we calculate $\rho_{{}_{E}}(r,t)$ as a function of $\delta$ for a fixed $N=2048$ (we set $b=4$) for different values of time after the system has attained equilibrium. The evolution of $\rho_{E}(r,t)$ at times $t=1000,1200,1400$ for different values of $\delta$ is shown in Fig. 8. Note that, quite surprisingly, the speed of propagation is slower for $\delta=0.0$ ($\rho_{{}_{E}}(r,t)$ is quite localized around $x=0$) than for, say, $\delta=5.0$ (for which the long-range model approaches the nearest neighbor FPU model).

The spatiotemporal distributions $\rho_{{}_{E}}(r,t)$ for different times $t$ can be collapsed by rescaling the curves as $t^{\gamma}\rho_{{}_{E}}(rt^{-\gamma},t)$ where $\gamma<1/2,\gamma=1/2,\gamma>1/2,$ and $\gamma=1$ imply sub-diffusion, normal diffusion, super-diffusion and ballistic spreading respectively. The exponent $\gamma$ is related to the exponent $\beta$, governing the temporal variation of the mean square deviation (MSD) $\langle\Delta r^{2}(t)\rangle_{{}_{E}}\sim t^{\beta}$, as $\beta=2\gamma$. Thus when heat transport is ballistic $\beta=2$.

From Fig. 9(a) we find that for $\gamma=1.0$ one obtains an excellent collapse of the rescaled functions $\rho_{{}_{E}}(r,t)$ shown in the inset for different times $t=800,1000,1200,1400$. Thus $\beta=2\gamma=2.0$ and therefore one should expect a ballistic spread of energy. In Fig. 9(b) we exhibit the function $\rho_{{}_{E}}(r=0,t)$ that demonstrate that such a scaling (with $\gamma=1$) persists up to the largest times $t=2000$ that we have simulated. We also compute the MSD of the excess energy distribution NJPSpread

[TABLE]

The variation of MSD with time $t$ is particularly useful to quantify the speed and identify the nature of energy diffusion. This is shown in Fig. 9(c) for the LR-FPU model at $\delta=2.0$. As can be seen, we obtain $\langle\Delta r^{2}(t)\rangle_{{}_{E}}\sim t^{\beta}$ with $\beta=2.0$ indicating ballistic propagation. All these results strongly suggest that the ballistic-like heat transport (marked by the nearly linear divergence of the conductivity with system size $N$) as was seen in Fig. 3 for $\delta=2.0$ is a physical property of our model. For other values of $\delta$ depicted in Fig. 8, we estimate that $\gamma$ is in the range $3/5<\gamma<3/4$, thus implying super-diffusion and therefore violation of Fourier’s law in nonequilibrium, as we expect from Fig. 3(a). One can also compute the exponents $\beta$ and relate it to the $\alpha$ using the so called connection relations Rel1 ; Rel2 to verify if they hold. A recent study however argued that such a connection between thermal transport and energy diffusion may not exist in general for anomalous transport NoConn .

V Conclusion

In summary, we have studied thermal transport properties of the Fermi-Past-Ulam model in the presence of long-ranged interactions that decay with (shortest) distance between the lattice sites as a power-law $|i-j|^{-\delta}$. By attaching two Langevin heat baths at the two ends of the open chain we compute the heat current, conductivity, and the temperature profiles of the long-ranged FPU system for different values of the exponent $\delta\geq 0$ using large-scale nonequilibrium molecular dynamics simulations. We find that $\kappa$ exhibits an interesting non-monotonic $\delta$ dependence with a maximum at $\delta=2.0$. From the system size dependence of $\kappa$, we obtain a power-law divergence $\kappa\sim N^{\alpha}$, where $1/3\lesssim\alpha\lesssim 1$, for all values of $\delta$ (violation of Fourier’s law); for $\delta=2.0$ we obtain $\alpha\approx 1$ and the temperature profile has a negligible (wrong) slope.

In order to explain these results, we next look into the equilibrium dynamics of the long-ranged model (by removing the heat baths). We compute the largest Lyapunov exponent $\lambda_{max}$ as a function of $\delta$ and obtain a nonmonotonic dependence that enabled us to understand the $\delta$ variation of the conductivity $\kappa$, both qualitatively and, to some extent, quantitatively. The near ballistic divergence of $\kappa$ at $\delta=2.0$ is explained by the suppression of chaos for the same $\delta$ value in the equilibrium simulations. Thus, although the range of interaction changes monotonically with the long-range parameter $\delta$, it is intriguing to see non-monotonicity in all the quantities we have computed namely, $\kappa$, $\alpha$, and $\lambda_{max}$.

From the computation of the spatiotemporal excess energy correlation function, we confirmed the almost linear divergence of $\kappa$ for $\delta=2.0$, which corresponds to a ballistic energy diffusion in the equilibrium model. Thus for $\delta=2.0$, from both the nonequilibrium and equilibrium simulations, we find that the nonlinear long-ranged interacting FPU model behaves similar to a nearly integrable system NIM . This is quite exciting, since integrable nonlinear long-ranged models are not encountered very frequently (one such example is the many-body Calogero-Moser system CM which is nonlinear, long-ranged, and exactly integrable) and therefore this needs to be explored further. At this point, it is also quite tempting to recall the connection between the Fermi-Pasta-Ulam model and the nonlinear yet integrable Toda model Toda3 , both with nearest-neighbor interactions.

Note that we cannot completely rule out finite size effects in the results presented here which therefore may change as $N$ is increased further (for larger $N>16384$, performing numerical calculations become computationally impractical). Finite size effects have plagued several works even for the short-ranged FPU model in the past and in recent times FiniteSize . Nevertheless, we believe that the results reported here, besides their academic interest, will also be important from the experimental perspective where one deals with finite sized systems at finite temperatures.

There are a number of possible open questions that we believe are important for this model (or similar models) but have not been addressed in this paper. This includes analyzing the temperature dependences of conductivity, the generality of these results by studying other classes of interaction potentials, the presence of disorder and on-site potentials, the existence of discrete breathers that often affect heat transport DB3 , extension to higher dimensions etc. We plan to address these issues in the future. Of course any analytical result along these lines will be extremely desirable to obtain. We hope that the results presented here will encourage further investigations of the heat transport properties of long-range interacting systems in the future.

Acknowledgments: The author gratefully acknowledges fruitful discussions with A. Dhar. This work has been funded by the John Templeton Foundation.