Multivariate fluctuation relations for currents

This paper is devoted to multivariate fluctuation relations for all the currents flowing across an open system in contact with several reservoirs at different temperatures and chemical potentials, or driven by time-independent external mechanical forces. After some transient behavior, the open system is supposed to reach a nonequilibrium steady state that is controlled by the thermodynamic and mechanical forces, called the affinities. The time-reversal symmetry of the underlying Hamiltonian dynamics implies symmetry relations among the statistical properties of the fluctuating currents, depending on the values of the affinities. These multivariate fluctuation relations are not only compatible with the second law of thermodynamics, but they also imply remarkable relations between the linear or nonlinear response coefficients and the cumulants of the fluctuating currents. These relations include the Onsager and Casimir reciprocity relations, as well as their generalizations beyond linear response. Methods to deduce multivariate fluctuation relations are presented for classical, stochastic and quantum systems. In this way, multivariate fluctuation relations are obtained for energy or particle transport in the effusion of an ideal gas, heat transport in Hamiltonian systems coupled by Langevin stochastic forces to heat reservoirs, driven Brownian motion of an electrically charged particle subjected to an external magnetic field, and quantum electron transport in multi-terminal mesoscopic circuits where the link to the scattering approach is established.


Introduction
If driving systems out of equilibrium requires free energy supply, controlling their dynamics is based on the distribution of this supply into separate channels and the transduction of free energy between its different forms: chemical, electrical, mechanical or else. During recent decades, this preoccupation has shifted from the macroscale down to micro-and nanometric systems such as Brownian particles driven in optical traps [1][2][3], self-propelled Janus particles [4,5], enzymes and biomolecular motors [6][7][8][9][10][11][12], molecular machines [13,14], catalytic particles [15] or semiconducting electronic circuits [16][17][18][19]. These small systems are in contact with heat, chemical or electrical reservoirs. They may be driven away from equilibrium by external mechanical forces, or thermodynamic forces given by differences in temperature or chemical potentials between several reservoirs in contact with the small system. The reservoirs are supposed to be large enough for their thermodynamic conditions to remain constant in time. All these possible driving parameters are called the affinities and they induce heat or particle currents flowing across the system [20][21][22][23][24]. In the case of chemical reactions, the currents are defined by the fluxes of molecules transforming from reactants into products under the drive of the chemical affinities.
For nonequilibrium systems close to equilibrium, the currents are proportional to the affinities by the linear response coefficients, which obey the Onsager-Casimir reciprocity relations as the consequence of microreversibility, as well as the Green-Kubo formulae [24][25][26][27]. Since the 1950s, such results have been extended to the nonlinear response properties [28][29][30][31][32][33]. Remarkable relationships have been discovered by Bochkov and Kuzovlev for closed and open systems ruled by time-reversal symmetric Hamiltonian dynamics and driven out of equilibrium by mechanical forces [30][31][32]. More recently, close connections with nonequilibrium thermodynamics and the second law have been established thanks to the so-called fluctuation relations [34][35][36][37][38][39][40][41][42][43][44][45][46][47]. Such relations were first obtained for a unique fluctuating variable globally characterizing the nonequilibrium driving of the system. In this context, an important issue is that systems are driven out of equilibrium not only by mechanical forces, but also by chemical or thermal affinities, which are of statistical origin. These latter affinities, which are very common in natural phenomena, are induced by coupling the system of interest with several reservoirs at different chemical potentials or temperatures. Indeed, a system in contact with a single reservoir could not be driven out of equilibrium without mechanical forces. In the case of molecular motors, the chemical affinities are due to differences of free energies between reactants and products mixed together at specific concentrations in the solution surrounding the motor. The description of these systems should combine very different formulations to take into account, on the one hand, the mechanical forces and, on the other hand, the chemical affinities of statistical origin. Nevertheless, the challenge can be overcome and multivariate fluctuation relations have been proved for nonequilibrium systems driven by chemical or thermal affinities and performing free energy transduction [48][49][50][51][52][53][54][55][56].
In this direction, many questions remain open about systems driven in nonequilibrium steady states by combined mechanical, chemical, electrical or thermal affinities, especially, if they are moreover subjected to an external magnetic field. Such systems are the stage of thermoelectric, galvanomagnetic or thermomagnetic effects, which have been investigated mainly in linear regimes for lack of general methods beyond. Furthermore, fundamental questions arise concerning the origins of thermodynamic time asymmetry in such systems, because it is still unclear at which level the time-reversal symmetry is broken in the description of nonequilibrium processes.
The purpose of the present paper is to address these issues from the unifying viewpoint provided by the multivariate fluctuation relations and their consequences exposed in section 2. For this purpose, methods are developed in order to establish multivariate fluctuation relations in different kinds of classical, stochastic and quantum systems. A key issue is to define the affinities of nonmechanical origins. In this regard, several reservoirs should be coupled to the system and the nonmechanical affinities can thus be defined with respect to a reference reservoir.
For the effusion process of an ideal gas through a small hole in a thin wall separating two reservoirs, a multivariate fluctuation relation is here deduced in section 3 directly from the underlying Hamiltonian dynamics. This deduction shows that such relations hold because of time-reversal symmetry breaking at the statistical level of description, but without requiring the approximations at the basis of Markovian master equations.
Nevertheless, many natural processes such as driven Brownian motion or Hamiltonian systems in contact with heat baths by Langevin forces are very well described in terms of stochastic processes. For such systems, multivariate fluctuation relations are here established thanks to the method of the time-evolution operator that is modified to include parameters counting heat and particle transports [40,57]. In section 4, the modification is carried out by using the path-integral theory of Machlup and Onsager [58] and a multivariate fluctuation relation is proved for all the heat currents flowing between the heat baths across an arbitrary Hamiltonian system. In section 5, this method is extended and a multivariate fluctuation relation is obtained for driven Brownian motion in an external magnetic field, allowing us to consider galvanomagnetic effects in nonlinear regimes.
Remarkably, multivariate fluctuation relations can also be established in the same unifying framework for open quantum systems in nonequilibrium steady states and subjected to an external magnetic field. In the case of quantum electron transport in multi-terminal circuits, it is here shown that the two-time measurement set up [59,60] is closely related to Landauer's scattering approach [61][62][63]. In particular, the cumulant generating function of the full counting statistics obtained in the two-time measurement set up is here shown to be equivalently given by the Levitov-Lesovik formula [64] also in the presence of an external magnetic field. For these quantum systems, the multivariate fluctuation relation is directly deduced from the time-reversal symmetry of the scattering matrix.
In this way, different kinds of open systems driven in nonequilibrium steady states by either mechanical and nonmechanical affinities are here shown to obey multivariate fluctuation relations as the consequence of microreversibility. Therefore, the Onsager-Casimir reciprocity relations and Green-Kubo formulae can be generalized from the linear to the nonlinear response properties, as explained section 2. In the presence of an external magnetic field, these properties concern the galvanomagnetic and thermomagnetic effects in nonlinear regimes. Conclusions are drawn in section 7.

The multivariate fluctuation relation
An open system in contact with several reservoirs at different temperatures and chemical potentials is crossed by heat and particle currents (see figure 1). The system is at equilibrium if the temperature and the chemical potentials are uniform in all the reservoirs. Otherwise, the global system is out of equilibrium and a nonequilibrium steady state may establish itself after some transient behavior. At the macroscale, the steady state is characterized by the average values of the currents across the open system. These currents depend on the thermodynamic forces or affinities A including the thermal affinities: the chemical affinities: between the reservoirs j = 1, 2, . . . , r − 1 at the temperatures {T j } and the chemical potentials {µ j p } for the different particle species p = 1, 2, . . . , s, and the reference reservoir j = r Figure 1. Schematic representation of (a) an open system S in contact with four reservoirs; (b) two interacting open systems S 1 and S 2 , each in contact with two reservoirs. If s particle species are transported between every reservoir besides energy, the system (a) is controlled by d = 3s + 3 affinities and the system (b) by d = 2s + 3 affinities. [20][21][22][23][24]. k B is Boltzmann's constant. The nonequilibrium conditions are thus fixed by d = (s + 1)(r − 1) different affinities in a globally coupled system such as the one in figure 1(a).
If the system is composed of c channels where the s particle species are separately conserved and if the channels interact with each other so that the total energy is globally conserved as in figure 1(b), there are c different reference chemical potentials for each one of the s species but a unique reference temperature, so that the nonequilibrium steady state of this system is controlled by d = s(r − c) + r − 1 affinities. We notice that the thermodynamic forces (1) and (2) may have to be supplemented by external mechanical forces if the system is subjected to them. The system reaches the steady state of thermodynamic equilibrium if all these affinities are vanishing. At the mesoscale, the currents are fluctuating and their statistical properties are described by a stationary probability distribution if the system is in a steady state. The random variables of interest are the quantities of matter and energy Q = { Q i } that are transferred between the reservoirs during a time interval [0, t]. There are as many such quantities as there are affinities. The fluctuating currents are thus defined as J = Q/t. If P A denotes the probability distribution describing the steady state corresponding to the affinities A, the multivariate fluctuation relation for the currents reads which compares the opposite fluctuations of the currents. This result has been proved in particular for Markovian stochastic jump processes [51,55] using the Hill-Schnakenberg cycle decomposition of the graph associated with the Markovian process [6,65]. It has been extended to semi-Markov processes in particular in [53]. At equilibrium where A = 0, we recover the principle of detailed balancing, according to which opposite fluctuations are equiprobable. This principle results from the fact that equilibrium states are described by Gibbsian probability distributions that are functions of the Hamiltonian. Since the Hamiltonian is invariant under time reversal (in the absence of external magnetic field), the Gibbsian probability distributions are also invariant. Accordingly, every phase-space trajectory over a time interval [0, t] has the same probability as its reversal, hence the principle of detailed balancing.
Out of equilibrium where A = 0, the time-reversal symmetry is broken at the statistical level of description in terms of the probability distribution P A and the probabilities of opposite fluctuations are thus no longer equal. Nevertheless, the multivariate fluctuation relation (3) holds thanks to the microreversibility of the underlying microscopic dynamics. Since the ratio of probabilities is increasing or decreasing exponentially with time, one of both fluctuations soon dominates over the opposite and a directionality appears in the open system. This directionality is controlled by the affinities and can thus be reversed by changing the sign of the affinities.
Remarkably, the relation (3) shows that the current fluctuations are larger in the direction of the affinities, which thus control the directionality. Therefore, the knowledge of the affinities in the frame defined by the reservoirs determines the direction of the most probable current fluctuations [66,67]. In this regard, the question arises whether it is possible to guess the direction of the affinities A imposed by the reservoirs from the sole observation of a current fluctuation J in some macroscopic steady state. Using Bayesian inference [45,68], the likelihood of the hypotheses that the observed current fluctuation J goes either forward (+) or backward (−) with respect to the direction A is given by Tossing a fair coin on either hypothesis amounts to taking P(+) = P(−) = 1 2 , while P(J|±) = P A (±J). According to the multivariate fluctuation relation (3), the likelihood that the observed current fluctuation J indeed goes forward with respect to the direction of the affinities A is thus given by where θ(x) is Heaviside's function. In the limit of a long positive time, the likelihood reaches the unit probability that A · J 0, if the system is out of equilibrium with A = 0. In contrast, if the affinities vanish A = 0, the likelihood remains at the value P(+|J) 1 2 , which confirms the absence of directionality at equilibrium. The reasoning is compatible with the overall timereversal symmetry because a similar result holds in the limit t → −∞ for A · J replaced by −A · J. As aforementioned, the question is pertinent only if the reservoirs and their affinities remain unknown to the observer.

The multivariate fluctuation relation and the second law of thermodynamics
The multivariate fluctuation relation (3) implies that the thermodynamic entropy production is always non-negative. In consistency with nonequilibrium thermodynamics at the macroscale, the entropy production should be identified with the sum of the affinities multiplied by the average values of the currents in units of Boltzmann's constant [24]. In a nonequilibrium steady state, these average values are obtained in the long-time limit as if the process is ergodic with respect to the probability distribution P A . According to the multivariate fluctuation relation (3), the entropy production is given by in terms of the Kullback-Leibler divergence between the distributions P A (J) and P A (−J), which is always non-negative, hence the agreement with the second law of thermodynamics. In this regard, the fluctuation relation shows that the thermodynamic entropy production is a measure of the breaking of the time-reversal symmetry by the nonequilibrium probability distribution P A (J). A further remark is that the inverse of the entropy production in units of Boltzmann's constant gives an estimation of the characteristic time taken by the likelihood (5) to approach asymptotically the unit probability: This characteristic time becomes tiny if the system is macroscopic and driven far from equilibrium by increasing the affinities A.
Since free energy dissipation generates entropy production, its non-negativity can be used to define thermodynamic efficiencies for different possible processes transducing free energy into a specific form. For instance, in order to drive a particular current J k A in the direction opposite to its associated affinity A k , energy should be supplied by using the other currents. In this case, A k J k A < 0 and a thermodynamic efficiency can be defined according to which is always less or equal to the unit value according to the second law of thermodynamics (7) as a consequence of the multivariate fluctuation relation (3).

Consequences for the response properties
An alternative expression of the fluctuation relation (3) can be obtained in terms of the generating function of the statistical cumulants defined as where λ are the so-called counting parameters. The full counting statistics of the currents is known to be completely characterized by this function. The theory of large deviations provides the methods to determine the long-time behavior of the probability distribution P A (J) from the knowledge of the cumulant generating function [69]. An important remark is that this function may be piecewise defined with respect to the counting parameters λ ∈ R d , if the probability distribution P A (J) presents exponential tails for some specific currents. In the following, the cumulant generating function is supposed to be indefinitely differentiable at λ = 0, as well as at The average values of the currents, their diffusivities, as well as their higher cumulants are given by taking the successive derivatives of the generating function (9) with respect to the counting parameters: . . .
On the other hand, the average value of some current can be expanded in powers of the affinities as which defines the linear and nonlinear response coefficients: . . .
Remarkably, these properties are interrelated thanks to the multivariate fluctuation relation (3), as the following reasoning shows. Inserting the fluctuation relation (3) in the definition (9) of the cumulant generating function yields the symmetry relation Taking successive derivatives of this relation with respect to the counting parameters and the affinities, the cumulants (11)-(13), etc and the response coefficients (15)- (17), etc are found to be interrelated. Indeed, taking the derivative of equation (18) with respect to the counting parameter λ i at λ = 0 gives two equivalent expressions for the corresponding average current: Expanding the second expression in powers of the counting parameters, which are here set equal to the affinities λ = A, the average current can be expressed only in terms of the cumulants (11)-(13), etc as which shows that the response coefficients are completely determined by the cumulants characterizing the fluctuations around the nonequilibrium or equilibrium steady state [31,32]. Now, the response coefficients are obtained by expanding the cumulants in powers of the affinities. Accordingly, the symmetry relation (18) implies the Green-Kubo formulae and the Onsager reciprocity relations for the linear response coefficients [31,36] With higher derivatives, generalizations of these relations to the higher cumulants and the nonlinear response coefficients can be deduced for systems driven away from equilibrium by mechanical forces [31,32] and chemical or thermal affinities [48,50,52]. In particular, the nonlinear response coefficients (16) and (17) are related to the second and third cumulants according to . . .
which are generalizations of the Green-Kubo formulae. A generalization of Onsager reciprocity relations is given by the total symmetry of the following fourth-order tensor: which is obtained by considering all the fourth derivatives of the generating function (9) with respect to the affinities and the counting parameters [48,50,52]. We notice that the thirdcumulant tensor (12) is vanishing at equilibrium, C i jk (0) = 0. Similar relations exist at higher orders as well [31,52]. They are the consequences of the underlying microreversibility. As discussed in sections 5 and 6, similar results generalizing the Casimir-Onsager reciprocity relations can also be obtained for systems in an external magnetic field.
The average value of this random variable is equal to the entropy production in units of Boltzmann's constant. Its probability distribution is given by where δ(x) denotes the Dirac delta distribution.
Using the multivariate fluctuation relation (3), we deduce the univariate relation for ζ : The random variable ζ is related to the quantity called the action functional in [40].

Univariate fluctuation relation for a sole current.
In general, the multivariate fluctuation relation (3) holds for all the currents flowing across an open system and it does not imply the validity of similar relations for a subset of currents. Let us consider a system with two currents, in which case the fluctuation relation (3) is bivariate and reads We may wonder [70] if there exist specific conditions under which the univariate fluctuation relation is satisfied for the marginal distribution of the current J 1 defined as and for an effective affinityÃ 1 that may differ from the actual affinity A 1 . It turns out that such univariate fluctuation relations may indeed hold, but under the following specific conditions.

Tight coupling between the currents.
The tight-coupling condition is defined by requiring that both currents remain proportional to each other during their random time evolution: which is only possible under special circumstances encountered for instance in molecular motors [8]. Inserting the condition (32) in the bivariate fluctuation relation (29), we get the univariate fluctuation relation (30) with the effective affinitỹ which is associated with the coupled currents. In this case, the entropy production (7) reduces to so that the sole affinityÃ 1 drives the system out of equilibrium. Under the tight-coupling condition (32), the thermodynamic efficiency is directly determined by the affinities The condition of tight coupling between mechanics and chemistry is encountered in the rotary motor F 1 -ATPase if it is subjected to a small enough external torque. Otherwise, the torque decouples the mechanical and chemical currents and the motor enters in the regime of loose coupling [7,8].

Separation of time scales.
In other circumstances, the transition rates for one current may be much higher than for the other current, e.g. | J 2 | | J 1 |, in such a way that the assumption holds for some A 1 (A 1 , A 2 ) in the long-time limit. Inserting the bivariate fluctuation relation (29) in the assumption (36), we get the univariate fluctuation relation (30) for the marginal (31) and the effective affinityÃ 1 , which depends on both affinities A 1 and A 2 .
Using Jensen's inequality e X e X here with X =Ã 1 J 1 t − (A 1 J 1 + A 2 J 2 ) t and noting that e X 1 as the consequence of the assumption (36), we find that the entropy production has the following lower bound: in terms of the effective affinityÃ 1 of the univariate fluctuation relation (30) [71,72]. Accordingly, the thermodynamic efficiency is here limited to a value lower than unity: This situation is encountered in mesoscopic electronic devices where the current J 1 in quantum dots is driven by the Coulomb drag of a large current J 2 in a quantum point contact capacitively coupled to the quantum dots [71,72]. We notice that the bound (38) is attained for the tightcoupling condition, which is thus stronger. These results show that the multivariate fluctuation relation is very useful to characterize the coupling between different currents and the efficiencies of free energy transduction.

Connection to other types of large-deviation relationships
Generally speaking, the multivariate fluctuation relation (3) characterizes the large-deviation properties of the current fluctuations. Further quantities of interest can be observed during the time evolution of the system. This is the case for Brownian particles driven by an optical trap where the position of the Brownian particle can be recorded with a certain sampling time, giving the statistics of the random histories followed by the system [2,3]. Such observations provide much more details on the time evolution of the system than the measurement of the fluctuating currents. From the viewpoint that every observable corresponds to some degree of coarse graining, the fluctuating currents are coarser than the positions of the particles sampled at a finite frequency. The finer the graining, the higher the rate of data accumulation in the recording of the time evolution of the system. The quantitative theory of data accumulation rates has been developed, in particular, for the study of chaotic dynamics [73]. In this context, large-deviation relationships have also been discovered relating chaotic and transport properties [74][75][76][77].
Considering microreversibility, the data accumulation rate of the typical histories can be compared with the data accumulation rates of their reversals. Out of equilibrium, both rates differ because the time-reversal symmetry is broken by the probability distribution P A of nonequilibrium steady states. A remarkable result is that the difference between both is the Kullback-Leibler divergence between the probabilities of the histories and their reversals, which is again related to the thermodynamic entropy production, but at a finer degree of coarse graining than in equation (7) [78,79]. The dependence on the resolution used for coarse graining has been studied in the case of driven Brownian particles and noisy electric RC circuits [2,3].
In the following, we shall see how multivariate fluctuation relations can be established for different classes of systems.

Multivariate fluctuation relation for heat and particle currents in effusion
In this section, our purpose is to show that the multivariate fluctuation relation can be directly deduced from the microscopic Hamiltonian dynamics without using a stochastic master equation. This demonstration emphasizes the importance of the breaking of timereversal symmetry by the nonequilibrium probability weights given to the trajectories that are the solutions of Hamilton's equations, although the equations themselves remain perfectly symmetric. Furthermore, the demonstration shows that there is no need to use the Markovian stochastic master equation to deduce the fluctuation relation [82]. The multivariate fluctuation relation is thus more fundamental than previously thought.

The invariant probability measure for an ideal gas
With this aim, let us consider the effusion of an ideal gas through a small hole in a thin wall separating the gas between two reservoirs at different temperatures and densities (see figure 2). The ideal gas is composed of independent monoatomic particles of mass m moving in free flight and undergoing elastic collision on the thin wall W . The free flights are the solutions of Hamilton's equations for the single-particle position r ∈ R 3 and momentum p ∈ R 3 with the Hamiltonian function given by H = p 2 /(2m). Moreover, the free flights are eventually interrupted by elastic collisions on the thin wall W reflecting the z-component of the momentum according to the collision rule of specular reflection on the wall W . The motion of every particle is thus ruled by a Hamiltonian flow t in the single-particle phase space }. This single-particle Hamiltonian flow is symplectic, obeys Liouville's theorem, and is symmetric under the time-reversal transformation (r, p) = (r, −p): • t • = −t . In the left-hand reservoir, the particles arrive from z = −∞ with velocities distributed according to a Maxwell-Boltzmann distribution at the temperature T L and the density n L . In the right-hand reservoir, the particles arrive from z = +∞ with velocities distributed according to a Maxwell-Boltzmann distribution at the temperature T R and the density n R . The single-particle distribution function has thus the Maxwell-Boltzmann form where the particle density n C and the temperature T C take the values associated with the trajectory C to which the phase-space point (r, p) belongs and corresponding to the domain from which the trajectory is coming. There are three types of possible trajectories: (i) the The key point is that the trajectories {C 2 } crossing the hole move to a reservoir at a different temperature and density than where they come from. Therefore, the single-particle distribution function (39) may take different values for these trajectories than for their corresponding reversals { C 2 }. Accordingly, we have that The distribution function is invariant under the Hamiltonian flow, , p), if the reservoirs have different temperatures or densities: T L = T R or n L = n R . For monoatomic particles, the chemical potentials are related to the densities by Because of the infinite spatial extension of the reservoirs, the ideal gas contains an infinite number of particles and its microstate is given by the point In this infinite-particle system, the Hamiltonian flow Γ t = Φ t (Γ 0 ) is also time-reversal symmetric: where Θ is the time-reversal transformation acting on the microstate (42) by reversing all the momenta: p i → −p i for i = 1, 2, 3, . . . , N , . . .. For the infinite-particle system, an invariant probability measure µ can be constructed as a Poisson suspension on the basis of the invariant distribution function (39) [80]. For any sixdimensional domain D ⊂ M of the single-particle phase space M, we consider the random events for which there are N particles in the domain D. The average number of particles in the phasespace domain D is given by dr dp (45) in terms of the single-particle distribution function (39). The number N of particles in the phasespace domain D is a random variable of Poisson distribution Moreover, random events in which disjoint phase-space domains have given particle numbers are statistically independent so that Both equations (46) and (47) define the probability distribution of the so-called Poisson suspension [80]. This probability measure is invariant under the time evolution of the Hamiltonian flow Φ t : for any random event A. This is the consequence of the stationarity of the single-particle distribution function (39), f [ t (r, p)] = f (r, p), which implies the invariance ν t D = ν (D) of the measure (45). The flow Φ t in the infinite phase space M and the invariant probability measure of the Poisson suspension defines the infinite-particle dynamical system ( t , M, µ), which is known to be ergodic and mixing [80]. Although the flow Φ t has the microreversibility (43), the dynamical system is not symmetric under time reversal because as already pointed out in [75]. The time-reversal symmetry is thus broken at the statistical level of description in terms of the stationary probability distribution µ under nonequilibrium conditions, which is essential for the multivariate fluctuation relation, as shown here below.

Exact multivariate fluctuation relation for the effusion of an ideal gas
In the previous framework, the counting statistics of the particle and energy currents can be directly established. The method is similar as for the sole particle current [81], but more involved since the groups of particles with different energies should be separated. The transfers of energy and particles from the left-hand to the right-hand reservoirs during the time interval [0, t] are given by where dN ε+ (respectively, dN ε− ) denotes the number of particles with a positive (respectively, negative) z-component v z of their velocity v = p/m, a kinetic energy such that H = 1 2 mv 2 ∈ [ε, ε + dε], and going through the hole in the wall during the time interval [0, t]. In the phase space, these particles belong to the domain D ε+ . In position space, this domain corresponds to the cylinder of basis given by the hole of area σ , of height equal to v z t, and of volume σ v z t (see figure 2). The average number of these particles is given by where we used equation (39) for an orbit C = L, on which the condition (40) holds. Similarly, the average number of particles of kinetic energy H = 1 2 mv 2 ∈ [ε, ε + dε] flowing from the righthand to the left-hand reservoir during the time interval [0, t] is given by The ratio of these average numbers satisfies the relation in terms of the thermal and chemical affinities here given by Now, the cumulant generating function for the transport of energy and particles between the reservoir is defined as in terms of the affinities A = (A E , A N ), the counting parameters λ = (λ E , λ N ), and the transfers (50) and (51). The axis of energy 0 ε < ∞ is divided into small intervals [ε, ε + dε] of width dε. In every one of these intervals, the numbers of particles transferred from the lefthand to the right-hand reservoir or vice versa constitute Poisson random variables, n ε+ and n ε− , of average values given by equations (52) and (53), respectively. Therefore, the moment generating function can be calculated as follows: in terms of the multiple Poisson distribution where ε and denote the sum and product over all the intervals [ε, ε + dε] partitioning the energy axis. Performing the summation over the random numbers {n ε+ , n ε− }, we get Using equations (52) and (53) and taking the limit dε → 0, the generating function (57) becomes Integrating over the energy ε, the cumulant generating function for heat and particle transport in effusion is finally obtained as as it should [82,83]. Transforming equation (62) with the ratio (54), the multivariate fluctuation relation is established for effusion in the form: in terms of the affinities (55) and (56). Accordingly, the consequences presented in the previous section, in particular concerning the response properties, apply to effusion [83]. We notice that the demonstration has been directly carried out from the microscopic Hamiltonian dynamics and the invariant probability distribution defining the infinite-particle dynamical system without using any stochastic master equation. The demonstration clearly shows that the fluctuation relation characterizes time-reversal symmetry breaking by the invariant probability distribution under nonequilibrium conditions, while microreversibility is always satisfied by the underlying Hamiltonian dynamics.
If the wall between both reservoirs had a thickness and the hole a billiard-like geometry with a possible chaotic repeller, the scattering of particles by the wall should be studied to determine which trajectories are going through the hole and which are back-scattered, and the result would be similar.
In the case where the gas is not ideal but diluted or rarefied, the effusion process can be described in terms of the fluctuating Boltzmann equation, for which a multivariate fluctuation relation has also been established [84].

The system and its time evolution
Heat transport is a fascinating topic that has recently focused much interest in connection with fluctuation relations [85][86][87][88][89][90][91][92][93][94][95][96][97]. Here, our purpose is to show that multivariate fluctuation relations hold for a Hamiltonian system coupled to heat reservoirs, as schematically depicted in figure 3. We suppose that the Hamiltonian system contains N particles of masses m i , positions r i ∈ R 3 , and momenta p i ∈ R 3 with i = 1, 2, . . . , N . The total energy of the isolated system is given by the Hamiltonian function Moreover, the particles j = 1, 2, . . . , r are coupled to so many heat reservoirs of temperatures {T j } r j=1 as Brownian particles of friction coefficients {ζ j } r j=1 . Accordingly, Newton's equations ruling their motion include a friction force proportional to their velocity, −ζ jṙ j , and a Langevin random force f j (t) given by a corresponding Gaussian white noise characterized by for j, j = 1, 2, . . . , r , where 1 is the 3 × 3 unit matrix. Newton's equations thus read as for j = 1, 2, . . . , r , and for i = j. If ρ = ρ(R, P) denotes the probability density in the phase space (R, P) = {r i , p i } N i=1 ∈ R 6N of the system, the Fokker-Planck equation ruling the Langevin stochastic process (67) and (68) is given by in terms of the Poisson bracket and the operatorŝ

The evolution operator with the counting parameters
The Fokker-Planck operator (69) is modified to include the counting parameters λ = {λ j } r −1 j=1 of the heat transfers from r − 1 heat reservoirs to the Hamiltonian system and the reference reservoir j = r . For this purpose, we use the theory by Machlup and Onsager [58]. The probability that the Gaussian white noises of the r Langevin forces (66) would take specific values during the time interval [0, t] gives the contribution of all the r reservoirs to the probability of the path R(t) = {r i (t)} N i=1 followed by the system: We notice that the path probability also includes the Liouvillian contribution from the Hamiltonian dynamics, which is deterministic and preserves phase-space volumes. The probability weight given in equation (72) corresponds to the operators (71) in the Fokker-Planck equation (69). Now, we consider the moment generating function of the heat currents in a steady state where Q = { Q j } r −1 j=1 denote the quantities of heat transferred from the reservoirs to the Hamiltonian system during the time interval [0, t].
Using Stratonovich's integral [86], the heat transfers can be expressed in terms of the instantaneous heat currents as In order to make the usual correspondence with the partial differential operator, Itô's integral is required [98], which is given by Combining with the path probability (72), we get Provided that the path probability (72) corresponds to the Fokker-Planck operator (69), the modified path probability (76) corresponds to the operator whereL 0 is the Fokker-Planck operator (69). The sum has been restricted to r − 1 terms because the reference reservoir j = r does not need a counting parameter since heat transfers are measured with respect to it. Moreover, the identity ∂ ∂p j · (p j ρ) = p j · ∂ρ ∂p j + 3ρ has been used, which explains that the last term comes with a plus sign in place of the minus sign in the last term of equation (76).
Using equation (69), the modified operator (77) can be written in the form whereK j,λ j denotes the operator depending on the counting parameter λ j in equation (77). For the following, the key point is that the leading eigenvalue of the modified operator gives the cumulant generating function aŝ Indeed, the cumulant generating function is the long-time decay rate of the moment generating function (73), but the time evolution of this latter is ruled by the modified operator, whereupon its leading eigenvalue gives the cumulant generating function. Although the modified operator is not Hermitian, the Perron-Frobenius theorem guarantees that its adjoint has the same leading eigenvalue [99].

The time-reversal symmetry and the multivariate fluctuation relation
Since the modified operator and its adjoint share the same leading eigenvalue, the fluctuation relation (18) results from the symmetry relation where (R, P) = (R, −P) is the time-reversal transformation and β r = (k B T r ) −1 is the inverse temperature of the reference reservoir. This symmetry relation is proved as follows. Firstly, the Hamiltonian part of the modified operator (78) satisfies the symmetry since the Liouvillian operator {H, · } is anti-Hermitian. The equation (81) because the Hamiltonian function depends on the momentum p r as H = p 2 r 2m r + · · ·. Thirdly, we verify that with the thermal affinity for every j = 1, 2, . . . , r − 1. Consequently, the symmetry relation (80) holds for the modified operator and its leading eigenvalue has the symmetry (18), which proves the multivariate fluctuation relation for heat transport.

Consequences on heat conduction properties
Microreversibility and the fluctuation relation (18) have fundamental consequences on the heat conduction properties of mesoscopic devices [100]. Thermal rectifiers made of carbon and boron nitride nanotubes [101] or quantum dots [102] are studied in the new field of phononics [103]. Rectification in thermal diodes is based on the nonlinear response of these devices to the thermal affinity A = (k B T 2 ) −1 − (k B T 1 ) −1 . Thanks to the fluctuation relation (18), the nonlinear as well as the linear response coefficients are given in terms of the fluctuation properties of the heat current in nonequilibrium steady states. For a thermal diode, there is a single heat current J A . Accordingly, the second and third cumulants (11) and (12) can be expanded in powers of the single affinity A as where D n = d n D dA n (0) and C n = d n C dA n (0). The coefficients D 0 and C 0 = 0 characterize the fluctuations at equilibrium. Inserting these expansions in equation (20), the heat current is given by so that all the response coefficients turn out to be determined by the statistical cumulants.
Rectification is due to the terms that are even powers of the affinity. The first of them has the coefficient D 1 , which characterizes the sensitivity of the diffusivity of the current fluctuations with respect to some perturbation A away from equilibrium. This coefficient gives the main contribution to the rectification ratio Surprisingly, thermal rectification appears as a feature of heat current fluctuations under nonequilibrium conditions.

Equation of motion and Fokker-Planck equation
In this section, we consider an electrically charged particle of mass m undergoing Brownian motion in a spatially periodic energy potential V (r), in a uniform external magnetic field B, and driven out of equilibrium by a uniform external electric field E. This particle is for instance a charged Brownian particle moving in an optical lattice or a charged ion in a crystal. The potential has thus the spatial periodicity of a crystallographic group: V (r) = V (r + a) for all the vectors a ∈ L of the lattice L. The position r = (x, y, z) of this particle obeys a stochastic differential equation of Langevin type: where e is the electric charge of the particle, ζ its friction coefficient and f(t) a Langevin random force such that with i, j = x, y, z and the temperature T of the medium where the particle moves. This model has been formulated for the study of orbital diamagnetism [104]. More recently, it has been used to investigate Jarzynski's identity [105] and Crooks' fluctuation relation [38,39] for a charged Brownian particle in a time-dependent potential or magnetic field [106,107]. The issue of geometric magnetism has also been studied with such a model [108]. Here, we are interested by the multivariate fluctuation relation in nonequilibrium steady states. The potential V (r) is thus time independent, as well as the external electric and magnetic fields.
In the absence of external electric field E = 0, the particle is at equilibrium in the medium at the temperature T . If the electric field is switched on E = 0, the particle undergoes a biased Brownian motion with a mean drift velocity to be determined. The affinity controlling the nonequilibrium drive of this system is defined as The motion of the particle admits a Hamiltonian description in terms of the momentum and the Hamiltonian function Under the time-reversal transformation (r, p) = (r, −p), the external magnetic field should also be reversed to get the symmetry of the Hamiltonian function In these terms, the equation of motion (89) becomes If ρ(r, p) denotes the probability density in phase space, the Fokker-Planck equation corresponding to the stochastic differential equation reads which defines the Fokker-Planck operatorL A,B .

The evolution operator with the counting parameters
Under nonequilibrium conditions, the drift velocity of the Brownian particle represents the current of this transport process. Integrating the velocityṙ over a time interval [0, t], we find that Jt = r t = r t − r 0 . Therefore, the cumulant generating function is defined as in terms of the counting parameters λ = (λ x , λ y , λ z ).
Since the integrated current is given in terms of the basic observable that is the position r, the modification of the evolution operator to include the counting parameter is here given byL in terms of the Fokker-Planck operator (96). Since the Poisson bracket {·, ·} involves a derivative with respect to the position, the modified operator readŝ for an arbitrary differentiable function . In order to emphasize the fact that the Brownian particle is driven out of equilibrium by the external electric field E, it is interesting to write the modified operator in the form

The time-reversal symmetry and the multivariate fluctuation relation
In the presence of an external magnetic field B = 0, but zero electric field E = 0, the Fokker-Planck operator (96) obeys the following time-reversal symmetry: with the inverse temperature β = (k B T ) −1 . This identity can be verified by direct calculation. We notice that the magnetic field is reversed in the left-hand side of this identity, but it is not the case for adjoint operator in the right-hand side. Now, if the affinity (91) is switched on A = 0, the modified operator satisfies We notice that the Hamiltonian function in the Boltzmann factor continues to be taken with zero electric field. The reason is that the Boltzmann factor defines the reference equilibrium state so that the electric field should be equal to zero otherwise the Brownian particle would be driven out of equilibrium. The symmetry (102) is established by focusing on the operator for which we verify that In the first line, we have used the fact that time reversal changes the sign of the momentum and the magnetic field. At the second line, the chain rule of derivatives is used. At the third line, the point is that ∂ p H E,B does not depend on the electric field E and the definition (91) of the affinity is used. The fourth line is finally obtained because ∂ † p = −∂ p . This proves the symmetry (102) of the modified operator.
Again, the cumulant generating function is obtained as the leading eigenvalue of the modified operator or its adjoint aŝ given the fact that the generating function is real. Using the symmetry (102), we find that the eigenfunctions of the operator and its adjoint are related bỹ while the cumulant generating function obeys the symmetry relation which corresponds to the following multivariate fluctuation relation: as can be checked using the definition (97) of the cumulant generating function.

Consequences on galvanomagnetic properties
The multivariate fluctuation relation implies a hierarchy of time-reversal symmetry relations between the response coefficients and the cumulants of the random drift of the Brownian particle, here, in the presence of the external magnetic field. The average drift velocity is in general a nonlinear function of the affinities, which vanishes at equilibrium. The response coefficients are introduced by expanding the average drift velocity in powers of the affinities as with i, j, k = x, y, z. Taking second derivatives of the symmetry relation (108), the Casimir-Onsager reciprocity relations are recovered We notice that the Hall effect is identified by decomposing the linear response tensor (111) into its symmetric and antisymmetric parts as with The symmetric part is even in the magnetic field, while the antisymmetric part is odd. In three-dimensional space, the antisymmetric tensor can be written in terms of an axial vector characterizing the Hall effect. Moreover, the antisymmetric part behaves linearly in weak magnetic field: [109]. Remarkably, the symmetry relation (108) also implies similar relations valid beyond linear response. In particular, the following relations hold for the nonlinear response coefficients at second order in the perturbations with respect to equilibrium [60]: If i = j = k, we infer from equation (116) that unidirectional response coefficients are given by The other coefficients remain linked together by equations (115)- (116). Taking the third derivatives of the symmetry relation (108), we can also deduce that, at equilibrium where A = 0, the third cumulants (12) are odd with respect to the magnetic field C i jk (0; B) = −C i jk (0; −B). (118) These cumulants characterize the magnetic-field asymmetry of the fluctuations [110,111]. In the absence of magnetic field B = 0, this magnetic-field asymmetry disappears because the third cumulants vanish at equilibrium, C i jk (0; 0) = 0. Accordingly, these cumulants are proportional to the external magnetic field, as it was the case for the antisymmetric part of the linear response tensor (114): Furthermore, the third cumulants are related to the responses of the diffusivities (11) with respect to the affinities according to Taking the difference between equation (116) for ±B and using equation (120), we can infer that the third cumulants are also given by which could be verified by direct experimental measurements. Generalizations have also been deduced at higher orders [60].
These results show that the galvanomagnetic nonlinear response properties obey as fundamental relationships as the Casimir-Onsager reciprocity relations and the Green-Kubo formulae are for linear response.

Multivariate fluctuation relation for quantum electron transport in multi-terminal circuits
In this section, we turn to quantum transport in an open system in contact with several reservoirs, as depicted in figure 1. Advances in microelectronics have led to the fabrication of submicrometric semiconducting circuits where electron transport is essentially ballistic. These circuits are composed of several elements such as quantum dots and quantum point contacts. The coherence of electron transport between the reservoirs can be described thanks to the scattering approach developed by Landauer, Büttiker, Imry and others [61][62][63].
If early interest was focused on the quantization of the conductance itself, recent efforts have been devoted to the noise properties and to the full counting statistics [62]. A formula has been obtained by Levitov and Lesovik for the cumulant generating function of full counting statistics within the scattering approach [64]. More recently, the connection has been established with the fluctuation theorem and the issue of time-reversal symmetry [44,46,50,59,60,112].
Thanks to the formalism of quantum mechanics, multivariate fluctuation relations have already been established for all the energy and particle currents flowing across an open quantum system. Here, our purpose is to show the equivalence between the theory developed in [60] and the Levitov-Lesovik formula obtained in the scattering approach [64]. Furthermore, the time-reversal symmetry of the generating function is shown to result directly from the known symmetry of the scattering matrix in the presence of an external magnetic field [61].

Quantum theory of full counting statistics
The counting statistics of the energy and particle currents is established in a scheme with an initial and a final quantum measurements. The system evolves in time under the effect of the interaction coupling the open system to the reservoirs during a time interval [0, t]. Before the initial time t = 0 when the interaction is switched on, the energies and particle numbers can be measured in the reservoirs during the semi-infinite lapse of time −∞ < t < 0. After the time interval [0, t], the interaction is switched off and the energies and particle numbers of the reservoirs can again be measured during the semi-infinite lapse of time t < t < +∞. The differences between the final and the initial values of the energies and particle numbers give the corresponding transfers, which fluctuate over the many measurements performed to obtain the statistics.
During the time interval [0, t], the Hamiltonian operator is given bŷ whereĤ j is the Hamiltonian of the jth reservoir andV the interaction. In the presence of an external magnetic field B, the Hamiltonian operator (122) has the symmetrŷ under the time-reversal transformationˆ . The quantum states evolve in time under the unitary evolution operator as The evolution operator has the symmetrŷ under time reversal. The reservoirs contain different species of particles p = 1, 2, . . . , s. The particle number operators {N j p } are observables that commute among themselves and with the reservoir Hamiltonian operators {Ĥ j }. The total particle numberŝ with p = 1, 2, . . . , s are constants of motion together with the total Hamiltonian (122). Here, we assume that the borders between the reservoirs are located in the middle of the system itself [113].
Before the interaction is switched on, every reservoir is supposed to be in a grand-canonical equilibrium state. The reservoirs may have different temperatures {T j } and chemical potentials {µ j p }. The initial density matrix of the total system is thus given bŷ where j is the grand-canonical partition function of the jth reservoir and β j = (k B T j ) −1 its inverse temperature.
The initial quantum measurement is performed before the time t = 0 on the energies and particle numbers of the decoupled reservoirs: The outcome of this measurement is the quantum state | k , which appears with the probability A subsequent quantum measurement is performed after the time evolution (124) during the time interval [0, t] and after the interaction between the reservoirs is switched off: Its outcome is the other quantum state | l , which appears after the initial state | k and thus with the probability Accordingly, the energies and particle numbers that are transferred between the reservoirs during the time interval [0, t] are defined as The probability density to observe these transfers is therefore given by P t ( ε j , n j p ; B) = kl j p δ ε j − ε jl + ε jk δ n j p − n j pl + n j pk The moment generating function is defined as where {ξ j , η j p } are the counting parameters. Using the expression (136), the moment generating function becomes whereĤ This set up with a preparation period −∞ < t < 0 and a post-interaction measurement for t < t < +∞ is reminiscent of the scattering approach [114,115]. In this approach, the interaction evolves the total wavefunction from the incoming to the outgoing state. The former is asymptotic towards the past and the latter towards the future. This similarity allows us to prove here below the equivalence between the counting statistics obtained in the two-measurement set up and in the scattering approach for independent particles.

Connection to the scattering approach for independent electrons
We consider the transport of independent electrons in a circuit as the one of figure 1(a). For independent electrons, the Hamiltonian and other operators are quadratic functions of the annihilation and creation operators {ĉ q ,ĉ † q } of an electron in every one-electron state {φ q } of the system:X whereX denotes a many-body operator andx = (x qq ) the corresponding one-body operator.
Klich has shown that the trace of products of exponential functions of many-body operators can be expressed as appropriate determinants involving the one-body operators [116,117]. For fermions, Klich's formula reads tr eX eŶ = tr e (x) e (ŷ) = det 1 + ex eŷ .
Now, the moment generating function (138) has the form which commute, [X ,Ŷ ] = 0. According to Klich's formula (142), the moment generating function becomes Using the property that the determinant of a product of operators is equal to the product of determinants, the generating function can be rewritten as with the operator describing the Fermi-Dirac distributions in the reservoirs. In the following, we consider a single species of fermions s = 1.
The connection to the scattering approach can be established by introducing the unitary scattering operator defined in the long-time limit bŷ S = lim t→∞ e iĥ 0 t/2 e −iĥt e iĥ 0 t/2 (149) in terms of the full Hamiltonianĥ and the noninteracting Hamiltonianĥ 0 = j h j and with h = 1 [114,115]. Therefore, we have that e −ŷ t = e iĥt e −ŷ e −iĥt e iĥ 0 t/2Ŝ † e −ŷŜ e −iĥ 0 t/2 (150) over a long enough time interval. Since the noninteracting Hamiltonianĥ 0 = j h j commutes with the operator (148), the moment generating function can be written as To make sense, the determinant should be taken over an appropriate discrete set of oneelectron states forming a quasi continuum in the long-time limit. Indeed, the left-hand side of equation (151) concerns the many-body system although the expression in the right-hand side concerns single independent electrons flowing one by one across the system. Every one-body operator can be decomposed on the eigenstates of the noninteracting Hamiltonianĥ 0 . This is the case for the scattering operator: whereŜ(ε) is the scattering matrix. This is a r × r unitary matrix if there is a single open channel in every reservoir, which is the simplifying assumption we use in the following. Similarly, the operator (148) decomposes aŝ in terms of the matrix with the Fermi-Dirac distributions of the reservoirs: Moreover, the counting parameters appear in the matrix whereξ = (ξ j δ j j ) andη = (η j δ j j ) are the diagonal r × r matrices containing the counting parameters on their diagonal. Accordingly, the moment generating function (151) becomes a product over all the relevant single-electron states {ε, σ } as where σ is the spin label. Using the aforementioned property of determinants, the unitarity of scattering matrix, and the commutativity of the diagonal Hermitian operatorsf (ε) and e εξ +η , the expression (157) can be shown to be real, as it should. A crucial property is that the r × r determinants in equation (157) are invariant under the transformations e εξ +η → e εξ +η−χ(ε)1 (158) for any function χ (ε) of the energy ε and where1 is the identity matrix. In [60], the function χ(ε) = 1 r r j=1 (ε ξ j + η j ) was used, which corresponds to counting the currents with respect to a reference taken by averaging over the r reservoirs. In the following, we consider the function χ(ε) = ε ξ r + η r , which corresponds to using the r th reservoir as the reference. Consequently, the generating function (157) only depends on the differences of the counting parameters with respect to those of the reference reservoir, as it should to properly describe nonequilibrium steady states. We can thus define the 2(r − 1) counting parameters: with j = 1, 2, . . . , r − 1, corresponding to the r − 1 independent energy and electron currents in the circuit. These currents are induced by the 2(r − 1) thermal and chemical affinities given by for j = 1, 2, . . . , r − 1.
The cumulant generating function is introduced as where A = {A j E , A j N } r −1 j=1 and λ = {λ j E , λ j N } r −1 j=1 . Taking the logarithm of the function (157), the product over the relevant states {ε, σ } becomes a sum of logarithms. For a process lasting over a time interval [0, t], the spacing between the relevant energies {ε} is equal to ε = 2πh/t. Since these relevant energies form a quasi continuum, the sum is replaced by an integral in the infinite-time limit. Moreover, the spin orientation takes the two values σ = ± for electrons. For spin one-half particles such as electrons, we thus find that with the factor 2 s due to the electron spin [62]. The cumulant generating function is finally obtained as for independent electrons. This result is equivalent to the Levitov-Lesovik formula [64], as can be checked for two-terminal circuits [118]. In equation (165), the counting parameters λ only appear in the operators e εξ +η . The scattering matrixŜ(ε) depends on the external magnetic field B. The Fermi-Dirac distributionsf (ε) depend on the temperatures and the chemical potentials of the r reservoirs or, equivalently, on the r − 1 affinities A together with the temperature T r and the chemical potential µ r of the reference reservoir. The Fermi-Dirac distribution may also depend on the magnitude B of the external magnetic field if this latter extends to the reservoirs, but not on its sign if there is no spontaneous magnetization in the reservoirs.

The time-reversal symmetry and the multivariate fluctuation relation
For spinless particles or for each spin component in systems without spin-orbit interaction, the time-reversal symmetry implies that the scattering matrix in the presence of an external magnetic field B obeyŝ where T denotes the transpose [61,63]. In order to deduce the multivariate fluctuation relation in the form (18), we use the property the fact that the determinants in equation (157) are real, and the symmetry (166). We thus have the following identities: After using the cyclic property (167), equation (168) results from the fact that the determinant is real and equation (169) from the time-reversal symmetry (166) of the scattering matrix. Finally, equation (170) is the consequence of the identitŷ whereβ = (β j δ j j ) andμ = (µ j δ j j ) are the diagonal r × r matrices containing the inverse temperatures and the chemical potentials on their diagonal. Therefore, equation (170)

Consequences on the galvanomagnetic and thermomagnetic properties
The multivariate fluctuation relation (172), which finds its origin in the fundamental symmetry (123), has consequences on the nonlinear galvanomagnetic and thermomagnetic properties of mesoscopic devices performing free energy transduction between multiple electron and heat reservoirs [100]. These consequences are the fundamental relations (115)- (121), which generalize the Casimir-Onsager reciprocity relations (111) from the linear to the nonlinear response and fluctuation properties. Here, the indices refer to the different independent affinities and currents between the reservoirs (instead of the three spatial dimensions in section 5). Till now, such relations have been experimentally investigated for the single-current case in a circuit with an Aharonov-Bohm ring in a external magnetic field at low temperature [17,18]. However, equations (115) and (116) and equations (118)-(121) constitute predictions for the coefficients characterizing the coupling between the multiple electrical and thermal currents flowing across mesoscopic devices. These fundamental results concern not only electronic devices such as quantum dots at low temperature, but also transistors and other devices sensitive to an external magnetic field at room temperature. In such devices, noise manifests itself as in RC circuits and experimental methods are available to measure independently the response coefficients and the statistical cumulants in order to verify the predictions of the multivariate fluctuation relation [119].

Conclusions
In this paper, multivariate fluctuation relations are demonstrated for all the currents flowing across different kinds of open systems. A unified presentation is given for classical, stochastic and quantum systems. They are maintained in nonequilibrium steady states by thermodynamic forces called affinities, which are defined in terms of the difference of temperatures and chemical potentials between several reservoirs, or by mechanical forces. The multiple affinities controlling the nonequilibrium steady state induce so many currents. They may be coupled together, allowing free energy transduction in the open system. Multivariate fluctuation relations find their origin in the microreversibility of the underlying Hamiltonian dynamics. Yet, the time-reversal symmetry is broken by the probability distribution describing the nonequilibrium steady state, which is reminiscent of symmetry-breaking phenomena at equilibrium [120,121]. The non-negativity of the entropy production can be directly deduced from the multivariate fluctuation relation, which is thus compatible with the second law of thermodynamics. The multivariate fluctuation relation implies the Onsager or Casimir reciprocity relations between the linear response coefficients, as well as the Green-Kubo formulae between these coefficients and the second cumulants of the fluctuating currents, i.e. the current diffusivities. Most remarkably, the multivariate fluctuation relation also implies generalizations of these relations between the nonlinear response coefficients, higherorder coefficients and their responses to perturbations with respect to equilibrium. Under certain conditions presented in section 2.4, the multivariate fluctuation relation may reduce to univariate fluctuation relations.
Different methods are here presented to obtain the multivariate fluctuation relation for specific open systems.
For the effusion of an ideal gas through a small hole in a thin wall separating two reservoirs at different temperatures and chemical potentials, the multivariate fluctuation relation can be directly deduced from the Hamiltonian dynamics of the independent particles, using the classical scattering approach [75]. The invariant probability distribution describing the nonequilibrium steady state is given by a many-particle Poisson distribution based on the Maxwell-Boltzmann single-particle distribution. The cumulant generating function of the currents can be calculated from this invariant probability distribution and its symmetry under time reversal can be established, depending on the affinities of the reservoirs.
For stochastic processes ruled by a Fokker-Planck master equation, the multivariate fluctuation relation for the currents can be obtained from a time-reversal symmetry of the Fokker-Planck evolution operator after its modification to include the counting parameters associated with the energy or particle currents between the reservoirs. Accordingly, the leading eigenvalue of this modified evolution operator gives the cumulant generating function of the currents. The symmetry of the modified evolution operator implies the symmetry of the cumulant generating function, hence the multivariate fluctuation relation. The advantage of this method is that the symmetry can be established by analyzing separately the different terms composing the modified evolution operator and corresponding to different processes contributing to the overall time evolution [40,57,84].
This method is applied to heat transport in Hamiltonian systems coupled to several heat reservoirs by Langevin stochastic forces. As shown in section 4, the evolution operator can be modified by using the Machlup-Onsager path integral defining the generating function of the heat currents. A fundamental symmetry is obtained for the so-modified stochastic evolution operator, which established the multivariate fluctuation relation. This result has fundamental consequences on the nonlinear properties of heat conduction in thermal diodes or other mesoscopic devices. In particular, a formula is deduced showing that thermal rectification is determined at leading order by the sensitivity to nonequilibrium perturbations of the diffusivity characterizing the heat current fluctuations.
In section 5, a charged Brownian particle driven in a spatially periodic potential by an external electric field and subjected to a magnetic field is considered with a similar method. Here, microreversibility should be conceived together with the reversal of the external magnetic field. Consequently, the multivariate fluctuation relation involves two invariant probability distributions differing by the sign of the magnetic field. In this respect, the Casimir-Onsager reciprocity relations are recovered for the linear response coefficients and the multivariate fluctuation relation also yields their generalizations to higher orders, which allows us to make predictions on the nonlinear galvanomagnetic properties of this system. For open quantum systems, the multivariate fluctuation relation can be deduced by considering the transfers of energy and particles during the time interval starting and ending by a quantum measurement of the energy and particle numbers in the reservoirs. In section 6, an equivalence is established between the two-measurement set up and the quantum scattering approach. For independent electrons moving in multi-terminal mesoscopic circuits, the Levitov-Lesovik formula [64] is shown to be equivalent to the formula obtained in [60]. Moreover, the multivariate fluctuation relation is directly proved from the time-reversal symmetry of the scattering matrix in the presence of an external magnetic field. In this way, the relations between the nonlinear response coefficients and the statistical cumulants of the current fluctuations are shown to hold also for the nonlinear galvanomagnetic and thermomagnetic properties.
The unified presentation shows the generality of the multivariate fluctuation relations in classical, stochastic and quantum open systems. These relations allow us to understand much better than previously the subtle interplay between the microreversibility of the underlying Hamiltonian dynamics of the particles composing matter and the thermodynamic time asymmetry appearing at the statistical level of description. The perspectives given by the multivariate fluctuation relations extend from physics to chemistry and biology. In chemistry, multivariate fluctuation relations apply to nonequilibrium reactions in gases, liquids or at surfaces [122][123][124]. They are fundamental to our understanding of autonomous molecular machines, self-propelled Janus particles or catalysts, where chemical energy powers mechanical motion or pattern oscillations [4,5,15]. In biology, multivariate fluctuation relations provide a framework to understand free energy transduction in molecular motors and their thermodynamic efficiencies in regimes of loose or tight coupling [7,8,49]. More fundamentally, these new advances shed a new light on the thermodynamic origins of directionality and dynamical order in natural systems.