\section{Introduction}
\label{secintroduction}

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 \cite{BustamanteETAL2005,Ritort2006}. 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 \textit{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
\footnote{For a comprehensive review  see the report by Seifert Ref.\cite{Seifert2012} and the references therein. In particular, for the case of stochastic Markovian dynamics with jumps between a finite number of configurations see Ref.\cite{HarrisSchutz2007}.}.
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
\footnote{Short introductions which  point out the role of entropy  are  to be found \textit{e.g.} in Refs.\cite{Mallick2009,CornuBauer2013I,VanDenBroeck2013Varenna}.}.
 In particular  the class of systems with a finite number of discrete degrees of freedom has provided firmly established fluctuation relations \cite{HarrisSchutz2007}.

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
\cite{ChouMallickZia2011,Mallick2015}). On the other hand Ising spin chains (with nearest-neighbor ferromagnetic interactions) where all spins are flipped by one of two thermostats.
\vskip 0.3 cm

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 \cite{Glauber1963} 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
\cite{BedeauxEtAl1970,Felderhof1971-A}.
In particular, Felderhof \cite{Felderhof1971-A,Felderhof1971-B} 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 \cite{RaczZia1994,MobiliaZiaSchmittmann2004,LavrentovichZia2010,Lavrentovich2012}. 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 \cite{FaragoPitard2007,FaragoPitard2008}. 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 \cite{CornuBauer2013}. 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
\footnote{In 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 \cite{Cornu2017a}.
}.
The explicit  stationary probability distributions of microscopic configurations have also been obtained  for other archetypal models: the asymmetric exclusion process \cite{Mallick2015} and several variants of the zero-range process \cite{EvansHanney2005,EvansWaclaw2014}. 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\cite{Mallick2015}.
\vskip 0.3cm

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$, $0$, 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.)
\vskip0.3cm

The key to solvability is
the observation
\footnote{This observation goes back at least to Garrido {\it et al.}
\cite{GarridoEtAl1987}, 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 \cite{Felderhof1971-A,Felderhof1971-B} 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)$
\footnote{We 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,  (\ref{FRoneBis})-(\ref{FRoneS}) and (\ref{symmfNz})-(\ref{INfluctrel}), 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,
\begin{equation}
{\sf A} = \bar{\nu}_1\bar{\nu}_2(1-\gamma_1\gamma_2), \qquad
{\sf B} = \bar{\nu}_1\bar{\nu}_2(\gamma_2-\gamma_1),
\label{xABintro}
\end{equation}
where $\gamma_a = \tanh 2\beta_aE$ for $a=1,2$
\footnote{They are the same as the $A$ and $B$ of
Ref.\,\cite{CornuBauer2013}, 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 \cite{CornuBauer2013}
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
\footnote{
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$.}.
\vskip0.3cm

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 (\ref{xABintro}) 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 \ref{secdynIsing} we define the Ising model between two thermostats.
In section \ref{secdirect} we discuss  the instantaneous energy current,
whose average $\langle j\rangle$ per site we determine by elementary means.
In section \ref{secinjected} 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 \ref{secaverages} 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 (\ref{FRoneBis})-(\ref{FRoneS}) which are retrieved from general considerations.
In section \ref{Qstat} 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 \ref{Physicaleffects}
 we study physical effects in  various regimes of the thermodynamic and kinetic parameters for a finite chain.  In  section \ref{secInfiniteSize} we consider a large size chain at very low effective temperature: from the study of some divergent coefficients performed  in Appendix \ref{Details} we exhibit the first correction to the leading $N$-behavior of the cumulants. In section \ref{secconclusion} we briefly conclude.

\section{Ising model coupled to two thermostats}
\label{secdynIsing}

We consider a chain of Ising spins $s_n=\pm1$, 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
\begin{equation}
H(s)=-E\sum_{n=1}^{N} s_ns_{n+1}\,,
\label{defH}
\end{equation}
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
\cite{KadanoffSwift1968} 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
\begin{equation}
|P(\tau)\rangle = \sum_{s}P(s;\tau)|s\rangle.
\label{defketP}
\end{equation}
Since the classical discrete variables $s_n$ all commute, the Ising model has no dynamics of itself.
In 1963 Glauber \cite{Glauber1963}
 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
\begin{equation}
w_n(s;\beta_1)=\frac{1}{2}\bar{\nu}_1[1-\tfrac{1}{2}\gamma_1 s_n(s_{n-1}+s_{n+1})],
\label{defwn}
\end{equation}
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
\begin{equation}
\frac{\mbox{d}}{\mbox{d} \tau} |P(\tau)\rangle = \bar{\nu}_1{M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_1)|P(\tau)\rangle
\label{defmeq}
\end{equation}
with a ``master operator'' ${M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_1)$
whose expression is originally due to Felderhof \cite{Felderhof1971-A,Felderhof1971-B},
\begin{equation}
{M_{\scriptscriptstyle{\text{\mdseries 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],
\label{defMbeta}
\end{equation}
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
\begin{equation}
|P_{\rm eq}(\beta_1)\rangle = \rho_{\rm eq}(\beta_1)|1\rangle,
\qquad |1\rangle\equiv \sum_s|s\rangle,
\label{statstate}
\end{equation}
in which we have
\begin{equation}
\rho_{\rm eq}(\beta_1)= \frac{\mbox{e}^{-\beta_1 {\cal H}}}{Z(\beta_1)}\,, \qquad
{\cal H} =-E\sum_{n=1}^N \sigma_n^z\sigma_{n+1}^z\,, \qquad
Z(\beta_1)=\mbox{Tr}\,\mbox{e}^{-\beta_1 {\cal H}}.
\label{defcalP}
\end{equation}
We remark that $H(s)$ in Eq.\,(\ref{defH}) is an eigenvalue of ${\cal H}$.

By means of fermionization the operator ${M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_1)$ may be completely
diagonalized and all its eigenvectors determined \cite{Felderhof1971-A,Felderhof1971-B}.
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$\,,
\begin{equation}
M=
\bar{\nu}_1{M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_1) + \bar{\nu}_2{M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_2).
\label{defMbis}
\end{equation}
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 (\ref{defMbis}) can be rewritten as
\begin{equation}
M ={M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta),
\label{exM}
\end{equation}
where $\beta$ represents an effective
temperature intermediate between $\beta_1$ and $\beta_2$ given by
\begin{equation}
\tanh 2\beta E = \bar{\nu}_1\tanh 2\beta_1 E + \bar{\nu}_2\tanh 2\beta_2 E.
\label{defbeta}
\end{equation}
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
\footnote{The same observation was made by Cornu and Bauer
\cite{CornuBauer2013} 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.

\section{Energy current between the thermostats}
\label{secdirect}

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_ns_{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
\begin{equation}
f_{\rm al,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],
\qquad N\geq 2.
\label{deffalop}
\end{equation}

We consider the action on this system by the operator $\bar{\nu}_1{M_{\scriptscriptstyle{\text{\mdseries 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
\begin{equation}
w_{\rm al,op}(\beta_1) = \tfrac{1}{2}\bar{\nu}_1(1 \mp \tanh 2\beta_1 E),
\label{defwalop}
\end{equation}
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
\begin{eqnarray}
\langle j_1\rangle &=& f_{\rm al}\,w_{\rm al}(\beta_1)
             - f_{\rm op}\,w_{\rm op}(\beta_1) \nonumber\\[2mm]
&& = \frac{\bar{\nu}_1}{2}
\left[ \frac{\zeta+\zeta^{N-1}}{1+\zeta^N}
 - \tfrac{1}{2}(1+\zeta^2)\frac{1+\zeta^{N-2}}{1+\zeta^N}
\,\tanh 2\beta_1 E \right].
\label{xj1av}
\end{eqnarray}
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{\mdseries th}}}}(\beta_2)$. From (\ref{defbeta}) and (\ref{xj1av}) 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
\begin{equation}
\langle j \rangle = \tfrac{1}{4}\bar{\nu}_1\bar{\nu}_2
(1+\zeta^2)\frac{1+\zeta^{N-2}}{1+\zeta^N}
\big[ \tanh2\beta_2 E - \tanh2\beta_1 E \big].
\label{resjav}
\end{equation}
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
\begin{equation}
\langle {\cal Q} \rangle \simeq \langle j\rangle N\tau, \qquad \tau\to\infty,
\label{Javdirect}
\end{equation}
and $\langle j\rangle$ given by (\ref{resjav}).
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 (\ref{resjav}) 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_BT^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 (\ref{defbeta}) we then have
\begin{equation}
\delta\beta_1=\bar{\nu}_2\delta\beta_{12}, \qquad
\delta\beta_2=-\bar{\nu}_1\delta\beta_{12}.
\label{xdeltabetai}
\end{equation}
Calling the linearized current $\delta j$, we obtain from (\ref{resjav})
\begin{equation}
\langle \delta j \rangle = \lambda_{\rm T}\delta T, \qquad
\lambda_{\rm 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_{\rm B}.
\label{resjavlin}
\end{equation}
The heat conduction coefficient $\lambda_{\rm T}$
tends to zero in both limits $\beta\to 0$  and $\beta\to\infty$, with $E$ fixed.

\section{Extended master operator: definition and diagonalization}
\label{secinjected}

\subsection{Extended master operator $\widehat{\cal M}$}
\label{secextendedME}

Each spin reversal is due to either
${M_{\scriptscriptstyle{\text{\mdseries th}}}}(\beta_1)$ or ${M_{\scriptscriptstyle{\text{\mdseries 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{\mdseries th}}}}(\beta_1)$
and ${M_{\scriptscriptstyle{\text{\mdseries 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)
\label{sumAP}
$
and
the initial condition
\begin{equation}
 P(s,\vec{Q}; 0)=\delta_{\vec{Q},\vec0} \, P(s; 0).
\end{equation}
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
 \begin{eqnarray}
\frac{ \mbox{d} P(s,\vec{Q}; \tau)}{\mbox{d} \tau}&=& - \left[\sum_{a=1,2} \sum_{n=1}^N w_n(s;\beta_a)\right] P(s,\vec{Q}; \tau)
\nonumber
\\ &+&\sum_{n=1}^N w_n(s^{n};\beta_1) P(s^{n}, Q_1+ \Delta Q_n(s), Q_2; \tau)
\nonumber
\\ &+&  \sum_{n=1}^N w_n(s^{n};\beta_2) P(s^{n}, Q_1, Q_2 + \Delta Q_n(s); \tau)
\label{BalanceP}
\end{eqnarray}
By analogy with the representation (\ref{defketP}) of $P(s;\tau)$, we
represent the probability $P(s, \vec{Q};\tau)$ by the time dependent ket
\begin{equation}
|P(\vec{Q}; \tau)\rangle = \sum_s P(s,\vec{Q};\tau) | s\rangle.
\label{defketPbis}
\end{equation}
We consider the Fourier transformed ket
\begin{equation}
|\widehat{P}(\vec{p}; \tau)\rangle =
\sum_{\vec{Q}}\,\, \mbox{e}^{{\rm i} \vec{p}\cdot \vec{Q}} | P(\vec{Q};\tau)\rangle,
\label{defketPF}
\end{equation}
where $\vec{p}=(p_1,p_2)$ with $-\pi<p_1,p_2\leq\pi$.
Upon taking the Fourier transform of
the balance equation
(\ref{BalanceP}) we get the evolution equation for the ket (\ref{defketPF}),
\begin{equation}
\frac{ \mbox{d} |\widehat{P}(\vec{p}; \tau)\rangle}{\mbox{d} \tau}= \widehat{\cal M}(\vec{p}) |\widehat{P}(\vec{p}; \tau)\rangle,
\label{evolTFP}
\end{equation}
in which
\begin{equation}
\widehat{\cal M}(\vec{p})= \bar{\nu}_1 \widehat{{\cal M}}_{\rm th}(p_1;\beta_1)+ \bar{\nu}_2
\widehat{{\cal M}}_{\rm th}(p_2;\beta_2)
\label{whcalMsum}
\end{equation}
where, by analogy with (\ref{defMbeta}),
\begin{equation}
\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].
\label{defTFMbeta}
\end{equation}
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$, $0$, 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 (\ref{whcalMsum}) may be rewritten as
\begin{eqnarray}
\widehat{\cal M}(\vec{p}) &=& \tfrac{1}{2}\sum_{n=1}^N
\Big[ \tfrac{1}{2}(1+C)\sigma_n^x
-\tfrac{1}{2}D\,\sigma_n^x\sigma_n^z(\sigma_{n-1}^z+\sigma_{n+1}^z)
\nonumber\\
&& \phantom{XXX}-\tfrac{1}{2}(1-C)\sigma_n^x\sigma_{n-1}^z\sigma_{n+1}^z
-1
+\tfrac{1}{2}\gamma\sigma_n^z(\sigma_{n-1}^z+\sigma_{n+1}^z) \Big],
\phantom{xxx}
\label{exMp}
\end{eqnarray}
in which
\begin{eqnarray}
C(\vec{p})&=&\bar{\nu}_1 [\cos p_1 -{\rm i}\gamma_1\sin p_1]+\bar{\nu}_2[\cos p_2 -{\rm i}\gamma_2\sin p_2],
\nonumber\\
D(\vec{p})&=& \bar{\nu}_1 [\gamma_1\cos p_1 -{\rm i}\sin p_1]+\bar{\nu}_2[\gamma_2\cos p_2 -{\rm i}\sin p_2]\,.
\label{exAB}
\end{eqnarray}
These coefficients are real when $p_1$ and $p_2$ are pure imaginary.

\subsection{Symmetrizing the master operator}
\label{secsymm}

We apply to $\widehat{\cal M}(\vec{p})$ a similarity transformation
and define
\begin{equation}
\widetilde{\cal M}(\vec{p})
= \rho_{\rm eq}^{-\frac{1}{2}}({\beta_*})\widehat{\cal M}(\vec{p})\,
\rho_{\rm eq}^{\frac{1}{2}}({\beta_*}),
\label{dPsymm}
\end{equation}
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 (\ref{dPsymm}) is \cite{Felderhof1971-A,Felderhof1971-B}
\begin{eqnarray}
\tilde{\sigma}_n^x ({\beta_*})&\equiv&
\rho_{\rm eq}^{-\frac{1}{2}}({\beta_*}) \sigma_n^x
\rho_{\rm eq}^{ \frac{1}{2}}({\beta_*})
\nonumber\\[2mm]
&=&\sigma_n^x\Big[\cosh^2{\beta_*} E +
\sigma_{n-1}^z\sigma_{n+1}^z\sinh^2{\beta_*} E
\nonumber\\[2mm]
&& \phantom{\big[}
+ \sigma_n^z(\sigma_{n-1}^z+\sigma_{n+1}^z)
\sinh{\beta_*} E\cosh{\beta_*} E \Big],
\label{deftsigma}
\end{eqnarray}
which is easily derived.
The result is that $\widehat{\cal M}(\vec{p})$ of equation (\ref{exMp})
becomes an expression $\widetilde{\cal M}(\vec{p})$ which is of the same form
as (\ref{exMp}) but with $C$ and $D$ of equation (\ref{exAB})
replaced with $\widetilde{C}$ and $\widetilde{D}$,
respectively, where
\begin{eqnarray}
\widetilde{C}(\vec{p},{\beta_*}) &=&
C(\vec{p})\cosh 2{\beta_*} E - D(\vec{p})\sinh 2{\beta_*} E,
\nonumber\\[2mm]
\widetilde{D}(\vec{p},{\beta_*}) &=&
C(\vec{p})\sinh 2{\beta_*} E - D(\vec{p})\cosh 2{\beta_*} E.
\label{extAtB}
\end{eqnarray}
We now choose ${\beta_*}$ such that the coefficient $\widetilde{D}({\beta_*})$
of the non-Hermitian term vanishes.
This amounts to taking
\begin{equation}
\tanh 2{\beta_*}(\vec{p})E = \frac{D(\vec{p})}{C(\vec{p})}
\label{eqbetastar}
\end{equation}
where $C(\vec{p})$ and $D(\vec{p})$ are given by (\ref{exAB}).
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
\begin{equation}
\widetilde{\cal M}(\vec{p})=\tfrac{1}{2}\sum_{n=1}^N \Big[
\tfrac{1}{2}(1+C_*)\sigma_n^x
-\tfrac{1}{2}(1-C_*)\sigma_n^x\sigma_{n-1}^z\sigma_{n+1}^z
-1+\tfrac{1}{2}\gamma(\sigma_{n-1}^z\sigma_n^z+\sigma_n^z\sigma_{n+1}^z)
\Big]
\label{Mpsigma}
\end{equation}
in which $C_*$ is given by
\begin{equation}
C_*(\vec{p}) \equiv \widetilde{C}(\vec{p},{\beta_*}(\vec{p})).
\label{exAstarp}
\end{equation}
All $\vec{p}$ dependence of $\widehat{\cal M}(\vec{p})$ is seen to enter through
the single coefficient $C_*(\vec{p})$.

After substituting (\ref{exAB}) in (\ref{extAtB})
and (\ref{extAtB}) in (\ref{exAstarp})
we find that this quantity may be written as
\begin{equation}
C_*^2(\vec{p})=1-\gamma^2 + \Theta(\vec{p})
\label{xAstar}
\end{equation}
where
\begin{equation}
\Theta(\vec{p})=2{\sf A}\left[\cos(p_1-p_2) -1 \right]
+ 2{\rm i}{\sf B} \sin(p_1-p_2)
\label{defTheta}
\end{equation}
with ${\sf A}$ and ${\sf B}$ given by (\ref{xABintro}) in the Introduction.
These coefficients
will appear again in our final
results in section \ref{Qstat}.

\subsection{Transformation to fermion operators}
\label{secfermion}

We define fermionic quasi-particles by means of the Jordan-Wigner \cite{JordanWigner1928} transformation
\begin{eqnarray}
c_n^\dagger &=& \tfrac{1}{2}
\left[ \,\prod_{j=1}^{n-1}\sigma_j^x \right](\sigma_n^z+{\rm i}\sigma_n^y),
\nonumber\\[2mm]
c_n         &=& \tfrac{1}{2}
\left[ \,\prod_{j=1}^{n-1}\sigma_j^x \right](\sigma_n^z-{\rm i}\sigma_n^y),
\qquad n=1,2,\ldots,N.
\label{defcdagc}
\end{eqnarray}
The vacuum state of these $c$-particles is the state $|1\rangle$ defined in
(\ref{statstate}).
It is now straightforward to express $\widetilde{\cal M}(\vec{p})$ in terms of these
fermion operators. We find from (\ref{Mpsigma})
\begin{eqnarray}
\widetilde{{\cal M}}(\vec{p}) &=& -\tfrac{1}{4}N(1-C_*)
            -\tfrac{1}{2}(1+C_*)\sum_{n=1}^N c^\dag_n c_n \nonumber\\
&& +\tfrac{1}{2}\gamma\sum_{n=1}^N(c^\dag_n-c_n)(c^\dag_{n+1}-c_{n+1})
\nonumber\\
&&
-\tfrac{1}{4}(1-C_*)\sum_{n=1}^N(c^\dag_n-c_n)(c^\dag_{n+2}-c_{n+2})
\label{xMc}
\end{eqnarray}
with the understanding that
the creation and annihilation operators whose indices exceed $N$ are
defined by
\begin{eqnarray}
c^\dag_{N+m} &=& -c^\dag_m(-1)^{{\cal N}}, \nonumber\\[2mm]
c_{N+m} &=& -c_m(-1)^{{\cal N}}, \qquad m=1,2,
\label{defcNpm}
\end{eqnarray}
in which ${\cal N}=\sum_{n=1}^N c^\dag_n c_n$ is the operator for the total number of
quasi-particles.

\subsection{Diagonalizing in terms of fermion operators}

For convenience we hence restrict ourselves to even $N$.
We define fermion operators $\eta^\dag_q$ and $\eta_q$ by
\begin{eqnarray}
c^\dag_n &=& N^{-1/2}\sum_q \mbox{e}^{-{\rm i}qn} \eta^\dag_q\,,  \nonumber\\[2mm]
c_n      &=& N^{-1/2}\sum_q \mbox{e}^{ {\rm i}qn} \eta_q\,, \qquad n=1,...,N,
\label{defetaq}
\end{eqnarray}
where the wavenumber $q$ runs through the $N$ values
\begin{equation}
q=\pm\frac{\pi}{N},\pm\frac{3\pi}{N},...,\pm\frac{(N-1)\pi}{N}\,.
\label{defrangeq}
\end{equation}
Equation (\ref{defetaq}) is easily inverted to find the $\eta^\dag_q$ and $\eta_q$
in terms of the $c^\dag_n$ and $c_n$.
This equation guarantees the periodicity conditions (\ref{defcNpm})
in the subspace where ${\cal N}$ is {\it even}. In that subspace
equation (\ref{defetaq})
may also be used in (\ref{xMc}) for $n=N+1$ and $n=N+2$.
Obviously the $c$ vacuum $|1\rangle$
is also the $\eta$ vacuum.

Applying transformation (\ref{defetaq}) to (\ref{xMc}) we get
\begin{equation}
\widetilde{\cal M}(\vec{p}) = -\tfrac{1}{2}N -\tfrac{1}{2}\sum_q
\left[C_q(\etadqq + \etadmqmq - 1) -{\rm i}D_q(\etadqdmq + \etaqmq) \right],
\label{xMeta}
\end{equation}
valid in the subspace with an even number ${\cal N}$ of $c$-particles.
\footnote{In the subspace with an {\it odd\,} number of $c$-particles
  $\widetilde{\cal M}(\vec{p})$ takes a slightly different form, as discussed in detail in
  references \cite{Felderhof1971-A,Felderhof1971-B,HilhorstETAL1972,Hilhorst1973,Hilhorst1974}
We will not need that form in this work.},
and where the coefficients $C_q$ and $D_q$ are given by
\begin{eqnarray}
C_q(\vec{p}) &=& \tfrac{1}{2}(1+C_*) -\gamma\cos q +\tfrac{1}{2}(1-C_*)\cos 2q,
\nonumber\\[2mm]
D_q(\vec{p}) &=& \gamma\sin q -\tfrac{1}{2}(1-C_*)\sin 2q,
\label{defCD}
\end{eqnarray}
where the $\vec{p}$ dependence comes in through the $q$ independent coefficient
$C_*(\vec{p})$ defined in (\ref{exAstarp}).
Extending the approach of Ref.\,\cite{Felderhof1971-A,Felderhof1971-B}
to nonzero $\vec{p}$ we define angles $\chi_q$ (that are in general complex) by
\begin{equation}
\cos\chi_q(\vec{p}) = \frac{C_q}{\sqrt{C_q^2+D_q^2}}\,, \qquad
\sin\chi_q(\vec{p}) = \frac{D_q}{\sqrt{C_q^2+D_q^2}}\,,
\label{defchi}
\end{equation}
and perform in the space of the pair
$\{\eta_q,\eta^\dag_{-q}\}$  a Bogoliubov-Valatin
\cite{Bogoliubov1958,Valatin1958} operator rotation
\begin{eqnarray}
\xi_q(\vec{p}) &=& \phantom{-}
(\cos\tfrac{1}{2}\chi_q)\,\eta_q - {\rm i}(\sin\tfrac{1}{2}\chi_q)\,\eta^\dag_{-q}\,,
\nonumber\\[2mm]
\xi^\dag_{-q}(\vec{p}) &=&
- {\rm i}(\sin\tfrac{1}{2}\chi_q)\,\eta_q + (\cos\tfrac{1}{2}\chi_q)\,\eta^\dag_{-q}\,.
\label{defxi}
\end{eqnarray}
It is useful to note that $\chi_{-q}=-\chi_q$. Upon using (\ref{defxi})
to transform (\ref{xMeta}) to $\xi$ operators
we find the diagonal form
\begin{equation}
\widetilde{\cal M}(\vec{p}) = -\mu_* -\sum_q\mu_q\xi^\dag_q\xi_q
\label{x1Mxi}
\end{equation}
where
\begin{equation}
\mu_q(\vec{p}) = \sqrt{C_q^2+D_q^2}
\label{defmu}
\end{equation}
and
\begin{equation}
\mu_*(\vec{p}) = \frac{1}{2}\sum_{q}(1-\mu_q).
\label{defmustar}
\end{equation}
From (\ref{defCD}) and (\ref{defmu}) it is easily seen that $\mu_q=\mu_{-q}$.
Upon combining both equations we get for $\mu_q$ the
explicit expression
\begin{equation}
\mu_q^2 = (\gamma-\cos q)^2 + C_*^2(\vec{p})\sin^2 q,
\label{x1muq}
\end{equation}
with $C_*(\vec{p})$ given by (\ref{xAstar}).
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.\,\cite{Felderhof1971-A,Felderhof1971-B}, namely $\mu_q(\vec{0})=1-\gamma\cos q$.

A different way, useful for later, to write the eigenvalue $\mu_q(\vec{p})$ is
\begin{equation}
\mu_q^2 = (1-\gamma\cos q)^2 + \Theta(\vec{p})\sin^2 q,
\label{x2muq}
\end{equation}
where $\Theta$ has been defined in (\ref{defTheta}).
We observe for later use that
\begin{equation}
\mu_*(\vec{0})=0, \qquad \Theta(\vec{0})=0.
\label{relvecpzero}
\end{equation}
It is convenient to rewrite the diagonalized form
(\ref{x1Mxi}) of the master operator as
\begin{equation}
\widetilde{\cal M}(\vec{p}) = -\tfrac{1}{2} N - \sum_{q>0} \mu_q (\xidqq + \xidmqmq - 1),
\label{x2Mxi}
\end{equation}
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 (\ref{defrangeq}).

\section{Joint probability distribution of the  time-integrated energy currents}
\label{secaverages}

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

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 (\ref{defketPbis})
this probability distribution is
given by $P(\vec{Q};\tau)=\sum_s \, \langle s |P(\vec{Q};\tau)\rangle$  and upon inverting (\ref{defketPF}) we find
\begin{equation}
P(\vec{Q};\tau) = \int_{-\pi}^{\pi}\frac{\mbox{d} p_1}{2\pi}\int_{-\pi}^{\pi}\frac{\mbox{d} p_2}{2\pi}\,\mbox{e}^{-{\rm i}\vec{p}\cdot\vec{Q}}\,\,\widehat{P}(\vec{p};\tau),
\label{PInverseTF}
\end{equation}
where $\widehat{P}(\vec{p};\tau)=\sum_s \langle s |\widehat{P}(\vec{p};\tau)\rangle$.
The evolution equation (\ref{evolTFP}) may be formally solved as
\begin{equation}
|\widehat{P}(\vec{p};\tau)\rangle = \mbox{e}^{{\widehat{\cal M}(\vec{p})} \tau} |\widehat{P}(\vec{p};0)\rangle,
\end{equation}
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},\vec0} |P_{\rm eq}(\beta_0))$. With the definitions
(\ref{statstate}) the  probability $\widehat{P}(\vec{p};\tau)$  reads
\begin{equation}
\widehat{P}(\vec{p};\tau)=\langle 1|\mbox{e}^{{\widehat{\cal M}(\vec{p})} \tau} \rho_{\rm eq}(\beta_0)|1\rangle\,.
\label{PInverseTFbis}
\end{equation}
This expression takes advantage of the fact that
${\widehat{\cal M}}$ is block diagonal in the subspaces of fixed $\vec{p}$.

\subsection{Rewriting $P(\vec{Q};\tau)$}
\label{secrewrPi}

According to the complete diagonalization performed in section \ref{secinjected}
 the matrix element in the Fourier transform (\ref{PInverseTFbis}) is an expectation value in the $\eta$ vacuum and
 can be rewritten as
\begin{equation}
P(\vec{p};\tau)=
\langle 1|\rho_{\rm eq}^{1/2}({\beta_*})\,\mbox{e}^{\widetilde{\cal M}(\vec{p})\tau}
\rho_{\rm eq}^{-1/2}({\beta_*})\rho_{\rm eq}(\beta_0)|1\rangle,
\label{x1hPi}
\end{equation}
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
\begin{equation}
|1\rangle = 2^{N/2}\otimes_{q>0}|0_q0_{-q}\rangle,
\label{defqvac}
\end{equation}
where $|0_q0_{-q}\rangle$ is the state in which the
quasi-particles of wavenumbers $\pm q$ are absent,
\begin{equation}
\eta_q|0_q0_{-q}\rangle=0, \qquad \eta_{-q}|0_q0_{-q}\rangle=0,
\label{relqvac}
\end{equation}
and $\langle 0_q0_{-q}|0_q0_{-q}\rangle=1$.

It is useful to rewrite the time evolution operator (\ref{x2Mxi}) as
\begin{equation}
{\widetilde{\cal M}}(\vec{p}) = -\tfrac{1}{2} N - \sum_{q>0} \mu_{q}\mathbb{X}_q
\label{Mrewrite}
\end{equation}
with
\begin{equation}
\mathbb{X}_q = \xidqq + \xidmqmq - 1,
\label{defX}
\end{equation}
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
\begin{eqnarray}
{\cal H} &=& -2E\sum_{q>0} \mathbb{H}_q\,,
\nonumber\\[2mm]
\mathbb{H}_q &=& \mathbb{A}_q\cos q + {\rm i}\mathbb{B}_q \sin q,
\label{Hrewrite}
\end{eqnarray}
where
\begin{eqnarray}
\mathbb{A}_q &=& \etadqq + \etadmqmq - 1,  \nonumber\\[2mm]
\mathbb{B}_q &=& \etadqdmq + \etaqmq\,,  \nonumber\\[2mm]
\mathbb{D}_q &=& \etadqq\etadmqmq\,,
\label{defAB}
\end{eqnarray}
where we included $\mathbb{D}_q$ for later reference.
We now use the fact that $\mathbb{X}_q$
[in view of relations (\ref{defX}) and (\ref{defxi})]
and $\mathbb{H}_q$ [in view of (\ref{Hrewrite}) and (\ref{defAB})]
are quadratic in the $\eta$ operators and that therefore
\begin{equation}
[\mathbb{X}_q,\mathbb{X}_{q^\prime}]=
[\mathbb{X}_q,\mathbb{H}_{q^\prime}]=
[\mathbb{H}_q,\mathbb{H}_{q^\prime}]=0,
\qquad q\neq q^\prime.
\label{relcommut}
\end{equation}
Upon using (\ref{defqvac}) in (\ref{x1hPi})
we may factorize  $\widehat{P}(\vec{p};\tau)$ according to
\begin{equation}
{\widehat{P}}(\vec{p};\tau)= \frac{2^N}{Z(\beta_0)}\,\mbox{e}^{-\frac{1}{2}N\tau}
\prod_{q>0}\Pi_q(\vec{p};\tau)
\label{factPi}
\end{equation}
in which
\begin{equation}
\Pi_q(\vec{p};\tau)= \langle 0_q0_{-q}|\mbox{e}^{{\beta_*} E\,\mathbb{H}_q}
\,\mbox{e}^{-\mu_q\mathbb{X}_q\tau}
\,\mbox{e}^{(2\beta_0-{\beta_*})E\,\mathbb{H}_q}|0_q0_{-q}\rangle,
\label{defMqp}
\end{equation}
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_q0_{-q}\rangle$
defined above and the two-particle state $|1_q1_{-q}\rangle$ defined by
\begin{equation}
|1_q1_{-q}\rangle=\eta^\dag_q\eta^\dag_{-q}|0_q0_{-q}\rangle.
\label{def1q}
\end{equation}

To make the action of $\exp(K\mathbb{H}_q)$ more explicit we expand the
exponential using the relations
\begin{equation}
\mathbb{H}_q^2 = -\mathbb{A}_q+2\mathbb{D}_q, \qquad
\mathbb{H}_q^3 = \mathbb{H}_q, \qquad  \mathbb{H}_q^4=\mathbb{H}_q^2,
\label{relH}
\end{equation}
which are easily checked.
One then obtains
\begin{eqnarray}
\mbox{e}^{K\mathbb{H}_q} &=& 1 + \mathbb{H}_q\sinh K + \mathbb{H}_q^2(\cosh K -1)
\nonumber\\[2mm]
&=& d_0(K) + d_1(K)(\etadqq + \etadmqmq) + d_2(K)(\etadqdmq + \etaqmq) \nonumber\\[2mm]
&& + d_4(K)\etadqq\etadmqmq
\label{expH}
\end{eqnarray}
in which
\begin{eqnarray}
d_0(K) &=& \cosh K - \cos q\sinh K,      \nonumber\\[2mm]
d_1(K) &=& 1 - \cosh K + \cos q \sinh K, \nonumber\\[2mm]
d_2(K) &=& {\rm i}\sin q \sinh K,        \nonumber\\[2mm]
d_4(K) &=& 2(\cosh K - 1).
\label{defdi}
\end{eqnarray}
For two specific choices of $K$ we will use below the notation
\begin{equation}
b_i=d_i([2\beta_0-{\beta_*}]E), \qquad c_i=d_i({\beta_*} E), \qquad i=0,1,2,4.
\label{defbc}
\end{equation}
We have to similarly expand $\exp(-\mu_q\mathbb{X}_q\tau)$ and obtain along the
same lines
\begin{eqnarray}
\mbox{e}^{-\mu_q\mathbb{X}_q\tau} &=& 1 -\mathbb{X}_q\sinh\mu_q\tau
+ \mathbb{X}_q^2(\cosh\mu_q\tau - 1)
\nonumber\\[2mm]
&=&  \mbox{e}^{\mu_q\tau} - (\mbox{e}^{\mu_q\tau}-1)(\xidqq + \xidmqmq)
+ 2(\cosh\mu_q\tau - 1)\xidqq\xidmqmq.
\nonumber\\
&&
\label{exp1M}
\end{eqnarray}
With the aid of the relations (\ref{defxi})
we can turn this into an expansion of the form
\begin{eqnarray}
\mbox{e}^{-\mu_q\mathbb{X}_q\tau} &=& a_0 + a_1(\etadqq + \etadmqmq) \nonumber\\[2mm]
&& + a_2(\etadqdmq + \etaqmq) + a_4\etadqq\etadmqmq.
\label{expM}
\end{eqnarray}
After a fair amount of algebra
one finds for the coefficients $a_i$ the expressions
\begin{eqnarray}
a_0 &=&     \cosh\mu_qt + \cos\chi_q \sinh\mu_qt,
\nonumber\\[2mm]
a_1 &=& 1 - \cosh\mu_qt - \cos\chi_q \sinh\mu_qt, \nonumber\\[2mm]
a_2 &=& {\rm i}\sin\chi_q \sinh\mu_qt, \nonumber\\[2mm]
a_4 &=& 2(\cosh\mu_qt - 1).
\label{defai}
\end{eqnarray}
We substitute now expansions (\ref{expH})
for $K={\beta_*} E$ or $K=(2\beta_0-{\beta_*})E$ and (\ref{expM}) in
(\ref{defMqp}).
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 (\ref{defMqp}) as
\begin{eqnarray}
\Pi_q(\vec{p};\tau) &=& \la0_q0_{-q}|[c_0+c_2\eta_q\eta_{-q}] \nonumber\\[2mm]
&& \times[a_0+a_1(\etadqq + \etadmqmq)+a_2(\etadqdmq + \etaqmq)+a_4\etadqq\etadmqmq] \nonumber\\[2mm]
&& \times[b_0+b_2\eta^\dag_q\eta^\dag_{-q}]|0_q0_{-q}\rangle \nonumber\\[2mm]
&=& a_0b_0c_0 - a_2(b_0c_2+b_2c_0) - (a_0+2a_1+a_4)b_2c_2.
\label{x1Mpq}
\end{eqnarray}
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 (\ref{x1Mpq})
we may cast the $\Pi_q(\vec{p};\tau)$ in the form
\begin{equation}
\Pi_q(\vec{p};\tau) = S_q(\beta_0)\left[\cosh\mu_q(\vec{p})  \tau +  \frac{ T_q(\vec{p};\beta_0)}{S_q(\beta_0)}\frac{\sinh\mu_q(\vec{p})\tau}{\mu_q(\vec{p})} \right],
\label{x2Mpq}
\end{equation}
where $\mu_q$ is given by (\ref{x2muq}) while
\begin{eqnarray}
S_q (\beta_0) &=& \cosh(2\beta_0E) -  \sinh(2\beta_0E) \cos q\,, \nonumber\\[2mm]
T_q(\vec{p}; \beta_0) &=& S_q (\beta_0)(1-\gamma \cos q )+ U(\vec{p};\beta_0) \sin^2 q   \,
\label{defSTUV}
\end{eqnarray}
with
\begin{equation}
U(\vec{p};\beta_0)=\bar{\nu}_1 \, u(p_1; \beta_0-\beta_1) + \bar{\nu}_2 \,u(p_2; \beta_0-\beta_2)
\end{equation}
and
\begin{equation}
u(p_a; \beta)=\cosh(2\beta E + {\rm i} p_a) -\cosh(2\beta E).
\end{equation}
Combining (\ref{factPi}) and (\ref{x2Mpq}) we get
\begin{equation}
\widehat{P}(\vec{p};\tau)=\frac{2^N}{Z(\beta_0)}\,\mbox{e}^{-\frac{1}{2}N\tau}\,
\left[ \prod_{q>0}S_q(\beta_0) \right]
\prod_{q>0}\left[ \cosh\mu_q \tau + \frac{T_q(\vec{p};\beta_0)}{S_q(\beta_0) \mu_q(\vec{p})}\sinh\mu_q \tau\right].
\label{x2hPi}
\end{equation}
where the partition function defined in (\ref{defcalP}) reads
\begin{equation}
Z(\beta_0)=2^N \left[ \cosh^N \beta_0 E + \sinh^N \beta_0 E \right].
\end{equation}
It is easy to verify the relation
\begin{equation}
\frac{2^N}{Z(\beta_0)}\,\prod_{q>0}S_q (\beta_0)= 1.
\label{normcond}
\end{equation}
Using it in (\ref{x2hPi}) and substituting (\ref{x2hPi}) in
(\ref{PInverseTF}) we finally obtain
\begin{eqnarray}
P(\vec{Q};\tau) &=& \mbox{e}^{-\frac{1}{2}N\tau}\int_{-\pi}^{\pi}\frac{\mbox{d} p_1}{2\pi}\int_{-\pi}^{\pi}\frac{\mbox{d} p_2}{2\pi}\,
\mbox{e}^{-{\rm i}\vec{p}\cdot\vec{Q}}\,
\prod_{q>0}\left[ \cosh\mu_q \tau + \frac{T_q(\vec{p};\beta_0)}{S_q(\beta_0) \mu_q(\vec{p})}\sinh\mu_q \tau\right].
\nonumber\\
&&
\label{resfinPi}
\end{eqnarray}
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

\begin{equation}
\frac{T_q(\vec{0};\beta_0)}{\mu_q(\vec{0})} = S_q(\beta_0),
\end{equation}
which together with (\ref{x2hPi}), (\ref{normcond}), and (\ref{relvecpzero})
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$.

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

One can check  on the explicit expression (\ref{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
\begin{equation}
\frac{P(Q_1,Q_2;\tau)}{P(-Q_1,-Q_2;\tau)}=\mbox{e}^{-4E[(\beta_1-\beta_0)Q_1+(\beta_2-\beta_0)Q_2]}.
\label{FRoneBis}
\end{equation}

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',\vec{Q'}} \langle s,\vec{Q}| {\cal M}_\textrm{ext}| s',\vec{Q'}\rangle P(s',\vec{Q'}; \tau)$ and whose explicit expression is derived from the balance equation
(\ref{BalanceP}). 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  (\ref{BalanceP}), 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  (\ref{defwn}). 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$
\begin{equation}
\frac{w_n^{(a)} \left(s, \Delta Q^{(1)}_n(s), \Delta Q^{(2)}_n(s)\right)}
{w_n^{(a)} \left(s_n, -\Delta Q^{(1)}_n(s), -\Delta Q^{(2)}_n(s)\right)}
=\mbox{e}^{-\beta_a \Delta Q^{(a)}_n(s)}
\end{equation}
This symmetry can be considered as an extension of the
 so-called  generalized detailed balance
\footnote{Several terminologies can be found in the literature :  ''local'' detailed balance \cite{LebowitzSpohn1999,Seifert2012},  "generalized"  detailed balance \cite{Derrida2007} or "modified" detailed balance \cite{Kurchan1998,CornuBauer2013I}.}
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 \cite{CornuBauer2013I,BauerCornu2015} 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 (\ref{FRoneBis}) 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_aQ_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_1Q_1+\beta_2Q_2)$. As in Ref.\cite{CornuBauer2013I,BauerCornu2015}  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
\begin{equation}
\Delta_{\rm exch}^{\rm excs, \beta_0} S=4E\left[ (\beta_1-\beta_0)Q_1 +(\beta_2-\beta_0)Q_2\right].
\end{equation}
Hence the fluctuation relation (\ref{FRoneBis}) 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
\begin{equation}
P(\Delta_{\rm exch}^{\rm excs, \beta_0}S ;\tau)=\mbox{e}^{-\Delta_{\rm exch}^{\rm excs, \beta_0} S \,}P(-\Delta_{\rm exch}^{\rm excs, \beta_0}S ;\tau).
\label{FRoneS}
\end{equation}

\section{Statistics of the time-integrated energy current}
\label{Qstat}

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

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
\begin{equation}
{\cal Q}=\tfrac{1}{2}\left(Q_1-Q_2\right),
\label{defQ}
\end{equation}
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 (\ref{defQ})
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 (\ref{resfinPi})  and the preceding definitions we then get
\begin{equation}
P({\cal Q};\tau)=
\int_{-\pi}^{\pi}\frac{\mbox{d} p}{2\pi}\,\mbox{e}^{-2{\rm i}p{\cal Q}}\,{\widehat{P}^\star}(p;\tau),
\label{defhPip}
\end{equation}
in which
\begin{equation}
{\widehat{P}^\star}(p;\tau)= \mbox{e}^{-\tfrac{1}{2}N\tau}
\prod_{q>0}\left[ \cosh \mu^\star_q(p) \tau
  + \frac{T^\star_q(p;\beta_0)}{S_q(\beta_0) \mu^\star_q(p)} \sinh \mu^\star_q(p) \tau \right].
\label{x1PQ}
\end{equation}
We observe that  $P({\cal Q};\tau)$  given by (\ref{defhPip}) and (\ref{x1PQ}) 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
\begin{equation}
\langle \mbox{e}^{\lambda {\cal Q}} \rangle =\widehat{P}^\star\left(- \frac{{\rm i} \lambda}{2};\tau\right).
\label{ExpMGF}
\end{equation}

\subsection{Cumulants of ${\cal Q}$ in the long-time limit}
\label{xcumulants}

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
 \begin{equation}
Ng_N(\lambda)\equiv\lim_{\tau\to \infty} \frac{1}{\tau} \ln
\langle \mbox{e}^{\lambda {\cal Q}}\rangle
\label{defgN}
\end{equation}
The cumulants of interest are the values of the derivatives of $g_N(\lambda)$ with respect to $\lambda$ taken at $\lambda=0$,
\begin{equation}
\frac{\langle {\cal Q}^n\rangle_{\rm c}}{N\tau} =
\frac{\mbox{d}^n g_N(\lambda)}  {\mbox{d} \lambda^n} \Bigg|_{\lambda=0}.
\label{cumulantsDef}
\end{equation}

According to (\ref{ExpMGF}) and (\ref{defgN}) together with the explicit value (\ref{x1PQ}) of $\widehat{P}^\star(-{\rm i} \lambda/2;\tau)$ we get
\begin{equation}
g_N(\lambda)=
-\frac{1}{2}+\frac{1}{N}\sum_{q>0}\mu^\star_q\left(- \frac{{\rm i}\lambda}{2}\right)
\label{expgN1}
\end{equation}
The expression for  $\mu_q^\star(p)=\mu_q(p,-p)$ is given by (\ref{x2muq}) where $\Theta(\vec{p})$ is defined in (\ref{defTheta}). We set $\theta (\lambda)=\Theta(-{\rm i} \lambda/2,{\rm i} \lambda/2)$ and get
\begin{equation}
g_N(\lambda) =
-\frac{1}{2}+\frac{1}{N}\sum_{q>0} \sqrt{(1-\gamma\cos q)^2 + \theta(\lambda) \sin^2 q}
\label{expgN2}
\end{equation}
with
\begin{equation}
\theta(\lambda)=2{\sf A}[\cosh\lambda -1] + 2{\sf B}\sinh\lambda,
\label{recallThetap}
\end{equation}
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 \cite{PopkovSchutz2011}. 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 (\ref{expgN1}) 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
(\ref{x1hPi}) 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 (\ref{Mrewrite}) and (\ref{defX}) lead to (\ref{expgN1}).

Eventually the  cumulants of the time-integrated energy current per site and  unit of time in the long-time limit  are given by
(\ref{cumulantsDef}) and (\ref{expgN2})
\begin{eqnarray}
\lim_{\tau\to\infty}\frac{\langle {\cal Q}\rangle}{N\tau} &=&
\tfrac{1}{2}{\sf B}\Sigma_1(N,\gamma),
\label{xfirstcum}
\\[2mm]
\lim_{\tau\to\infty}\frac{\langle {\cal Q}^2\rangle_{\rm c}}{N\tau} &=&
\tfrac{1}{2}\big[ {\sf A} \Sigma_1(N,\gamma)
- {\sf B}^2 \Sigma_2(N,\gamma) \big],
\label{xsecondcum}\\[2mm]
\lim_{\tau\to\infty}\frac{\langle {\cal Q}^3\rangle_{\rm c}}{N\tau} &=&
\tfrac{1}{2}\big[ {\sf B} \Sigma_1(N,\gamma)
- 3{\sf A}{\sf B} \Sigma_2(N,\gamma)
+ 3{\sf B}^3 \Sigma_3(N,\gamma) \big],
\label{xthirdcum}\\[2mm]
\lim_{\tau\to\infty}\frac{\langle {\cal Q}^4\rangle_{\rm c}}{N\tau} &=&
\tfrac{1}{2}\big[ {\sf A} \Sigma_1(N,\gamma)
- (3{\sf A}^2+4{\sf B}^2) \Sigma_2(N,\gamma)
+ 18{\sf A}{\sf B}^2 \Sigma_3(N,\gamma) \nonumber\\[2mm]
&& \phantom{\tfrac{1}{2}[}
-15{\sf B}^4 \Sigma_4(N,\gamma)\big].
\label{xfourthcum}
\end{eqnarray}
where we have introduced
\begin{equation}
\Sigma_n(N,\gamma) =
\frac{2}{N}\sum_{q>0}\,\,\frac{\sin^{2n} q}{(1-\gamma\cos q)^{2n-1}}\,.
\label{defSigman}
\end{equation}
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 (\ref{defbeta}).The $\Sigma_n(N,\gamma)$ are monotonically increasing with $\gamma$.

Expressions for higher order cumulants may be derived by increasing
algebraic effort.
Expression (\ref{xfirstcum}) for the time-averaged energy current has to
coincide with Eqs.\,(\ref{resjav})-(\ref{Javdirect})
of section \ref{secdirect}. Upon inserting the explicit expressions
for both one obtains the identity
\begin{equation}
\Sigma_1(N,\gamma)  =\frac{2}{N}\sum_{\ell=1}^{N/2}
\frac{\sin^2 \frac{(2\ell-1)\pi}{N}}
{1-\gamma\cos\frac{(2\ell-1)\pi}{N}}
\,=\,\frac{(1+\zeta^2)(1+\zeta^{N-2})}{2(1+\zeta^N)}\,, \qquad
N=2,4,6\ldots,
\label{ident1}
\end{equation}
where we recall that $\zeta=\tanh\beta E$ while
$\gamma=\tanh 2\beta E$.
Eq.\,(\ref{ident1}) 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$.
\vskip 0.3cm

The expressions for cumulants, of which the first four
are given in (\ref{xfirstcum})-(\ref{xfourthcum}), 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.
\footnote{Recall 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 \ref{secInfiniteSize}.

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
\cite{CornuBauer2013}.

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 (\ref{xfirstcum})-(\ref{xfourthcum}) are rewritten in
 dimensionful time $t=\tau/(\nu_1+\nu_2)$
 it appears that   the cumulants  \textit{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 \cite{CornuBauer2013}. (In other words, in dimensionful time, when $N=2$ the scaled generating function  for cumulants \textit{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.

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

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
\begin{equation}
{\cal Q}=\bar{\jmath} N\tau,
\label{defq}
\end{equation}
According to definition (\ref{defQ}) 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
\footnote{For a precise discussion of this definition see Ref.\cite{CornuBauer2013I} section 5 and Appendix E.}
\begin{equation}
\lim_{\varepsilon\to 0}\lim_{\tau\to \infty} -\frac{1}{N\tau}\ln P\left(\frac{{\cal Q}}{N\tau} \in [\bar{\jmath}-\varepsilon,\bar{\jmath}+\varepsilon]  \right)
= {\cal I}_N(\bar{\jmath}).
\label{defIN}
\end{equation}
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 (\ref{defIN}) by considering $P({\cal Q};\tau)$ and  applying the saddle point method   to its inverse Fourier transform representation (\ref{defhPip})-(\ref{x1PQ}) 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.\cite{CornuBauer2013}).

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  (\ref{defgN}). 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  (\ref{defgN}), exists
and  is differentiable for all real $\lambda$, as shown by its expression (\ref{expgN2})-(\ref{recallThetap}). Thus
$g_N(\lambda)$ satisfies the hypothesis of the simplified version of the
Gärtner-Ellis theorem (see, {\it e.g.,}
\cite{Touchette2009,DemboZeitouni1998}).
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,
\begin{equation}
 {\cal I}_N(\bar{\jmath}) =\max_{\lambda \in \mathbb{R}}\{ \lambda \bar{\jmath} - g_N(\lambda)  \}.
 \label{LFTransform}
\end{equation}
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,
\begin{equation}
{\cal I}_N(\bar{\jmath}) = \bar{\jmath} \lambda_{\bar{\jmath}} - g_N(\lambda_{\bar{\jmath}}),
\label{xINt}
\end{equation}
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 (\ref{expgN2}) 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.\cite{CornuBauer2013}.
\vspace{3mm}

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 (\ref{expgN2})-(\ref{recallThetap}) has the symmetry property
\begin{equation}
g_N(\lambda)=g_N\big( -\lambda-4(\beta_2-\beta_1)E \big).
\label{symmfNz}
\end{equation}
As a consequence ${\cal I}_N(\bar{\jmath})$ obeys the fluctuation relation
\begin{equation}
{\cal I}_N(\bar{\jmath})={\cal I}_N(-\bar{\jmath})+4(\beta_2-\beta_1)E\bar{\jmath}.
\label{INfluctrel}
\end{equation}
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 \cite{BodineauDerrida2007}.
We notice that the symmetry property (\ref{symmfNz}) determines the value of the large deviation function for $\bar{\jmath}=0$. Indeed expression (\ref{LFTransform}) for the large deviation function together with (\ref{symmfNz}) 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 (\ref{INfluctrel}) for ${\cal Q}$ can be derived from the finite time fluctuation relation (\ref{FRoneBis}) for the couple of variables $Q_1$ and $Q_2$.

\subsection{Infinite size chain at finite effective temperature}
\label{InfiniteChainFiniteT}

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 (\ref{expgN2})-(\ref{recallThetap}) reads
 \begin{eqnarray}
&&\lim_{N\to\infty} g_N(\lambda; \gamma)=-\frac{1}{2}
\\  \nonumber
&& +  \frac{1}{2\pi} \int_{0}^{\pi} \mbox{d} q\,\sqrt{(1-\gamma\cos q)^2+ 2\left[{\sf A}(\cosh\lambda-1) + {\sf B}\sinh\lambda\right] \sin^2q}.
 \label{limgNTer}
 \end{eqnarray}
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 (\ref{xfirstcum}) -(\ref{xfourthcum}) 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$
 \begin{equation}
\lim_{N\to\infty} \Sigma_n(N,\gamma) = \frac{1}{\pi}
\int_{0}^{\pi} \mbox{d} q\,\,\frac{\sin^{2n}q}{(1-\gamma\cos q)^{2n-1}}
\label{limSigma}
\end{equation}
is  finite for all $n\geq 1$.

\section{Various physical effects}
\label{Physicaleffects}

\subsection{Kinetic effects}
\label{seckinetic}

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 (\ref{expgN2}) and (\ref{recallThetap})
in a power series in $\nu_2/\nu_1$ we find that
\begin{equation}
(\nu_1+\nu_2)g_N(\lambda) = \frac{\nu_2}{2}\left\{
\left[p_+ \mbox{e}^{\lambda} + p_- \mbox{e}^{-\lambda} -(p_++p_-)\right]
\Sigma_1(N,\gamma_1)
+ {\cal O}\left(\frac{\nu_2}{\nu_1}\right)\right\},
\label{frw}
\end{equation}
 in which
\begin{equation}
p_+=\tfrac{1}{2}(1-\gamma_1)(1+\gamma_2),
\qquad
p_-=\tfrac{1}{2}(1+\gamma_1)(1-\gamma_2).
\label{xppm}
\end{equation}
The argument $\gamma_1$ of the function $\Sigma_{N,1}$ in (\ref{frw})
is the leading order
term of the expansion of $\gamma$ for small $\nu_2/\nu_1$.

The leading order term in (\ref{frw})
is in fact the scaled generating function
for the cumulants of a {\it 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$ \cite{CoxMiller1965,vanKampen1992}.
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 $\textrm{f}$ and $\textrm{s}$ denote the fast and slow thermostats, respectively, then the scaled generating function given in (\ref{frw}) takes the generic form
\begin{equation}
(\nu_1+\nu_2)g_N(\lambda) = \frac{\nu_\textrm{s}}{2}\left\{
\left[p_+ \mbox{e}^{\lambda} + p_- \mbox{e}^{-\lambda} -(p_++p_-)\right]
\Sigma_1(N,\gamma_\textrm{f})
+ {\cal O}\left(\frac{\nu_\textrm{s}}{\nu_\textrm{f}}\right)\right\},
\label{frwBis}
\end{equation}
From the generic relation (\ref{cumulantsDef}) cumulants read to leading order in $\nu_\textrm{s}/\nu_\textrm{f}$ \begin{eqnarray}
\nonumber
\lim_{t\to \infty}\frac{\langle {\cal Q}^{2m-1}\rangle_{\rm c}}{Nt} &=&
\frac{\nu_\textrm{s} }{2}\left[(\gamma_2-\gamma_1)\Sigma_1(N,\gamma_\textrm{f})
+ {\cal O}\left(\frac{\nu_\textrm{s}}{\nu_\textrm{f}}\right)\right],
\\[2mm]
\label{CumRW}
\lim_{t\to \infty}\frac{\langle {\cal Q}^{2m}\rangle_{\rm c}}{Nt} &=& \frac{\nu_\textrm{s}}{2} \left[(1-\gamma_1\gamma_2)\Sigma_1(N,\gamma_\textrm{f})
+ {\cal O}\left(\frac{\nu_\textrm{s}}{\nu_\textrm{f}}\right)\right],
\end{eqnarray}
for $m=1,2,\ldots$.
The latter expressions, with $t$ in place of $\tau$, can be retrieved from our expressions  (\ref{xfirstcum})-(\ref{xfourthcum}) by multiplying them by  $\tau/t=\nu_1 +\nu_2$ and expanding them to leading order in $\nu_\textrm{s}/\nu_\textrm{f}$.

\subsection{One thermostat  at zero temperature}
\label{secZeroTtwo}

Dissipation towards a thermostat at zero temperature was studied
by Farago and Pitard \cite{FaragoPitard2007,FaragoPitard2008}
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 (\ref{xABintro})
implies that ${\sf A} = {\sf B} = 1-\gamma_1$. Combined with
(\ref{recallThetap}) 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 (\ref{expgN2}), we get that when $\gamma_2=1$
\begin{equation}
g_N(\lambda)=
-\frac{1}{2}+\frac{1}{N}\sum_{q>0} \sqrt{(1-\gamma\cos q)^2 + 2\bar{\nu}_1\bar{\nu}_2(1-\gamma_1)\left[\mbox{e}^{\lambda}-1\right] \sin^2 q}.
\label{expgN2ZeroTtwo}
\end{equation}
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
 \begin{equation}
 {\cal I}_N(\bar{\jmath})=\infty, \qquad \bar{\jmath}<0.
 \end{equation}
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 \ref{xcumulants} 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 (\ref{xfirstcum}) and
(\ref{xsecondcum}) now read
\begin{eqnarray}
\lim_{\gamma_2\to 1}\lim_{t\to \infty}
\frac{\langle Q\rangle}{{N t}}&=& \frac{\nu_1+\nu_2}{2}\bar{\nu}_1\bar{\nu}_2(1-\gamma_1)\Sigma_{N,1}(\gamma),
\label{dissipationcum}
\\[2mm]
\nonumber
\lim_{\gamma_2\to 1}\lim_{t\to \infty}
\frac{\langle {\cal Q}^2 \rangle_{\rm c}}{{N t}} &=&
 \frac{\nu_1+\nu_2}{2}\bar{\nu}_1\bar{\nu}_2(1-\gamma_1)
\left[\Sigma_{N,1}(\gamma) -\bar{\nu}_1\bar{\nu}_2(1-\gamma_1)\Sigma_2(N,\gamma) \right],
\end{eqnarray}
where the inverse temperature $\beta$ in the special case $\gamma_2=1$ is  determined from
$\tanh2\beta E =\gamma =1-\bar{\nu}_1(1-\gamma_1)$.

\subsection{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.\,(\ref{frw})  becomes
\begin{equation}
(\nu_1+\nu_2)g_N(\lambda) = \frac{\nu_\textrm{s}}{2}\left\{
(1-\gamma_1)[ \mbox{e}^{\lambda} -1]\
\Sigma_1(N,\gamma_\textrm{f})
+ {\cal O}\left(\frac{\nu_\textrm{s}}{\nu_\textrm{f}}\right)\right\}.
\label{frw2}
\end{equation}
 This is the generating function for a {\it 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}$,
\begin{equation}
\lim_{t\to \infty}\frac{\langle Q^{n}\rangle_{\rm c}}{Nt}= \frac{1}{2}\nu_\textrm{s}
\left[(1-\gamma_1)\Sigma_1(N,\gamma_\textrm{f})
+ {\cal O}\left(\frac{\nu_\textrm{s}}{\nu_\textrm{f}}\right)\right]
\label{nu2ll1Tzero}
\end{equation}
for $n=1,2,\ldots$.
By comparing (\ref{dissipationcum}) and (\ref{nu2ll1Tzero}) one sees that the limits $T_2\to 0$ and $\nu_\textrm{s}\ll \nu_\textrm{f}$ commute.

\section{Large size and low effective temperature}
\label{secInfiniteSize}

\subsection{Parameters at low effective temperature}

We now consider the regime where $N\gg1$ and $0<1-\gamma\ll1$. According to the relation $\gamma=\gamma_2-\nu_1(\gamma_2-\gamma_1)$ the condition  $0<1-\gamma\ll1$  corresponds to
\begin{equation}
0\leq 1-\gamma_2\ll1
\end{equation}
 while
\begin{equation}
 0<\gamma_2-\gamma_1\ll 1\quad\textrm{and/or}\quad 0<\bar{\nu}_1 \ll 1.
\end{equation}
We notice that  in the case $\gamma_1=\gamma_2$  and $0\leq 1-\gamma_2\ll1$ 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 (\ref{xABintro}), as
\begin{equation}
{\sf A}=(1-\gamma)\, \bar{\nu}_1\bar{\nu}_2 \,{\sf a},
\qquad
{\sf B}=(1-\gamma)\, \bar{\nu}_1\bar{\nu}_2\, {\sf b}\, ,
\label{defab}
\end{equation}
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$.
\vskip0.3cm

For the sake of conciseness, from now on we denote the long time cumulants per site and unit of time
\begin{equation}
\kappa^{(n)}(N,\gamma)=\frac{1}{N}\lim_{t\to\infty}\frac{\langle {\cal Q}^n\rangle_{\rm c}}{(\nu_1+\nu_2) t}.
\label{defkappatilde}
\end{equation}
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 (\ref{CumRW}) in the case of two thermostats whose kinetic parameters are of  different orders of magnitude. The cumulant of order $n$ can be written as
\begin{equation}
\kappa^{(n)}_\textrm{RW}(N,\gamma)=(1-\gamma)\,\,\Sigma_1(N,\gamma)\,\,
 {\sf k}_n,
\label{DefKappaRW}
\end{equation}
with the definition
\begin{equation}
{\sf k}_n= \frac{1}{2} \bar{\nu}_1\bar{\nu}_2 \left[\frac{1+(-1)^n}{2}{{\sf a}} + \frac{1-(-1)^n}{2}{{\sf b}}\right],
\label{defc1}
\end{equation}
where ${\sf a}$ and ${\sf b}$ are defined in (\ref{defab}).
As illustrated by the expressions (\ref{xfirstcum})-(\ref{xfourthcum}) 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
 \begin{equation}
\kappa^{(1)}(N,\gamma)=\kappa^{(1)}_\textrm{RW}(N,\gamma),
\label{rescaleKappa1}
\end{equation}
while for $n\geq 2$
\begin{equation}
\kappa^{(n)}(N,\gamma)=\kappa^{(n)}_\textrm{RW}(N,\gamma) + \Delta \kappa^{(n)}(N,\gamma)
\label{decomp}
\end{equation}
where the deviation $ \Delta \kappa^{(n)}(N,\gamma)$ from the random walk process reads
\begin{equation}
\Delta \kappa^{(n)}(N,\gamma)=  \sum_{p=2}^{n}(1-\gamma)^p \,\Sigma_p(N,\gamma)  (\bar{\nu}_1\bar{\nu}_2 )^p{\sf c}_p^{(n)}({\sf a}, {\sf b}),
\label{rescaleKappan}
\end{equation}
In (\ref{rescaleKappan}) 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 (\ref{cumulantsDef}) for  every  cumulant  per site and unit of time and the expression (\ref{expgN2})-(\ref{recallThetap}) for their  generating function.

\subsection{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 (\ref{expgN2})-(\ref{recallThetap}) 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\to1} {\cal I}_N(\bar{\jmath}; \gamma)$.

The random walk contribution to the cumulant $\kappa^{(n)}$ is defined in (\ref{DefKappaRW}).
According to the explicit expression (\ref{ident1}) 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
\begin{equation}
\lim_{\gamma\to 1} \frac{\kappa^{(n)}_\textrm{RW}(N,\gamma)}{1-\gamma}= {\sf k}_n,
\label{FirstOrderBehaviorRW}
\end{equation}
where ${\sf k}_n$ is defined in (\ref{defc1}).

We now turn to  the ${\cal Q}$ process.  By virtue of (\ref{rescaleKappa1}) 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
(\ref{FirstOrderBehaviorRW}).
Besides, for all $n\geq 2$ the coefficient $\Sigma_n(N,\gamma)$ defined in (\ref{defSigman}), remains finite when $\gamma=1$ at $N$ fixed. Thus, according to the expression (\ref{rescaleKappan}), 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$,
\begin{equation}
\Delta \kappa^{(n)}(N,\gamma)\underset{\gamma\to 1}={\cal O}\left((1-\gamma)^2\right).
\end{equation}
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 (\ref{FirstOrderBehaviorRW}) it reads
\begin{equation}
\lim_{\gamma\to 1} \frac{\kappa^{(n)}(N,\gamma)}{1-\gamma}= {\sf k}_n.
\label{FirstOrderBehavior}
\end{equation}

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 (\ref{rescaleKappa1}), 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 (\ref{DefKappaRW}) for the random walk cumulant, it reads
 \begin{equation}
 \Delta \kappa^{(1)}= (1-\gamma) \left[\Sigma_1(N,\gamma)-1\right]{\sf k}_1.
 \label{defDeltaKappa1}
\end{equation}
It can be rewritten in terms of a single finite sum
\begin{equation}
 \Delta \kappa^{(1)}= (1-\gamma)^2 \,{\sf k}_1\sum_{\ell =1}^{N/2}\Delta s_{1,\ell}(N,\gamma)
\end{equation}
where the increment $\Delta s_{1,\ell}(N,\gamma)$ is written in  (\ref{defDeltas1l}). This finite sum indeed converges when $\gamma\to 1$, and
 \begin{equation}
 \Delta \kappa^{(1)}\underset{\gamma\to 1}={\cal O}\left((1-\gamma)^2\right).
\end{equation}

\subsection{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 (\ref{InfiniteChainFiniteT}),  the large deviation function exists and  all cumulants are finite.
When $\gamma\to1$,  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\vert_{\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 (\ref{DefKappaRW})
\begin{equation}
\lim_{N\to\infty,\gamma\to 1} \frac{\kappa^{(n)}_\textrm{RW}(N,\gamma)}{1-\gamma}
= {\sf k}_n.
\label{kappaRWinfinitechain}
\end{equation}
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 (\ref{ident1}), 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 (\ref{rescaleKappa1}). 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 (\ref{rescaleKappa1})-(\ref{decomp}), 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  (\ref{decomp})-(\ref{rescaleKappan}).
For $N$ finite  $\Delta \kappa^{(n)}(N,\gamma)$ vanishes as $(1-\gamma)^2$  when $\gamma\to1$. 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 (\ref{defSigman}) and as can also be seen in the integral representation (\ref{limSigma}) 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
\begin{equation}
\lim_{N\to\infty,\gamma\to 1} \frac{\kappa^{(n)}(N,\gamma)}{1-\gamma}
= {\sf k}_n.
\label{kappainfinitechain}
\end{equation}

 Now we turn to the correction to this leading $1-\gamma$ term. In the regime where $N\gg1$ and $0<1-\gamma\ll1$, for every cumulant $\kappa^{(n)}(N,\gamma)$ the correction to  the leading $(1-\gamma)$-term (\ref{FirstOrderBehavior})
is the sum of two contributions arising from  (\ref{rescaleKappa1})-(\ref{decomp}):  on the one hand, the  correction (\ref{defDeltaKappa1}) to the leading $(1-\gamma)$-term in the first moment of the random walk defined in (\ref{DefKappaRW}),
and on the other hand,  the leading behavior  of  the deviation $\Delta \kappa^{(n)}$ defined in (\ref{rescaleKappan}).
In the double limit $N\gg1$ and $0<1-\gamma\ll1$,  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 \ref{Details}, 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
\begin{eqnarray}
\nonumber
&& \textrm{scaling regime $[\textrm{I}]$:} \quad N\sqrt{1-\gamma)} \to  +\infty
\\
\nonumber
&&\textrm{scaling regime $[\textrm{II}]$:} \quad N\sqrt{1-\gamma)} = \rho /\sqrt{2} \quad \textrm{with} \quad
0\leq \rho<\infty.
\label{defregimescaling}
\end{eqnarray}
Eventually in  scaling regime $[\textrm{I}]$, the cumulants behave as
\begin{equation}
\kappa^{(n)}(N,\gamma)
\underset{\textrm{scl} \, [\textrm{I}]}{\sim}
\left[ 1-\gamma +  (1-\gamma)^{3/2}\,
\bar{\nu}_1\bar{\nu}_2\,F_n^{[\textrm{I}]}(\bar{\nu}_1\bar{\nu}_2; {\sf a} ,{\sf b}) \right]
{\sf k}_n
\label{corr1rescaleKappan1}
\end{equation}
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}]$
\begin{equation}
\kappa^{(n)}(N,\gamma)
\underset{\textrm{scl} \, [\textrm{II}]}{\sim}
\left[ 1-\gamma  +  (1-\gamma)^{3/2}\,
\bar{\nu}_1\bar{\nu}_2\,  \rho\,F_n^{[\textrm{II}]}(\rho ,\bar{\nu}_1\bar{\nu}_2; {\sf a} ,{\sf b})\right]{\sf k}_n .
\label{corr1rescaleKappan2}
\end{equation}
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\to0$.
We notice that the factor $\bar{\nu}_1\bar{\nu}_2$ in (\ref{corr1rescaleKappan1}) and
(\ref{corr1rescaleKappan2})  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 (\ref{DefKappaRW}).

\subsection{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 \cite{Glauber1963} 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
\begin{equation}
\lim_{t\to\infty} \frac{\langle{\cal Q}^n\rangle_{\rm c}}{t/t_\textrm{rel}}=\lim_{t\to\infty} \frac{\langle {\cal Q}^n\rangle_{\rm c}}{(1-\gamma) t},
\end{equation}
which are referred to as  ``long time rescaled cumulants''  in the following.
According to (\ref{kappainfinitechain})  all rescaled cumulants for the whole chain scale as the system size $N$ in the low-temperature regime,
\begin{equation}
\lim_{t\to\infty}\frac{\langle {\cal Q}^n\rangle_{\rm c}}{t / t_\textrm{rel}} \underset{\xi\gg 1 \, , \,N\gg 1}{\quad \sim \quad} N \,{\sf k}_n
\label{limKdominant}
\end{equation}
where the finite coefficient  $k_n$ is  a random walk cumulant given by (\ref{defc1}):  for $m\geq 1$
\begin{eqnarray}
{\sf k}_{2m-1}&=&\frac{1}{2} \bar{\nu}_1\bar{\nu}_2 \, {\sf b}
\\
\nonumber
{\sf k}_{2m}&=&\frac{1}{2} \bar{\nu}_1\bar{\nu}_2 \, {\sf a}
\end{eqnarray}
where ${\sf a}$ and ${\sf b}$ are defined in (\ref{defab}).

Second, the correlation length $\xi$ is defined from  the  correlation $\langle s_ns_{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=\tanh2\beta E$.
The equilibrium correlation $\langle s_ns_{n+r}\rangle$ in the  Ising chain with finite size $N$ reads $\langle s_ns_{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_ns_{n+r}\rangle = \zeta^r$ at any distance $r$. Therefore the dimensionless  correlation length $\xi(\beta)$ in the system is
\begin{equation}
\xi (\beta)= [-\ln \tanh \beta E]^{-1}
\label{defxicorr}
\end{equation}
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
\begin{equation}
\sqrt{2(1-\gamma)}=\frac{1}{\xi} +{\cal O}\left(\frac{1}{\xi^2}\right).
\label{gammaxi}
\end{equation}
\vskip0.3cm

In the low effective temperature regime, at leading order  all rescaled cumulants for the whole chain scale as the system size $N$ (see (\ref{limKdominant})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)^2N\gg1$. According to (\ref{gammaxi}) 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 (\ref{corr1rescaleKappan1}) for the whole chain can be rewritten for $n\geq 1$ as
\begin{equation}
  \lim_{t\to\infty}\frac{\langle {\cal Q}^n\rangle_{\rm c}}{t / t_\textrm{rel}}
 \underset{1 \ll \xi \ll N}{\sim}  \left[ N + \frac{N}{\xi} f_n^{[\textrm{I}]}\right] {\sf k}_n
\end{equation}
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}$,
\begin{equation}
\frac{N}{\xi}  \underset{\textrm{scl} \, [\textrm{I}]}{\sim}N_\textrm{dw}
\quad \textrm{with}  \quad 1\ll N_\textrm{dw} \ll N
\end{equation}
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 (\ref{gammaxi}),
\begin{equation}
\frac{N}{\xi}  \underset{\textrm{scl} \, [\textrm{II}]}{\sim} \rho < \infty.
\end{equation}
Then the behavior (\ref{corr1rescaleKappan2}) of the rescaled cumulants for the whole chain can be rewritten for $n\geq 1$ as
\begin{equation}
 \lim_{t\to\infty}\frac{\langle {\cal Q}^n\rangle_{\rm c}}{t / t_\textrm{rel}}
\underset{\textrm{scl} \, [\textrm{II}]}{\sim}  \left[ N + \rho^2 f_n^{[\textrm{II}]}(\rho)\right] {\sf k}_n
\end{equation}
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
\cite{DerridaDoucotRoche2004},  on a ring with periodic boundary conditions
 \cite{AppertDerridaETAL2008},
 for a one-dimensional hard particle gas on a ring or with open boundary conditions
\cite{BrunetETAL2010}
and for the one-dimensional lattice gas model  $ABC$
in the vicinity of a phase transition \cite{GerschenfeldDerrida2011}.
In Ref.\cite{MallickProlhac2009}  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.

\section{Conclusion}
\label{secconclusion}

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 \cite{Glauber1963} 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.\cite{FaragoPitard2007,FaragoPitard2008,PopkovSchutz2011}, 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.\cite{CornuBauer2013}. 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 $\textrm{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.