Typical and rare fluctuations in nonlinear driven diffusive systems with dissipation

Typical and rare fluctuations in nonlinear driven diffusive systems with dissipation Pablo I. Hurtado∗ Instituto Carlos I de F´ısica Te´orica y Computacional, and Departamento de Electromagnetismo y F´ısica de la Materia, Universidad de Granada, Granada 18071, Spain A. Lasanta†and A. Prados‡ F´ısica Te´orica, Universidad de Sevilla, Apdo. de Correos 1065, Sevilla 41080, Spain (Dated: February 27, 2013) We consider fluctuations of the dissipated energy in nonlinear driven diffusive systems subject to bulk dissipation and boundary driving. With this aim, we extend the recently-introduced macroscopic fluctuation theory to nonlinear driven dissipative media, starting from the fluctuating hydrodynamic equations describing the system mesoscopic evolution. Interestingly, the action associated to a path in mesoscopic phase-space, from which large-deviation functions for macroscopic observables can be derived, has the same simple form as in non-dissipative systems. This is a consequence of the quasi-elasticity of microscopic dynamics, required in order to have a nontrivial competition between diffusion and dissipation at the mesoscale. Euler-Lagrange equations for the optimal density and current fields that sustain an arbitrary dissipation fluctuation are also derived. A perturbative solution thereof shows that the probability distribution of small fluctuations is always gaussian, as expected from the central limit theorem. On the other hand, strong separation from the gaussian behavior is observed for large fluctuations, with a distribution which shows no negative branch, thus violating the Gallavotti-Cohen fluctuation theorem as expected from the irreversibility of the dynamics. The dissipation large-deviation function exhibits simple and general scaling forms for weakly and strongly dissipative systems, with large fluctuations favored in the former case but heavily supressed in the latter. We apply our results to a general class of diffusive lattice models for which dissipation, nonlinear diffusion and driving are the key ingredients. The theoretical predictions are compared to extensive numerical simulations of the microscopic models, and excellent agreement is found. Interestingly, the large-deviation function is in some cases non-convex beyond some dissipation. These results show that a suitable generalization of macroscopic fluctuation theory is capable of describing in detail the fluctuating behavior of nonlinear driven dissipative media. I. INTRODUCTION Fluctuations are inherent to many physical phenomena, reflecting the hectic microscopic dynamics at macroscopic scales. In spite of their apparent random origin, essential physical information is encoded therein [1]. A classical example is the fluctuation-dissipation theorem, which relates the linear response of a system to an external perturbation to the fluctuation properties of the system in thermal equilibrium [2, 3]. More recently, the investigation of general properties of fluctuations in nonequilibrium steady states is opening new paths for understanding physics far from equilibrium [4– 14]. The study of fluctuation statistics of macroscopic observables provides an alternative path to obtain thermodynamic potentials, a complementary approach to the usual ensemble description. This observation, valid both in equilibrium [1] and nonequilibrium [7, 9], is most relevant in the latter case because no general bottom-up approach, connecting microscopic dynamics to macroscopic nonequilibrium properties, has been found yet. In this way, the large deviation function (LDF) controlling the ∗ph[email protected] †alasan[email protected] ‡[email protected] statistics of these fluctuations may play in nonequilibrium statistical mechanics a role similar to the equilibrium free energy [4, 10]. A central point for this emerging paradigm is the identification of the relevant macroscopic observables characterizing the out of equilibrium behavior of the system at hand. The system dynamics often conserves locally some magnitude (a density of particles, energy, momentum, charge, etc.), and the essential nonequilibrium observable is thus the current or flux sustained by the system when subject to boundary-induced gradients or external fields. Therefore, the understanding of current statistics in terms of the microscopic dynamics represents one of the main problems of nonequilibrium statistical mechanics, triggering an intense research effort in recent years. In this context, key results are the Gallavotti-Cohen fluctuation theorem [5], which relates the probability of observing a given current fluctuation ~ Jwith the probability of the reversed event −~ J, or the recently-introduced Isometric Fluctuation Relation [14], which relates the probability of any pair of isometric current fluctuations ( ~ J, ~ J0), with |~ J|=|~ J0|. These intriguing “symmetries” appear as a consequence of the invariance under time reversal of the underlying microscopic dynamics. These and other recent results [5–14] are however restricted to nonequilibrium “conservative” systems characterized by currents. On the other hand, many nonequilibrium systems are inherently dissipative, that is, they need a continuous inarXiv:1302.6544v1 [cond-mat.stat-mech] 26 Feb 2013 2 put of energy in order to reach a steady state. In this class of systems, the relevant macroscopic observable is not only the current: the dissipated energy is also expected to play a main role. These systems include granular media [16], dissipative biophysical systems [18], turbulent fluids [19], active matter [20], chemical reactions [21], population dynamics [22], etc. In general, this class comprises all sort of reaction-diffusion systems where dissipation, diffusion and driving are the main physical mechanisms. Fluctuations in dissipative media have been much less investigated, most probably because their physics is more complicated as a result of the irreversibility of their microscopic dynamics. In principle, most of the results for nonequilibrium steady states we have referred to in the previous paragraph are not applicable to dissipative systems, since they stem from the reversibility of the underlying microscopic dynamics in the conservative case. Therefore, a question naturally arises as to whether it is possible to extend some of these ideas to dissipative media. One of the main goals of the present paper is to give a (partial) answer to this question. In this work, we analyze both typical and rare fluctuations in nonlinear driven diffusive systems with dissipation. This is done by combining a suitable generalization of macroscopic fluctuation theory (MFT) [7] to the realm of dissipative media, and extensive numerical simulations of a particular albeit broad class of microscopic models. Our starting point is a general fluctuating balance equation for the (energy) density, with a drift term proportional to the spatial derivative of the current and a sink term. This mesoscopic description is expected to be valid for many driven dissipative media over a certain “hydrodynamic” time scale, much larger than the one characteristic of the microscopic dynamics. Over the fast (microscopic) time scale, the system forgets the initial conditions and relaxes to a local equilibrium state in which all the properties of the system become functionals of a few “hydrodynamic” fields, here the density, the current and the dissipation. Afterwards, over the much slower hydrodynamic time scale, the system eventually approaches the steady state following the mesoscopic balance equation. We focus on the fluctuations of the system in this nonequilibrium steady state, in which the dissipation and the injection of energy balance each other. By using this fluctuating hydrodynamic picture together with a path integral formulation, we derive a general form for the action associated to a history of the density, current and dissipation fields (that is, a path in mesoscopic phase space). Remarkably, this action takes the same form as in conservative nonequilibrium systems [7–14], simplifying the analysis in the dissipative case. This is both an important and a surprising result, which stems from the quasi-elastic character of the underlying microscopic dynamics in the large system size limit. This quasi-elasticity is necessary in order to have a balanced competition between diffusion and dissipation at the mesoscopic level. From the derived action functional, and using the recently-introduced additivity conjecture [8, 11–13], a general form for the LDF of the dissipated energy is derived, with a “Lagrangian” including second order derivatives. Therefrom, we derive the Euler- Lagrange equation (a fourth-order differential equation) for the optimal fields responsible of an arbitrary fluctuation. This Lagrangian variational problem can be mapped onto an equivalent Hamiltonian problem (four coupled first-order differential equations) which turns out to simplify the analysis. We use this Hamiltonian picture to analyze in detail three different limits, namely small fluctuations around the average for arbitrary dissipation coefficient, and the whole spectrum of fluctuations (typical and rare) for weakly- and strongly-dissipative systems. The statistics of typical (that is, small) fluctuations is gaussian as expected from the central limit theorem. However, strong separation from gaussian behavior is observed for rare fluctuations, with a distribution which shows no negative branch, thus violating the Gallavotti-Cohen fluctuation theorem as otherwise expected from the irreversible character of microscopic dynamics. We study in general the weakly-dissipative system limit using a singular perturbation expansion. This yields a simple scaling form for the dissipation LDF, showing that large dissipation fluctuations are favored in this weakly-dissipative limit, with a LDF which extends over a broad regime and decays slowly in the far positive tail. On the other hand, a different perturbative analysis in the strongly-dissipative system limit can be carried based on the formation of boundary energy layers in this limit for all fluctuations, which effectively decouples the system in two almost-independent parts. This analysis shows that large dissipation fluctuations are heavily suppressed in this limit, as opposed to the weakly-dissipative regime result. We apply this theoretical scheme to a general class of d-dimensional dissipative lattice models with stochastic microscopic dynamics, for which the hydrodynamic fluctuating picture used above as starting point can be demonstrated in the large system size limit [39] (we will focus here in one dimension for simplicity). In these models there is one particle at each lattice site, characterized by its energy. Dynamics is stochastic and proceeds via collisions between nearest neighbors, at a rate which depends on the energy of the colliding pair. In a collision, a certain fraction of the pair energy is dissipated, and the remaining energy is randomly distributed within the pair. This mechanism gives rise to a nonlinear competition between diffusion and dissipation in the macroscopic limit, provided that the microscopic dissipation coefficient scales adequately with the system size. This class of models represents at a coarse-grained level the physics of many reaction-diffusion systems of technological as well as theoretical interest. In particular, when the colliding pair is chosen completely at random, independently of the value of its energy, the Kipnis-Marchioro-Presutti (KMP) model [23] for heat conduction is recovered in the conservative case. The KMP model plays a main role in nonequilibrium statistical physics as a touchstone to test 3 theoretical advances [5, 6, 8, 11–15, 23]. Our general class of models contains the essential ingredients characterizing most dissipative media, namely: (i) diffusive dynamics, (ii) bulk dissipation, and (iii) boundary injection. The chances are that our results remain valid for more complex dissipative media described at the mesoscopic level by a similar evolution equation. Here we report analytical and simulation results for the statistics of the dissipated energy in this general class of models using both standard simulations and an advanced Monte Carlo method [27]. The latter allows the sampling of the tails of the distribution, and implies simulating a large number of clones of the system. The plan of the paper is as follows. Section II describes a suitable generalization of macroscopic fluctuation theory to nonlinear driven dissipative media. The largedeviation statistics of the dissipated energy, a central observable in this type of systems, is investigated here within both the Lagrangian and Hamiltonian equivalent frameworks. Section III is devoted to the detailed study of different asymptotic behaviors within the Hamiltonian formulation, which turns out to simplify the analysis. In section IV we define a general class of microscopic lattice models whose stochastic dynamics is dissipative [25, 39]. The theoretical framework developed in the previous sections is applied to this family of models in section V, and the LDF for the dissipated energy is explicitly worked out. The analytical predictions are compared to extensive numerical simulations of the microscopic models, and a very good agreement is found. A summary of the main results of the paper, together with a physical discussion thereof, is given in sec. VI. Finally, the appendix deals with some technical details that, for the sake of clarity, we have preferred to omit in the main text. II. MACROSCOPIC FLUCTUATION THEORY FOR DRIVEN DISSIPATIVE SYSTEMS In this work, we will analyze a general class of systems whose dynamics at the mesoscale is described by the following fluctuating evolution equation ∂tρ(x, t) = −∂xj(x, t) + d(x, t).(2.1) We focus here in one dimension for simplicity, but our analysis can be carried out in an equivalent manner in ddimensions. In eq. (2.1), ρ(x, t), j(x, t) and d(x, t) are the density, current and dissipation fields, respectively, and t and x∈[−1/2,1/2] are the macroscopic time and space variables, obtained after a diffusive scaling limit such that x= ˜x/L and t=˜ t/L2, with ˜xand ˜ tthe microscopic space and time variables and Lthe system length. These coarse-grained spatial and temporal scales emerge from a suitable continuum limit of the underlying microscopic dynamics [39]. The current field is a fluctuating quantity, and can be written as j(x, t) = −D(ρ)∂xρ(x, t) + ξ(x, t).(2.2) The first term is Fourier’s law, where D(ρ) is the diffusivity (which might be a nonlinear function of the local density), and ξ(x, t) is the current noise that is gaussian and white, hξ(x, t)i= 0,hξ(x, t)ξ(x0, t0)i=σ(ρ) Lδ(x−x0)δ(t−t0), (2.3) with σ(ρ) being the so-called mobility. This gaussian fluctuating field is expected to emerge for most situations in the appropriate mesoscopic limit as a result of a central limit theorem: although microscopic interactions for a given model can be highly complicated, the ensuing fluctuations of the slow hydrodynamic fields result from the sum of an enormous amount of random events at the microscale which give rise to gaussian statistics, with an amplitude of the order of L−1/2, in the mesoscopic regime in which eq. (2.1) emerges. On the other hand, the dissipation field d(x, t) is d(x, t) = −νR(ρ(x, t)),(2.4) where νis the macroscopic dissipation coefficient, and R(ρ) is a certain function of the density ρ. For the calculations which follow throughout this section, it is useful to introduce a new variable y, such that y=R(ρ),(2.5a) d(x, t) = −νy(x, t).(2.5b) The dissipation field is present at the mesoscopic level because the microscopic stochastic dynamics of the models of interest dissipates some energy, that is, we have the equivalent of a microscopic restitution coefficient α, so that the amount of dissipated energy is proportional to 1 −α. The macroscopic dissipation coefficient νis thus proportional to 1 −α. Note however that there is no noise term in eq. (2.4), so the local fluctuations of the dissipation field are enslaved to those of the density ρ(x, t). The physical reason for this behavior is that the microscopic dynamics must be quasi-elastic in order to ensure that dissipation and diffusion take place over the same time scale in the thermodynamic limit. Typically, 1−αmust scale as L−2that is the order of magnitude of the diffusive term in a system of length L[25, 39]. The boundary conditions for eq (2.1) depend on the physical situation of interest. For instance, we may consider that the system is kept in contact with two thermal reservoirs at x=±1/2, at the same temperature T, so ρ(±1/2, t) = T. In that case, the system eventually reaches a steady state in the long time limit, for which the injection of energy through the boundaries and the dissipation balance each other. The stationary average (macroscopic) solution of (2.1) verifies j0 av(x) + νR(ρav(x)) = 0, jav(x) = −D(ρav(x))ρ0 av(x), (2.6) where the prime indicates spatial derivative. The first equation in (2.6) follows from (2.1), and the second one 4 is Fourier’s law for the averages. Equivalently, a closed second-order equation for ρmay be written, d dx [D(ρav)ρ0 av] = νR(ρav),(2.7) with the boundary conditions ρav(±1/2) = T. Equations (2.6) and (2.7) can be also written for the variable y introduced in eq. (2.5a), j0 av(x) + νyav(x)=0, jav(x) = −ˆ D(yav(x))y0 av(x), (2.8) with ˆ D(y) = dy dρ−1 D(ρ),(2.9) since jav(x, t) = −ˆ D(yav)∂xy(x, t).(2.10) Thus ˆ Dis an “effective” diffusivity: it is the factor multiplying the spatial gradient when writing Fourier’s equation in terms of the new variable y. Equation (2.8) can also be summarized in a second order differential equation for yav, hˆ D(yav)y0 avi0=νyav, yav(±1/2) = R(T).(2.11) Interestingly, it can be shown (see below) that ˆ Dis constant, independent of y, whenever y=R(ρ) depends algebraically in ρ, a case we will study in detail in section V. This observation considerably simplifies the subsequent analysis. The probability of observing a history {ρ(x, t), j(x, t)}τ 0 of duration τfor the density and current fields, starting from a given initial state, can be written now as a path integral over all the possible realizations of the current noise {ξ(x, t)}τ 0, weighted by its gaussian measure, and restricted to those realizations compatible with eq. (2.1) at every point of space and time [12]. This probability hence obeys a large deviation principle of the form [4, 7– 13, 25] P({ρ, j}τ 0)∼exp (+LIτ[ρ, j]) ,(2.12) with a rate functional [7, 9] Iτ[ρ, j] = −Zτ 0 dt Z1/2 −1/2 dx [j+D(ρ)∂xρ]2 2σ(ρ)(2.13) with ρ(x, t) and j(x, t) coupled via the balance equation (2.1), and the dissipation d(x, t) given in terms of ρ(x, t) by (2.4). Equation (2.13) expresses the gaussian nature of the local current fluctuations around its average (Fourier’s law) behavior. The functional in (2.13) is the same as in the conservative case (that is, with no bulk dissipation), due to the quasi-elasticity of the microscopic dynamics, which makes the current noise be the FIG. 1. (Color online) Convergence of the space&timeintegrated dissipation to its ensemble value for many different realizations, and sketch of the probability concentration as time increases, associated with the large deviation principle, eq. (2.15). only relevant one in the hydrodynamic description, see discussion in section IV [39]. We focus now on the fluctuations of the dissipated energy, integrated over space and time d=−1 τZτ 0 dt Z1/2 −1/2 dx d(x, t) (2.14) =ν τZτ 0 dt Z1/2 −1/2 dx R(ρ(x, t)) >0, where we have introduced a minus sign for the sake of convenience, in order to make dpositive. As discussed above, this is a fundamental observable to understand the statistical physics of driven dissipative media. The probability of such a fluctuation Pτ(d) scales in the longtime limit as Pτ(d)∼exp [+τL G(d)] , G(d) = 1 τmax ρ,j Iτ[ρ, j]. (2.15) This defines a new large deviation principle for d, see Fig. 1, such that G(d) is obtained from Iτ[ρ, j] via a saddlepoint calculation for long times (that is, it follows from the contraction of the original rate function Iτ[10]). The optimal fields ρ0(x, t;d), j0(x, t;d) which are the solution of the variational problem (2.15) must be consistent with the prescribed value of the dissipated energy din (2.14), and are also related by the balance equation (2.1), supplemented with (2.4) and the appropriate boundary conditions. These optimal fields can be interpreted as the ones adopted by the system to sustain a long-time fluctuation of the space&time-integrated dissipation d. For the sake of simplicity, we have not explicitly introduced in our notation the parametric dependence of the LDF G(d) and the associated optimal profiles on the boundary temperature T, though this should be borne in mind for latter reference. 5 A. The constrained variational problem We now assume that these optimal profiles do not depend on time. In conservative systems, this conjecture has been shown [7] to be equivalent to the additivity principle recently introduced to study current fluctuations in diffusive media [8]. The validity of this additivity scenario has been recently confirmed in extensive numerical simulations for a broad interval of fluctuations [11, 14], though it may eventually break down for extreme fluctuations via a dynamic phase transition [30, 31]. As we will see below, the applicability of this generalization of the additivity conjecture to dissipative systems is well supported by numerical evidence. Under this simplifying hypothesis, the fluctuating balance equation (2.1) reduces to j0(x) + νy(x)=0, y(x) = −j0(x)/ν, (2.16) making use of the variable ydefined in eq. (2.5a). Moreover, we can integrate over time in the definition (2.14) of the integrated dissipation d, d=νZ1/2 −1/2 dx y(x),(2.17a) or, equivalently, d=−Z1/2 −1/2 dx j0(x) = j(−1/2) −j(1/2) >0.(2.17b) In this way, by using the additivity hypothesis we can eliminate ρ(x) and write G(d) in terms of only one variable as G(d) = −min j(x)S[j],with S[j] = Z1/2 −1/2 dx L(j, j0j00), (2.18a) L(j, j0, j00) = [j−ˆ D(−j0/ν)j00 ν]2 2ˆσ(−j0/ν),(2.18b) where ˆ Dis the effective diffusivity defined in eq. (2.9), and ˆσis the mobility, defined in eq. (2.3), both written in terms of y=−j0/ν. The function L(j, j0, j00) is a generalized Lagrangian with dependence on first and also second order derivatives, see the appendix. We have to find the optimal current profile j0(x;d), that is, the solution of the variational problem (2.18), with the constraint that the integrated dissipation dhas a definite value, as given by (2.17). Therefore we must use the Lagrange multiplier procedure [32, 33], that is, look for an extremum of Sλ[j] = S[j]−λZ1/2 −1/2 dx (j0+d) (2.19) =Z1/2 −1/2 dx Lλ(j, j0, j00) (2.20) where Lλ(j, j0, j00) = L(j, j0, j00)−λ(j0+d),(2.21) with λbeing the Lagrange multiplier. The extremum of Sλfollows from two conditions: (i) δSλ= 0, and (ii) ∂Sλ/∂λ = 0. The first condition implies d2 dx2∂Lλ ∂j00 −d dx ∂Lλ ∂j0+∂Lλ ∂j = 0,(2.22) which is the Euler-Lagrange equation for a Lagrangian Lλcontaining second order derivatives (see the appendix). Condition (ii) leads to the constraint on the integrated dissipation, given by Eq. (2.17). The boundary conditions for the Euler-Lagrange equation are j0(±1/2) = −νR(T), pλj(±1/2) = 0.(2.23) We have introduced the generalized momentum pjconjugate of j, for the new Lagrangian Lλ, as pλj =∂Lλ ∂j0−d dx ∂Lλ ∂j00 .(2.24) These boundary conditions arise from (1) the values of the density at the boundaries, which are prescribed, ρ(±1/2) = T, and (2) the condition δSλ= 0, which provides the additional needed conditions when there are not enough values of the variables fixed at the boundaries (see the appendix, and also [33]). B. Mapping the constraint to boundary conditions Taking into account the relation between Lλand L, Eq. (2.21), the generalized momentum pλj verifies pλj =pj−λ, (2.25) where pjis the generalized momentum for the Lagrangian L, that is, pj=∂L ∂j0−d dx ∂L ∂j00 .(2.26) Moreover, the Euler-Lagrange equation (2.22) implies that d2 dx2∂L ∂j00 −d dx ∂L ∂j0+∂L ∂j = 0,(2.27) that is, we also obtain the Euler-Lagrange equation corresponding to the original Lagrangian L. The definition of pjin eq. (2.26) guarantees that p0 j=∂L/∂j, as in the case of the usual variational problem with a Lagrangian including only first-order derivatives. Now, the boundary conditions can be written as j0(±1/2) = −νR(T), pj(±1/2) = λ, (2.28) 6 which follow from Eqs. (2.23) and (2.25). The above result imply that our constrained variational problem can be mapped onto a unconstrained variational problem with the original Lagrangian L, its associated Euler- Lagrange equation (2.27) and the boundary conditions (2.28). The unknown value λfor the generalized momentum pjat the boundaries must be determined by imposing the prescribed value of the integrated dissipation, as given by Eq. (2.17), that is, λ=λ(d). In particular, λ= 0 is equivalent to imposing no restrictions on the integrated dissipation, so that we should recover (as we will see later) the average profiles and dissipation in this case. In this sense, a non-zero value of λ=pj(±1/2) is a measure of the departure from the average hydrodynamic behaviour. On physical grounds, we expect the corresponding optimal density profile to be an even function of x, because of the symmetry of our system around the center x= 0. In fact, the Euler-Lagrange equation (2.27) admits solutions with well-defined parity. Since the Lagrangian has the symmetry property L(−j, j0,−j00) = L(j, j0, j00), Eq. (2.27) has solutions with jbeing an odd function of x, which implies that y(and therefore ρ) is an even function of x. From now on, we will restrict ourselves to these symmetric solutions of the variational problem. Thus, Eq. (2.17) reduces to d= 2j(−1/2; d) = −2j(1/2; d),(2.29) so the boundary conditions for the Euler-Lagrange equation boil down to j0(±1/2; d) = −νR(T), j(−1/2; d) = −j(1/2; d) = d/2, (2.30) i.e. much simpler than Eq. (2.28). This follows from pj being an even function of xfor the solutions with welldefined parity we are considering. In this way, symmetry considerations lead to the simpler boundary conditions (2.30), which in turn allow us to get rid of the Lagrange multiplier λ. In summary, we have mapped our original variational problem with the subsidiary condition that the dissipation has a given value to an unconstrained variational problem, with the original Lagrangian L(j, j0, j00) and prescribed values of jand j0at the boundaries. Once the optimal current profile is obtained, the optimal density profile can be calculated from the balance equation (2.16). Of course, the density profile so obtained obeys the boundary conditions ρ(±1/2; d) = T. It must be stressed that the Euler-Lagrange equation (2.27) with boundary conditions (2.30) gives the correct solution to the constrained variational problem when the optimal profiles have a well-defined parity. Nevertheless, one cannot rule out the existence of symmetry-breaking solutions without well-defined parity, since in general a variational problem may have multiple solutions [32]. In that case, one must solve the more complex variational problem comprising the Euler-Lagrange equation (2.27) with the boundary conditions (2.28), where the Lagrange multiplier λ=λ(d) is determined by imposing the constraint (2.17). We note however that numerical evidence strongly supports the validity of symmetric solutions (see below). The LDF G(d) depends on dand Tthrough the boundary conditions (recall that its T-dependence has been omitted in notation for simplicity). The derivation of the Euler-Lagrange equation (2.27) (see the appendix) shows that, for the solutions with well-defined parity, δG(d) = pj(1/2)δd + 2pj0(1/2)νdR(T) dT δT, (2.31) where pjis the generalized momentum conjugate of j, defined in Eq. (2.26), and pj0is the generalized momentum conjugate of j0, pj0=∂L ∂j00 ,(2.32) which is an odd function of xfor the solutions with welldefined parity. In this way Eq. (2.31) offers a geometric interpretation for the values of the generalized momenta at the boundaries, as they are directly related to the partial derivatives of the LDF, ∂G ∂d =pj(1/2),∂G ∂T = 2νdR(T) dT pj0(1/2).(2.33) C. A Hamiltonian formulation of the problem We will not write the detailed form of the general fourth-order differential equation (2.27) for the optimal profile j(x;d), since it is not particularly illuminating. Instead we now write a set of four equivalent first order differential equations arising in the equivalent “Hamiltonian” description. In the following, we sketch the procedure to introduce the Hamiltonian for a Lagrangian with higher-order derivatives [32, 33], adapted to the present case. As the Euler-Lagrange equation is a fourth-order differential equation, we should have two canonical coordinates and their two corresponding canonical momenta. The first canonical coordinate is the current j, and we choose the second one to be y, which is proportional to j0, as given by Eq. (2.16). This choice is suggested by the structure of the Lagrangian in eq. (2.18). It is worth recalling that the density profile can be directly obtained from yby making use of its definition, Eq. (2.5a). Next, we introduce the canonical momenta pyand pjconjugate to yand j, respectively. The momentum pjhas been defined in Eq. (2.26), and pyis given by py≡ −ν∂L ∂j00 ,(2.34) which follows from the definition of pj0, Eq. (2.32), and Eq. (2.16) for y. The Hamiltonian is then introduced in the usual way, H ≡ y0py+j0pj−L ≡ y0py−νypj−L.(2.35) 7 After some algebra, we get H=1 2Q(y)p2 y−ˆ D−1(y)jpy−νypj,(2.36a) Q(y)≡ˆσ(y) ˆ D2(y),(2.36b) where we have defined the auxiliary function Q(y), with Q(y)>0 for all y. We have also made use of eq. (2.16), of the expression for pywhich follows from its definition (2.34), py=ˆ D(y)j+ˆ D(y)y0 ˆσ(y),(2.37) and of the Lagrangian L=ˆσ(y)p2 y 2ˆ D2(y)=1 2Q(y)p2 y,(2.38) written in terms of the canonical variables, with the aid of eq. (2.37). As usual, His a function of (y, j, py, pj), which satisfy the following set of four ”canonical” firstorder differential equations, y0=∂H ∂py =Q(y)py−ˆ D−1(y)j, (2.39a) j0=∂H ∂pj =−νy, (2.39b) p0 y=−∂H ∂y =−dQ(y) dy p2 y 2+dˆ D−1(y) dy jpy+νpj,(2.39c) p0 j=−∂H ∂j =ˆ D−1(y)py(2.39d) that are equivalent to the fourth-order Euler-Lagrange equation (2.27). Note that, as is usual in physics and in order not to clutter our formulae, we have dropped the subindex 0 for the optimal profiles, which are now solutions of the above canonical equations; the same notation is used for the canonical variables in the Hamiltonian and for the solutions of Hamilton’s equations. On the other hand, as the Hamiltonian does not depend explicitly on x, it is a first integral of the system (2.39): H= const. over any of its solutions. This property may be used to simplify the integration of the system. In general, for a given value of the dissipation d, we have to solve the system of equations (2.39) with the boundary conditions y(±1/2) = R(T), j(−1/2) = −j(1/2) = d/2,(2.40) which follow from Eqs. (2.30) and (2.16). Again, we have to look for solutions of Eq. (2.39) with well-defined parity, that is, yand jare even and odd functions of x, respectively (and, therefore, pyis odd and pjeven). The solution of these canonical equations is then inserted into the expression of the LDF G(d), which can be written in terms of the canonical variables as G(d) = −Z1/2 −1/2 dx L=−1 2Z1/2 −1/2 dx Q(y)p2 y,(2.41) by combining eqs. (2.18) and (2.38). In this way, we obtain the LDF for an arbitrary value of the integrated dissipation dwithin the Hamiltonian formulation of the variational problem. Equation (2.41) shows clearly that the most probable (average) profiles correspond to a solution with py= 0 for all x, for which G(d) vanishes. By substituting py= 0 in Eqs. (2.39c)-(2.39d), we also have that pj= 0 for all x. Moreover, Eqs. (2.39a) and (2.39b) simplify to eq. (2.8), that is, the average profiles are reobtained. Therefore, there is always a solution of the canonical equations (2.39) with identically vanishing canonical momenta, which corresponds to the average solution of the hydrodynamic equation (2.8) [40]. These average hydrodynamic profiles {ρav, jav}lead to the average value of the integrated dissipation dav =νZ1/2 −1/2 dx R(ρav) = νZ1/2 −1/2 dx yav(x) = 2jav(−1/2).(2.42) This discussion is consistent with the one below Eq. (2.28), which was done within the framework of the equivalent Lagrangian description. Fluctuations involve non-zero values for the canonical momenta, whose magnitude is then a measure of the departure from the average behaviour (d−dav)/dav. III. ANALYSIS OF THE LDF IN SOME LIMITING CASES In the following subsections, we further analyze the form of the LDF in certain limits of interest, for which some general results can be obtained. First, we focus on the behavior of G(d) for small fluctuations around the average, where a quadratic shape of the LDF is expected (corresponding to gaussian fluctuations). We then analyze the limit of weakly-dissipative systems, ν1, for which an adequate perturbative expansion allows us to obtain a non-trivial and interesting scaling form for the LDF. Finally, we consider the opposite limit of stronglydissipative systems, ν1, for which a different scaling for the LDF is found. A. Small fluctuations around the average As the average behavior corresponds to the particular solution of the canonical equations corresponding to vanishing momenta pj= 0, pρ= 0, small fluctuations can be 8 thus analyzed by assuming that the canonical momenta are small. Let us define the dimensionles parameter =d−dav dav (3.1) to measure the separation from the average integrated dissipation dav. As we have just discussed, the canonical momenta vanish for = 0. We write y=yav +∆y, j =jav +∆j, py=∆py, pj=∆pj, (3.2) and linearize Eqs. (2.39) around the average solution, that is, we only retain terms linear in . Then, ∆y0=Q(yav)∆py−ˆ D−1(yav)∆j−jav dˆ D−1(yav) dyav ∆y, (3.3a) ∆j0=−ν∆y, (3.3b) ∆p0 y=jav dˆ D−1(yav) dyav ∆py+ν∆pj,(3.3c) ∆p0 j=ˆ D−1(yav)∆py.(3.3d) The boundary conditions for these equations are ∆y(±1/2) = 0, ∆j(−1/2) = −∆j(1/2) = dav/2. The solution of this system of equations must be inserted in the large deviation function (2.18). Using the expression (2.41) for the LDF, it follows that G(d)∼ −2 2Z1/2 −1/2 dx Q(yav)∆p2 y,(3.4) for small fluctuations of the dissipation around the average. Taking into account (3.1) and the large deviation principle (2.15), eq. (3.4) means that the probability of such small fluctuations of the integrated dissipation dis approximately gaussian, Pτ(d)||1 ∝exp −Lτ (d−dav)2 2d2 avΛ2 ν(3.5) with Λ2 νgiven by Λ2 ν= Z1/2 −1/2 dx Q(yav)∆p2 y!−1 .(3.6) In this way, the gaussian estimation for the standard deviation of the dissipation, by comparing (3.5) to (3.6), is given by χ≡davΛν/√τL. In order to make a more detailed study of the LDF, concrete functional dependences of the diffusivity D, the mobility σand the dissipation R on the density ρmust be considered. This is done in the following sections of the paper, where we will consider a broad family of models for which the transport coefficients can be explicitly obtained. On the other hand, it is important to notice that gaussian statistics is only expected for small fluctuations around the average dissipation. In general, the solution of the variational problem given by the integration of eq.(2.39), when inserted into (2.41), will give rise to non-gaussian statistics (that is, a non-quadratic dependence of the LDF) for an arbitrary fluctuation of the dissipated energy d. B. Weakly-dissipative systems, ν1 We proceed now by analysing the canonical equations (2.39) in the limit ν1. Unsurprisingly, a regular perturbation expansion in powers of νbreaks down, since it is not possible to impose the necessary boundary conditions for the current. This singularity of the elastic limit was to be expected on a physical basis, as it is not possible to obtain the behavior of weakly dissipative systems (ν1) as a correction around the conservative case ν= 0, for which ρ(x) = Tand j(x) = 0. Therefore, a singular perturbation analysis should be done, looking for a suitable rescaling of the variables for ν1. Equation (2.8) for the averages implies that yav =R(T) + O(ν), jav =−νR(T)x+O(ν2), ν 1. (3.7) The average current vanishes linearly in νin the limit ν→0+, as expected. Moreover, the average dissipation, obtained by combining eqs. (2.42) and (3.7), is given by dav =νR(T) + O(ν2).(3.8) Therefore, it is sensible to propose the following rescaling of variables j(x) = νψ(x), pj(x) = Πψ(x) ν,(3.9) which is consistent with the canonical equations (2.39), since ψ0=1 νj0=1 ν ∂H ∂pj =∂H ∂Πψ ,(3.10a) Π0 ψ=νp0 j=−ν∂H ∂j =−∂H ∂ψ ,(3.10b) with the same Hamiltonian H. In other words, eq. (3.9) defines a “canonical transformation” from the pair of canonical conjugate variables {j, pj}to {ψ, Πψ}, a transformation that heals the singular behavior in the ν→0+ limit. The Hamiltonian can be now written as H=1 2Q(y)p2 y−yΠψ−νˆ D(y)−1ψpy(3.11) in the rescaled variables. Notice that the transformation introduced is essential to obtain the correct “dominant balance” [42] to the lowest order. In particular, before the rescaling, the term proportional to ypjwas of the order of νand the term proportional to jpywas of the order of unity; after the rescaling the orders of magnitude 9 are interchanged, the term proportional to yΠψis of the order of unity while the term proportional to ψpyis of the order of ν. We now start from the zero-th order rescaled Hamiltonian by putting ν= 0 in eq. (3.11), H0=1 2Q(y)p2 y−yΠψ.(3.12) from which we we arrive at y0=∂H0 ∂py =Q(y)py,(3.13a) p0 y=−∂H0 ∂y =−1 2 dQ(y) dy p2 y+ Πψ,(3.13b) ψ0=∂H0 ∂Πψ =−y , (3.13c) Π0 ψ=−∂H0 ∂ψ = 0 .(3.13d) In order not to clutter our formulas, we do not introduce a different notation for the canonical variables, although the approximate canonical equations (with H0) are different from the exact ones (with H). We have only to remember that our results are valid to the lowest order in ν. The canonical equations (3.13) have to be solved with the boundary conditions y(±1/2) = R(T), ψ(−1/2) = −ψ(1/2) = ∆/2, (3.14) where ∆ = d ν=R(T)d dav (3.15) is assumed to be of the order of unity, that is, d=O(ν) or d/dav =O(1). Thus, our rescaling allows us to obtain a solution for the optimal profiles for the density ρ, by inverting the relation y=R(ρ), and the current j= νψ, for integrated dissipations dvery different from its average value dav. It is worth noticing that ψis a cyclic variable and its conjugate momentum is thus constant, Πψ≡Πψ0= const., see Eq. (3.14); this fact allows us to obtain a closed first order differential equation for y(x) in the ν1 limit. Moreover, by recalling Eq. (2.33), we have that Πψ0=∂G ∂∆,(3.16) which gives the physical interpretation of this first integral of the approximate canonical equations: it is the partial derivative of the LDF with respect to the rescaled dissipation. The Hamiltonian H0is also constant, since it does not depend explicitly on x, and combining (3.12) and (3.13a), y02= 2Q(y)(H0+yΠψ0), y(±1/2) = R(T).(3.17) Once this is solved, the rescaled current ψcan be obtained from (3.13c) ψ0=−y, ψ(−1/2) = −ψ(1/2) = R(T)d 2dav ,(3.18) so that the two constants H0and Πψ0will be given in terms of the temperature Tand d/dav. There are no more constants to be adjusted in the solution of eqs. (3.17) and (3.18) due to the parity properties of (y, ψ): yis and even function of xand ψis an odd function of xin the interval [−1/2,1/2]. Of course, the average profiles yav(x) and jav(x) are reobtained from the canonical equations by putting Πψ= 0 and py= 0 therein. Equation (3.12) implies that H0= 0 over the average profiles. The simple form of the differential equation (3.17) allows us to infer some of the properties of the optimal profile y(x) associated to a given dissipation fluctuation in the limit of weakly-dissipative systems. First, notice that in general the solution of eq. (3.17) will be non-monotonic, exhibiting extrema in the interval x∈[−1 2,1 2]. Moreover, taking into account that the function Q(y) is positive defined, it follows that the profile at the extrema will take an unique value y0≡ − H0 Πψ0 (3.19) Note that Πψ06= 0 for d6=dav and, moreover, it must have a different sign that H0, i.e. sgn(Πψ0)6= sgn(H0), since y(x)>0∀x. Therefore, the optimal profile y(x) can only have a single extremum (minimum or maximum) [41], that is located at x= 0 because of symmetry reasons. By rewriting eq. (3.17) as y02= 2Q(y)H01−y y0(3.20) we conclude that the constant H0and y0−y(x) must have the same sign ∀x∈[−1 2,1 2]. Thus, for H0>0 the profile y(x) has a single maximum, y(x)> y(±1/2) = R(T)∀x, and thus d>dav. On the other hand, H0<0 implies a single minimum, y(x)< R(T)∀xand d < dav. All these properties are confirmed below for particular examples, both analytically and numerically. Interestingly, the leading behavior for the LDF can be also easily obtained in terms of the first integrals H0and Πψ0. In fact G(d)∼ −1 2Z1/2 −1/2 dx Q(y)p2 y=−1 2Z1/2 −1/2 dx y02 Q(y), (3.21) and making use of eq. (3.17), G(d)∼ − H0+ Πψ0Z1/2 −1/2 dx y(x)! =−(H0+ Πψ0∆) = −H01−∆ y0.(3.22) Thus, the remaining task consists in writing the constants H0and y0(or equivalently H0and Πψ0) in terms of the integrated dissipation dand the temperature at the boundaries T. Once this is done, the LDF follows from the simple expression given by eq. (3.22). Furthermore, we may obtain bounds for the profile extremum y0 16 ∆j0=−ν∆y, (5.19b) ∆p0 y=ν∆pj,(5.19c) ∆p0 j=1 ˆ D∆py,(5.19d) which particularizes eq. (3.3) for our family of models. The boundary conditions are ∆y(±1/2) = 0, ∆j(−1/2) = −∆j(1/2) = dav/2. The solution of this system must be inserted into eq. (3.6), which gives the variance of the gaussian distribution, χ2≡d2 avΛ2 ν/Lτ. The canonical momentum ∆pyis directly obtained by integrating eqs. (5.19c) and (5.19d), ∆py=Ksinh xrν ˆ D,(5.20) since ∆pymust be an odd function of xas a consequence of ybeing even. The constant Kis to be determined with the aid of the boundary conditions, but this can only be done after solving eqs. (5.19a) and (5.19b), having previously inserted (5.20) into them. Substitution of eq. (5.20) into eq. (3.6) gives Λ2 ν=K2 Z1/2 −1/2 dx Q(yav) sinh2xrν ˆ D!−1 .(5.21) In order to evaluate the integral, eq. (5.14) for Q(y) must be used, Q(y)∝yγ, with the parameter γbeing a function of β, 1 < γ ≤2. We now analyze the simplest choice β= 0, that is, γ= 2, that corresponds to the dissipative version of the KMP model introduced in [25]. In this case, the calculation is straightforward and yields Λ2 ν=sinh(2√2ν)−2√2ν 4√2νsinh2(√2ν).(5.22) Interestingly, Λ2 ν∼1/3 independent of νin the limit of weakly-dissipative systems ν1. This can be understood as a reminiscence of the scaling of G(d) derived in section III. In fact, eq. (3.34) tells us that, for γ= 2, G(d) is just a function of d/dav. In the gaussian approximation, this implies the convergence of Λ2 νto a constant value in the quasi-elastic limit as ν→0+. On the other hand, Λ2 ν∼(2√2ν)−1for ν1, which is consistent with the suppression of dissipation fluctuations previously found in the strongly inelastic regime, as expressed by the general scaling of the LDF given by Eq. (3.48). The same qualitative observations apply to other values of β, though the calculation is more convoluted. We have tested the above predictions in standard Monte Carlo simulations of the dissipative KMP model described in this section, for the particular case β= 0. Figure 3 shows the probability density function (pdf) for the dissipated energy, integrated over the whole system and over a long time τfor many different values of the 10-3 10-2 10-1 100101102103 ν 10-3 10-2 10-1 100 101 102 dav, Λν 2 ~ν1/2 ~ν-1/2 Λν 2~1/3 dav~ν FIG. 4. (Color online) Measured average dissipation and its variance as a function of νfor β= 0. The solid line corresponds to the theoretical prediction for dav, eq. (5.18), while the dashed line is the gaussian estimation of the dissipation variance parameter, Λ2 ν, see eq. (5.22). The agreement is excellent in all cases. Notice in particular the scaling with νof both observables in the weakly- and strongly-dissipative system limits. macroscopic dissipation coefficient ν∈[10−3,103]. In order to minimize finite-size effects in the measurements, we performed simulations for systems with increasing size as νgrows, L∝`−1 ν, in such a way that the number of lattice sites per unit typical length is constant and large enough so we are within the hydrodynamic regime. Furthermore, the integration time τ=O(1) for the continuous, diffusive, timescale over which the hydrodynamic predictions should hold [45]. Standard Monte Carlo simulations do not allow us to sample the tails of the distribution, but they are useful to study the typical fluctuations around the average we are interested in here (e.g., a regime of 5 standard deviations around the average). Figure 3 shows that, when plotted against the reduced variable z≡(d−dav)/χ, the distribution Pτ(z) follows approximately a normal distribution for typical fluctuations. Moreover, all curves for different νcollapse in this regime. However, even at this standard simulation level, it becomes apparent that the tails of the distribution (corresponding to moderate dissipation fluctuations) deviate from gaussian behavior, see inset in Fig. 3, showing asymmetric tails and breaking the collapse to gaussian behavior observed for small fluctuations. As we will show below, the analysis of the dissipation LDF shows that the large fluctuations statistics is far from gaussian. In order to further check our theory, we have also compared the measured average dissipation and its variance with the analytical results above, as a function of the macroscopic dissipation coefficient ν, varying in a range which covers 6 orders of magnitude. Again, we see in Fig. 4 that the agreement is excellent in all cases. In par- 17 ticular, the average dissipation grows as ν(resp. ν1/2) in the weakly (resp. strongly) dissipative system limit, while the variance remains constant for ν1 but decays as ν−1/2for ν1. Remarkably, the gaussian approximation for the variance turns out to be an excellent estimator of the empirical dissipation variance. For Lτ 1, large fluctuations of the dissipation are very rare and most of the probability concentrates in a region of width proportional to (Lτ)−1/2around the average value, a regime described by the gaussian approximation. B. Complete fluctuation spectrum for the integrated dissipation We now investigate the whole spectrum of fluctuations (both typical and rare) of the integrated dissipation. Thus, we need to evaluate the LDF G(d) for arbitrary values of d, in general not close to its average value dav, both analytically and numerically. Exploring in standard simulations the tails of the dissipation distribution associated to the nontrivial structure of G(d) is an daunting task, since LDFs involve by definition exponentially-unlikely rare events, see eq. (2.12). This has been corroborated in Fig. 3, where the dissipation distribution has been measured directly but we are unable to gather enough statistics in the tails of the pdf to obtain clear-cut results in the non-gaussian regime. A recent series of works have addressed this issue, developing an efficient method to measure directly LDFs in many particle systems [27–29]. The method is based on a modification of the dynamics so that the rare events responsible of the large deviation are no longer rare [27], and it has been developed for discrete- [27] and continuous-time Markov dynamics [28]. For a recent review, which also discusses Hamiltonian systems, see ref. [29]. The method yields the Legendre-Fenchel transform of the dissipation LDF, which is usually defined as µ(s) = maxd[G(d) + sd] [10, 46]. In particular, if UC0Cis the transition rate from configuration Cto C0of the associated stochastic process, the modified dynamics is defined as ˜ UC0C(s) = UC0Cexp(s dC0C), where dC0Cis the energy dissipated in the elementary transition C→C0. It can be then shown [11, 27–29] that the natural logarithm of the largest eigenvalue of matrix ˜ U(s) gives µ(s), which in turn can be Legendre-transformed back to obtain a Monte Carlo estimate of G(d). The method of refs. [27–29] thus provides a way to measure µ(s) by evolving a large number Mof copies or clones of the system using the modified dynamics ˜ U(s). This method is exact in the limit M→ ∞, but in practice we are able to simulate a large but finite population of clones, typically M∈[103,104]. This introduces additional finite-size effects related to the population of clones which must be considered with care, see [13] for further discussion along this line. The numerical results for the LDF in the following sections have been obtained using these advanced Monte Carlo techniques. 1. Weakly-dissipative systems, ν1 We now focus our attention on the analysis of LDF of the integrated dissipation for weakly dissipative systems, in which ν1. In the general framework developed in section II, we found a scaling property for G(d), as given by eq. (3.34), dav νγ−2 G(d) = −"1−d/dav Y0(e H)#e H,(5.23) where γ= (2 + β)/(1 + β), Y0(e H) is determined by Eq. (3.29), and the constant e Hdepends only on the ratio d/dav, as given by eq. (3.31). For the sake of concreteness, let us consider now the simplest case β= 0, corresponding to the dissipative KMP model introduced in [25]. Equation (3.26) for the rescaled density profile now reads Y0(x)2= 8 e HY21−Y Y0, Y (±1/2) = 1,(5.24) which can be explicitly integrated, with the solution Y(x, e H) = Y0sech2(xp2e H), Y0= cosh2se H 2,(5.25) where we have already used that Y(x) must be an even function of x. The rescaled current profile Ψ(x) introduced in (3.27) is Ψ(x, e H) = −cosh2qe H 2 p2e H tanh(xp2e H).(5.26) The optimal profiles for the density and the current can be now readily written by combining the previous two equations with Eqs. (3.35)-(3.36), yielding ρ(x) = TY (x) = Tcosh2se H 2sech2(xp2e H),(5.27a) j(x) = davΨ(x) = −dav cosh2qe H 2 p2e H tanh(xp2e H) (5.27b) in terms of e H=e H(d). We have taken into account that y≡ρfor β= 0. Note that the curves ρ(x)/T =Y(x, e H) for different values of νplotted as a function of xonly depend on the relative dissipation d/dav. Now, eq. (3.31) implies that d dav = 2Ψ(−1/2) = sinh p2e H p2e H ,(dav ∼νT ),(5.28) which gives the constant e Himplicitly in terms of d/dav. Finally, particularizing eq. (5.23) for the case γ= 2 we 18 0123456 d/dav -8 -6 -4 -2 0 G(d) ν=0.01 ν=0.1 0123456 d/dav -4 -3 -2 -1 0 G(d) β=1 β=0.5 β=0 β=1 β=0.5 β=0 FIG. 5. (Color online) Scaling of the dissipation LDF in the quasi-elastic limit (ν1) for N= 50, T= 1 and varying β for two different values of ν, namely ν= 0.01 (filled symbols) and ν= 0.1 (open symbols). The solid lines are the MFT predictions in each case. Curves have been shifted vertically for convenience, G(dav) = 0, ∀ν, β. For the case β= 0, the simulation curves are plotted for d<dI, with dIbeing the inflection point at which G(d) changes convexity in the limit ν1 (see the text and also Fig. 6). Inset: Comparison of the theoretical G(d) for different β, where it is clear that increasing βfavors larger dissipation fluctuations. are analyzing (that is, β= 0), we obtain G(d) = p2e Htanh se H 2−e H.(5.29) Note that eq. (5.28) for e H(d) requires some careful analysis. From the general discussion in Section III B, we have that e H>0 for d > dav, and thus p2e His a real number, while e H<0 for d<dav and p2e His imaginary. The latter case poses no problem for G(d), which is always real-valued. In fact, if we write pe H=iq|e H| we arrive at G(d) = −q2|e H|tan q| e H| 2+|e H| for e H<0. In the limit as d→0, we have that e H → −π2/2, and thus G(d)→ −∞ as expected on a physical basis. Equation (5.29) gives a simple scaling form for G(d), independent of ν, for the linear (β= 0) dissipative KMP model in the low-dissipation limit ν1 [25]. As anticipated by (5.23), the curve of G(d) vs the relative dissipation d/dav is independent of νin this quasi-elastic regime. This scaling is fully confirmed in Fig. 5, in which we plot G(d) for different, small values of ν∈[10−2,10−1] measured in simulations of the dissipative KMP model using the advanced Monte Carlo technique described at the beginning of this subsection. In particular, the agreement between theory and simulations is excellent in the broad fluctuation regime that we could measure (see below). The dissipation LDF is highly skewed with a fast decrease for fluctuations d<dav and no negative branch, so fluctu- -100 -80 -60 -40 -20 0 s dav -25 -20 -15 -10 -5 0 5 µ(s) ν=0.01 ν=0.1 0 0.5 1 s dav 0 0.5 1 µ(s) 0 2 4 6 8 10 d/dav -1 -0.5 0 G(n)(d) G G' G'' FIG. 6. (Color online) Scaling plot of the Legendre transform of the dissipation LDF, µ(s) = maxd[G(d) + sd], in the quasielastic limit ν1 for N= 50, T= 1, β= 0 and two different values of ν, namely ν= 0.01 (circles) and ν= 0.1 (triangles). The solid line is the MFT prediction, see eq. (5.34). Notice that µ(s) is defined up to a threshold value sI= 0.878458/dav, beyond which the Legendre-Fenchel transform diverges. The top-right inset shows a zoom around the threshold sI. This is related to the existence of an inflection point in G(d) for dI= 2.27672dav, i.e. a point at which G00(dI) = 0, see middleleft inset, beyond which the dissipation LDF is non-convex, see discussion in main text. ation theorem-type relations linking the probabilities of a given integrated dissipation dand the inverse event −ddo not hold [5, 6]. This was of course expected from the lack of microreversibility, a basic tenet for the fluctuation theorem to apply [36]. The limit e H  1 corresponds to large dissipation fluctuations, where G(d)≈ −1 2[ln(d/dav)]2, that is, a very slow decay which shows that such large fluctuations are far more probable than expected within gaussian statistics (∼ −3 2(d/dav)2). In fact, such slow decay implies the presence of an inflection point in G(d): there is a value dIsuch that G00(dI) = 0. The convexity of G(d) changes at d=dI,G00(d)<0 for d<dI while G00(d)>0 for d > dI. Specifically, Eqs. (5.28) and (5.29) imply that dI/dav = 2.27672 (see middleleft inset in Fig. 6). The complete measurement of non-convex LDFs in computer simulations is a challenge which remains unsolved. The reason is that the advanced Monte Carlo method described above to directly measure large-deviation functions in simulations is based on the Legendre-Fenchel transform for the LDF of interest, which is not well-behaved in regimes where the LDF is non-convex [10]. To better understand this issue, recall that the Legendre-Fenchel transform of the dissipation LDF can be written as µ(s) = max d[G(d) + sd] = G[d∗(s)] + s d∗(s),(5.30) 19 where d∗(s) is solution of the equation ∂G(d) ∂d =−s . (5.31) Note that, mathematically, µ(s) is the Legendre-Fenchel transform of −G(d), because the Legendre-Fenchel transform is defined for convex functions [10]. The partial derivative of Gwith respect to dis related to the first integral of Hamilton equations Πψ0, see Eq. (3.16), which in turn can be obtained from eqs. (3.37) and (5.25), yielding Πψ0=ν∂G ∂d =−e H Tsech2se H 2.(5.32) Equivalently, s=−∂G ∂d =e H νT sech2se H 2.(5.33) In this way, making use of Eqs. (5.28), (5.29) and (5.33), the Legendre transform of the dissipation LDF can be written as µ(s) = 2p2e Htanh se H 2−e H,(5.34) in terms of e H, which is obtained implicitly as a function of sfrom Eq. (5.33). Note that the scaling of G(d) with d/dav, see Eq. (5.23), implies a similar collapse for µ(s) when plotted as a function of s dav. Eq. (5.31) has a single solution d∗(s) for s < 0 and hence poses no problem. On the other hand, due to the existence of an inflection point, G0(d) exhibits a minimum at dI, increasingly smoothly afterward to reach asymptoticaly zero in the limit d→ ∞, see middle-left inset in Fig. 6. Therefore, for s > 0 there exist two solutions d∗ 1(s)≤ dI≤d∗ 2(s) for eq. (5.31), but only the first one maximizes eq. (5.30). This means that we cannot obtain G(d) by inverse Legendre-transforming µ(s) for dissipations above the inflection point dI= 2.27672 dav. In fact, µ(s) is defined up to a critical sI, such that sI= 0.87845/dav (the slope of −G(d) at the inflection point), beyond which µ(s) diverges. This can be seen by noticing the main properties of µ(s), namely ∂µ ∂s =d , ∂2µ ∂s2=−∂2G ∂d2−1 .(5.35) Therefore, µhas a singularity at the value of the slope sI corresponding to the inflection point dI, where ∂2µ/∂s2 diverges. The transition to non-convex behavior thus implies that we can only measure the statistics of rare dissipation fluctuations up to dIusing the cloning algorithm [27–29]. Fig. 6 shows a comparison between the measured µ(s) for two different values of ν1 and the theoretical expectation, up to the critical sI. The agreement is excellent in all cases, and the collapse of µ(s) -0.4 -0.2 0 0.2 0.4 x 0 1 2 3 ρ0(x;d) d/dav=0.30 d/dav=0.37 d/dav=0.46 d/dav=0.54 d/dav=0.78 d/dav=1.01 d/dav=1.20 d/dav=1.46 d/dav=1.73 d/dav=1.97 d/dav=2.22 FIG. 7. (Color online) Top: Optimal energy profiles for varying d/dav and β= 0, measured for ν= 10−3(symbols) and ν= 10−2(dashed lines), and MFT predictions (solid lines). Agreement is very good in all cases. Bottom: MFT prediction for the optimal density profiles for varying d/dav. when plotted against s dav is confirmed. The challenge remains to devise computational techniques capable of exploring rare-event statistics even in regimes where the associated LDF is non-convex. We may solve in a similar way the MFT for the integrated dissipation for arbitrary values of the exponent β, though mathematical expressions are far more convoluted that in the illustrative case β= 0 described above. Fig. 5 also shows the dissipation LDF for other exponents β > 0, as well as the results of numerical experiments in these cases. Qualitatively, the results are equivalent to those discussed above, with a ν-independent scaling form of the LDF in the ν1 limit which goes rapidly to zero as d→0 and has a relatively fat tail for ddav. This tail changes convexity (based on a numerical analysis) at a large dissipation dIwhich increases with (a) νfor fixed β(b) βfor fixed ν. For β=0, we have dI/dav ≃2.8 for ν= 1, while dI/dav >6 for ν= 10. On the other hand, for β= 1, G00(d)<0 in the considered region and no inflection point therein. For the intermediate value, β= 0.5, the positive values of G00(d) are so small that 20 FIG. 8. (Color online) Collapse of the optimal energy profiles measured for N= 50 and T= 1 as a function of the relative dissipation d/dav for different values of ν1, namely ν= 10−2(dotted lines) and ν= 10−1(dashed lines), for β= 0 (top, green), β= 0.5 (middle, red) and β= 1 (bottom, blue). The larger β, the less pronounced the central overshoot is for d > dav. Solid lines correspond to MFT predictions. we have chosen not to eliminate the points behind the numerical inflection point, dI/dav ≃3.2 for ν1 and dI/dav ≃5.1 for ν= 1, although they roughly coincide with the values at which the theoretical and the simulation curves begin to separate. Furthermore, comparison with numerical results is excellent in all cases. Interestingly, see inset in Fig. 5, increasing βresults in a broader dissipation LDF, meaning that large dissipation fluctuations are enhanced as βgrows away from the linear case β= 0. We have also measured the typical energy profile associated to a given dissipation fluctuation for the case β= 0, see top panel in Fig. 7, finding also very good agreement with the macroscopic fluctuating theory developed in this paper. Remarkably, optimal profiles for varying ν1 also collapse for constant d/dav (all the simulations have been done with the same value of the energy density at the boundaries T= 1), as predicted by eq. (5.27a). Furthermore, profiles exhibit the x↔ −x symmetry conjectured in section II B in all cases, with a single extremum which can be minimum or maximum depending on the value of the relative dissipation d/dav, a property which was deduced from the general formalism in section III B. Interestingly, profiles associated to dissipation fluctuations above the average exhibit an energy overshoot in the bulk. This observation suggests that the mechanism responsible for large dissipation fluctuations consists in a continued over-injection of energy from the boundary bath, which is transported to and stored in the bulk before being dissipated. The same qualitative observations and good agreement between theory and simulations is observed for other values of the exponent β > 0, see Fig. 8. Notice in particular the nice collapse of optimal profiles for different values of ν1 0123456 d/dav -12 -10 -8 -6 -4 -2 0 G(d) ν=10-2 ν=10-1 ν=1 ν=10 β=1.0 β=0.5 β=0.0 FIG. 9. (Color online) Dissipation LDF for N= 50, T= 1 and varying β= 0,0.5,1.0 and ν∈[10−2,10]. Curves for β= 0.5 and 0 have been shifted vertically for convenience (recall that G(dav) = 0 ∀ν, β), so that β= 1, 0.5 and 0 correspond to top (blue), medium (red) and bottom (green). The MFT predictions are plotted with lines: solid for β= 1, dashed for β= 0.5, and dotted for β= 0. As in Fig. 5, for a fixed ν increasing βresults in larger dissipation fluctuations. but equal relative dissipation. An interesting observation is that optimal density profiles are less pronounced the larger de nonlinearity exponent βis, see Fig. 8. This gives a plausible explanation of the widening of G(d) as βincreases: for the same value of d/dav and increasing β, the associated optimal profile is closer to the hydrodynamic solution the larger βis, and hence this fluctuation cost decreases, having a larger associated probability. 2. Arbitrary dissipation coefficient ν For arbitrary values of ν&1 no general scaling function can be derived in principle for G(d). For each particular case, the whole variational problem, eqs. (5.11)- (5.16), must be solved, which is often analytically intractable. In order to further advance, we resort now to a numerical evaluation of the optimal profiles, which are used in turn to compute the dissipation LDF. Fig. 9 shows the theoretical predictions for G(d) for increasing, non-perturbative values of ν, together with numerical results from simulations, for different values of β. As for the weakly-dissipative system limit previously discussed, the agreement between theory and measurements in Fig. 9 is quite good. We attribute the observed differences between theory and simulation to finite size effects in the latter, which are more apparent for large νas compared to the weakly-dissipative system limit ν1, compare with Fig. 5, see also [39]. Such strong finite-size effects are expected since the natural length scale associated to a given νis `ν=qˆ D/ν. As follows from Eq. (5.7) and the associated discussion, `νdecreases as νgrows 21 FIG. 10. (Color online) Top: Optimal energy profiles as a function of the relative dissipation measured for ν= 10, N= 50 and T= 1 for the particular case β= 0. Thick (green) lines correspond to measurements while thin (pink) lines are MFT predictions. Bottom: Measured optimal energy profiles for ν= 10, N= 50 and T= 1, and varying values of the nonlinearity exponent β. For a given relative dissipation, energy localization around thermal baths decreases as βincreases. so larger system sizes are needed to observe convergence to the macroscopic limit. In addition, finite-size effects related to the number of clones Mused for the sampling become an issue in this limit [13, 31]. In any case, the sharpening of G(d) as νincreases for any βshows that large dissipation fluctuations are strongly suppressed in this regime, as was argued for ν1 on quite general grounds in Sec. III C). In this strongly-dissipative system limit ν1 the scale `ν→0, and the system decouples effectively into two independent boundary shells. Thus, the energy is concentrated around the boundary baths, a picture which agrees again with the analysis of Sec. III C. This behaviour is evidenced by the optimal energy profiles for a given dmeasured for ν= 10, see Fig. 10, in contrast to the behavior observed for ν1, see Figs. 7-8. The agreement of the observed profiles with MFT predictions is rather good, taking into account the non-negligible finite-size effects affecting these measurements. Bottom panel in Fig. 10 shows the measured energy profiles as a function of the relative dissipation and for different values of the nonlinearity exponent β. From this figure, it is clear that for a given relative dissipation, energy localization around thermal baths decreases as βincreases. This suggests again that, as in the ν1 limit, the probability of a fixed relative dissipation fluctuation d/dav, increases as βgrows, giving rise to a broadening of G(d) with β. VI. SUMMARY AND CONCLUSIONS In this paper we have developed a general theoretical framework for calculating the probability of large deviations for the dissipated energy in a general class of nonlinear driven diffusive systems with dissipation. Our starting point is a mesoscopic fluctuating hydrodynamic theory for the energy density in terms if a few slow hydrodynamic fields, that is, a fluctuating reaction-diffusion equation with a drift term compatible with Fourier’s law and a sink term which can be written in terms of the local energy density. The validity of this hydrodynamic description can be demonstrated for a large family of stochastic microscopic models [39], but it is expected to describe the coarse-grained physics of many real systems sharing the same main ingredients, namely: (i) nolinear diffusive dynamics, (ii) bulk dissipation, and (iii) boundary driving. From this fluctuating hydrodynamic description, and using a standard path integral formulation of the problem, we can write the probability of a path in mesoscopic phase space, that is, the space spanned by the slow hydrodynamic fields. Interestingly, the action associated to this path, from which large-deviation functions for macroscopic observables can be derived, has the same simple form as in non-dissipative systems. This is a consequence of the quasi-elasticity of microscopic dynamics, required in order to have a nontrivial competition between diffusion and dissipation at the mesoscale [39]. We use the derived action functional to investigate the large deviation function of the dissipated energy. The energy dissipated in a non-conserving diffusive system is, together with the energy current, the relevant macroscopic observable characterizing nonequilibrium behavior. A simple and powerful additivity conjecture simplifies the resulting variational problem for the dissipation LDF, from which we arrive at Euler-Lagrange equations for the optimal density and current fields that sustain an arbitrary dissipation fluctuation. A Hamiltonian reformulation of this variational problem greatly simplifies the calculations, allowing us to analyze the general theory in certain interesting limits. A perturbative solution thereof shows that the probability distribution of small (that is, typical) fluctuations of the dissipated energy is always gaussian, as expected from the central limit theorem. Moreover, a general expression for the variance of the distribution in the gaussian approximation has been derived which compares nicely with numerical results. On the other hand, strong separation from the gaussian behavior is expected for large dissipation fluctuations, with a distribution which shows no negative branch, thus violating the Gallavotti-Cohen fluctuation theorem as expected from the irreversibility of the dynamics. Furthermore, the dissipation LDF exhibits simple and general scaling forms in the weakly- and strongly-dissipative system limits, which can be analyzed in general without knowing the explicit solution of the canonical equations. We apply our results to a general class of diffusive lattice models for which dissipation, nonlinear diffusion 22 and driving are the key ingredients. The theoretical predictions, which can be explicitely worked out in certain cases, are compared to extensive numerical simulations of the microscopic models (which cover both typical fluctuations and rare events), and excellent agreement is found in all cases. In particular, the simple scaling for the dissipation large-deviation function in the weakly-dissipative system limit is fully confirmed for different values of the nonlinearity exponent β, exhibiting non-convex behavior for large enough fluctuations. Interestingly, in this limit ν1 energy profiles associated to large dissipation fluctuations exhibit an overshoot in the bulk resulting from an excess energy injection from boundary baths. On the other hand, in the strongly-dissipative system limit ν1 the typical lengthscale goes to zero and the system decouples into two almost-independent boundary shells, giving rise to a different scaling form for the LDF and a strong suppression of the dissipation fluctuations in this regime. Recently, a similar hydrodynamic theory has been developed to study large fluctuations in a particular class of driven dissipative media [37], but its predictions do not compare well with numerical results of the dissipative lattice models here studied. The reason is that Ref. [37] studies systems with two competing dynamics, one conservative and another nonconservative, thus resulting in independent fluctuations for the density and dissipation fields. In our theory, as is the case in many driven dissipative systems, dissipation is linked to the collision process, and hence dissipation fluctuations are enslaved to density profile deviations. In fact, both theories coincide in the limit where the optimal dissipation profile is given in terms of the optimal density field. In summary, our results show that a suitable generalization of macroscopic fluctuation theory [7] is capable of describing in detail the fluctuating behavior of general nonlinear driven dissipative media. In this scheme, the dissipation LDF follows from a variational problem whose solution also gives the optimal profiles that the system has to sustain to achieve the considered fluctuation. The proposed framework is very general, as MFT is based only on (a) the knowledge of the conservation laws governing a system, which allow to write down the balance equations for the fluctuating fields, and (b) a few transport coefficients appearing in these fluctuating balance equations. This opens the door to further general results in the nonequilibrium statistical physics of dissipative media. In particular, it would be interesting to explore the existence of phase transitions and spontaneous symmetry breaking at the fluctuating level therein, in a way similar to the phenomenon reported in conservative systems [30, 31]. Moreover, as the relevant magnitudes characterizing nonequilibrium behavior in dissipative systems are both the dissipated energy and the current, it would be worth analyzing the joint fluctuations of these two observables within the MFT approach. ACKNOWLEDGMENTS We acknowledge financial support from Spanish Ministerio de Ciencia e Innovaci´on projects FIS2011-24460 and FIS2009-08451, EU-FEDER funds, and Junta de Andaluc´ıa projects P07-FQM02725 and P09-FQM4682. Appendix: Variational problem with a Lagrangian including second-order derivatives Let us analyze a variational problem in which the “action” is defined as the integral of a “Lagrangian” with second order derivatives, that is S[j] = Zx2 x1 dx L(j, j0, j00).(A.1) The action S[j] is a functional of the profile j(x) in the fixed interval x1≤x≤x2. The variational problem arises when one looks for the “optimal” profile j(x) for which the functional S[j] is a extremum. For the sake of concreteness, let us consider a problem similar to the one analyzed in this paper: we are interested in calculating Gdefined as G=−min j(x)S[j].(A.2) Then, we consider the variation δSof the functional when a given profile j(x) is slightly changed to j(x) + δj(x), δS=Zx2 x1 dx ∂L ∂j δj +∂L ∂j0δj0+∂L ∂j00 δj00.(A.3) Now, we take into account that δj0=d dxδj, δj00 =d2 dx2δj (A.4) in order to integrate by parts (i) once the term proportional to δj0(ii) twice the term proportional to δj00. We arrive thus at δS=∂L ∂j0−d dx ∂L ∂j00 δj +∂L ∂j00 δj0x2 x1 +Zx2 x1 dx ∂L ∂j −d dx ∂L ∂j0+d2 dx2∂L ∂j00 ,(A.5) 23 where [f]x2 x1=f(x2)−f(x1). By analogy with the case of the usual Lagrangian with only first-order derivatives, we introduce the generalized momenta as pj=∂L ∂j0−d dx ∂L ∂j00 , pj0=∂L ∂j00 .(A.6) In this way, the boundary term has the usual form and Eq. (A.5) can be rewritten as δS= [pjδj +pj0δj0]x2 x1+Zx2 x1 dx ∂L ∂j −d dx ∂L ∂j0+d2 dx2∂L ∂j00 δj. (A.7) The extremum condition is δS = 0. If the values of j and j0are prescribed at the boundaries, both δj and δj0 vanish at x1,2and the boundary term vanishes. Then, as δj is arbitrary for x1< x < x2, the “optimal” profile solution of the variational problem verifies the Euler- Lagrange equation d2 dx2∂L ∂j00 −d dx ∂L ∂j0+∂L ∂j = 0,(A.8) which is a fourth-order differential equation. Interestingly, Eq. (A.8) can be written as dpj/dx =∂L/∂j that is formally identical to the usual Euler-Lagrange equation for Lagrangians with only first-order derivatives. The boundary conditions for the Euler-Lagrange equation are the prescribed values of jand j0at the boundaries; four conditions for the fourth-order differential equation. However, in physical problems there are sometimes less prescribed quantities at the boundaries than necessary. In that case, as pointed out by Lanczos [33], the extremum condition δS= 0 provides the “missing” boundary conditions. For instance, if we only have fixed values of j0at the boundaries (as in the LDF problem we have dealt with in the main text), δj0(x1) = δj0(x2) = 0 but δj(x1) and δj(x2) are free parameters. Equation (A.7) still implies the Euler-Lagrange equation but also that pj(x1) = pj(x2)=0.(A.9) The generalized momentum conjugate of the variable that is not fixed at the boundary must vanish: the solution of the variational problem verifies then the Euler- Lagrange equation (A.8) with the prescribed values of j0 at the boundaries and the “extra” conditions provided by Eq. (A.9). In this way, we obtain the four conditions needed to determine completely the solution of the Euler-Lagrange equation. The function Gdefined in Eq. (A.2) depends on the values of jand j0at the boundaries. Making use of Eq. (A.7), and taking into account that the optimal profile j(x) verifies the Euler-Lagrange equation, we get δG =−pj,2δj2−pj0,2δj0 2+pj,1δj1+pj0,1δj0 1.(A.10) We have introduced the notation pj,i =pj(xi), δji= δj(xi), i= 1,2 and so on. Equation (A.10) implies that pj,2=−∂G ∂j2 , pj0,2=−∂G ∂j0 2 , pj,1=∂G ∂j1 , pj0,1=∂G ∂j0 1 . (A.11) Equation (2.33) of the main paper is the particularization of this result for the case (i) x2=−x1= 1/2, (ii) solutions of the Euler-Lagrange equation with well-defined parity, in which pj,2=pj,1,pj0,2=−pj0,1, and (iii) the boundary conditions of Eq. (2.30). [1] L. D. Landau and E. M. Lifshitz, Statistical Physics 3rd edition, Course of Theoretical Physics Vol. 5 (Pergamon Press, Oxford, 1980). [2] H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951). [3] R. Kubo, Reports on Progress in Physics 29, 255 (1966). [4] R.S. Ellis, Entropy, Large Deviations and Statistical Mechanics (Springer, New York, 1985). [5] G. Gallavotti and E.G.D. Cohen, Phys. Rev. Lett. 74, 2694 (1995). [6] J.L. Lebowitz and H. Spohn, J. Stat. Phys. 95, 333 (1999). [7] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim, Phys. Rev. Lett. 87, 040601 (2001); Phys. Rev. Lett. 94, 030601 (2005); J. Stat. Mech. P07014 (2007); J. Stat. Phys. 135, 857 (2009). [8] T. Bodineau and B. Derrida, Phys. Rev. Lett. 92, 180601 (2004) [9] B. Derrida, J. Stat. Mech. P07023 (2007). [10] H. Touchette, Phys. Rep. 478,1 (2009). [11] P.I. Hurtado and P.L. Garrido, Phys. Rev. Lett. 102, 250601 (2009); Phys. Rev. E 81, 041102 (2010). [12] C. P´erez-Espigares, J.J. del Pozo, P.L. Garrido and P.I. Hurtado, AIP Conf. Proc. 1332, 204 (2011). [13] P.I. Hurtado and P.L. Garrido, J. Stat. Mech. P02032 (2009). [14] P.I. Hurtado, C. P´erez-Espigares, J.J. del Pozo and P.L. Garrido, Proc. Natl. Acad. Sci. USA 108, 7704 (2011). [15] P. L. Krapivsky and B. Meerson, Phys. Rev. E 86, 031106 (2012). 24 [16] J. J. Brey, M. I. Garc´ıa de Soria, P. Maynar and M. J. Ruiz-Montero, Phys. Rev. Lett. 94, 098001 (2005). [17] J.-F. Millithaler, E. Alfinito and L. Reggiani, IEEE Conference Proceedings, vol. 21st International Conference on Noise and Fluctuations (ICNF), 417 (2011), doi: 10.1109/ICNF.2011.5994358 [18] J. Murray, Mathematical Biology, Springer Verlag (2002). [19] S. T. Bramwell, P. C. W. Holdsworth and J.-F. Pinton, Nature 396, 552 (1998). [20] S Ramaswamy, Annu. Rev. Condens. Matt. Phys. 1323 (2010); M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao and R. Aditi Simha, arXiv:1207.2929v1. [21] I. Epstein and J. Pojman, An Introduction to Nonlinear Chemical Dynamics, Oxford University Press, Oxford (1998); A. De Wit, Adv. Chem. Phys. 109, 435 (1999); P. Gaspard, J. Chem. Phys. 120, 8898 (2004). [22] R.S. Cantrell and C. Corner, Spatial Ecology via Reaction-Diffusion Equations, Wiley Series (2003). [23] C. Kipnis, C. Marchioro and E. Presutti, J. Stat. Phys. 27, 65 (1982). [24] D. Levanony and D. Levine, Phys. Rev. E 73, 055102(R) (2006). [25] A. Prados, A. Lasanta, and P.I. Hurtado, Phys. Rev. Lett. 107, 140601 (2011). [26] T. P¨oschel and S. Luding eds., Granular Gases, Lecture Notes in Physics vol.564 (Springer-Verlag, Berlin, 2001). [27] C. Giardin`a, J. Kurchan and L. Peliti, Phys. Rev. Lett. 96, 120603 (2006). [28] V. Lecomte and J. Tailleur, J. Stat. Mech. (2007) P03004; J. Tailleur and V. Lecomte, AIP Conf. Proc. 1091, 212 (2009). [29] C. Giardina, J. Kurchan, V. Lecomte, J. Tailleur, J. Stat. Phys. 145, 787 (2011). [30] T. Bodineau and B. Derrida, Phys. Rev. E 72, 066110 (2005). [31] P.I. Hurtado and P.L. Garrido, Phys. Rev. Lett. 107, 180601 (2011). [32] I. M. Gelfand and S. V. Fomin, Calculus of Variations (Dover, New York, 2000). [33] C. Lanczos, The Variational Principles of Mechanics (Dover, New York, 1986). [34] G.Eyink, J. L. Lebowitz, and H. Spohn, Commun. Math. Phys. 132, 253 (1990); 140, 119 (1991). [35] J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. E 62, 5339 (2000). [36] A. Puglisi, P. Visco, A. Barrat, E. Trizac, and F. van Wijland, Phys. Rev. Lett. 95, 110202 (2005). [37] T. Bodineau and M. Lagouge, J. Stat. Phys. 139, 201 (2010) [38] P.I. Hurtado and P.L. Krapivsky, Phys. Rev. E 85, 060103(R) (2012). [39] A. Prados, A. Lasanta and P. I. Hurtado, Phys. Rev. E 86, 031134 (2012). [40] A mathematically similar situation was found in Ref. [15] for the fluctuations of the current in conservative systems, where a “standard” variational problem with a Lagrangian with only first-order derivatives of the current j arose. Thus, in the canonical description, there was only one canonical momentum pj, and the particular solution of the canonical equations with pj= 0 also corresponded to the average behaviour. [41] Multiple maxima (minima), if present, must be separated by intermediate minima (maxima), but this is not possible as y(x) takes the same value y0at all these extrema. The possibility of an infinite (continuum) set of extrema is ruled out because of the smoothness of Q(y) for the typical systems of interest. [42] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (New York, Springer, 1999). [43] J. P. Boon, J. F. Lutsko, and C. Lutsko, Phys. Rev. E 85, 021126 (2012). [44] Γ(z) is the gamma function, Γ(z) = R∞ 0dt tz−1e−t, with the property Γ(z+ 1) = zΓ(z) [45] In the simulations, time is measured in Monte Carlo steps. A Monte Carlo step comprises Ncollision events of the microscopic dynamics introduced in Sec. IV. On the other hand, the unit of our continuous timescale t comprises N2Monte Carlo steps [39]. [46] The parameter sappearing in the Legendre transform of G(d) is conjugated to the dissipation via the equation s= −∂G/∂d, and turns out to be just minus the Lagrange multiplier introduced in Sec. II A, s=−λ.