# Graph-combinatorial approach for large deviations of Markov chains

Giorgio Carugno<EMAIL_ADDRESS>Department of Mathematics, King’s
College London, Strand, London WC2R 2LS, UK  Pierpaolo Vivo
<EMAIL_ADDRESS>Department of Mathematics, King’s College London,
Strand, London WC2R 2LS, UK  Francesco Coghi<EMAIL_ADDRESS>Nordita,
KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg
12, SE-106 91 Stockholm, Sweden

###### Abstract

We consider discrete-time Markov chains and study large deviations of the pair
empirical occupation measure, which is useful to compute fluctuations of pure-
additive and jump-type observables. We provide an exact expression for the
finite-time moment generating function, which is split in cycles and paths
contributions, and scaled cumulant generating function of the pair empirical
occupation measure via a graph-combinatorial approach. The expression obtained
allows us to give a physical interpretation of interaction and entropic terms,
and of the Lagrange multipliers, and may serve as a starting point for sub-
leading asymptotics. We illustrate the use of the method for a simple two-
state Markov chain.

##  I Fluctuations for discrete-time Markov chains in the large deviations
regime

Markov chains are widely-used stochastic models of in and out-of-equilibrium
physical systems. We consider a discrete-time ergodic Markov chain
$X=\left(X_{\ell}\right)_{\ell=1}^{n+1}=(X_{1},X_{2},\dots,X_{n+1})$ evolving
in a finite discrete state space $\Gamma$ of $N$ states according to the
(irreducible and aperiodic) transition matrix $\Pi$. The matrix $\Pi$
characterises the probability of going from a state $X_{\ell}=i$ at time
$\ell$ to a state $X_{\ell+1}=j$ at time $\ell+1$. We will use the index
$\ell$ to refer to time and the indices $i$ and $j$ to refer to general states
of the state space.

In this setting, one and two-point observables having the general form

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}f(X_{\ell},X_{\ell+1})\,,$ (1)

where $f$ is any function that may depend both on the starting and landing
state, are of fundamental importance to characterise the typical and
fluctuating behaviour of the associated physical systems. [Notice that by
taking $f(i,j)=g(i)$, $C_{n}$ in (1) can also cover the case of purely time-
additive observables111This marginalisation only works for discrete-time
processes..] Just to give an example, the observable in (1) can represent the
number of transitions over $[1,n]$ in a particular subset of the state space
[1], obtained by fixing $f=\mathbf{1}_{\Delta}$, with $\Delta$ the
characteristic function of the subset. Furthermore, in certain contexts,
$C_{n}$ can also express heat [2], two-point correlation functions, activities
[3, 4, 5, 6], particle and energy currents [7], efficiency [8, 9, 10, 11, 12],
entropy production [13, 14, 6], and many others.

To study fluctuations of $C_{n}$, the probabilistic theory of large deviations
and, in particular, the Donsker–Varadhan approach may be used as it offers
analytical and numerical methods to calculate the large deviation (or rate)
function

$I(c)=-\lim_{n\rightarrow\infty}\frac{1}{n}\ln\mathbb{P}(C_{n}=c)$ (2)

characterising the time-leading exponential behaviour – provided there is one
– of the probability distribution $\mathbb{P}(C_{n}=c)$ [15, 16, 17, 18, 19,
1, 20]. The rate function in (2) is always positive and measures the extent of
the fluctuations of $C_{n}$ around its typical value $c^{*}$, which, for
ergodic Markov chains, is the unique zero of $I$ [16, 17, 18]. The existence
of the rate function $I$ is referred to as the validity of a large deviation
principle for the observable $C_{n}$ and can be seen as an extension of the
weak law of large numbers as it provides information on the speed –
exponential in $n$ – of convergence of $C_{n}$ to $c^{*}$.

In the context of Markov chains, there are several ways to compute the rate
function $I$. It is known that, by means of spectral large deviation
techniques (see, for instance, [21, 22, 23, 24, 5, 25]), one could calculate
the scaled cumulant generating function (SCGF)

$\Psi_{N}(s)\coloneqq\lim_{n\rightarrow\infty}\frac{1}{n}\ln\mathbb{E}\left[e^{nsC_{n}}\right]\,,$
(3)

where $s$ is the Lagrange (or tilting, in the large deviation jargon)
parameter dual to $C_{n}=c$. The SCGF $\Psi$ represents the leading
exponential behaviour of the moment generating function, associated with the
observable in (1). To obtain the rate function $I$ in (2), it would then be
enough to Legendre–Fenchel transform the SCGF, provided it be a differentiable
function – a result known as Gärtner–Ellis theorem [16, 17, 18]. Although
these methods serve well to the scope, variational techniques can also be
employed and one may derive the rate function $I$ by solving a variational
problem [3, 26, 20]. The advantage of employing variational methods is, at
least, twofold. In case of non-analytically solvable problems, variational
methods offer ways to bound the true rate function (see, for instance, [27,
28]) and, at the same time, alternative numerical techniques – inherited from
the fields of optimization theory and PDEs – are available to compute it [29].

In the case considered here, it is known that all the information on the
fluctuations of one and two-point observables can be obtained by studying the
pair empirical occupation measure

$L^{(2)}_{n}(i,j)=\frac{1}{n}\sum_{\ell=1}^{n}\delta_{X_{\ell},i}\delta_{X_{\ell+1},j}\hskip
28.45274pt\forall i,j\in\Gamma\,,$ (4)

as the value of $C_{n}$ can be deduced via the formula

$C_{n}=\sum_{i,j=1}^{N}f(i,j)L_{n}^{(2)}(i,j)\,.$ (5)

Interestingly, the long-time behaviour of (4), denoted by
$\rho=(\rho_{ij})_{i,j=1}^{N}$, can be interpreted as the amount of time that
the Markov chain $X$ spends transiting from a state $i$ to a state $j$ of
$\Gamma$ [16, 17, 18].

The pair empirical occupation measure of (4) is known to satisfy a large
deviation principle of the form

$\mathbb{P}\left(L_{n}^{(2)}=\nu\right)=e^{-nH[\nu]+o(n)}\,,$ (6)

with rate function

$H[\nu]=\sum_{i,j}\nu_{ij}\ln\left(\frac{\nu_{ij}}{\mu_{i}\Pi_{ij}}\right)\,,$
(7)

where $\nu=(\nu_{ij})_{i,j=1}^{N}$ belongs to the set of probability measures
satisfying two constraints: the global balance on the state space, i.e.,
$\sum_{j}\nu_{ij}=\sum_{j}\nu_{ji}$, such that the sum of probability density
currents flowing in and out of an arbitrary state $i$ is conserved, and the
normalisation $\sum_{i,j}\nu_{ij}=1$ (with $\sum_{j}\nu_{ij}=\mu_{i}$ the
occupation measure). The rate function $H$ in (7) is known to be finite,
continuous, and convex for densities $\nu$ that satisfy the global balance on
the state space, featuring minimum and zero for $\nu=\rho$ [16]. Here, it is
interesting to notice that the rate function $I$ associated with $C_{n}$, can
be obtained variationally by solving the following contraction222This term is
commonly used in the large deviation theory jargon when referring to the
solution of a variational problem passing from a higher-up level in the large
deviation hierarchy to a lower one. (minimisation) problem

$I(c)=\inf_{\begin{subarray}{c}\nu:\\\
c=\sum_{i,j}f(i,j)\nu_{ij}\end{subarray}}H[\nu]\,,$ (8)

where the constraint appearing beneath the $\inf$ symbol is the formula (5),
which selects $c$, the fluctuation of interest for the observable $C_{n}$ in
(1).

The functional $H$ in (7) is thus a key ingredient for the variational study
of fluctuations in discrete-time Markov chains and, as mentioned, it plays a
pivotal role in statistical mechanics as many interesting dynamical
observables arising in physics have the two-point form in (1).

The form in (7) has been derived with various methods. Among these, the
exponential tilting procedure combined with the Radon–Nikodym change of
measure [16] (see [26] for continuous-time processes) holds a leading position
as it offers a simple and straightforward way to tackle the calculation,
provided that the form of the rate function for the i.i.d. process (or any
other useful process) is known. We will review and discuss this method in
Section II.

Although simple and well suited to large deviation estimates, the exponential
tilting procedure does not allow for the calculation of $o(n)$-exponential
sub-leading terms in the probability distribution of the pair empirical
occupation measure (4). In the probability and applied statistics literature,
however, exact combinatorial derivations that work at finite time can be
found. These may lead to the evaluation of sub-leading order terms that,
although not significant in the large deviation regime, would be important if
one wanted to study transition regimes. The first combinatorial result goes
back to [30], later on reviewed in [31], and more recently recalled in [32].
Another graph-combinatorial derivation for the probability distribution of the
pair empirical occupation measure was proposed in [33] and later on extended
in [34]. More recently, [35] provided an explicit—although not fully
rigorous—expression of subleading terms in (6), and constructed a Gauge theory
for typical fluctuations of $C_{n}$ around its expected value.

In the main Section III of our paper we use similar arguments to provide an
alternative, exact expression for the moment generating function of the pair
empirical occupation measure. We make use of notation and terminology that are
more familiar to the theoretical physics audience, and show—in line with
previous literature [17]—that our expression for the SCGF, akin to a Helmoltz
(canonical) free energy, allows us to give a straightforward physical
interpretation of all the terms and of the Lagrange multipliers that fix the
necessary constraints. Furthermore, we establish a direct link with spectral
methods and show an alternative variational formulation of the so-called
driven process [36, 37, 1, 20] (the Markov process responsible for the
creation of fluctuations in the large-deviation regime). In section IV we show
explicitly in a general two-state model the equivalence of our approach and
the standard spectral techniques to compute the moment generating function at
finite time $n$.

##  II Pair empirical measure rate functional

In this Section we show how the rate functional $H$ in (7) can be derived via
the exponential-tilting method. We start by writing the path-probability
definition

$\displaystyle\mathbb{P}\left(L_{n}^{(2)}=\frac{T}{n}\right)\coloneqq\mathbb{P}\left(L_{n}^{(2)}(i,j)=\frac{t_{ij}}{n}\;\forall
i,j\in\Gamma\right)$
$\displaystyle=\sum_{X_{1},X_{2},\ldots,X_{n+1}}\mathbb{P}(X_{1},X_{2},\cdots,X_{n+1})\delta_{L_{n}^{(2)},T/n}$
(9)
$\displaystyle=\sum_{X_{1},X_{2},\dots,X_{n+1}}\mathbb{P}(X_{1})\Pi_{X_{1},X_{2}}\cdots\Pi_{X_{n},X_{n+1}}\delta_{L_{n}^{(2)},T/n}\,,$
(10)

where $t_{ij}$ represents the number of jumps that the Markov chain $X$ makes
between nodes $i$ and $j$, and in (10) we make use of the Markov property. We
also notice in (9) that we can interpret the set of $t_{ij}$s as the elements
of a matrix $T$, which will be a central object in the rest of this work.

We now introduce a new i.i.d. process
$X^{\prime}=\left(X^{\prime}_{\ell}\right)_{\ell=1}^{n+1}=(X^{\prime}_{1},X^{\prime}_{2},\dots,X^{\prime}_{n+1})$
based on the probability distribution $\zeta=\left(\zeta_{i}\right)_{i=1}^{N}$
on the state space and with its own pair empirical occupation measure that,
with abuse of notation, have the same form of (4). A large deviation principle
for the pair empirical measure of $X^{\prime}$ is known to hold (see, for
instance, Chapter 9 of [38] or Section II.2 of [16]) with rate functional

$H_{\text{i.i.d.}}[\nu]=\sum_{i,j}\nu_{ij}\ln\left(\frac{\nu_{ij}}{\mu_{i}\zeta_{j}}\right)\,.$
(11)

Consequently, we multiply and divide in the summation of (10) by the path-
probability $\mathbb{P}^{\prime}(L_{n}^{(2)}=T/n)$ of this i.i.d. process and
then introduce an exponential function as follows

$\displaystyle\mathbb{P}\left(L_{n}^{(2)}=\frac{T}{n}\right)$
$\displaystyle=\sum_{X_{1},X_{2},\dots,X_{n+1}}\mathbb{P}(X_{1})\Pi_{X_{1},X_{2}}\dots\Pi_{X_{n},X_{n+1}}\frac{\mathbb{P^{\prime}}(X_{1})\mathbb{P^{\prime}}(X_{2})\dots\mathbb{P^{\prime}}(X_{n+1})}{\mathbb{P^{\prime}}(X_{1})\mathbb{P^{\prime}}(X_{2})\dots\mathbb{P^{\prime}}(X_{n+1})}\delta_{L_{n}^{(2)},T/n}$
(12)
$\displaystyle=\sum_{X_{1},X_{2},\dots,X_{n+1}}\frac{\mathbb{P}(X_{1})}{\mathbb{P^{\prime}}(X_{1})}e^{\sum_{\ell=1}^{n}\left[\ln\Pi_{X_{\ell},X_{\ell+1}}-\ln\zeta_{X_{\ell+1}}\right]}\mathbb{P^{\prime}}(X_{1})\mathbb{P^{\prime}}(X_{2})\dots\mathbb{P^{\prime}}(X_{n+1})\delta_{L_{n}^{(2)},T/n}\,.$
(13)

The derivation continues by observing, in the exponential function, the
equality

$\sum_{\ell=1}^{n}\left[\ln\Pi_{X_{\ell},X_{\ell+1}}-\ln\zeta_{X_{\ell+1}}\right]=n\sum_{i,j=1}^{N}L_{n}^{(2)}(i,j)\left[\ln\Pi_{ij}-\ln\zeta_{j}\right]\,,$
(14)

obtained by using the definition of the pair empirical measure (7). Hence, we
get

$\mathbb{P}\left(L_{n}^{(2)}=\frac{T}{n}\right)=\sum_{X_{1},X_{2},\dots,X_{n+1}}\frac{\mathbb{P}(X_{1})}{\mathbb{P^{\prime}}(X_{1})}e^{n\sum_{i,j=1}^{N}L_{n}^{(2)}(i,j)\left[\ln\Pi_{ij}-\ln\zeta_{j}\right]}\mathbb{P^{\prime}}(X_{1})\mathbb{P^{\prime}}(X_{2})\dots\mathbb{P^{\prime}}(X_{n+1})\delta_{L_{n}^{(2)},T/n}\,.$
(15)

Eventually, by taking minus the logarithm of the probability $\mathbb{P}$,
dividing by $n$, and taking the limit $n\rightarrow\infty$ we get

$\displaystyle-\lim_{n\rightarrow\infty}\frac{1}{n}\ln\mathbb{P}\left(L_{n}^{(2)}=\nu\right)$
$\displaystyle=\sum_{i,j=1}^{N}\nu_{ij}\ln\frac{\zeta_{j}}{\Pi_{ij}}+\lim_{n\rightarrow\infty}\frac{1}{n}\ln\left(\sum_{X_{1},X_{2},\dots,X_{n+1}}\frac{\mathbb{P}(X_{1})}{\mathbb{P^{\prime}}(X_{1})}\mathbb{P^{\prime}}(X_{1})\mathbb{P^{\prime}}(X_{2})\dots\mathbb{P^{\prime}}(X_{n+1})\delta_{L_{n}^{(2)},\nu}\right)$
$\displaystyle=\sum_{i,j=1}^{N}\nu_{ij}\left(\ln\frac{\zeta_{j}}{\Pi_{ij}}-\ln\left(\frac{\nu_{ij}}{\mu_{i}\zeta_{j}}\right)\right)\,,$
(16)

where the matrix $\nu$ is defined as

$\nu=\frac{T}{n}\,.$ (17)

In the derivation of (16), we make use of the fact that
$L_{n}^{(2)}\rightarrow\nu$, and also that for the probability
$\mathbb{P}^{\prime}(L_{n}^{(2)}=\nu)$ a large deviation principle holds with
rate functional (11). The last formula obtained in (16) is exactly (7). We
remark that it is only because of the long-time limit that we can get rid of
the boundary term $\mathbb{P}(X_{1})/\mathbb{P}^{\prime}(X_{1})$ in (16) and
thus get the form of the rate functional for the pair empirical occupation
measure of the Markov process $X$. We also notice that, although extremely
useful, the use of an i.i.d. process with its pair empirical rate functional
is not strictly necessary for the purpose of the proof. Indeed, if the
asymptotics of the pair empirical probability of another process were known
and easy to handle, we could have tilted the path probability measure of the
Markov process in (12) with respect to it and we would have obtained the same
result. For further details on this and on how to best use the tilting method
we refer to [1].

The derivation presented in this Section makes use of methods that are well
known in the large deviation community. Nevertheless, for a more rigorous
proof of the large deviation principle for the pair empirical measure (4)
having rate functional (7) – which focuses on lower and upper bounds over
closed and open sets – we refer the reader to [38, 16, 18].

The derivation presented in this Section takes into consideration only leading
order terms in $n$ and, furthermore, lacks some physical interpretations of
the form of the rate functional $\eqref{eq:LevPair}$. The finite $n$
behaviour, captured by subleading terms in (6), is in general much harder to
study than the large deviations regime. For the continuous-time setting, in
[39] the authors use matrix product states to study finite-time large
fluctuations of one-dimensional lattice models. For discrete Markov Chains,
estimates and bounds for subleading terms in (6) are known in the literature
[40, 41], and derived by using spectral methods. In [35] the author proposes a
characterisation of subleading terms using graph-combinatorial arguments.
Using a similar approach as [35], we provide an exact formula for the moment
generating function valid for any finite $n$, a first step towards an
alternative derivation of the subleading terms in (6).

##  III Graph-combinatorial approach

In this Section, we present an alternative derivation of the rate functional
associated with the pair empirical occupation measure in (4). The proposed
derivation moves the focus from the probability distribution $\mathbb{P}$ and
rate functional $H$ in (4) to the moment generating function $Z_{N,n}$ and
SCGF

$\lambda_{N}[s]\coloneqq\lim_{n\rightarrow\infty}\frac{1}{n}\ln
Z_{N,n}[s]=\lim_{n\rightarrow\infty}\frac{1}{n}\ln\mathbb{E}\left[e^{ns\cdot
L_{n}^{(2)}}\right]\,$ (18)

where, with abuse of notation with respect to (3),
$s=\left(s_{ij}\right)_{i,j=1}^{N}$ is now a set of Lagrange parameters. This
paradigm shift is equivalent to a change of ensemble in statistical mechanics
[1]. Instead of working with the probability distribution
$\mathbb{P}(L_{n}^{(2)}=T/n)$ at a fixed $t$, we introduce Lagrange parameters
$s_{ij}$s that fix the $t_{ij}$s only on average, and thus work with a moment
generating function. The equilibrium statistical mechanics analogue would be a
change from the microcanonical ensemble, where the energy is fixed, to the
canonical ensemble, where only the average energy is fixed by the Lagrange
parameter $\beta$, the inverse temperature.

In this canonical framework, thanks to Markovianity and ergodicity, it is
known [17, 18] that we can map the large deviation problem to a spectral one.
This is because the SCGF can be calculated as the logarithm of the dominant
eigenvalue of the so-called tilted matrix
$\Pi_{s}=\left(\Pi_{s}\right)_{ij}\forall i,j\in\Gamma$, which has the form

$\left(\Pi_{s}\right)_{ij}=\Pi_{ij}e^{s_{ij}}.$ (19)

Noticeably, thanks to a graph-combinatorial mapping [33, 34], we can derive an
exact expression for the moment generating function $Z_{N,n}$ at finite $N$
and $n$. In principle, the exact form $Z_{N,n}$ allows one to evaluate sub-
leading terms (in $n$) that cannot be calculated within a purely large
deviation approach as that of Section II. Historically, graph-combinatorial
arguments similar to those used in this work have been proposed for cyclic
Markov chains by Dawson and Good in [33], and later on extended for general
Markovian paths by Goodman in [34]. The derivation that follows explains in
the details, with a theoretical-physics approach, a similar graph-
combinatorial calculation but moves the focus onto the moment generating
function $Z_{N,n}$ of the pair empirical occupation measure. This allows us to
naturally give a physical interpretation of the interaction and entropic terms
in the SCGF $\lambda_{N}$ (18).

###  III.1 An alternative expression for the moment generating function

The graph-combinatorial approach is based on the representation of the state
space connectivity as a graph $\mathbf{G}$ with associated adjacency matrix
$A$ – see Fig. 2(a). This has elements $A_{ij}=1$ if state $j$ can directly be
reached from $i$, and $0$ otherwise. In this context, we will refer to states
also as nodes or vertices. The transition matrix $\Pi$ of the Markov chain
$X$, in turn, embeds in its elements the connectivity of the state space as
$\Pi_{ij}=A_{ij}p_{ij}$ with $p_{ij}$ the jump probability between $i$ and
$j$.

The moment generating function of the probability
$\mathbb{P}(L_{n}^{(2)}=T/n)$ is

$Z_{N,n}(s)=\sum_{X_{1},\dots,X_{n+1}}\mathbb{P}(X_{1})\prod_{\ell=1}^{n}\Pi_{X_{\ell},X_{\ell+1}}e^{\sum_{ij}s_{ij}\delta_{X_{\ell},i}\delta_{X_{\ell+1},j}}\
,$ (20)

where $\mathbb{P}(X_{1})$ indicates the probability distribution of our
process at initial time $n=1$ and $s=(s_{ij})_{i,j=1}^{N}$ indicates the set
of tilting parameters.

The specific form of the distribution $\mathbb{P}(X_{1})$ will play a role for
finite time behaviour or sub-leading asymptotics, but it will not matter in
the large deviation regime – it only amounts to a boundary term – provided
that the graph $\mathbf{G}$ is connected. For convenience, we choose
$\mathbb{P}(X_{1})=\delta_{X_{1},1}$, viz. the starting node is fixed to be
node $1$.

The core idea of our work is to perform a change of variables: we transform
the sum over all states $X_{\ell}$, for $\ell\in\\{1,\dots,n+1\\}$, to a sum
over variables $t_{ij}\in\\{0,\dots,n\\}$, for nodes $i,j\in\mathbf{G}$. The
new variable $t_{ij}$, as in the previous Section II, is the number of times
the Markov chain jumps from state $i$ to state $j$, in particular $t_{ij}$ can
be different from zero only if there is an edge in $\mathbf{G}$ between nodes
$i$ and $j$.

For a matrix $T$ to represent the number of jumps of a chain of states
$(X_{1},X_{2},\dots,X_{n+1})$ the following constraints must be satisfied: (i)
the total number of jumps is equal to the total length of the chain minus one,
$\sum_{ij}t_{ij}=n$, (ii) all jumps can be temporally arranged like domino
tiles $(1,X_{2}),(X_{2},X_{3}),\dots,(X_{n},X_{n+1})$ reflecting the fact that
if at time $\ell$ the Markov chain jumps to state $i$, then at time $\ell+1$
it has to start from state $i$. Constraints (i) and (ii) do not make the
change of variables one to one – there can be many instances of the Markov
chain that correspond to the same set of $t_{ij}$s. In fact, the variables
$t_{ij}$ do not carry any information regarding the temporal order of the
jumps. In other words, given an instance of $T$ we have to count in how many
ways we can order the jumps as $(1,X_{2}),(X_{2},X_{3}),\dots,(X_{n},X_{n+1})$
so that $\sum_{\ell=1}^{n}\delta_{X_{\ell},i}\delta_{X_{\ell+1},j}=t_{ij}$ and
$\\{X_{1},X_{2},\dots,X_{n+1}\\}$ realises a walk in $\mathbf{G}$: we call
this number $\Theta_{T}$. Hence, we can express $Z_{N,n}$ as

$Z_{N,n}(s)=\sum_{t_{11}=0}^{n}\dots\sum_{t_{ij}=0}^{n}\dots\sum_{t_{NN}=0}^{n}\delta_{\sum_{ij}t_{ij},n}\
\Theta_{T}\prod_{i,j}\left(\Pi_{ij}^{t_{ij}}e^{s_{ij}t_{ij}}\right)\ .$ (21)

Figure 1: Flowchart that summarizes the computation of $\Theta_{T}$. From top-
left, we start from an instance of the matrix of jumps $T$. We first check if
the total number of jumps is $n$ and if the multi-graph $\mathbf{M}_{T}$
associated with $T$ is connected. Then, we proceed by checking if Kirchhoff
law is satisfied. There are two possible positive scenarios: i) it is
satisfied in every node – cycle scenario; ii) it is satisfied in every node
except node $1$ and another node which we call $F$: in the former node there
is one more outgoing link, while in the latter there is one more incoming link
– path scenario. In either of these cases, we can relate $\Theta_{T}$ to the
number of T-eulerian cycles (or paths) in $\mathbf{W}_{T}$, which can be
computed knowing the number of eulerian cycles (or paths) in $\mathbf{M}_{T}$.

We now face the problem of computing $\Theta_{T}$. We notice that for many
instances of $T$, this number is simply zero: this is because the
aforementioned domino-like constraint (ii) imposes stringent conditions on the
form of $T$. First of all, the set of edges $(i,j)$, for which $t_{ij}>0$,
together with the union of all their extremes $i$ and $j$ must form a
connected graph. This is because the Markov chain starting from a node $i$ can
only hop to neighbours of $i$ according to the connectivity of $G$.
Mathematically, this condition is equivalent to requiring that the dimension
of the kernel of the Laplacian $L=D_{\mathrm{in}}-T$ is $1$ [42], where
$D_{\mathrm{in}}$ is a diagonal matrix with elements
$(D_{\mathrm{in}})_{ii}=\sum_{j}t_{ji}$. Second, the number of times a Markov
chain jumps towards a state $i$ have to be related to the number of jumps
starting from that state $i$, a phenomenon analogous to the Kirchhoff law in
electric circuits that encodes the global balance of the dynamics. We hereby
distinguish two possible scenarios in which these conditions on $T$ are
satisfied. In the first one, for every state the incoming flux and outgoing
flux are equal, that is $\sum_{j}t_{ij}=\sum_{j}t_{ji}$: this situation
corresponds to a Markov chain starting and ending in the same node, and we
will refer to this as the cycle scenario. In the second one, for all but two
states the incoming and outgoing fluxes are equal. The two special states are
the initial, that we set to $1$ choosing $\mathbb{P}(X_{1})=\delta_{X_{1},1}$,
and the final, $F$, for which one must have $\sum_{j}t_{1j}=1+\sum_{j}t_{j1}$
and $1+\sum_{j}t_{Fj}=\sum_{j}t_{jF}$: we will refer to this as the path
scenario. This leads to a natural way to express

$\displaystyle\Theta_{T}=\Theta_{T}^{\text{path}}\sum_{F\neq
1}\left(\prod\limits_{i\neq 1,F}^{N}\delta_{\sum_{j=1}^{N}t_{ij},\
\sum_{j=1}^{N}t_{ji}}\right)\delta_{\sum_{j=1}^{N}t_{1j}+1,\
\sum_{j=1}^{N}t_{j1}}\delta_{\sum_{j=1}^{N}t_{Fj},\ \sum_{j=1}^{N}t_{jF}+1}+$
$\displaystyle+\Theta_{T}^{\text{cycle}}\prod\limits_{i=1}^{N}\delta_{\sum_{j=1}^{N}t_{ij},\
\sum_{j=1}^{N}t_{ji}}\left(1-\delta_{\sum_{j=1}^{N}t_{1j},0}\right)\ ,$ (22)

where $\Theta_{T}^{\text{path}}$ ($\Theta_{T}^{\text{cycle}}$) is the number
of distinct permutations of the set of $t_{ij}$s in the path (cycle) scenario
which give a realisation of a walk in $\mathbf{G}$ and the deltas enforce
Kirchhoff law. The factor $1-\delta_{\sum_{j=1}^{N}t_{1j},0}$ ensures that the
cycle will pass at least once from node $1$: this condition is required
because the starting node is node $1$. We will show that it is not necessary
to enforce connectedness explicitly because the expressions for
$\Theta_{T}^{\text{path}}$ and $\Theta_{T}^{\text{cycle}}$ are automatically
zero when $T$ is not connected.

We now note that we can interpret the matrix $T$ as the adjacency matrix of a
directed multi-graph $\mathbf{M}_{T}$ with $t_{ij}$ directed links, having
unitary weight, between nodes $i$ and node $j$ – see Fig. 2(b). A directed
multi-graph is a collection of nodes and directed links, in which multiple
links between two nodes are permitted. We refer to the collection of links
between two nodes as a multi-link. As a preliminary step in the computation of
$\Theta_{T}^{\text{path}}$ ($\Theta_{T}^{\text{cycle}}$), we consider a
related combinatorial problem, that is counting how many paths there are on
$\mathbf{M}_{T}$ that start in $1$ and end in $F$ (cycles that start in $1$)
and pass through every link exactly once. We can interpret this as the number
of non-distinct ways we can arrange the jumps like domino tiles that respect
the matrix of jumps $T$. This number overestimates $\Theta_{T}^{\text{path}}$
(respectively $\Theta_{T}^{\text{cycle}}$). To see this, we can consider a
multi-link in $\mathbf{M}_{T}$ having at least two links $l_{1}$ and $l_{2}$.
Given a path (or cycle) that passes through every link in $\mathbf{M}_{T}$, we
can, for instance, generate another distinct one by swapping the order in
which we visit $l_{1}$ and $l_{2}$. This new path (cycle) will not contribute
to $\Theta_{T}^{\text{path}}$ ($\Theta_{T}^{\text{cycle}}$), as the time-
ordered jumps $(1,X_{2}),(X_{2},X_{3}),\dots,(X_{n},X_{n+1})$ are unaffected
by the swap. Nonetheless, this calculation is a useful starting point as we
can compute this number using results available in the literature [43], and we
will show how to correct this overcounting later on.

(a)(b)(c)

Figure 2: (a) State-space connectivity $\mathbf{G}$, un-directed and un-
weighted, with adjacency matrix $A$; (b) directed un-weighted multi-graph
$\mathbf{M}_{T}$, with adjacency matrix $T$, a multi-link from $i$ to $j$ is
composed by $t_{ij}$ links; (c) directed weighted graph $\mathbf{W}_{T}$, with
adjacency matrix $T$, the boldness of links is proportional to the integer
weights $t_{ij}$.

###  III.2 Computation of $\Theta_{T}^{\text{path}}$ and
$\Theta_{T}^{\text{cycle}}$

So far we have described the key steps that underlie our approach, which are
also summarised in the flowchart in Fig. 1. In the following instead we will
provide the details of the calculation. For this reason, some definitions will
be useful and we collect them in this paragraph. An eulerian multi-graph is a
multi-graph for which, at every node $i$, in-degree and out-degree are the
same, viz. $k_{i}^{in}=k_{i}^{out}$. Noticing that in $\mathbf{M}_{T}$ we have
$k_{i}^{in}=\sum_{j}t_{ji}$ and $k_{i}^{out}=\sum_{j}t_{ij}$, it follows from
Kirchhoff law that $\mathbf{M}_{T}$ is either an eulerian multi-graph (for the
cycle scenario) or close to an eulerian multi-graph (for the path scenario),
in the sense that only the initial and final nodes do not satisfy
$k_{in}=k_{out}$. An eulerian cycle (path) on a multi-graph $\mathbf{M}_{T}$,
as already mentioned in Section II, is a cycle (path) that passes through
every link exactly once. We denote the number of eulerian cycles (paths) with
$ec(\mathbf{M}_{T})$ ($ep(\mathbf{M}_{T})$). Furthermore, we will indicate by
$ec(\mathbf{M}_{T}|\chi)$ ($ep(\mathbf{M}_{T}|\chi)$) the number of eulerian
cycles (paths) given some specified condition $\chi$, that in our case will be
a combination of the starting node $1$, the final node $F$, the starting edge
$e_{1}$ and the final edge $e_{F}$.

In the literature on the topic, a result is known for the number
$ec(\mathbf{M}_{T}|e_{1})$ of eulerian cycles of an eulerian multi-graph
$\mathbf{M}_{T}$ with a fixed starting edge. This goes by the name of BEST
theorem [43, 44, 45] and reads

$ec(\mathbf{M}_{T}|e_{1})=\Omega_{w}(\mathbf{M}_{T})\prod_{i=1}^{N}(k^{in}_{i}-1)!,$
(23)

where $\Omega_{w}(\mathbf{M}_{T})$ is the number of arborescences, i.e.,
spanning trees rooted in a node $w$ such that there exists a unique path from
every vertex of $\mathbf{M}_{T}$ to $w$. We note that $\Omega_{w}$ does not
depend on the choice of root $w$ when $\mathbf{M}_{T}$ is an eulerian multi-
graph, so that $\Omega_{w}(\mathbf{M}_{T})=\Omega(\mathbf{M}_{T})$ [44, 46].
Similarly, the r.h.s. of (23) does not show any explicit dependence on the
starting edge $e_{1}$ because of the inherent symmetry in $\mathbf{M}_{T}$. An
explicit expression for $\Omega(\mathbf{M}_{T})$ is given by

$\Omega(\mathbf{M}_{T})=\det(L_{w})\ ,$ (24)

where det is the determinant operator and $L_{w}$ is a submatrix of the
Laplacian of the eulerian multi-graph $\mathbf{M}_{T}$ obtained by removing
(any) $w$-th row and column, a result known in the literature as Tutte’s
theorem or Matrix tree theorem.

In the following, we first consider the path scenario. In this case, the
multi-graph $\mathbf{M}_{T}$ is not eulerian, but we can make it so simply by
adding a link $e_{F}$ from $F$ to $1$. We refer to this modified graph as
$\tilde{\mathbf{M}}_{T}$. Using BEST theorem we have

$ec(\tilde{\mathbf{M}}_{T}|e_{1})=\Omega(\tilde{\mathbf{M}}_{T})\prod_{i\neq
1}^{N}\left(\sum_{j=1}^{N}t_{ji}-1\right)!\left(\sum_{j=1}^{N}t_{j1}\right)\
,$ (25)

where we use the fact that the in-degree of node $1$ is
$\sum_{j=1}^{N}t_{j1}+1$ in $\tilde{\mathbf{M}}_{T}$. The number of eulerian
cycles starting from node $1$ is related to $ec(\tilde{\mathbf{M}}_{T}|e_{1})$
by
$ec(\tilde{\mathbf{M}}_{T}|1)=ec(\tilde{\mathbf{M}}_{T}|e_{1})\left(\sum_{j}t_{j1}+1\right)$.
Furthermore, $ep(\mathbf{M}_{T}|1,F)$ is equal to the number of eulerian
cycles in $\tilde{\mathbf{M}}_{T}$ starting in $1$ and ending with the link we
added to construct it, viz. $ec(\tilde{\mathbf{M}}_{T}|1,e_{F})$. This number
can be computed by considering an eulerian cycle in $\mathbf{M}_{T}$ as a
collection of loops passing through node $1$. The number of these loops is
given by the in-degree of node $1$, so that we have
$ep(\mathbf{M}_{T}|1,F)=ec(\tilde{\mathbf{M}}_{T}|1)/(\sum_{j}t_{j1}+1)$. All
these considerations put together give

$ep(\mathbf{M}_{T}|1,F)=\Omega_{1}(\mathbf{M}_{T})\prod_{i\neq
1}^{N}\left(\sum_{j=1}^{N}t_{ji}-1\right)!\left(\sum_{j=1}^{N}t_{j1}\right)\ $
(26)

where we used
$\Omega(\tilde{\mathbf{M}}_{T})=\text{det}(L_{1})=\Omega_{1}(\mathbf{M}_{T})$
with $L_{1}$ the cofactor of the graph Laplacian $L$ obtained by removing the
first row and column. We note that while $\Omega(\tilde{\mathbf{M}}_{T})$ does
not depend on the choice of the root since $\tilde{\mathbf{M}}_{T}$ is
eulerian, $\Omega_{1}(\mathbf{M}_{T})$ does, because $\mathbf{M}_{T}$ is not
eulerian.

We now consider the cycle scenario. In this case, since $\mathbf{M}_{T}$ is
already an eulerian graph, we can readily express $ec(\mathbf{M}_{T},1)$ as

$ec(\mathbf{M}_{T}|1)=ec(\mathbf{M}_{T}|e_{1})\sum_{j=1}^{N}t_{j1}=\Omega(\mathbf{M}_{T})\prod_{i\neq
1}^{N}\left(\sum_{j=1}^{N}t_{ji}-1\right)!\left(\sum_{j=1}^{N}t_{j1}\right)\
.$ (27)

As previously argued, $ep(\mathbf{M}_{T}|1,F)$ ($ec(\mathbf{M}_{T}|1)$)
overestimates $\Theta_{T}^{\text{path}}$ ($\Theta_{T}^{\text{cycle}}$). To
correct this, one must consider all the links belonging to a given multi-link
as totally equivalent. This boils down to considering a weighted graph
$\mathbf{W}_{T}$ (see 2(c)) in place of the multi-graph $\mathbf{M}_{T}$. The
weighted graph $\mathbf{W}_{T}$ has adjacency matrix $T$ and directed links
(e.g., between nodes $i$ and $j$) obtained by merging all the multi-links
(between $i$ and $j$) in $\mathbf{M}_{T}$ together. In analogy with [45], we
define the notion of $T$-eulerian cycle (path) as a cycle (path) that passes
through every link $(i,j)$ a number $t_{ij}$ of times. With an abuse of
notation, we denote the number of $T$-eulerian cycles (paths) by
$ec(\mathbf{W}_{T}|\chi)$ ($ep(\mathbf{W}_{T}|\chi)$), as it will be clear by
the graph we are considering whether we are referring to eulerian or
$T$-eulerian cycles (paths). Crucially, in the cycle scenario
$\Theta_{T}^{\text{cycle}}$ is equal to the number of $T$-eulerian cycles
starting from node $1$, i.e., $ec(\mathbf{W}_{T}|1)$ in $\mathbf{W}_{T}$,
whereas in the path scenario $\Theta_{T}^{\text{path}}$ is equal to the number
of $T$-eulerian paths from $1$ to $F$, i.e., $ep(\mathbf{W}_{T}|1,F)$ in
$\mathbf{W}_{T}$. The combinatorial factor connecting $ep(\mathbf{W}_{T}|1,F)$
($ec(\mathbf{W}_{T}|1)$ ) to $ep(\mathbf{M}_{T}|1,F)$ ($ec(\mathbf{M}_{T}|1)$)
is simply the number of permutations of links in a multi-link for every multi-
link in $\mathbf{M}_{T}$

$\displaystyle\begin{split}&ep(\mathbf{M}_{T}|1,F)=ep(\mathbf{W}_{T}|1,F)\prod_{i,j=1}^{N}t_{ij}!\\\
&ec(\mathbf{M}_{T}|1)=ec(\mathbf{W}_{T}|1)\prod_{i,j=1}^{N}t_{ij}!\
.\end{split}$ (28)

This allows us to write explicit expressions for $\Theta_{T}^{\text{path}}$
and $\Theta_{T}^{\text{cycle}}$

$\Theta_{T}^{\text{path}}=\Theta_{T}^{\text{cycle}}=\det(L_{1})\prod_{i=1}^{N}\frac{(\sum_{j=1}^{N}t_{ji}-1)!}{\prod_{j=1}^{N}t_{ij}!}\sum_{k=1}^{N}t_{k1}\
,$ (29)

where we recall that $L_{1}$ is – in both the cycle and the path scenarios –
the submatrix of the graph Laplacian $L$ obtained by removing the first row
and column. We remark that, although the expressions for
$\Theta_{T}^{\text{path}}$ and $\Theta_{T}^{\text{cycle}}$ in (29) are
formally the same, the variable $T$ is of different nature in the path and
cycle scenario as it satisfies different sets of constraints. With this
expression, we can write the moment generating function explicitly as

$\displaystyle\begin{split}Z_{N,n}(s)=&\sum_{t_{11}=0}^{n}\dots\sum_{t_{ij}=0}^{n}\dots\sum_{t_{NN}=0}^{n}\delta_{\sum_{ij}t_{ij},n}\
\prod_{i,j}\left(\Pi_{ij}^{t_{ij}}e^{s_{ij}t_{ij}}\right)\det(L_{1})\sum_{j=1}^{N}t_{j1}\prod_{i=1}^{N}\frac{(\sum_{j=1}^{N}t_{ji}-1)!}{\prod_{j=1}^{N}t_{ij}!}\\\
&\Bigg{(}\sum\limits_{F\neq 1}\left(\prod\limits_{i\neq
1,F}^{N}\delta_{\sum_{j=1}^{N}t_{ij},\
\sum_{j=1}^{N}t_{ji}}\right)\delta_{\sum_{j=1}^{N}t_{1j}+1,\
\sum_{j=1}^{N}t_{j1}}\delta_{\sum_{j=1}^{N}t_{Fj},\
\sum_{j=1}^{N}t_{jF}+1}+\\\
&\left(\prod\limits_{i=1}^{N}\delta_{\sum_{j=1}^{N}t_{ij},\
\sum_{j=1}^{N}t_{ji}}\right)\left(1-\delta_{\sum_{j}t_{1j},0}\right)\Bigg{)}\
.\end{split}$ (30)

We note that the factor $\det(L_{1})$ kills configurations of $T$ that have
the null-space dimension of the Laplacian greater than $1$. This ensures that
we only consider graphs $\mathbf{M}_{T}$ – equivalently, $\mathbf{W}_{T}$ –
that are connected, as it is known in the literature that the dimension of the
null-space of the graph Laplacian is the number of connected components of a
graph [42].  Remarkably, in equation (30) the contributions for paths and
cycles are split, giving an interesting physical perspective. In general, this
difference is more pronounced when $n$ is small, in particular when the walker
has not explored the full state space. In the limit of large $n$,
contributions relative to paths and cycles are comparable and share the same
asymptotics, as we show in the next Section.

Compared to the spectral method to compute the moment generating function
[17], which requires the computation of all eigenvalues and eigenvectors of an
$N\times N$ matrix, our formula is computationally favourable when $n$ is
small and $N$ is large. If $n$ is large, instead, the spectral method is
numerically more efficient.

###  III.3 Long-time asymptotics

Expression (30) is valid for every finite $n$, and can be used to derive the
large $n$ limit and, in principle, finite $n$ corrections. In the following,
we focus on the large deviation regime, which corresponds to taking $n$ to be
much greater than the longest relaxation time of the system $\tau(N)$,
$n\gg\tau(N)$. In this limit it is useful to rescale time-additive variables
with $n$ as in (17), as we can approximate the sums over $t_{11},\dots,t_{NN}$
with integrals

$\frac{1}{n^{|E_{\mathbf{W}_{T}}|}}\sum_{t_{11}}\dots\sum_{t_{NN}}\
\to\;\prod_{i,j}\int_{0}^{1}\\!\mathrm{d}\nu_{ij}\ ,$ (31)

where $|E_{\mathbf{W}_{T}}|$ is the number of directed edges in the weighted
graph $\mathbf{W}_{T}$ and $\nu_{ij}$s are defined as in (17). In the r.h.s.
of (31) and in the following, by $\sum_{ij}$ and $\prod_{ij}$ we mean sums and
products over $(i,j)$ such that $(i,j)$ is a directed link in
$\mathbf{W}_{T}$. To leading order in $n$ we obtain the following asymptotic
expressions

$\displaystyle\prod_{i,j}\left(\Pi_{ij}^{t_{ij}}e^{s_{ij}t_{ij}}\right)\
\to\;e^{n\sum_{i=1}^{N}\sum_{j=1}^{N}\left(s_{ij}+\log\Pi_{ij}\right)\
\nu_{ij}}$ (32)
$\displaystyle\left(\sum_{j=1}^{N}t_{j1}\right)\left(\prod_{i=1}^{N}\frac{(\sum_{j=1}^{N}t_{ji}-1)!}{\prod_{j=1}^{N}t_{ij}!}\right)\
\to\;e^{n\sum_{i=1}^{N}\sum_{j=1}^{N}\nu_{ij}(\log(\sum_{k=1}^{N}\nu_{ik})-\log(\nu_{ij}))}$
(33) $\displaystyle\delta_{\sum_{ij}t_{ij},n}\
\to\;\delta_{\sum_{ij}\nu_{ij},1}$ (34) $\displaystyle
n^{|E_{\mathbf{W}_{T}}|}\ \to\;e^{|E_{\mathbf{W}_{T}}|\ \log n}\ .$ (35)

The Kirchhoff constraints tend to the same form for large $n$, giving
explicitly

$\displaystyle\sum_{F\neq 1}\left(\prod_{i\neq
1,F}\delta_{\sum_{j=1}^{N}t_{ji},\sum_{j=1}^{N}t_{ij}}\right)\delta_{\sum_{j=1}^{N}t_{j1}+1,\sum_{j=1}^{N}t_{1j}}\delta_{\sum_{j=1}^{N}t_{jF},1+\sum_{j=1}^{N}t_{Fj}}\
\to\;(N-1)\prod_{i=1}^{N}\delta_{\sum_{j=1}^{N}\nu_{ji},\sum_{j=1}^{N}\nu_{ij}}$
(36)
$\displaystyle\prod_{i=1}^{N}\delta_{\sum_{j=1}^{N}t_{ji},\sum_{j=1}^{N}t_{ij}}\to\;\prod_{i=1}^{N}\delta_{\sum_{j=1}^{N}\nu_{ji},\sum_{j=1}^{N}\nu_{ij}}\
.$ (37)

We also notice that

$\text{det}(L_{1})=n^{N-1}\det\left(\frac{L_{1}}{n}\right)=e^{(N-1)\log
n+\text{Tr}\left[\log\left(\frac{L_{1}}{n}\right)\right]}\ \to\;e^{(N-1)\log
n}\ ,$ (38)

where we use the fact that the determinant is multi-linear in the rows and
that each element in $L_{1}$ is proportional to $n$ by construction, so that
$\text{Tr}\left[\log\left(\frac{L_{1}}{n}\right)\right]$ is finite for large
$n$. Remarkably $\text{det}(L_{1})$ becomes sub-leading in the large $n$
limit, while it may be an interesting term to study the finite-time transient
behaviour of the Markov chain. Finally, the term $\delta_{\sum_{j}t_{1j},0}$
present in the factor $1-\delta_{\sum_{j}t_{1j},0}$ becomes negligible for
large $n$.

Putting all together, we obtain to exponential leading order in $n$

$Z_{N,n}(s)\approx\int_{0}^{1}\dots\int_{0}^{1}\\!\left(\prod_{i,j}\mathrm{d}\nu_{ij}\right)e^{n\left[\sum_{ij}\nu_{ij}(\log(\sum_{k}\nu_{ik})-\log(\nu_{ij}))+\sum_{ij}\left(s_{ij}+\log\Pi_{ij}\right)\
\nu_{ij}\right]}\left(\prod_{i}\delta_{\sum_{j}\nu_{ji},\sum_{j}\nu_{ij}}\right)\delta_{\sum_{ij}\nu_{ij},1}\
.$ (39)

We note that the integrand in (39) can be brought to the form
$e^{n\lambda_{N}[\nu]}$, with the following definitions:

$\displaystyle\lambda_{N}$
$\displaystyle[\nu]=\lambda_{1}[\nu]+\lambda_{2}[\nu]+\lambda_{3}[\nu]+\lambda_{4}[\nu]$
(40) $\displaystyle\lambda_{1}$
$\displaystyle[\nu]=\sum_{i=1}^{N}\sum_{j=1}^{N}\nu_{ij}\left(\log\left(\sum_{k=1}^{N}\nu_{ik}\right)-\log(\nu_{ij})\right)$
(41) $\displaystyle\lambda_{2}$
$\displaystyle[\nu]=\sum_{i=1}^{N}\sum_{j=1}^{N}\log(\Pi_{ij})\ \nu_{ij}$ (42)
$\displaystyle\lambda_{3}$
$\displaystyle[\nu]=\sum_{i=1}^{N}\sum_{j=1}^{N}s_{ij}\ \nu_{ij}$ (43)
$\displaystyle\lambda_{4}$
$\displaystyle[\nu]=\epsilon\left(\sum_{i=1}^{N}\sum_{j=1}^{N}\nu_{ij}-1\right)+\sum_{i=1}^{N}\eta_{i}\left(\sum_{j=1}^{N}\nu_{ij}-\sum_{j=1}^{N}\nu_{ji}\right)\
,$ (44)

where $\epsilon$ and $\eta_{i}$ are Lagrange multipliers fixing the respective
constraints. In (40) each term has a clear physical interpretation:
$\lambda_{1}$, in (41), is the geometric – viz. related to the connectivity of
the graph $\mathbf{G}$ – entropy of a random walk on a graph with nodes and
links contributions, akin to the entropy of a free particle; $\lambda_{2}$, in
(42), is the entropy due to the dynamics, encoded in the transition matrix;
$\lambda_{3}$, in (43), is the tilting potential necessary to drive the system
towards a fluctuation of the pair empirical occupation measure; finally,
$\lambda_{4}$ in (44), enforces the normalisation and Kirchhoff-law (global
balance).

We can calculate the leading order in $n$ of (39) via a saddle-point
approximation, arriving at

$Z_{N,n}(s)\approx e^{n\lambda_{N}[\nu^{*}]}\ ,$ (45)

where $\nu^{*}=\text{argmin}_{\nu,\epsilon,\eta}\ \lambda_{N}[\nu]$ and
$\nu^{*}$ are the minimisers of $\lambda_{N}[\nu]$ with respect to the set of
$\nu_{ij}$s, $\eta_{i}$s and $\epsilon$. From the Euler–Lagrange equations for
critical points of (40), we find the following implicit expression for
$\nu_{ij}^{*}$:

$\nu_{ij}^{*}=\left(\Pi_{s}\right)_{ij}\left(\frac{e^{-\eta_{j}}}{e^{-\epsilon}e^{-\eta_{i}}}\right)\sum_{k=1}^{N}\nu_{jk}^{*}\
,$ (46)

where the tilted matrix introduced in (19) appears. From (46) we can write
self-consistent conditions for $\epsilon$ and $\eta_{i}$ as follows:

$\displaystyle\sum\limits_{j}\left(\Pi_{s}\right)_{ij}e^{-\eta_{j}}$
$\displaystyle=$ $\displaystyle e^{-\epsilon}e^{-\eta_{i}}$ (47)
$\displaystyle\sum\limits_{i}\left(\Pi_{s}\right)_{ij}\frac{\sum_{k}\nu_{ik}^{*}}{e^{-\eta_{i}}}$
$\displaystyle=$ $\displaystyle
e^{-\epsilon}\frac{\sum_{k}\nu_{jk}^{*}}{e^{-\eta_{j}}}\ ,$ (48)

which reveal that $e^{-\epsilon}$ is an eigenvalue of the tilted matrix
$\Pi_{s}$ with right eigenvector components $r_{i}=e^{-\eta_{i}}$ and left
eigenvector components $l_{j}=\sum_{k}\nu_{jk}^{*}/e^{-\eta_{j}}$.
Substituting (47) into (40) we get

$\lambda_{N}[\nu^{*}]=-\epsilon\ ,$ (49)

and, in particular, since $\lambda_{N}[\nu^{*}]$ is a maximum, $e^{\epsilon}$
is the dominant eigenvalue of (19). The same conclusion can be reached by
noticing that the left and right eigenvector elements in (47) and (48) are all
positive, which is true only for the dominant eigenvalue. These arguments
provide a direct link with spectral methods. In particular, (49) provides an
expression for the logarithm of the dominant eigenvalue of the tilted matrix.

Remarkably, this approach also provides an alternative expression for the so-
called driven (or effective) process. This is a modified Markov chain that
explains how specific fluctuations are created in time [36, 37, 1, 20]; under
certain conditions, it is equivalent to the original Markov chain conditioned
to visiting the fluctuation of interest. Useful spectral and variational
expressions of the driven process already appeared in the papers just
mentioned. Here, we offer another explicit variational representation valid
for discrete-time Markov chains. In agreement with [20], the minimisers
$\nu^{*}=\left\\{\nu^{*}_{ij}\right\\}$ of the action functional (40)
characterise the driven process transition matrix with components

$\tilde{\Pi}_{ij}=\frac{\nu^{*}_{ij}}{\sum_{k=1}^{N}\nu^{*}_{ik}}\ .$ (50)

This last expression offers an alternative way to physically study and
simulate the appearance of fluctuations and rare events in discrete-time
Markov chain models.

Concluding, in (40) we have obtained $\lambda_{N}$, the SCGF associated with
the probability distribution of the pair empirical occupation measure in (4).
To get the rate functional (2) we only need to Legendre–Fenchel transform the
SCGF in (40), i.e.,

$\sup_{s}\left(\sum_{i=1}^{N}\sum_{j=1}^{N}s_{ij}\nu^{*}_{ij}-\lambda_{N}[\nu^{*}]\right)=\sup_{s}\left(\lambda_{3}[\nu^{*}]-\lambda_{N}[\nu^{*}]\right)=-\lambda_{1}[\nu^{*}]-\lambda_{2}[\nu^{*}]-\lambda_{4}[\nu^{*}]=H[\nu^{*}]\
,$ (51)

where in the last step we recognise the pair empirical rate functional (with
the necessary constraints – mentioned and understood in (2) – fixed by the
Lagrange multipliers in $\lambda_{4}$).

Assuming that one is interested in studying large fluctuations of an
observable of the form (1), we remark that the associated SCGF can be obtained
simply replacing $\lambda_{3}[\nu]$ in (43) with

$\lambda_{3}[\nu]=s\sum_{i,j=1}^{N}f(i,j)\nu_{ij}\ ,$ (52)

where $s$ is the tilting parameter conjugated to $C_{n}$. For instance, in
physics applications, it is often of interest to consider the empirical
current $\mathbb{J}_{n}(i,j)=L_{n}^{(2)}(i,j)-L_{n}^{(2)}(j,i)$, viz. the
antisymmetric part of the pair-empirical occupation measure in (4), or again
the occupation measure itself $L_{n}(i)=\sum_{j=1}^{N}L_{n}^{(2)}(i,j)$. The
empirical current is an important observable as it allows us to estimate how
far a system lies from equilibrium, whereas the occupation measure gives an
estimate of the time spent by the system in each state of the state space.

##  IV Two-state model

In this Section, in order give a more pedagogical understanding of how one
could use (30) to derive leading, i.e., the SCGF in (40), and finite $n$
behaviour, we compare our method with the more standard spectral approach on a
simple two-state Markov chain. We show that the two methods give equivalent
results and propose a physical interpretation of all terms appearing in the
SCGF. We consider a general two-state Markov chain, whose transition matrix
$\Pi$ reads

$\Pi=\begin{pmatrix}1-p&p\\\ q&1-q\end{pmatrix}\ ,$ (53)

with $p$ and $q$ between $0$ and $1$. We choose to observe the flux between
node $1$ and node $2$, that is

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}\delta_{X_{\ell},1}\delta_{X_{\ell+1},2}=\frac{t_{12}}{n}\
.$ (54)

The long-time behaviour of $C_{n}$ is given by
$\lim_{n\rightarrow\infty}C_{n}=\frac{pq}{p+q}\eqqcolon c^{*}$. Intuitively,
when $t_{12}$ is large, the Markov chain jumps frequently from $1$ to $2$ and
from $2$ to $1$; instead, when $t_{12}$ is small the chain spends most of the
time jumping from $1$ to $1$ and/or from $2$ to $2$. This situation is
reminiscent of a particle in a double well potential immersed in a thermal
bath, where temperature—that is, the strength of noise—regulates the frequency
of jumps between the two minima. In this two-state model, the tilting
parameter $s$ plays a role analogous to the temperature.

###  IV.1 Spectral approach

The moment generating function can be computed using spectral methods. We
start from (20) (restricted to the case of the observable (54)), i.e.,

$Z_{N,n}(s)=\sum_{X_{1},\dots,X_{n+1}}\mathbb{P}(X_{1})\prod_{\ell=1}^{n}\Pi_{X_{\ell},X_{\ell+1}}e^{s\delta_{X_{\ell},1}\delta_{X_{\ell+1},2}}\
,$ (55)

which can be cast in the form

$Z_{N,n}(s)=\langle\mathbb{P}_{1}|\left(\Pi_{s}\right)^{n}|1\rangle\ ,$ (56)

where $\langle\mathbb{P}_{1}|=(1,0)$ is the vector of initial probabilities,
$|1\rangle=(1,1)$ and $\Pi_{s}$ is the tilted matrix, viz. (19) restricted to
the case at hand, which reads

$\Pi_{s}=\begin{pmatrix}1-p&pe^{s}\\\ q&1-q\end{pmatrix}\ .$ (57)

We can use the spectral decomposition of $\Pi_{s}$ to get

$Z_{2,n}(s)=\langle P_{0}|\left(|r^{+}\rangle\langle
l^{+}|\Lambda_{+}^{n}+|r^{-}\rangle\langle
l^{-}|\Lambda_{-}^{n}\right)|1\rangle=r_{1}^{+}\Lambda_{+}^{n}(l_{1}^{+}+l_{2}^{+})+r_{1}^{-}\Lambda_{-}^{n}(l_{1}^{-}+l_{2}^{-})\
,$ (58)

where $\Lambda_{\pm}$ are the eigenvalues of $\Pi_{s}$ and $l^{\pm}$,
$r^{\pm}$ the corresponding left and right eigenvectors, respectively. We
notice that—for the spectral decomposition of $\Pi_{s}$ to be valid—left and
right eigenvectors have to be bi-orthonormal.

Figure 3: From top left to bottom right: Cycles and Paths contributions to the
moment generating function, their sum (Comb.) obtained using the graph-
combinatorial approach, in comparison with the moment generating function
obtained via the spectral decomposition (Spect.) for increasing values of $n$.
Comb. and Spect. curves fully overlap and the green-ish color shown is
obtained by the combination of blue (Comb.) and yellow (Spect.). Furthermore,
as expected, the Cycles contribution to the moment generating function is
smaller for $s>0$ and larger for $s<0$ than the Paths one. To generate these
plots, we have used $p=q=1/2$.

By computing the eigenvalues and eigenvectors of $\Pi_{s}$ explicitly, we
arrive at

$\displaystyle\begin{split}&Z_{2,n}(s)=\frac{1}{2^{n+1}\sqrt{(p-q)^{2}+4pqe^{s}}}\times\\\
&\hskip
28.45274pt\times\Bigg{\\{}\left((1-2e^{s})p-q\right)\left[\left(2-p-q-\sqrt{(p-q)^{2}+4pqe^{s}}\right)^{n}-\left(2-p-q+\sqrt{(p-q)^{2}+4pqe^{s}}\right)^{n}\right]\\\
&\hskip
28.45274pt\left[\sqrt{(p-q)^{2}+4pqe^{s}}\left(2-p-q-\sqrt{(p-q)^{2}+4pqe^{s}}\right)^{n}+\left(2-p-q+\sqrt{(p-q)^{2}+4pqe^{s}}\right)^{n}\right]\Bigg{\\}}\
.\end{split}$ (59)

###  IV.2 Graph-combinatorial approach

The moment generating function can also be computed using (30). Remarkably,
this other formulation highlights two different contributions coming from
cycles and paths traveled starting from state $1$ of the state space. These
can explicitly be written as

$\displaystyle Z_{2,n}^{\text{Cycles}}(s)$
$\displaystyle=(1-p)^{n}+\sum_{t_{11}=0}^{n}\sum_{t_{12}=1}^{n}\sum_{t_{22}=0}^{n}\delta_{t_{11}+2t_{12}+t_{22},n}\left(pqe^{s}\right)^{t_{12}}(1-p)^{t_{11}}(1-q)^{t_{22}}{t_{11}+t_{12}\choose
t_{11}}{t_{22}+t_{12}-1\choose t_{22}}$ (60) $\displaystyle
Z_{2,n}^{\text{Paths}}(s)$
$\displaystyle=\sum_{t_{11}=0}^{n}\sum_{t_{12}=1}^{n}\sum_{t_{22}=0}^{n}\delta_{t_{11}+2t_{12}+t_{22}-1,n}\frac{\left(pqe^{s}\right)^{t_{12}}}{q}(1-p)^{t_{11}}(1-q)^{t_{22}}{t_{11}+t_{12}-1\choose
t_{11}}{t_{22}+t_{12}-1\choose t_{22}}\ ,$ (61)

where in $Z_{2,n}^{\text{Cycles}}(s)$ we made explicit the cycle contribution
coming from staying for $n$ consecutive steps on state $1$, and in
$Z_{2,n}^{\text{Paths}}(s)$ the counting needs to start from $t_{12}=1$
because to have a meaningful path contribution the Markov chain needs to hop
at least once from state $1$ to state $2$. We remark that in this simple model
$\det(L_{1})=t_{12}$ if $t_{12}\neq 0$ (in such a case the term is absorbed in
(60) and (61) by the binomial coefficients), while when $t_{12}=0$ the
Laplacian is a $1\times 1$ matrix: $L_{1}$ is thus an empty matrix, and we
take its determinant to be $1$ for consistency.

We can find an explicit expression for (60) and (61) analytically. We replace
the delta functions appearing by their contour integral representations

$\delta_{i,j}=\frac{1}{2\pi i}\oint_{|z|=1}z^{i-j-1}\mathrm{d}z\ .$ (62)

After making the substitution, we notice that the integrands in (60) and (61)
are analytic functions everywhere except in $0$. This allows us to deform the
integration contour to a circle of radius $\epsilon\ll 1$. The reason for this
is to avoid spurious poles in the following steps.

We now let all the sums run up to $\infty$. This procedure is allowed as
higher order terms in the sums do not affect the residue in $0$. The infinite
sums can be explicitly evaluated and by doing so we get

$\displaystyle Z_{2,n}^{\text{Cycles}}(s)$
$\displaystyle=(1-p)^{n}+\frac{1}{(e^{s}pq-pq+q+p-1)}\frac{1}{2\pi
i}\oint_{|z|=\epsilon}\frac{1}{z^{n-1}}\frac{e^{s}pq}{(z(1-p)-1)(z-z^{*}_{1})(z-z^{*}_{2})}$
(63) $\displaystyle Z_{2,n}^{\text{Paths}}(s)$
$\displaystyle=-\frac{1}{(e^{s}pq-pq+q+p-1)}\frac{1}{2\pi
i}\oint_{|z|=\epsilon}\frac{1}{z^{n}}\frac{e^{s}p}{(z-z^{*}_{1})(z-z^{*}_{2})}\
,$ (64)

where

$\displaystyle z^{*}_{1}$
$\displaystyle=\frac{2}{2-p-q-\sqrt{(p-q)^{2}+4e^{s}pq}}$ (65) $\displaystyle
z^{*}_{2}$ $\displaystyle=\frac{2}{2-p-q+\sqrt{(p-q)^{2}+4e^{s}pq}}\ .$ (66)

Notice that $z^{*}_{1}$ and $z^{*}_{2}$ are exactly the inverse of the
eigenvalues found with spectral methods. We remark that the integrands in (63)
and (64) have acquired new singularities, in the form of simple poles at
$z^{*}_{1}$, $z^{*}_{2}$ and $1/(1-p)$: these poles are unphysical, in the
sense that their residue should not be considered when computing the contour
integrals.

Figure 4: Ratios $\frac{Z_{2,n}^{\text{Cycles}}}{Z_{2,n}^{\text{Paths}}}$ as a
function of $n$ for five different values of $s$, which are, from top to
bottom: cyan, $s=-1.0$; magenta, $s=-0.5$; black, $s=0$; red, $s=0.5$; orange,
$s=1.0$

We can express $Z_{2,n}^{\text{Cycles}}(s)$ and $Z_{2,n}^{\text{Paths}}(s)$ as

$\displaystyle Z_{2,n}^{\text{Cycles}}(s)$
$\displaystyle=(1-p)^{n}+\frac{1}{(e^{s}pq-
pq+q+p-1)}\text{Res}_{z=0}\left(\frac{1}{z^{n-1}}\frac{e^{s}pq}{(z(1-p)-1)(z-z^{*}_{1})(z-z^{*}_{2})}\right)$
(67) $\displaystyle Z_{2,n}^{\text{Paths}}(s)$
$\displaystyle=-\frac{1}{(e^{s}pq-
pq+q+p-1)}\text{Res}_{z=0}\left(\frac{1}{z^{n}}\frac{e^{s}p}{(z-z^{*}_{1})(z-z^{*}_{2})}\right)\
.$ (68)

Computing the residues we find

$\displaystyle Z_{2,n}^{\text{Cycles}}(s)$ $\displaystyle=(1-p)^{n}+$ (69)
$\displaystyle+\frac{(e^{s}pq(z_{2}^{*})^{-n}(z_{1}^{*}((z_{2}^{*}/z_{1}^{*})^{n}-(z_{2}^{*}-pz_{2}^{*})^{n})+z_{2}^{*}(-1+(z_{2}^{*}-pz_{2}^{*})^{n}+(-1+p)z_{1}^{*}(-1+(z_{2}^{*}/z_{1}^{*})^{n}))))}{(e^{s}pq-
pq+p+q-1)((1+(-1+p)z_{1}^{*})(z_{1}^{*}-z_{2}^{*})(1+(-1+p)z_{2}^{*}))}$ (70)
$\displaystyle Z_{2,n}^{\text{Paths}}(s)$
$\displaystyle=\frac{(e^{s}p(z_{2}^{*})^{-n}(-1+(z_{2}^{*}/z_{1}^{*})^{n}))}{(e^{s}pq-
pq+q+p-1)(z_{1}^{*}-z_{2}^{*})}\ .$ (71)

By summing these two contributions and replacing $z^{*}_{1}$ and $z^{*}_{2}$
from (65) and (66), we obtain exactly (59).

In Fig. 3 we show the functions $Z_{2,n}^{\text{Cycles}}$ and
$Z_{2,n}^{\text{Paths}}$ and compare them with the moment generating function
previously obtained via spectral methods.

Evidently, the moment generating function obtained by summing up cycles and
paths contributions completely matches the moment generating function obtained
with spectral methods, as the two curves are indistinguishable. An advantage
of the graph-combinatorial approach with respect to the spectral calculation
is the possibility to split the contributions coming from cycles and paths. As
expected for the simple model investigated, cycles contribute less to the
moment generating function for $s>0$ with respect to paths, and viceversa for
$s<0$. The reason for this is that in the path scenario the Markov chain has
to jump at least once from $1$ to $2$, contributing to $C_{n}$. The larger the
$n$, the less pronounced is this effect. We also show in Fig. 4 the ratio
$Z_{2,n}^{\text{Cycles}}/Z_{2,n}^{\text{Paths}}$ for a few fixed values of the
tilting parameter $s$ as a function of time $n$. Noticeably, the ratios become
constant for $n$ big enough, supporting the fact that both
$Z_{2,n}^{\text{Cycles}}$ and $Z_{2,n}^{\text{Paths}}$ share the same
asymptotics for large $n$ and differ only by a constant prefactor that is a
function of $s$.

###  IV.3 Large deviation regime

We now investigate fluctuations in the large-$n$ limit computing the SCGF
$\lambda(s)$. We compare the spectral and the variational formulae, to
highlight the benefits of both approaches.

Using the spectral approach, the logarithm of the dominant eigenvalue is the
SCGF (59) and reads

$\lambda(s)=\log\frac{2-p-q+\sqrt{(p-q)^{2}+4pqe^{s}}}{2}\ .$ (72)

We can arrive at the same result by minimising (40). Noticing that Kirchhoff
law reduces to $\nu_{12}=\nu_{21}$, the action functional reduces to

Figure 5: Top-left panel: Plot of the minimizers $\nu^{*}_{11}$,
$\nu^{*}_{12}$, $\nu^{*}_{22}$ of the form (74), (75), (76) for $p=q=0.5$ as a
function of $s$. We notice that the curve of $\nu^{*}_{11}$ coincides with
that of $\nu^{*}_{22}$, due to the symmetry of $p$ and $q$. Top-right panel:
Plot of all the contributions to the SCGF as defined in (42), (43), (44),
(77), (78) of the two-state model and their sum as a function of $s$ for
$p=q=0.5$. For this choice of parameters, we notice that the curve of
$\lambda_{1,\text{states}}$ coincides with that of $\lambda_{2}$, and
interestingly they do not depend on $s$. Bottom-left panel: Plot of the
minimizers $\nu^{*}_{11}$, $\nu^{*}_{12}$, $\nu^{*}_{22}$ of the form (74),
(75), (76) for $p=0.5$, $q=0.9$ as a function of $s$. Bottom-right panel: Plot
of all the contributions to the SCGF as defined in (42), (43), (44), (77),
(78) of the two-state model and their sum as a function of $s$ for $p=0.5$ and
$q=0.9$.

$\displaystyle\begin{split}&\lambda[\nu]=(\nu_{11}+\nu_{12})\log(\nu_{11}+\nu_{12})+(\nu_{12}+\nu_{22})\log(\nu_{12}+\nu_{22})-\log{\nu_{11}}-\log{\nu_{22}}\\\
&+\nu_{11}\log(1-p)+\nu_{12}(\log p+\log
q)+\nu_{22}\log(1-q)+s\nu_{12}+\epsilon(\nu_{11}+2\nu_{12}+\nu_{22}-1)\
.\end{split}$ (73)

The determinant of the system of equations satisfied by the minimum of
(LABEL:eq:action2states), which is linear in $\nu$, must be $0$ to have non-
trivial solutions. This condition gives an equation for $e^{-\epsilon}$
($\epsilon$ is the Lagrange multiplier fixing the normalisation condition)
whose solution gives—through (49)—expression (72). This last can be replaced
in the form of the minimisers obtained by solving the linear system, which
read

$\displaystyle\nu^{*}_{11}$
$\displaystyle=\frac{(1-q)^{-1}e^{-\epsilon}-1}{(1-q)^{-1}e^{-\epsilon}+2((1-p)^{-1}e^{-\epsilon}-1)((1-q)^{-1}e^{-\epsilon}-1)+(1-p)^{-1}e^{-\epsilon}-2}$
(74) $\displaystyle\nu^{*}_{12}$
$\displaystyle=\frac{((1-p)^{-1}e^{-\epsilon}-1)((1-q)^{-1}e^{-\epsilon}-1)}{(1-q)^{-1}e^{-\epsilon}+2((1-p)^{-1}e^{-\epsilon}-1)((1-q)^{-1}e^{-\epsilon}-1)+(1-p)^{-1}e^{-\epsilon}-2}$
(75) $\displaystyle\nu^{*}_{22}$
$\displaystyle=\frac{(1-p)^{-1}e^{-\epsilon}-1}{(1-q)^{-1}e^{-\epsilon}+2((1-p)^{-1}e^{-\epsilon}-1)((1-q)^{-1}e^{-\epsilon}-1)+(1-p)^{-1}e^{-\epsilon}-2}\
,$ (76)

to get their explicit form as a function of $p$, $q$, and the tilting
parameter $s$.

In the top-left and bottom-left panels of Fig. 5 we plot the minimisers (74),
(75), and (76) as a function of the tilting parameter $s$. For the top-left
case, we use $p=q=0.5$, while for the bottom case $p=0.5$ and $q=0.9$. We
notice that when $p=q$, as in the top-left panel, $\nu^{*}_{11}=\nu^{*}_{22}$
identically. This reflects a permutation symmetry of the system: when $p=q$,
switching states $1$ and $2$ does not affect the transition matrix. In this
case, the Markov chain smoothly transitions between two regimes: for $s\ll 0$,
the chain spends half of the time in node $1$ and half on $2$; for $s\gg 0$,
the chain spends all the time jumping from state $1$ to state $2$ and back.
When $p\neq q$, instead, for $s<0$ the system smoothly transitions to a
localised state, where the Markov chain is mostly located on $1$ (resp. $2$)
if $p<q$ (resp. $p>q$). Interestingly, we notice that for $p<q$ the maximum of
$\nu^{*}_{22}$ occurs at a finite and negative value of $s$.

In the top-right and bottom-right panels of Fig. 5 we plot each contribution
to the SCGF obtained with our approach alongside their sum. The SCGF is in
perfect agreement with the one obtained using spectral methods. Furthermore,
our approach allows us to understand the magnitude of each physical term. We
split $\lambda_{1}(s)$, as defined in (41), into two terms as follows:

$\displaystyle\lambda_{1,\text{states}}=\sum_{i=1}^{2}\sum_{j=1}^{2}\nu_{ij}\log\left(\sum_{k=1}^{2}\nu_{ik}\right)$
(77)
$\displaystyle\lambda_{1,\text{links}}=-\sum_{i=1}^{2}\sum_{j=1}^{2}\nu_{ij}\log(\nu_{ij})\
,$ (78)

and plot them separately. In the top-right panel we used $p=q=0.5$, while in
the bottom-right panel we used $p=0.5$ and $q=0.9$. In both cases, when $|s|$
is large we notice that $\lambda_{1,\text{states}}$ and
$\lambda_{1,\text{links}}$ balance each other, and their sum is close to zero.
This is because in both cases, the Markov chain spends most of the time in
just a fraction of the available links. For $s\gg 0$, the dominant
contribution in both cases is due to the tilting term $\lambda_{3}(s)$. For
$s\ll 0$, we notice a striking difference: when $p=q$, both
$\lambda_{1,\text{states}}$ and $\lambda_{1,\text{links}}$ tend to a finite
value. This is because the chain still visits both node $1$ and node $2$.
Instead, in the case $p\neq q$, $\lambda_{1,\text{states}}$ and
$\lambda_{1,\text{links}}$ tend to $0$. This is because of the aforementioned
localisation behaviour. In both cases, since the tilting term $\lambda_{3}(s)$
becomes negligible, the SCGF $\lambda$ is well approximated by $\lambda_{2}$,
the dynamical entropy.

##  V Conclusion

In this work we propose a way to study the large deviation regime of
fluctuations of two-point observables of a discrete-time Markov chain.
Adopting graph-combinatorial arguments similar to those in [33, 34, 35], we
show how to calculate the finite-time moment generating function and the
scaled cumulant generating function, objects that have a clear interpretation
in the framework of statistical physics as they correspond, respectively, to
the canonical partition function and Helmoltz free energy. In particular, all
terms of the Helmoltz free energy have a clear physical meaning—see (40) and
following discussion. We establish a direct and explicit link with spectral
methods, as the Lagrange multipliers in (40) can be shown to be the dominant
eigenvalue and right eigenvector of the tilted matrix—see (47). Furthermore,
from the minimisers $\nu^{*}$ we show how to compute in a simple way the
occupation measure on the nodes and the driven process.

We illustrate the benefits of our method in a general two-state model, for
which we can compute analytically both the moment generating function and the
SCGF. We show plots where we highlight the new information accessible with our
method: in particular, we compare the different contributions of paths and
cycles to the moment generating function. For the large deviation regime,
analysing the minimisers $\nu^{*}$ as well as all the terms in our formula for
the SCGF, we find an interesting localisation behaviour of the Markov chain
when the two-state model is not symmetric.

Remarkably, the finite-time expression for the moment generating function
could be used as the starting point for future investigations on the role of
sub-leading terms in the fluctuations of observables, for which to our
knowledge not much is known. A remarkable contribution in this direction is
[39], where authors use matrix products states to characterise fluctuations at
finite time. An interesting avenue for future research would be to try to
apply our methods in the continuous-time setting.

Furthermore, once we fix the state-space connectivity and probability weights,
it would be interesting to understand the interplay between the long-time
limit and the large number of states limit. In the framework of Markov chains
satisfying detailed balance, this approach could, in principle, be adopted to
investigate transient behaviour and metastability in rough energy landscapes,
a problem relevant to many areas in statistical physics [47, 48]. More
generally, Markov chains that satisfy global balance but not detailed balance
are a paradigmatic model for out of equilibrium phenomena. In this context,
understanding finite-time behaviour is challenging—see [49] for applications
to biology. Out of equilibrium steady states are directly accessible in the
large deviations framework [19, 50] and are of interest to many communities.

Finally, we remark that the Helmoltz free energy associated with the pair
empirical occupation measure is a powerful tool to investigate dynamical phase
transitions in fluctuations of one and two-point observables. For instance, in
[21, 24] the authors show evidence of a localisation phase transition in
random walks on random graphs. In an upcoming work, we intend to investigate
toy models where this phenomenon can be analytically characterised using the
approach outlined in this paper.

## Acknowledgments

GC and FC are thankful to Gianmichele Di Matteo for insightful discussions and
to Mayank Shreshtha for having designed Fig. 2. FC is grateful to Hugo
Touchette for pointing to interesting literature in the topic and for the
hospitality in Stellenbosch (South Africa) during the writing stage of the
manuscript. GC is supported by the EPSRC Centre for Doctoral Training in
Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES,
EP/L015854/1).

## References

  * [1] R. Chetrite and H. Touchette, “Nonequilibrium Markov processes conditioned on large deviations,” Annales Henri Poincare, vol. 16, no. 9, pp. 2005–2057, 2015.
  * [2] K. Sekimoto, Stochastic Energetics, vol. 799 of Lecture Notes in Physics.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.
  * [3] C. Maes and K. Netovcný, “Canonical structure of dynamical fluctuations in mesoscopic nonequilibrium steady states,” Europhysics Letters, vol. 82, no. 3, 2008.
  * [4] M. Baiesi, C. Maes, and B. Wynants, “Fluctuations and response of nonequilibrium states,” Physical Review Letters, vol. 103, no. 1, p. 010602, 2009.
  * [5] R. Gutierrez and C. Perez-Espigares, “Generalized optimal paths and weight distributions revealed through the large deviations of random walks on networks,” Physical Review E, vol. 103, no. 2, p. 022319, 2021.
  * [6] A. Dechant, “Minimum entropy production, detailed balance and Wasserstein distance for continuous-time Markov processes,” arXiv:2110.01141, 2021\.
  * [7] B. Derrida, “Non-equilibrium steady states: Fluctuations and large deviations of the density and of the current,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2007, no. 7, p. P07023, 2007.
  * [8] G. Verley, M. Esposito, T. Willaert, and C. Van Den Broeck, “The unlikely Carnot efficiency,” Nature Communications, vol. 5, no. 1, pp. 1–5, 2014\.
  * [9] G. Verley, T. Willaert, C. Van Den Broeck, and M. Esposito, “Universal theory of efficiency fluctuations,” Physical Review E, vol. 90, no. 5, p. 052145, 2014.
  * [10] T. R. Gingrich, G. M. Rotskoff, S. Vaikuntanathan, and P. L. Geissler, “Efficiency and large deviations in time-asymmetric stochastic heat engines,” New Journal of Physics, vol. 16, no. 10, p. 102003, 2014.
  * [11] S. K. Manikandan, L. Dabelow, R. Eichhorn, and S. Krishnamurthy, “Efficiency fluctuations in microscopic machines,” Physical Review Letters, vol. 122, no. 14, p. 140601, 2019.
  * [12] F. Coghi and R. J. Harris, “A large deviation perspective on ratio observables in reset processes: robustness of rate functions,” Journal of Statistical Physics, vol. 179, no. 1, pp. 131–154, 2020.
  * [13] J. L. Lebowitz and H. Spohn, “A Gallavotti–Cohen-type symmetry in the large deviation functional for stochastic dynamics,” Journal of Statistical Physics 1999 95:1, vol. 95, no. 1, pp. 333–365, 1999.
  * [14] J. Mehl, T. Speck, and U. Seifert, “Large deviation function for entropy production in driven one-dimensional systems,” Physical Review E, vol. 78, no. 1, p. 011123, 2008.
  * [15] J.-D. Deuschel and D. W. Stroock, Large deviations.  Academic Press, 1989.
  * [16] F. den Hollander, Large Deviations.  American Mathematical Society, 2000.
  * [17] H. Touchette, “The large deviation approach to statistical mechanics,” Physics Reports, vol. 478, no. 1-3, pp. 1–69, 2009.
  * [18] A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications, vol. 38 of Stochastic Modelling and Applied Probability.  Springer Berlin Heidelberg, 2010.
  * [19] H. Touchette and R. J. Harris, “Large deviation approach to nonequilibrium systems,” in Nonequilibrium Statistical Physics of Small Systems: Fluctuation Relations and Beyond, vol. 6, pp. 335–360, Weinheim: Wiley-VCH, 2013\.
  * [20] R. Chetrite and H. Touchette, “Variational and optimal control representations of conditioned and driven processes,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2015, no. 12, p. P12001, 2015\.
  * [21] C. De Bacco, A. Guggiola, R. Kühn, and P. Paga, “Rare events statistics of random walks on networks: localisation and other dynamical phase transitions,” Journal of Physics A: Mathematical and Theoretical, vol. 49, no. 18, p. 184003, 2016.
  * [22] P. Tsobgni Nyawo and H. Touchette, “Large deviations of the current for driven periodic diffusions,” Physical Review E, vol. 94, no. 3, p. 032101, 2016.
  * [23] S. Whitelam, “Large deviations in the presence of cooperativity and slow dynamics,” Physical Review E, vol. 97, p. 62109, 2018.
  * [24] F. Coghi, J. Morand, and H. Touchette, “Large deviations of random walks on random graphs,” Physical Review E, vol. 99, no. 2, p. 022137, 2019.
  * [25] S. Whitelam and D. Jacobson, “Varied phenomenology of models displaying dynamical large-deviation singularities,” Physical Review E, vol. 103, no. 3, p. 032152, 2021.
  * [26] A. C. Barato and R. Chetrite, “A formal view on level 2.5 large deviations and fluctuation relations,” Journal of Statistical Physics, vol. 160, no. 5, pp. 1154–1172, 2015.
  * [27] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, “Dissipation bounds all steady-state current fluctuations,” Physical Review Letters, vol. 116, no. 12, p. 120601, 2016.
  * [28] F. Coghi, R. Chetrite, and H. Touchette, “Role of current fluctuations in nonreversible samplers,” Physical Review E, vol. 103, no. 6, p. 062142, 2021.
  * [29] J. Hoppenau, D. Nickelsen, and A. Engel, “Level 2 and level 2.5 large deviation functionals for systems with and without detailed balance,” New Journal of Physics, vol. 18, no. 8, p. 083010, 2016.
  * [30] P. Whittle, “Some distribution and moment formulae for the Markov chain,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 17, no. 2, pp. 235–242, 1955.
  * [31] P. Billingsley, “Statistical methods in Markov chains,” The Annales of Mathematical Statistics, vol. 32, no. 1, pp. 12–40, 1961.
  * [32] I. Csiszár, T. M. Cover, and B. S. Choi, “Conditional Limit theorems under Markov conditioning,” IEEE Transactions on Information Theory, vol. 33, no. 6, pp. 788–801, 1987.
  * [33] R. Dawson and I. J. Good, “Exact Markov probabilities from oriented linear graphs,” The Annales of Mathematical Statistics, vol. 28, no. 4, pp. 946–956, 1957.
  * [34] L. A. Goodman, “Exact probabilities and asymptotic relationships for some statistics from $m$-th order Markov chains,” The Annals of Mathematical Statistics, vol. 29, no. 2, pp. 476–490, 1958.
  * [35] M. Polettini, “BEST statistics of Markovian fluxes: a tale of Eulerian tours and Fermionic ghosts,” Journal of Physics A: Mathematical and Theoretical, vol. 48, p. 365005, aug 2015.
  * [36] R. L. Jack and P. Sollich, “Large deviations and ensembles of trajectories in stochastic models,” Progress of Theoretical Physics Supplement, vol. 184, pp. 304–317, mar 2010.
  * [37] R. Chetrite and H. Touchette, “Nonequilibrium microcanonical and canonical ensembles and their equivalence,” Physical Review Letters, vol. 111, no. 12, p. 120601, 2013.
  * [38] R. S. Ellis, Entropy, Large Deviations, and Statistical Mechanics.  Classics in Mathematics, Springer New York, 1985.
  * [39] L. Causer, M. C. Bañuls, and J. P. Garrahan, “Finite time large deviations via matrix product states,” Physical Review Letters, vol. 128, p. 090605, mar 2022.
  * [40] I. Kontoyiannis and S. P. Meyn, “Spectral theory and limit theorems for geometrically ergodic markov processes,” The Annals of Applied Probability, vol. 13, no. 1, pp. 304–362, 2003.
  * [41] I. Kontoyiannis and S. Meyn, “Large deviations asymptotics and the spectral theory of multiplicatively regular markov processes,” Electronic Journal of Probability, vol. 10, pp. 61–123, 2005.
  * [42] M. E. J. Newman, Networks: An Introduction.  Oxford University Press, 2010.
  * [43] van T. Aardenne-Ehrenfest and de N.G. Bruijn, “Circuits and trees in oriented linear graphs,” Simon Stevin, vol. 28, pp. 203–217, 1951.
  * [44] W. T. Tutte and C. A. B. Smith, “On unicursal paths in a network of degree 4,” The American Mathematical Monthly, vol. 48, no. 4, p. 233, 1941.
  * [45] M. Farrell and L. Levine, “Multi-Eulerian tours of directed graphs,” The Electronic Journal of Combinatorics, vol. 22, 2015.
  * [46] M. Rubey, Counting spanning trees.  PhD thesis, Universität Wien, 2000.
  * [47] A. Baronchelli, A. Barrat, and R. Pastor-Satorras, “Glass transition and random walks on complex energy landscapes.,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 80 2 Pt 1, p. 20102, 2009\.
  * [48] F. Arceri, F. P. Landes, L. Berthier, and G. Biroli, “Glasses and aging: A statistical mechanics perspective,” 2020.
  * [49] C. Jia, M. Qian, and D. Jiang, “Overshoot in biological systems modelled by Markov chains: a non-equilibrium dynamic phenomenon.,” IET systems biology, vol. 8, no. 4, pp. 138–145, 2014.
  * [50] R. Cofré, L. Videla, and F. Rosas, “An introduction to the non-equilibrium steady states of maximum entropy spike trains,” Entropy, vol. 21, no. 9, 2019.