{'fraction_non_alphanumeric': 0.048821289427923666, 'fraction_numerical': 0.012678567497963946, 'mean_word_length': 4.008488590012127, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heat fluctuations of Brownian oscillators in nonstationary processes: fluctuation theorem and condensation transition

We study analytically the probability distribution of the heat released by an ensemble of harmonic oscillators to the thermal bath, in the nonequilibrium relaxation process following a temperature quench. We focus on the asymmetry properties of the heat distribution in the nonstationary dynamics, in order to study the forms taken by the Fluctuation Theorem as the number of degrees of freedom is varied. After analysing in great detail the cases of one and two oscillators, we consider the limit of a large number of oscillators, where the behavior of fluctuations is enriched by a condensation transition with a nontrivial phase diagram, characterized by reentrant behavior. Numerical simulations confirm our analytical findings. We also discuss and highlight how concepts borrowed from the study of fluctuations in equilibrium under symmetry breaking conditions [Gaspard, J. Stat. Mech. P08021 (2012)] turn out to be quite useful in understanding the deviations from the standard Fluctuation Theorem.


I. INTRODUCTION
At the level of the thermodynamical description of out of equilibrium transformations, the irreversibility of macroscopic processes is encoded in the second law expressed through strict inequalities satisfied by various thermodynamic observables, like entropy and free energies. At the more refined level of the statistical mechanical description of the same processes, the irreversibility ought to manifest as an asymmetry in the probability distributions of the fluctuations of these same quantities. The central and yet unsolved problem is to find general methods, comparable to those available in equilibrium, to characterize these probability distributions. A major advance in this direction has been made in the last two decades with the development of Fluctuations Theorems (FT), which allow to constrain the form of the distributions with a certain degree of universality and in conditions arbitrarily far from equilibrium [1]. This theoretical approach turned out to be very useful, both in experimental and numerical studies, in the characterization of several nonequilibrium systems, such as colloidal particles in harmonic traps [2][3][4][5], vibrated granular media [6][7][8], models of coupled Langevin equations [9][10][11], driven stochastic Lorentz gases [12,13], and active matter [14], just to name a few examples.
An important development in the field has been achieved with the understanding of the impact of symmetry and symmetry breaking on the FT of currents in stationary diffusive systems [15]. More recently, it has * Electronic address: andrea.crisanti@uniroma1.it † Electronic address: alessandro.sarracino@roma1.infn.it ‡ Electronic address: mrc.zannetti@gmail.com been shown that relations in the form of FT are of general occurrence also in the probability distributions of order parameters in equilibrium statistical mechanics, in connection with symmetry breaking [16]. Here we show that these equilibrium results do feedback in the understanding of the nonequilibrium FT, shedding new light on the conditions under which the FT in the standard form (namely, with a linear asymmetry function) is expected to hold [17][18][19][20]. In particular, we study the fluctuations of the heat exchanged with the environment by one or more Brownian oscillators in an interval of time [t w , t], during the relaxation following the instantaneous quench from high to low temperature. This protocol induces a nonstationary dynamics that is noninvariant under time translation and therefore the heat exchanged along a given trajectory explicitly depends on the two times (t w , t). The study of the FT in similar processes has been addressed in [21][22][23]. Here, we derive exact expressions for the heat probability density functions, for one, two and a large number of independent oscillators. Our analysis shows how, in the case of more than one oscillator, the heat probability obeys an FT which deviates from the standard form and whose physical meaning can be rationalized by resorting to the analogy with the equilibrium problem in the presence of a nonuniform external field.
In the case of a large number of degrees of freedom, arising, for instance, in the normal mode decomposition of an extended system, we present a computation based on the steepest descent method, which allows us to obtain an accurate description of the large deviation function of the exchanged heat. Here, we will use large deviation theory with a large number of degrees of freedom (and not for long time intervals) [24,25]. This kind of problem was considered previously in Refs. [26][27][28], showing that fluctuations undergo a condensation transition. Briefly, by this is meant that fluctuations in a multi-component system do condense if there exists a critical threshold above which the fluctuation is feeded by just one of the components (or degrees of freedom). Here, we analyse in detail the crossover region between the normal and the condensed phase, providing an explicit expression for the large deviation function. We also reconstruct numerically the phase diagram at finite temperature, in the parameter space (t w , τ = t − t w ), which shows a nontrivial reentrant behavior.
The paper is organized as follows. In Section II we present the general framework within which various FT forms are derived and we argue on physical grounds for the particular form we adopt. We also recall in some detail how an FT arises in statics as a consequence of symmetry breaking, highlighting its relevance for the dynamical problem. Section III is devoted to the simplest case of a single oscillator. We compute exactly the probability distribution of heat fluctuations, we introduce the basic concept of time-dependent effective temperature and we derive the FT in the "Gallavotti-Cohen" form [29]. In Section IV we consider the case of two oscillators, and we analyse the modifications arising in the FT form due to the presence of more than one degree of freedom. In Section V we consider the case of a large number of oscillators. We discuss the conditions leading to condensation and we map out numerically the phase diagram in the parameter space. The consequences on the FT are analysed. Finally, conclusions are drawn in Section VI.

II. GENERAL SETUP
There exist many variants of FT whose derivation [1,30,31] can be unified into a single master theorem. Let us first outline the general setting which applies to all fluctuation problems, in or out of equilibrium. Consider a sample space Ω with elements σ. Let µ(σ) and µ ′ (σ) be two arbitrary probability measures over Ω and let W(σ) be defined by where the * operation denotes an involutory transformation of Ω onto itself, i.e. (σ * ) * = σ. In the following this will be taken as the representation of the inversion element of the Z 2 = (I, * ) group, where I is the identity operator. Introducing an arbitrary function F (σ), from the above relation there follows the identity where · and · ′ denote expectations with respect to µ(σ) and µ ′ (σ), respectively. Taking for F (σ) the characteristic function θ M (σ|M ) of a certain observable M(σ), that is from Eq. (2) follows where P ′ and P are the probabilities of the events in the arguments induced by µ ′ (σ) and µ(σ), respectively, and in the right hand side there appears the expectation with respect to µ(σ), conditioned to M(σ) = M . Defining the * transformation on the set of random variables by M(σ) → M * (σ) = M(σ * ), which is also an involution, the above relation can be recast in the form with It appears, then, that Eq. (5) is the transposition to the level of the observable M of the underlying basic relation (1). In particular, regarding W as the bias which is necessary to apply on µ in order to construct µ ′ and K as the analogous quantity relating P and P ′ , Eq. (6) tells how these biases are related one to the other. As long as µ(σ) and µ ′ (σ) are arbitrary, the above result does not have predictive power. It becomes the master FT when µ ′ (σ) is taken with a definite relation to µ(σ), which constrains the form of P (M ) if the right hand side of Eq. (5) is accessible without having to actually compute the expectation, possibly via symmetry arguments. In the simplest case µ ′ (σ) is taken to be the same as µ(σ). Then, from Eq. (1) one has showing that W is odd W(σ * ) = −W(σ) and characterizes the asymmetry of the probability measure under inversion. Now, the question is to what extent this asymmetry is preserved or distorted as the description is moved up from the microscopic level to the higher one of random variables. The answer from Eq. (5) is given by Taking M to have definite parity, if it is even no information about the symmetry of P is obtained and Eq. (8) yields the conditional integral FT Conversely, if M is odd, one has whose usefulness depends on the possibility of assessing the form of K(M ), which is called the asymmetry function (AF). In particular, if M ∝ W, there follows an FT of the Gallavotti-Cohen [29] form with a linear AF If, instead, W and M are not simply related, the meaning of K(M ) in general is not immediately transparent.
For an overview of the variety of the different FT forms arising in the general case see Ref. [1].

A. FT and symmetry breaking in equilibrium
The latter remarks are well clarified in the equilibrium context, used in Ref. [16] to investigate the relation between FT and symmetry breaking. Let Ω and σ be the system's phase space and configurations, respectively. For definiteness, let σ = (s 1 , ..., s N ), with s i = ±1, be a spin configuration of a magnetic system on the lattice in the presence of an external site dependent field B = {B i }, whose equilibrium state is described by the probability measure where µ 0 (σ) is symmetric under spin inversion σ * = (−s 1 , ..., −s N ), while the exponential term breaks explicitly the Z 2 symmetry. We are interested in the fluctuations of the global magnetization M(σ) = i s i . From Eq. (12) follows and which takes a simple form only if the external field is uniform B i = B, ∀i, yielding an FT with the linear AF [16] K(M ) = 2βBM.
Instead, if B is not uniform, deviations from the FT arise.
In order to take a closer look, let the number of spins to become large and consider the case of the ideal paramagnet, in which µ 0 (σ) in Eq. (12) is the uniform measure Z −1 = [ i 2 cosh(βB i )] −1 . Then, by a straightforward saddle point computation one obtains the large deviation principle where m = M/N is the magnetization per spin and the large deviation function is given by Here, is the Helmotz free energy density, which depends on the field configuration B, we have defined B * (m) = {B * i (m)} with and x * (m) is obtained by solving with respect to x the equation of state Consequently, β −1 x * (m) is the shift to be applied to the external field on each site in order to produce m as the average magnetization per spin. From the definition (10) follows Averaging Eq. (19) over i, we can write with and Eq. (21) can be put in the form where we have defined In the particular case of the uniform external field and is the Legendre transform of f (B), which is even under m reversal, since B * (m) is odd. Hence, in this case Eq. (24) reproduces the result (15). In the nonuniform case g(m) is not the Legendre transform of f (B) and in general does not have a definite parity, leading to a nonlinear AF function. We will come back on this point in Section IV. What the above exercise shows is that the FT, in the sense of a linear AF, holds as long as the macrovariable M, whose fluctuations are considered, is conjugate to the symmetry breaking field. Instead, if M is not a conjugate variable, as it is the case with a site dependent B, the FT in the form (15) does not hold.

B. FT out of equilibrium
Let us, next, consider the nonequilibrium context. Assuming stochastic evolution, take for Ω the space of stochastic trajectories and for σ an individual trajectory. Then, σ * stands for the time reversed trajectory, while µ(σ) and µ ′ (σ) are the probability measures associated with two different evolutions, whose relation is specified from case to case. Here, as in the equilibrium problem, we shall be concerned with µ ′ (σ) = µ(σ), which in the dynamical context arises when there are no time dependent external parameters and the system evolves in contact with a single thermal reservoir at the final temperature T . Then, only heat is exchanged with the environment and one has [1] where P 0 is the initial probability distribution, x 0 and x t are the initial and final entries in the trajectory σ, β = 1/T and Q(σ) is the heat exchanged along the trajectory, which we take as negative if released to the environment.
Using the above form of W and the definition (6), the heat AF is given by whose understanding requires some clue on the role of the boundary terms. This can be gained from the work of Puglisi et al. [17]. In this paper, as anticipated in the Introduction, we shall be interested in the heat exchanged by a system of Brownian oscillators with the environment in an interval of time [t w , t]. Eventually, this will lead to recast the above equation in the form where k are single oscillators labels and Q k the heat exchanged by each one of them. The above expression is clearly analogous to Eq. (14), where the role of the nonuniform external field B is played by the set of affinities {∆β k }, which are k-dependent differences of inverse temperatures. The correspondence between the two problems helps to understand the deviations or modifications of the FT in terms of a collection of nonuniform degrees of freedom. It should be emphasized that the choice of taking µ ′ = µ, and therefore P ′ = P , is dictated by the particular physical setting of interest, since we consider the relaxation following a temperature quench and we want to compare the probability of exchanging the heat Q with that of exchanging −Q, in the same quench process, that is without time reversal. We focus on the asymmetry of the heat distribution in the given process, as it was done, for instance, in the experimental work of Gomez-Solano et al. [32].

III. BROWNIAN OSCILLATOR
The equation of motion for the single overdamped Brownian oscillator is of the Langevin typė where ω is the frequency and η is the white noise, modeling the interaction with the thermal bath at the temperature T , with expectations The Boltzmann constant will be taken k B = 1 throughout. Initially the system is in equilibrium at the temperature T 0 , with the position probability distribution where H(x) = 1 2 ωx 2 is the energy of the oscillator. Instantaneous cooling (quenches) or heating processes are realized by putting, at the time t = 0, the system in contact with the thermal bath at the temperature T < T 0 or T > T 0 , respectively. In the following, we shall be mainly interested in the case of the temperature quench.

A. Fluctuations of exchanged heat
Let us focus on the fluctuations of the heat exchanged by the oscillator with the thermal bath in the time interval (t w ≥ 0, t > t w ) after the temperature step. Since no work can be carried out on or by the system, due to ω constant, the heat exchanged in a single realization of the dynamical evolution coincides with the energy difference which is positive if heat is absorbed from the bath and negative if it is released to the bath. Then, the probability of exchanging the amount Q of heat is given by where P (x, t; x w , t w ) is the joint probability of the two events (x, t) and (x w , t w ), given by where G(τ ) = e −ωτ is the response function dependent on the time difference is the position variance at the time t w , whose initial value is given by ν 0 = T 0 /ω. The position probability distribution at the generic time t, obtained by integrating the above quantity over x w , maintains the equilibrium form (34) where β eff (t) is the time dependent inverse effective temperature defined by the equipartition-like statement [32] which yields For this quantity, which will play an important role in the following, we shall use the short hand notation and similarly for the inverse temperature β w = T −1 w or β t = T −1 t . Returning to Eq. (36) and introducing the integral representation of the δ function we obtain Carrying out the x and x w integrations and rotating the λ integration from the imaginary to the real axis, this can be rewritten as where with after using the notation which will be adopted from now on. Let us briefly comment on the parameters a and b introduced above. The asymmetry between λ + and λ − is controlled by b, while a controls the symmetric part. The sign of b depends on that of ∆T , implying λ + > λ − in the quench process and λ + < λ − in the heating process. The most asymmetrical situation is obtained for a = 0, which is realized either setting T = 0, or T 0 = 0 and t w = 0, yielding Instead, the symmetrical situation with is obtained when the system is equilibrated (b = 0) for ∆T = 0 or for t w → ∞.
The integral in Eq. (44) is evaluated by closing the contour either on the upper half plane or on the lower half plane around the branch cuts running along the imaginary axis from ±iλ ∓ up to ±i∞ (see Fig. 1), depending on Q < 0 or Q > 0, and obtaining where while K 0 is the modified Bessel function of the second kind.
In the limit t w → ∞ one has λ + = λ − = β, recovering the equilibrium result which was derived in Refs. [4,33] for a Brownian particle optically trapped in a stationary harmonic potential. The complementary case of a Brownian oscillator in the strongly underdamped limit has been studied in the recent work of Salazar and Lira [34], where a similar result for the heat distribution is derived. The dependence on t w of the asymmetry of the distribution, for a quench to a small but finite temperature T = 0.1, is illustrated in Fig. 2, where the analytical expression (51) is compared with the numerical simulation of the process (31).

B. Decomposition into cooling and heating contributions
The two quantities λ ± turn out to be the basic building blocks in all what follows. The physical meaning can be readily understood by rewriting P (Q) as the convolution of the two distributions arising in the purely cooling and in the purely heating process. Defining and carrying out the integral, which involves only one of the branch points, one finds and from which follows where Q ± are the average heats exchanged in the processes in which heat can only be absorbed or released. Then, the probability (44) can be rewritten as which reproduces the result of Eq. (51) and the total average heat exchanged is clearly given by The physical conditions for P (Q) = P − (Q) or P (Q) = P + (Q) are obtained putting the oscillator in contact with the thermal reservoir at T = 0, or by arranging a purely heating process with T 0 = 0, T > 0 and t w = 0. In the first case, one has T w = T 0 G 2 w , a = 0, b = (T w − T t ) from which follows, according to Eqs. (48) and (49), and, therefore, according to Eq. (55) (60) Similarly, in the second case, with T 0 = 0 and t w = 0, one (61)

C. FT and time reversal symmetry breaking
In this subsection we discuss the symmetry properties of the heat distribution. From Eq. (51) follows immediately the FT in the standard form A similar result was derived in Ref. [34] with t w = 0 and ∆β = β 0 − β, while in Ref. [32] there appears ∆β = β w −β t , which holds true only in the limit of large τ when the system is time decorrelated and β t ∼ β. Hence, as anticipated in section II, the AF is linear since the role of the initial probability P 0 now is played by P (x, t w ), which after Eqs. (35) and (38) yields Here, we want to highlight the connection between the above result and the breaking of the time reversal symmetry, by adapting to the present context the approach of Gaspard [16] sketched in Section II. Taking for Ω the set of ordered pairs σ = (x w , x), the time reversal symmetry operation is represented by the involution σ * = (x, x w ) which exchanges the order of x w and x. The task is to rewrite the joint probability Eq. (37) in the form of Eq. (12), that is as the product of a Z 2 symmetric factor µ 0 (σ) times an exponential where there appears the heat Q in the role of the observable explicitly breaking the Z 2 symmetry. In order to do this, let us first cast the joint probability Eq. (37) in the form where the "action" A(σ) can be read out from Eq. (37) as and the normalization factor is given by Z = 2πββwω 2 (1−G 2 ) . Decomposing A(σ) into the sum of an even and an odd part with and where Q(σ) is defined by Eq. (35), the form (12) of µ(σ) is obtained where is Z 2 invariant. Therefore, the invariance under time reversal in µ(σ) is explicitly broken by the exponential factor involving Q(σ) in Eq. (70). Once this is established, the FT in the form of Eq. (62) is recovered straightforwardly. The following remarks are in order: 1. the affinity ∆β driving the heat flow plays the role of the external field causing the explicit breaking of the Z 2 invariance.
and, comparing with Eq. (51), one obtains 3. From Eq. (44) it is straightforward to derive the moment generating function implying that the cumulant generating function ψ(ϑ) = ln e ϑQ satisfies the same symmetry which was found for a particle coupled to two thermostats at different temperatures [9].

IV. TWO BROWNIAN OSCILLATORS
Let us, next, see how the overall heat fluctuations are modified when there are internal degrees of freedom. We consider first the case of two uncoupled oscillators with frequencies ω 1 and ω 2 and, in the next section, we consider the limit of a large number of oscillators.
As before, initially the system is in equilibrium at the temperature T 0 and at the time t = 0 is put in contact with the reservoir at the temperature T . Denoting by Q 1 and Q 2 the amounts of heat exchanged by each oscillator, the probability that the system as a whole exchanges the heat Q is given by where P i (Q i ), with i = 1, 2, are the probabilities pertaining to each component. Inserting the expression (51) and changing the integration variable from Q 1 to y = Q 1 − Q/2, this can be put in the form analogous to Eq. (51) where and with Here, λ +,i are defined in Eq. (45), the index i = 1, 2 denoting the oscillator frequencies, ω 1 and ω 2 , respectively and ∆β i is defined by Eqs. (39) and (52) with frequency ω i . The resulting AF is given by Hence, if the two oscillators are equal (ω 1 = ω 2 ) the second term disappears, since R(Q) becomes an even function, and the linear result (63) is recovered. If they are unequal, R(Q) is not even under the change of sign of Q and the second term in the above equation gives a nonlinear contribution. The modifications introduced by this second contribution are illustrated by Fig. 3 obtained with parameters ω 1 = 0.5 and ω 2 = 0.2, T 0 = 10, T = 1 at times t = 12, t w = 10. Around Q = 0 one observes a linear behavior with slope ∆β (see dashed line in the inset), while for large values of Q the AF shows a nonmonotonic behavior. In order to expose the connection between the presence of internal structure and the equilibrium problem mentioned in the Introduction, let us go back to the microscopic measure. With two oscillators the phase space enlarges to the set of pairs Ω = {σ 1 , σ 2 } and the corresponding action, due to the absence of coupling, is additive The two terms in the right hand side are given by Eq. (66), keeping into account that the frequencies ω 1 and ω 2 are different, that is Carrying out the decomposition into even and odd parts, as in Eq. (67), the measure can be written as where is even under * : is odd. Hence, W(σ 1 , σ 2 ) = −2O(σ 1 , σ 2 ) and from Eq. (10) there follows with The formal similarity with Eq. (14) is evident, with the affinities ∆β i playing the role of the site dependent external field {B i } and the partial heat Q i that of the local order parameter s i . Accordingly, the above expression simplifies in the particular case akin to that of the uniform external field, that is for ω 1 = ω 2 or for t w = 0, yielding ∆β i = β 0 − β for all i.

V. EXTENDED SYSTEM
In this section we consider the case of a large number of oscillators. The aim is to analyse important qualitatively new features arising in the behavior of fluctuations when the size of the system becomes large. As it is usually the case for large systems, fluctuations of extensive quantities, such as heat, obey a large deviation principle [24,25]. The non trivial feature, unexpected for a linear system, is that in the large deviation function there appears a singularity corresponding to a condensation transition [26][27][28]. This means, as briefly anticipated in the Introduction, that there exists a critical threshold Q c of the exchanged heat, such that fluctuations above threshold are due to a single oscillator. We give a careful treatment of the saddle point computation involved in the computation of the large deviation function and we study the phase diagram of the transition. Mathematical details are collected in Appendix A.
We consider an extended system of linear size L and volume V = L d , where d is the space dimensionality. We suppose that this system is linear and that the normal modes decomposition with periodic boundary conditions gives rise to a set of harmonic components for the low lying modes with frequencies obeying the dispersion relation [35] ω(k) = k α + ω 0 , where ω 0 > 0, k = 2π L n, n i = 0, ±1, ±2, ... and α > 0. Just to fix ideas, a dispersion relation of this type arises in the Gaussian model or Van Hove theory of critical phenomena [36]. The k-space volume per mode is (2π) d /V . Assuming that there exists an underlying lattice structure with lattice spacing a 0 , the first Brillouin zone is bounded by Λ = 2π/a 0 and the total number of modes is given by N = V /a d 0 . Due to modes independence, the probability of a generic microscopic state at the time t w is given by the product measure where, according to Eq. (40), is the inverse effective temperature of the mode k and H k (x k ) = (1/2)ω(k)x 2 k is the single mode energy. Hence, Eq. (64) generalizes to is the heat exchanged by the mode k, and the sum is restricted to the first Brillouin zone. Inserting into Eq. (29) we recover Eq. (30) Comparing with Eq. (14), we see that the differences ∆β k = β k,w − β play the role of the site dependent external field in the equilibrium problem, as remarked in section IV, and that the AF linear form is recovered when these differences become independent of k, as for instance for t w = 0, or when the oscillators are identical, i.e. α = 0 (see Appendix B). Let us, next, see what form takes the AF by letting N to become large. Generalizing Eq. (76) and keeping V finite, the probability that the system exchanges the amount of heat Q is given by where λ ± (k) are defined in Eq. (45), evaluated for ω(k) as in Eq. (90). Denoting by λ + 0 and λ − 0 the upper and lower edges of the two branches of the spectrum, that is the contour of integration can be deformed to an arbitrary line Γ on the complex plane, provided it crosses the real axis within the gap [−λ − 0 , λ + 0 ]. This allows to rewrite the above integral as where q = Q/V is the density of the exchanged heat and The large V behavior of P (Q, V ) can be obtained by using the steepest descent method to estimate the integral in Eq. (97) [37].
with dµ(k) = Υ d k d−1 dk and Υ d = [2 d−1 π d/2 Γ(d/2)] −1 coming from the angular integration, is an uniform approximation to the function L(z; V ) as V ≫ 1. Here, Γ(x) denotes the Euler gamma function. Denoting by φ(z; q) = zq + L(z) and by primes derivatives with respect to z, it is straightforward to verify that Im φ ′ (z; q) ∼ Im z, implying that the stationary point z * of φ(z; q) is on the real axis. Thus, setting z = x + iy, the coordinate x * of the saddle point is obtained solving the equation φ ′ (x * ; q) = 0, which explicitly reads The function −L ′ (x) is monotonically increasing in the interval x ∈ [−λ − 0 , λ + 0 ], so Eq. (100) admits solution on the condition that q lies between the limiting values: with q − c < 0 < q + c . In this case the steepest descent path in the neighborhood of the saddle point is parallel to the imaginary axis and a straightforward calculation leads to the asymptotic result for V ≫ 1 If both q ± c diverge this result is valid for arbitrary finite values of q because x * is always inside the gap [−λ − 0 , λ + 0 ], and far from the edges. If, instead, one or both q ± c are finite, the saddle point x * may reach the boundary of the gap [−λ − 0 , λ + 0 ], and leave it for q < q − c or q > q + c . When this occurs the asymptotic approximation in Eq. (102) is no more valid.
The boundary values q ± c are finite or infinite depending on whether the singularities of the integrand function in L ′ (x) are integrable or not. More specifically, let us denote by k ± 0 the wave vectors at which λ ± (k) attain the minimum value, that is λ ± (k ± 0 ) = λ ± 0 . Then, in the neighborhoods of k ± 0 , where C is a positive constant, n = 1, 2, . . . and α > 0. Therefore, q ± c diverge if k ± 0 > 0 or if k ± 0 = 0 and d ≤ α, while they are finite if k ± 0 = 0 and d > α. Which is the case depends on the parameters of the quench T 0 , T, τ, t w . So, in general, this manifold of parameters is partitioned into phases distinguished by q ± c being finite or infinite. For what follows, and before discussing how the steepest descent calculation must be modified, it is useful to get some insight into the meaning of the two different cases. Recalling that λ + (k) and −λ − (k) are the inverse average heat absorbed or released, it is evident from Eq. (100) that if q coincides with the average heat density q , then x * ( q ) = 0, implying −L ′ (0) = q . Consequently, if q differs from q then x * (q) = 0 and acquire the meaning of the inverse average heat absorbed or released in new conditions, such that q would be the new average total heat density, exactly as the {B * i } of section II have been recognized to be the shifted external fields necessary to render the fluctuation m equal to the average magnetization per spin. In other words, a deviation of q from the average q brings in a shift by x * (q) of the inverse heats exchanged by the single modes. Accordingly, when x * (q) approaches the edges of the gap the above defined λ ± * 0 vanish, which means that the corresponding edge modes give an infinite contribution to the exchanged heat. This may happen either because the fluctuation q itself is infinite, or because the contributions of all the modes, other than the edge ones, can sum up at most to the finite amounts q ± c . In the latter case, if |q| goes above the thresholds |q ± c |, in order to make up for the finite difference |q| − |q ± c | a spike contribution must come from the edge modes at k ± 0 . This phenomenon is the condensation of fluctuations in the edge mode, since as anticipated in the Introduction and at the beginning of this section, the entire amount of the fluctuation above the critical threshold comes from a single degree of freedom. The mechanism of the transition can be recognized to be the same as that of Bose-Eintein condensation [38]. More in general, condensation of fluctuations has been found in widely different contexts such as information theory [39], finance [40] and statistical mechanics encompassing both equilibirum and out of equilibrium situations [41].
When one or both q ± c are finite, the steepest descent calculation must be modified. This case will be dealt with in the next subsection. Returning to the case in which Eq. (102) holds, neglecting subdominant terms and taking the ratio P (−Q)/P (Q), one obtains the following explicit form of the AF which is clearly reminiscent of Eq. (21) in sect. II. In order to expose the analogy, let us subtract equations (104) one from the other and take the average over k, obtaining where and similar expressions for λ * ± (q). Inserting into Eq. (105), we obtain where Hence, comparing with Eqs. (24) and (81), we can recognize the same structure: the prefactor λ − − λ + of the linear contribution (the same appearing in the single oscillator case) plays the role of the external average field B, while the nonlinear term [Ψ(−q) − Ψ(q)], arises from "free energy" contributions. As remarked above, the presence of this latter term is due to the k-dependence of ∆β k , which plays the same role as the i-dependence of B i .

A. Condensation transition
On physical grounds one can argue that condensation at finite q − c can occur only in cooling experiments, where the final temperature is lower than the initial one, while condensation at finite q + c can occur only in heating experiment. We do not have a proof of this, but numerical analysis of the conditions for condensation confirm this conjecture.
Let us thus concentrate on the case with q + c = ∞ and q − c finite, which occurs for k − 0 = 0 and d > α. The analysis of the opposite case in which q + c is finite and q − c = −∞, (or both finite, if such a case exists) is straightforward. When q − c is finite, Eq. (100) does not admit a solution with x * ∈ [−λ − 0 , λ + 0 ] whenever q < q − c . The problem is well known from the theory of Bose-Einstein condensation of an ideal gas of bosons [38] or from the "sticking" of the saddle point to a singularity in the solution of the spherical model of ferromagnetism [42]. It arises because when q < q − c = −L ′ (λ − 0 ) the steepest descent path cannot pass through the stationary point x * < −λ − 0 of φ(z; q) = zq + L(z) and must traverse the real axis at the gap edge z = −λ − 0 . Since φ ′ (−λ − 0 ; q) = q − q − c is negative for q < q − c , the steepest descent path at z = −λ − 0 bends toward the negative real axis forming the cusp characteristic of saddle point sticking at the gap edge (see Fig. 4). The large V behavior of P (Q, V ) is dominated by the neighborhood of the gap edge because φ(−λ − 0 ; q) < φ(x * ; q) leading to the asymptotic behavior for V ≫ 1: valid for q < q − c . As a matter of fact, this expression and Eq. (102) valid for q > q − c , hold for x * not to close to the gap edge −λ − 0 . That is, respectively, in the condensed phase and normal phase for q not too close to q − c . The analysis of the asymptotic behavior of P (Q, V ) as V ≫ 1 for all values of q, including close to q − c , requires some care. Details can be found in Appendix A.
Adding and subtracting λ − 0 q − c in the exponent of Eq. (110), P (Q; V ) in the condensed phase far from the critical point can be written as where Q − c = V q − c . Comparing the first term with the expression (55) for the distribution in the purely cooling process we have where, up to normalization factors, This means that negative fluctuations below the critical lower threshold Q − c condense into the cooling contribution of the k = 0 mode for the exceeding part (Q − Q − c ), while the contribution of all the other modes is locked onto Q − c . The comparison of the saddle point estimates of P (Q, V ), in the normal and in the condensed phase, with the "exact" numerical computation is illustrated in Fig. 5.

B. Phase diagram
In order to establish the "phase" structure, it is necessary to analyse the behavior of k ± 0 > 0. Let us first recall the results of Ref. [28], where the problem was analysed for the quench to T = 0. This is the simplest case because, as specified in section III B, heat can only be released and we have to deal only with the λ − (k) branch of the spectrum. Differentiating with respect to k the expression for λ − (k), given by the first line of Eq. (49), we get where C(k, τ, t w ) = b −2 ∆T e −2ω( k)(tw+τ ) αk α−1 is a positive quantity, while has the sign of (k α − E), with and vanishes at k α = E. Therefore, imposing the condition E(τ, t w ) = 0, in the (τ, t w ) plane there remains defined the critical line, given by such that above it k − 0 = 0, while below k − 0 > 0. Thus, k − 0 or, equivalently, q − c act as order parameters and below the critical line the system is in the normal phase, corresponding to q − c = −∞, while above it is in the condensed phase, corresponding to q − c finite, as illustrated in Fig. 6 (black line).
If T > 0, we must take into account both branches λ ± (k) and keep track of the two order parameters q ± c . The analytical search of the critical lines turns out to be quite complicated. So, we have looked for the absolute minima of λ ± (k) numerically. The resulting phase diagram for q − c is depicted in Fig. 6 for a few values of T , ranging from very low to almost equal to T 0 . When the system is equilibrated (T = T 0 ), there is no condensed phase. Fig. 6 shows that upon lowering T there appears a condensed phase which, starting from the far right, pronges toward the left eventually filling the entire region above the T = 0 critical line. The prominent qualitative difference between the T = 0 and the T > 0 phase diagrams for q − c is that in the latter case the transition driven by an increasing t w , for fixed τ , manifests reentrant behavior. The same computation for q + c does not show the existence of a condensed phase, namely q + c = ∞ all over the explored (τ, t w ) plane.
The reentrant behavior of the phase diagram can be interpreted as follows. Let us consider a quench at finite temperature, and let us fix τ . Then, for small enough t w , all oscillators are out of equilibrium and can exchange an arbitrary amount of heat with the bath. Therefore no condensation can take place in this region. Upon increasing t w , all oscillators do equilibrate but the slowest one, corresponding to the mode k = 0, that can account for large exchange of heat (above the threshold), leading to the condensation of fluctuations. Eventually, for large t w , all oscillators are equilibrated and again the normal phase is recovered.

VI. CONCLUSIONS
In this paper we have studied the fluctuations of the heat exchanged with the environment by a system of oscillators relaxing after a temperature quench. Focusing on the FT, we have investigated the relation between the deviations from linearity of the AF and the internal structure of the system, building on the analogy with the behavior of fluctuations in equilibrium when a symmetry is broken by nonuniform external perturbations. We have first analysed in great detail the case of a single oscillator, reducing it to the convolution of the two elementary processes in which heat can be only released or only absorbed. The inverse average heats released (λ − ) and absorbed (λ + ) in these processes turn out to be the basic objects underlying the behavior of the quantities of interest. In the one-oscillator case the AF is linear with slope given by the difference ∆β = λ − − λ + and, in the framework of the above mentioned analogy with the equilibrium problem, this quantity acts like an external field explicitly breaking the time reversal symmetry. By adding a second oscillator it starts to surface the role of the dishomogeneity of the external perturbation in determining deviations from linearity in the AF. It should be pointed out that in the dynamical problem the lack of homogeneity is due to the existence of different relaxation rates, related to the presence of degrees of freedom evolving on different time scales. This produces a temporal dishomogeneity, which induces a differentiation of the affinities ∆β 1 and ∆β 2 , analogous to the spatial heterogeneity generated by the external field {B i }.
This picture emerges most clearly in the case of a large number of oscillators. In this case, the system, although diagonal, presents nontrivial features, the most notable of which is the possible condensation of fluctuations [26][27][28]. Here, we have presented a detailed analytical study of the large deviation function of the heat distribution, with an accurate treatment of the crossover from the normal to the condensed phase, pointing out some interesting technical points involved in the application of the steepest descent method in the presence of a transition. In addition, we have mapped out the phase diagram of the condensation transition, discovering an unexpected reentrant behavior, in the (t w , τ ) plane.
Our analysis allowed us to establish a close correspondence with the structure of the AF for the paramagnet in equilibrium under the action of a nonuniform external field. In the equilibrium case it seems reasonable to make the statement that the linearity of the AF depends on whether the observable of interest M(σ) is conjugate to the external perturbation, the implication of which being that deviations from the FT are due to the lack of conjugation. Clearly, the interesting and challenging issue is to understand whether this way of looking at the AF and its deviations from linearity can be extended to the nonequilibrium case. The study we have presented in this paper is a first step in this direction, paving the way to futures investigations within the context of more complex interacting systems, showing slow relaxation and aging phenomena.

Appendix A
In this Appendix we derive the uniform asymptotic expansions of P (Q, V ) as V ≫ 1 valid for all values of q = Q/V , including at the condensation transition. As in the main text we shall assume q + c = ∞ and q − c finite so that the steepest descent path hits the gap lower edge at the finite value q = q − c . Close to the edge z = −λ − 0 the dangerous boundary term (1/V ) ln[λ − 0 +z] is no more negligible, thus the uniform approximation L(z), Eq. (99), to L(z; V ) as V ≫ 1 breaks up and shuld be replaced by (98). Therefore, to construct an uniform asymptotic approximation to P (Q; V ) as V ≫ 1 valid for all values of q the function L(z) must be regularised. Hence, in Eq. (A1) the function L(z) is replaced by where δ = O(1/V 1/d ) is an infrared cut-off [43]. Without loss of generality we can take δ = V −1/d because any proportionality constant can be absorbed into the O(1/V ) corrections in Eq. (A1). Notice that the boundary value of q separating the two phases becomes: Using λ + (k) − λ − 0 ∼ Ck α as k ≪ 1, the finite volume correction to the critical point reads in the large V limit: where Υ d = [2 d−1 π d/2 Γ(d/2)] −1 and Γ(x) is the Euler's gamma function. Replacing L(z; V ) into Eq. (97), neglecting the contribution from the O(1/V ) terms and taking z = −λ − 0 +t to move the end point of the cut (−∞; −λ − 0 ] on the negative real axis to t = 0, leads to: where φ(z; q) = zq + L δ (z). As discussed in the main text, the stationary point changes sign at q = q − δ , and is negative if q < q − δ and positive if q > q − δ . Far from the critical point, i.e., |q −q − δ | = O(1) as V ≫ 1, the steepest descent calculation is straightforward. If q > q − δ the saddle point lies on the positive t-axis and the steepest descent path can pass through it. Near t * the steepest descent path is parallel to the imaginary axis and as V ≫ 1, cfr. Eq. (102). In the opposite case q < q − δ the stationary point t * is negative and the steepest descent path cannot pass through it. Thus it must traverse the real axis at the gap edge t = 0. Since φ ′ δ (−λ − 0 ; q) is negative for q < q − δ , the steepest descent path bends toward the negative real axis at t = 0 forming the cusp characteristic of saddle point sticking at the gap edge (see Fig.  4). The asymptotic behavior for V ≫ 1 is dominated by the gap edge because φ(−λ − 0 ; q) < φ(−λ − 0 + t * ; q). In the neighborhood of t = 0 the steepest descent path is parallel to the negative real axis, hence using (A6), as V ≫ 1, cfr. Eq. (110). When |q − q − δ | ≪ 1 as V ≫ 1, q is very close to q − δ and the stationary point is at |t * | ≪ 1. In this region we can expand φ(−λ − 0 + t; q) in powers of t: The first three terms of the expansion must be retained; subsequent terms just give corrections to the leading asymptotic expansion and can be neglected.
Using the expansion (A9) the stationary point is at: .
The point t = 0 lies on the steepest ascent path issuing from the saddle point at t * and hence it will dominate the integral as V ≫ 1 whenever t * = O(1) as V ≫ 1. In this case only the paths Γ b contribute and a straightforward calculation leads to Eq. (A8).
To study the behaviour in the critical region V |φ ′′ 0 |t * = O(1) as V ≫ 1, it is convenient to notice that the steepest descent path Γ a + Γ b can be deformed into the imaginary axis x = 0 without changing the value of the integral. Taking t = iy/ V |φ ′′ 0 | the integral becomes: where ξ is defined in Eq. (A13), which with Eq. (A20) leads to Eq. (A16). However, notice that now ξ < 0 because we are in the condensed phase. Using the asymptotic expansion D −1/2 (−|ξ|) ∼ e ξ 2 /4 / |ξ| valid for −ξ ≫ 1 Eq. (A16) reduces to and matches asymptotically Eq. (A8). Thus, from Eq. (A18), the uniform asymptotic approximation to Eq. (A5) as V ≫ 1 for q ≤ q − δ reads Notice that close to q − δ , inside the critical region, P (Q; V ) is given by Eq. (A16) regardless of q being larger or smaller than q − δ . The only difference is the sign of ξ. As a consequence P (Q; V ) is regular across the transition between the two phases. The singularity shows up only in the strict V → ∞ limit.