Open Markov Chains: Cumulant Dynamics, Fluctuations and Correlations

Abstract

In this work we propose a model for open Markov chains that can be interpreted as a system of non-interacting particles evolving according to the rules of a Markov chain. The number of particles in the system is not constant, because we allow the particles to arrive or leave the state space according to prescribed protocols. We describe this system by looking at the population of particles on every state by establishing the rules of time-evolution of the distribution of particles. We show that it is possible to describe the distribution of particles over the state space through the corresponding moment generating function. This description is given through the dynamics ruling the behavior of such a moment generating function and we prove that the system is able to attain the stationarity under some conditions. We also show that it is possible to describe the dynamics of the two first cumulants of the distribution of particles, which in some way is a simpler technique to obtain useful information of the open Markov chain for practical purposes. Finally we also study the behavior of the time-dependent correlation functions of the number of particles present in the system. We give some simple examples of open chains that either, can be fully described through the moment generating function or partially described through the exact solution of the cumulant dynamics.

1. Introduction

Markov chains are discrete-time models for stochastic evolution, widely used to model systems in physics [1,2], chemistry [2], biology [3,4] as well as social sciences [5,6,7]. Roughly speaking, a Markov chain consists of a sequence of random variables $\{X_t\in \mathcal{S}:t\in {\mathbb{N}}_0\}$ taking values from a (finite or countable) set $\mathcal{S}$, called state space. The jump of the random variable from one state to another in one time step occurs with a prescribed probability, and the probabilities of all the possible jumps are collected in a matrix called Markov matrix which is a stochastic matrix. A natural way to interpret a Markov chain comes from physics; we can think of the random variable $X_t$ as the position at time t of given particle, and this particle moves on the discrete space $\mathcal{S}$. Thus, the stationary probability vector $\pi$ (provided it exists) is interpreted from a point of view of ensembles: if we have a collection of N non-interacting particles moving according to the rules of the Markov chain, then, the stationary distribution of particles on $\mathcal{S}$ is $N\pi$ if N is large enough. This point of view clearly shows that a Markov chain is a closed system, since there is no inflow of particles to $\mathcal{S}$ nor outflow of particles from $\mathcal{S}$.

In this paper we shall be concerned with the case in which we allow the particles to arrive or leave the state space according to a prescribed protocol. To be precise let us consider a state $j\in \mathcal{S}$. On one hand, at every time step we allow a certain number of particles already present in the state j to leave this state to the “outside”. On the other hand we also allow a certain number of particles (from the “outside”) to arrive at the state j. Both, the number of incoming particles and the number of outgoing particles, are modeled as random variables (or sequences of random variables) with a distribution given a priori. Our main goal in this paper is to describe the population of particles on the state space as well as its fluctuations. We are particularly interested in the behavior of both, the space correlations and time correlations for several possible scenarios for the incoming and outgoing protocols.

At this point it is convenient to mention some works related to our model. In the current literature one can found several models for open stochastic processes. Among the earlier results on models of open Markov populations we can found those due to Bartholomew [7] and Gani [8]. Since then, several authors have obtained results in different classes of models for the evolution of open Markov populations. For instance, Guerreiro and Mexia have studied a class of models called stochastic vortices model [9] as an alternative approach to estimate the so-called long run distribution for a bonus malus system [10,11]. Guerreiro and colleges have also studied a class of open of population models subjected to periodical re-classifications [12,13,14]. The case of open Markov chains with Poisson recruitment has been dealt in [15], and remarkable extensions of these models to continuous time has been studied in [16,17,18,19,20,21,22,23]. Open Markov chain schemes fed by a second order stationary and non-stationary processes have also been studied in [24] where the authors consider that the inflow of new population elements is modeled by a time series coming from a second order stationary process, i.e., a stationary process with a deterministic bias. Additionally, open Markov chains with not independent inflow processes has been considered in [25]. All these different models of stochastic open systems have been useful for several important applications. These models have been used to study consumption credit portfolios [26], open automobile portfolios [27,28], hospital planning [19], enrollments and degrees awarded in universities [8], manpower models [20,21,29] as well as traffic flow by means of random networks [30]. Open Markov processes have also been used to understand in general several properties of bonus-malus systems [9,10,11,27]. Particularly Floriani et al. [30] focus their study to the case where the Markov chain has some absorbing states (in which the mass accumulates) and the mass is supplied by either, a constant source or a periodic source. In contrast, in our model the number of incoming particles at every time step is not constant but a sequence of random variables not necessarily independent and identically distributed (i.i.d.). Moreover, instead of modeling the outflow by absorbing states, we define a protocol of outgoing particles, which allows every particle to leave the chain with a state-dependent probability. Another work which is worth mentioning here is the one of Pollard and co-workers [22,31,32]. They consider a class of open Markov process in which the state space is discrete and the time continuous. They assume that the incoming and outgoing fluxes are regulated by means of a set of special states of $\mathcal{S}$ (which they call “boundary”). One of the main differences of our approach with respect to the one proposed by Pollard is that they suppose that the elements of the boundary has a distribution prescribed a priori (which may be even time-dependent). In contrast, our approach considers every state as a source or sink of particles. In this sense our model can be thought of as a Markov chain in contact with a reservoir of particles. Thus, according to a prescribed protocol the particles go from the reservoir to the chain and vice versa, the particles go from the chain to the reservoir with a prescribed protocol. This way of modeling the source of particle is, in some way, similar to the grand canonical ensemble in thermodynamics, in which a system is in contact with a reservoir allowing particles to be interchanged.

2. A Model for Open Markov Chains

The main idea behind our model for an “open” Markov chain is that we allow the particles to arrive and leave the state space $\mathcal{S}$. We have mentioned that the particle can enter the state space according to a protocol which is modeled as a sequence of random variables. Such a sequence of random variables determines the number of particles arriving at certain state every time step. On the other hand, the particle can leave the state space depending on which state they are. The most natural way to model this situation is by defining, for every state, a given probability with which the particle leaves such a state towards the reservoir. This probability must satisfy a compatibility condition, consisting in the fact that a particle in a given state has only two options (i) jump to any other state in $\mathcal{S}$ or (ii) jump to the reservoir. The sum of all these probabilities should be one, in order for the “jump” to be well defined. Notice that due to the compatibility condition, we have that the protocol of outgoing particles is completely determined by means of a non-negative matrix $\mathbf{Q}$ with a spectral radius strictly less than one. This is because the “missing probability” in $\mathbf{Q}$ (necessary for $\mathbf{Q}$ to be a stochastic matrix) is interpreted as “jump probabilities” towards the outside.

Definition 1

(Open Markov chain). Let $\mathcal{S}$ be a finite set, whose cardinality is denoted by $\#\mathcal{S}=S$, and let $\mathbf{Q}:\mathcal{S}\times \mathcal{S}\to [0,1]\subset \mathbb{R}$ be an irreducible and aperiodic matrix with spectral radius strictly less than one. Let $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ be a sequence of random vectors taking values in ${\mathbb{N}}_0^S$. We say that $\left(\mathcal{S},\mathbf{Q},\{{\mathbf{J}}^t:t\in \mathbb{N}\}\right)$ is an open Markov chain with state space $\mathcal{S}$, jump matrix $\mathbf{Q}$ and incoming protocol $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$. Now let $\{{\mathbf{N}}^t:t\in \mathbb{R}\}$ be a sequence of random vectors taking values in ${\mathbb{N}}_0^S$. Such a sequence is defined as follows. Given the initial random vector ${\mathbf{N}}^0$ with a given distribution, we define ${\mathbf{N}}^t$ recursively as

$${\mathbf{N}}^{t+1}:={\mathbf{J}}^t+{\mathbf{R}}^t,$$

(1)

where ${\mathbf{R}}^t$ is a random vector taking values in ${\mathbb{N}}_0^S$, whose components are given by,

$$R_j^t:=\sum _{i=1}^S{B}_{i,j}^t.$$

(2)

The random variables $B_{i,j}^t$ are defined such that the $(S+1)$-dimensional vector ${\mathbf{A}}_i^t$, with components

$${\left({\mathbf{A}}_i^t\right)}_j=\left\{\begin{array}{cc}\hfill B_{i,j}^t& \hfill 1\le j\le S\\ \hfill 1-{\sum }_{j=1}^S{B}_{i,j}^t& \hfill j=S+1,\end{array}\right.$$

has multinomial distribution, i.e., ${\mathbf{A}}_i^t\sim Multinomial({\mathbf{z}}_i,N_i^t)$. Additionally we assume that ${\mathbf{A}}_i^t$ and ${\mathbf{A}}_i^s$ are independent if $t\ne s$ and that ${\mathbf{A}}_i^t$ and ${\mathbf{A}}_j^t$ are independent if $i\ne j$. Here $\{{\mathbf{z}}_i:1\le i\le S\}$ is a set of probability vectors defined as

$${\left({\mathbf{z}}_i\right)}_j=\left\{\begin{array}{cc}\hfill q_{i,j}& \hfill 1\le j\le S\\ \hfill e_i& \hfill j=S+1,\end{array}\right.$$

where $q_{i,j}$ is the $(i,j)$th component of $\mathbf{Q}$ and $e_i$ is defined as,

$$e_i:=1-\sum _{j=1}^S{q}_{i,j}.$$

We say that ${\mathbf{N}}^t$ is the distribution over the state space at time t with initial condition ${\mathbf{N}}^0$.

Notice that Equation (1) establishes the evolution of the number of particles. This equation states that the number of particles at time $t+1$ is the number of particles having arrived as the state space (represented by ${\mathbf{J}}^t$) plus the number of particles having remained in the state space (which is represented by ${\mathbf{R}}^t$). Observe that the random vector ${\mathbf{R}}^t$ can be seen as a sum of independent random vectors ${\mathbf{B}}_i^t={\left(B_{i,j}\right)}_{j=1}^S$ representing the quantity of particles departing from state i towards other states. Notice also that the random vector ${\mathbf{A}}_i^t$ is the “enlarged” version of the vector ${\mathbf{B}}_i^t$, because the vector ${\mathbf{A}}_i^t$ represents the quantity of particles departing from state i to other states and to the outside. It is clear that, due to conservation of particles during the process of redistribution, the vector ${\mathbf{A}}_i^t$ should have multinomial distribution (its components actually sum up $N_i^t$ as we can see from the definition above). Hence we have that ${\mathbf{N}}^{t+1}$ depends on ${\mathbf{N}}^t$ through the random variable ${\mathbf{R}}^t$ which is the responsible of the redistribution of particles among the internal states.

From the above definition we can appreciate that the quantity of foremost interest is ${\mathbf{N}}^t$, the distribution of particles on the state space at time t. Our goal is then to provide a way for determining the probability function of ${\mathbf{N}}^t$ and to determine whether or not the process $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$ is able to reach a stationary distribution. We will see that $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$ is actually a Markov process and the main goal is to determine its properties. Before establishing a result on this direction let us give some examples of open Markov chains.

Example 1.

Let us consider the most simple case in which the system has only one state. In this case the matrix $\mathbf{Q}$ consists of a single number q (a $1\times 1$ matrix), which should be non-negative and strictly less than one, i.e., $0\le q<1$. The most simple case for the incoming protocol consists of a sequence of constant random variables all these taking the same value, which we call $J_0\in {\mathbb{N}}_0$. This means that the number of particles arriving at every time-step is $J_0$. Every particle arriving at the unique available state has only two options jump to the outside or remain in its state. The probability of remaining in the state is q and the probability of jumping to the outside is $1-q$. This simple example for open Markov chain can be illustrated as a graph with only one vertex and three edges as shown in Figure 1. Notice that there are two special edges, one establishing the incoming protocol (labeled by $J_0$) and one defining the outgoing protocol (labeled by $1-q$).

Example 2.

A less trivial example is provided by giving a matrix larger than $1\times 1$. Let us consider for example the matrix given by

$$\mathbf{Q}=\left[\begin{array}{ccc}0& 1/2& 1/4\\ 1/4& 1/4& 0\\ 1/4& 1/2& 1/4\end{array}\right]$$

(3)

Notice that the above matrix is not stochastic, because some rows does not add to one, but less than one. The latter means that not all the states allow the particles to leave to the outside. As we can see the first row adds to $3/4$, which means that every particle on the state 1 has a probability $1/4$ to go out to the reservoir. The second row adds to $1/2$, this means that a given particle on the state 2 has a probability $1/2$ to go out to the outside. Finally, the third row adds to one, meaning that a particle on the state 3 can only jump to the other states 1 or 2 or remains in its current state, but it cannot leave the state space to go to the outside.

Now assume that the incoming protocol $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ is a set of i.i.d. random vectors, i.e., the protocol is time-independent. Then the number of incoming particles at every time-step can be considered as independent realizations of a single random vector, which we denote by $\mathbf{J}$. To be more precise, we can chose in particular this random vector as ${\mathbf{J}}^{\mathbf{t}}=(J_1^t,0,J_3^t)$, with $J_1^t$ and $J_3^t$ two independent random variables with Bernoulli distribution with parameters $p_1$ and $p_3$ respectively. If, for instance, the parameters of every Bernoulli distribution were given by $p_1=0.1$ and $p_3=0.6$, then this would mean that the average number of particles arriving at the vertex 1 is lower than the average number of particles arriving at the vertex 3. Figure 2 shows the graphical representation of this open Markov chain. Then, the number of particles ${\mathbf{N}}^t=(N_1^t,N_3^t,N_3^t)$ evolves according to the rule

$${\mathbf{N}}^{t+1}={\mathbf{J}}^t+{\mathbf{R}}^t,$$

(4)

where ${\mathbf{R}}^t$ is the redistribution random vector given in Equation (2). The above vectorial equation can be written coordinate-wise as follows,

$$\begin{array}{ccc}\hfill N_1^{t+1}& =& J_1^t+R_1^t,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill N_2^{t+1}& =& R_2^t,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill N_3^{t+1}& =& J_3^t+R_2^t.\hfill \end{array}$$

(7)

We should keep in mind that the random vectors $R_i^t$, for $1\le i\le 3$, are conditioned to the values of ${\mathbf{N}}^t$, i.e., the distribution depends on the specific realization of ${\mathbf{N}}^t$, which is a manifestation of the Markovian hypothesis.

At this point it is interesting to compare the characteristics of the proposed model for open Markov chains with the continuous-time Markov process introduced in Ref. [22]. At first sight we would expect that the model described here might match with a discrete-time version of the formalism proposed in [22]. However this is not the case in general, which can be seen from the very definitions of both models. The main difference between both approaches is that the model in Ref. [22] considers a state space V with special states members of a boundary B. The role of the boundary is to provide the system with inflow and outflow of particles. In the remaining of the states, elements of V\B, occurs the dynamics of the system. In the approach proposed here, the state space $\mathcal{S}$ can, in principle, be identified with the latter, i.e., $\mathcal{S}=V$\B. However, we should recall that we admit the that every state of $\mathcal{S}$ is open, i.e., particles from the outside might arrive at any state of $\mathcal{S}$ and particles within the system might leave $\mathcal{S}$ from any state. This means that we should also interpret the boundaryB as the whole state space $\mathcal{S}$, but this is impossible since we identified $\mathcal{S}$ with V\B. This shows that both approaches does not math each other in general, although it would be interesting if they might coincide in particular cases.

3. The Evolution of the Particle Distribution

3.1. Evolution of the Moment Generating Function

In this section we will establish the evolution equation for the process ${\mathbf{N}}^t$, which, as we have anticipated, is a Markov process. This fact can actually be appreciated in Equation (1), where we have indicated that the number of particles at time $t+1$ is uniquely determined by the number of incoming particles and the redistribution of particles that were present at time t. Clearly the last condition states the Markov property for the process $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$. In order to obtain the stochastic matrix governing the behavior of ${\mathbf{N}}^t$, let us define $p_t\left(\mathbf{n}\right)$ as the probability vector associated to ${\mathbf{N}}^t$, i.e.,

$$p_t\left(\mathbf{n}\right):=\mathbb{P}({\mathbf{N}}^t=\mathbf{n}).$$

(8)

We will refer to $p_t\left(\mathbf{n}\right)$ as the distribution of particles over the state space or, to simplify, distribution over the state space Let us consider the probability vector at time $t+1$. Notice that, due the fact that ${\mathbf{N}}^{t+1}$ depends only on ${\mathbf{N}}^t$, which is a consequence of the relation (1), as we have mentioned above, it is clear that

$$\begin{array}{ccc}\hfill p_{t+1}\left(\mathbf{n}\right)& =& \mathbb{P}({\mathbf{N}}^{t+1}=\mathbf{n})\hfill \\ & =& \sum _{\mathbf{k}\in {\mathbb{N}}_0^S}\mathbb{P}({\mathbf{N}}^t=\mathbf{k})\mathbb{P}({\mathbf{N}}^{t+1}=\mathbf{n}|{\mathbf{N}}^t=\mathbf{k}).\hfill \end{array}$$

(9)

Observe that the last expression is a consequence of the Markov property of the process which was assumed in Equation (1). Notice that Equation (9) establishes the evolution equation for $p_t$, which can be written as

$$\begin{array}{ccc}\hfill p_{t+1}\left(\mathbf{n}\right)& =& \sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right)\mathbb{P}({\mathbf{N}}^{t+1}=\mathbf{n}|{\mathbf{N}}^t=\mathbf{k}).\hfill \end{array}$$

(10)

Now let us call $K(\mathbf{k},\mathbf{n})$ the conditional probability appearing in the above equation, i.e.,

$$K(\mathbf{k},\mathbf{n}):=\mathbb{P}({\mathbf{N}}^{t+1}=\mathbf{n}|{\mathbf{N}}^t=\mathbf{k}),$$

(11)

and observe that this quantity can be rewritten as follows

$$K(\mathbf{k},\mathbf{n})=\mathbb{P}({\mathbf{J}}^t+{\mathbf{R}}^t=\mathbf{n}).$$

(12)

The dependence on $\mathbf{k}$ in the above expression is implicit in the random vector ${\mathbf{R}}^t$, since the redistribution of particles depends on the number of particles on every state at time t, which is indeed given by $\mathbf{k}$. Thus, the function $K:{\mathbb{N}}_0^S\times {\mathbb{N}}_0^S\to [0,1]\subset \mathbb{R}$ can be thought as the stochastic matrix corresponding to the process $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$, since the evolution equation for the probability vector $p_t\left(\mathbf{n}\right)$ is given in terms of K as follows,

$$p_{t+1}\left(\mathbf{n}\right)=\sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right)K(\mathbf{k},\mathbf{n}).$$

(13)

In order to solve Equation (13) for $p_t\left(\mathbf{n}\right)$ it is necessary to make some assumptions on the nature of the random vectors ${\mathbf{J}}^t$ and ${\mathbf{R}}^t$ for all t. First of all, it is natural to assume that ${\mathbf{J}}^t$ and ${\mathbf{R}}^t$ are independent. This assumption actually means that the number of particles incoming to the state space does not have any influence on the redistribution of particles already present in the chain. This implies that the joint probability for the random vectors ${\mathbf{J}}^t$ and ${\mathbf{R}}^t$ can be factorized as the product of its corresponding probability vectors, i.e.,

$$\mathbb{P}({\mathbf{J}}^t=\mathbf{j};{\mathbf{R}}^t=\mathbf{r})=\mathbb{P}({\mathbf{J}}^t=\mathbf{j})\mathbb{P}({\mathbf{R}}^t=\mathbf{r}).$$

The above equality allows us to express the kernel $K(\mathbf{k},\mathbf{n})$ as,

$$\begin{array}{ccc}\hfill K(\mathbf{k},\mathbf{n})& :=& \mathbb{P}({\mathbf{J}}^t+{\mathbf{R}}^t=\mathbf{n})\hfill \\ & =& \sum _{\mathbf{j}+\mathbf{r}=\mathbf{n}}\mathbb{P}({\mathbf{J}}^t=\mathbf{j})\mathbb{P}({\mathbf{R}}^t=\mathbf{r}).\hfill \end{array}$$

(14)

Next, we will solve for $p_t$ by using the well-known technique of the moment generating function (m.g.f.). To this end let us introduce some notation. Let ${\mathcal{G}}_t:{\mathbb{R}}^S\times {\mathbb{R}}^S\to \mathbb{R}$ be the m.g.f. of ${\mathbf{N}}^t$, which is defined as

$${\mathcal{G}}_t\left(\mathit{\alpha }\right):=\mathbb{E}\left[e^{{\mathbf{N}}^t}{\mathit{\alpha }}^T\right]=\sum _{\mathbf{n}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{n}\right)e^{\mathbf{n}{\mathit{\alpha }}^T}.$$

(15)

At this point it is important to describe our convention for vectors in ${\mathbb{R}}^S$. First of all we should emphasize that we interpret the vectors $\mathbf{n}$, $\mathit{\alpha }$, ${\mathbf{N}}^t$, etc., as row vectors (i.e., matrices of size $1\times S$). Thus, the superscript T means, as usual, matrix transposition implying that the vector ${\mathit{\alpha }}^T$ is a column vector (a matrix of size $S\times 1$). Thus, within this convention, the product $\mathbf{n}{\mathit{\alpha }}^T$ should be understood in the sense of the usual matrix product, which in this case results in a single number.

Analogously, we also define the moment generating functions for ${\mathbf{J}}^t$ and ${\mathbf{R}}^t$ as follows,

$$\begin{array}{ccc}\hfill {\mathcal{F}}_t\left(\mathit{\alpha }\right)& :=& \mathbb{E}\left[e^{{\mathbf{J}}^t{\mathit{\alpha }}^T}\right]=\sum _{\mathbf{j}\in {\mathbb{N}}_0^S}\mathbb{P}({\mathbf{J}}^t=\mathbf{j})e^{\mathbf{j}{\mathit{\alpha }}^T},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \mathcal{H}\left(\mathit{\alpha }\right)& :=& \mathbb{E}\left[e^{{\mathbf{R}}^t{\mathit{\alpha }}^T}\right]=\sum _{\mathbf{r}\in {\mathbb{N}}_0^S}\mathbb{P}({\mathbf{R}}^t=\mathbf{r})e^{\mathbf{r}{\mathit{\alpha }}^T}.\hfill \end{array}$$

(17)

Notice that we omitted the superscript t in the m.g.f. for ${\mathbf{R}}^t$, since two random vectors, say for example ${\mathbf{R}}^t$ and ${\mathbf{R}}^s$, are independent and have the same distribution, and consequently, share the same moment generating function.

Next, our objective is to provide a recurrence relation for ${\mathcal{G}}_t$ using the evolution Equation (13). Thus, let us consider the m.g.f. for ${\mathbf{N}}^{t+1}$, and note that

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)& =& \sum _{\mathbf{n}\in {\mathbb{N}}_0^S}{p}_{t+1}\left(\mathbf{n}\right)e^{\mathbf{n}{\mathit{\alpha }}^T}\hfill \\ & =& \sum _{\mathbf{n}\in {\mathbb{N}}_0^S}\sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right)K(\mathbf{k},\mathbf{n})e^{\mathbf{n}{\mathit{\alpha }}^T}\hfill \\ & =& \sum _{\mathbf{n}\in {\mathbb{N}}_0^S}\sum _{\mathbf{k}\in {\mathbb{N}}_0^S}\sum _{\mathbf{j}+\mathbf{r}=\mathbf{n}}{p}_t\left(\mathbf{k}\right)\mathbb{P}({\mathbf{J}}^t=\mathbf{j})\mathbb{P}({\mathbf{R}}^t=\mathbf{r})e^{\mathbf{n}{\mathit{\alpha }}^T},\hfill \end{array}$$

(18)

where we have used the form of K given in Equation (14). We can appreciate that the summation over $\mathbf{n}$ together with the summation over the restriction $\mathbf{j}+\mathbf{r}=\mathbf{n}$ results in a double sum over the “indices” $\mathbf{j}$ and $\mathbf{r}$ without restrictions, i.e., we obtain two sums over independent indices. This observation allows us to write

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)& =& \sum _{\mathbf{k}\in {\mathbb{N}}_0^S}\sum _{\mathbf{j}\in {\mathbb{N}}_0^S}\sum _{\mathbf{r}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right)\mathbb{P}({\mathbf{J}}^t=\mathbf{j})\mathbb{P}({\mathbf{R}}^t=\mathbf{r})e^{(\mathbf{j}+\mathbf{r}){\mathit{\alpha }}^T}.\hfill \end{array}$$

(19)

In the above Equation (19) we can observe that the summation over $\mathbf{j}$ and $\mathbf{r}$ results in the m.g.f. for ${\mathbf{J}}^t$ and ${\mathbf{R}}^t$ respectively. This implies that

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)& =& \sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right){\mathcal{F}}_t\left(\mathit{\alpha }\right)\mathcal{H}\left(\mathit{\alpha }\right).\hfill \end{array}$$

(20)

In Appendix A we show that $\mathcal{H}\left(\mathit{\alpha }\right)$ can be written as

$$\mathcal{H}\left(\mathit{\alpha }\right)=e^{\mathbf{k}{\mathbf{H}}^T\left(\mathit{\alpha }\right)},$$

(21)

where the function $\mathbf{H}:{\mathbb{R}}^S\to {\mathbb{R}}^S$ is defined as follows. If $H_i\left(\mathit{\alpha }\right)={\left(\mathbf{H}\left(\mathit{\alpha }\right)\right)}_i$ denotes the ith component of $\mathbf{H}$ we have that

$$H_i\left(\mathit{\alpha }\right):=log\left(e_i+\sum _{j=1}^S{q}_{i,j}{e}^{{\mathit{\alpha }}_j}\right).$$

(22)

Observe that Equation (21) shows explicitly the dependence on $\mathbf{k}$ of the m.g.f. for ${\mathbf{R}}^t$. Then, if we substitute the relation (21) into (20), we obtain,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)& =& \sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right){\mathcal{F}}_t\left(\mathit{\alpha }\right)e^{\mathbf{k}{\mathbf{H}}^T\left(\mathit{\alpha }\right)}.\hfill \end{array}$$

(23)

We can easily see that the summation over $\mathbf{k}$ results in the m.g.f. for ${\mathbf{N}}^t$. Thus,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)& =& {\mathcal{F}}_t\left(\mathit{\alpha }\right){\mathcal{G}}_t\left(\mathbf{H}\left(\mathit{\alpha }\right)\right).\hfill \end{array}$$

(24)

Equation (24) is a recurrence relation governing the time-dependence of the m.g.f. for ${\mathbf{N}}^t$. This equation can be formally solved to obtain

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& {\mathcal{G}}_0\left({\mathbf{H}}^{\left(t\right)}\left(\mathit{\alpha }\right)\right)\prod _{r=0}^{t-1}{\mathcal{F}}_{t-r}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right),\hfill \end{array}$$

(25)

where ${\mathcal{G}}_0$ stand for the m.g.f for ${\mathbf{N}}^0$ (the initial distribution on the state space). We should emphasize that the superscript notation ${\mathbf{H}}^{\left(r\right)}$ stands for the rth iteration of the function $\mathbf{H}$, i.e., ${\mathbf{H}}^{\left(r\right)}:=\mathbf{H}\circ \mathbf{H}\circ \dots \mathbf{H}$, r times.

Now, let us assume that the process $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ is a sequence of identically distributed random vectors (not necessarily independent). In this case, the m.g.f. for ${\mathbf{J}}^t$ is the same for all t, consequently the formal solution for ${\mathcal{G}}_t$ can be expressed as,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& {\mathcal{G}}_0\left({\mathbf{H}}^{\left(t\right)}\left(\mathit{\alpha }\right)\right)\prod _{r=0}^{t-1}\mathcal{F}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right).\hfill \end{array}$$

(26)

This result, together with the fact that ${\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\to \mathbf{0}$ as $r\to \infty$ in an open neighborhood around $\mathit{\alpha }=\mathbf{0}$ (see Appendix A for a proof), implies that, for the case in which the random vectors ${\mathbf{J}}^t$ are identically distributed for all t, the process $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$ admits an stationary solution. This is because the m.g.f. ${\mathcal{G}}_t$ attains a limit when $t\to \infty$. Such a limit can be written as,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{\infty }\mathcal{F}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right).\hfill \end{array}$$

(27)

whenever the infinite product exist. If this is the case, the m.g.f. ${\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)$ is additionally a solution for the evolution Equation (24). This means that ${\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)$ corresponds to a m.g.f. of a distribution $p_{\mathrm{stat}}\left(\mathbf{n}\right)$ over the state space, which is invariant under the dynamics (13).

Example 3.

Let us consider the open chain consisting of only one vertex given in Example 1. In this case we will assume that the incoming protocol $\{J^t:t\in \mathbb{N}\}$ consists of a sequence of i.i.d. random variables having a Bernoulli distribution with parameter p. Since all $J^t$ have identical distribution, then we have only one m.g.f. $\mathcal{F}$ characterizing them. This function is given by,

$$\mathcal{F}\left(\mathit{\alpha }\right)=1-p+p{e}^{\mathit{\alpha }}.$$

On the other hand, due to the fact that $\mathbf{Q}$ is a $1\times 1$ matrix, we have that $\mathbf{H}$ is a real-valued function depending on one variable, α, given by,

$$H\left(\mathit{\alpha }\right)=log(1-q+q{e}^{\mathit{\alpha }}).$$

where we denoted by q the unique element of $\mathbf{Q}$. In this case, the rth iteration of H can be exactly calculated. A straightforward calculation shows that,

$$H^{\left(r\right)}\left(\mathit{\alpha }\right)=log\left(1-q^r+q^r{e}^{\mathit{\alpha }}\right).$$

Notice that $H^{\left(r\right)}\left(\mathit{\alpha }\right)\to 0$ as $r\to \infty$, as we have anticipated above. Hence, in this case, the stationary m.g.f. is given by

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{\infty }\mathcal{F}\left(H^{\left(r\right)}\left(\mathit{\alpha }\right)\right)=\prod _{r=0}^{\infty }\left(1-p+p\left(1-q^r+q^r{e}^{\mathit{\alpha }}\right)\right)\hfill \\ & =& \prod _{r=0}^{\infty }\left(1-p{q}^r+p{q}^r{e}^{\mathit{\alpha }}\right).\hfill \end{array}$$

(28)

The above result shows that the stationary distribution can be interpreted as the convolution of an infinite sequence of Bernoulli distributions with parameters $p{q}^r$. This particularly means that the random variable $N^t$, when it has reached the stationarity can be written as a sum of Bernoulli random variables $X_r$ (with parameter $p{q}^r$),

$$N^t=\sum _{r=0}^{\infty }{X}_r.$$

The above result allows us to compute, for example, the expected value $\mathbb{E}\left[N^t\right]$ and the variance,

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N^t\right]& =& \sum _{r=0}^{\infty }\mathbb{E}\left[X_r\right]=\sum _{r=0}^{\infty }p{q}^r=\frac{p}{1-q},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \mathrm{Var}\left(N^t\right)& =& \sum _{r=0}^{\infty }\mathrm{Var}\left(X_r\right)=\sum _{r=0}^{\infty }p{q}^r(1-p{q}^r)=\frac{p}{1-q}+\frac{p^2}{1-q^2}.\hfill \end{array}$$

(30)

We should observe that the fact that the expected value $\mathbb{E}\left[N^t\right]=p/(1-q)$ is finite implies that an equilibration is attained between the number of incoming particles and the number of outgoing particles. It is clear that the mean number of arriving particles per unit time is p, while, the number of leaving particles per unit time is $(1-q)\mathbb{E}\left[N^t\right]$. Once $N^t$ has reached the stationarity, we have an equality between these quantities, giving the result stated above. Although we have obtained the mean number of particles by using the argument of equilibration, the same line of reasoning cannot be applied to the variance. We have thus provided a way to compute the fluctuations in the number of particles that are present in the vertex.

Example 4.

Next, we will deal with another example that is closely related to the problem of the random walk on a ring, which well-known system in physics. A random walk on a ring consist of a particle moving randomly on ${\mathcal{S}}_L:=\mathbb{Z}/L\mathbb{Z}$ for some $L\in \mathbb{N}$. To open this system we allow the particles leave the state space ${\mathcal{S}}_L$ from a given site $i\in {\mathcal{S}}_L$. We will also allow the particles to arrive at the state space from a certain site $j\in {\mathcal{S}}_L$. For simplicity we will consider the case $i=j=0\in {\mathcal{S}}_L$ for arbitrary $L>1$. A schematic representation of the open chain with this state space is depicted in Figure 3. We will assume that once a particle is on a given state $i\in {\mathcal{S}}_L$, it can jump to $i+1$ with probability q or it can jump to $i-1$ with probability $1-q$ for all $i\ne 0$. If the particle is at the state $0\in {\mathcal{S}}_L$ it can jump to the state 1 with probability a and it can jump to the state $L-1$ with probability b. Clearly we have to assume that $0<q<1$ and $0<a+b<1$. Since we are assuming that the particles can only enter to the state space through the state 0, we can see that take the protocol of incoming particles can be written as

$$\begin{array}{ccc}\hfill \mathbf{J}& =& (J_0,0,0,\dots ,0),\hfill \end{array}$$

(31)

where $J_1$ is taken as a Bernoulli random variable with parameter $p\in [0,1]$.

Thus the system that we propose can be seen as a system of random walkers that arrive at the ring and leave the ring from the state 0. Now our goal is to obtain the time-dependent generating function of the distribution of particle over the sate space. As he have shown, the latter is given by

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& {\mathcal{G}}_0\left({\mathbf{H}}^{\left(t\right)}\left(\mathit{\alpha }\right)\right)\prod _{r=0}^{t-1}{\mathcal{F}}_{t-r}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right),\hfill \end{array}$$

(32)

Assuming that at $t=0$ there are no particles on the ring, it is clear that ${\mathcal{G}}_0\left(\mathit{\alpha }\right)=1$. On the other hand, since we assumed that the protocol of incoming particles is $\mathbf{J}=(J_1,0,0,\dots ,0)$, it is clear that the corresponding m.g.f $\mathcal{F}$ is given by

$$\begin{array}{c}\hfill \mathcal{F}\left(\mathit{\alpha }\right)=(1-p)+p{e}^{{\mathit{\alpha }}_0},\end{array}$$

(33)

where ${\mathit{\alpha }}_0$ is the 0th component of α. Thus we have that the m.g.f. of the distribution of particles over the state space can then be written as

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{t-1}\left[(1-p)+p{e}^{{\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right)}_0}\right]=\prod _{r=0}^{t-1}\left[(1-p)+p{e}^{H_0\left({\mathbf{H}}^{(r-1)}\left(\mathit{\alpha }\right)\right)}\right]\hfill \end{array}$$

(34)

Using the definition of $H_i\left(\mathit{\alpha }\right)$ given in Equation (22),

$$H_i\left(\mathit{\alpha }\right):=log\left(e_i+\sum _{j=1}^S{q}_{i,j}{e}^{{\mathit{\alpha }}_j}\right).$$

(35)

it is easy to verify that $e^{H_i\left({\mathbf{H}}^{(r-1)}\left(\mathit{\alpha }\right)\right)}$ can be written as,

$$e^{H_i\left({\mathbf{H}}^{(r-1)}\left(\mathit{\alpha }\right)\right)}=\sum _{k=0}^{r-1}{\left({\mathbf{Q}}^k{\mathbf{e}}^T\right)}_i+\sum _{j=0}^{L-1}{\left({\mathbf{Q}}^r\right)}_{i,j}{e}^{{\mathit{\alpha }}_j},$$

(36)

for all $0\le j\le L-1$. For the sake of simplicity let us adopt the following short-hand notation

$$\begin{array}{c}\hfill {\tilde{p}}_j\left(r\right):={\left({\mathbf{Q}}^r\right)}_{0,j}.\end{array}$$

(37)

for $0\le i\le L-1$. Then notice that due to the fact that $\mathbf{H}\left(\mathbf{0}\right)=\mathbf{0}$ we have the following identity

$$\sum _{k=0}^{r-1}{\left({\mathbf{Q}}^k{\mathbf{e}}^T\right)}_i+\sum _{j=0}^{L-1}{\left({\mathbf{Q}}^r\right)}_{i,j}=1.$$

The latter allows us to rewrite the first term of the right hand side of (36) (with $i=0$) as,

$$\sum _{k=0}^{r-1}{\left({\mathbf{Q}}^k{\mathbf{e}}^T\right)}_0=1-\sum _{j=0}^{L-1}{\tilde{p}}_j\left(r\right).$$

(38)

Then, from Equations (36)–(38), we arrive at the following identity,

$$e^{H_0\left({\mathbf{H}}^{(r-1)}\left(\mathit{\alpha }\right)\right)}=1-\sum _{j=0}^{L-1}{\tilde{p}}_j\left(r\right)+\sum _{j=0}^{L-1}{\tilde{p}}_j\left(r\right)e^{{\mathit{\alpha }}_j}.$$

(39)

Therefore it is clear that the time-dependent m.g.f. of the distribution of particles over the state space can be written as

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{t-1}\left[(1-p)+p\left(1-\sum _{j=0}^{L-1}{\tilde{p}}_j\left(r\right)+\sum _{j=0}^{L-1}{\tilde{p}}_j\left(r\right)e^{{\mathit{\alpha }}_j}\right)\right]\hfill \\ & =& \prod _{r=0}^{t-1}\left(1-\sum _{j=0}^{L-1}{p}_j\left(r\right)+\sum _{j=0}^{L-1}{p}_j\left(r\right)e^{{\mathit{\alpha }}_j}\right),\hfill \end{array}$$

(40)

where we denoted by $p_j\left(k\right)$ the product $p{\tilde{p}}_j\left(k\right)$, i.e.,

$$p_j\left(k\right):=p{\left({\mathbf{Q}}^r\right)}_{0,j}.$$

(41)

At this point it is important to stress that the quantity $p_j\left(r\right)$ is such that $0\le p_j\left(r\right)\le 1$ for all $0\le j\le L-1$ and all $r\in {\mathbb{N}}_0$. Next, we can observe from Equation (41) that

$$0\le \sum _{j=0}^{L-1}{p}_j\left(r\right)\le 1.$$

These properties of allows to interpret the function appearing in the product in Equation (40),

$${\mathcal{K}}^{\left(r\right)}\left(\mathit{\alpha }\right):=1-\sum _{j=0}^{L-1}{p}_j\left(r\right)+\sum _{j=0}^{L-1}{p}_j\left(r\right)e^{{\mathit{\alpha }}_j}$$

(42)

as the a m.g.f. of certain random vector that will be denoted by ${\mathbf{X}}_r$. This random vector has the following characteristics. First of all it is clear that ${\mathbf{X}}_r\in {\mathbb{R}}^L$ and that it can take at most $L+1$ possible values. Specifically,

$${\mathbf{X}}_r={\widehat{\mathbf{e}}}_j,with probability{p}_j\left(r\right),and{\mathbf{X}}_r=\mathbf{0},with probability1-\sum _{j=0}^{L-1}{p}_j\left(r\right).$$

(43)

where ${\widehat{\mathbf{e}}}_j:=(0,0,\dots ,0,1,0,\dots ,0,0)$ is the unitary vector along the jth direction.

This interpretation of the function ${\mathcal{K}}^{\left(r\right)}\left(\mathit{\alpha }\right)$ allows, in turn, interpret the m.g.f. of the distribution of particles over the state space as a convolution of the collection of independent random vectors ${\mathbf{X}}_r$,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{t-1}{\mathcal{K}}^{\left(r\right)}\left(\mathit{\alpha }\right).\hfill \end{array}$$

(44)

The latter means that the time-dependent distribution of particles over the state space ${\mathbf{N}}^t$ can formally be written as a sum of independent random vectors

$$\begin{array}{c}\hfill {\mathbf{N}}^t=\sum _{r=0}^{t-1}{\mathbf{X}}_r.\end{array}$$

(45)

This last expression is useful for several purposes. For instance, the first cumulant can be straightforwardly computed from the last expression,

$$\mathbb{E}\left[{\mathbf{N}}^t\right]=\sum _{r=0}^{t-1}\mathbb{E}\left[{\mathbf{X}}_r\right]=\sum _{r=0}^{t-1}\sum _{j=0}^{L-1}{p}_j\left(r\right),$$

(46)

or equivalently, in terms of the stochastic matrix and the parameter of the protocol of incoming particles, we have that

$$\mathbb{E}\left[{\mathbf{N}}^t\right]=\sum _{r=0}^{t-1}\sum _{j=0}^{L-1}p{\left({\mathbf{Q}}^r\right)}_{0,j}.$$

(47)

This shows that the properties of the time-dependent evolution of the distribution of particles can be obtained by means of the technique of the moment generating function for this class of open systems.

3.2. Cumulant Dynamics

Up to now we have obtained two main results, a recurrence relation of the time-dependent m.g.f. of ${\mathbf{N}}^t$ and a formal expression for the m.g.f. of ${\mathbf{N}}^t$ when the system has attained the stationarity. Now we will focus in two quantities which are of special interest, namely, the two first cumulants for ${\mathbf{N}}^t$, and how do these quantities evolve in time. Let us start by obtaining the first cumulant of ${\mathbf{N}}^t$. Notice that the first cumulant (which coincides with the first moment) can be obtained by taking the first derivative of the m.g.f. ${\mathcal{G}}_t\left(\mathit{\alpha }\right)$. Let us denote by ${\mathit{\mu }}_t$ the first moment of ${\mathbf{N}}^t$, i.e.,

$${\mathit{\mu }}_t:=\mathbb{E}\left[{\mathbf{N}}^t\right],$$

(48)

which we will refer to as the mean distribution of particles over the state space at time t, or simply, the mean distribution at time t. Now let us notice that,

$${\left({\mathit{\mu }}_t\right)}_i=\frac{\partial {\mathcal{G}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i}{|}_{\mathit{\alpha }=\mathbf{0}},$$

(49)

where ${\left({\mathit{\mu }}_t\right)}_i$ is the ith component of ${\mathit{\mu }}_t$ and ${\mathit{\alpha }}_i$ is the ith component of $\mathit{\alpha }$. Observe that the ith component of the mean distribution at time $t+1$ can be obtained from Equation (24), giving

$$\begin{array}{ccc}\hfill {\left({\mathit{\mu }}_{t+1}\right)}_i& =& \frac{\partial {\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i}{|}_{\mathit{\alpha }=\mathbf{0}}\hfill \\ & =& \frac{\partial {\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i}{\mathcal{G}}_t\left(\mathbf{H}\left(\mathit{\alpha }\right)\right){|}_{\mathit{\alpha }=\mathbf{0}}+{\mathcal{F}}_t\left(\mathit{\alpha }\right)\sum _{k=1}^S\frac{\partial {\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_k}\frac{\partial H_k}{\partial {\mathit{\alpha }}_i}{|}_{\mathit{\alpha }=\mathbf{0}}.\hfill \end{array}$$

(50)

Notice that $\mathbf{H}(\mathit{\alpha }=\mathbf{0})=\mathbf{0}$ and any moment generating function evaluated at $\mathbf{0}$ is one. Hence we have,

$$\begin{array}{ccc}\hfill {\left({\mathit{\mu }}_{t+1}\right)}_i& =& {\left({\mathit{ϵ}}_t\right)}_i+\sum _{k=1}^S{\left({\mathit{\mu }}_t\right)}_k{q}_{k,i},\hfill \end{array}$$

(51)

where ${\mathit{ϵ}}_t$ stands for the expected value of ${\mathbf{J}}^t$, which can be obtained by means of the first derivative of ${\mathcal{F}}_t$, i.e.,

$$\frac{\partial {\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i}{|}_{\mathit{\alpha }=\mathbf{0}}={\left({\mathit{ϵ}}_t\right)}_i.$$

In Equation (51) we make use of the fact that

$$\frac{\partial H_k}{\partial {\mathit{\alpha }}_i}{|}_{\mathit{\alpha }=\mathbf{0}}={\left(\mathbf{Q}\right)}_{k,i}=q_{k,i},$$

which is proved in Appendix A. Thus, it is clear that Equation (51) can be written as,

$${\mathit{\mu }}_{t+1}={\mathit{ϵ}}_t+{\mathit{\mu }}_t\mathbf{Q}.$$

(52)

The above expression states the dynamics for the evolution of the mean distribution ${\mathit{\mu }}_t$ in time. This evolution has two components, one involving the internal dynamics (which is given by the term ${\mathit{\mu }}_t\mathbf{Q}$ giving the internal redistribution of particles) and other one involving the external dynamics (which is given by the time-dependent mean number of incoming particles).

Now let us explore the behavior of the second cumulant of ${\mathbf{N}}^t$. The second cumulant corresponds to the variance matrix $\mathrm{Var}\left({\mathbf{N}}^t\right)$ which we will denote by ${\mathbf{\Sigma }}_t$. Notice that this matrix has entries given by

$${\left({\mathbf{\Sigma }}_t\right)}_{i,j}:={\left(\mathrm{Var}\left({\mathbf{N}}^t\right)\right)}_{i,j}=\mathbb{E}\left[N_i^t{N}_j^t\right]-{\left({\mathit{\mu }}_t\right)}_i{\left({\mathit{\mu }}_t\right)}_j.$$

(53)

Next we will use the dynamics of the m.g.f. of ${\mathbf{N}}^t$ to obtain a recurrence for the expected value $\mathbb{E}\left[N_i^t{N}_j^t\right]$. It is clear that

$$\mathbb{E}\left[N_i^{t+1}{N}_j^{t+1}\right]=\frac{{\partial }^2{\mathcal{G}}_{t+1}\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}{|}_{\mathit{\alpha }=\mathbf{0}}.$$

(54)

The above expression together with the evolution Equation (24) leads to

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i^{t+1}{N}_j^{t+1}\right]& =& \frac{{\partial }^2}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}\left({\mathcal{F}}_t\left(\mathit{\alpha }\right){\mathcal{G}}_t\left(\mathbf{H}\left(\mathit{\alpha }\right)\right)\right){|}_{\mathit{\alpha }=\mathbf{0}}\hfill \\ & =& \frac{\partial }{\partial {\mathit{\alpha }}_i}\left(\frac{\partial {\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_j}{\mathcal{G}}_t\left(\mathbf{H}\left(\mathit{\alpha }\right)\right)+{\mathcal{F}}_t\left(\mathit{\alpha }\right)\sum _{k=1}^S\frac{\partial {\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_k}\frac{\partial H_k}{\partial {\mathit{\alpha }}_j}\right){|}_{\mathit{\alpha }=\mathbf{0}}\hfill \\ & =& (\frac{{\partial }^2{\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}{\mathcal{G}}_t\left(\mathbf{H}\left(\mathit{\alpha }\right)\right)+\frac{\partial {\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_j}\sum _{k=1}^S\frac{\partial {\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_k}\frac{\partial H_k}{\partial {\mathit{\alpha }}_i}\hfill \\ & +& \frac{\partial {\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i}\sum _{k=1}^S\frac{\partial {\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_k}\frac{\partial H_k}{\partial {\mathit{\alpha }}_j}+{\mathcal{F}}_t\left(\mathit{\alpha }\right)\sum _{k=1}^S\sum _{l=1}^S\frac{{\partial }^2{\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_l\partial H_k}\frac{\partial H_k}{\partial {\mathit{\alpha }}_i}\frac{\partial H_l}{\partial {\mathit{\alpha }}_j}\hfill \\ & +& {\mathcal{F}}_t\left(\mathit{\alpha }\right)\sum _{k=1}^S\frac{\partial {\mathcal{G}}_t\left(\mathbf{H}\right)}{\partial H_k}\frac{{\partial }^2{H}_k}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}){|}_{\mathit{\alpha }=\mathbf{0}}.\hfill \end{array}$$

(55)

At this point it is necessary to introduce some notations. First let us denote by $\mathbf{\Delta }$ the variance matrix of the incoming flux, i.e.,

$${\mathbf{\Delta }}_t:=\mathrm{Var}\left({\mathbf{J}}^t\right).$$

(56)

The above quantity can be obtained through the second derivative of the m.g.f. ${\mathcal{F}}_t\left(\mathit{\alpha }\right)$ as follows,

$${\left({\mathbf{\Delta }}_t\right)}_{i,j}=\mathbb{E}\left[J_i^t{J}_j^t\right]-\mathbb{E}\left[J_i^t\right]\mathbb{E}\left[J_j^t\right]=\frac{{\partial }^2{\mathcal{F}}_t\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}-{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j.$$

(57)

Thus, performing the evaluation of Equation (55) at $\mathit{\alpha }=\mathbf{0}$ we obtain,

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i^{t+1}{N}_j^{t+1}\right]& =& {\left({\mathbf{\Delta }}_t\right)}_{i,j}+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j+{\left({\mathit{ϵ}}_t\right)}_j\sum _{k=1}^S{\left({\mathit{\mu }}_t\right)}_k{q}_{k,i}+{\left({\mathit{ϵ}}_t\right)}_i\sum _{k=1}^S{\left({\mathit{\mu }}_t\right)}_k{q}_{k,j}\hfill \\ & +& \sum _{k=1}^S\sum _{l=1}^S\mathbb{E}\left[N_k^t{N}_l^t\right]q_{k,i}{q}_{l,j}+\sum _{k=1}^S{\left({\mathit{\mu }}_t\right)}_k\left(q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}\right),\hfill \end{array}$$

(58)

where we used the fact that (see Appendix A),

$$\frac{{\partial }^2{H}_k\left(\mathit{\alpha }\right)}{\partial {\mathit{\alpha }}_i\partial {\mathit{\alpha }}_j}{|}_{\mathit{\alpha }=\mathbf{0}}=q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}.$$

(59)

Let us simplify Equation (58) by noticing that the summations can be written as matrix products,

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i^{t+1}{N}_j^{t+1}\right]& =& {\left({\mathbf{\Delta }}_t\right)}_{i,j}+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j+{\left({\mathit{ϵ}}_t\right)}_j{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j+{\left({\mathbf{\Lambda }}_t\right)}_{i,j}\hfill \\ & +& \sum _{k=1}^S\sum _{l=1}^S\mathbb{E}\left[N_k^t{N}_l^t\right]q_{k,i}{q}_{l,j},\hfill \end{array}$$

(60)

where we defined the matrix ${\mathbf{\Lambda }}_t$ as

$${\left({\mathbf{\Lambda }}_t\right)}_{i,j}:=\sum _{k=1}^S{\left({\mathit{\mu }}_t\right)}_k\left(q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}\right).$$

(61)

Equation (60) allows us to obtain the variance matrix $\mathrm{Var}\left({\mathbf{N}}^{t+1}\right)$,

$$\begin{array}{ccc}\hfill \mathrm{Var}\left({\mathbf{N}}^{t+1}\right)& =& \mathbb{E}\left[N_i^{t+1}{N}_j^{t+1}\right]-{\left({\mathit{\mu }}_{t+1}\right)}_i{\left({\mathit{\mu }}_{t+1}\right)}_j\hfill \\ & =& {\left({\mathbf{\Delta }}_t\right)}_{i,j}+{\left({\mathbf{\Lambda }}_t\right)}_{i,j}+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j+{\left({\mathit{ϵ}}_t\right)}_j{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j+\hfill \\ & +& \sum _{k=1}^S\sum _{l=1}^S\left(\mathbb{E}\left[N_k^t{N}_l^t\right]-{\left({\mathit{\mu }}_t\right)}_k{\left({\mathit{\mu }}_t\right)}_l\right)q_{k,i}{q}_{l,j}+\sum _{k=1}^S\sum _{l=1}^S{\left({\mathit{\mu }}_t\right)}_k{\left({\mathit{\mu }}_t\right)}_l{q}_{k,i}{q}_{l,j}-{\left({\mathit{\mu }}_{t+1}\right)}_i{\left({\mathit{\mu }}_{t+1}\right)}_j=\hfill \\ & =& {\left({\mathbf{\Delta }}_t\right)}_{i,j}+{\left({\mathbf{\Lambda }}_t\right)}_{i,j}+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j+{\left({\mathit{ϵ}}_t\right)}_j{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j\hfill \\ & +& \sum _{k=1}^S\sum _{l=1}^S{q}_{k,i}{\left(\mathrm{Var}\left({\mathbf{N}}^t\right)\right)}_{k,l}{q}_{l,j}+{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j-{\left({\mathit{\mu }}_{t+1}\right)}_i{\left({\mathit{\mu }}_{t+1}\right)}_j.\hfill \end{array}$$

(62)

Rearranging terms in the above expression and denoting by ${\mathbf{\Sigma }}_t$ the matrix variance $\mathrm{Var}\left({\mathbf{N}}^t\right)$, we obtain

$$\begin{array}{ccc}\hfill {\left({\mathbf{\Sigma }}_{t+1}\right)}_{i,j}& =& {({\mathbf{\Delta }}_t+{\mathbf{\Lambda }}_t)}_{i,j}+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{ϵ}}_t\right)}_j+{\left({\mathit{ϵ}}_t\right)}_j{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i+{\left({\mathit{ϵ}}_t\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j+{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_i{\left({\mathit{\mu }}_t\mathbf{Q}\right)}_j\hfill \\ & +& \sum _{k=1}^S\sum _{l=1}^S{q}_{k,i}{\left({\mathbf{\Sigma }}_t\right)}_{k,l}{q}_{l,j}-{\left({\mathit{\mu }}_{t+1}\right)}_i{\left({\mathit{\mu }}_{t+1}\right)}_j\hfill \\ & =& {({\mathbf{\Delta }}_t+{\mathbf{\Lambda }}_t)}_{i,j}+{\left({\mathit{ϵ}}_t+{\mathit{\mu }}_t\mathbf{Q}\right)}_i{\left({\mathit{ϵ}}_t+{\mathit{\mu }}_t\mathbf{Q}\right)}_j+{\left({\mathbf{Q}}^T{\mathbf{\Sigma }}_t\mathbf{Q}\right)}_{i,j}-{\left({\mathit{\mu }}_{t+1}\right)}_i{\left({\mathit{\mu }}_{t+1}\right)}_j.\hfill \\ \hfill \end{array}$$

(63)

Finally, observing that ${\mathit{ϵ}}_t+{\mathit{\mu }}_t\mathbf{Q}={\mathit{\mu }}_{t+1}$, it is clear that the variance matrix satisfy the evolution equation,

$${\mathbf{\Sigma }}_{t+1}={\mathbf{\Delta }}_t+{\mathbf{\Lambda }}_t+{\mathbf{Q}}^T{\mathbf{\Sigma }}_t\mathbf{Q}.$$

(64)

Equations (52) and (64) govern the dynamics of the first and second cumulants and are valid even when the incoming flux has a time-dependent distribution. In the case where the protocol of incoming particles $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ is a stationary process (for which the two first cumulants are time-independent) we have that the system can reach the stationarity. Particularly we have that the dynamics Equations (52) and (64) have stationary solutions (proved in Appendix A) given by,

$$\begin{array}{ccc}\hfill \overline{\mathit{\mu }}& =& \mathit{ϵ}{(\mathbf{1}-\mathbf{Q})}^{-1}\hfill \end{array}$$

(65)

$$\begin{array}{ccc}\hfill \overline{\mathbf{\Sigma }}& =& \sum _{k=0}^{\infty }{\left({\mathbf{Q}}^T\right)}^k(\mathbf{\Delta }+\overline{\mathbf{\Lambda }}){\mathbf{Q}}^k,\hfill \end{array}$$

(66)

where $\mathit{ϵ}$ and $\mathbf{\Delta }$ are the mean vector and the variance matrix of the stationary process $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ and $\overline{\mathit{\mu }}$ and $\overline{\mathbf{\Sigma }}$ denote the mean distribution $\mathbb{E}\left[{\mathbf{N}}^t\right]$ and the variance matrix $\mathrm{Var}\left({\mathbf{N}}^t\right)$ when the process $\{{\mathbf{N}}^t:t\in {\mathbb{N}}_0\}$ has reached the stationarity. We also defined $\overline{\mathbf{\Lambda }}$ as

$${\left(\overline{\mathbf{\Lambda }}\right)}_{i,j}=\sum _{k=0}^S{\left(\overline{\mathit{\mu }}\right)}_k\left(q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}\right).$$

(67)

Example 5.

Let us consider a three states open chain with jump matrix given by

$$\mathbf{Q}=\left[\begin{array}{ccc}0& q& q\\ q& 0& q\\ q& q& 0\end{array}\right]$$

(68)

with p a parameter restricted to take values in the interval $0<q<1/2$. Let us also assume that the protocol of incoming particles $\{{\mathbf{J}}^t=(J_1^t,J_2^t,J_3^t):t\in \mathbb{N}\}$ is a stationary process whose joint probability distribution at time t can be written as

$$\mathbb{P}(J_1^t=j_1;J_2^t=j_2;J_3^t=j_3)=:f_{1,2}(j_1,j_2)f_3\left(j_3\right).$$

Particularly, we chose $f_{1,2}$ and $f_3$ as

$$\begin{array}{ccc}\hfill f_{1,2}(j_1,j_2)& =& \left\{\begin{array}{ccc}\hfill p/2& if& (j_1,j_2)=(1,1)\hfill \\ \hfill p/2& if& (j_1,j_2)=(0,0)\hfill \\ \hfill (1-p)/2& if& (j_1,j_2)=(1,0)\hfill \\ \hfill (1-p)/2& if& (j_1,j_2)=(0,1)\hfill \\ \hfill 0& if& \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(69)

$$\begin{array}{ccc}\hfill f_3\left(j_3\right)& =& \left\{\begin{array}{ccc}\hfill 1/2& if& j_3=0\hfill \\ \hfill 1/2& if& j_3=1\hfill \\ \hfill 0& if& \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(70)

Notice that each random variable $J_i^t$ can only take the values 0 or 1. Moreover, the above choice for the joint probability distribution for the random vector $(J_1^t,J_2^t,J_3^t)$ implies that the random variables $J_1^t$ and $J_2^t$ are dependent and that $J_3^t$ is independent. It is easy to see that all these random variable have, separately, a Bernoulli distribution with parameter $1/2$, i.e., at every time-step one particle arrives at every node with probability $1/2$ and no particles arrive at every node with probability $1/2$. In Figure 4 we show a graphical representation of the open chain.

A straightforward calculation shows that the first and second moments of ${\mathbf{J}}^t$ are given by

$$\begin{array}{ccc}\hfill \mathbb{E}\left[{\mathbf{J}}^t\right]& =& \mathit{ϵ}=(1/2,1/2,1/2),\hfill \end{array}$$

(71)

$$\begin{array}{ccc}\hfill \mathrm{Var}\left({\mathbf{J}}^t\right)& =& \mathbf{\Delta }=\frac{1}{4}\left[\begin{array}{ccc}1& 2p-1& 0\\ 2p-1& 1& 0\\ 0& 0& 1\end{array}\right].\hfill \end{array}$$

(72)

Now we calculate the mean stationary distribution of particles, $\overline{\mathit{\mu }}$, and the stationary variance matrix, $\overline{\mathbf{\Sigma }}$. The calculation of $\overline{\mathit{\mu }}$ is straightforward. We only need to obtain the inverse of the matrix $\mathbf{1}-\mathbf{Q}$, which is given by

$${(\mathbf{1}-\mathbf{Q})}^{-1}=\frac{1}{1-3q^2-2q^3}\left[\begin{array}{ccc}1-q^2& q+q^2& q+q^2\\ q+q^2& 1-q^2& q+q^2\\ q+q^2& q+q^2& 1-q^2\end{array}\right].$$

With the above result we can see that the mean stationary distribution $\overline{\mathit{\mu }}$ is given by

$$\overline{\mathit{\mu }}=\mathit{ϵ}{(\mathbf{1}-\mathbf{Q})}^{-1}=\left(\frac{1}{2-4q},\frac{1}{2-4q},\frac{1}{2-4q}\right).$$

The last result implies that the particles on the state space are uniformly distributed when the system has reached the stationarity. Moreover, we can also appreciate that the number of particles on every state diverges as the parameter p tends to $1/2$. This divergence is actually a consequence of the fact that the system does not allow that the particles leave the state space when $p=1/2$. Thus, for such a parameter value, the system is still receiving particles but it does not allow that the particles escape to the outside, thus increasing indefinitely the number of particles inside the system.

Now let us compute the stationary variance matrix $\overline{\mathbf{\Sigma }}$. Recall that this quantity can be obtained by means of the formula,

$$\overline{\mathbf{\Sigma }}=\sum _{k=0}^{\infty }{\left({\mathbf{Q}}^T\right)}^k(\mathbf{\Delta }+\overline{\mathbf{\Lambda }}){\mathbf{Q}}^k.$$

(73)

First let us obtain the explicit form of $\overline{\mathbf{\Lambda }}$. According to Equation (67) we have that

$$\begin{array}{ccc}\hfill {\left(\overline{\mathbf{\Lambda }}\right)}_{i,j}& =& \sum _{k=1}^S{\left(\overline{\mathit{\mu }}\right)}_k\left(q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}\right)\hfill \\ & =& \frac{1}{2-4q}\sum _{k=1}^3\left(q_{k,i}{\delta }_{i,j}-q_{k,i}{q}_{k,j}\right)\hfill \\ & =& \frac{q}{1-2q}{\delta }_{i,j}-\frac{1}{2-4q}\sum _{k=1}^S{q}_{k,i}{q}_{k,j}.\hfill \end{array}$$

Thus, we have

$$\overline{\mathbf{\Lambda }}=\frac{1}{2-4q}\left[\begin{array}{ccc}2q(1-q)& -q^2& -q^2\\ -q^2& 2q(1-q)& -q^2\\ -q^2& -q^2& 2q(1-q)\end{array}\right].$$

Next we need to compute the nth power of the matrix $\mathbf{Q}$. It is not hard to prove that

$${\mathbf{Q}}^k=\frac{q^k}{3}\left[\begin{array}{ccc}{2}^k+2{(-1)}^k& 2^k-{(-1)}^k& 2^k-{(-1)}^k\\ 2^k-{(-1)}^k& 2^k+2{(-1)}^k& 2^k-{(-1)}^k\\ 2^k-{(-1)}^k& 2^k-{(-1)}^k& 2^k+2{(-1)}^k\end{array}\right].$$

(74)

This result shows that for this particular case

$${\left({\mathbf{Q}}^T\right)}^k={\mathbf{Q}}^k,$$

because $\mathbf{Q}$ is itself a symmetric matrix.

It is also not hard to see that the term

$${\mathbf{Q}}^k(\mathbf{\Lambda }+\mathbf{\Delta }){\left({\mathbf{Q}}^T\right)}^k$$

is a matrix whose components are all exponential (or linear combinations of exponentials) in the variable k. This observation allows us to see that the infinite summation (73) can be exactly computed. Then, we can obtain a closed expression for $\overline{\mathbf{\Sigma }}$ by using symbolic calculations performed in the software Mathematica. Thus we obtain,

$$\begin{array}{ccc}\hfill \overline{\mathbf{\Sigma }}& =& \left[\begin{array}{ccc}\frac{-8q^5+(8p-2)q^4+4q^3+3q^2+4q+1}{4\left(8q^6-6q^4-3q^2+1\right)}& \frac{2q^4+q^2+p\left(-8q^4-4q^2+2\right)-1}{4\left(8q^6-6q^4-3q^2+1\right)}& -\frac{q^2\left(q^2-p+1\right)}{2\left(8q^6-6q^4-3q^2+1\right)}\\ \frac{2q^4+q^2+p\left(-8q^4-4q^2+2\right)-1}{4\left(8q^6-6q^4-3q^2+1\right)}& \frac{-8q^5+(8p-2)q^4+4q^3+3q^2+4q+1}{4\left(8q^6-6q^4-3q^2+1\right)}& -\frac{q^2\left(q^2-p+1\right)}{2\left(8q^6-6q^4-3q^2+1\right)}\\ -\frac{q^2\left(q^2-p+1\right)}{2\left(8q^6-6q^4-3q^2+1\right)}& -\frac{q^2\left(q^2-p+1\right)}{2\left(8q^6-6q^4-3q^2+1\right)}& \frac{-8q^5+(6-8p)q^4+4q^3+(4p+1)q^2+4q+1}{4\left(8q^6-6q^4-3q^2+1\right)}\end{array}\right]\hfill \end{array}$$

Now, let us define the space correlation functions. We will denote by ${\kappa }_{i,j}$ the correlation function between the random variables $N_i^t$ and $N_j^t$ as follows,

$$\begin{array}{ccc}\hfill {\kappa }_{i,j}& :=& \mathrm{Corr}(N_i^t,N_j^t)=\frac{{\overline{\mathbf{\Sigma }}}_{i,j}}{\sqrt{{\overline{\mathbf{\Sigma }}}_{i,i}{\overline{\mathbf{\Sigma }}}_{j,j}}},\hfill \end{array}$$

(75)

for all $i,j\in \mathcal{S}$. It is not hard to see that the correlation functions for this example are given by,

$$\begin{array}{ccc}\hfill {\kappa }_{1,2}& =& \frac{p\left(-8q^4-4q^2+2\right)+2q^4+q^2-1}{\sqrt{{\left((8p-2)q^4-8q^5+4q^3+3q^2+4q+1\right)}^2}},\hfill \end{array}$$

(76)

$$\begin{array}{ccc}\hfill {\kappa }_{1,3}& =& -\frac{2q^2\left(-p+q^2+1\right)}{\sqrt{\left((6-8p)q^4+(4p+1)q^2-8q^5+4q^3+4q+1\right)}}\frac{1}{\sqrt{\left((8p-2)q^4-8q^5+4q^3+3q^2+4q+1\right)}},\hfill \end{array}$$

(77)

$$\begin{array}{ccc}\hfill {\kappa }_{2,3}& =& \frac{-2q^2\left(-p+q^2+1\right)}{\sqrt{\left((6-8p)q^4+(4p+1)q^2-8q^5+4q^3+4q+1\right)}}\frac{1}{\sqrt{\left((8p-2)q^4-8q^5+4q^3+3q^2+4q+1\right)}}.\hfill \end{array}$$

(78)

In order to test the results we obtained for this example we performed numerical simulations. We have simulated the dynamics of the open chain for several parameter values during a time of $5\times {10}^5$ steps. Then, from the time series obtained we estimated the correlation functions ${\kappa }_{1,2}$, ${\kappa }_{1,3}$ and ${\kappa }_{2,3}$ and we compare them against the theoretical prediction given in Equations (76)–(78) (note that for our example ${\kappa }_{1,3}={\kappa }_{2,3}$). In Figure 5a,b we show the correlation functions obtained by numerical simulations and computed from Equations (76) and (78). We take the parameter value $q=0.25$ and plotted ${\kappa }_{1,2}$ (solid line, filled circles) and ${\kappa }_{2,3}$ (dashed line, filled squares) as a function of p (Figure 5a). The same graph is done but using the parameter value $q=0.45$ (Figure 5b). As we can see, our the results obtained from numerical simulations are consistent with the formulas that we obtained theoretically.

3.3. Distribution of Particles Leaving the State Space

Up to now we have given an expression for the m.g.f. for the number of particles in the state space. Since the system is open, at every time step there is a number of particles arriving to the system, which is determined by the random vector ${\mathbf{J}}^t$. The total number of particles per step arriving to the system is then

$$I_t=\sum _{i=1}^s{J}_i^t.$$

When the system reaches the stationarity the mean number of particles within the system attains a constant value, meaning that there is an equilibration between the number of incoming and outgoing particles. This is a consequence of the fact that we assumed that the jump matrix is irreducible and aperiodic.

To be precise, if we denote by $O_t$ the total number of particles leaving the state space, then, under stationarity we should have that

$$\mathbb{E}\left[I_t\right]=\mathbb{E}\left[O_t\right],$$

a fact that can be inferred by the “conservation” of the mean number of particles at stationarity. Our goal here is to go beyond the above expression, we would like to characterize how the random variable $O_t$ evolves in time and how much it is influenced by the incoming number of particles and the “jumping” rules of the Markov chain. To this end, let us define some quantities which will allow to describe the random variable $O_t$ explicitly.

Definition 2.

Let $\left(\mathcal{S},\mathbf{Q},\{{\mathbf{J}}^t:t\in \mathbb{N}\}\right)$ be an open Markov chain with state space $\mathcal{S}$, jump matrix $\mathbf{Q}$ (with components ${\left(\mathbf{Q}\right)}_{i,j}=q_{i,j}$) and incoming protocol $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$. Let $e_i$ be defined as the escape probability from the state i, i.e., the probability with which a particle in the state i leaves the system to the outside,

$$e_i:=1-\sum _{j=1}^S{q}_{i,j}.$$

(79)

Next, let ${\mathbf{U}}^t=(U_1^t,U_2^t,\dots ,U_S^t)$ be a random vector whose components have binomial distribution as follows,

$$U_i^t\sim \mathrm{Binom}(N_i^t,e_i)$$

where ${\mathbf{N}}^t$ is time-dependent distribution over the state space. Then we say the $U_i^t$ is the number of particles leaving the state i to the outside at time t. The total number of particles $O_t$ leaving the system is then the random variable given by

$$O_t:=\sum _{i=1}^S{U}_i^t.$$

(80)

Finally, we define the vector $\mathbf{e}$, with components $e_i={\left(\mathbf{e}\right)}_i$, which will be referred to as the escape probability vector of the chain. Additionally let $\mathbf{E}$ be a diagonal matrix with components ${\left(\mathbf{E}\right)}_{i,j}=E_{i,j}$ defined as

$$E_{i,j}:=e_i{\delta }_{i,j},$$

(81)

which will be referred to as the escape probability matrix of the chain.

Our main goal is now to characterize the random vector ${\mathbf{U}}^t$ giving the number of particles leaving the chain. We should notice that the distribution of ${\mathbf{U}}^t$ depends on ${\mathbf{N}}^t$ which is also a random vector. This implies that, given the value of ${\mathbf{N}}^t$, we can specify the conditional distribution for ${\mathbf{U}}^t$, i.e.,

$$T(\mathbf{n};\mathbf{m}):=\mathbb{P}\left({\mathbf{U}}^t=\mathbf{m}|{\mathbf{N}}^t=\mathbf{n}\right),$$

(82)

which is given by

$$T(\mathbf{n};\mathbf{m})=\prod _{i=1}^S\frac{n_i!}{(n_i-m_i)!m_i!}{e}_i^{m_i}{(1-e_i)}^{n_i-m_i}.$$

(83)

Once we know $T(\mathbf{n};\mathbf{m})$, we can write an expression for the probability distribution for the random vector ${\mathbf{U}}^t$,

$$\mathbb{P}\left({\mathbf{U}}^t=\mathbf{m}\right)=\sum _{\mathbf{n}\in {\mathbb{N}}_0^S}\mathbb{P}\left({\mathbf{U}}^t=\mathbf{m}|{\mathbf{N}}^t=\mathbf{n}\right)\mathbb{P}\left({\mathbf{N}}^t=\mathbf{n}\right).$$

(84)

If we denote by $r_t\left(\mathbf{m}\right)$ the probability distribution of the random vector ${\mathbf{U}}^t$, i.e.,

$$r_t\left(\mathbf{n}\right):=\mathbb{P}\left({\mathbf{U}}^t=\mathbf{m}\right)$$

then it is clear that Equation (84) can be rewritten as

$$r_t\left(\mathbf{m}\right)=\sum _{n\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{n}\right)T(\mathbf{n};\mathbf{m}).$$

(85)

The next step consists in obtaining the moment generating function of ${\mathbf{U}}^t$. This is because, as we saw above, the expression for the m.g.f. of ${\mathbf{N}}_t$ can be written explicitly. Thus, let ${\mathcal{R}}_t:{\mathbb{R}}^S\to \mathbb{R}$ be the m.g.f. of ${\mathbf{U}}^t$,

$${\mathcal{R}}_t\left(\mathit{\alpha }\right):=\mathbb{E}\left[e^{{\mathbf{U}}^t{\mathit{\alpha }}^T}\right]=\sum _{\mathbf{m}\in {\mathbb{N}}_0^S}{r}_t\left(\mathbf{m}\right)e^{\mathbf{m}{\mathit{\alpha }}^T}.$$

(86)

Thus we can see that, using Expression (85), ${\mathbf{R}}_t$ can be written as follows

$${\mathcal{R}}_t\left(\mathit{\alpha }\right)=\sum _{\mathbf{m}\in {\mathbb{N}}_0^S}\sum _{n\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{n}\right)T(\mathbf{n};\mathbf{m})e^{\mathbf{m}{\mathit{\alpha }}^T}.$$

(87)

Since $T(\mathbf{n};\mathbf{m})$ is a product of binomial distributions, it is clear that the sum over $\mathbf{m}$ results in the product of moment generating functions of binomial random variables,

$$\sum _{m\in {\mathbb{N}}_0^S}T(\mathbf{n};\mathbf{m})e^{\mathbf{m}{\mathit{\alpha }}^T}=\prod _{i=1}^S{\left(1-e_i+e_i{e}^{{\mathit{\alpha }}_i}\right)}^{n_i}.$$

(88)

Moreover, if we define the function $\mathbf{C}=(C_1,C_2,\dots ,C_S):{\mathbb{R}}^S\to {\mathbb{R}}^S$ as follows,

$$C_i\left(\mathit{\alpha }\right)=C_i\left({\mathit{\alpha }}_i\right)=log\left(1-e_i+e_i{e}^{{\mathit{\alpha }}_i}\right),\mathrm{for}1\le i\le S,$$

(89)

it is clear that we can write

$$\sum _{m\in {\mathbb{N}}_0^S}T(\mathbf{n};\mathbf{m})e^{\mathbf{m}{\mathit{\alpha }}^T}=exp\left(\mathbf{n}{\mathbf{C}}^T\left(\mathit{\alpha }\right)\right).$$

(90)

The above expression, together with Equation (87), gives us

$${\mathcal{R}}_t\left(\mathit{\alpha }\right)={\mathcal{G}}_t\left(\mathbf{C}\left(\mathit{\alpha }\right)\right).$$

(91)

The last relation states that the moment generating function of ${\mathbf{U}}^t$ can be obtained by means of the moment generating function of ${\mathbf{N}}^t$, which is mediated by the transformation $\mathbf{C}$.

Observe that ${\mathcal{R}}_t\left(\mathit{\alpha }\right)$ allows us to obtain the first and the second moment of number of leaving particles. Taking the first derivative (gradient) of ${\mathcal{R}}_t\left(\mathit{\alpha }\right)$ and evaluating it in $\mathit{\alpha }=\mathbf{0}$ we obtain,

$$\mathbb{E}\left[{\mathbf{U}}^t\right]={\mathit{\mu }}_t\mathbf{E}=({\left({\mathit{\mu }}_t\right)}_1{e}_1,{\left({\mathit{\mu }}_t\right)}_2{e}_2,\dots ,{\left({\mathit{\mu }}_t\right)}_S{e}_S).$$

(92)

The above result allows us to calculate the mean number of leaving particles $\mathbb{E}\left[O_t\right]$,

$$\mathbb{E}\left[O_t\right]=\sum _{i=1}^S\mathbb{E}\left[U_i^t\right]=\sum _{i=1}^S{\left({\mathit{\mu }}_t\right)}_i{e}_i={\mathit{\mu }}_t{\mathbf{e}}^T.$$

(93)

Now, when the system reaches stationarity we have to replace ${\mathit{\mu }}_t$ by $\overline{\mathit{\mu }}$,

$${\left(\mathbb{E}\left[O_t\right]\right)}_{\mathrm{stat}}=\sum _{i=1}^S{\left(\overline{\mathit{\mu }}\right)}_i{e}_i.$$

Thus, from the definition of $e_i$ we obtain

$$\begin{array}{ccc}\hfill {\left(\mathbb{E}\left[O_t\right]\right)}_{\mathrm{stat}}& =& \sum _{i=1}^S{\left(\overline{\mathit{\mu }}\right)}_i\left(1-\sum _{j=1}^S{q}_{i,j}\right)=\sum _{i=1}^S{\left(\overline{\mathit{\mu }}\right)}_i-\sum _{j=1}^S\sum _{i=1}^S{\left(\overline{\mathit{\mu }}\right)}_i{q}_{i,j}\hfill \\ & =& \sum _{i=1}^S{\left(\overline{\mathit{\mu }}\right)}_i-\sum _{j=1}^S{\left(\overline{\mathit{\mu }}\mathbf{Q}\right)}_j=\sum _{i=1}^S{\left(\overline{\mathit{\mu }}-\overline{\mathit{\mu }}\mathbf{Q}\right)}_i,\hfill \end{array}$$

(94)

and recalling that $\overline{\mathit{\mu }}$ satisfy the equation

$$\overline{\mathit{\mu }}-\overline{\mathit{\mu }}\mathbf{Q}=\mathit{ϵ},$$

we can observe that

$$\begin{array}{ccc}\hfill {\left(\mathbb{E}\left[O_t\right]\right)}_{\mathrm{stat}}& =& \sum _{i=1}^S{\left(\mathit{ϵ}\right)}_i=\sum _{i=1}^S{\left(\mathbb{E}\left[{\mathbf{J}}^t\right]\right)}_i={\left(\mathbb{E}\left[I_t\right]\right)}_{\mathrm{stat}},\hfill \end{array}$$

(95)

which is the relation we have anticipated by invoking the particle number conservation principle.

The variance matrix $\mathrm{Var}\left({\mathbf{U}}^t\right)$ can also be obtained by means of the second derivative of ${\mathcal{R}}_t\left(\mathit{\alpha }\right)$. A calculation achieved in Appendix A shows that

$$\mathrm{Var}\left({\mathbf{U}}^t\right)=\mathbf{E}{\mathbf{\Sigma }}_t\mathbf{E}+{\mathbf{D}}_t,$$

(96)

where $\mathbf{D}$ is a diagonal matrix with components

$${\left({\mathbf{D}}_t\right)}_{i,j}:={\left({\mathit{\mu }}_t\right)}_i{e}_i(1-e_i){\delta }_{i,j},$$

where ${\delta }_{i,j}$ is the well-known Krönecker delta. The variance matrix ${\mathbf{U}}^t$ allows us in turn to obtain a closed formula for the variance of the total number of leaving particles,

$$\begin{array}{ccc}\hfill \mathrm{Var}\left(O_t\right)& =& \sum _{i=1}^S\sum _{j=1}^S{\left(\mathrm{Var}\left({\mathbf{U}}^t\right)\right)}_{i,j}=\sum _{i=1}^S\sum _{j=1}^S{e}_i{\left({\mathbf{\Sigma }}_t\right)}_{i,j}{e}_j+\sum _{i=1}^S{\left({\mathit{\mu }}_t\right)}_i(1-e_i)e_i\hfill \\ & =& \mathbf{e}{\mathbf{\Sigma }}_t{\mathbf{e}}^T+{\mathit{\mu }}_t\left(\mathbf{1}-\mathbf{E}\right){\mathbf{e}}^T,\hfill \end{array}$$

(97)

where $\mathbf{1}$ stands for the $S\times S$ identity matrix. We observe that the variance matrix of the number of outgoing particles is not the same as the variance matrix of the incoming particles nor the variance matrix of the particles in the system. Thus we have that the fluctuations and correlations are modulated by the internal dynamics of the system. This is important because measuring the correlations and fluctuations of the outgoing flux gives information on the internal dynamics of the system.

Example 6.

Let us consider again the one-vertex model given in Example 1. We have seen that the m.g.f. of ${\mathbf{N}}^t$ is given by

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{\infty }\left(1-p{q}^r+p{q}^r{e}^{\mathit{\alpha }}\right).\hfill \end{array}$$

(98)

Notice that the probability escape “vector” consists of a single number, given by $e_0=1-q$ (recall that the jump matrix is also a single number). Thus we have that the transformation $C:\mathbb{R}\to \mathbb{R}$ can be written as

$$C\left(\mathit{\alpha }\right)=log\left(1-e_0+e_0{e}^{\mathit{\alpha }}\right)=log\left(q+(1-q)e^{\mathit{\alpha }}\right).$$

Now we can obtain the m.g.f. of $U^t$, the number of particles leaving the system,

$${\mathcal{R}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)={\mathcal{G}}_{\mathrm{stat}}\left(C\left(\mathit{\alpha }\right)\right)=\prod _{r=0}^{\infty }\left(1-p{q}^r+p{q}^r{e}^{C\left(\mathit{\alpha }\right)}\right).$$

Notice that $e^{C\left(\mathit{\alpha }\right)}$ can be written as $e^{C\left(\mathit{\alpha }\right)}=q+(1-q)e^{\mathit{\alpha }}$, thus we have,

$$\begin{array}{ccc}\hfill {\mathcal{R}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{\infty }\left(1-p{q}^r+p{q}^r\left(q+(1-q)e^{\mathit{\alpha }}\right)\right).\hfill \\ & =& \prod _{r=0}^{\infty }\left(1-p{q}^r(1-q)+p{q}^r(1-q)e^{\mathit{\alpha }}\right).\hfill \end{array}$$

The last expression establishes that the distribution of $U^t$ (at stationarity) can be seen as an infinite convolution of Bernoulli distributions, with parameters $p{q}^r(1-q)$ for $r\in {\mathbb{N}}_0$. Thus, the random variable $U^t$ can be written as an infinite sum of i.i.d. random variables $Y_r$ (with the above-mentioned Bernoulli distribution),

$$U^t=\sum _{r=0}^{\infty }{Y}_r.$$

whenever $U^t$ has reached stationarity. Particularly the mean number of leaving particles as well as its variance can be exactly determined,

$$\begin{array}{ccc}\hfill \mathbb{E}\left[U^t\right]& =& \sum _{r=0}^{\infty }\mathbb{E}\left[Y^r\right]=\sum _{r=0}^{\infty }p{q}^r(1-q)=p.\hfill \\ \hfill \mathrm{Var}\left(U^t\right)& =& \sum _{r=0}^{\infty }\mathrm{Var}\left(Y_r\right)=\sum _{r=0}^{\infty }p{q}^r(1-q)(1-p{q}^r(1-q))\hfill \end{array}$$

(99)

$$\begin{array}{ccc}& =& p-\frac{p^2{(1-q)}^2}{1-q^2}.\hfill \end{array}$$

(100)

We should emphasize that the above expressions can be obtained through the formulas for $\mathbb{E}\left[O_t\right]$ and $\mathrm{Var}\left(O_t\right)$, given in Equations (93) and (97). Since $S=1$ (because we have only one state) we have that

$$\begin{array}{ccc}\hfill \mathbb{E}\left[O^t\right]& =& \mathbb{E}\left[N^t\right]e_0=\frac{p}{1-q}{e}_0=p,\hfill \\ \hfill \mathrm{Var}\left(O^t\right)& =& \mathrm{Var}\left(N^t\right)e_0^2+\mathbb{E}\left[N^t\right]e_0(1-e_0).\hfill \\ & =& \left(\frac{p}{1-q}+\frac{p^2}{1-q^2}\right){(1-q)}^2+\left(\frac{p}{1-q}\right)q(1-q)\hfill \\ & =& p+\frac{p^2{(1-q)}^2}{1-q^2}.\hfill \end{array}$$

(101)

where we used the expressions for $\mathbb{E}\left[N^t\right]$ and $\mathrm{Var}\left(N^t\right)$ given in Equations (29) and (30).

4. The Influence of Incoming Particles on Time-Correlations

Time Correlations for the Open Markov Chain

Up to now we have seen that it is possible to find an explicit expression for the time-dependent distribution on the state space. This evolution is fully characterized by the two first cumulants, the mean distribution over the state space ${\mathit{\mu }}_t$ and the variance matrix ${\mathbf{\Sigma }}_t$. It is important to emphasize that we made no assumptions on the time correlations of the sequence of random vectors $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$. This is because to obtain the distribution over the state space, ${\mathbf{N}}^t$, for a given time t, it was enough to know the number of incoming particles at time t. On the other hand, if we would like to compute the two-times correlation function for certain observable we necessarily have to known the number of incoming particles a two different times. This information unavoidably will be related to the two-times covariance matrix of the process $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$, i.e., the covariance matrix between the random vector ${\mathbf{J}}^t$ and ${\mathbf{J}}^{t+s}$ for $s,t\in \mathbb{N}$. Our goal in this section is to obtain an expression for the covariance $C_{i,j}(t,t+s)$ between the ith coordinate of ${\mathbf{N}}^t$ and the jth coordinate of ${\mathbf{N}}^{t+s}$, i.e.,

$$C_{i,j}(t,t+s):=\mathbb{E}\left[N_i^t{N}_j^{t+s}\right]-{\left({\mathit{\mu }}^t\right)}_i{\left({\mathit{\mu }}^{t+s}\right)}_j,$$

(102)

where ${\overline{\mathit{\mu }}}_i$ is the ith coordinate of the mean stationary distribution.

In order to compute the expected value $\mathbb{E}\left[N_i^t{N}_j^{t+s}\right]$ it is necessary to have an expression for the stationary joint distribution $P_{t,t+s}(\mathbf{n},\mathbf{m})$. This quantity is defined as,

$$P_{t,t+s}(\mathbf{n},\mathbf{m}):=\mathbb{P}\left({\mathbf{N}}^t=\mathbf{n};{\mathbf{N}}^{t+s}=\mathbf{m}\right),$$

(103)

Our goal here is to establish a method to obtain the joint distribution $P_t(\mathbf{n},\mathbf{m})$. Actually, we will first determine an expression for a more general quantity. Let $P({\mathbf{n}}_0,{\mathbf{n}}_1,{\mathbf{n}}_s)$ denote the joint probability function of the random vectors ${\mathbf{N}}^t,{\mathbf{N}}^{t+1},\dots ,{\mathbf{N}}^{t+s}$ (Notice that, to be strict, the probability function P depends on $t,t+1,\dots ,t+s$, and we should denote this dependence explicitly by using subscripts, i.e., $P=P_{t,t+1,\dots ,t+1}$. However, we will not use such a notation by the sake of simplicity in further calculations. The same convention will be adopted for other “multiple-times” joint probability functions or its corresponding moment generating functions.),

$$P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s):=\mathbb{P}\left({\mathbf{N}}^t={\mathbf{n}}_0;{\mathbf{N}}^{t+1}={\mathbf{n}}_1;\dots ;{\mathbf{N}}^{t+s}={\mathbf{n}}_s\right).$$

(104)

First of all, let us introduce some notation that will be useful to perform further calculations. Let $f({\mathbf{j}}_0,{\mathbf{j}}_1,\dots ,{\mathbf{j}}_{s-1})$ be the joint probability function of the random vectors ${\mathbf{J}}^t,{\mathbf{J}}^{t+1},\dots ,{\mathbf{J}}^{t+s-1}$, i.e.,

$$f({\mathbf{j}}_0,{\mathbf{j}}_1,\dots ,{\mathbf{j}}_{s-1}):=\mathbb{P}\left({\mathbf{J}}^t={\mathbf{j}}_0;{\mathbf{J}}^{t+1}={\mathbf{j}}_1;\dots ,{\mathbf{J}}^{t+s-1}={\mathbf{j}}_{s-1}\right).$$

(105)

Let us also denote by $h(\mathbf{r};\mathbf{k})$ the probability function of the random vector ${\mathbf{R}}^t$, which, as we saw in Section 3, depends on the value taken by the random vector ${\mathbf{N}}^t$ (a value which we denote by $\mathbf{k}$ in the probability function h).

With the above-introduced notation it is possible to write the joint probability function, given in Equation (104), in terms of the probability functions f and h,

$$\begin{array}{ccc}\hfill P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s)& =& \mathbb{P}\left({\mathbf{N}}^t={\mathbf{n}}_0;{\mathbf{N}}^{t+1}={\mathbf{n}}_1;\dots ;{\mathbf{N}}^{t+s}={\mathbf{n}}_s\right)\hfill \\ \hfill & =& \mathbb{P}\left({\mathbf{N}}^{t+1}={\mathbf{n}}_1;{\mathbf{N}}^{t+2}={\mathbf{n}}_2;\dots ;{\mathbf{N}}^{t+s}={\mathbf{n}}_s|{\mathbf{N}}^t={\mathbf{n}}_0\right)\mathbb{P}({\mathbf{N}}^t={\mathbf{n}}_0)\hfill \\ \hfill & =& \mathbb{P}\left({\mathbf{J}}^t+{\mathbf{R}}^t={\mathbf{n}}_1;\dots ;{\mathbf{J}}^{t+s-1}+{\mathbf{R}}^{t+s-1}={\mathbf{n}}_s|{\mathbf{N}}^t={\mathbf{n}}_0\right)p_t\left({\mathbf{n}}_0\right).\hfill \end{array}$$

Notice that the random vector ${\mathbf{R}}^{t+j}$ depends on ${\mathbf{n}}_j$ for $1\le j\le s-1$, values which are given a priori. Thus, the random vectors ${\mathbf{R}}^t,{\mathbf{R}}^{t+1},\dots ,{\mathbf{R}}^{t+s-1}$ are all independent (because the values taken by the random vectors ${\mathbf{N}}^{t+j}$ for $0\le j\le s$ are all fixed), which allows us to write

$$\begin{array}{ccc}\hfill P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s)& =& \sum _{{\mathbf{j}}_0+{\mathbf{r}}_0={\mathbf{n}}_1}\sum _{{\mathbf{j}}_1+{\mathbf{r}}_1={\mathbf{n}}_2}\dots \sum _{{\mathbf{j}}_{s-1}+{\mathbf{r}}_{s-1}={\mathbf{n}}_s}\mathbb{P}\left({\mathbf{J}}^t={\mathbf{j}}_0;{\mathbf{J}}^{t+1}={\mathbf{j}}_1;\dots ;{\mathbf{J}}^{t+s-1}={\mathbf{j}}_{s-1}\right)\hfill \\ \hfill & \times & \mathbb{P}\left({\mathbf{R}}^t={\mathbf{r}}_0\right)\mathbb{P}\left({\mathbf{R}}^{t+1}={\mathbf{r}}_1\right)\dots \mathbb{P}\left({\mathbf{R}}^{t+s-1}={\mathbf{r}}_{s-1}\right)p_t\left({\mathbf{n}}_0\right).\hfill \end{array}$$

In terms of the probability functions defined above we have,

$$\begin{array}{ccc}\hfill P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s)& =& \sum _{{\mathbf{j}}_0+{\mathbf{r}}_0={\mathbf{n}}_1}\sum _{{\mathbf{j}}_1+{\mathbf{r}}_1={\mathbf{n}}_2}\dots \sum _{{\mathbf{j}}_{s-1}+{\mathbf{r}}_{s-1}={\mathbf{n}}_s}f({\mathbf{j}}_0,{\mathbf{j}}_1,\dots ,{\mathbf{j}}_{s-1})h({\mathbf{r}}_0;{\mathbf{n}}_0)h({\mathbf{r}}_1;{\mathbf{n}}_1)\dots h({\mathbf{r}}_{s-1};{\mathbf{n}}_{s-1})p_t\left({\mathbf{n}}_0\right).\hfill \end{array}$$

The above expression states that the joint distribution $P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s)$ can be written in terms of the stationary distribution $p_{\mathrm{stat}}$ and the probability functions of the random vector ${\mathbf{R}}^t$ and the joint distribution of the random vectors ${\mathbf{J}}^t,{\mathbf{J}}^{t+1},\dots ,{\mathbf{J}}^{t+s-1}$, distributions that are given a priori. Once knowing the joint distribution $P({\mathbf{n}}_0,{\mathbf{n}}_1,\dots ,{\mathbf{n}}_s)$, we can compute the two-times joint distribution $P_{t,t+s}(\mathbf{n},\mathbf{m})$,

$$P_{t,t+s}(\mathbf{n},\mathbf{m})=\sum _{{\mathbf{n}}_1\in {\mathbb{N}}_0}\sum _{{\mathbf{n}}_2\in {\mathbb{N}}_0}\dots \sum _{{\mathbf{n}}_{s-1}\in {\mathbb{N}}_0}P(\mathbf{n},{\mathbf{n}}_1,\dots ,{\mathbf{n}}_{s-1},\mathbf{m}).$$

(106)

The above expression can be used to obtain the moment generating function of $({\mathbf{N}}^t,{\mathbf{N}}^{t+s})$ and then the corresponding two-times covariance matrix $C_{i,j}(t,t+s)={\left(\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})\right)}_{i,j}$. Those calculations are performed in Appendix A, here we only write down the result,

$$\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})={\mathbf{\Sigma }}_t{\mathbf{Q}}^s.$$

(107)

We should emphasize that Equation (107) is valid even if the system has not necessarily reached stationarity. If we assume that the system has attained the stationarity (which means that ${\mathbf{\Sigma }}_t$ no longer depend on time), we obtain

$$\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})=\overline{\mathbf{\Sigma }}{\mathbf{Q}}^s.$$

(108)

Due to stationarity, it is clear that the covariance matrix depends only on s, the difference between the times t and $t+s$.

Example 7.

Let us consider the system introduced in Example 5. We should notice that we were able to obtain an exact expression for the stationary variance matrix $\overline{\mathbf{\Sigma }}$. Thus, computing the covariance matrix involves only the product of two matrices. The resulting expression for the covariance matrix is too long to write down here. Thus, instead of giving explicitly the expression for the covariance matrix, we will show the theoretically computed time-dependent correlation functions defined as

$${\tilde{C}}_{i,j}\left(s\right):=\frac{{\left(\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})\right)}_{i,j}}{\sqrt{{\overline{\mathbf{\Sigma }}}_{i,i}{\overline{\mathbf{\Sigma }}}_{j,j}}}$$

(109)

In Figure 6 we show the behavior of the correlation functions ${\tilde{C}}_{1,1}\left(s\right)$ and ${\tilde{C}}_{1,2}\left(s\right)$ for the three-states open Markov chain studied in Example 5. For the parameter values we display the theoretically computed correlation function using Equations (108) and (109). The figure also shows the correlation functions obtained by means of numerical simulations of the system. The total time-steps performed to obtain the data from simulations was $5\times {10}^5$. We appreciate that the theoretical results agree with the simulations showing the consistency of our results.

5. Summary of Main Results

In this section we summarize the main results reported here. First of all, the proposed model for open Markov chain is described through the distribution of particles over the state space,

$$p_t\left(\mathbf{n}\right):=\mathbb{P}({\mathbf{N}}^t=\mathbf{n}),$$

(110)

which is ruled by the evolution equation

$$p_{t+1}\left(\mathbf{n}\right)=\sum _{\mathbf{k}\in {\mathbb{N}}_0^S}{p}_t\left(\mathbf{k}\right)K(\mathbf{k},\mathbf{n}),$$

(111)

where $K(\mathbf{k},\mathbf{n})=\mathbb{P}({\mathbf{J}}^t+{\mathbf{R}}^t=\mathbf{n})$. We proved that the above evolution equation can be formally solved using the moment generating function formalism. The time-dependent m.g.f. ${\mathcal{G}}_t\left(\mathit{\alpha }\right)=\mathbb{E}\left[e^{{\mathbf{N}}^t{\mathit{\alpha }}^T}\right]$ of $p_t\left(\mathbf{n}\right)$ can be written explicitly as

$$\begin{array}{ccc}\hfill {\mathcal{G}}_t\left(\mathit{\alpha }\right)& =& {\mathcal{G}}_0\left({\mathbf{H}}^{\left(t\right)}\left(\mathit{\alpha }\right)\right)\prod _{r=0}^{t-1}{\mathcal{F}}_{t-r}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right),\hfill \end{array}$$

(112)

where ${\mathcal{G}}_0$ stand for the m.g.f for ${\mathbf{N}}^0$ (the initial distribution over the state space) and ${\mathbf{H}}^{\left(r\right)}$ stands for the rth iterate of $\mathbf{H}$, i.e., ${\mathbf{H}}^{\left(r\right)}:=\mathbf{H}\circ \mathbf{H}\circ \dots \mathbf{H}$, r times, assuming stationarity. We also showed that, under some mild conditions, the system attains stationarity as t goes to infinity. In such a case the stationary m.g.f. ${\mathcal{G}}_{\mathrm{stat}}$ is given by,

$$\begin{array}{ccc}\hfill {\mathcal{G}}_{\mathrm{stat}}\left(\mathit{\alpha }\right)& =& \prod _{r=0}^{\infty }\mathcal{F}\left({\mathbf{H}}^{\left(r\right)}\left(\mathit{\alpha }\right)\right).\hfill \end{array}$$

(113)

We showed through a couple of examples that the above expression was useful to describe the stationary state as an infinite convolution of certain random variables. This approach allow, for instance, to compute the first cumulants. However we also developed a couple of formulas for the dynamics of the two first cumulants in terms of the variables of the internal dynamics and the protocols of outgoing particles. These dynamic equations are

$$\begin{array}{ccc}\hfill {\mathit{\mu }}_{t+1}& =& {\mathit{ϵ}}_t+{\mathit{\mu }}_t\mathbf{Q},\hfill \end{array}$$

(114)

$$\begin{array}{ccc}\hfill {\mathbf{\Sigma }}_{t+1}& =& {\mathbf{\Delta }}_t+{\mathbf{\Lambda }}_t+{\mathbf{Q}}^T{\mathbf{\Sigma }}_t\mathbf{Q},\hfill \end{array}$$

(115)

where ${\mathit{ϵ}}_t$ and $\mathbf{\Delta }$ are the mean vector and the variance matrix of the process $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$, i.e., ${\mathit{ϵ}}_t:=\mathbb{E}\left[{\mathbf{J}}^t\right]$ and ${\mathbf{\Delta }}_t:=\mathrm{Var}\left({\mathbf{J}}^t\right)$. Besides, ${\mathit{\mu }}_t$ and ${\mathbf{\Sigma }}_t$ correspond to the mean and the variance matrix of the number of particles in the state space respectively. The stationary version of Equations (114) and (115) can be written as

$$\begin{array}{ccc}\hfill \overline{\mathit{\mu }}& =& \mathit{ϵ}+\overline{\mathit{\mu }}\mathbf{Q},\hfill \end{array}$$

(116)

$$\begin{array}{ccc}\hfill \overline{\mathbf{\Sigma }}& =& \mathbf{\Delta }+\overline{\mathbf{\Lambda }}+{\mathbf{Q}}^T\overline{\mathbf{\Sigma }}\mathbf{Q},\hfill \end{array}$$

(117)

whose formal solution is

$$\begin{array}{ccc}\hfill \overline{\mathit{\mu }}& =& \mathit{ϵ}{(\mathbf{1}-\mathbf{Q})}^{-1},\hfill \end{array}$$

(118)

$$\begin{array}{ccc}\hfill \overline{\mathbf{\Sigma }}& =& \sum _{k=0}^{\infty }{\left({\mathbf{Q}}^T\right)}^k(\mathbf{\Delta }+\overline{\mathbf{\Lambda }}){\mathbf{Q}}^k.\hfill \end{array}$$

(119)

Clearly, to guarantee the existence of the stationary state it is necessary to assume that the process $\{{\mathbf{J}}^t:t\in \mathbb{N}\}$ be stationary, which implies that ${\mathit{ϵ}}_t$ and ${\mathbf{\Delta }}_t$ no longer depends on time. Following this formalism of moment generating functions we were also able to compute the dynamics of the m.g.f. ${\mathcal{R}}_t\left(\mathit{\alpha }\right)$ of the outgoing number of particles ${\mathbf{O}}^t$. Indeed we showed that

$${\mathcal{R}}_t\left(\mathit{\alpha }\right)={\mathcal{G}}_t\left(\mathbf{C}\left(\mathit{\alpha }\right)\right).$$

(120)

The above formula allowed us to prove that in the stationarity the number of outgoing particles equals the number of incoming particles, which can be interpreted as a kind of conservation law of the number of particles. Moreover, we were also able to obtain the behavior of the fluctuations of ${\mathbf{O}}^t$, i.e., we obtained the variance matrix of ${\mathbf{O}}^t$,

$$\begin{array}{ccc}\hfill \mathrm{Var}\left(O_t\right)& =& \mathbf{e}{\mathbf{\Sigma }}_t{\mathbf{e}}^T+{\mathit{\mu }}_t\left(\mathbf{1}-\mathbf{E}\right){\mathbf{e}}^T.\hfill \end{array}$$

(121)

Finally we also proved that the two-times covariance matrix $C_{i,j}\left(s\right)$ defined as

$$C_{i,j}(t,t+s):=\mathbb{E}\left[N_i^t{N}_j^{t+s}\right]-{\left({\mathit{\mu }}^t\right)}_i{\left({\mathit{\mu }}^{t+s}\right)}_j,$$

(122)

evolves according to

$$\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})={\mathbf{\Sigma }}_t{\mathbf{Q}}^s.$$

(123)

We should emphasize that Equation (123) is valid even if the system has not necessarily reached stationarity. If we assume that the system has attained the stationarity we obtain,

$$\mathrm{Cov}({\mathbf{N}}^t,{\mathbf{N}}^{t+s})=\overline{\mathbf{\Sigma }}{\mathbf{Q}}^s.$$

(124)

6. Conclusions

We have introduced a simple model for open Markov chains by interpreting the state space of a usual Markov chain as physical “sites” where non-interacting particles can be placed and moving throughout it according to “jumping rules” given by a kind of stochastic matrix. The conditions for the chain to be open are given as a protocol of incoming particles, defined by a discrete-time stochastic process, and by a protocol of outgoing particles, which is implicitly defined by the condition that the “stochastic matrix” (called here a jump matrix) has a spectral radius strictly less than one. These conditions establish the rules by means of which the particles arrive and leave the state space to the outside. We have shown that this model can be treated by means of the moment generating function technique, allowing us to obtain, in a closed form, the moment generating function of the distribution of particles over the state space. We have also shown that the system can be partially described by the dynamics of the two first cumulants of the distribution of particles over the state space. Actually, we have given closed formulas for the two first cumulants when the system is able to reach the stationarity. We have also studied how the correlations in the incoming protocol of particles are processed by the open chain. We have obtained closed formulas allowing to compute the two-times covariance matrix for the random vector defined as the number of particles on the states. Our main result is that the stationary two-times covariance matrix does not depend on the correlations of the particles arriving at the state space. This means that the stationary correlation functions essentially behaves as a closed Markov chain, i.e., that the correlations vanishes exponentially in time. The non-stationary correlations might probably content some information on the correlations of the incoming particles, but it would be necessary a more exhaustive study in this direction for a better understanding of such a process.