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

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.


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 [1], the response of equilibrium systems may be entirely described by its dissipation, or entropy production [2], 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 [3,4]. 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 [5] or larger oscillating devices used in gravitational wave detectors [6]. 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 [7,8,9,10,11,12,13,14,15,16,17,18,19] were recently complemented by specific works on thermal response [20,21,22,23,24,25,26,27,28]. In particular, a linear response approach first developed for the case of mechanical perturbations [17,18,19] was recently extended to include the response to temperature variations [24,25,26,27,28,29]. 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 [24]. 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 [30,31], 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 [32,33]. 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 fluctuationresponse 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. The same extension seems not possible for systems putting in contact diffusing and deterministic degrees of freedom, namely the FPU system driven by Langevin thermostats [34].
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 equationṡ 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 = 1 to a temperature T i is achieved by a Andersen thermostat [32,33]. The velocity v i is updated at some instants with a value v extracted from an equilibrium one-dimensional Maxwell-Boltzmann distribution where the Boltzmann constant is set k B = 1. As in Markovian jump processes, instants of updates are separated by time intervals ∆t extracted from an exponential distribution where τ is a parameter of the simulation. We first study the response to the perturbation of a single temperature activated at time 0 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 has a path-weight denoted by P[ω] and related averages are written as . . . . The starting point of a trajectory is from an arbitrary initial distribution ρ({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 ρ but evolve with modified For an observable O(t) by definition we have a thermal susceptibility While O could be a functional O[ω] of the whole trajectory, we keep the notation O(t) to highlight the evaluation of such observable after a time t from the formal activation of the perturbation θ. Let us express the mean perturbed value of O by where ω 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 which is an average in the unperturbed system, namely 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 P θ /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 ω in a time span t with n i updates of v i , we use the index α for velocity updates and denote by v α+ i the velocity right after the jump α, that is the velocity extracted from the equilibrium distribution in the Andersen scheme. With this notation we have The expansion in the last line allow us to keep the first order in θ and write the susceptibility as the unperturbed correlation where We see that it is a correlation between the observable and a sum of kineticminus-bath temperatures. In order to gain more physical insight, we may split the functional Γ [ω] in two terms with well-defined time parity, We call v α− i the velocity of the oscillator the instant right before the update with index α. In the time-reversed trajectory, that jump would thus be a transition from It is easy to find that We see that S[ω] 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[ω] is not as straightforward. We may call K[ω] the dynamical activity [35,36,37,38] or frenesy [19,39] 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 with In equilibrium at temperature T i = T , for state observables O(t) (i.e., here we do not consider functionals O[ω]), we expect the combination 1 2 (C − + C + ) to turn into C − (the susceptibility should be the time integral of a response function from Kubo's formula), namely 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 to prove that our approach recovers the standard structure in equilibrium, we just need to show that Γ [ω]O(t) 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Γ [ω] changes its jump variables to those (v α+ i ) after a jump. From these considerations we get 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 ∈ I and thermostats are independent stochastic processes, where now jumps have indices α i that refer to their respective velocity index.

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 The presence of κ 3 = 0 yields an asymmetric potential and κ 4 > 0 makes it well concave for r → ∞. We selected   In the following 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 length 1 L ≡ N −1 i=1 r i = x N − x 1 and the total potential energy U . In the first case, the susceptibility χ L (t) quantifies the thermal response of the system size, which is always an expansion because κ 3 < 0. The susceptibility χ 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 θ = 0.01 for perturbations. We recover the correct equi-  The susceptibilities computed without perturbing the system via equations (12)- (14) match those determined with a small perturbation, see Fig.1. Moreover, at equilibrium, one notes that C + = 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 χ = 2C − .
Out of equilibrium, we do not find anymore C + = 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)- (14) 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)- (14) 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 The extrapolation to t → ∞ is obtained by a linear fit of the data as a function of 1/t, see  (2) for T = 0.5. The nonequilibrium cases suggest that χ s L 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 χ s L = 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 χ s L = 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.

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 [26,27,29]. 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 [34]. 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.