Glauber's Ising chain between two thermostats

TL;DR

This paper analyzes the energy current in a one-dimensional Ising model with spins connected to two thermostats at different temperatures, using fermionization to compute cumulants and explore various physical regimes.

Contribution

It provides an exact calculation of the cumulant generating function for the energy current in the Glauber dynamics of the Ising chain with two thermostats, revealing new phenomena and scaling behaviors.

Findings

Explicit formulas for cumulants up to fourth order.

Identification of regimes with kinetic mean-field effects and zero-temperature dissipation.

Linear growth of cumulants with system size in the zero-temperature limit.

Abstract

We consider a one-dimensional Ising model each of whose $N$ spins is in contact with two thermostats of distinct temperatures $T_1$ and $T_2$. Under Glauber dynamics the stationary state happens to coincide with the equilibrium state at an effective intermediate temperature $T(T_1,T_2)$. The system nevertheless carries a nontrivial energy current between the thermostats. By means of the fermionization technique, for a chain initially in equilibrium at an arbitrary temperature $T_0$ we calculate the Fourier transform of the probability $P({\cal Q};\tau)$ for the time-integrated energy current ${\cal Q}$ during a finite time interval $\tau$. In the long time limit we determine the corresponding generating function for the cumulants per site and unit of time $\langle{\cal Q}^n\rangle_{\rm c}/(N\tau)$ and explicitly give those with $n=1,2,3,4.$ We exhibit various phenomena in specific regimes: kinetic mean-field effects when one thermostat flips any spin less often than the other one, as well as dissipation towards a thermostat at zero temperature. Moreover, when the system size $N$ goes to infinity while the effective temperature $T$ vanishes, the cumulants of ${\cal Q}$ per unit of time grow linearly with $N$ and are equal to those of a random walk process. In two adequate scaling regimes involving $T$ and $N$ we exhibit the dependence of the first correction upon the ratio of the spin-spin correlation length $\xi(T)$ and the size $N$.

Equations

$$\mathsf{A}={\bar{\nu }}_1{\bar{\nu }}_2(1-{\gamma }_1{\gamma }_2),\mathsf{B}={\bar{\nu }}_1{\bar{\nu }}_2({\gamma }_2-{\gamma }_1),$$

$$H(s)=-E\sum _{n=1}^N{s}_n{s}_{n+1},$$

$$|P(\tau )⟩=\sum _{s}P(s;\tau )|s⟩.$$

$$w_n(s;{\beta }_1)=\frac{1}{2}{\bar{\nu }}_1[1-\frac{1}{2}{\gamma }_1{s}_n(s_{n-1}+s_{n+1})],$$

$$\frac{\text{mbox}d}{\text{mbox}d\tau }|P(\tau )⟩={\bar{\nu }}_1{M}_{\text{th}}({\beta }_1)|P(\tau )⟩$$

$${M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})=\tfrac{1}{2}\sum_{n=1}^{N}(\sigma_{n}^{x}-1)\big{[}1-\tfrac{1}{2}\gamma_{1}\,\sigma_{n}^{z}(\sigma_{n-1}^{z}+\sigma_{n+1}^{z})\big{]},$$

$$|P_{\mathrm{eq}}({\beta }_1)⟩={\rho }_{\mathrm{eq}}({\beta }_1)|1⟩,|1⟩\equiv \sum _s|s⟩,$$

$${\rho }_{\mathrm{eq}}({\beta }_1)=\frac{\text{mbox}{e}^{-{\beta }_1\mathcal{H}}}{Z({\beta }_1)},\mathcal{H}=-E\sum _{n=1}^N{\sigma }_n^z{\sigma }_{n+1}^z,Z({\beta }_1)=\text{mbox}Tr\text{mbox}{e}^{-{\beta }_1\mathcal{H}}.$$

$$M={\bar{\nu }}_1{M}_{\text{th}}({\beta }_1)+{\bar{\nu }}_2{M}_{\text{th}}({\beta }_2).$$

$$M=M_{\text{th}}(\beta ),$$

$$\tanh 2\beta E={\bar{\nu }}_1\tanh 2{\beta }_{1}E+{\bar{\nu }}_2\tanh 2{\beta }_{2}E.$$

$$f_{\mathrm{al},\mathrm{op}}=\frac{1}{4}\left[1\pm 2\frac{\zeta +{\zeta }^{N-1}}{1+{\zeta }^N}+\frac{{\zeta }^2+{\zeta }^{N-2}}{1+{\zeta }^N}\right],N\ge 2.$$

$$w_{\mathrm{al},\mathrm{op}}({\beta }_1)=\frac{1}{2}{\bar{\nu }}_1(1\mp \tanh 2{\beta }_{1}E),$$

$$⟨j_1⟩$$

$$\langle j\rangle=\tfrac{1}{4}\bar{\nu}_{1}\bar{\nu}_{2}(1+\zeta^{2})\frac{1+\zeta^{N-2}}{1+\zeta^{N}}\big{[}\tanh 2\beta_{2}E-\tanh 2\beta_{1}E\big{]}.$$

$$⟨\mathcal{Q}⟩\simeq ⟨j⟩N\tau ,\tau \to \infty ,$$

$$\delta {\beta }_1={\bar{\nu }}_2\delta {\beta }_{12},\delta {\beta }_2=-{\bar{\nu }}_1\delta {\beta }_{12}.$$

$$⟨\delta j⟩={\lambda }_{\mathrm{T}}\delta T,{\lambda }_{\mathrm{T}}={\bar{\nu }}_1{\bar{\nu }}_2\frac{(1-{\zeta }^2{)}^2(1+{\zeta }^{N-2})}{2(1+{\zeta }^2)(1+{\zeta }^N)}{\beta }^{2}E{k}_{\mathrm{B}}.$$

$$P(s,\vec{Q};0)={\delta }_{\vec{Q},\vec{0}}P(s;0).$$

$$\frac{\text{mbox}dP(s,\vec{Q};\tau )}{\text{mbox}d\tau }$$

$$|P(\vec{Q};\tau )⟩=\sum _{s}P(s,\vec{Q};\tau )|s⟩.$$

$$|\hat{P}(\vec{p};\tau )⟩=\sum _{\vec{Q}}\text{mbox}{e}^{\mathrm{i}\vec{p}\cdot \vec{Q}}|P(\vec{Q};\tau )⟩,$$

$$\frac{\text{mbox}d|\hat{P}(\vec{p};\tau )⟩}{\text{mbox}d\tau }=\hat{\mathcal{M}}(\vec{p})|\hat{P}(\vec{p};\tau )⟩,$$

$$\hat{\mathcal{M}}(\vec{p})={\bar{\nu }}_1{\hat{\mathcal{M}}}_{\mathrm{th}}(p_1;{\beta }_1)+{\bar{\nu }}_2{\hat{\mathcal{M}}}_{\mathrm{th}}(p_2;{\beta }_2)$$

$$\widehat{{\cal M}}_{\rm th}(p_{a};\beta_{a})=\tfrac{1}{2}\sum_{n=1}^{N}\big{(}\sigma_{n}^{x}\mbox{e}^{-\tfrac{1}{2}{\rm i}p_{a}\sigma_{n}^{z}(\sigma_{n-1}^{z}+\sigma_{n+1}^{z})}-1\big{)}\big{[}1-\tfrac{1}{2}\gamma_{a}\,\sigma_{n}^{z}(\sigma_{n-1}^{z}+\sigma_{n+1}^{z})\big{]}.$$

$$\hat{\mathcal{M}}(\vec{p})$$

$$C(\vec{p})$$

$$D(\vec{p})$$

$$\tilde{\mathcal{M}}(\vec{p})={\rho }_{\mathrm{eq}}^{-\frac{1}{2}}({\beta }_{\ast })\hat{\mathcal{M}}(\vec{p}){\rho }_{\mathrm{eq}}^{\frac{1}{2}}({\beta }_{\ast }),$$

$${\tilde{\sigma }}_n^x({\beta }_{\ast })$$

$$\tilde{C}(\vec{p},{\beta }_{\ast })$$

Full text

Keywords: driven Ising model, time-integrated energy flux, large deviation function, exact solution

LPT Orsay 16/xx

1. 1 Introduction

2. 2 Ising model coupled to two thermostats

3. 3 Energy current between the thermostats

4. 4 Extended master operator: definition and diagonalization

5. 4.1 Extended master operator $\widehat{\cal M}$

6. 4.2 Symmetrizing the master operator

7. 4.3 Transformation to fermion operators

8. 4.4 Diagonalizing in terms of fermion operators

9. 5 Joint probability distribution of the time-integrated energy currents

10. 5.1 Joint probability distribution $P(\vec{Q};\tau)$ of the time-integrated energy currents

11. 5.2 Rewriting $P(\vec{Q};\tau)$

12. 5.3 Finite time fluctuation relation for $P(\vec{Q};\tau)$

13. 6 Statistics of the time-integrated energy current

14. 6.1 Distribution $P({\cal Q};\tau)$ of the time-integrated energy current

15. 6.2 Cumulants of ${\cal Q}$ in the long-time limit

16. 6.3 Large deviation function of the time-integrated current ${\cal Q}/\tau$

17. 6.4 Infinite size chain at finite effective temperature

18. 7 Various physical effects

19. 7.1 Kinetic effects

20. 7.2 One thermostat at zero temperature

21. 7.3 Kinetic effects when colder thermostat is at zero temperature

22. 8 Large size and low effective temperature

23. 8.1 Parameters at low effective temperature

24. 8.2 Finite chain at zero effective temperature

25. 8.3 Infinite size chain at low effective temperature

26. 8.4 Interpretation of the scaling regimes

27. 9 Conclusion

28. A Behavior of coefficients $\Sigma_{n}(N,\gamma)$

1 Introduction

Since a few decades statistics of the currents that characterize an out-of equilibrium state have been intensely studied both experimentally and theoretically. Indeed the fluctuations of these currents in small systems are non-negligible with respect to their mean value, and they now can be investigated at nano scale thanks to very fast technological improvements [1, 2]. Meanwhile, the theory of stochastic thermodynamics has been developed and the large fluctuations of time-integrated currents in out-of-equilibrium systems have been shown to obey generic fluctuation relations. The latter have been derived under various hypotheses about the microscopic dynamics: deterministic or stochastic with either discrete or continuous degrees of freedom 111For a comprehensive review see the report by Seifert Ref.[3] and the references therein. In particular, for the case of stochastic Markovian dynamics with jumps between a finite number of configurations see Ref.[4].. These fluctuation relations for time-integrated currents quantify how the second law of thermodynamics, valid for mean currents, is modified at the scale of fluctuations; they are linked in some way to the fluctuations of the time-integrated entropy production rate in the system 222Short introductions which point out the role of entropy are to be found e.g. in Refs.[5, 6, 7].. In particular the class of systems with a finite number of discrete degrees of freedom has provided firmly established fluctuation relations [4].

Besides these generic fluctuation relations based on symmetry arguments, solvable models have provided better insight into more detailed statistical properties of non-equilibrium stationary states (NESS). This is most valuable in the absence of any equivalent of the equilibrium Gibbs ensemble theory for the description of NESS. In particular two paradigmatic kinetic models where a stationary current of particles or energy quanta flows from one reservoir to another have been widely investigated under various forms. On the one hand one-dimensional systems of particles endowed with a simple exclusion process and non-equilibrium open boundary conditions; such models describe particle exchange between two reservoirs connected to both ends of the system and which have different chemical potentials (see reviews [8, 9]). On the other hand Ising spin chains (with nearest-neighbor ferromagnetic interactions) where all spins are flipped by one of two thermostats.

In this paper we will introduce and study analytically a particuiar version of an Ising chain coupled to two thermostats. We begin by briefly recalling a few exact analytic results about kinetic Ising models.

In 1963 Glauber [10] endowed the Ising spin chain with a stochastic dynamics in order to describe the relaxation of this chain to its canonical equilibrium, which is determined only by the Ising energy and a given temperature $T$. A spin flip is interpreted as an energy exchange with a thermostat at temperature $T$. A single spin is flipped at a time, and the corresponding Markovian process is described by a master equation in spin configuration space.

The relaxation to the canonical equilibrium is ensured by the choice of the transition rates made by Glauber: these are the simplest ones that obey the detailed balance with the canonical configuration probability. The solution to the full description of the approach to equilibrium in this kinetic model was made in successive steps. First Glauber determined the evolution of the average magnetization and spin-spin correlations, and studied the linear response to an applied magnetic field. In the early 1970’s higher order correlation functions were studied [11, 12]. In particular, Felderhof [12, 13] was the first one to apply the fermionization technique to the Glauber model and showed that the master equation is fully solvable: that is, for a system of $N$ spins the $2^{N}$ eigenvalues and eigenvectors of the Markov matrix were all found exactly.

Later kinetic models for the Ising chain have been introduced in order to investigate the non-equilibrium stationary state (NESS) sustained by this Ising chain when the spins are flipped by two thermostats at different temperatures.

Exact results about the stationary probability distribution of the spin configurations have been obtained through the determination of mean instantaneous quantities in various models [14, 15, 16, 17]. Analytical expressions for the large deviation function of the time-integrated energy current in the non-equilibrium stationary state (NESS) have been obtained for simpler models [18, 19]. The complete description of the time-integrated energy currents has been obtained for a model where thermal contact between two thermostats is ensured by the interaction inside a set of independent Ising spin pairs, where each thermostat flips only one spin in the pair according to the corresponding Glauber dynamics [20]. The explicit joint probability of the cumulative heats received from each thermostat at any time and the analytical expression for the large deviation function of the time-integrated heat transfer from one thermostat to the other were obtained 333In the case of interacting Ising spin pairs one can obtain a partial description of the energy transfer from one thermostat to the other: the generating function for the long time cumulants per unit of time can be calculated analytically [21]. . The explicit stationary probability distributions of microscopic configurations have also been obtained for other archetypal models: the asymmetric exclusion process [9] and several variants of the zero-range process [22, 23]. The generating function for the cumulants of the time-integrated particle current have been obtained by sophisticated methods for various models endowed with an simple exclusion process[9].

In this work we study the Ising chain with a ferromagnetic nearest-neighbor coupling $E$, a finite number $N$ of spins, and periodic boundary conditions. The chain is coupled to two thermostats at temperatures $T_{1}$ and $T_{2}$ in the simplest of all possible ways: each spin may be reversed by either thermostat according to Glauber transition rates with inverse time constants (inverse time scales of random jumps) $\nu_{1}$ and $\nu_{2}$, respectively. These are kinetic parameters which depend on the microscopic dynamics of the system, as opposed to the thermodynamic parameters $T_{1}$ and $T_{2}$ of the energy reservoirs. The amount of energy received by the chain for each spin flip is equal to $-E$, [math], or $+E$. In the following all energies will be expressed as multiples of $4E$. We will take $T_{1}>T_{2}$ throughout this work. We rescale the physical time $t$ as $\tau=(\nu_{1}+\nu_{2})t$ and the kinetic parameters as $\bar{\nu}_{a}=\nu_{a}/(\nu_{1}+\nu_{2})$, where $a=1,2$.

We are interested in the joint probability $P(Q_{1},Q_{2};\tau)$ for the stochastic energy amounts $Q_{1}$ and $Q_{2}$ received by the Ising chain from the thermostats during a given time $\tau$. Then the probability for the time-integrated energy current ${\cal Q}$ (or net total energy that has flowed) from thermostat $1$ to thermostat $2$ during time $\tau$ is obtained as the marginal probability for the variable ${\cal Q}=\tfrac{1}{2}(Q_{1}-Q_{2})$. (We recall that $Q_{a}$ (with $a=1,2$) is an integer.)

The key to solvability is the observation 444This observation goes back at least to Garrido et al. [24], whose focus is however different from ours. that the sum of two Glauber rates at temperatures $T_{1}$ and $T_{2}$ is a Glauber rate with an effective kinetic parameter $\nu_{1}+\nu_{2}$ and at an intermediate temperature $T$ which is function of $T_{1}$, $T_{2}$ and $\bar{\nu}_{1}=\nu_{1}/(\nu_{1}+\nu_{2})$. As a consequence, on the one hand, the transition rates obey the canonical detailed balance and in a finite time the Ising spin chain reaches its stationary state where the probability for a spin configuration is the Boltzmann-Gibbs weight at the effective temperature $T(T_{1},T_{2},\bar{\nu}_{1})$. Then the net instantaneous energy current on each site has a zero mean, $\langle j\rangle=0$, but the contribution to this mean current from each thermostat does not vanish, $\langle j_{1}\rangle=-\langle j_{2}\rangle\neq 0$.

In order to deal with the extended Markov matrix which governs the evolution of the Fourier transform of the joint probability $P(s,Q_{1},Q_{2};\tau)$ for spin configurations $s$ and exchanged quantities $Q_{1}$ and $Q_{2}$, we extend the original method introduced by Felderhof [12, 13] for the Markov matrix of the probability $P(s;\tau)$ of the spin configurations during the relaxation to equilibrium for an Ising chain coupled to a single thermostat. This extended method yields all eigenvalues and eigenvectors of the extended master equation. It allows us to calculate the Fourier transform of $P(Q_{1},Q_{2};\tau)$ and $P({\cal Q};\tau)$ 555We use the same symbol $P$ for various different probabilities; the meaning will always be clear. at any time $\tau$. The system fulfills the hypotheses of various generic fluctuation relations, (5.32)-(5.35) and (6.20)-(6.21), which are indeed satisfied by the explicit expressions for the involved quantities.

From the expression for the Fourier transform of the probability $P({\cal Q};\tau)$ of the time-integrated energy current ${\cal Q}=\tfrac{1}{2}(Q_{1}-Q_{2})$, we obtain the explicit expression of the generating function for the infinite time limit of the cumulants of ${\cal Q}$ per site and unit of time, to be denoted as $\langle{\cal Q}^{n}\rangle_{\rm c}/N\tau$. The $n$th cumulant (per site and unit of time) of interest, $\lim_{\tau\to\infty}\langle{\cal Q}^{n}\rangle_{\rm c}/(N\tau)$, appears to be a $n$th degree polynomial in two variables ${\sf A}$ and ${\sf B}$ that are combinations of the thermodynamic and kinetic parameters,

[TABLE]

where $\gamma_{a}=\tanh 2\beta_{a}E$ for $a=1,2$ 666They are the same as the $A$ and $B$ of Ref. [20], except that our ${\sf B}$ has a minus sign compared to $B$, due to an inversion of the roles of the two thermostats.. These polynomials have coefficients $\Sigma_{n}(N,\gamma)$ which depend on the system size $N$ and the inverse effective temperature $\beta=(1/2E)\operatorname{artanh}\gamma$. They generalize the constant-coefficient polynomials that appeared in work by Cornu and Bauer [20] for a model where each thermostat flips only the spin on a given site. Although their model is different from the present chain with $N=2$, its various symmetries render its energetics identical to that of the present $N=2$ system 777 Properties of their model that are invariant by a global spin flip are equivalent to the properties of our system that are left-right invariant along the chain with $N=2$..

The explicit solution for the long time cumulants per site and unit of time allows one to investigate several physical effects beyond the generic symmetry relations. Indeed kinetic and dissipation effects specific to various regimes of the thermodynamic and kinetic parameters can be investigated. They are summarized in the conclusion.

Moreover size effects generated by the interaction between spins can be controlled. The model makes sense only if the effective temperature $\beta$ is finite ($\gamma\neq 1$). Then the large deviation function exists in the infinite size limit and all long time cumulants per unit of time for the whole chain, $\lim_{\tau\to\infty}\langle{\cal Q}^{n}\rangle_{\rm c}/\tau$, are proportional to the size $N$ of the chain at leading order in $N$. In the double limit where the effective temperature $1/\beta$ goes to zero while the size $N$ goes to infinity, all these cumulants are proportional to $(1-\gamma)N$ at leading order in $N$ and $1-\gamma$. We notice that the factor $(1-\gamma)$ disappears if one considers the rescaled cumulants per unit of time when the unit of time is the magnetization relaxation time $\tau_{\textrm{rel}}$, which is equal to $[(\nu_{1}+\nu_{2})(1-\gamma)]^{-1}$. In this double limit the variables ${\sf A}$ and ${\sf B}$ defined in (1.1) vanish as $1-\gamma$ while the coefficients $\lim_{N\to\infty}\Sigma_{n}(N,\gamma)$ with $n\geq 2$ diverge. As a consequence, the leading behavior of the rescaled cumulants per unit of time is a random walk contribution of order $N$, whereas the first correction to it is not of order zero in $N$ when $1-\gamma\to 0$. In fact one has to consider two scaling regimes where the increase of $N\gg 1$ is related to the decrease of $1-\gamma\ll 1$; we exhibit how the first correction in the cumulants depends upon the ratio of the spin-spin correlation length $\xi(T)$ and the size $N$.

This paper is set up as follows. In section 2 we define the Ising model between two thermostats. In section 3 we discuss the instantaneous energy current, whose average $\langle j\rangle$ per site we determine by elementary means. In section 4 we define and diagonalize the master operator in the extended space of spin configurations and energies $Q_{1}$ and $Q_{2}$ received by the spin chain from both thermostats during a time interval $\tau$, and in section 5 we determine the Fourier transform of the joint probability $P(Q_{1},Q_{2};\tau)$. We check that the explicit expression of $P(Q_{1},Q_{2};\tau)$ in the present model does satisfy the fluctuation relations (5.32)-(5.35) which are retrieved from general considerations. In section 6 we obtain the Fourier transform of the probability $P({\cal Q};\tau)$ of the time-integrated energy current ${\cal Q}$ from one thermostat to the other during a time $\tau$. We determine the cumulants per site and unit of time of ${\cal Q}$ in the long-time limit and discuss their structure. In section 7 we study physical effects in various regimes of the thermodynamic and kinetic parameters for a finite chain. In section 8 we consider a large size chain at very low effective temperature: from the study of some divergent coefficients performed in Appendix A we exhibit the first correction to the leading $N$-behavior of the cumulants. In section 9 we briefly conclude.

2 Ising model coupled to two thermostats

We consider a chain of Ising spins $s_{n}=\pm 1$, where $n=1,2,\ldots,N$ and $N\geq 2$ is an arbitrary integer. A configuration $s=(s_{1},s_{2},\ldots,s_{N})$ of the Ising model has an energy $H(s)$ given by

[TABLE]

where we adopt the periodic boundary condition $s_{N+n}=s_{n}$. We will be concerned with time dependent probability $P(s;\tau)$ in configuration space.

In a formalism that goes back at least to Kadanoff and Swift [25] we associate with each $s$ a ket $|s\rangle=\otimes_{n=1}^{N}|s_{n}\rangle$. A probability $P(s;\tau)$ is then represented by a time dependent ket

[TABLE]

Since the classical discrete variables $s_{n}$ all commute, the Ising model has no dynamics of itself. In 1963 Glauber [10] stipulated that when the system is in contact with a thermostat at temperature $T_{1}$, then in a configuration $s$ the spin $s_{n}$ on the $n$th lattice site may reverse its state with a transition rate given in dimensionless time $\tau=(\nu_{1}+\nu_{2})t$ (where $\nu_{a}$ is an inverse time) by

[TABLE]

where $\bar{\nu}_{1}=\nu_{1}/(\nu_{1}+\nu_{2})$ is an inverse time, $\gamma_{1}=\tanh 2\beta_{1}E$, and $\beta_{1}=1/k_{\rm B}T_{1}$ is the inverse temperature. The ket $|P(\tau)\rangle$ then evolves according to the master equation

[TABLE]

with a “master operator” ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})$ whose expression is originally due to Felderhof [12, 13],

[TABLE]

in which $\sigma_{n}^{z}$ and $\sigma_{n}^{x}$ are the usual Pauli spin operators defined by $\sigma_{n}^{z}|s_{n}\rangle=s_{n}|s_{n}\rangle$ and $\sigma_{n}^{x}|s_{n}\rangle=|-s_{n}\rangle$. The master equation is easily shown to have the unique stationary state

[TABLE]

in which we have

[TABLE]

We remark that $H(s)$ in Eq. (2.1) is an eigenvalue of ${\cal H}$.

By means of fermionization the operator ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})$ may be completely diagonalized and all its eigenvectors determined [12, 13]. That means that, in principle, this problem is fully understood. Recent renewal of interest in kinetic Ising models, as mentioned in the introduction, is due to the development of the study of non-equilibrium stationary state systems. With this perspective in mind we will here couple the same system to two thermostats at inverse temperatures $\beta_{1}$ and $\beta_{2}$ and acting with rates $\nu_{1}$ and $\nu_{2}$, respectively. The total operator describing the system, denoted by $M$, then becomes a weighted sum of the Glauber operators at inverse temperatures $\beta_{1}$ and $\beta_{2}$ ,

[TABLE]

with $\bar{\nu}_{1}+\bar{\nu}_{2}=1$. In this work we study this model in detail.

Normally a system in contact with two reservoirs in different equilibrium states will tend to a stationary state. Usually the precise properties of such a state are not easy to determine. In the present case a simplification occurs since the operator $M$ of equation (2.8) can be rewritten as

[TABLE]

where $\beta$ represents an effective temperature intermediate between $\beta_{1}$ and $\beta_{2}$ given by

[TABLE]

We will employ below the abbreviations $\gamma=\tanh 2\beta{E}$ and $\gamma_{a}=\tanh 2\beta_{a}{E}$ for $a=1,2$.

It follows that the stationary state in this case actually happens to be equal to the equilibrium state at the effective temperature 888The same observation was made by Cornu and Bauer [20] for their two-spin system with only two energy levels.. This does not mean that we immediately know the answers to the questions raised above considering the energy injection and dissipation. It means, however, that they can be calculated, which is what we do in this work.

3 Energy current between the thermostats

We consider the system in its stationary state, that is, in the equilibrium state at inverse temperature $\beta$. The reversal of a spin involves an energy change only if the two neighbors of that spin are mutually parallel. Let $f_{\rm al}$ be the fraction of all spins that have their two neighbors mutually parallel and aligned to it, and $f_{\rm op}$ the fraction of those having them mutually parallel and opposite to it. The indicator function for a spin $s_{n}$ aligned with (opposite to) both of its neighbors is $\tfrac{1}{4}(1\pm s_{n-1}s_{n})(1\pm s_{n}s_{n+1})$. Ensemble averaging this by standard methods which lead to the result $\langle s_{n}s_{n+r}\rangle=[\zeta^{r}+\zeta^{N-r}]/[1+\zeta^{N}]$, with $\zeta=\tanh\beta E$, we obtain for the periodic Ising chain

[TABLE]

We consider the action on this system by the operator $\bar{\nu}_{1}{M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})$. The spins of the two classes $f_{\rm al}$ and $f_{\rm op}$ are reversed with transition rates expressed in the dimensionless time $\tau=(\nu_{1}+\nu_{2})t$ as

[TABLE]

respectively. (The minus sign corresponds to $w_{\rm al}$.) Let $\langle j_{1}\rangle$ be the net average instantaneous energy current per unit of chain length from thermostat $1$ into the system. Expressed in units of $4E$ it reads

[TABLE]

A similar expression holds for the net average current $\langle j_{2}\rangle$ from thermostat $2$ into the system under the action of $\bar{\nu}_{2}{M_{\scriptscriptstyle{\text{th}}}}(\beta_{2})$. From (2.10) and (3.3) together with the relation $\tanh 2\beta E=\zeta^{2}/(1+\zeta^{2})$, we get that $\langle j_{1}\rangle+\langle j_{2}\rangle=0$ : in a stationary state the finite system cannot accumulate energy. Then $\langle j\rangle=\langle j_{1}\rangle=-\langle j_{2}\rangle$ represents the net average energy current per site (= unit of chain length) that traverses the system from thermostat $1$ to thermostat $2$. The most elegant expression for this quantity is obtained by remembering that $\bar{\nu}_{1}+\bar{\nu}_{2}=1$ and writing it as $\langle j\rangle=\bar{\nu}_{2}\langle j_{1}\rangle-\bar{\nu}_{1}\langle j_{2}\rangle$ with the result

[TABLE]

This is our ‘direct’ result for the average instantaneous energy current density, valid in a finite periodic chain. Let ${\cal Q}$ stand for the net total energy (i.e. time-integrated energy current), expressed in units of $4E$, that during a time interval $[0,\tau]$ passes through the system from thermostat 1 to thermostat 2. We will let $\bar{\jmath}\equiv{\cal Q}/N\tau$ stand for the dimensionless integrated current per site and per unit of time. In the long-time limit $\langle\bar{\jmath}\rangle=\langle j\rangle$ and $\langle{\cal Q}\rangle$ diverges with the time $\tau$ as

[TABLE]

and $\langle j\rangle$ given by (3.4). There is no such simple method to calculate the higher order moments $\langle{\cal Q}^{n}\rangle$ for $n\geq 2$. The work of this paper will lead us to expressions for the cumulants $\langle{\cal Q}^{n}\rangle_{\rm c}$. It will confirm equation (3.4) as a particular case.

It is of some interest to consider the linearization in temperature around the equilibrium state where $\beta_{1}=\beta_{2}=\beta$. Let $\beta_{a}=\beta+\delta\beta_{a}$ for $a=1,2$ and let us set $\delta\beta_{12}=\delta\beta_{1}-\delta\beta_{2}=-\delta T/k_{B}T^{2}$, where $T=1/k_{\scriptscriptstyle B}\beta$ ($k_{\scriptscriptstyle B}$ Boltzmann constant) and the infinitesimal temperature difference is $\delta T=T_{1}-T_{2}$. Because of the relation (2.10) we then have

[TABLE]

Calling the linearized current $\delta j$, we obtain from (3.4)

[TABLE]

The heat conduction coefficient $\lambda_{\rm T}$ tends to zero in both limits $\beta\to 0$ and $\beta\to\infty$, with $E$ fixed.

4 Extended master operator: definition and diagonalization

4.1 Extended master operator $\widehat{\cal M}$

Each spin reversal is due to either ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})$ or ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{2})$, and each spin reversal involves the injection or the release of a quantum of energy equal to 0 or to $\pm 4E$. Let the integers $Q_{1}$ and $Q_{2}$ denote the total energy, measured in units of $4E$, furnished to the system by the operators ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{1})$ and ${M_{\scriptscriptstyle{\text{th}}}}(\beta_{2})$, respectively, in a time interval of duration $\tau$. For $T_{1}>T_{2}$ both $Q_{1}$ and $-Q_{2}$ will have positive expectation values. We will write $\vec{Q}=(Q_{1},Q_{2})$. We are interested in the joint probability distribution $P(s,\vec{Q};\tau)$, which satisfies $\sum_{\vec{Q}}P(s,\vec{Q};\tau)=P(s;\tau)$ and the initial condition

[TABLE]

Let $s^{n}$ denote the configuration obtained from $s$ by flipping the spin at site $n$, and let $\Delta Q_{n}(s)$ denote the increment in either $Q_{1}$ or $Q_{2}$ associated with the jump from $s$ to $s^{n}$, that is, $\Delta Q_{n}(s)=\tfrac{1}{2}s_{n}(s_{n-1}+s_{n+1})$. (For the reversed spin flip at site $n$, namely the jump from $s^{n}$ to $s$, the increment in either $Q_{1}$ or $Q_{2}$ is $\Delta Q_{n}(s^{n})=-\Delta Q_{n}(s)$.) The probability $P(s,\vec{Q};\tau)$ then obeys the balance equation

[TABLE]

By analogy with the representation (2.2) of $P(s;\tau)$, we represent the probability $P(s,\vec{Q};\tau)$ by the time dependent ket

[TABLE]

We consider the Fourier transformed ket

[TABLE]

where $\vec{p}=(p_{1},p_{2})$ with $-\pi<p_{1},p_{2}\leq\pi$. Upon taking the Fourier transform of the balance equation (4.2) we get the evolution equation for the ket (4.4),

[TABLE]

in which

[TABLE]

where, by analogy with (2.5),

[TABLE]

In this expression the operator $O_{n}\equiv\tfrac{1}{2}\sigma_{n}^{z}(\sigma_{n-1}^{z}+\sigma_{n+1}^{z})$, whose eigenvalues are $1$, [math], and $-1$, has the properties $O_{n}^{2k}=O_{n}^{2}=\tfrac{1}{2}(1+\sigma_{n-1}^{z}\sigma_{n+1}^{z})$ for $k\geq 1$ and $O_{n}^{2k+1}=O_{n}$ for $k\geq 0$. Hence $\mbox{e}^{-{\rm i}p_{a}O_{n}}=1-({\rm i}\sin p_{a})O_{n}+(\cos p_{a}-1)O_{n}^{2}$. As a consequence expression (4.6) may be rewritten as

[TABLE]

in which

[TABLE]

These coefficients are real when $p_{1}$ and $p_{2}$ are pure imaginary.

4.2 Symmetrizing the master operator

We apply to $\widehat{\cal M}(\vec{p})$ a similarity transformation and define

[TABLE]

with a ${\beta_{*}}(\vec{p})$ left to be determined in such a way that $\widetilde{\cal M}(\vec{p})$ be Hermitian. The only nontrivial relation needed to find an explicit expression for (4.10) is [12, 13]

[TABLE]

which is easily derived. The result is that $\widehat{\cal M}(\vec{p})$ of equation (4.8) becomes an expression $\widetilde{\cal M}(\vec{p})$ which is of the same form as (4.8) but with $C$ and $D$ of equation (4.9) replaced with $\widetilde{C}$ and $\widetilde{D}$, respectively, where

[TABLE]

We now choose ${\beta_{*}}$ such that the coefficient $\widetilde{D}({\beta_{*}})$ of the non-Hermitian term vanishes. This amounts to taking

[TABLE]

where $C(\vec{p})$ and $D(\vec{p})$ are given by (4.9). We see that ${\beta_{*}}(\vec{0})=\beta$ and that ${\beta_{*}}(\vec{p})$ is real when $p_{1}$ and $p_{2}$ are pure imaginary. As a result the symmetrized operator $\widetilde{\cal M}(\vec{p})$ takes the form

[TABLE]

in which $C_{*}$ is given by

[TABLE]

All $\vec{p}$ dependence of $\widehat{\cal M}(\vec{p})$ is seen to enter through the single coefficient $C_{*}(\vec{p})$.

After substituting (4.9) in (4.12) and (4.12) in (4.15) we find that this quantity may be written as

[TABLE]

where

[TABLE]

with ${\sf A}$ and ${\sf B}$ given by (1.1) in the Introduction. These coefficients will appear again in our final results in section 6.

4.3 Transformation to fermion operators

We define fermionic quasi-particles by means of the Jordan-Wigner [26] transformation

[TABLE]

The vacuum state of these $c$-particles is the state $|1\rangle$ defined in (2.6). It is now straightforward to express $\widetilde{\cal M}(\vec{p})$ in terms of these fermion operators. We find from (4.14)

[TABLE]

with the understanding that the creation and annihilation operators whose indices exceed $N$ are defined by

[TABLE]

in which ${\cal N}=\sum_{n=1}^{N}c^{\dagger}_{n}c_{n}$ is the operator for the total number of quasi-particles.

4.4 Diagonalizing in terms of fermion operators

For convenience we hence restrict ourselves to even $N$. We define fermion operators $\eta^{\dagger}_{q}$ and $\eta_{q}$ by

[TABLE]

where the wavenumber $q$ runs through the $N$ values

[TABLE]

Equation (4.21) is easily inverted to find the $\eta^{\dagger}_{q}$ and $\eta_{q}$ in terms of the $c^{\dagger}_{n}$ and $c_{n}$. This equation guarantees the periodicity conditions (4.20) in the subspace where ${\cal N}$ is even. In that subspace equation (4.21) may also be used in (4.19) for $n=N+1$ and $n=N+2$. Obviously the $c$ vacuum $|1\rangle$ is also the $\eta$ vacuum.

Applying transformation (4.21) to (4.19) we get

[TABLE]

valid in the subspace with an even number ${\cal N}$ of $c$-particles. 999In the subspace with an *odd * number of $c$-particles $\widetilde{\cal M}(\vec{p})$ takes a slightly different form, as discussed in detail in references [12, 13, 27, 28, 29] We will not need that form in this work., and where the coefficients $C_{q}$ and $D_{q}$ are given by

[TABLE]

where the $\vec{p}$ dependence comes in through the $q$ independent coefficient $C_{*}(\vec{p})$ defined in (4.15). Extending the approach of Ref. [12, 13] to nonzero $\vec{p}$ we define angles $\chi_{q}$ (that are in general complex) by

[TABLE]

and perform in the space of the pair $\{\eta_{q},\eta^{\dagger}_{-q}\}$ a Bogoliubov-Valatin [30, 31] operator rotation

[TABLE]

It is useful to note that $\chi_{-q}=-\chi_{q}$. Upon using (4.26) to transform (4.23) to $\xi$ operators we find the diagonal form

[TABLE]

where

[TABLE]

and

[TABLE]

From (4.24) and (4.28) it is easily seen that $\mu_{q}=\mu_{-q}$. Upon combining both equations we get for $\mu_{q}$ the explicit expression

[TABLE]

with $C_{*}(\vec{p})$ given by (4.16). We note that the generally complex quantity $C_{*}$ does not depend on $q$ and that the $\vec{p}$ dependence of this diagonalization process comes in only through $C_{*}(\vec{p})$. For $\vec{p}=\vec{0}$ our results for $\mu_{q}(\vec{p})$ reduces to that of Ref. [12, 13], namely $\mu_{q}(\vec{0})=1-\gamma\cos q$.

A different way, useful for later, to write the eigenvalue $\mu_{q}(\vec{p})$ is

[TABLE]

where $\Theta$ has been defined in (4.17). We observe for later use that

[TABLE]

It is convenient to rewrite the diagonalized form (4.27) of the master operator as

[TABLE]

where the symmetry property $\mu_{q}=\mu_{-q}$ has been employed and where, here and hence, ‘$q>0$’ refers to the $\tfrac{1}{2}N$ positive values of $q$ among those given in (4.22).

5 Joint probability distribution of the time-integrated energy currents

5.1 Joint probability distribution $P(\vec{Q};\tau)$ of the time-integrated energy currents

Let $P(\vec{Q};\tau)$ be the probability that at time $\tau$ the time-integrated energies furnished by the thermostats $1$ and $2$ to the system, counted in units of $4E$, have the values $Q_{1}$ and $Q_{2}$ respectively. Then according to (4.3) this probability distribution is given by $P(\vec{Q};\tau)=\sum_{s}\,\langle s|P(\vec{Q};\tau)\rangle$ and upon inverting (4.4) we find

[TABLE]

where $\widehat{P}(\vec{p};\tau)=\sum_{s}\langle s|\widehat{P}(\vec{p};\tau)\rangle$. The evolution equation (4.5) may be formally solved as

[TABLE]

where $|\widehat{P}(\vec{p};0)$ is the Fourier transform of the initial state $P(\vec{Q};0)\rangle$. Our protocol will be to take for the initial configuration the equilibrium state at an arbitrary inverse temperature $\beta_{0}$. Moreover, since at time $\tau=0$ no energy exchange has taken place yet we choose this probability concentrated in $\vec{Q}=0$, that is, $|P(\vec{Q};0))=\delta_{\vec{Q},\vec{0}}|P_{\rm eq}(\beta_{0}))$. With the definitions (2.6) the probability $\widehat{P}(\vec{p};\tau)$ reads

[TABLE]

This expression takes advantage of the fact that ${\widehat{\cal M}}$ is block diagonal in the subspaces of fixed $\vec{p}$.

5.2 Rewriting $P(\vec{Q};\tau)$

According to the complete diagonalization performed in section 4 the matrix element in the Fourier transform (5.3) is an expectation value in the $\eta$ vacuum and can be rewritten as

[TABLE]

where we have passed to the symmetrized operator $\widetilde{\cal M}(\vec{p})$ and $|1\rangle$ denotes the $\eta$-vacuum. We decompose the $\eta$-vacuum as

[TABLE]

where $|0_{q}0_{-q}\rangle$ is the state in which the quasi-particles of wavenumbers $\pm q$ are absent,

[TABLE]

and $\langle 0_{q}0_{-q}|0_{q}0_{-q}\rangle=1$.

It is useful to rewrite the time evolution operator (4.33) as

[TABLE]

with

[TABLE]

where we have not indicated explicitly the $\vec{p}$ dependence of the $\mathbb{X}_{q}$, $\xi_{q}^{\dagger}$, and $\xi_{q}^{\dagger}$ operators. Furthermore we may express the Hamiltonian in terms of fermion operators, which yields

[TABLE]

where

[TABLE]

where we included $\mathbb{D}_{q}$ for later reference. We now use the fact that $\mathbb{X}_{q}$ [in view of relations (5.8) and (4.26)] and $\mathbb{H}_{q}$ [in view of (5.9) and (5.10)] are quadratic in the $\eta$ operators and that therefore

[TABLE]

Upon using (5.5) in (5.4) we may factorize $\widehat{P}(\vec{p};\tau)$ according to

[TABLE]

in which

[TABLE]

Since $\exp(K\mathbb{H}_{q})$ (for $K={\beta_{*}}E$ or $K=(2\beta_{0}-{\beta_{*}})E$) and $\exp(-\mu_{q}\mathbb{X}_{q}\tau)$ are both quadratic in the fermion operators, they act in the two-dimensional space spanned by the vacuum $|0_{q}0_{-q}\rangle$ defined above and the two-particle state $|1_{q}1_{-q}\rangle$ defined by

[TABLE]

To make the action of $\exp(K\mathbb{H}_{q})$ more explicit we expand the exponential using the relations

[TABLE]

which are easily checked. One then obtains

[TABLE]

in which

[TABLE]

For two specific choices of $K$ we will use below the notation

[TABLE]

We have to similarly expand $\exp(-\mu_{q}\mathbb{X}_{q}\tau)$ and obtain along the same lines

[TABLE]

With the aid of the relations (4.26) we can turn this into an expansion of the form

[TABLE]

After a fair amount of algebra one finds for the coefficients $a_{i}$ the expressions

[TABLE]

We substitute now expansions (5.16) for $K={\beta_{*}}E$ or $K=(2\beta_{0}-{\beta_{*}})E$ and (5.20) in (5.13). Taking into account that creation (annihilation) operators acting to the left (to the right) on the $\eta$-vacuum give zero, we may suppress the corresponding terms in the expansions and can write (5.13) as

[TABLE]

In the last line each term correspond to a sequence of creations and annihilations as one reads the first line from the right to the left, starting from and ending up in the vacuum. We may now substitute the values of the $a_{i}$, $b_{i}$, and $c_{i}$ found above.

After some algebra applied to (5.22) we may cast the $\Pi_{q}(\vec{p};\tau)$ in the form

[TABLE]

where $\mu_{q}$ is given by (4.31) while

[TABLE]

with

[TABLE]

and

[TABLE]

Combining (5.12) and (5.23) we get

[TABLE]

where the partition function defined in (2.7) reads

[TABLE]

It is easy to verify the relation

[TABLE]

Using it in (5.27) and substituting (5.27) in (5.1) we finally obtain

[TABLE]

This expression depends on the initial inverse temperature $\beta_{0}$ through the ratio $T_{q}(\vec{p};\beta_{0})/S_{q}(\beta_{0})$. It is possible to show with the aid of considerable algebra that

[TABLE]

which together with (5.27), (5.29), and (4.32) for $\mu_{*}(\vec{0})$ implies that $\widehat{P}(\vec{0};\tau)=1$, equivalent to the normalization condition $\sum_{\vec{Q}}P(\vec{Q};\tau)=1$.

5.3 Finite time fluctuation relation for $P(\vec{Q};\tau)$

One can check on the explicit expression (LABEL:resfinPi) that $P(Q_{1},Q_{2};\tau)$ obeys a finite time fluctuation relation: by virtue of the relation $\ln\sqrt{({\sf A}+{\sf B})/({\sf A}-{\sf B})}=2(\beta_{2}-\beta_{1})E$, the ratio of probabilities for opposite values of the couple $(Q_{1},Q_{2})$ is given at any time by

[TABLE]

In fact this relation relies on two key properties. First one can define the extended transition rates associated with the extended master operator ${\cal M}_{\textrm{ext}}$ such that $\mbox{d}P(s,\vec{Q};\tau)/\mbox{d}\tau=\sum_{s^{\prime},\vec{Q^{\prime}}}\langle s,\vec{Q}|{\cal M}_{\textrm{ext}}|s^{\prime},\vec{Q^{\prime}}\rangle P(s^{\prime},\vec{Q^{\prime}};\tau)$ and whose explicit expression is derived from the balance equation (4.2). These extended transition rates are defined between two triplets, each of which involves a spin configuration together with the two energies received from thermostats since the beginning of the considered history of the system: when the system is in spin configuration $s$ and the spin at site $n$ is flipped by thermostat $a$ (with $a=1,2$) they read $w_{n}^{(a)}\left(s,\Delta Q^{(1)}_{n}(s),\Delta Q^{(2)}_{n}(s)\right)=w_{n}(s;\beta_{a})$ where the expression for $\Delta Q^{(a)}_{n}(s)=\Delta Q_{n}(s)$ is given before (4.2), that for the other thermostat $b$ is $\Delta Q^{(b)}_{n}(s)=0$, and the transition rate between spin configuration $w_{n}(s;\beta_{a})$ is given in (2.3). These extended transition rates have the symmetry property obeyed by the transition rates $w_{n}(s;\beta_{a})$ for the two reversed transitions $s\to s_{n}$ and $s_{n}\to s$

[TABLE]

This symmetry can be considered as an extension of the so-called generalized detailed balance 101010Several terminologies can be found in the literature : ”local” detailed balance [32, 3], ”generalized” detailed balance [33] or ”modified” detailed balance [34, 6]. which involves only the transition rates between two configurations, and where the values of $\Delta Q^{(1)}_{n}(s)$ and $\Delta Q^{(2)}_{n}(s)$ are determined solely by the transition $s\to s_{n}$ (which is the case when a spin at a given site can be flipped by only one thermostat). The second key property arises from the considered protocol and the specificity of the stationary configuration probability in the model. Indeed the initial spin configuration distribution is the stationary configuration probability at the effective inverse temperature $\beta_{0}$ (with a nonvanishing mean current from thermostat $1$ to thermostat $2$). Besides in the present model the latter configuration probability is the canonical equilibrium distribution at inverse temperature $\beta_{0}$. The two key properties altogether allow one to apply usual arguments for the derivation of fluctuation relations. In the present case the precise argument is a mere transposition of that to be found for instance in [6, 35] where a transition between two spin configurations is caused by only one thermostat.

As in the case where the generalized detailed balance is met by the mere transition rates between spin configurations, the property (5.32) can be interpreted in terms of some time-integrated entropy variation as follows. When a thermostat at inverse temperature $\beta_{a}$ gives an energy $4EQ_{a}$ to the system, its entropy variation is $\Delta S_{a}=-4E\beta_{a}Q_{a}$. The exchange contribution $\Delta_{\rm exch}^{\beta_{1},\beta_{2}}S$ to the entropy variation of the system is defined as $-(\Delta S_{1}+\Delta S_{2})$, namely the opposite of the sum of the entropy variations of the two thermostats ; hence $\Delta_{\rm exch}^{\beta_{1},\beta_{2}}S=4E(\beta_{1}Q_{1}+\beta_{2}Q_{2})$. As in Ref.[6, 35] we introduce the excess exchange entropy variation of the system $\Delta_{\rm exch}^{\rm excs,\beta_{0}}S$ which is defined as the difference between the exchange entropy variation under the non-equilibrium external constraint $\beta_{1}\neq\beta_{2}$ and its value under the equilibrium condition $\beta_{1}=\beta_{2}=\beta_{0}$ and for the same values of the energies $4EQ_{1}$ and $4EQ_{2}$ received by the system: $\Delta_{\rm exch}^{\rm excs,\beta_{0}}S=\Delta_{\rm exch}^{\beta_{1},\beta_{2}}S-\Delta_{\rm exch}^{\beta_{0},\beta_{0}}S$. It also reads

[TABLE]

Hence the fluctuation relation (5.32) can be rewritten as $P(Q_{1},Q_{2};\tau)=\\ \exp[-\Delta_{\rm exch}^{\rm excs,\beta_{0}}S]P(-Q_{1},-Q_{2};\tau)$. As a consequence, the probability of the excess exchange entropy variation obeys a finite time fluctuation relation which takes the “universal” form

[TABLE]

6 Statistics of the time-integrated energy current

6.1 Distribution $P({\cal Q};\tau)$ of the time-integrated energy current

We now restrict our interest to the time-integrated current that during a time interval $[0,\tau]$ has traversed the system. It is defined as

[TABLE]

which may be integer or half-integer. It measures, in units $4E$, half the energy furnished to the system by thermostat 1 plus half the energy extracted from it by thermostat 2. Since for long times no energy can accumulate in the system, this quantity is, in the long time limit, equal to the time-integrated energy current. The particular definition (6.1) is motivated by the fact that it is antisymmetric under exchange of the two thermostats, which makes subsequent calculations easier.

Let $P({\cal Q};\tau)$ be the probability of ${\cal Q}$ at time $\tau$. This marginal probability of $P(\vec{Q};\tau)$ is obtained as $P({\cal Q};\tau)=\sum_{Q_{1},Q_{2}}\delta_{Q_{1}-Q_{2},2{\cal Q}}\,P(\vec{Q};\tau)$. We will from here on, for any $\vec{p}$  dependent quantity $X(\vec{p})$, employ the notation ${X}(p,-p)\equiv X^{\star}(p)$. From (LABEL:resfinPi) and the preceding definitions we then get

[TABLE]

in which

[TABLE]

We observe that $P({\cal Q};\tau)$ given by (6.2) and (6.3) still depends on the initial state parameter $\beta_{0}$. For the choice $\beta_{0}=\beta$ the system is in a stationary state for all $\tau\geq 0$; for $\beta_{0}\neq\beta$ it will asymptotically tend to that state.

We take advantage of the analyticity in $p$ of the integrand $\widehat{P}^{\star}(p;\tau)$ to point out that the moment generating function of $P({\cal Q};\tau)$, defined as $\langle\mbox{e}^{\lambda{\cal Q}}\rangle\equiv\sum_{{\cal Q}}\,\mbox{e}^{\lambda{{\cal Q}}}P({{\cal Q}};\tau)$, exists for all real $\lambda$ and is given by

[TABLE]

6.2 Cumulants of ${\cal Q}$ in the long-time limit

In the long-time limit the cumulants per site and unit of time $\langle{\cal Q}^{n}\rangle_{\rm c}/N\tau$ of the time-integrated energy current per site and unit of time are obtained from the scaled cumulant generating function $g_{N}(\lambda)$ defined as

[TABLE]

The cumulants of interest are the values of the derivatives of $g_{N}(\lambda)$ with respect to $\lambda$ taken at $\lambda=0$,

[TABLE]

According to (6.4) and (6.5) together with the explicit value (6.3) of $\widehat{P}^{\star}(-{\rm i}\lambda/2;\tau)$ we get

[TABLE]

The expression for $\mu_{q}^{\star}(p)=\mu_{q}(p,-p)$ is given by (4.31) where $\Theta(\vec{p})$ is defined in (4.17). We set $\theta(\lambda)=\Theta(-{\rm i}\lambda/2,{\rm i}\lambda/2)$ and get

[TABLE]

with

[TABLE]

in which ${\sf A}=\bar{\nu}_{1}\bar{\nu}_{2}(1-\gamma_{1}\gamma_{2})$ and ${\sf B}=\bar{\nu}_{1}\bar{\nu}_{2}(\gamma_{2}-\gamma_{1})$ depend only on the kinetic and thermodynamic parameters of the model. We notice that the expression of the scaled generating function for the time-integrated current of energy has a form similar to that for various currents of interest in the case of a system of diffusing particles with pair creation and annihilation [36]. This is due to a connection between the model considered by these authors and an Ising spin chain with Glauber dynamics.

We notice that expression (6.7) for $g_{N}(\lambda)$ can be obtained without knowing the explicit expression of the moment generating function $\langle\mbox{e}^{\lambda{\cal Q}}\rangle$ at any time $\tau$. Indeed, the evolution of $\langle\mbox{e}^{\lambda{\cal Q}}\rangle$ is Markovian, as shown by (5.4) with $p_{1}=-{\rm i}\lambda$ and $p_{2}={\rm i}\lambda$. Hence $g_{N}(\lambda)$ is the largest eigenvalue of $\widetilde{\cal M}(-{\rm i}\lambda,{\rm i}\lambda)$ and the operator expressions (5.7) and (5.8) lead to (6.7).

Eventually the cumulants of the time-integrated energy current per site and unit of time in the long-time limit are given by (6.6) and (6.8)

[TABLE]

where we have introduced

[TABLE]

We have indicated explicitly the dependence of $\Sigma_{n}(N,\gamma)$ on the size $N$ and the effective intermediate inverse temperature $\beta$ of the stationary state; we recall that $\beta$ is defined in terms of the parameter $\gamma$ through (2.10).The $\Sigma_{n}(N,\gamma)$ are monotonically increasing with $\gamma$.

Expressions for higher order cumulants may be derived by increasing algebraic effort. Expression (6.10) for the time-averaged energy current has to coincide with Eqs. (3.4)-(3.5) of section 3. Upon inserting the explicit expressions for both one obtains the identity

[TABLE]

where we recall that $\zeta=\tanh\beta E$ while $\gamma=\tanh 2\beta E$. Eq. (6.15) may be checked by explicit calculation. It shows that $\frac{1}{2}\leq\Sigma_{1}(N,\gamma)\leq 1$. We have not found similarly simple expressions for the $\Sigma_{n}(N,\gamma)$ with $n\geq 2$.

The expressions for cumulants, of which the first four are given in (6.10)-(6.13), have an interesting structure. The $n$th cumulant is an $n$th degree polynomial in the two variables ${\sf A}$ and ${\sf B}$ with coefficients $\Sigma_{1}(N,\gamma),\ldots,\Sigma_{n}(N,\gamma)$. The variables ${\sf A}$ and ${\sf B}$ depend on both thermostat temperatures $T_{1}$ and $T_{2}$ but are independent of the system size $N$. On the contrary the coefficients $\Sigma_{n}(N,\gamma)$ vary with the system size $N$, but depend only on the intermediate effective temperature $T$ and not on $T_{1}$ and $T_{2}$ separately. 111111Recall that $T$ depends also on the kinetic parameters $\nu_{1}$ and $\nu_{2}$ . We will analyze the $\Sigma_{n}(N,\gamma)$ in detail in the limit of large $N$ and low effective temperature $T$ in section 8.

When the two thermostats have equal temperatures, $T_{1}=T_{2}$, one has ${\sf B}=0$. Then only the even cumulants are nonzero, as must be the case when one considers the energy transfer between two thermostats at the same temperature. The even cumulants with $n\geq 4$ do not vanish: when $T_{1}=T_{2}$ the distribution of ${\cal Q}$ is an even but non-Gaussian function [20].

For a two-spin system ($N=2$ and $q=\pm\pi/2$) we have that $\Sigma_{n}(2,\gamma)=1$ for all $n$ and $\gamma$, and when expressions (6.10)-(6.13) are rewritten in dimensionful time $t=\tau/(\nu_{1}+\nu_{2})$ it appears that the cumulants per lattice site $\tfrac{1}{2}\langle{\cal Q}^{n}\rangle_{\rm c}/t$, with ${\cal Q}=\tfrac{1}{2}(Q_{1}-Q_{2})$, are equal to the cumulants $\langle Q_{1}^{n}\rangle_{\rm c}/t$ for a pair in the model considered by Cornu and Bauer [20]. (In other words, in dimensionful time, when $N=2$ the scaled generating function for cumulants per lattice site, $(\nu_{1}+\nu_{2})g_{2}(\lambda)$, is equal to the the scaled cumulant generating function $\alpha(\lambda)$ for the pair of model.) Indeed, in their model where each spin is reversed by only one thermostat, an increment in ${\cal Q}$ is invariant under global flip of the two spins in the initial configuration of a transition, while in the present model where each spin is reversed by both thermostats, an increment in ${\cal Q}$ is invariant under the left-right exchange of the two spins in the $N=2$ chain.

6.3 Large deviation function of the time-integrated current ${\cal Q}/\tau$

The energy ${\cal Q}$ which goes through the system from thermostat 1 to thermostat 2 during a given time $\tau$ is determined by the whole history of the successive changes of spin configurations. We consider the time-integrated current per site and per unit of time (in multiples of $4E$), $\bar{\jmath}$, defined by

[TABLE]

According to definition (6.1) this variable takes the discrete values $\bar{\jmath}_{m}=m/(2N\tau)$, where $m$ is an integer. As time increases the number of discrete values $\bar{\jmath}_{m}$ in a given interval $[\bar{\jmath}-\varepsilon,\bar{\jmath}+\varepsilon]$ (with $\varepsilon>0$) becomes larger and larger. Then the variable $\bar{\jmath}$ is said to satisfy a large deviation principle if there exists a function ${\cal I}_{N}(\bar{\jmath})$ such that 121212For a precise discussion of this definition see Ref.[6] section 5 and Appendix E.

[TABLE]

The limit ${\cal I}_{N}(\bar{\jmath})$ is the so called large deviation function of the current $\bar{\jmath}$. It vanishes for the most probable value of $\bar{\jmath}$, namely when $\bar{\jmath}$ is equal to $\lim_{\tau\to\infty}\langle{\cal Q}\rangle/N\tau$. This value coincides with the mean instantaneous current per site in the stationary state, $\lim_{\tau\to\infty}\langle{\cal Q}\rangle/N\tau=\langle j\rangle$.

One might try to evaluate the large deviation function ${\cal I}_{N}(\bar{\jmath})$ from the definition (6.17) by considering $P({\cal Q};\tau)$ and applying the saddle point method to its inverse Fourier transform representation (6.2)-(6.3) rewritten as an integral on the unit circle by setting $z=\mbox{e}^{{\rm i}p}$. However this method is mathematically tricky because of the singularities in the complex $z$-plane. This is exemplified by the explicit calculation of the leading behavior of $P(\bar{\jmath}_{2};\tau)$ for the time-integrated current $\bar{\jmath}_{2}$ received from thermostat $2$ in the case of a two spin model. We point out that the limit $\tau\to\infty$ must be taken under the condition that $\bar{\jmath}_{2}\tau$ takes only integer values (see section 6 of Ref.[20]).

A far simpler method relies on the Gärtner-Ellis theorem, which ensures that, under weak hypotheses which are fulfilled in the generic case, the expression for ${\cal I}_{N}(\bar{\jmath})$ can be derived from the sole knowledge of the scaled cumulant generating function $g_{N}(\lambda)$ defined in (6.5). For a Markovian process the determination of $g_{N}(\lambda)$ is reduced to the calculation of the largest eigenvalue of the operator that governs the evolution of the generating function $\langle\mbox{e}^{\lambda{\cal Q}}\rangle=\widehat{P}^{\star}(-{\rm i}\lambda/2;\tau)$. In the present case the scaled cumulant generating function $g_{N}(\lambda)$, defined in (6.5), exists and is differentiable for all real $\lambda$, as shown by its expression (6.8)-(6.9). Thus $g_{N}(\lambda)$ satisfies the hypothesis of the simplified version of the Gärtner-Ellis theorem (see, e.g., [37, 38]). This version guarantees that the large deviation function ${\cal I}_{N}(\bar{\jmath})$ of the time-integrated energy current per site exists and can be calculated as the Legendre-Fenchel transform of $g_{N}(\lambda)$, that is,

[TABLE]

Moreover in the present case $g_{N}(\lambda)$ is strictly convex and continuously differentiable for all real $\lambda$. As a consequence the maximum in the definition of the Legendre-Fenchel transform may be calculated with the aid of the Legendre transform,

[TABLE]

where $\lambda_{\bar{\jmath}}$ is the solution of the extremum equation $\mbox{d}g_{N}(\lambda)/\mbox{d}\lambda=\bar{\jmath}$.

In the present case this extremum cannot be solved analytically except for the case $N=2$. Indeed, when $N=2$ only one wave number $q=1/\pi$ is involved in the expression (6.8) for $g_{N}(\lambda)$ and the corresponding expression $g_{2}(\lambda)$ happens to coincide with the scaled cumulant generating function $g(\lambda)$ for another two spin model considered by Cornu and Bauer. Various explicit expressions of ${\cal I}_{2}(\bar{\jmath})$, together with some properties, can be found in section 6.1.2 of Ref.[20].

We point out that ${\cal I}_{N}(\bar{\jmath})$ obeys a generic fluctuation relation which relies on the ratio of transition rates for two reversed jumps of configurations. It can be retrieved for the explicit solution of the paper in various ways. First, since $\ln\sqrt{({\sf A}+{\sf B})/({\sf A}-{\sf B})}=2(\beta_{2}-\beta_{1})E$, the scaled generating function $g_{N}(\lambda)$ given by (6.8)-(6.9) has the symmetry property

[TABLE]

As a consequence ${\cal I}_{N}(\bar{\jmath})$ obeys the fluctuation relation

[TABLE]

This relation also appears for a system of particles moving along a line between two thermostats at different temperatures and endowed with the kinetics of a simple exclusion process [39]. We notice that the symmetry property (6.20) determines the value of the large deviation function for $\bar{\jmath}=0$. Indeed expression (6.18) for the large deviation function together with (6.20) implies that for zero current the minimum is located at the point of symmetry of $g_{N}(\lambda)$, namely $\lambda_{0}=-2(\beta_{2}-\beta_{1})E$. As a consequence {\cal I}_{N}(0)=-g_{N}\big{(}-2(\beta_{2}-\beta_{1})E\big{)}. We notice that, since the system as a finite number of energy levels the long time fluctuation relation (6.21) for ${\cal Q}$ can be derived from the finite time fluctuation relation (5.32) for the couple of variables $Q_{1}$ and $Q_{2}$.

6.4 Infinite size chain at finite effective temperature

When the system size goes to infinity at finite effective temperature ($N\to\infty$ with $\gamma<1$), the limit of the generating function $g_{N}(\lambda;\gamma)$ given by (6.8)-(6.9) reads

[TABLE]

The function $\lim_{N\to\infty}g_{N}(\lambda;\gamma)$, as well as its first derivative with respect to $\lambda$, are well defined for all real values of $\lambda$. Therefore, according to Gärtner-Ellis theorem, when $N$ goes to infinity, there exists a large deviation function ${\cal I}(\bar{\jmath};\gamma)$ given by ${\cal I}(\bar{\jmath};\gamma)=\lim_{N\to\infty}{\cal I}_{N}(\bar{\jmath};\gamma)$.

Moreover, not only the first long time cumulant per site and unit of time $\lim_{\tau\to\infty}\langle{\cal Q}\rangle/N\tau$, but also all other cumulants with $n\geq 2$ remain finite in the limit of infinite size at finite effective temperature. Indeed, when $\gamma<1$, all derivatives of $\lim_{N\to\infty}g_{N}(\lambda;\gamma)$ with respect to $\lambda$ have a finite value at $\lambda=0$ in this limit. The fact that all cumulants per site and unit of time remain finite in this limit can be also retrieved from the structure of the cumulants exhibited by the expressions (6.10) -(6.13) for cumulants of order $n=1,2,3,4$). Indeed the $n$th cumulant is a polynomial of order $n$ in the variables ${\sf A}$ and ${\sf B}$ with coefficients proportional to the $\Sigma_{p}(N,\gamma)$ with $p\leq n$. The finite values of ${\sf A}$ and ${\sf B}$ are independent of $N$ while if $\gamma<1$

[TABLE]

is finite for all $n\geq 1$.

7 Various physical effects

7.1 Kinetic effects

We call ‘kinetic’ those effects that are related to the kinetic parameters $\bar{\nu}_{1}$ and $\bar{\nu}_{2}$ governing the mean frequencies of the spin flips by each thermostat. It is of interest to consider, at arbitrary fixed temperatures $T_{1}$ and $T_{2}$ , the condition $\nu_{2}/\nu_{1}\ll 1$. That is, the colder thermostat flips any spin more slowly than the hotter one. We restore in the discussion below the dimensionful physical time variable $t=\tau/(\nu_{1}+\nu_{2})$. Upon expanding $(\nu_{1}+\nu_{2})g_{N}$ as given by (6.8) and (6.9) in a power series in $\nu_{2}/\nu_{1}$ we find that

[TABLE]

in which

[TABLE]

The argument $\gamma_{1}$ of the function $\Sigma_{N,1}$ in (7.1) is the leading order term of the expansion of $\gamma$ for small $\nu_{2}/\nu_{1}$.

The leading order term in (7.1) is in fact the scaled generating function for the cumulants of a *biased random walk * with step rates $p_{+}\Sigma_{1}(N,\gamma_{1})$ to the right and $p_{-}\Sigma_{1}(N,\gamma_{1})$ to the left, and dimensionful kinetic parameter $\nu_{2}$ [40, 41]. The corresponding formulae in the case where $\nu_{1}\ll\nu_{2}$ are obtained by exchanging $\nu_{1}$ and $\nu_{2}$ and replacing $\gamma_{1}$ by $\gamma_{2}$. In other words, if the indices f and s denote the fast and slow thermostats, respectively, then the scaled generating function given in (7.1) takes the generic form

[TABLE]

From the generic relation (6.6) cumulants read to leading order in $\nu_{\textrm{s}}/\nu_{\textrm{f}}$

[TABLE]

for $m=1,2,\ldots$. The latter expressions, with $t$ in place of $\tau$, can be retrieved from our expressions (6.10)-(6.13) by multiplying them by $\tau/t=\nu_{1}+\nu_{2}$ and expanding them to leading order in $\nu_{\textrm{s}}/\nu_{\textrm{f}}$.

7.2 One thermostat at zero temperature

Dissipation towards a thermostat at zero temperature was studied by Farago and Pitard [18, 19] for an Ising chain in which the energy is injected at a single site. We consider here the corresponding limit for the present model.

Let thermostat 2 have $T_{2}=0$ while we keep $T_{1}>0$. Consequently $\gamma_{2}=1$, which for ${\sf A}$ and ${\sf B}$ given by (1.1) implies that ${\sf A}={\sf B}=1-\gamma_{1}$. Combined with (6.9) this yields ${\theta}(\lambda)=2\bar{\nu}_{1}\bar{\nu}_{2}(1-\gamma_{1})[\mbox{e}^{\lambda}-1].$ When the latter expression is substituted in (6.8), we get that when $\gamma_{2}=1$

[TABLE]

The function $g_{N}(\lambda)$ is now monotonous increasing on the whole real $\lambda$ axis. It follows that the saddle point equation $\mbox{d}g_{N}(\lambda)/\mbox{d}\lambda=\bar{\jmath}$ has no solution for $\bar{\jmath}<0$, which may be restated as

[TABLE]

This expresses the strict impossibility for the energy to flow from the thermostat at $T_{2}=0$ to the one at finite temperature $T_{1}>0$.

The calculation of the cumulants in section 6.2 nevertheless remains valid and their expressions now simplify. The cumulants now become polynomials in $\bar{\nu}_{1}\bar{\nu}_{2}(1-\gamma_{1})$. For instance, the first two cumulants (6.10) and (6.11) now read

[TABLE]

where the inverse temperature $\beta$ in the special case $\gamma_{2}=1$ is determined from $\tanh 2\beta E=\gamma=1-\bar{\nu}_{1}(1-\gamma_{1})$.

7.3 Kinetic effects when colder thermostat is at zero temperature

When the colder thermostat is at zero temperature, $T_{2}=0$, and one thermostat is faster than the other, the scaled generating function given by Eq. (7.1) becomes

[TABLE]

This is the generating function for a Poisson process. As is well known, its cumulants are all equal, and indeed we find, to leading order in $\nu_{\textrm{s}}/\nu_{\textrm{f}}$,

[TABLE]

for $n=1,2,\ldots$. By comparing (7.7) and (7.9) one sees that the limits $T_{2}\to 0$ and $\nu_{\textrm{s}}\ll\nu_{\textrm{f}}$ commute.

8 Large size and low effective temperature

8.1 Parameters at low effective temperature

We now consider the regime where $N\gg 1$ and $0<1-\gamma\ll 1$. According to the relation $\gamma=\gamma_{2}-\nu_{1}(\gamma_{2}-\gamma_{1})$ the condition $0<1-\gamma\ll 1$ corresponds to

[TABLE]

while

[TABLE]

We notice that in the case $\gamma_{1}=\gamma_{2}$ and $0\leq 1-\gamma_{2}\ll 1$ the stationary state would correspond to an equilibrium state at very low temperature.

In view of later analysis we rewrite ${\sf A}$ and ${\sf B}$, defined in (1.1), as

[TABLE]

where ${\sf a}=(1-\gamma_{1}\gamma_{2})/(1-\gamma)$ and ${\sf b}=(\gamma_{2}-\gamma_{1})/(1-\gamma)$. The model is defined for $\bar{\nu}_{1}\bar{\nu}_{2}\neq 0$ and the non-equilibrium condition reads $\gamma_{1}<\gamma_{2}$. As a result the identity $\gamma=\gamma_{2}-\bar{\nu}_{1}(\gamma_{2}-\gamma_{1})$ entails the hierarchy $\gamma_{1}\gamma_{2}\leq\gamma_{1}<\gamma<\gamma_{2}\leq 1$, and $0<{\sf b}\leq{\sf a}<1$.

For the sake of conciseness, from now on we denote the long time cumulants per site and unit of time

[TABLE]

The cumulants can be conveniently split into two contributions : a random walk process with the same first two cumulants as for the ${\cal Q}$ process and a deviation from it. The cumulants $\kappa^{(n)}$ for the random walk are denoted by $\kappa^{(n)}_{\textrm{RW}}$. All even (odd) cumulants take the same value, as exemplified by (7.4) in the case of two thermostats whose kinetic parameters are of different orders of magnitude. The cumulant of order $n$ can be written as

[TABLE]

with the definition

[TABLE]

where ${\sf a}$ and ${\sf b}$ are defined in (8.3). As illustrated by the expressions (6.10)-(6.13) for the first four cumulants, the generic expression of the cumulants $\kappa^{(n)}$ are related to those for the corresponding random walk as follows. The first cumulant of the ${\cal Q}$ process can be reduced to the random wall contribution

[TABLE]

while for $n\geq 2$

[TABLE]

where the deviation $\Delta\kappa^{(n)}(N,\gamma)$ from the random walk process reads

[TABLE]

In (8.9) the factor ${\sf c}_{p}^{(n)}({\sf a},{\sf b})$ is a linear combination of terms ${\sf a}^{q}{\sf b}^{p-q}$, with $q=0,\ldots,p$, where the numerical coefficients depend on the order $n$ of the cumulant; it is determined from the definition (6.6) for every cumulant per site and unit of time and the expression (6.8)-(6.9) for their generating function.

8.2 Finite chain at zero effective temperature

For a finite size chain, the limit of zero effective temperature for the scaled cumulant generating function, $\lim_{\gamma\to 1}g_{N}(\lambda;\gamma)$, is a finite sum given by (6.8)-(6.9) with $\gamma$ equal to one. This function and all its derivative with respect to $\lambda$ are well defined for all real values of $\lambda$. As a consequence, when $\gamma\to 1$ all cumulants are finite and the large deviation function exists and is given by ${\cal I}_{N}(\bar{\jmath};1)=\lim_{\gamma\to 1}{\cal I}_{N}(\bar{\jmath};\gamma)$.

The random walk contribution to the cumulant $\kappa^{(n)}$ is defined in (8.5). According to the explicit expression (6.15) for $\Sigma_{1}(N,\gamma)$, its value at $\gamma=1$ is merely $\Sigma_{1}(N,1)=1$ for all $N$. Therefore in the limit $\gamma\to 1$ the random walk contribution $\kappa^{(n)}_{\textrm{RW}}(N,\gamma)$ vanishes as $1-\gamma$. More precisely

[TABLE]

where ${\sf k}_{n}$ is defined in (8.6).

We now turn to the ${\cal Q}$ process. By virtue of (8.7) its first moment coincides with the first moment of the corresponding random walk, $\kappa^{(1)}(N,\gamma)=\kappa^{(1)}_{\textrm{RW}}(N,\gamma)$, and its leading behavior is the leading behavior of $\kappa^{(1)}_{\textrm{RW}}(N,\gamma)$ given by (8.10). Besides, for all $n\geq 2$ the coefficient $\Sigma_{n}(N,\gamma)$ defined in (6.14), remains finite when $\gamma=1$ at $N$ fixed. Thus, according to the expression (8.9), the deviation $\Delta\kappa^{(n)}(N,\gamma)$ of a cumulant from the corresponding random walk expression vanishes as $(1-\gamma)^{2}$ when $\gamma\to 1$,

[TABLE]

Eventually the leading $(1-\gamma)$-term in the cumulant $\kappa^{(n)}$ is equal to the $(1-\gamma)$-term in the corresponding random walk contribution $\kappa^{(n)}_{\textrm{RW}}$. By virtue of (8.10) it reads

[TABLE]

We point out that the deviation of the first moment $\kappa^{(1)}(N,\gamma)$ from its leading contribution of order $1-\gamma$, denoted by $\Delta\kappa^{(1)}(N,\gamma)$, vanishes as $(1-\gamma)^{2}$, as it is the case for the deviation $\Delta\kappa^{(n)}(N,\gamma)$ of every higher order cumulant from the corresponding random walk cumulant. Indeed, by virtue of (8.7), the difference $\Delta\kappa^{(1)}(N,\gamma)$ is the difference between the first moment of the random walk and its leading $(1-\gamma)$ term, and according to the expression (8.5) for the random walk cumulant, it reads

[TABLE]

It can be rewritten in terms of a single finite sum

[TABLE]

where the increment $\Delta s_{1,\ell}(N,\gamma)$ is written in (A.2). This finite sum indeed converges when $\gamma\to 1$, and

[TABLE]

8.3 Infinite size chain at low effective temperature

In the case of an infinite size chain at finite effective temperature ($\gamma<1$), as discussed in subsection (6.4), the large deviation function exists and all cumulants are finite. When $\gamma\to 1$, the scaled generating function for the cumulants still exists and it is differentiable for all $\lambda$. As a consequence, the large deviation exists and is given by ${\cal I}(\bar{\jmath})=\left.\lim_{N\to\infty}{\cal I}_{N}(\bar{\jmath};\gamma)\right|_{\gamma=1}$, while the first cumulant remains finite.

First we consider the double limit $N\to\infty$ and $\gamma\to 1$ for the random walk process. By virtue of the definition (8.5)

[TABLE]

where the notation for the limit is meant to emphasize the commutativity of the limits $N\to\infty$ and $\gamma\to 1$ for the leading $(1-\gamma)$-term in every random walk cumulant. Indeed, according to the expression (6.15), on the one hand $\lim_{N\to\infty}\Sigma_{1}(N,\gamma)=[1+(\tanh\beta E)^{2}]/2$ and $\lim_{\gamma\to 1}\lim_{N\to\infty}\Sigma_{1}(N,\gamma)=1$, while, on the other hand, for all $N$ $\lim_{\gamma\to 1}\Sigma_{1}(N,\gamma)=1$ and $\lim_{N\to\infty}\lim_{\gamma\to 1}\Sigma_{1}(N,\gamma)=1$.

For the infinite chain (as for the finite chain) the first cumulant $\kappa^{(1)}(\infty,\gamma)$ coincides with the first cumulant of the corresponding random walk $\kappa^{(1)}_{\textrm{RW}}(\infty,\gamma)$ according to (8.7). Therefore the first cumulant in the double limit $N\to\infty$ and $\gamma\to 1$ also vanishes as $(1-\gamma){\sf k}_{1}$. According to the decomposition (8.7)-(8.8), the deviation of $\kappa^{(n)}(N,\gamma)$ from the random walk process, $\Delta\kappa^{(n)}(N,\gamma)$, is a linear combination of terms $(1-\gamma)^{p}\Sigma_{p}(N,\gamma)$ with $2\leq p\leq n$ given by (8.8)-(8.9). For $N$ finite $\Delta\kappa^{(n)}(N,\gamma)$ vanishes as $(1-\gamma)^{2}$ when $\gamma\to 1$. In the double limit $N\to\infty$ and $\gamma\to 1$ every $\Sigma_{p}(N,\gamma)$ with $2\leq p$ diverges according to its definition (6.14) and as can also be seen in the integral representation (6.23) for $\Sigma_{n}(\infty,\gamma)$. Therefore the limits $\gamma\to 1$ and $N\to\infty$ cannot be taken independently of each other for the calculation of $\Delta\kappa^{(n)}$. However $\Delta\kappa^{(n)}(N,\gamma)$ is expected to decay more slowly than $(1-\gamma)^{2}$ but still faster than $1-\gamma$ in an adequate scaling for $N$ and $1-\gamma$. Eventually, in the case of the infinite chain the cumulants vanish as $1-\gamma$, when $\gamma\to 1$ with the same coefficient as the random walk contribution

[TABLE]

Now we turn to the correction to this leading $1-\gamma$ term. In the regime where $N\gg 1$ and $0<1-\gamma\ll 1$, for every cumulant $\kappa^{(n)}(N,\gamma)$ the correction to the leading $(1-\gamma)$-term (8.12) is the sum of two contributions arising from (8.7)-(8.8): on the one hand, the correction (8.13) to the leading $(1-\gamma)$-term in the first moment of the random walk defined in (8.5), and on the other hand, the leading behavior of the deviation $\Delta\kappa^{(n)}$ defined in (8.9). In the double limit $N\gg 1$ and $0<1-\gamma\ll 1$, the sum $\left[1-\gamma\right]^{-1}\left[\Sigma_{1}(N,\gamma)-1\right]$, diverges as well as the $\Sigma_{n}(N,\gamma)$ for $n\geq 2$. Indeed, as detailed in Appendix A, and in the double limit $N\to\infty$ and $\gamma\to 1$ these sums diverge. These divergences can be controlled in two scaling regimes which compare the increasing rates of $N$ and $[\gamma-1]^{-1}$, namely

[TABLE]

Eventually in scaling regime $[\textrm{I}]$, the cumulants behave as

[TABLE]

where $F_{n}^{[\textrm{I}]}(\bar{\nu}_{1}\bar{\nu}_{2};{\sf a},{\sf b})$ does not vanish when $\bar{\nu}_{1}\bar{\nu}_{2}\to 0$, while in scaling regime $[\textrm{II}]$

[TABLE]

where $F_{n}^{[\textrm{II}]}(\rho,\bar{\nu}_{1}\bar{\nu}_{2};{\sf a},{\sf b})$ does not vanish when $\rho=0$ or $\bar{\nu}_{1}\bar{\nu}_{2}\to 0$. We notice that the factor $\bar{\nu}_{1}\bar{\nu}_{2}$ in (8.18) and (8.19) ensures that, when either $\nu_{1}\ll\nu_{2}$ or $\nu_{1}\gg\nu_{2}$ even in the scaling regimes $[\textrm{I}]$ and $[\textrm{II}]$, the leading behavior of the cumulants $\kappa^{(n)}$ per site and unit of time is still given by that of the corresponding random walk defined in (8.5).

8.4 Interpretation of the scaling regimes

Previous results can be interpreted by introducing two physical quantities: the relaxation time to the stationary state and the spin correlation length.

First we recall that in the present model the dynamics for the spin configurations of the system can be seen as a Glauber dynamics with effective kinetic parameter $\nu_{1}+\nu_{2}$ and effective inverse temperature $\beta$. Hence the stationary distribution of spin configurations is the canonical equilibrium distribution at inverse temperature $\beta$, and the evolution of the spin configurations from an initial probability distribution to the stationary one is that of a relaxation to equilibrium. It has been shown by Glauber [10] that in the course of this relaxation the magnetization of the whole chain decays exponentially to its stationary value over the time scale $t_{\textrm{rel}}=[(1-\gamma)(\nu_{1}+\nu_{2})]^{-1}$. In other words, $t_{\textrm{rel}}$ is the relaxation time to the stationary state, a characteristic of which is that the mean magnetization is constant. In the limit $\gamma\to 1$, $\nu_{1}+\nu_{2}$ remains finite and the relaxation time $t_{\textrm{rel}}$ goes to infinity. Therefore it is convenient to consider the long time cumulants per site and per unit of magnetization relaxation time, namely

[TABLE]

which are referred to as “long time rescaled cumulants” in the following. According to (8.17) all rescaled cumulants for the whole chain scale as the system size $N$ in the low-temperature regime,

[TABLE]

where the finite coefficient $k_{n}$ is a random walk cumulant given by (8.6): for $m\geq 1$

[TABLE]

where ${\sf a}$ and ${\sf b}$ are defined in (8.3).

Second, the correlation length $\xi$ is defined from the correlation $\langle s_{n}s_{n+r}\rangle$ between spins at sites $n$ and $n+r$ in the infinite size chain (limit $N\to\infty$ at fixed $\beta$) when the distance $r$ is large. In the present model, the stationary state for spin configurations is the equilibrium state at the effective inverse temperature $\beta$ defined from $\gamma$ by $\gamma=\tanh 2\beta E$. The equilibrium correlation $\langle s_{n}s_{n+r}\rangle$ in the Ising chain with finite size $N$ reads $\langle s_{n}s_{n+r}\rangle=(\zeta^{r}+\zeta^{N-r})/(1+\zeta^{N})$ with $\zeta=\tanh\beta E$. When $N\to\infty$ at fixed $\beta$, it takes the form $\langle s_{n}s_{n+r}\rangle=\zeta^{r}$ at any distance $r$. Therefore the dimensionless correlation length $\xi(\beta)$ in the system is

[TABLE]

In the low effective temperature regime $[\xi(\gamma)]^{-1}\sim 2\mbox{e}^{-2\beta E}\left[1+{\cal O}\left(\mbox{e}^{-2\beta E}\right)\right]$ while $1-\gamma\sim 2\mbox{e}^{-4\beta E}\left[1+{\cal O}\left(e^{-4\beta E}\right)\right]$. Therefore

[TABLE]

In the low effective temperature regime, at leading order all rescaled cumulants for the whole chain scale as the system size $N$ (see (8.21)but the behavior of the subleading term depends on the scaling regime for $N$ and $\beta$.

In scaling regime $[\textrm{I}]$, the size $N$ grows much faster than $\mbox{e}^{2\beta E}$ so that $(1-\gamma)^{2}N\gg 1$. According to (8.24) the latter condition implies that in scaling regime $[\textrm{I}]$, when the temperature decreases the correlation length $\xi(\beta)$ increases but the size $N$ of the chain grows even much faster: $(N/\xi)\gg 1$. Then the scaling behavior (8.18) for the whole chain can be rewritten for $n\geq 1$ as

[TABLE]

where $f_{n}^{[\textrm{I}]}=(1/\sqrt{2})(\bar{\nu}_{1}\bar{\nu}_{2})F_{n}^{[\textrm{I}]}(\bar{\nu}_{1}\bar{\nu}_{2};{\sf a},{\sf b})$. The correlation length $\xi$ may be viewed as the typical size of domains of parallel spins. Thus $N/\xi$ is the typical number of domains with parallel spins or equivalently the number of domain walls, $N_{\textrm{dw}}$,

[TABLE]

Eventually, any rescaled cumulant of the whole chain grows linearly with the number of sites $N$, and in scaling $[\textrm{I}]$ the correction to this leading $N$-behavior scales as the number of parallel spins domains $N_{\textrm{dw}}$.

In scaling regime $[\textrm{II}]$, the size $N$ grows as the temperature goes to zero in such a way that $(1-\gamma)N^{2}=\rho^{2}/2$ with $\rho$ fixed and finite, namely by virtue of (8.24),

[TABLE]

Then the behavior (8.19) of the rescaled cumulants for the whole chain can be rewritten for $n\geq 1$ as

[TABLE]

where $f_{n}^{[\textrm{II}]}(\rho)=(1/\sqrt{2})(\bar{\nu}_{1}\bar{\nu}_{2})F_{n}^{[\textrm{II}]}/(\bar{\nu}_{1}\bar{\nu}_{2};{\sf a},{\sf b})$. In the limit where $\rho\to 0$ the function $f_{n}^{[\textrm{II}]}(\rho)$ goes to a non-vanishing value. The latter regime corresponds to the equilibrium at inverse temperature $\beta$ in the limit of very low temperature. Eventually, in scaling $[\textrm{II}]$ the correction to the leading $N$-behavior of every rescaled cumulant of the whole chain is a finite contribution.

The size dependence of the cumulants of particle currents has been investigated for various exclusion processes: the one-dimensional symmetric simple exclusion process with open boundaries [42], on a ring with periodic boundary conditions [43], for a one-dimensional hard particle gas on a ring or with open boundary conditions [44] and for the one-dimensional lattice gas model $ABC$ in the vicinity of a phase transition [45]. In Ref.[46] the weakly asymmetric exclusion process on a ring has been considered in a scaling regime where the parameter which drives the system out of equilibrium tends to zero as the inverse of the system size; all cumulants of a current are calculated at both leading order and next-to-leading order in the size of the system.

9 Conclusion

The one-dimensional Ising chain has been for a very long time a laboratory for developing the methods of statistical physics. In this work we have contributed to that enterprise. We have considered the $N$-spin cyclic chain with each spin coupled to two thermostats at distinct temperatures $T_{1}$ and $T_{2}$ and a dynamics that generalizes the Glauber [10] model. There appears, as expected, an energy current from the hotter to the colder thermostat. Our fermionization method is a direct extension of the method introduced by Felderhof for the evolution of the spin probability distribution in an Ising chain with Glauber dynamics. It has allowed us to obtain the full spectrum of eigenvalues and eigenvectors of a master equation acting in the product space of the spin configurations and two “counters” that keep track of the net energy furnished by each individual thermostats.

In other words, we have calculated the statistics of the total time-integrated energy current ${\cal Q}$ between the thermostats after a given time interval $\tau$. We found an explicit expression for the probability distribution $P({\cal Q};\tau)$ (at arbitrary finite $N$) at any time $\tau$. In the long time limit we exhibit the generating function for the long time cumulants per site and unit of time $\lim_{\tau\to\infty}\langle{\cal Q}^{n}\rangle_{\rm c}$ for the transferred energy ${\cal Q}$. Their expressions can be determined at any order $n$. We notice that, since the evolution of the joint probability $P(s,{\cal Q};\tau)$ where $s$ is the spin configuration is Markovian, the corresponding generating function is equal to the largest eigenvalue of the matrix that governs the evolution of the Laplace transform of $P(s,{\cal Q};\tau)$ with respect to the variable ${\cal Q}$. Indeed, in models solved by fermionic techniques such as those in Refs.[18, 19, 36], the large deviation of the time-integrated current $X$ of interest is obtained as the largest eigenvalue of that matrix. However in these works the Laplace transform of $P(X;\tau)$, which describes the full statistics, is not exhibited.

The explicit solution for the long time cumulants per site and unit of time has allowed us to investigate effects specific to various regimes of the thermodynamic and kinetic parameters. The main effects are the following. When thermostat 2 is at zero temperature, the current from thermostat 1 to thermostat 2 cannot have negative fluctuations and the large deviation function is non-zero only for positive time-integrated currents : there is pure dissipation towards the zero temperature bath. When one thermostat is very slow with respect to the other one, the generating function for the long time cumulants of ${\cal Q}$ per site and unit of time becomes that of a biased random walk: all odd (even) cumulants are equal to the same value. In this asymmetric random walk the effective kinetic parameter is that of the slower thermostat. This effect has already been exhibited in the two spin model of Ref.[20]. In the present model with $N\geq 2$ spins the sole coefficient due to $N$-body effects that does contribute to the asymmetric random walk cumulants is $\Sigma_{1}(N,\gamma_{\textrm{f}})$, where the index f refers to the slower thermostat: the $N$-body effects involve only the inverse temperature of the faster thermostat. If the colder thermostat is at zero temperature, the generating function for the long time cumulants per site and unit of time becomes that of a Poisson process with an effective kinetic parameter equal to that of the slower thermostat: the random biased walk is confined to positive values of ${\cal Q}$.

In this work we have dealt only with global quantities. However, our results allow for the calculation, in principle, of any quantity related to the energy currents, and in particular energy current-current correlation functions at different points in space and time. This is the subject of ongoing investigation.

Appendix A Behavior of coefficients $\Sigma_{n}(N,\gamma)$

In order to investigate the leading behavior of the correction $\Sigma_{1}(N,\gamma)-1$, where $1=\lim_{N\to\infty,\,\gamma\to 1}\Sigma_{1}(N,\gamma)$, as well as the divergence of $\Sigma_{n}(N,\gamma)$ for $n\geq 2$ in the double limit $N\gg 1$ and $0<1-\gamma\ll 1$, we consider the following finite sums of interest. First we use the property $1=\Sigma_{1}(N,1)$ in order to rewrite the correction $\Sigma_{1}(N,\gamma)-1$ as a single sum

[TABLE]

where

[TABLE]

and the discrete variable $q_{\ell}=(2\ell-1)\pi/N$ varies between $q_{1}=\pi/N$ and $q_{N/2}=\pi[1-1/N]$. Similarly the definition (6.14) can be rewritten as $\Sigma_{n}(N,\gamma)=\sum_{\ell=1}^{N/2}s_{n,\ell}(N,\gamma)$ with $n\geq 2$ and

[TABLE]

In the double limit where $N\to\infty$ and $\gamma\to 1$ the increments defined in (A.2) and (A.3) have the following behavior,

[TABLE]

and

[TABLE]

with the denominator

[TABLE]

At this point one has to distinguish two scaling regimes of parameters.

The scaling regime $[\textrm{I}]$ corresponds to $(1-\gamma)N^{2}\gg 1$. Then we rewrite the denominator $D_{\ell}(N,\gamma)$ as

[TABLE]

with $q_{\ell}^{\star}=(2\ell-1)\pi/[\sqrt{1-\gamma}N]$. Hence, from the definition (A.1), we get that when $(1-\gamma)N^{2}\to\infty$

[TABLE]

In the same limit, for $n\geq 2$ the expression $\sqrt{2(1-\gamma)}^{(2n-3)}\sum_{\ell=1}^{N/2}s_{n,\ell}^{\star}(N,\gamma)$ tends to a constant denoted as $\sigma_{n}^{[\textrm{I}]}$ and

[TABLE]

with

[TABLE]

We notice that, if $N\to\infty$ at $\gamma<1$ fixed, then $\Sigma_{n}(\infty,\gamma)$ is given by (6.23) and it diverges as $1/(\sqrt{1-\gamma})^{2n-3}$ as $\gamma\to 1$ with the same behavior as that given in (A.9). In other words, the result from the successive limits $N\to\infty$ and then $1-\gamma\ll 1$ leads to the same divergence in $1-\gamma$ as if one considers scaling regime $[\textrm{I}]$ where $(1-\gamma)N^{2}\to\infty$.

The scaling regime $[\textrm{II}]$ corresponds to $(1-\gamma)N^{2}=\tfrac{1}{2}\rho^{2}$ with $\rho$ fixed. Then the denominator $D_{\ell}$ defined in (A.6) is conveniently rewritten as

[TABLE]

In the scaling regime $[\textrm{II}]$ the series $[(1-\gamma)N]^{-1}\sum_{\ell=1}^{N/2}\Delta s_{1,\ell}^{\star}(N,\gamma)$ tends a constant denoted as $2\times C_{\textrm{RW}}(\rho)$. Then, from the definition (A.1), we get that

[TABLE]

By virtue of the relation $(1-\gamma)N^{2}=\tfrac{1}{2}\rho^{2}$, the latter behavior can be rewritten as

[TABLE]

with

[TABLE]

In the same scaling $N^{-(2n-3)}\sum_{\ell=1}^{N/2}s_{n,\ell}^{\star}(N,\gamma)$ tends to a constant denoted as $2\times C_{n}(\rho)$,

[TABLE]

By virtue of the relation $(1-\gamma)N^{2}=\tfrac{1}{2}\rho^{2}$, the latter behavior can be rewritten as

[TABLE]

with

[TABLE]

We notice that, if the limit $\gamma\to 1$ is taken at $N$ fixed, then the behavior of $\Sigma_{n}(N,1)$ at large $N$ is given by that of a sum where the $\ell^{\textrm{th}}$ increment $s_{n,\ell}^{\star}(N,1)$ has the denominator $D_{\ell}(N,1)=1/(2N^{2})(2\ell-1)^{2}\pi^{2}$. Then $\Sigma_{n}(N,1))$ behaves as $2N^{2n-3}C_{n}(0)$ where the constant $C_{n}(0)$ happens to be the value of $C_{n}(\rho)$ (A.17) taken at $\rho=0$. In other words, the result from the successive limits $\gamma\to 1$ and then $N\gg 1$ coincides with the behavior (A.15) of $\Sigma_{n}(N,\gamma)$ in scaling $[\textrm{II}]$. In other words, the divergence in $N$ of $\Sigma_{n}(N,\gamma)$ when the limit $\gamma\to 1$ is taken first is the same as in the scaling regime $[\textrm{II}]$ where $\rho=N\sqrt{2(1-\gamma)}$ is fixed and then sent to zero.