Non-linear Fluctuating Hydrodynamics with many conserved fields: the case of a 3D anharmonic chain

TL;DR

This paper extends non-linear fluctuating hydrodynamics to a 3D anharmonic chain with multiple conserved fields, analyzing mode behaviors and comparing theoretical predictions with molecular dynamics simulations at high temperatures.

Contribution

It develops a 3D model within NLFHT, identifying distinct diffusive and superdiffusive behaviors of modes, and provides detailed correlation analyses.

Findings

Transverse modes are diffusive.

Heat and longitudinal modes are superdiffusive.

Good agreement between theory and simulations at high temperatures.

Abstract

We propose a model for a chain vibrating in 3 dimensions (3D), with first neighbors anharmonic interatomic potential which depends on their distance, and subjected to an external tension. In the framework of the non-linear fluctuating hydrodynamic theory (NLFHT), that was successfully applied to 1D chains, we obtain a heat mode, two longitudinal and four transverse sound modes. We compute their spatio-temporal correlations comparing the theoretical results with molecular dynamics simulations, finding a good agreement for high temperatures. We find that the transverse sound modes behave diffusively, meanwhile the heat and longitudinal sound modes superdiffusively, exploring their possible scaling functions and characteristic exponents.

Figures

Tables

${\varphi }_1$	Heat mode
${\varphi }_2$	Longitudinal mode	$x+$
${\varphi }_3$	Longitudinal mode	$x-$
${\varphi }_4$	Transverse mode	$y+$
${\varphi }_5$	Transverse mode	$y-$
${\varphi }_6$	Transverse mode	$z+$
${\varphi }_7$	Transverse mode	$z-$

Table 2. Table 1: Theoretical and numerical velocities of longitudinal and transverse sound mode correlations.

Set	$c_x$ Theo	$c_x$ Num	$c_{⟂}$ Theo	$c_{⟂}$ Num
A	1.12923	1.17 $\pm$ 0.01	0.608417	0.61 $\pm$ 0.01
B	1.24087	1.37 $\pm$ 0.01	0.661516	0.66 $\pm$ 0.01
C	2.05064	2.07 $\pm$ 0.02	1.40817	1.41 $\pm$ 0.01
D	2.25594	2.30 $\pm$ 0.02	1.45634	1.46 $\pm$ 0.01

Table 3. Table 2: Scaling exponents for the heat, longitudinal and transverse sound modes.

Set	${\delta }_1$	${\delta }_2$	${\delta }_4$
A)	0.74 $\pm$ 0.02	0.70 $\pm$ 0.02	0.46 $\pm$ 0.01
B)	0.78 $\pm$ 0.02	0.70 $\pm$ 0.02	0.46 $\pm$ 0.01
C)	0.60 $\pm$ 0.02	0.69 $\pm$ 0.02	0.50 $\pm$ 0.01
D)	0.62 $\pm$ 0.02	0.69 $\pm$ 0.02	0.50 $\pm$ 0.01

Equations

$$\mathcal{H}(\{{\mathbf{R}}_n,{\mathbf{P}}_n\})=$$

$$\frac{1}{2m}\sum _{n=1}^N|{\mathbf{P}}_n{|}^2+\sum _{n=1}^{N}V(|{\mathbf{R}}_n-{\mathbf{R}}_{n-1}|)-\mathbf{F}\cdot {\mathbf{R}}_N.$$

$${\mathbf{f}}_n=-V^{\prime }(r_n)\frac{{\mathbf{r}}_n}{r_n},$$

$$\frac{\text{d}{\mathbf{r}}_n}{\text{d}t}$$

$$\frac{\text{d}{\mathbf{P}}_n}{\text{d}t}$$

$${ϵ}_n=\frac{1}{2m}|{\mathbf{P}}_n{|}^2+V(r_n).$$

$$\frac{\text{d}{ϵ}_n}{\text{d}t}=\frac{1}{m}{\mathbf{f}}_n\cdot {\mathbf{P}}_{n-1}-\frac{1}{m}{\mathbf{f}}_{n+1}\cdot {\mathbf{P}}_n,$$

$$\mathcal{P}(\{{\mathbf{r}}_n,{\mathbf{P}}_n\})\text{d}\Gamma =$$

$$\frac{1}{Z}\prod _{n=1}^N\exp \left(-\frac{\beta }{2m}|{\mathbf{P}}_n{|}^2\right)\mathrm{d}{\mathbf{P}}_n\exp \left(-\beta V(r_n)+\beta \mathbf{F}\cdot {\mathbf{r}}_n\right)\mathrm{d}{\mathbf{r}}_n.$$

$$Z(\beta ,\mathbf{F})={\left(\frac{2\pi m}{\beta }\right)}^{3N/2}\zeta (\beta ,\mathbf{F}{)}^N,$$

$$S_k(\beta ,F)$$

$$C_k(\beta ,F)$$

$$\zeta$$

$$⟨x⟩$$

$$⟨ϵ⟩$$

$$⟨x^2⟩$$

$$⟨y^2⟩$$

$$⟨P_x^2⟩$$

$$⟨{ϵ}^2⟩$$

$$⟨x⟩=\ell (\beta ,F)\frac{F_x}{F},⟨y⟩=\ell (\beta ,F)\frac{F_y}{F},$$

$$⟨z⟩=\ell (\beta ,F)\frac{F_z}{F},$$

$$⟨P_x⟩=m{V}_{0x},⟨P_y⟩=m{V}_{0y},$$

$$⟨P_z⟩=m{V}_{0z},$$

$$⟨ϵ⟩=\frac{1}{2}m|{\mathbf{V}}_0{|}^2+e(\beta ,F).$$

$$g_1(n,t)=x_n(t),g_2(n,t)=y_n(t),g_3(n,t)=z_n(t),$$

$$g_4(n,t)=P_{x,n}(t),g_5(n,t)=P_{y,n}(t),g_6(n,t)=P_{z,n}(t),$$

$$g_7(n,t)={ϵ}_n(t),$$

$$J_1(n,t)=-\frac{1}{m}{P}_{x,n}(t),J_2(n,t)=-\frac{1}{m}{P}_{y,n}(t),$$

$$J_3(n,t)=-\frac{1}{m}{P}_{z,n}(t),J_4(n,t)=f_{x,n+1}(t),$$

$$J_5(n,t)=f_{y,n+1}(t),J_6(n,t)=f_{z,n+1}(t),$$

$$J_7(n,t)=\frac{1}{m}{\mathbf{f}}_{n+1}(t)\cdot {\mathbf{P}}_n(t).$$

$$\frac{\text{d}}{\text{d}t}{g}_{\alpha }(n,t)+J_{\alpha }(n,t)-J_{\alpha }(n-1,t)=0,$$

$$⟨h⟩(n,t)=\sum _{n^{\prime }=1}^N{h}_n\Theta (n^{\prime }-n),$$

$$\frac{\partial }{\partial t}{\mathcal{G}}_{\alpha }(n,t)+\frac{\partial }{\partial n}{\mathcal{J}}_{\alpha }(n,t)=0,$$

$${\mathcal{G}}_1$$

$${\mathcal{G}}_4$$

$${\mathcal{G}}_7$$

$$F_x$$

$$V_{0x}$$

$$e$$

$${\mathcal{J}}_1$$

Full text

In the framework of the non-linear fluctuating hydrodynamic theory (NLFHT), that was successfully applied to 1D chains, we obtain a heat mode, two longitudinal and four transverse sound modes. We compute their spatio-temporal correlations comparing the theoretical results with molecular dynamics simulations, finding a good agreement for high temperatures. We find that the transverse sound modes behave diffusively, meanwhile the heat and longitudinal sound modes superdiffusively, exploring their possible scaling functions and characteristic exponents.

I Introduction

Anomalous transport have been observed in a wide variety of systems revealing that this phenomenon is ubiquitous in nature. Its emergence can be found in processes taking place in plasmas, glassy materials, porous media, polymers and flexible filaments in viscoelastic media, biological cells, heat conduction, chemical reaction-diffusion, epidemic spreading, just to mention some libro1 ; libro2 ; cells .

Low dimensional systems, such as polymers, nanowires, and nanoribbons, generally present anomalous thermal transport properties, in other words Fourier’s Law is not fulfilled. This anomalous behaviour have been studied theoretically, numerically and also experimentally dhar2 ; li .

Among different theoretical approaches to study anomalous thermal transport, it has been recently proposed a nonlinear extension of fluctuating hydrodynamics. It is a solid scheme that can be applied to general dynamics as long as the interactions are local and translationally invariant, and there are locally conserved fields. Under these conditions it is possible to compute steady state average currents, that are functions of the conserved fields Spohn2014 .

Models for a chain of atoms with anharmonic interactions have been considered to study thermal transport under this approach. However, in these models the atoms were allowed to vibrate only in one dimension and periodic boundary conditions were usually considered EvenOlla2014 ; Spohn2014 ; sponh2 ; dhar ; Spohn2015 . In more realistic systems the interaction depends on the relative position between atoms, so vibrations in different directions are coupled. Moreover, experimental conditions usually require setups that are suspended or clamped by their ends and subject to external stresses.

In this context, we extend recent works on anomalous thermal transport in low dimensional systems nos1 ; nos2 in the framework of the NLFHT to the case of stressed anharmonic chains with 3D motion. The theoretical results are checked by molecular dynamics simulations.

In section II, we present the model and the equations of motion, calculating average values in the canonical ensemble. In Section III, we apply NLFHT to our model transforming the elongation, momentum and energy (fields) of the particles to a normal mode basis, characterizing their correlations. We also discuss the scaling functions in the framework of mode-coupling theory (MCT). In section IV we compute through molecular dynamics simulations the evolution and spatio-temporal correlations of the fields, comparing with the theoretical predictions. Finally, in Section V we present some brief conclusions.

II Model and global equilibrium

The model consists of a chain of $N+1$ identical particles of mass $m$, labeled with index $0\leq n\leq N$. The first particle $(n=0)$ is fixed in the coordinate origin. The other $N$ particles can move in the three $x$-$y$-$z$ directions and are subjected to a nearest neighbor potential $V(r)$, that depends only on the distance between them. At this point, we consider a general potential with the condition to have one minimum at a finite equilibrium distance $r_{0}$. Moreover, the last particle $n=N$ is subjected to an external constant force ${\bf F}$, which provides a tension along the chain.

In terms of positions ${\bf R}_{n}$ and momenta ${\bf P}_{n}$ of the particles, the Hamiltonian is

[TABLE]

Due to the nearest neighbor interactions, the equations of motion given by this Hamiltonian can be rewritten in a simpler way in terms of the inter-particle coordinates ${\bf r}_{n}={\bf R}_{n}-{\bf R}_{n-1}$, and the force that particle $(n-1)$ does on particle $n$

[TABLE]

with $r_{n}=|{\bf r}_{n}|$, finally giving

[TABLE]

These equations are valid for $1\leq n\leq N$, reminding that ${\bf P}_{0}=0$, and defining ${\bf f}_{N+1}=-{\bf F}$.

The energy per site is defined as

[TABLE]

This is the kinetic energy of particle $n$ plus the potential energy between particle $(n-1)$ and $n$. We stress that this is an arbitrary definition, because the potential energy of a bond is shared by two particles. The derivative with respect to time of this site energy is

[TABLE]

which corresponds to a local conservation of energy. In the limit $N\rightarrow\infty$, equilibrium statistics of the system tend to the same results both in the microcanonical or canonical ensembles. In the canonical ensemble, the probability density to find the system in a particular configuration is ${\cal P}(\{{\bf R}_{n},{\bf P}_{n}\})\textnormal{d}\Gamma=Z^{-1}\exp(-\beta{\cal H})\textnormal{d}\Gamma$, with $\beta=1/(k_{\mathrm{B}}T)$ the inverse temperature, $Z$ the canonical partition function, and $\textnormal{d}\Gamma=\prod_{n}\mathrm{d}{\bf R}_{n}\mathrm{d}{\bf P}_{n}$ the differential of phase space. The average value of any quantity is computed integrating it over the phase space with this probability density. Although the canonical variables are ${\bf R}_{n}$ and ${\bf P}_{n}$, we can again take advantage of the structure of the Hamiltonian, and change the space coordinates from ${\bf R}_{n}$ to the ${\bf r}_{n}$. This transformation has a Jacobian equals to one. Moreover, ${\bf R}_{N}=\sum_{n=1}^{N}{\bf r}_{n}$, and the probability density factorizes completely to

[TABLE]

We see explicitly that each component of momentum is a random Gaussian variable. The elongations ${\bf r}_{n}$ for different bonds are also decorrelated, although the three cartesian components of each one are correlated through the term $V(r_{n})$. Integrating this probability density we obtain the partition function

[TABLE]

with $\zeta(\beta,{\bf F})=\int\exp\left[-\beta(V(r)-{\bf F}\cdot{\bf r})\right]\mathrm{d}{\bf r}$. To compute this function and other mean values on the canonical ensemble, it is useful to change to spherical coordinates, putting the polar axis in the direction of the external force. In this case the integrals in the angular coordinates are easily performed. To simplify the explicit expressions it is also useful to define the following radial integrals

[TABLE]

These integrals converge only for potentials growing faster than linear, i.e. $V(r)/r\rightarrow\infty$ for $r\rightarrow\infty$. This condition corresponds to interparticle forces at large elongations being stronger than the applied external force. In terms of these radial integrals, the partition function, the mean values of the elongation, and the energy are expressed as

[TABLE]

The functions $\ell$ and $e$ are of relevance for the following theoretical calculations. To pass to the microcanonical ensemble, in the limit $N\rightarrow\infty$, one can invert the relations to obtain $\beta(\ell,e)$ and $F(\ell,e)$.

The mean values of momentum and force are $\langle P_{x}\rangle=\langle P_{y}\rangle=\langle P_{z}\rangle=0$, $\langle f_{x}\rangle=-F$, and $\langle f_{y}\rangle=\langle f_{z}\rangle=0$. It is also possible to compute the second moments, obtaining

[TABLE]

We can generalize this ensemble to a situation where there is a center of mass velocity ${\bf V}_{0}$ (it corresponds to a constant velocity of the first particle $n=0$ originally fixed), and an arbitrary direction of the force. With these changes, the mean values become

[TABLE]

where $\ell$ and $e$ are the same previously defined functions.

III Nonlinear fluctuating hydrodynamic theory

III.1 Mesoscopic fields and currents

Following the theoretical developments in Spohn2014 applied to anharmonic chains, we extend it to our model for a stressed chain vibrating in 3D. We first define seven microscopic fields

[TABLE]

and their corresponding microscopic currents

[TABLE]

By these definitions, the equations of motion (3) and (4) and conservation of energy (6) can be written in the compact form

[TABLE]

which are the Euler equations. The index $n$ is discrete, and the equations are valid for $n=1$ up to $n=(N-1)$, being the first particle ($n=0$) fixed, and the last particle ($n=N$) is subjected to the external force ${\bf F}$.

To connect these fields to the results of the generalized canonical ensemble, we define for any arbitrary microscopic quantity $h_{n}(t)$ a coarse-graining or mesoscopic average

[TABLE]

with $\Theta(x)$ a smoothing function with the following properties: positive $\Theta(x)\geq 0\ \forall x$, normalized $\int\Theta(x)\textnormal{d}x=1$, and with a finite variance $\sigma^{2}=\int x^{2}\Theta(x)\textnormal{d}x$ such that $1\ll\sigma\ll N$, typically a Gaussian function. Now $n$ becomes a continuous variable along the chain, instead of a discrete index. This mesoscopic average is valid far from the ends of the chain $(3\sigma\lesssim n\lesssim N-3\sigma)$.

This coarse-graining, in a statistical sense, is related to a local thermal equilibrium of the chain around $n$. The averaged fields and currents become smooth functions, not varying too much in the scale of few sites, finally providing the (continuous) Euler equations

[TABLE]

where ${\cal G}_{\alpha}=\langle g_{\alpha}\rangle$ and ${\cal J}_{\alpha}=\langle J_{\alpha}\rangle$ are the smoothed fields and currents, respectively. A local microcanonical equilibrium can be mapped to the general canonical ensemble defined by the seven parameters ${\bf V}_{0},{\bf F}$ and $\beta$ (local center of mass velocity, local tension and local temperature), which can smoothly vary over the chain, not far from the global equilibrium, giving for the smooth fields

[TABLE]

One can invert these relations to obtain the seven parameters of the ensemble in terms of the fields

[TABLE]

where now $\ell=\sqrt{{\cal G}_{1}^{2}+{\cal G}_{2}^{2}+{\cal G}_{3}^{2}}$ is also a function of the smooth fields. Temperature $\beta$ and force modulus $F$ are defined implicitly in terms of $\ell$ and $e$. This inversion allows to express the local mesoscopic currents in terms of the averaged fields:

[TABLE]

For the last current ${\cal J}_{7}$, we have used that the microscopic force, which only depends on spatial coordinates, is decorralated from the microscopic momentum, as in the general canonical ensemble.

We see explicitly that the currents are non-linear functions of the fields (except for the first three components). Considering that the local equilibrium is not far from a global equilibrium, the currents can be expanded up to second order.

The global equilibrium is a uniform state over the chain, where the fields are

[TABLE]

This state can be defined by a global temperature $\beta_{0}$ and external tension $F_{0}$ in the $x$ direction, therefore $\ell_{0}=\ell(\beta_{0},F_{0})$ and $e_{0}=e(\beta_{0},F_{0})$. On the other hand, we can constrain the last particle to move in a plane $y$-$z$ at distance $L_{0}$ from the first particle, with a fixed total energy of the chain $E_{0}$, then $\ell_{0}=L_{0}/N$ and $e_{0}=E_{0}/N$. In the thermodynamic limit $N\rightarrow\infty$ both conditions would give the same results.

Around this global equilibrium, we expand the currents up to second order

[TABLE]

with the Jacobian and Hessians matrices

[TABLE]

[TABLE]

and $u_{\alpha}(n,t)={\cal G}_{\alpha}(n,t)-{\cal G}_{0\alpha}$, the fluctuations of the fields around the global equilibrium, which are typically small.

The matrix $A$ can be computed explicitly from Eqs. (25), in terms of the derivatives $F_{\ell}=\partial F/\partial\ell$, and $F_{e}=\partial F/\partial e$, giving

[TABLE]

Although the force $F$ is an implicit function of $\ell$ and $e$, its derivatives can be computed in terms of derivatives of $\ell$ and $e$ with respect to $\beta$ and $F$ Spohn2014 , which in turn can be expressed in terms of the second moments of the fields. In the same way, the Hessian matrices $H^{\alpha}$ are written in terms of first and second derivatives of the external force, needing up to third moments of the fields. Their full explicit expressions are in the Appendix.

We are interested in the spatio-temporal correlations between the fluctuation of the fields

[TABLE]

This average is done over different reference sites $n_{0}$ far from the borders, and for different reference times $t_{0}$. For $n=0$ and $t=0$, this auto-correlation matrix can be computed explicitly from the second moments computed in the canonical ensemble:

[TABLE]

with

[TABLE]

where for simplicity $\delta x=x-\langle x\rangle$, and $\delta V=V-\langle V\rangle$.

With the previous expansion of the currents, the Euler equations up to second order in the fields are

[TABLE]

In these equations all fields are coupled. Nevertheless, up to linear order, the equations can be decoupled by diagonalizing the matrix $A$, which is guaranteed by the property $CA=CA^{\text{T}}$, being $C=C_{\alpha\beta}(0,0)$, which was checked explicitly in our model. This provides seven normal modes with eigenvalues

[TABLE]

where

[TABLE]

and

[TABLE]

The solutions of the Euler equations up to linear order are travelling waves with velocities $c_{i}$. The mode with eigenvalue zero is called the heat mode. There are two longitudinal sound modes (each travelling in opposite directions along the chain), and four transverse sound modes, with degenerate sound velocities $c_{\perp}$. Explicit expressions of the eigenvectors are in the Appendix. In this mode representation, the original $u_{\alpha}$ fields are transformed to

[TABLE]

where $R$ is the basis change matrix composed by the left eigenvectors of $A$ as rows. We order the modes in the following way:

In this new basis the spatio-temporal correlations are:

[TABLE]

where the average is computed in the same way as it was explained for the original $u_{\alpha}$ fields.

III.2 Noise, diffusion and spatio-temporal correlations

After considering the coarse-graining smoothing of the elongation, momentum and energy fields, and their corresponding currents, it is necessary to analyze the short range fluctuations to characterize the broadeaning of correlations in time.

On top of the smooth fluctuations along the chain, there are short range fluctuations that can be considered random. One can model these fluctuations adding diffusion and noise to Eq. (35)

[TABLE]

with $u_{\alpha}$ indicating the fields subjected to the stochastic fluctuations $\xi_{\beta}(n,t)$ which are random white noise components whose correlations are

[TABLE]

$B_{\alpha\beta}$ is the noise strength matrix, which is diagonal sponh2 . The noise does not affect the first three elongation components, and due to the symmetry of the model in the transverse directions the matrix $B$ has the following structure

[TABLE]

The symmetric diffusion matrix $D_{\alpha\beta}$ is related to the noise matrix and the correlations by the generalized fluctuation-dissipation theorem

[TABLE]

From this relation we obtain the diffusion matrix explicit and uniquely (see Appendix for its full expression).

In the mode basis, Eq. (41) transforms to:

[TABLE]

where $\tilde{D}=RDR^{-1}$ and $\tilde{B}=RB$. In the new basis the relation Eq. (44) transforms to $\tilde{D}+\tilde{D}^{\text{T}}=\tilde{B}\tilde{B}^{\text{T}}$. The matrices $G^{\alpha}$ are obtained form the $H^{\alpha}$ matrices as:

[TABLE]

As the modes travel with their own sound velocities, the crossed terms for the correlations ${\cal S}_{\alpha,\beta}(n,t)$ for $\alpha\neq\beta$ tends to zero for long times. In our model the only exceptions could be ${\cal S}_{4,6}$ and ${\cal S}_{5,7}$, because the sound velocities in the directions $y$ and $z$ are degenerated.

Given also the symmetry between modes $\phi_{2}$, $\phi_{4}$, and $\phi_{6}$ (travelling along the chain in positive direction) with modes $\phi_{3}$, $\phi_{5}$, and $\phi_{7}$ (travelling in the negative direction), the relevant correlations to study are ${\cal S}_{1,1}$, ${\cal S}_{2,2}$ and ${\cal S}_{4,4}$.

For long time it is expected the following general scalings for the correlations:

[TABLE]

The longitudinal sound mode has a self-coupling quadratic term $G^{2}_{22}$ that does not vanish. Therefore, the corresponding Eq. (41) has the structure of a noisy Burgers equation whose solution is the scaling Kardar-Parisi-Zhang (KPZ) function, with a characteristic exponent $\delta_{2}=2/3$.

For the heat and transverse sound modes, the self-coupling terms $G^{1}_{11}$ and $G^{4}_{44}$ vanish. Therefore, to compute their auto-correlations ${\cal S}_{1,1}$ and ${\cal S}_{4,4}$, respectively, it is necessary to take into account subleading corrections. This can be done within the mode-coupling theory. In this approach the dynamical evolution of the correlations is given by integro-differential equations involving the following memory kernel Spohn2014 :

[TABLE]

For long times, due to the spread of the correlations at different velocities, the leading contributions to this kernel come from the diagonal terms, plus eventual contributions from the degenerate transverse modes.

For the transverse sound mode, due to the symmetries of the $G^{4}$ matrix, the kernel in Eq. (50) vanishes. Under these conditions the equation for this mode reduces to a normal diffusion equation, with a Gaussian solution with a characteristic exponent $\delta_{4}=1/2$.

For the heat mode, due to the structure of the matrix $G^{1}$, the kernel has contributions from all the other six sound modes. If one considers only the coupling with the two longitudinal sound modes, whose correlation function approaches to a KPZ with exponent 2/3, one expects for the heat mode correlation a Levy function with exponent $\delta_{1}=3/5$. On the other hand, considering only a coupling with the four transverse sound modes, whose correlations are Gaussian functions, one would expect a Levy function with exponent $\delta_{1}=2/3$ for the heat mode correlation Spohn2015 .

Our model is more complex because it mixes both situations, so the exponent $\delta_{1}$ depends on the particular parameters, through the $G^{1}$ matrix elements. Nevertheless, the previous analysis indicates a superdiffusive behavior $\delta_{1}>1/2$.

IV Molecular Dynamics simulations

We apply the previous theoretical results to the case of a Fermi-Pasta-Ulam potential

[TABLE]

which can be considered as a general expansion of an arbitrary potential around its minimum. The quartic term $k_{4}$ should be bigger than zero, to guarantee the convergence of all integrals needed for the canonical averages. We consider $k_{2}=1$, $k_{3}=-5$ and $k_{4}=16$. These values guarantee the existence of only one minimum at $r=r_{0}$. On the other hand, $k_{3}$ tunes the asymmetry of the potential which is negative in our case, to model a repulsive force at short distances.

We take $r_{0}$ as unit of length and $k_{2}r_{0}^{2}$ as unit of energy and temperature.

Depending on temperature and on the applied tension the chain presents two different configurations. When the mean elongation $\ell$ is less than one, the chain is soft and the system has a huge number of available configurations, making difficult to explore numerically the full phase space. On the other hand, if $\ell>1$, the chain is more tight and the available phase space for the system is drastically reduced.

We consider four sets of parameters that are representative of regimes for low and high temperatures with and without strain:

A) $T=0.1$ and $\ell=1$ ,

B) $T=0.1$ and $\ell=1.1$ ,

C) $T=1$ and $\ell=1$ ,

D) $T=1$ and $\ell=1.1$ .

For each set we compute analytically the $R$ matrix in order to obtain the seven corresponding fields given by Eq. (39) as a function of time, from the elongation, momentum and energy.

To compute numerically the spatio-temporal correlations of these fields it is necessary to do a statistical average starting from a thermal equilibrated initial condition. For this purpose, we first integrate numerically the stochastic extension of Eqs. (3) and (4) at the desired temperature, using a Runge-Kutta algorithm. After the equilibration was achieved, these equations are integrated in the microcanonical ensemble using a Velocity-Verlet algorithm.

The last stage of the numerical integration allows to calculate the normal mode fields. Doing statistical average in time, we obtain the corresponding spatio-temporal correlations.

As the sound modes come in pairs travelling in opposite directions along the chain, and there is a degeneration in the transverse directions, we only show the correlations for the heat mode $\phi_{1}$, the longitudinal sound mode $\phi_{2}$ and the transverse sound mode $\phi_{4}$.

In Fig. 2 we plot the mode correlations ${\cal S}_{1,1}$, ${\cal S}_{2,2}$ and ${\cal S}_{4,4}$ for the set of parameters C) and D) as a function of the site and for different times.

We observe that in all cases the peaks broaden and flatten for increasing times. The sound modes propagates along the chain with constant velocities while the heat mode remains at its initial position.

Also, for the heat peak, we observe the presence of two small lateral peaks, which corresponds to the coupling with the transverse sound mode due to the nonlinear terms in the hydrodynamic expansion. It is expected that at longer times these modes are decoupled.

From the position of the peaks, we compute numerically the velocities of the sound modes (see Table 1), observing a good agreement with the theoretical predictions for all sets of parameters.

We obtain numerically the characteristic exponents $\delta_{1}$, $\delta_{2}$ and $\delta_{4}$ in Eqs. (49), from the broadening rate for each peak, which is the same as the flattening rate, within the numerical errors. The results are presented in Table 2.

We find for the heat mode and for the longitudinal sound mode that their exponents are in all cases larger that $1/2$ corresponding to a superdiffusive behavior. On the other hand for the transverse sound mode the exponent is close to $1/2$, corresponding to a normal diffusion.

For sets C and D, the numerical exponents agree with the theoretical predictions, while for sets A and B the exponents deviate, specially for the heat mode. These deviations at low temperatures are expected because the system explores regions around the minimum of the potential, where the non-linear terms are negligible.

On the other hand, at high temperatures the system is highly anharmonic, so NLFHT and MCT shoud work fine.

For the transverse sound mode, the normal diffusive behaviour gives a Gaussian scaling function, with an exponent $\delta_{4}=1/2$. For the longitudinal sound mode, as it was discussed previously, we expect a KPZ scaling function, with an exponent $\delta_{2}=2/3$, which is in near the numerical value.

For the heat mode, we obtain characteristic exponents $\delta_{1}$ in the expected range between $3/5$ and $2/3$, but closer to $3/5$.

We show in Fig. 3, for different times, the shifted and rescaled correlations peaks by their expected velocities and characteristic exponents. We also superimpose the expected corresponding scaling functions. It is observed a good agreement for the transverse and longitudinal sound modes peaks with a Gaussian and a KPZ function, respectively. For the heat mode we propose a Levy function with an exponent $5/3$, which is closer to our numerical estimations. The deviations observed at the tails may indicate that the long time limit was not yet achieved.

Numerically, we observe for high temperature, that the tension does not modify the characteristic exponents of the sound modes (within the statistical errors), nor their scaling functions (see Table 2). However, for the sound mode there is a slight dependence of its exponent with the tension.

V Conclusions

In this paper we studied an anharmonic chain vibrating in 3 dimensions, with first neighbors interatomic potential that depends on the distance, and subjected to an external tension. This is a more realistic model that can be applied to e.g. suspended nanowires and polymers.

The 3D motion and the tension gives to our model a complex behavior in comparison with previous 1D anharmonic chain models. We applied the NLFHT where the new features of our model result in the emergence of two longitudinal and four transverse sound modes, added to a single heat mode. For longitudinal sound modes, the self-coupling term does not vanish, implying for its spatio-temporal correlation a superdiffusive behavior with a characteristic exponent of 2/3.

On the other hand, for a transverse sound mode, the self-coupling term vanishes, and it was necessary to apply MCT. Even within this theory, this mode does not couple to any others, resulting in a diffusive behavior for its spatio-temporal correlation.

Finally, for the heat mode, the self-coupling term also vanishes. MCT shows that this mode couples to all the sound modes, infering a more complex scaling function with characteristic exponent in between 3/5 and 2/3.

As it is expected, we have checked numerically that this theory does not work at low temperatures, where the particles feels an almost harmonic potential.

For high temperatures, the numerical simulations are in a reasonable good agreement with the theoretical predictions for the sound velocities and the exponents. For the heat mode the exponent is closer to 3/5, but it seems to slightly depend on the tension, which affects the couplings with the sound modes. The role of the tension, as well as the boundary conditions, deserve further investigation.

These results for the correlations of the fields can be applied to obtain the thermal conductivity of the chain through the Green-Kubo formula, suggesting an anomalous thermal conduction.

Acknowledgments

The authors want to thank for the funding project PIO-CONICET 14420140100013CO.

Appendix

For a general derivative with respect to the microcanonical variables $\ell$ and $e$ of a quantity $h$, we have

[TABLE]

Additionally, if $h$ is an average on the canonical ensemble, we have the following rules for derivation with respect to $\beta$ and $F$

[TABLE]

These rules allow to compute the derivatives of the external force, needed for the $A$ matrix, which are obtained from the following expressions:

[TABLE]

Using that $F=F(\ell(\beta,F),e(\beta,F))$ and $\beta=\beta(\ell(\beta,F),e(\beta,F))$ it is possible to determine that

[TABLE]

with

[TABLE]

Replacing Eqs. (A5) we obtain

[TABLE]

We can also compute the derivatives of the correlators of order two

[TABLE]

in terms of correlators of order three. These results allow to compute the second derivatives of the external force.

For the $A$ matrix, the right eigenvectors are

[TABLE]

with $Z_{i}$ normalization constants. As columns, they form the matrix $R^{-1}$. The left eigenvectors of $A$ are

[TABLE]

which as rows, they form the matrix $R$. The constants $Z_{i}$ can be determined uniquely (up to a sign) imposing that the correlation matrix in the mode basis is the identity $(RCR^{\textnormal{T}}=I)$

[TABLE]

The Hessian matrices are explicitly

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The diffusion matrix has the following structure:

[TABLE]

with

[TABLE]