Thermal response of a Fermi-Pasta-Ulam chain with Andersen thermostats

TL;DR

This paper develops a fluctuation-response relation for the thermal response of Fermi-Pasta-Ulam chains with Andersen thermostats, enabling better predictions of thermal properties out of equilibrium.

Contribution

It introduces a novel fluctuation-response relation for inertial heat conducting systems driven out of equilibrium using Andersen thermostats.

Findings

The relation predicts thermal expansion and heat capacity more accurately than small perturbation methods.

Simulations confirm improved susceptibility estimates with the new relation.

The approach extends linear response theory to nonequilibrium inertial systems.

Abstract

The linear response to temperature variations is well characterised for equilibrium systems but a similar theory is not available, for example, for inertial heat conducting systems, whose paradigm is the Fermi-Pasta-Ulam (FPU) model driven by two different boundary temperatures. For models of inertial systems out of equilibrium, including relaxing systems, we show that Andersen thermostats are a natural tool for studying the thermal response. We derive a fluctuation-response relation that allows to predict thermal expansion coefficients or the heat capacitance in nonequilibrium regimes. Simulations of the FPU chain of oscillators suggest that estimates of susceptibilities obtained with our relation are better than those obtained via a small perturbation.

Equations

$${\dot{x}}_i$$

$$m_i{\dot{v}}_i$$

$${\rho }_{eq}(v)=\sqrt{\frac{m_i}{2\pi T_i}}{e}^{-\frac{m{v}^2}{2T_i}}$$

$$P(\Delta t)\propto e^{-\frac{\Delta t}{\tau }}$$

$$T_i\to T_i+\theta .$$

$${\chi }_{\mathcal{O}}(t)\equiv {\frac{\partial }{\partial \theta }⟨\mathcal{O}(t){⟩}_{\theta }|}_{\theta =0}$$

$$⟨\mathcal{O}(t){⟩}_{\theta }=\sum _{\omega }{\mathcal{P}}_{\theta }[\omega ]\mathcal{O}(t),$$

$$⟨\mathcal{O}(t){⟩}_{\theta }=\sum _{\omega }\mathcal{P}[\omega ]\frac{{\mathcal{P}}_{\theta }}{\mathcal{P}}[\omega ]\mathcal{O}(t)$$

$$⟨\mathcal{O}(t){⟩}_{\theta }=⟨\frac{{\mathcal{P}}_{\theta }}{\mathcal{P}}[\omega ]\mathcal{O}(t)⟩$$

$$\frac{{\mathcal{P}}_{\theta }}{\mathcal{P}}[\omega ]$$

$$=\prod _{\alpha }^{n_i}\sqrt{\frac{T_i}{T_i+\theta }}{e}^{\frac{m_i(v_i^{\alpha +}{)}^2}{2}\frac{\theta }{T_i(T_i+\theta )}}$$

$$=1+\sum _{\alpha }^{n_i}\frac{m_i(v_i^{\alpha +}{)}^2-T_i}{2T_i^2}\theta +\mathcal{O}({\theta }^2)$$

$${\chi }_{\mathcal{O}}(t)$$

$$\equiv \frac{1}{T_i}⟨\Gamma [\omega ]\mathcal{O}(t)⟩.$$

$$\Gamma [\omega ]=\frac{S[\omega ]-K[\omega ]}{2}$$

$$S[{\omega }^R]$$

$$K[{\omega }^R]$$

$$S[\omega ]$$

$$K[\omega ]$$

$${\chi }_{\mathcal{O}}(t)$$

$$\equiv {\mathcal{C}}^{-}(t)+{\mathcal{C}}^{+}(t)$$

$${\mathcal{C}}^{-}(t)$$

$${\mathcal{C}}^{+}(t)$$

$${\chi }_O(t)=\frac{1}{T}⟨S[\omega ]\mathcal{O}(t){⟩}_{\mathrm{eq}}$$

$$K[\omega ]$$

$$\text{with}\bar{\Gamma }[\omega ]$$

$${⟨\sum _{\alpha }[(v_i^{\alpha -}{)}^2-T_i]\mathcal{O}(t)⟩}_{\mathrm{eq}}$$

$${⟨\sum _{\alpha }[(v_i^{\alpha +}{)}^2-T_i]\mathcal{O}(0)⟩}_{\mathrm{eq}}$$

$${\chi }_{\mathcal{O}}(t)$$

$$U$$

$$U_{\mathrm{FPU}}(r)$$

$$r_i$$

$${\chi }_L^s\equiv \underset{t\to \infty }{\lim }{\chi }_L(t)$$

Full text

1 Introduction

The thermal response is usually associated with coefficients such as the specific heat or the thermal expansion coefficient. In equilibrium these coefficients may be computed by means of fluctuation-response formulas, in which the response to temperature variations is related to the natural fluctuations of the systems. For instance, the specific heat is the variance of the energy in equilibrium. As encoded in the Kubo formula tod92 , the response of equilibrium systems may be entirely described by its dissipation, or entropy production bai13 , in the transient toward a new equilibrium after the perturbation.

The thermal response of systems out of equilibrium is less understood. For example, the equilibrium theory cannot be used to fully determine how glasses undergoing a relaxation do respond to an increase of temperature bar99 ; may06 ; mag10 . Another instance of nonequilibrium regime is a system carrying a heat flow from a hot to a cold reservoir, say a cantilever heated by a laser at one end gei17 or larger oscillating devices used in gravitational wave detectors con13 . A change in the laser intensity, translated into a change of one boundary temperature, would lead to forms of thermal response that, again, cannot be described within the framework of equilibrium systems.

Studies on nonequilibrium linear response cug94 ; nak08 ; che08z ; vil09 ; spe09 ; sei10 ; pro09 ; ver11a ; lip07 ; alt16 ; bai09 ; bai09b ; bai10 were recently complemented by specific works on thermal response che09b ; bok11 ; bra15 ; pro16 ; bai14 ; yol16 ; fal16 ; fal16b ; bai16 . In particular, a linear response approach first developed for the case of mechanical perturbations bai09 ; bai09b ; bai10 was recently extended to include the response to temperature variations bai14 ; yol16 ; fal16 ; fal16b ; bai16 ; yol17 . Those works dealt with overdamped stochastic systems and provided fluctuation-response relations based on standard integrals, thus regularizing singular terms present in the earlier contribution bai14 . However, currently we do not have similar results for inertial systems. This means that the basic models of heat conduction, such as the classical driven Fermi-Pasta-Ulam (FPU) chain of coupled oscillators lep03 ; dha08 , cannot be described with a linear response theory for temperature variations. In fact, to our knowledge, despite the extensive amount of studies focusing on FPU models, there is no theory describing their thermal response to a change of one of the two boundary temperatures.

We derive a thermal fluctuation-response relation for (nonequilibrium) inertial systems driven by Andersen thermostats and80 ; frenkel . These thermostats are a standard tool for equilibrium molecular dynamics simulations in the canonical ensemble. We show that, as a key advantage, the Andersen thermostats do not bring all mathematical problems of Langevin heat baths, yet they bring to expressions displaying the physically relevant quantities, such as the entropy production in the reservoirs. The fluctuation-response relation we obtain is based on the stochastic nature of Andersen thermostats, which allows us to apply a well-established scheme for determining response functions and susceptibilities from ratios of path-weights. As a numerical example, we show two forms of susceptibility of the FPU system to a variation of one boundary temperature.

2 Thermal susceptibility of systems driven by Andersen thermostats

We consider systems evolving via Newton’s equations

[TABLE]

where $x_{i}$, $v_{i}$, $m_{i}$, and $F_{i}$ are, respectively, positions, velocities, masses and forces.

A thermalization of a degree of freedom $i$ to a temperature $T_{i}$ is achieved by a Andersen thermostat and80 ; frenkel . The velocity $v_{i}$ is updated at some instants with a value $v$ extracted from an equilibrium one-dimensional Maxwell-Boltzmann distribution

[TABLE]

where the Boltzmann constant is set $k_{B}=1$. As in Markovian jump processes, instants of updates are separated by time intervals $\Delta t$ extracted from an exponential distribution

[TABLE]

where $\tau$ is a parameter of the simulation.

We first study the response to the perturbation of a single temperature

[TABLE]

activated at time [math] and persistent till time $t$. A simultaneous perturbation of several temperatures would be simply be a linear combination of the following basic results, as shown at the end of this section. A trajectory $\omega=\{x_{i}(s),v_{i}(s)\}$, $s\in[0,t]$ has a path-weight denoted by ${\mathcal{P}}[\omega]$ and related averages are written as $\langle\ldots\rangle$. The starting point of a trajectory is from an arbitrary initial distribution $\rho(\{x_{i}(0),v_{i}(0)\})$ which is not recalled explicitly in the notation (in general it can also be non-steady state distribution). The notation for perturbed trajectories, which start from the same $\rho$ but evolve with modified $T_{i}+\theta$, is turned to ${\mathcal{P}}_{\theta}[\omega]$ and $\langle\ldots\rangle_{\theta}$.

For an observable ${\mathcal{O}}(t)$ by definition we have a thermal susceptibility

[TABLE]

While ${\mathcal{O}}$ could be a functional ${\mathcal{O}}[\omega]$ of the whole trajectory, we keep the notation ${\mathcal{O}}(t)$ to highlight the evaluation of such observable after a time $t$ from the formal activation of the perturbation $\theta$. Let us express the mean perturbed value of ${\mathcal{O}}$ by

[TABLE]

where $\sum_{\omega}$ is a simplified notation for indicating a formal inclusion in the statistics of all (uncountable) trajectories of the system. To write this susceptibility as a function of unperturbed correlation functions, we rewrite (5) as

[TABLE]

which is an average in the unperturbed system, namely

[TABLE]

To compute the path weights ratio, note that the velocity updates are stochastic only at instants when they are replaced by values extracted from Maxwell-Boltzmann distributions; the deterministic evolution in between these updates is the same in the perturbed and unperturbed dynamics. Thus the path-weight ratio ${\mathcal{P}}_{\theta}/{\mathcal{P}}$ can be written as a product over all velocity jumps of the velocity probabilities (perturbed and unperturbed) at every jump, which gives factors different from $1$ only for the degree of freedom $i$ driven by the perturbed reservoir. In a trajectory $\omega$ in a time span $t$ with $n_{i}$ updates of $v_{i}$, we use the index $\alpha$ for velocity updates and denote by $v_{i}^{\alpha+}$ the velocity right after the jump $\alpha$, that is the velocity extracted from the equilibrium distribution in the Andersen scheme. With this notation we have

[TABLE]

The expansion in the last line allow us to keep the first order in $\theta$ and write the susceptibility as the unperturbed correlation

[TABLE]

where $\Gamma[\omega]\equiv\frac{1}{2T_{i}}\sum_{\alpha=1}^{n_{i}}[m_{i}(v_{i}^{\alpha+})^{2}-T_{i}]$. We see that it is a correlation between the observable and a sum of kinetic-minus-bath temperatures. In order to gain more physical insight, we may split the functional $\Gamma[\omega]$ in two terms with well-defined time parity,

[TABLE]

with a time-antisymmetric component $S[\omega]$ and a time symmetric one, $K[\omega]$. This symmetries are understood by considering a time-reversed trajectory $\omega^{R}$, which corresponds to the evolution backward in time, with variables $\{x_{i},-v_{i}\}$, from the final state of $\omega$. By definition,

[TABLE]

We call $v_{i}^{\alpha-}$ the velocity of mass $i$ at the instant right before the update with index $\alpha$. In the time-reversed trajectory, that jump would thus be a transition from $-v_{i}^{\alpha+}$ to $-v_{i}^{\alpha-}$. It is easy to find that

[TABLE]

We see that $S[\omega]$ is a sum of kinetic energy gains (of the system) during interactions with the Andersen thermostat, all divided by temperature. This is the entropy production in the heat bath. The interpretation of $K[\omega]$ is not as straightforward. We may call $K[\omega]$ the dynamical activity lec05 ; mer05 ; gar09 ; ful13 or frenesy bai10 ; bas15 of the systems: in our model it contains all deviations between the bath temperature and the kinetic temperature “during” a jump (i.e., the average of the two values before and after the jump). It can be seen as a measure of how agitated the system is in comparison to what it should be on average according to $T_{i}$. With (12) and (13) one can reinterpret the susceptibility (9) as a sum of two correlations

[TABLE]

with

[TABLE]

In equilibrium at temperature $T_{i}=T$, for state observables ${\mathcal{O}}(t)={\mathcal{O}}(\{x_{i}(t),v_{i}(t)\})$ (thus excluding functionals ${\mathcal{O}}[\omega]$ such as integrated currents), we expect the combination $\frac{1}{2}(\mathcal{C^{-}}+\mathcal{C^{+}})$ to turn into $\mathcal{C^{-}}$ (the susceptibility should be the time integral of a response function from Kubo’s formula), namely

[TABLE]

which is the equilibrium correlation between observable and entropy produced into the environment (times the $1/T$ coming from the definition of thermal susceptibility). Since

[TABLE]

to prove that our approach recovers the standard structure in equilibrium, we just need to show that $\langle\bar{\Gamma}[\omega]O(t)\rangle_{\rm eq}=0$. In equilibrium we can exploit the time-translational invariance of correlation functions as well as their invariance for time reversal. Using time reversal on $\bar{\Gamma}[\omega]$ changes its jump variables to those $(v_{i}^{\alpha+})$ after a jump. From these considerations we get

[TABLE]

where the last term equals zero because velocities extracted by the Andersen procedure do not correlate with an observable at the beginning of the trajectory.

Before discussing some numerical results, we conclude this section by noting that the generalization to the perturbation of many degrees of freedom is trivial. If these have indices $i\in I$ and thermostats are independent stochastic processes,

[TABLE]

where now jumps have indices $\alpha_{i}$ that refer to their respective velocity index.

3 Numerical results for the FPU model

As a toy model for heat conduction, we consider a FPU chain of coupled oscillators ordered from $i=1$ to $i=N$, each one with with mass $m=1$. Forces are determined by the FPU interparticle quartic potential

[TABLE]

The presence of $\kappa_{3}\neq 0$ yields an asymmetric potential and $\kappa_{4}>0$ makes it well concave for $r\to\infty$. We selected $\kappa_{2}=27$, $\kappa_{3}=-63$ and $\kappa_{4}=48$, and $N=10$ oscillators. Thermostats are applied only at the two boundaries, with temperatures $T_{1}$ and $T_{N}$. These temperatures are in natural units, i.e. we continue using $k_{B}=1$. We use the same value of the typical time $\tau=1$ between velocity renewals in both Andersen thermostats. This time scale is of the same order of the oscillation times in the potential wells: indeed, we use temperatures in the range $T=0.1\div 0.5$, corresponding to thermal velocities $v_{T}\approx\sqrt{T}\approx 0.3\div 0.7$ and oscillation ranges $r_{o}\approx 0.4\div 1$ which are covered in timescales $\approx r_{o}/v_{T}$ of order $1$. Simulations are run with a much shorter time step $dt=0.005$ and data collection starts after a conservatively long relaxation that ensures having reached the stationary state. This is also guaranteed by the small size ($N=10$) of the FPU system. We thus should not be dealing with the issues of thermalization in one-dimensional long FPU systems lep03 ; dha08 ; ben11a .

In the following we show that the typical response of our FPU system takes place in time scales at least one order of magnitude longer than $\tau$. In fact, we observe a slow, scale-free convergence $\sim 1/t$ to the asymptotic $t\to\infty$ limit of the static response. The thermal fluctuations introduced randomly by thermostats at the boundary velocities naturally propagate in the chain due to the FPU interaction between the oscillators. A heat flux is set up within the chain if the two thermal reservoirs are at two different temperatures $T_{1}\neq T_{N}$.

We perturb the temperature $T_{N}$, hence the perturbed degree of freedom in the notation of the previous section is $i=N$. The two observables we consider are the chain length111Since $x_{i}$’s should be interpreted as displacements from equilibrium positions, the absolute length of the system is not specified and $L$ is more appropriately interpreted as the difference between the length and that at zero temperature. $L\equiv\sum_{i=1}^{N-1}r_{i}=x_{N}-x_{1}$ and the total potential energy $U$. In the first case, the susceptibility $\chi_{L}(t)$ quantifies the thermal response of the system size, which is always an expansion because $\kappa_{3}<0$. The susceptibility $\chi_{U}(t)$ of the energy instead generalizes the concept of heat capacitance.

To check the algorithm, we first look at equilibrium results. We set $T_{1}=T_{N}=0.1$ and, for the direct evaluation of susceptibilities, we use the temperature step $\theta=0.01$ for perturbations. We recover the correct equilibrium Boltzmann distributions (not shown) for kinetic and potential energies.

The susceptibilities computed without perturbing the system via equations (12)-(2) match those determined with a small perturbation, see Fig.1. Moreover, at equilibrium, one notes that $\mathcal{C^{+}}=\mathcal{C^{-}}$ as expected. Therefore one may reconstruct the susceptibility from only the usual correlation between the entropy production and the observable, i.e., the Kubo formula $\chi=2\mathcal{C^{-}}$.

Out of equilibrium, we do not find anymore $\mathcal{C^{+}}=\mathcal{C^{-}}$, rather the symmetric and antisymmetric terms may diverge dramatically, as shown for $T_{1}=0.1$ and $T_{N}=0.5$ in Fig. 2. Yet, with (12)-(2) we are able to correctly compute the susceptibility also in a strong nonequilibrium regime and without perturbing the system. One may also note that the susceptibility computed with the fluctuation-response formula, both in equilibrium and out of equilibrium, is smoother and affected by smaller errors than the susceptibility computed by perturbing the system.

Since the numerical results obtained with (12)-(2) are better than those collected via a real perturbation of the system, we use them for computing the susceptibility of the chain length to a variation of $T_{N}$ in different regimes. We then derive a static susceptibility

[TABLE]

The extrapolation to $t\to\infty$ is obtained by noting that data seem to converge toward the asymptotic value with a gap decaying $\sim 1/t$. A linear fit of the data as a function of $1/t$ is shown in Fig. 3. Having sampled at three different equilibrium conditions we note that $\chi_{L}^{s}$ is not constant in this model, rather it is increasing with $T=T_{1}=T_{N}$: we get $\chi_{L}^{s}=0.524(3)$ for $T=0.1$, $\chi_{L}^{s}=0.601(2)$ for $T=0.3$, and $\chi_{L}^{s}=0.625(2)$ for $T=0.5$.

The nonequilibrium cases suggest that $\chi_{L}^{s}$ is not very sensible to the temperature on the unperturbed side ($T_{1}$) if it is lower than $T_{N}$ (for $T_{1}=0.1$ and $T_{N}=0.5$ we get again $\chi_{L}^{s}=0.625(2)$). This seems not to be the case for the opposite condition $T_{1}=0.5$ and $T_{N}=0.1$, in which the perturbation of $T_{N}$ leads to a $\chi_{L}^{s}=0.553(3)$ higher than that in equilibrium at $T_{1}=T_{N}=0.1$. This is a genuine nonequilibrium effect related to the heat flux coming from the side of $T_{1}>T_{N}$. Finally, note that both nonequilibrium cases are different from that in equilibrium at the average temperature $T=0.3$.

4 Conclusions

A thermal response relation can be easily derived for inertial systems driven by Andersen thermostats and it turns out to contain quite simple expressions. This contrasts with the difficulty of deriving a fluctuation-response relation for temperature variations in systems following a Langevin dynamics. A sound version of the latter (namely a theory devoid of finite time-step singularities) is currently available only for overdamped systems fal16 ; fal16b ; yol17 . Moreover, a linear response theory for systems with deterministic degrees of freedom coexisting with stochastic ones, as in the FPU model driven at the boundaries, seems impossible for Langevin systems priv-com-GF : by reguralizing path integrals as previously done fal16 ; fal16b ; yol17 , part of the perturbation is shifted to deterministic degrees of freedom, for which one cannot expand the path probability as normally done for stochastic degrees of freedom. This suggests that Andersen thermostats are an ideal tool for studying thermal susceptibilities, in particular in boundary-driven systems.

The example we analyzed, a FPU chain driven by two boundary temperatures, shows how to define the alike of heat capacitance and of thermal expansion susceptibility in a simple system carrying a heat flux. Numerical results suggest that the response to a change of a temperature is better estimated with our fluctuation-response relation than by actually perturbing the system (which is already not convenient in simulations, as it would also require more CPU time). Relaxing systems could be studied in the same framework.

Acknowledgements

We thank Gianmaria Falasco for useful discussions.

Both authors wrote this paper and contributed with analytical calculations and with analyses of numerical data. FDA implemented the simulations.