# From random point processes to hierarchical Cavity Master Equations for
stochastic dynamics of disordered systems in Random Graphs: Ising models and
epidemics

D. Machado  Group of Complex Systems and Statistical Physics. Department of
Theoretical Physics, Physics Faculty, University of Havana, Cuba
<EMAIL_ADDRESS>R. Mulet  Group of Complex Systems and Statistical
Physics. Department of Theoretical Physics, Physics Faculty, University of
Havana, Cuba<EMAIL_ADDRESS>
###### Abstract

We start from the Theory of Random Point Processes to derive n-point coupled
master equations describing the continuous dynamics of discrete variables in
random graphs. These equations constitute a hierarchical set of approximations
that generalize and improve the Cavity Master Equation (CME) previously
obtained in CME-PRE . Our derivation clarifies some of the hypotheses and
approximations that originally lead to the CME, considered now as the first
order of a more general technique. We tested the new scheme in the dynamics of
three models defined over diluted graphs: the Ising ferromagnet, the Viana-
Bray spin-glass and the susceptible-infectious-susceptible model for
epidemics. In the first two, the new equations perform similarly to the best-
known approaches in the literature. In the latter, they outperform the well-
known Pair Quenched Mean-Field Approximation.

##  I Introduction

Emissions from a radioactive source, time series of electrical energy in a
nerve fiber, the instants of arrival of customers in a queue, and the flips of
a spin due to thermal activation are all examples of random processes that
occur in such short intervals of time that can be classified as punctual.
Sometimes the distinction between point processes and other stochastic
processes is hard to define. For example, any stochastic process in continuous
time in which the sample paths are step functions is associated with a point
process, also any process with a discrete state space where a time of entry
into a state occurs randomly.

The Theory of Random Point Processes (TRPP) vanKampen ; Daley2002RPP  provides
the formal and practical background to study and understand these and similar
random collections of point occurrences. Usually, it is a matter of taste of
the researchers, or the convenience for the research to attain, or not, to the
TRPP to approach a specific problem.

In this work we study the continuous dynamics of a system of discrete
interacting variables. The model systems of interest can be viewed as a
multivariate and multidimensional Random Point Process. Alternatively, the
dynamics of each variable can be considered as a Random Point Process itself.
We will show that starting from the TRPP it is possible to write hierarchical
master-like equations for any group of variables in the system.

The first attempt in this direction was proposed a few years ago through the
adaptation of dynamic message-passing equations from discrete-time to
continuous-time dynamics. It landed first in the field of epidemics
propagation Karrer2010SIR ; shrestha2015message  without a general formulation
but with good results. Later, in CME-PRE  was proposed a new closure for the
master equation that exploited the TRPP. The core of this method is a
differential equation for the cavity conditional probability densities: the
Cavity Master Equation (CME). This approach has proven to be very general as
it has been applied to several models in graphs with finite connectivity like
the Ising ferromagnet, the random field Ising model and the Viana-Bray spin-
glass model CME-PRE , the ferromagnetic $p$-spin under Glauber dynamics CME-
Pspin  and more recently also the dynamics of a focused search algorithm to
solve the three-SAT problem in a random graph CME-PRL . In this work we
generalize the CME derived in CME-PRE , providing master equations for the
probability densities of any group of connected variables.

We concentrate on the continuous-time dynamics of discrete-spin variables.
They can be described by a master equation for the probability density of the
states of the system Glauber63 ; coolen2005theory . But to fully solve this
master equation is a cumbersome task and results have been elusive, except for
special cases. For example, the Sherrington-Kirkpatrick (SK) model in its non-
symmetric version has received preferential attention among the fully-
connected family Crisanti88 ; Rieger89 , and the exact solution for parallel
and asynchronous dynamics of the dilute fully-asymmetric neural network model
dates back to 1987 derrida1987exactly .

We provide here a general alternative solution that improves the approach
presented in CME-PRE . Although it is devised for single instances, it does
not prevents its use to obtain global information. One can always numerically
average the results of the integration of the single-instance equations. On
the other hand, in CME-AvAndBP  the CME equations were use as the starting
point in the derivation of an average case description for the dynamics of the
Ising ferromagnet, in this case defined on Erdos-Renyi random graphs.

To compute averaged quantities, one of the greatest advances in the field is
the Dynamical Replica Theory (DRT), first introduced in coolen1996DRT . This
approach permits the derivation of average case equations for the
probabilities of some macroscopic observables in fully-connected coolen1996DRT
; laughton1996order  and diluted mozeika2008dynamical ; mozeika2009dynamical
graphs. In practice DRT successfully reduces the dimensionality of the system,
making the problem tractable. However, it also assumes that the microscopic
probability distribution function is a constant within a subspace with a
finite number of order parameters, something that has been put into question
on specific models nishimori1996dynamical ; delferraro2014pert , and is not
necessarily true for general non-equilibrium situations.

Another relevant approach to these problems is the scheme of analytical
approximations introduced in BarthelwalkSAT2003 ; SemerjianwalkSAT2003  for
the description of the dynamics of local search algorithms. In
Semerjian_dyn_2004  this methodology is revisited and consolidated, showing
also its equivalence to DRT under the assumption of replica symmetry. Due to
its successful combination of accuracy and simplicity, we will come back to
this scheme as a reference for evaluating our numerical results.

This work contributes to the settlement of the dynamic cavity method as a
general tool to study the continuous-time dynamics of discrete-spin models. A
brief description of the theory of RPP is given in subsection II. Section III
presents our hierarchical system of corrections to the original equations
derived in CME-PRE  and used in CME-Pspin ; CME-PRL ; CME-AvAndBP . They are
numerically tested at section IV. We do not intend either to compare these
directly with DRT, whose numerics we consider appreciably harder in this kind
of graphs, nor to give a better description of the dynamics of average
quantities. Instead, we use the equivalent but simpler approximation schemes
in Semerjian_dyn_2004  as a theoretical and numeric reference, and as a proxy
to the comparison with DRT.

Unlike the two approaches previously described, our method gives information
about local probability distributions and can be straightforwardly used to
study problems that are intrinsically out of equilibrium CME-PRL , i.e., where
detailed balance is not conserved. As an example, the last subsection of IV is
devoted to the exploration of susceptible-infectious-susceptible model for
epidemics. Finally, section V contains the conclusions of our work.

##  II Random Point Processes

To facilitate the reading of section III we take a brief tour into the Theory
of Random Point Processes, which we will use to parametrize probability
distributions of spin histories. After some basic definitions and for
completeness we will refer in subsection II.1 to the simple case of an
independent binary variable that randomly changes its state. There, we will
use the random point processes formalism to derive a master equation for the
variable’s dynamics, and we will show how to extend this to the case of many
interacting variables.

The core object of random point processes are spin trajectories or histories.
A specific spin history that starts with $\sigma(t_{0})=\sigma_{0}$ and ends
with $\sigma(t)=\sigma$ is parametrized by the number of spin flips, the time
in which they occur and the initial state of the system. The spin trajectory
is then nothing but a random point process vanKampen ; Daley2002RPP  and the
probability measure in this sample space may be denoted as:

$Q^{t}(X)=Q_{s}(t_{0},t_{1},...,t_{s},t\mid\sigma(t_{0})=\sigma_{0})$ (1)

which represents the probability density of having a trajectory with $s$ jumps
at $(t_{1},t_{1}+dt_{1}),...,(t_{s},t_{s}+dt_{s})$. When needed, we may write
$X$ as $X(t)$ to emphasize that the final time of the spin history is
precisely $t$. Seeking simplicity, in almost all equations we will not write
the conditioning of all dynamic quantities on a given initial condition (for
example, in (1) we have $\sigma(t_{0})=\sigma_{0}$).

The probability density $Q(X)$ fulfills the following normalization relation:

$1=\sum_{X}^{t}Q^{t}(X)=\sum_{s=0}^{\infty}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\,...\,\int_{t_{s-1}}^{t}dt_{s}\,Q_{s}(t_{0},t_{1},...,t_{s},t\mid\sigma(t_{0})=\sigma_{0})$
(2)

The sum $\sum_{X}^{t}$ in equation (2) goes over all histories that occur
between $t_{0}$ and $t$, starting with $\sigma(t_{0})=\sigma_{0}$. It is
explicitly written in the right hand side as a sum $\sum_{s=0}^{\infty}$ of
time integrals. Each term of that sum corresponds to histories that have
exactly $s$ jumps, and the integrals account for all the possible jumping
times $t_{1},t_{2},...,t_{s}$.

Instantaneous magnitudes can be written as marginals of $Q^{t}(X)$. For
example:

$P^{t}(\sigma)=\sum_{X\mid\sigma(t)=\sigma}^{t}Q^{t}(X)$ (3)

where the sum $\sum_{X\mid\sigma(t)=\sigma}^{t}$ goes over all the histories
$X(t)$ such that $\sigma(t)=\sigma$ and can be written similarly as
$\sum_{X}^{t}$ (see (2)). In what follows we will use the symbol
$\sum_{X\mid\sigma}^{t}$ as a shortening of
$\sum_{X\mid\sigma(t)=\sigma}^{t}$.

As we will show, some properties of the time evolution of magnitudes like
$P^{t}(\sigma)$ can be obtained by working directly over $Q^{t}(X)$.

###  II.1 From random point processes to master equations

Master equations are a common instrument in the study of continuous-time
dynamics. As a proof of the validity of the random point process approach, we
will show its equivalence with the formulation of master equations in two
examples: the time evolution of a single non-interacting variable and the
dynamics of N-interacting spins. Through the derivations, the reader will
become familiar with properties and procedures that will be used in section
III.

####  II.1.1 Single variable

The master equation for the time evolution a single spin without interactions
reads:

$\dfrac{dP^{t}(\sigma)}{dt}=-r(\sigma)P^{t}(\sigma)+r(-\sigma)P^{t}(\sigma)$
(4)

where $r(\sigma)$ is the transition rate between the states $\sigma$ and
$-\sigma$.

In order to derive equation (15) from the TRPP we will differentiate (3).
Differentiation in this context should be handled carefully since increasing
$t$ means we are changing the sample space itself. Instead of using calculus
rules to differentiate we use the definition of derivative as the limit of the
incremental ratio:

$\frac{dP^{t}(\sigma)}{dt}=\lim_{\Delta t\rightarrow 0}\dfrac{P^{t+\Delta
t}(\sigma)-P^{t}(\sigma)}{\Delta t}$ (5)

For the sake of the discussion in the next sections let’s write down the
calculations in detail.

Let’s say that the final state of the trajectory $X$ is $X(t)=\sigma$. Then
the probability density $Q^{t}(X)$ can take the form vanKampen ; Daley2002RPP
:

$Q^{t}(X)=r(\sigma_{0})e^{-r(\sigma_{0})(t_{1}-t_{0})}r(-\sigma_{0})e^{-r(-\sigma_{0})(t_{2}-t_{1})}\ldots
r(-\sigma)e^{-r(-\sigma)(t_{s}-t_{s-1})}e^{-r(\sigma)(t_{f}-t_{s})}$ (6)

The right hand side of equation (6) is the probability density of the waiting
times of the history $X:[t_{0},t]\rightarrow\\{1,-1\\}$, whose jumps occur at
$(t_{1},t_{2},...\,,t_{s})$, with rates $r(\sigma)$.

Now we will write $P^{t+\Delta t}(\sigma)$ as the marginalization of
$Q^{t+\Delta t}(X)$ and then expand to first order in $\Delta t$. Explicitly
writing that marginalization as in (2):

$\displaystyle P^{t+\Delta t}(\sigma)=\sum_{s}\int_{t_{0}}^{t+\Delta
t}dt_{1}\int_{t_{1}}^{t+\Delta t}dt_{2}\ldots\int_{t_{s-1}}^{t+\Delta
t}dt_{s}\,Q^{t+\Delta t}(X)$ (7)

Here, the sum on the right-hand side goes over all the values of the number of
jumps (denoted as $s$) such that $X(t)=\sigma$.

Now we must keep only the order $\Delta t$ term in (7). Let’s expand first
$Q^{t+\Delta t}(X)$ and leave the rest untouched. The equation (6) can be
rearranged as:

$Q^{t+\Delta
t}(X)=\Big{[}\prod_{l=1}^{s}r(\sigma(t_{l}))\Big{]}\exp\Big{\\{}-\int_{t_{0}}^{t+\Delta
t}r(\sigma(\tau))d\tau\Big{\\}}$ (8)

where the jumps in the history $X$ occur at the times $t_{l}$, with
$l=1,2,\ldots,s$.

In (8), the probability of having the last jump in the interval $[t,t+\Delta
t]$ is $r(-\sigma)\Delta t$ and, more generally, the probability of occurrence
of $n$ jumps is proportional to $(\Delta t)^{n}$. Therefore, when $\Delta t$
goes to zero, with probability $1$ there are no jumps in $[t,t+\Delta t]$ and
none of the times $t_{l}$ belongs to that interval.

We can then expand the exponential in (8) and write:

$\displaystyle Q^{t+\Delta t}(X)$ $\displaystyle=$
$\displaystyle\Big{[}\prod_{l=1}^{s}r(\sigma(t_{l}))\Big{]}\exp\Big{\\{}-\int_{t_{0}}^{t}r(\sigma(\tau))d\tau\Big{\\}}\big{[}1-r(\sigma)\Delta
t\big{]}+o(\Delta t)$ $\displaystyle Q^{t+\Delta t}(X)$ $\displaystyle=$
$\displaystyle Q^{t}(X)\big{[}1-r(\sigma)\Delta t\big{]}+o(\Delta t)$ (9)

On the other hand, the sum of iterated integrals in (7) is of order $(\Delta
t)^{s}$ by itself. However, for each integral we have the property
$\int_{t_{0}}^{t+\Delta t}d\tau=\int_{t_{0}}^{t}d\tau+\int_{t}^{t+\Delta
t}d\tau$, and we can write:

$\displaystyle\sum_{s}\int_{t_{0}}^{t+\Delta t}dt_{1}\int_{t_{1}}^{t+\Delta
t}dt_{2}\ldots\int_{t_{s-1}}^{t+\Delta t}dt_{s}$ $\displaystyle=$
$\displaystyle\sum_{s}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\ldots\int_{t_{s-1}}^{t}dt_{s}\;+$
(10)
$\displaystyle+\sum_{s}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\ldots\int_{t}^{t+\Delta
t}dt_{s}+o(\Delta t)$

The first term of the right hand side of (10) is an operator over the space of
the histories that occur in the interval $[t_{0},t]$. It can be safely applied
to the expansion (9) of the probability density $Q^{t+\Delta t}(X)$ to give a
contribution of order $O((\Delta t)^{0})\equiv O(1)$

$I_{0}=\sum_{s}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\ldots\int_{t_{s-1}}^{t}dt_{s}Q^{t}(X)=P^{t}(\sigma)$
(11)

and a contribution of order $O(\Delta t)$

$I_{1}=-\sum_{s}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\ldots\int_{t_{s-1}}^{t}dt_{s}Q^{t}(X)\,r(\sigma)\Delta
t=-P^{t}(\sigma)\,r(\sigma)\Delta t$ (12)

Together, $I_{0}+I_{1}$ represent the probability density of having
$\sigma(t)=\sigma$ and no jumps in $[t,t+\Delta t]$.

In the second sum of (10) we have an operator that acts over the space where
all jumps, except the last one, took place in $[t_{0},t]$, and the last jump
occurs in $[t,t+\Delta t]$. Therefore, we can obtain a second contribution of
order $O(\Delta t)$ by applying this operator to the probability density $Q$
of having only one jump in $[t,t+\Delta t]$. Remembering the parametrization
(8) and expanding in powers of $\Delta t$:

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\sum_{s}\int_{t_{0}}^{t}dt_{1}\ldots\int_{t_{s-2}}^{t}dt_{s-1}\Big{[}\prod_{l=1}^{s-1}r(\sigma(t_{l}))\Big{]}e^{-\int_{t_{0}}^{t}r(\sigma(\tau))d\tau}\int_{t}^{t+\Delta
t}dt_{s}\,r(\sigma(t_{s}))e^{-\int_{t}^{t+\Delta t}r(\sigma(\tau))d\tau}$
$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\Big{(}\sum_{s}\int_{t_{0}}^{t}dt_{1}\ldots\int_{t_{s-2}}^{t}dt_{s-1}\Big{[}\prod_{l=1}^{s-1}r(\sigma(t_{l}))\Big{]}e^{-\int_{t_{0}}^{t}r(\sigma(\tau))d\tau}\Big{)}\Big{(}r(-\sigma)\Delta
t\Big{)}+o(\Delta t)$ (13)

In (13) we have used that, in order to fulfill the condition
$\sigma(t)=\sigma$, the last jump must happen from $-\sigma$ to $\sigma$, with
rate $r(-\sigma)$. The sum inside parenthesis in the equation now goes over
all histories that end with $\sigma(t)=-\sigma$, and therefore is equal to
$P^{t}(-\sigma)$. Thus $I_{2}=P^{t}(-\sigma)\,r(-\sigma)\Delta t+o(\Delta t)$
represents the probability density of having $\sigma(t)=-\sigma$ and only one
jump in $[t,t+\Delta t]$.

Putting all this together:

$\displaystyle P^{t+\Delta t}(\sigma)$ $\displaystyle=$ $\displaystyle
I_{0}+I_{1}+I_{2}+o(\Delta t)$ $\displaystyle P^{t+\Delta t}(\sigma)$
$\displaystyle=$ $\displaystyle P^{t}(\sigma)-P^{t}(\sigma)\,r(\sigma)\Delta
t+P^{t}(-\sigma)\,r(-\sigma)\Delta t+o(\Delta t)$ (14)

By subtracting $P^{t}(\sigma)$ from both sides of (14), dividing by $\Delta t$
and taking the limit $\Delta t\rightarrow 0$ we obtain the desired result (4):

$\dfrac{dP^{t}(\sigma)}{dt}=-r(\sigma)P^{t}(\sigma)+r(-\sigma)P^{t}(-\sigma)$

####  II.1.2 Multiple variables

In this case, where we have the time evolution of the N-interacting spins
$\vec{\sigma}=\\{\sigma_{1},...,\sigma_{N}\\}$, the master equation is:

$\dfrac{dP^{t}(\vec{\sigma})}{dt}=-\sum_{i=1}^{N}\big{[}r_{i}(\vec{\sigma})P^{t}(\vec{\sigma})+r_{i}(F_{i}[\vec{\sigma}])P^{t}(F_{i}[\vec{\sigma}])\big{]}$
(15)

where $F_{i}[\vec{\sigma}]$ is an operator that transforms the state
$\vec{\sigma}=\\{\sigma_{1},...,\sigma_{i},...,\sigma_{N}\\}$ into the state
$F_{i}[\vec{\sigma}]=\\{\sigma_{1},...,-\sigma_{i},...,\sigma_{N}\\}$, and
$r_{i}(\vec{\sigma})$ is the transition rate between $\vec{\sigma}$ and
$F_{i}[\vec{\sigma}]$.

Similarly as in (3), the probability density $P^{t}(\vec{\sigma})$ is given
by:

$P^{t}(\vec{\sigma})=\sum_{\vec{X}\mid\vec{\sigma}}^{t}Q^{t}(\vec{X})$ (16)

where $\vec{X}=\\{X_{1},\ldots,X_{N}\\}$ is the vector of the histories of all
spins.

Now we must differentiate (16) similarly as with (3):

$\frac{dP^{t}(\vec{\sigma})}{dt}=\lim_{\Delta t\rightarrow
0}\dfrac{P^{t+\Delta t}(\vec{\sigma})-P^{t}(\vec{\sigma})}{\Delta t}$ (17)

The set of individual histories of $N$-interacting spins, $Q^{t}(\vec{X})$,
can be written as a product of probability densities
$\Phi_{a}^{t}(X_{a}\mid\vec{X}_{\setminus a})$ of the history $X_{a}$ with the
histories of all the other spins $\vec{X}_{\setminus a}$ fixed kipnis99 :

$Q^{t}(X_{1},X_{2},...,X_{N})=\prod_{i=1}^{N}\Phi_{i}^{t}(X_{i}\mid\vec{X}_{\setminus
i})$ (18)

Each $\Phi_{i}^{t}(X_{i}\mid\vec{X}_{\setminus i})$ can be parametrized like
in (8) vanKampen ; Daley2002RPP :

$\Phi_{i}^{t}(X_{i}|\vec{X}_{\setminus
i})=\Big{[}\prod_{l_{i}=1}^{s_{i}}r_{i}(\vec{\sigma}(t_{l_{i}}))\Big{]}\exp\Big{\\{}-\int_{t_{0}}^{t}r_{i}(\vec{\sigma}(\tau))d\tau\Big{\\}}$
(19)

Notice that when any $X_{j}$ with $j\in\partial i$ has a jump, the value of
$r_{i}(\vec{\sigma}(\tau))$ changes in the integral inside the exponential
function. As before, we will write $P^{t+\Delta t}(\vec{\sigma})$ as the
marginalization of $Q^{t+\Delta t}(\vec{X})$ and then expand to first order in
$\Delta t$. Explicitly writing that marginalization as in (7):

$\displaystyle P^{t+\Delta
t}(\vec{\sigma})=\sum_{s_{1}=0}^{\infty}\sum_{s_{2}=0}^{\infty}\\!\\!\ldots\\!\\!\sum_{s_{N}=0}^{\infty}\Big{[}\prod_{i=1}^{N}\int_{t_{0}}^{t+\Delta
t}\\!\\!\\!\\!\\!dt_{1}^{i}\int_{t_{1}^{i}}^{t+\Delta
t}\\!\\!\\!\\!\\!dt_{2}^{i}\ldots\int_{t_{s-1}^{i}}^{t+\Delta
t}\\!\\!\\!\\!\\!dt_{s}^{i}\Big{]}Q^{t+\Delta t}(X_{1},X_{2},\ldots,X_{N})$
(20)

To expand $Q^{t+\Delta t}(\vec{X})$ requires expanding each
$\Phi_{i}^{t+\Delta t}(X_{i}\mid\vec{X}_{\setminus i})$, which is almost the
same we did in (9):

$\Phi_{i}^{t+\Delta t}\big{(}X_{i}\mid\vec{X}_{\setminus
i}\big{)}=\Phi_{i}^{t}\big{(}X_{i}\mid\vec{X}_{\setminus
i}\big{)}\,\big{[}1-r_{i}(\vec{\sigma}(t))\,\Delta t\big{]}+o(\Delta t)$ (21)

Thus:

$\displaystyle Q^{t+\Delta t}(\vec{X})$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\Phi_{i}^{t+\Delta
t}\big{(}X_{i}\mid\vec{X}_{\setminus i}\big{)}$ $\displaystyle Q^{t+\Delta
t}(\vec{X})$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\Phi_{i}^{t}\big{(}X_{i}\mid\vec{X}_{\setminus
i}\big{)}\Big{[}1-\Delta t\sum_{k=1}^{N}r_{k}(\vec{\sigma}(t))\Big{]}+o(\Delta
t)$ $\displaystyle Q^{t+\Delta t}(\vec{X})$ $\displaystyle=$ $\displaystyle
Q^{t}(\vec{X})\Big{[}1-\Delta
t\sum_{k=1}^{N}r_{k}(\vec{\sigma}(t))\Big{]}+o(\Delta t)$ (22)

Similarly as with (11) and (12) we will use (22) to get two contributions of
order $O(1)$ and $O(\Delta t)$, respectively:

$\displaystyle I_{0}$ $\displaystyle=$
$\displaystyle\sum_{s_{1}=0}^{\infty}\\!\ldots\\!\sum_{s_{N}=0}^{\infty}\Big{[}\prod_{i=1}^{N}\int_{t_{0}}^{t}\\!\\!\\!dt_{1}^{i}\ldots\int_{t_{s-1}^{i}}^{t}\\!\\!\\!dt_{s}^{i}\Big{]}Q^{t}(\vec{X})=P^{t}(\vec{\sigma})$
(23) $\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\Big{(}\sum_{s_{1}=0}^{\infty}\\!\ldots\\!\sum_{s_{N}=0}^{\infty}\Big{[}\prod_{i=1}^{N}\int_{t_{0}}^{t}\\!\\!\\!dt_{1}^{i}\ldots\int_{t_{s-1}^{i}}^{t}\\!\\!\\!dt_{s}^{i}\Big{]}Q^{t}(\vec{X})\Big{)}\Big{(}-\Delta
t\sum_{k=1}^{N}r_{k}(\vec{\sigma}(t))\Big{)}$ $\displaystyle I_{1}$
$\displaystyle=$ $\displaystyle-P^{t}(\vec{\sigma})\,\Delta
t\,\sum_{k=1}^{N}r_{k}(\vec{\sigma}(t))$ (24)

As in (10) we can split the sum of iterated integrals in the right hand side
of (20) to get $O(\Delta t)$ contributions. The latter is an operator that
acts over the space of the histories where one and only one spin jumps in
$[t,t+\Delta t]$. he rest is a calculation essentially equal to the one in
(13):

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N}\,\sum_{s_{1}=0}^{\infty}\\!\ldots\\!\sum_{s_{N}=0}^{\infty}\Big{[}\prod_{i\neq
k}\int_{t_{0}}^{t}\\!\\!\\!dt_{1}^{i}\ldots\int_{t_{s-1}^{i}}^{t}\\!\\!\\!dt_{s}^{i}\Big{]}\int_{t_{0}}^{t}\\!\\!\\!dt_{1}^{k}\ldots\int_{t_{s-2}}^{t}\\!\\!\\!dt_{s-1}^{k}\Big{[}\prod_{i\neq
k}\Phi_{i}^{t}(X_{i}\mid\vec{X}_{\setminus i})\Big{]}\times$
$\displaystyle\times\Big{[}\prod_{l_{k}=1}^{s_{k}-1}r_{k}(\vec{\sigma}(t_{l_{k}}))\Big{]}e^{-\int_{t_{0}}^{t}r_{k}(\vec{\sigma}(\tau))d\tau}r_{k}\big{(}F_{k}[\vec{\sigma}(t)]\big{)}\Delta
t+o(\Delta t)$ $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N}P^{t}(F_{k}[\vec{\sigma}])\,r_{k}\big{(}F_{k}[\vec{\sigma}(t)]\big{)}\,\Delta
t$ (25)

Finally, the expansion of $P^{t+\Delta t}(\vec{\sigma})$ gives:

$\displaystyle P^{t+\Delta t}(\vec{\sigma})$ $\displaystyle=$ $\displaystyle
I_{0}+I_{1}+I_{2}+o(\Delta t)$ $\displaystyle P^{t+\Delta t}(\vec{\sigma})$
$\displaystyle=$ $\displaystyle P^{t}(\vec{\sigma})-\Delta
t\sum_{k=1}^{N}P^{t}(\vec{\sigma})r_{k}(\vec{\sigma})+\Delta
t\sum_{k=1}^{N}P^{t}(F_{k}[\vec{\sigma}])r_{k}(F_{k}[\vec{\sigma}])+o(\Delta
t)$ (26)

and using (17) we obtain the usual Master Equation vanKampen  for a set of
N-interacting spins (15):

$\dfrac{dP^{t}(\vec{\sigma})}{dt}=-\sum_{i=1}^{N}\big{[}r_{i}(\vec{\sigma})P^{t}(\vec{\sigma})+r_{i}(F_{i}[\vec{\sigma}])P^{t}(F_{i}[\vec{\sigma}])\big{]}$
(27)

Unfortunately, to solve (27), even numerically, is in general a very difficult
task because the densities $P^{t}(\vec{\sigma})$ are high dimensional objects.

##  III Hierarchical Cavity Master Equations

This section contains the main analytical contribution of our work. We will
exploit techniques similar to the ones presented above to write down a set of
closed differential equations for the stochastic dynamics of discrete
variables in a Random Graph. We generalize a closure presented in CME-PRE
substituting an uncontrolled approximation by new equations derived from first
principles through the Theory of Random Point Processes.

To simplify the reading this section is divided in five subsections. At III.1
we introduce the dynamic cavity method to study systems defined over tree-like
graphs. In subsection III.2 we use that formulation to write the local
probability densities sitting on any group of connected nodes in terms of
dynamic cavity messages. Then we derive the known local master equations for
those probability densities. In subsection III.3 we repeat the same procedure
to obtain analogous equations for cavity probability densities. A general
approximated method for closing these equations is presented and discussed in
subsection III.4. There, we show how these closed dynamic equations can be
organized in a system of hierarchical approximations.

###  III.1 Treelike architecture and cavity messages

In this subsection we will introduce the dynamic cavity messages for tree-like
graphs, where it is possible to separate the spins into two disconnected
networks just by removing a single connection between two nodes.

We first select a spin, say $i$, and rewrite equation (18) expanding the tree
around it and making use of its structure:

$Q^{t}(X_{1},\ldots,X_{N})=\Phi_{i}^{t}(X_{i}|X_{\partial
i})\prod_{k\in\partial i}\Big{[}\Phi_{k}^{t}(X_{k}|X_{\partial
k})\prod_{m\in\partial k\setminus i}\Big{(}\Phi_{k}^{t}(X_{m}|X_{\partial
m})\prod_{l\in\partial m\setminus k}\ldots\Big{)}\Big{]}$ (28)

In equation (28) we used the symbol $\partial i$ for the set of spins in the
neighborhood of $i$. Those are the only spins that directly interact with $i$
and therefore we adapted the previous notation of the $\Phi$ probability
densities. The symbol $X_{\partial i}$ represents the set of histories of the
nodes in $\partial i$. Furthermore, let $G_{k}^{(i)}$ be the sub-graph
expanded from the site $k$ after removing the link $(ik)$. We define as
$\\{X\\}_{ik}$ the set of histories of the spins included in $G_{k}^{(i)}$
except $X_{k}$ itself. With these definitions we express (28) as:

$Q^{t}(X_{1},\ldots,X_{N})=\Phi_{i}^{t}(X_{i}|X_{\partial
i})\prod_{k\in\partial i}M_{ki}^{t}(X_{i},X_{k},\\{X\\}_{ik})$ (29)

Here $M_{ki}^{t}$ is just a shortening for the expression inside brackets.
Marginalizing $Q$ on all histories except $X_{i},X_{\partial i}$ we get:

$Q^{t}(X_{i},X_{\partial i})=\Phi_{i}^{t}(X_{i}|X_{\partial
i})\prod_{k\in\partial i}\mu_{k\rightarrow(ki)}^{t}(X_{k}|X_{i})$ (30)

The new functions $\mu_{k\rightarrow(ki)}^{t}(X_{k}|X_{i})$, called dynamic
cavity messages, are the marginals:

$\mu_{k\rightarrow(ik)}^{t}(X_{k}|X_{i})=\sum_{\\{X\\}_{ik}}^{t}M_{ki}^{t}(X_{i},X_{k},\\{X\\}_{ik})$
(31)

and have the interpretation of the probability density of history $X_{k}$
given $X_{i}$ fixed.

If we take now two neighbors, $i$ and $j$, and make a similar reasoning we
conclude that their marginal probability density can be written as:

$Q^{t}(X_{i},X_{j})=\mu_{i\rightarrow(ij)}^{t}(X_{i}|X_{j})\mu_{j\rightarrow(ji)}^{t}(X_{j}|X_{i})$
(32)

These cavity messages can be parametrized as other similar dynamical
quantities (see equations (6), (19)):

$\displaystyle\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid
X_{j})=\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t_{0})\,e^{-\int_{t_{0}}^{t_{1}}\lambda_{i}^{\tau}d\tau}\times\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t_{1})e^{-\int_{t_{1}}^{t_{2}}\lambda_{i}^{\tau}d\tau}\times...$
(33)
$\displaystyle...\times\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t_{s_{i}})e^{-\int_{t_{s_{i}}}^{t}\lambda_{i}^{\tau}d\tau}$

only that now the jumps occur with unknown rates
$\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t)$.

Then, as in (21) we can expand $\Phi$ and $\mu$ to order $\Delta t$ (for the
latter, see Appendix A):

$\displaystyle\Phi_{i}^{t+\Delta t}\big{(}X_{i}\mid X_{\partial i}\big{)}$
$\displaystyle=$ $\displaystyle\Phi_{i}^{t}\big{(}X_{i}\mid X_{\partial
i}\big{)}\,\big{[}1-r_{i}(\sigma_{i}(t),\sigma_{\partial i}(t))\,\Delta
t\big{]}+o(\Delta t)$ (34) $\displaystyle\mu_{i\rightarrow(ij)}^{t+\Delta
t}(X_{i}|X_{j})$ $\displaystyle=$
$\displaystyle\mu_{i\rightarrow(ij)}^{t}(X_{i}|X_{j})\,[1-\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t)\Delta
t]+o(\Delta t)$ (35)

The equations (31)-(35) were already present in CME-PRE , but they will be
also at the basis of our derivation. They will be exploited below to write
down a more general set of equations than the ones presented in CME-PRE .

###  III.2 Equations for local probability densities

By eliminating a group of connected nodes from a tree-like graph, we always
divide it into several tree-like sub-graphs. We illustrate this in the top-
left panel of figure (1). There, we colored in gray the sub-trees obtained
after removing the white nodes in the center. As we learned in subsection
III.1, if we marginalize the full distribution $Q^{t}(\vec{X})$ over all the
histories in these gray sub-graphs, we obtain the probability density of the
set of histories corresponding to the connected group we have selected.

Figure 1: (Top-left panel): Selection of a connected set with $4$ nodes
(white) from a tree-like graph. The full graph is divided in $6$ tree-like
sub-graphs (gray nodes). (Top-right panel): Dividing a connected set with $8$
nodes in two groups: inner (white) and outer nodes (black). This is a
connected set of the first kind, because every outer node has only one
neighbor inside the set. (Bottom-left panel): Connected set with 7 nodes.
According to our classification the connected set is of the second kind
because the node $l$ has two neighbors, $i$ and $l_{2}$, which belong to the
set. (Bottom-right panel): Local weights and cavity messages present in the
probability densities related to the connected set shown in top-right panel.

Therefore, if we denote the central node in the top-left panel of figure (1)
by the letter $i$, and remembering the steps from (28) to (30), we have:

$\displaystyle Q^{t}(\vec{X}_{\text{connected set}})$ $\displaystyle=$
$\displaystyle\sum_{\vec{X}_{\text{gray nodes}}}^{t}Q^{t}(\vec{X})$
$\displaystyle Q^{t}(X_{i},X_{\partial i})$ $\displaystyle=$
$\displaystyle\Phi_{i}^{t}(X_{i}|X_{\partial i})\prod_{k\in\partial
i}\mu_{k\rightarrow(ki)}^{t}(X_{k}|X_{i})$ (36)

This procedure can be easily generalized to write any local probability
density of the histories of a connected set of nodes. We just have to
recognize inner and outer nodes: the first ones have all their neighbors
belonging to the connected set, and the second ones have at least one neighbor
outside the connected set. The top-right and the bottom-left panels of figure
(1) show examples of connected sets. In each case, we have a group of inner
nodes (in white), and a group of outer nodes (in black).

The probability density of a connected set is always a marginal of the full
probability density. However, in order to avoid a more involved analysis, we
will separate the connected sets into two categories. A connected set is _of
the first kind_ when all its outer nodes have only one neighbor inside the
set, and is _of the second kind_ when at least one outer node has more than
one neighbor inside the set. This distinction is important because in the
first case we can directly write the corresponding probability density as a
product of local weights $\Phi$ and dynamic cavity messages $\mu$. We must
simply include a local weight for each inner node and a cavity message for
each outer node.

For example, the probability density corresponding to the connected set in the
top-right panel of figure (1) is:

$\displaystyle Q^{t}(\vec{X}_{\text{connected set}})$ $\displaystyle=$
$\displaystyle\sum_{\vec{X}_{\text{gray nodes}}}^{t}Q^{t}(\vec{X})$
$\displaystyle Q^{t}(X_{i},X_{\partial i},X_{\partial k\setminus
i},X_{\partial l\setminus i})$ $\displaystyle=$
$\displaystyle\Phi_{i}^{t}(X_{i}|X_{\partial i})\Big{[}\prod_{m\in\partial
i\setminus j}\Phi_{m}^{t}(X_{m}|X_{\partial
m})\Big{]}\mu_{j\rightarrow(ji)}^{t}(X_{j}|X_{i})\times$ (37)
$\displaystyle\>\>\>\>\>\times\Big{[}\prod_{n\in\partial k\setminus
i}\mu_{n\rightarrow(nk)}^{t}(X_{n}|X_{k})\Big{]}\Big{[}\prod_{n\in\partial
l\setminus i}\mu_{n\rightarrow(nl)}^{t}(X_{n}|X_{l})\Big{]}$

where $j,k,l\in\partial i$. A representation of this product of $\Phi$ and
$\mu$ functions is shown in the bottom-right panel of (1).

Connected sets of the second kind do not require a different treatment.
Fortunately, we can always add nodes to convert these into sets of the first
kind. For example, we can transform the bottom-left panel of figure (1) into
the top-right panel just by adding the node $l_{1}$ to the connected set. Then
we can simply compute:

$Q^{t}(X_{i},X_{\partial i},X_{\partial k\setminus
i},X_{l_{2}})=\sum_{X_{l_{1}}}^{t}Q^{t}(X_{i},X_{\partial i},X_{\partial
k\setminus i},X_{l_{1}},X_{l_{2}})$ (38)

Thus, the probability density of any set of the second kind can be written as
a marginal of the probability density of a set of the first kind just like
(38) is written in terms of (37).

In order to derive a differential equation for the time evolution of any
instantaneous local probability density $P^{t}$, we should concentrate on
connected sets of the first kind. To fix ideas, let’s work with such a
connected set containing $n$ inner nodes and $m$ outer nodes. As each outer
node has one neighbor which is an inner node, we will have $n\geq m$.

Let’s construct two subsets of nodes: $\mathcal{O}$ and $\mathcal{I}$. First,
we fill $\mathcal{O}$ with all outer nodes in an arbitrary order:
$\mathcal{O}=\\{o_{1},o_{2},\ldots,o_{m}\\}$. Then we build the set
$\mathcal{I}=\\{i_{1},i_{2},\ldots,i_{m},i_{m+1},\ldots,i_{n}\\}$ in such a
way that each node $i_{k}\in\mathcal{I}$ is a neighbor of the node
$o_{k}\in\mathcal{O}$, with $k\leq m$.

The instantaneous probability density can be written as:

$P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})=\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$
(39)

where $\vec{\sigma}_{\mathcal{I}}$ and $\vec{\sigma}_{\mathcal{O}}$ are the
vectors of the instantaneous states of the inner and outer nodes,
respectively, and $\vec{X}_{\mathcal{I}}$ and $\vec{X}_{\mathcal{O}}$ are the
corresponding vectors of the histories of those nodes. On the other hand:

$Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})=\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow(o_{k}i_{k})}^{t}(X_{o_{k}}|X_{i_{k}})\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}(X_{i_{k}},X_{\partial
i_{k}})\Big{]}$ (40)

Now we need to expand (39) to order $\Delta t$. As in subsection II.1 we can
write an order $\Delta t$ expression for $Q^{t+\Delta
t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$. Remembering (34) and (35) we
have:

$\displaystyle Q^{t+\Delta t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$
$\displaystyle=$ $\displaystyle\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\,[1-\lambda_{o_{k}\rightarrow i_{k}}^{t}\Delta t]+o(\Delta
t)\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\,\big{[}1-r_{i_{k}}^{t}\,\Delta
t\big{]}+o(\Delta t)\Big{]}$ $\displaystyle Q^{t+\Delta
t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$ $\displaystyle=$
$\displaystyle Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})-\Delta
t\sum_{l=1}^{n}r_{i_{l}}^{t}\;\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}-$ (41)
$\displaystyle-\Delta t\sum_{l=1}^{m}\lambda_{o_{l}\rightarrow
i_{l}}^{t}\;\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}$

where we have shortened our notation as follows:

$\displaystyle\Phi_{i_{k}}^{t}(X_{i_{k}}|X_{\partial i_{k}})$
$\displaystyle\equiv$ $\displaystyle\Phi_{i_{k}}^{t}$
$\displaystyle\mu_{o_{k}\rightarrow(o_{k}i_{k})}^{t}(X_{o_{k}}|X_{i_{k}})$
$\displaystyle\equiv$ $\displaystyle\mu_{o_{k}\rightarrow i_{k}}^{t}$ (42)
$\displaystyle r_{i_{k}}(\sigma_{i_{k}}(t),\sigma_{\partial i_{k}}(t))$
$\displaystyle\equiv$ $\displaystyle r_{i_{k}}^{t}$
$\displaystyle\lambda_{o_{k}\rightarrow(o_{k}i_{k})}(X_{o_{k}},X_{i_{k}},t)$
$\displaystyle\equiv$ $\displaystyle\lambda_{o_{k}\rightarrow i_{k}}^{t}$

Similarly as before, there are three contributions to the expansion of (39) to
order $\Delta t$:

$\displaystyle I_{0}$ $\displaystyle=$
$\displaystyle\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})=P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$
(43) $\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle-\Delta
t\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\Big{[}\sum_{l=1}^{n}r_{i_{l}}^{t}+\sum_{l=1}^{m}\lambda_{o_{l}\rightarrow
i_{l}}^{t}\Big{]}$ (44) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\Delta
t\sum_{l=1}^{n}\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}]}^{t}\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>r_{i_{l}}^{t}+$
(45) $\displaystyle+\Delta
t\sum_{l=1}^{m}\sum_{\vec{X}_{\mathcal{O}}|F_{o_{l}}[\vec{\sigma}_{\mathcal{O}}]}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}$

In Appendix A (see equation (108)) we derive the identity:

$\lambda_{i\rightarrow j}^{t}\;\;\mu_{i\rightarrow
j}^{t}=\displaystyle\sum_{X_{\partial i\setminus
j}}^{t}r_{i}^{t}\,\Phi_{i}^{t}\prod_{k\in\partial i\setminus
j}\mu_{k\rightarrow i}^{t}$ (46)

which leads to:

$\displaystyle\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}$ $\displaystyle=$
$\displaystyle\Big{[}\prod_{\begin{subarray}{c}k=1\\\ k\neq
l\end{subarray}}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}\;\mu_{o_{l}\rightarrow i_{l}}^{t}$ (47)
$\displaystyle\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}$ $\displaystyle=$
$\displaystyle\Big{[}\prod_{\begin{subarray}{c}k=1\\\ k\neq
l\end{subarray}}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\sum_{X_{\partial
o_{l}\setminus
i_{l}}}^{t}r_{o_{l}}^{t}\,\Phi_{o_{l}}^{t}\Big{[}\prod_{k\in\partial
o_{l}\setminus i_{l}}\mu_{k\rightarrow o_{l}}^{t}\Big{]}$
$\displaystyle\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}$ $\displaystyle=$ $\displaystyle\sum_{X_{\partial o_{l}\setminus
i_{l}}}^{t}r_{o_{l}}^{t}\Big{[}\prod_{\begin{subarray}{c}k=1\\\ k\neq
l\end{subarray}}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\Phi_{o_{l}}^{t}\Big{[}\prod_{k\in\partial
o_{l}\setminus i_{l}}\mu_{k\rightarrow o_{l}}^{t}\Big{]}$ (48)

Now, what does the product $\Big{[}\prod_{\begin{subarray}{c}k=1\\\ k\neq
l\end{subarray}}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\Phi_{o_{l}}^{t}\Big{[}\prod_{k\in\partial
o_{l}\setminus i_{l}}\mu_{k\rightarrow o_{l}}^{t}\Big{]}$ in the right hand
side of (48) stands for? As in the example of figure (2), we took (47) and
removed the factor $\lambda_{o_{l}\rightarrow
i_{l}}^{t}\;\mu_{o_{l}\rightarrow i_{l}}^{t}$. Then, we included a local
weight, $\Phi_{o_{l}}^{t}$, and a group of cavity messages $\mu_{k\rightarrow
o_{l}}^{t}$. Therefore, in (48) we have all the local weights and cavity
messages that appear in the probability density of a bigger connected set,
which now includes all the neighbors of $o_{l}$.

Figure 2: This figure illustrates what happens when we apply the relation (46)
to equation (47) in the particular case of the connected set of the right
panels of figure (1). We took the outer node $k_{2}$ (see figure (1)) as
$o_{l}$. After including a new local weight $\Phi_{k_{2}}^{t}$ and two new
messages $\mu_{k_{21}\rightarrow k_{2}}^{t}$ and $\mu_{k_{21}\rightarrow
k_{2}}^{t}$, we end up having the probability density of a bigger connected
set that includes all neighbors of $k_{2}$.

Equation (48) transforms into:

$\displaystyle\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow
i_{k}}^{t}\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}\Big{]}\>\lambda_{o_{l}\rightarrow
i_{l}}^{t}$ $\displaystyle=$ $\displaystyle\sum_{X_{\partial o_{l}\setminus
i_{l}}}^{t}r_{o_{l}}^{t}\,Q^{t}(X_{\partial o_{l}\setminus
i_{l}},\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$ (49)

Thus, we can rewrite (44) and (45) to get:

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle-\Delta
t\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}\Big{[}\sum_{l=1}^{n}r_{i_{l}}^{t}Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})+\sum_{l=1}^{m}\sum_{X_{\partial
o_{l}\setminus i_{l}}}^{t}r_{o_{l}}^{t}\,Q(X_{\partial o_{l}\setminus
i_{l}},\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})\Big{]}$ $\displaystyle
I_{1}$ $\displaystyle=$ $\displaystyle-\Delta
t\Big{[}\sum_{l=1}^{n}r_{i_{l}}(\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})+\sum_{l=1}^{m}\sum_{\sigma_{\partial
o_{l}\setminus i_{l}}}r_{o_{l}}(\sigma_{o_{l}},\sigma_{\partial
o_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})\Big{]}$ (50)
$\displaystyle I_{2}$ $\displaystyle=$ $\displaystyle\Delta
t\sum_{l=1}^{n}\sum_{\vec{X}_{\mathcal{O}}|\vec{\sigma}_{\mathcal{O}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}]}^{t}r_{i_{l}}^{t}\,Q(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})+\Delta
t\sum_{l=1}^{m}\sum_{\vec{X}_{\mathcal{O}}|F_{o_{l}}[\vec{\sigma}_{\mathcal{O}}]}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}r_{o_{l}}^{t}\,Q(X_{\partial
o_{l}\setminus i_{l}},\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})$
$\displaystyle I_{2}$ $\displaystyle=$ $\displaystyle\Delta
t\Big{[}\sum_{l=1}^{n}r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}])+$
(51)
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>+\sum_{l=1}^{m}\sum_{\sigma_{\partial
o_{l}\setminus i_{l}}}r_{o_{l}}(-\sigma_{o_{l}},\sigma_{\partial
o_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},F_{o_{l}}[\vec{\sigma}_{\mathcal{O}}],\vec{\sigma}_{\mathcal{I}})\Big{]}$

where we have returned to the longer notation for the rates $r_{i}^{t}$ (see
equation (42)).

Putting (43), (50) and (51) together:

$\displaystyle P^{t+\Delta
t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ $\displaystyle=$
$\displaystyle I_{0}+I_{1}+I_{2}+o(\Delta t)$ $\displaystyle P^{t+\Delta
t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ $\displaystyle=$
$\displaystyle
P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})\\!-\\!\Delta
t\\!\sum_{l=1}^{n}\\!\Big{[}r_{i_{l}}(\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})\\!-\\!r_{i_{l}}(\\!-\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}])\Big{]}-$
(52) $\displaystyle-\Delta t\sum_{l=1}^{m}\sum_{\sigma_{\partial
o_{l}\setminus i_{l}}}r_{o_{l}}(\sigma_{o_{l}},\sigma_{\partial
o_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})+$
$\displaystyle+\Delta t\sum_{l=1}^{m}\sum_{\sigma_{\partial o_{l}\setminus
i_{l}}}r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},F_{o_{l}}[\vec{\sigma}_{\mathcal{O}}],\vec{\sigma}_{\mathcal{I}})+o(\Delta
t)$

And finally:

$\displaystyle\frac{dP^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})}{dt}$
$\displaystyle=$
$\displaystyle-\sum_{l=1}^{n}\Big{[}r_{i_{l}}(\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})-r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\vec{\sigma}_{\mathcal{O}},F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}])\Big{]}-$
(53) $\displaystyle-\sum_{l=1}^{m}\sum_{\sigma_{\partial o_{l}\setminus
i_{l}}}\Big{[}r_{o_{l}}(\sigma_{o_{l}},\sigma_{\partial
o_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},F_{o_{l}}[\vec{\sigma}_{\mathcal{O}}],\vec{\sigma}_{\mathcal{I}})\Big{]}$

Equation (53) is like a master equation for the combined set of variables
$\\{\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}}\\}$, but it has a
very peculiar structure. The first sum represents the contribution to the
derivative due to flipping rates of spins sitting at inner nodes, and the
second and third lines are the contribution related to outer nodes. Precisely
there we can notice that the time derivative of
$P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ depends on
probability densities defined over bigger connected sets:
$P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$. This means
that this is not a closed system of equations, and we need to complement it
with other relations in order to obtain the time dependence of its variables.

In this subsection we have presented a general method to obtain the exact
differential equation of any local probability density in a tree-like graph.
Although we have worked with connected sets of the first category (see figure
(1)), the equation for connected sets of the second category is just a
marginal of (53).

###  III.3 Equations for cavity probability densities

To proceed further, it is necessary to find a way to close the set of
equations (53). To do this, we need to introduce first a new set of master-
like equations. They will represent Cavity Master Equations to be solved
separately. With them we may write $P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ in terms of the
quantities $P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ (see
below), therefore closing the full system.

We start from the idea that, in a tree-like graph, any local probability
density $Q$ of the histories of a connected set can be written as:

$Q^{t}(\vec{X}_{\mathcal{O}},\vec{X}_{\mathcal{I}})=\Big{[}\prod_{k=1}^{m}\mu_{o_{k}\rightarrow(o_{k}i_{k})}^{t}(X_{o_{k}}|X_{i_{k}})\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}(X_{i_{k}},X_{\partial
i_{k}})\Big{]}$ (54)

or as a marginal of this kind of expressions. The densities $Q$ are defined
over the space of all possible configurations of the $N$-vector $\vec{X}$ of
nodes histories.

On the other hand, we may introduce a cavity probability density, $q$, defined
over a reduced space where we partially or completely fix some of the
histories. Thus, $Q$ and $q$ are fundamentally different.

Here, we will deal only with cavity probability densities which are defined
over connected sets. As in previous subsection, we can focus on connected sets
of the first kind without any loss in generality. The reduction of the
configuration space that gives birth to a cavity treatment occurs in this case
by fixing the history of one outer node. We can write one of those $q$
densities in terms of local weights and cavity messages from a probability
density $Q$ just by removing the cavity messages $\mu$ corresponding to the
fixed outer node. For example:

$q^{t}(\vec{X}_{\mathcal{O}\setminus
o_{j}},\vec{X}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})=\Big{[}\prod_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{m}\mu_{o_{k}\rightarrow(o_{k}i_{k})}^{t}(X_{o_{k}}|X_{i_{k}})\Big{]}\Big{[}\prod_{k=1}^{n}\Phi_{i_{k}}^{t}(X_{i_{k}},X_{\partial
i_{k}})\Big{]}$ (55)

We say that $q^{t}(\vec{X}_{\mathcal{O}\setminus
o_{j}},\vec{X}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})$ is defined in the
cavity where $X_{o_{j}}$ is fixed, which is represented by the symbol
$\big{|}\big{|}$. Figure (3) illustrates this idea by revisiting a particular
example shown in previous subsection.

By marginalizing (55) we get instantaneous cavity probability densities:

$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})=\sum_{\vec{X}_{\mathcal{O}\setminus
o_{j}}|\vec{\sigma}_{\mathcal{O}\setminus
o_{j}}}^{t}\sum_{\vec{X}_{\mathcal{I}}|\vec{\sigma}_{\mathcal{I}}}^{t}q^{t}(\vec{X}_{\mathcal{O}\setminus
o_{j}},\vec{X}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})$ (56)

whose time differentiation is analogous to what we did in previous subsection.
We will have $n$ terms like the ones in the first line of (53), and $m-1$
terms like the ones in the second and third lines of the same equation, only
that now everything is defined in the cavity where $X_{o_{j}}$ is fixed.

Figure 3: This figure illustrates the relation between probability densities
$Q$ (left panel) and cavity probability densities $q$ (right panel) defined
over connected sets. In this case $q$ is obtained by removing the cavity
message $\mu_{j\rightarrow i}^{t}$, and we say that $q$ is defined in the
cavity where $X_{j}$ is fixed.

Thus, proceeding as in III.2 it is easy to verify that:

$\displaystyle\frac{dp^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})}{dt}$
$\displaystyle=$
$\displaystyle-\sum_{l=1}^{n}\\!\Big{[}r_{i_{l}}(\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})-$ (57)
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i_{l}}(\\!-\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}]\>\big{|}\big{|}\>X_{o_{j}})\Big{]}-$
$\displaystyle-\sum_{l=1}^{m}\sum_{\sigma_{\partial o_{l}\setminus
i_{l}}}\Big{[}r_{o_{l}}(\sigma_{o_{l}},\sigma_{\partial
o_{l}})\,p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},F_{o_{l}}[\vec{\sigma}_{\mathcal{O}\setminus
o_{j}}],\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})\Big{]}$

As (53), equation (57) is exact in tree-like graphs. However, none of them can
be analytically or numerically solved. Equations (53) and (57) contain
expressions for the time derivatives of densities
$P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ and
$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>X_{o_{j}})$, but these
expressions are not closed on the same variables. They both depend on
probability densities defined over bigger sets of nodes.

###  III.4 Closure and hierarchical approximations

Although cavity probability densities are defined over configuration spaces
which are fundamentally different than the full configuration space of all
histories, we can use them to give a closure to equations like (53).

First off all, let’s make a Markovian approximation in equation (57) buy
substituting the full history $X_{o_{j}}$ by its final state $\sigma_{o_{j}}$:

$\displaystyle\frac{dp^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})}{dt}$
$\displaystyle=$
$\displaystyle-\sum_{l=1}^{n}\\!\Big{[}r_{i_{l}}(\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})-$ (58)
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i_{l}}(\\!-\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},F_{i_{l}}[\vec{\sigma}_{\mathcal{I}}]\mid\sigma_{o_{j}})\Big{]}-$
$\displaystyle-\sum_{l=1}^{m}\sum_{\sigma_{\partial o_{l}\setminus
i_{l}}}\Big{[}r_{o_{l}}(\sigma_{o_{l}},\sigma_{\partial
o_{l}})\,p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i_{l}}(-\sigma_{i_{l}},\sigma_{\partial
i_{l}})\,p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},F_{o_{l}}[\vec{\sigma}_{\mathcal{O}\setminus
o_{j}}],\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})\Big{]}$

Equation (58) gives the time derivative of a probability density
$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})$ that can
be defined over any connected set in the graph. Its first two lines represent
the contribution to the derivative due to flipping rates of spins sitting at
inner nodes, and the last two lines correspond to the contribution related to
outer nodes (except for $\sigma_{j}$). As in (53), this equation is not
closed. Densities like $p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})$ are
defined over connected sets which include more nodes than
$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})$. Again, we
need closure relations for (58).

####  III.4.1 First-order closures

It is useful at this time to re-visit the closure presented in CME-PRE ,
considered now as the simplest approximation we can make to equations like
(58). Indeed, the CME at CME-PRE  is written for the cavity densities
$p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$, each one involving only two
spins. Following (58), the equation for this variables reads:

$\displaystyle\frac{dp^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})}{dt}$
$\displaystyle=$ $\displaystyle-\\!\\!\sum_{\sigma_{\partial i\setminus
j}}\\!\\!\Big{\\{}r_{i}(\sigma_{i},\sigma_{\partial
i})\>p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})-r_{i}(-\sigma_{i},\sigma_{\partial
i})\>p^{t}(\sigma_{\partial i\setminus
j},-\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\Big{\\}}$ (59)

To close (59), we need to write $p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ in terms of
$p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$. In order to do so in a
transparent way, we make first a small detour.

Consider the simplest case of a stochastic process with three random variables
$A$, $B$ and $C$. There, an example of conditional cavity probability density
reads:

$p^{t}(A=a\mid
B=b\>\big{|}\big{|}\>C=c)=\frac{p^{t}(a,b\>\big{|}\big{|}\>c)}{p^{t}(b\>\big{|}\big{|}\>c)}$
(60)

where $a$, $b$ and $c$ are values of the variables $A$, $B$ and $C$.

Notice that we used two different symbols to represent conditional and cavity
relations: $\mid$ and $\big{|}\big{|}$, respectively. We intend to remark that
these are conceptually different. In this example the conventional conditional
probability density $P^{t}(a\mid b)$ is a measure of the set of configurations
with $A(t)=a$ restricted to the sub-space where $B(t)=b$. However, this
magnitude is defined in a stochastic process where $A$, $B$ and $C$ can take
any of its possible values in the full configuration space. On the other hand,
a cavity condition $p^{t}(a\big{|}\big{|}\>c)$ is a modification to our
stochastic problem whose realizations are now defined on a configuration space
where $C(t)=c$ in all cases. If that is not clear, a look back to (55) shows
that a cavity density is a result of removing the probabilistic weight that
corresponds to the fixed variable ($C$ in this case). The question about what
might be happening with $C$ at time $t$ is simply not asked.

It is in that context that we wrote $p^{t}(a\mid b\>\big{|}\big{|}\>c)$. This
density is a measure of the set of configurations with $A(t)=a$ restricted to
the sub-space where $B(t)=b$ of a stochastic process whose realizations always
have $C(t)=c$.

Now we can write the exact relation:

$p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})=p^{t}(\sigma_{\partial i\setminus
j}\mid\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\,p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$
(61)

In CME-PRE , this relation is approximated by:

$p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})=p^{t}(\sigma_{\partial i\setminus
j}\mid\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\,p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\approx\Big{[}\prod_{k\in\partial
i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}\sigma_{i})\Big{]}\,p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$
(62)

The substitution in (62) can be interpreted as a series of changes made to
$p^{t}(\sigma_{\partial i\setminus
j}\mid\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$. First, the cavity relation
with $\sigma_{j}$ is dropped to get $p^{t}(\sigma_{\partial i\setminus
j}\mid\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\approx p^{t}(\sigma_{\partial
i\setminus j}\mid\sigma_{i})$. Then, the conditional relation with
$\sigma_{i}$ is substituted by a cavity relation $p^{t}(\sigma_{\partial
i\setminus j}\mid\sigma_{i})\approx p^{t}(\sigma_{\partial i\setminus
j}\>\big{|}\big{|}\>\sigma_{i})$. These approximations increase their accuracy
when the correlations between $\sigma_{\partial i\setminus j}$, $\sigma_{j}$
and $\sigma_{i}$ decrease CME-PRE .

Completely dropping the dependence on $\sigma_{j}$ looks more drastic than
replacing $\mid\sigma_{i}$ by $\>\big{|}\big{|}\>\sigma_{i}$, but we must
notice that node $j$ is farther apart from any node in $\partial i\setminus j$
than node $i$, thus, $\sigma_{\partial i\setminus j}$ should be less
correlated with $\sigma_{j}$ than with $\sigma_{i}$. Finally, it is easy to
see that $p^{t}(\sigma_{\partial i\setminus j}\>\big{|}\big{|}\>\sigma_{i})$
can be exactly factorized to get (62), and we obtain the CME as presented in
CME-PRE .

$\displaystyle\frac{dp^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})}{dt}$
$\displaystyle=$ $\displaystyle-\sum_{\sigma_{\partial i\setminus
j}}\Big{\\{}r_{i}(\sigma_{i},\sigma_{\partial i})\Big{[}\prod_{k\in\partial
i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}\>\sigma_{i})\Big{]}p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})-$
(63)
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i}(-\sigma_{i},\sigma_{\partial
i})\Big{[}\prod_{k\in\partial i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}-\sigma_{i})\Big{]}p^{t}(-\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\Big{\\}}$

In CME-PRE ; CME-Pspin ; CME-PRL ; CME-AvAndBP  these kind of equations are
integrated together with:

$\displaystyle\frac{dP^{t}(\sigma_{i})}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{\sigma_{\partial
i}}\Big{\\{}r_{i}(\sigma_{i},\sigma_{\partial i})\,\Big{[}\prod_{k\in\partial
i}p^{t}(\sigma_{k}\>\big{|}\big{|}\sigma_{i})\Big{]}\,P^{t}(\sigma_{i})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i}(-\sigma_{i},\sigma_{\partial
i})\,\Big{[}\prod_{k\in\partial
i}p^{t}(\sigma_{k}\>\big{|}\big{|}-\sigma_{i})\Big{]}\,P^{t}(-\sigma_{i})\Big{\\}}$

where a similar factorization is used.

$P^{t}(\sigma_{\partial i},\sigma_{i})=P^{t}(\sigma_{\partial
i}\mid\sigma_{i})P^{t}(\sigma_{i})\approx p^{t}(\sigma_{\partial
i}\>\big{|}\big{|}\>\sigma_{i})\,P^{t}(\sigma_{i})=\Big{[}\prod_{k\in\partial
i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}\sigma_{i})\Big{]}\,P^{t}(\sigma_{i})$ (65)

Although these equations had proven their utility in previous works, we can go
beyond equation (63) by directly obtaining the time dependence of the
probability densities $p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$. These are defined over a larger
connected sets, and according to the general equation (58) its time derivative
can be written as:

$\displaystyle\frac{dp^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})}{dt}$ $\displaystyle=$
$\displaystyle-r_{i}(\sigma_{i},\sigma_{\partial i})\,p^{t}(\sigma_{\partial
i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})+r_{i}(-\sigma_{i},\sigma_{\partial
i})\,p^{t}(\sigma_{\partial i\setminus
j},-\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ $\displaystyle-\sum_{l\in\partial
i\setminus j}\sum_{\sigma_{\partial l\setminus
i}}\Big{\\{}r_{l}(\sigma_{l},\sigma_{\partial l})\,p^{t}(\sigma_{\partial
l\setminus i},\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{l}(-\sigma_{l},\sigma_{\partial
l})\,p^{t}(\sigma_{\partial l\setminus i},F_{l}[\sigma_{\partial i\setminus
j}],\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\Big{\\}}$

Again, we need to close these equations introducing some relation between the
probability densities $p^{t}(\sigma_{\partial l\setminus i},\sigma_{\partial
i\setminus j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ and
$p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$. We can always concentrate on
conditional cavity densities, which allow us to write:

$p^{t}(\sigma_{\partial l\setminus i},\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})=p^{t}(\sigma_{\partial l\setminus
i}\mid\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\,p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ (67)

As before, we can take $p^{t}(\sigma_{\partial l\setminus
i}\mid\sigma_{\partial i\setminus j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$
and drop the dependence on $\sigma_{j}$. As all the nodes in $\partial
i\setminus j$, except for $l$, are at the same distance from any node in
$\partial l\setminus i$ than node $j$, the corresponding spins will be in
average equally correlated with $\sigma_{\partial l\setminus i}$ than
$\sigma_{j}$. Let’s also drop the conventional conditional dependence on these
nodes to get the relation $p^{t}(\sigma_{\partial l\setminus
i}\mid\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\approx p^{t}(\sigma_{\partial
l\setminus i}\mid\sigma_{l},\sigma_{i})$. This procedure is now more accurate
than with (63), because we are neglecting correlations associated with nodes
which are farther apart in the graph.

Then, we replace the conventional conditional relation with $\sigma_{i}$ by a
cavity relation to obtain:

$p^{t}(\sigma_{\partial l\setminus i}\mid\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\approx p^{t}(\sigma_{\partial
l\setminus
i}\mid\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})=\frac{p^{t}(\sigma_{\partial
l\setminus i},\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})}{\sum_{\sigma_{\partial
l\setminus i}}p^{t}(\sigma_{\partial l\setminus
i},\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})}$ (68)

With this, we have closed (III.4.1) through approximations for conditional
cavity probability densities.

Equations (III.4.1) and (68) must also be combined with differential equations
for the probability densities $P^{t}$. In this case we can use, for example:

$\displaystyle\frac{dP^{t}(\sigma_{\partial i},\sigma_{i})}{dt}$
$\displaystyle=$ $\displaystyle-r_{i}(\sigma_{i},\sigma_{\partial
i})\,P^{t}(\sigma_{\partial i},\sigma_{i})+r_{i}(-\sigma_{i},\sigma_{\partial
i})\,P^{t}(\sigma_{\partial i},-\sigma_{i})$ $\displaystyle-\sum_{l\in\partial
i}\sum_{\sigma_{\partial l\setminus
i}}\Big{\\{}r_{l}(\sigma_{l},\sigma_{\partial l})\,p^{t}(\sigma_{\partial
l\setminus
i}\mid\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})\,P^{t}(\sigma_{\partial
i},\sigma_{i})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{l}(-\sigma_{l},\sigma_{\partial
l})\,p^{t}(\sigma_{\partial l\setminus
i}\mid-\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})\,P^{t}(F_{l}[\sigma_{\partial
i}],\sigma_{i})\Big{\\}}$

which involves essentially the same approximations than (68):

$P^{t}(\sigma_{\partial l\setminus i},\sigma_{\partial
i},\sigma_{i})=P^{t}(\sigma_{\partial l\setminus i}\mid\sigma_{\partial
i},\sigma_{i})\,P^{t}(\sigma_{\partial i},\sigma_{i})\approx
p^{t}(\sigma_{\partial l\setminus
i}\mid\sigma_{l}\>\big{|}\big{|}\sigma_{i}\>)\,P^{t}(\sigma_{\partial
i},\sigma_{i})$ (70)

The performance of the set of equations (III.4.1)-(70) that we call in the
following CME-2 will outperform the CME as presented in CME-PRE  (see
equations (63) and (III.4.1)), fact that we illustrate in section IV.

####  III.4.2 General closure

Although we gave an initial idea of how to use the dynamic cavity method to
obtain closed systems of differential equations for local probability
densities, our goal is to provide a general method that works for any equation
like (57). This requires some sense of structure and order that should
facilitate the reader’s comprehension. We will move among connected sets of
the first kind (see figure (1)) to establish a hierarchy using the concept of
distance within the graph.

We can find a path that connects every pair of nodes in (74). Furthermore, for
each pair there is always a path with minimum length, and we can define a
distance between nodes as the length of that minimal path. The distance
$d(k,l)$ between nodes $k$ and $l$, with $k\neq l$, is then a positive
integer.

From this, we can classify the connected set $\\{\mathcal{O}\setminus
o_{j},\mathcal{I},o_{j}\\}$ according to the distances between the node
$o_{j}$ (the origin) and all the outer nodes $\mathcal{O}\setminus o_{j}$.
When the maximum distance $\max_{o_{l}\in\mathcal{O}\setminus
o_{j}}\\{d(o_{j},o_{l})\\}$ is equal to the integer $Z$, we will say that the
connected set is in the class CS-Z. If all the outer nodes (except $o_{j}$)
are at the same distance from $o_{j}$, i.e., if
$d(o_{j},o_{l})=d(o_{j},o_{k})$ for all $o_{k},o_{l}\in\mathcal{O}\setminus
o_{j}$, we will say that the connected set is regular and use the notation
rCS-Z to refer to that specific set. Of course, every rCS-Z belongs to the
class CS-Z. Figure (4) illustrates this classification.

Figure 4: A cavity probability density
$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})$ is defined
over a general connected set. This figure illustrates the classification of
these sets according to the distance between the outer nodes and the origin
$o_{j}$. (Left panel): The maximum distance is $d=3$. Thus, the set belongs to
the class CS-3. In this example all distances are not equal. (Right panel):
Every outer node is at a distance $d=3$ from $o_{j}$. Thus, the set is rCS-3
and belongs to CS-3.

In this context can use the notation:

$p^{t}(\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})\equiv
p^{t}(\vec{\sigma}_{Z},\vec{\sigma}_{Z-1},\ldots,\vec{\sigma}_{2},\vec{\sigma}_{1}\>\big{|}\big{|}\>\sigma_{o_{j}})\equiv
p^{t}(\vec{\sigma}_{Z},\vec{\sigma}_{Z-1},\ldots,\vec{\sigma}_{2},\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
(71)

where $i_{j}$ is the inner node connected to $o_{j}$.

In (71), the vector $\vec{\sigma}_{z}$ ($1\leq z\leq Z$) contains the spins
$\sigma_{k}$ such that $d(k,o_{j})=z$, with the condition that
$k\in\mathcal{O}$ or $k\in\mathcal{I}$. The cavity probability density in (71)
is defined over a connected set that belongs to the class CS-Z, which we will
represent with the symbol $p^{t}(\text{{CS-Z}})$. In order to simplify our
language, we will simply say that all $p^{t}(\text{{CS-Z}})$ are _in_ CS-Z.

Using our new notation (71), the cavity probability densities in the last two
lines of (58), which are in CS-(Z+1), can be rewritten as:

$p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\>\sigma_{o_{j}})\equiv
p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}_{Z},\vec{\sigma}_{Z-1},\ldots,\vec{\sigma}_{2},\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
(72)

where $\vec{\sigma}_{Z+1}$ is exactly the same than $\sigma_{\partial
o_{l}\setminus i_{l}}$.

Notice that every $p^{t}(\text{{CS-Z}})$ is a marginal of a cavity probability
density defined over the corresponding regular set rCS-Z (see left panels of
(4)). In view of all this, what we need to write is a system of equations
like:

$\displaystyle\frac{dp^{t}(\text{{rCS-Z}})}{dt}=F\big{[}p^{t}(\text{{rCS-Z}})\big{]}$
(73)

which in our case is equivalent to find a relation
$p^{t}(\text{{CS-(Z+1)}})=f\big{[}p^{t}(\text{{CS-Z}})\big{]}$.

With this classification of connected sets according, essentially, to their
size, we are ready to find that $f[\cdot]$ function. Let’s write:

$p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})=\frac{p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})}{p^{t}(\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})}$
(74)

If the distance $Z+1$ is large compared to system’s correlation length, we can
assume that the fact of having fixed the value of $\sigma_{o_{j}}$ at time $t$
is irrelevant for the state $\vec{\sigma}_{Z+1}$ in the left hand side of
(74):

$p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})\approx
p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}_{Z},\ldots,\vec{\sigma}_{2},\sigma_{i_{j}})$
(75)

Following the same reasoning, we can drop all the spins in
$\\{\vec{\sigma}_{Z},\ldots,\vec{\sigma}_{2}\\}$ whose distance to the nodes
in $\vec{\sigma}_{Z+1}$ is larger than $Z$. We have already done this kind of
approximations when writing the equations for
$p^{t}(\sigma_{i}\>\big{|}\big{|}\sigma_{j})$ and $p^{t}(\sigma_{\partial
i\setminus j},\sigma_{i}\>\big{|}\big{|}\sigma_{j})$, here we are just
generalizing. The remaining spins
$\\{\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2},\sigma_{i_{j}}\\}$
are defined over a connected set of the class CS-Z, with origin in $i_{j}$.

Furthermore, when the spatial correlation between $\vec{\sigma}_{Z+1}$ and
$\sigma_{i_{j}}$ is low, we can make the approximation CME-PRE :

$p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2},\sigma_{i_{j}})\approx
p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})=\frac{p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}{p^{t}(\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}$
(76)

Notice that the distance between $i_{j}$ and the nodes related to
$\vec{\sigma}_{Z+1}$ is equal to $Z$, so we are assuming that the correlation
decays noticeably at that distance.

From (75) and (76) we have:

$p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})\approx\frac{p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}{\sum_{\vec{\sigma}_{Z+1}}p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}$
(77)

and the probability densities (72) become:

$\displaystyle
p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
$\displaystyle=$ $\displaystyle
p^{t}(\vec{\sigma}_{Z+1}\mid\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})\>p^{t}(\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
$\displaystyle
p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
$\displaystyle\approx$
$\displaystyle\frac{p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}{\sum_{\vec{\sigma}_{Z+1}}p^{t}(\vec{\sigma}_{Z+1},\vec{\sigma}^{\prime}_{Z},\ldots,\vec{\sigma}^{\prime}_{2}\>\big{|}\big{|}\>\sigma_{i_{j}})}\>p^{t}(\vec{\sigma}_{Z},\ldots,\sigma_{i_{j}}\>\big{|}\big{|}\>\sigma_{o_{j}})$
(78)

Thus, we have managed to write the function
$f\big{[}p^{t}(\text{{CS-Z}})\big{]}$ that gives any cavity probability
density in CS-(Z+1) in terms of several densities in CS-Z. This constitutes a
closure for (58), and we can now numerically solve these equations.

It is important to say that after solving for $p^{t}(\text{{CS-Z}})$, we can
compute all the cavity probability densities in CS-1, CS-2, …, CS-(Z-1) as its
marginals.

As for (53), we can write $P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ in terms of
$P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ and a cavity
probability density whose value we can obtain from the numerical integration
of (58) with the relation (78). We do:

$\displaystyle P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$
$\displaystyle=$ $\displaystyle P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}}\mid\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})\,P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$
$\displaystyle P^{t}(\sigma_{\partial o_{l}\setminus
i_{l}},\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$
$\displaystyle\approx$ $\displaystyle p^{t}(\sigma_{\partial o_{l}\setminus
i_{l}}\mid\vec{\sigma}_{\mathcal{O}\setminus
o_{j}},\vec{\sigma}_{\mathcal{I}}\>\big{|}\big{|}\sigma_{o_{j}})\,P^{t}(\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$
(79)

Here we changed the conditional probability density $P^{t}(\sigma_{\partial
o_{l}\setminus
i_{l}}\mid\vec{\sigma}_{\mathcal{O}},\vec{\sigma}_{\mathcal{I}})$ by a
conditional cavity density. Nevertheless, we introduced the cavity condition
by fixing only the spin $\sigma_{o_{j}}$. Again, this substitution should be
accurate when the distance between $o_{j}$ and the nodes in $\partial
o_{l}\setminus i_{l}$ is large enough compared with correlation length.
Finally, equations (58), (74), (78) and (79) give a closure to the
differential equation (53) and we are ready to numerically obtain the time
dependence of system’s observables.

It is important to say that our closure relations could be used in several
ways, but the key point will always be how to choose the distance $Z$ at which
we neglect correlations and substitute conventional conditional relations by
conditional cavity relations. This distance will be the parameter that defines
what we will call a _level of approximation_.

Summarizing, we define a level of approximation through the differential
equations for the cavity probability densities (58). If we write them for
densities $p^{t}(\text{{CS-Z}})$ we say that we are in the $Z^{\text{th}}$
level of approximation. Within this context, equations (63) and (III.4.1)
constitute a first level approximation, and we will say they are _in_ the
level CME-1. On the other hand, equations (III.4.1)-(70) are _in_ the second
level, which is CME-2.

##  IV Numeric results

The information about the nature of the interactions and the local dynamics
goes into the spin-flip rates $r_{i}(\sigma_{i},\sigma_{\partial i})$, which
of course depend on the graph’s structure. Here, we will work over two
families of diluted random graphs which are locally tree-like: Erdos-Renyi
erdHos1959random  and Random Regular Graphs.

In practice, once interactions are set, we select an initial condition for all
the probability densities ($P^{t_{0}}$ and $p^{t_{0}}$), and the integration
of the equations gives the full $P^{t}$ and $p^{t}$. In subsection IV.1 we
will explore the numerical differences obtained with first and second levels
of approximation, CME-1 and CME-2, into two well-known spin models from
statistical mechanics, the Ising ferromagnet, and the Viana-Bray spin-glass,
both defined over Erdos-Renyi graphs. In IV.2 we compare the CME-2 and CME-3
with the dynamical independent neighbor approximation presented in
Semerjian_dyn_2004 , for the dynamics of the Ising ferromagnet. Finally, at
IV.3 we also contrast the results of CME-2, with the ones of Pair Quenched
Mean-Field Approximation for the dynamics of the susceptible-infectious-
susceptible (SIS) model for epidemics cator2012second ; mata2013pair .

###  IV.1 First and second levels of approximation

One of the traditional forms that theorists choose for the spin flipping
probabilities per time unit is Glauber’s rule Glauber63 :

$r_{i}(\sigma_{i},\sigma_{\partial
i})=\frac{1}{2}\big{[}1-\sigma_{i}\tanh\big{(}\beta\sum_{k\in\partial
i}J_{ki}\sigma_{k}+\beta h_{i}\big{)}\big{]}$ (80)

As usual, the spin variables $\sigma_{i}$ can take the values $\sigma_{i}=\pm
1$. The interaction between each pair of connected spins
$(\sigma_{i},\sigma_{j})$ is encoded in the couplings parameters $J_{ij}=\pm
1$ and can be either satisfied ($J_{ij}\sigma_{i}\sigma_{j}=1$) or unsatisfied
($J_{ij}\sigma_{i}\sigma_{j}=-1$). On the other hand, the $h_{i}$ are local
fields that we will set here to zero. The parameter $\beta$ is the inverse of
the temperature $T$.

We can then write:

$r_{i}(\sigma_{i},\sigma_{\partial i})\equiv
r_{i}(c_{i},u_{i})=\frac{1}{2}\big{[}1-\tanh\big{(}\beta(2c_{i}-u_{i})\big{)}\big{]}$
(81)

where $c_{i}$ is spin’s connectivity and we can use the Kronecker’s delta to
write:

$u_{i}\equiv u_{i}(\sigma_{i},\sigma_{\partial i})=\sum_{k\in\partial
i}\delta_{J_{ki},\sigma_{k}\sigma_{i}}$ (82)

When the number $u_{i}$ of unsatisfied interactions between $i$ and its
neighbors is large, there is a high probability that $\sigma_{i}$ changes to
$-\sigma_{i}$ (see equation (81)).

Let’s assume we have an Erdos-Renyi graph with $N$ nodes, and in every node a
spin variable. The corresponding ferromagnetic Ising model is the result of
setting $J_{ij}=1$ for all the connected pairs $(i,j)$. With this definition
and with the rules (81) we can make Kinetic Monte Carlo simulations of the
system’s dynamics.

Figure (5) shows results of these simulations alongside the semi-analytical
output of our dynamic cavity method in the first (CME-1) and second (CME-2)
levels of approximation. Top-left panel contains the time evolution of
system’s magnetization (equation (83)), while top-right panel illustrates the
corresponding error for the local magnetization (equation (84)) between a
theoretical (TH) approach and the simulation (MC) :

$\displaystyle m(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{\sigma_{i}}\sigma_{i}\,P^{t}(\sigma_{i})\equiv\frac{1}{N}\sum_{i=1}^{N}m_{i}(t)$
(83) $\displaystyle\delta m^{\text{TH}}(t)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{N}\sum_{i=1}^{N}\big{(}m_{i}^{\text{TH}}(t)-m_{i}^{\text{MC}}(t)\big{)}^{2}}$
(84)

The behavior depicted in top-left panel is typical of ferromagnets. A
transition between two regimes occurs at some $T_{C}$ (when the average
connectivity is $c=3$ we have $T_{C}\approx 2.89$). Steady-state magnetization
is zero for high temperatures ($T>T_{C}$), and at low temperatures ($T<T_{C}$)
we have non-zero magnetization for all times.

Both levels of approximation, CME-1 and _CME-2_ , give a good qualitative and
quantitative description of the magnetization obtained from the simulations,
especially at the steady-state. The results of CME-2 are particularly
accurate, and this happens in the steady-state and also at the transient. As
top-right panel of (5) shows, the equations in the level CME-2 reproduce local
magnetizations $m_{i}^{\text{MC}}(t)$ also quite well, with errors of order
$10^{-3}-10^{-2}$.

Figure 5: Dynamics of the Ising ferromagnet according to Kinetic Monte Carlo
simulations (points), and the integration of equations in levels CME-1 (thick
lines) and CME-2 (dashed lines). In all cases an initially fully-magnetized
system evolves in time in contact with a heat bath at a given temperature.
Calculations were done for Erdos-Renyi graphs of size $N=5000$ with mean
connectivity $c=3$. Top-left panel: Time evolution of system’s magnetization
(see equation (83)). Top-right panel: Time evolution of the error (84) between
local magnetizations predicted by CME-1, 2 and the results of Kinetic Monte
Carlo simulations. Bottom-left panel: Time evolution of system’s energy
density (see equation (85)). Bottom-right panel: Time evolution of the error
(86) written for the expected values of local energy terms.

Bottom-left panel of (5) shows the time evolution of ferromagnet’s energy
density:

$e(t)=-\frac{1}{2N}\sum_{i\neq
j}\sum_{\sigma_{i}}\sum_{\sigma_{j}}P^{t}(\sigma_{i},\sigma_{j})J_{ij}\sigma_{i}\sigma_{j}=\frac{1}{2}\sum_{i\neq
j}e_{ij}(t)$ (85)

which is a measure of how many unsatisfied interactions are in the system.
Again, our dynamic cavity method gives very accurate results, even for the
expected values of local energy terms $e_{ij}(t)$. Bottom-right panel of (5)
shows the time evolution of local error:

$\delta e^{\text{TH}}(t)=\sqrt{\frac{1}{2M}\sum_{i\neq
j}\big{(}e_{ij}^{\text{TH}}(t)-e_{ij}^{\text{MC}}(t)\big{)}^{2}}$ (86)

where $M$ is the number of connected pairs $(i,j)$. This error remains of
order $10^{-3}-10^{-2}$ for all times.

Another way of choosing the couplings is to draw each one from the bi-modal
distribution $d(J_{ij})=1/2\delta(J_{ij}-1)+1/2\delta(J_{ij}+1)$. In this
model, also defined over Erdos-Renyi random graphs, we have a transition to a
spin-glass phase at a finite temperature $T_{SG}$. As correlations play an
important role when temperature decreases, we know that the approximations we
made in III.4 can’t be as accurate as with the ferromagnet.

Figure (6) shows results for this model in a graph with $c=3$, where
$T_{SG}\approx 1.52$. As expected, below $T_{SG}$ the description is not as
successful as for the ferromagnetic case. However, it still holds that CME-2
performs better than CME-1. In the top-right panel, we see that for $0\leq
t\leq 100$ the errors are on a scale of $10^{-1}$, which is one order higher
than what we saw in the ferromagnet. Furthermore, even with the approximation
CME-2, the error seems to monotonically increase with time.

However, our theoretical description of the simulations is noticeably improved
just by going up one level in our hierarchical approximations. Not only errors
for local magnetizations are appreciably smaller with the level CME-2 than
with CME-1. If we look at the bottom panels of the figure we see how this
second level of approximation already gives a good description of the average
$e(t)$ even at very low temperatures like $T=0.25\ll T_{SG}$. However, this
result might be misleading because the errors for the expected values of local
energy terms are still of order $10^{-1}$.

Figure 6: Dynamics of the Viana-Bray spin-glass according to Kinetic Monte
Carlo simulations (points), and the integration of equations in levels _CME-1_
(thick lines) and _CME-2_ (dashed lines). In all cases an initially fully-
magnetized system evolves in time in contact with a heat bath at a given
temperature. Calculations were done for Erdos-Renyi graphs of size $N=1000$
with mean connectivity $c=3$. Top-left panel: Time evolution of system’s
magnetization (see equation (83)). Top-right panel: Time evolution of the
error (84) between local magnetizations predicted by CME-1, 2 and the results
of Kinetic Monte Carlo simulations. Bottom-left panel: Time evolution of
system’s energy density (see equation (85)). Bottom-right panel: Time
evolution of the error (86) written for the expected values of local energy
terms.

Results in figures (5) and (6) indicate that the performance of the approach
used in CME-PRE ; CME-Pspin ; CME-PRL ; CME-AvAndBP  (in CME-1) is
significantly improved just by using the next level of approximation. It is
important to say that the integration of equations in the second level,
although very accurate, does not take high efforts considering the
capabilities of present-day personal computers, not to mention high-
performance computers.

###  IV.2 Comparison with other theoretical approach

The hierarchical system of approximations introduced for the study of
stochastic local search algorithms SemerjianwalkSAT2003 ; BarthelwalkSAT2003
and systematized in Semerjian_dyn_2004  combines theoretical simplicity with
numerical accuracy in a variety of traditional models from statistical
mechanics. Here, we will compare our theoretical approach with the Dynamical
Independent Neighbor Approximation (DINA in what follows), which constitutes
the second level of the scheme presented in Semerjian_dyn_2004 .

In the case of an Ising ferromagnet with Glauber rates (81) defined over
random regular graphs, where every node has $K$ neighbors, the DINA works with
the exact equation:

$\displaystyle\frac{d\hat{P}^{t}(\sigma,u)}{dt}$ $\displaystyle=$
$\displaystyle-r(K,u)\,\hat{P}^{t}(\sigma,u)+r(K,K-u)\,\hat{P}^{t}(-\sigma,K-u)-$
$\displaystyle-\sum_{u^{\prime}=0}^{K-1}(K-u^{\prime})\,r(K,u^{\prime})\,\hat{P}^{t}(\sigma,u^{\prime})\Big{[}\hat{P}^{t}(u\mid\sigma,\sigma,u^{\prime})-\hat{P}^{t}(u-1\mid\sigma,\sigma,u^{\prime})\Big{]}-$
$\displaystyle-\sum_{u^{\prime}=0}^{K-1}(u^{\prime}\\!+\\!1)\,r(K,u^{\prime}+1)\,\hat{P}^{t}(\sigma,u^{\prime}+1)\Big{[}\hat{P}^{t}(u-1\mid\sigma,-\sigma,u^{\prime})-\hat{P}^{t}(u\mid\sigma,-\sigma,u^{\prime})\Big{]}$

together with the closure relations Semerjian_dyn_2004 :

$\displaystyle\hat{P}^{t}(\hat{u}\mid\sigma,\sigma,\hat{u}^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{(K-\hat{u})\hat{P}^{t}(\sigma,\hat{u})}{\sum_{u^{\prime}}(K-u^{\prime})\,\hat{P}^{t}(\sigma,u^{\prime})}$
(88) $\displaystyle P^{t}(\hat{u}\mid\sigma,-\sigma,\hat{u}^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{(\hat{u}+1)P^{t}(\sigma,\hat{u}+1)}{\sum_{u^{\prime}}u^{\prime}\,\hat{P}^{t}(\sigma,u^{\prime})}$
(89)

Here, $\hat{P}^{t}(\sigma,u)$ is the probability density of having a node with
spin $\sigma$ and $u$ unsatisfied interactions with its neighbors, where the
variable $u$ is an integer between zero and $K$. This can be written in terms
of the densities $P^{t}(\sigma_{i},\sigma_{\partial i})$ as:

$\hat{P}^{t}(\sigma,u)=\lim_{N\rightarrow\infty}\ll\frac{1}{N}\sum_{i=1}^{N}\sum_{\sigma_{i}}\sum_{\sigma_{\partial
i}}P_{\xi_{K}(N)}^{t}(\sigma_{i},\sigma_{\partial
i})\,\delta_{\sigma,\sigma_{i}}\>\>\delta_{(\sum_{k\in\partial
i}\sigma_{k}),K-2u}\gg_{\xi_{K}(N)}$ (90)

where $\xi_{K}(N)$ is the ensemble of random regular graphs with connectivity
$K$ and size $N$, and $\delta_{x,y}$ is a Kronecker delta evaluated at
$(x,y)$. The symbol $\ll\cdot\gg_{\xi_{K}(N)}$ represents an average over this
ensemble, with proper probabilistic weights.

We can similarly define the probability density
$\hat{P}^{t}(\sigma,\sigma^{\prime},\hat{u},\hat{u}^{\prime})$ of having a
connected pair $(\sigma,\sigma^{\prime})$ with $\hat{u}$ and
$\hat{u}^{\prime}$ unsatisfied interactions with their other neighbors,
respectively. In this case, $\hat{u},\hat{u}^{\prime}$ are integers between
zero and $K-1$. The magnitude
$\hat{P}^{t}(\hat{u}\mid\sigma,\sigma^{\prime},\hat{u}^{\prime})$ that we have
in (IV.2) is the conditional probability density:

$\hat{P}^{t}(\hat{u}\mid\sigma,\sigma^{\prime},\hat{u}^{\prime})=\frac{\hat{P}^{t}(\sigma,\sigma^{\prime},\hat{u},\hat{u}^{\prime})}{\sum_{u^{\prime}}\hat{P}^{t}(\sigma,\sigma^{\prime},\hat{u},u^{\prime})}$
(91)

The closure relation used in DINA can be derived starting from the assumptions
of the dynamical replica approach under replica symmetry Semerjian_dyn_2004 .
This explains why DINA gives very accurate results in models where detailed
balance holds. In this case, it is more likely that the microscopic
probability distribution function is a constant within a subspace with a
finite number of order parameters. Indeed, we know that in those models the
condition is satisfied at least at equilibrium.

Figure 7: Dynamics of the Ising ferromagnet according to Kinetic Monte Carlo
simulations (points), the integration of the DINA (thin-dashed lines) and
equations in the levels _CME-1_ (thick lines),
_CME-2_ (dashed lines) and CME-3 (dashed-dotted lines). In all cases an
initially fully-magnetized system evolves in time in contact with a heat bath
at a given temperature. From top to bottom, we have $T=1.5,1.75,1.82,2.0$.
KMC’s calculations were done for a random regular graph of size $N=10000$ with
connectivity $K=3$. The main panel and the inserted graphic show the time
evolution of system’s magnetization (see equation (83)). The first one
compares the different levels of approximation of our dynamic cavity method
with the DINA, and the latter compares the DINA with the results of Kinetic
Monte Carlo simulations.

As can be seen, (IV.2) is an average over the graphs ensemble of the
differential equation for $P^{t}(\sigma_{i},\sigma_{\partial i})$, which has
the form (53). Therefore, closure relations (88) and (89) are equivalent to
our closure relation (70) in the sense that both are approximations for the
same kind of conditional probabilities. Actually, they both drop the
information about the number of unsatisfied relations $\hat{u}^{\prime}$. An
important difference is in the fact that the first ones are written for the
average case probability densities, while the latter is written for the
single-instance.

However, for our dynamic cavity method all sites become equivalent at random
regular graphs with homogeneous initial conditions. In that case, all
probability densities defined over connected sets with the same structure
follow the same differential equations. We can drop information about the
local structure of the graph. As we show at Appendix B, we can forget the
irrelevant details of, let’s say, $P^{t}(\sigma_{i},\sigma_{\partial i})$. All
the combinations of $\sigma_{i}$, $\sigma_{\partial i}$ with the same value of
$\sigma_{i}$ and the same number of unsatisfied interactions between this spin
and its neighborhood will have the same probability density
$P^{t}(\sigma=\sigma_{i},u=\sum_{k\in\partial
i}\delta_{\sigma_{k},-\sigma_{i}})$.

We are easily able to write average case equations (see Appendix B) for
several levels of approximation, as the density (90) can be directly written
as $\hat{P}^{t}(\sigma,u)=\binom{K}{u}P(\sigma,u)$. Here we present one of the
equations we used to get the time dependence of the densities at (90):

$\displaystyle\frac{d\hat{P}^{t}(\sigma,u)}{dt}$ $\displaystyle=$
$\displaystyle-r(K,u)\,\hat{P}^{t}(\sigma,u)+r(K,K-u)\,\hat{P}^{t}(-\sigma,K-u)-$
(92)
$\displaystyle-(K-u)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{P}^{t}(\sigma,u)+(u+1)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{P}^{t}(\sigma,u+1)-$
$\displaystyle-u\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{P}^{t}(\sigma,u)+(K-u+1)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{P}^{t}(\sigma,\hat{u}-1)$

where we have assumed that $1\leq u\leq K-1$ and used the averages
$\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma^{\prime}}=\sum_{u^{\prime}=0}^{K-1}r(K,u^{\prime})\,\hat{p}^{t}(u^{\prime}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})$.
The equations for $u=0,K$ are easily obtained from (92) by suppressing its
second or third line, respectively.

Equation (92) is analogous to (IV.2), not only because they both give the time
derivative of the same quantities, but also because they have similar
structures. The main difference is that we have made the approximation:

$\hat{P}^{t}(\sigma,u)\,\hat{P}^{t}(u^{\prime}\mid\sigma^{\prime},\sigma,u)\approx\hat{P}^{t}(\sigma,u)\,\hat{p}^{t}(u^{\prime}\mid\sigma^{\prime},\sigma)\approx\hat{P}^{t}(\sigma,u)\,\hat{p}^{t}(u^{\prime}\mid\sigma^{\prime}\>\big{|}\big{|}\>\sigma)$
(93)

The closure (93) works as a substitution drawn from our dynamic cavity method
for the conditional probabilities
$\hat{P}^{t}(u^{\prime}\mid\sigma^{\prime},\sigma,u)$. The densities
$\hat{p}^{t}(\hat{u}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})$ can be
obtained from our cavity treatment (see Appendix B). On the other hand,
closures (88) and (89) of DINA are the substitutions for
$\hat{P}^{t}(u^{\prime}\mid\sigma^{\prime},\sigma,u)$ that correspond to the
assumptions of DRT.

Main panel of figure (7) shows a comparison between approximations in the
first, second and third levels of our dynamic cavity method, and the DINA.
Each level is defined by the approximation we use to obtain the cavity
probability densities
$\hat{p}^{t}(\hat{u}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})$. As the
level increases, our method predicts a time dependence of the magnetization
which is closer to what is obtained from DINA. This means that our cavity
closure on the densities
$\hat{p}^{t}(\hat{u}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})$ becomes
similar to the DRT-like closures (88) and (89), which work very well in this
case.

###  IV.3 Susceptible-infectious-susceptible model for epidemics

The propagation of an epidemic is a relevant issue present in various
scenarios. Problems like the dissemination of a disease within a population
pastor2015epidemic  or the spreading of a computer virus or rumors over a
network pastor2001epidemic  are studied by numerous scientists all over the
world. We do not have to explain the significance that global COVID-19
pandemics has brought nowadays to the development of theoretical tools for
understanding epidemic outbreaks cuevas2020agent .

There are a vast variety of models for epidemic processes. Since the seminal
work by Kermack and McKendric kermack1927contribution  about the susceptible-
infectious-recovered (SIR) model, we are used to seeing compartment models
whose main idea is to divide the population into several groups. Each one of
those groups is assumed to be homogeneous, in the sense that all individuals
interact with rules that depend on which group they belong to. However, this
does not necessarily mean that all of the members of a compartment are
equivalent. The dynamics, for example, can be defined over a specific contact
network. Here, we will concentrate on the susceptible-infectious-susceptible
(SIS) model defined over random regular graphs.

Figure 8: Allowed transitions in SIS compartment model on networks.

The standard SIS uses two compartments (states)
$x_{i}=0\equiv\mbox{susceptible}$, or $x_{i}=1\equiv\mbox{infectious}$ and is
the simplest model for recurrent transmissible diseases. The epidemic is thus
a continuous-time stochastic process with only two admitted transitions. An
infectious node transmits the disease to each one of its neighbors with rate
$\beta$, and recovers with rate $\mu$, as represented in Figure 8. The ratio
$\lambda=\beta/\mu$, known as spreading rate, is commonly used as the control
parameter of the model.

Such stochastic process can be described by a master equation like (27), where
the relevant quantities are probability densities over discrete variables:
$P^{t}(\vec{x})$. This allows us to apply our dynamic cavity method here, and,
as said before, we only have to adapt the rates:

$r_{i}(x_{i},x_{\partial i})=\beta\,\delta_{x_{i},0}\sum_{k\in\partial
i}x_{k}+\mu\,\delta_{x_{i},1}$ (94)

The probability per time unit that $x_{i}$ changes its value,
$r_{i}(x_{i},x_{\partial i})$, is equal to $\mu$ when node $i$ is infectious
($x_{i}=1$), and is equal to $\beta$ times the number of infected neighbors of
$i$ when the node is occupied by a susceptible individual ($x_{i}=0$). This
corresponds with our description of the model.

The dynamics of SIS on random networks has motivated abundant literature. One
of the most successful theoretical approaches in this scenario is known as
quenched mean-field (QMF) theory vanmieghem2009QMF ; castellano2010epid ,
which allows calculating epidemics threshold by complementing the master
equation with some suitable factorization. The most accurate version of QMF is
Pair Quenched Mean-Field (PQMF), or pair-based mean-field. This approximation
considers pair correlations and has been intensively investigated lately
silva19PQMF ; silva20PQMF  with very good results.

Figure 9: Dynamics of the SIS model for epidemics according to Kinetic Monte
Carlo simulations (points), the integration of PQMF (dashed lines) and
equations in the first (CME-1, thick lines) and second (CME-2, thin lines)
levels of our dynamic cavity method. Calculations were done for $\beta=0.4$
and $\mu=0.25,0.4,0.5,0.6$ (that in the figure correspond to
$\lambda=1.6,1.0,0.8,0.67$). In all cases the infection started in one node of
a random regular graph of size $N=1000$ with connectivity $K=3$. Top panels
panel: Time evolution of system’s average
$\rho(t)=\frac{1}{N}\sum_{i}\rho_{i}(t)$. Bottom panels: Time evolution of the
local error (computed analogously to 84) between these semi-analytical methods
and Kinetic Monte Carlo simulations.

The PQMF equations read:

$\displaystyle\frac{d\rho_{i}}{dt}$ $\displaystyle=$
$\displaystyle-\mu\rho_{i}+\beta\sum_{j\in\partial i}(\rho_{j}-\psi_{ij})$
(95)

and

$\displaystyle\frac{d\psi_{ij}}{dt}$ $\displaystyle=$
$\displaystyle-2(\mu+\beta)\psi_{ij}+\beta(\rho_{i}+\rho_{j})+$ (96)
$\displaystyle+\beta\frac{\rho_{j}-\psi_{ij}}{1-\rho_{i}}\sum_{k\in\partial
i\setminus
j}(\rho_{k}-\psi_{ik})+\beta\frac{\rho_{i}-\psi_{ij}}{1-\rho_{j}}\sum_{k\in\partial
j\setminus i}(\rho_{k}-\psi_{jk})$

where $\rho_{i}(t)=\sum_{x_{i}}x_{i}\,P^{t}(x_{i})$ and
$\psi_{ij}(t)=\sum_{x_{i}}\sum_{x_{j}}x_{i}\,x_{j}\,P^{t}(x_{i},x_{j})$

An extensive comparison with equations in the first level of approximation of
our dynamic cavity method (see (63) and (III.4.1)) has been already carried
out at CME-Epid-1 . As it is shown there, PQMF gives more accurate predictions
for the system’s observables. Left panels of figure (9) illustrate this for
the particular case of an epidemic outbreak that begins with a single
infectious individual in a random regular graph of size $N=1000$ and
connectivity $K=3$. Although the steady-state seems to be similarly predicted,
PQMF performs indisputably better at the transient regime than the equations
in CME-1. The corresponding maximum local errors, depicted in the bottom-left
panel, are approximately four times bigger for the latter than for PQMF.

Right panels of the figure (9) compares Kinetic Monte Carlo simulations with
PQMF’s results and the numeric integration of our dynamic cavity equations at
the second level of approximation. Equations (III.4.1)-(70) give a very good
description of the simulated epidemic outbreak. The bottom-right panel shows
that its maximum local errors are of the same order as PQMF’s, and for
$\lambda\leq 1.0$ is clear that our cavity method gives even preciser
predictions for the steady-state probability densities. This is also
observable in the top-right panel of the figure, which illustrates the time
dependence of the average of $\rho_{i}$ over all sites. There we can also see
that PQMF does not capture the transient regime appreciably better than our
new approach.

Thus, the second level of our dynamic cavity method outperforms PQMF at least
in describing the propagation of epidemics through a random regular graph.
Nevertheless, we expect the same to happen in other families of random
networks, like Erdos-Renyi graphs. As PQMF is currently the state-of-the-art
among mean-field approximations for the SIS dynamics silva19PQMF ; silva20PQMF
, this is a very relevant result.

##  V Conclusions

We derived, using the Theory of Random Point Processes, a hierarchical scheme
of Cavity Master Equations directly applicable to the continuous-time dynamics
of systems with discrete variables. We described carefully the approximations
made at each level of the scheme and explained their significance, pointing to
the system’s correlation length as the parameter defining the accuracy of the
approximation. Our scheme clarifies some of the contents and approximations
made on the recent literature about what it has been called Cavity Master
Equation CME-PRE ; CME-Pspin ; CME-PRL ; CME-AvAndBP .

We performed numerical tests at different levels of the approximation in spin
and epidemic models. For spin models, the accuracy of the technique is
comparable with the ones obtained with well-known successful methods like
Dynamic Independent Neighbor Approximation and Dynamical Replica Theory, and
improves as we increase the level of approximation. Studying the susceptible-
infectious-susceptible model for epidemics we show that our equations give
better predictions for the stationary state of the epidemics than the widely
used Pair Quenched Mean-Field Approximation, and, depending on the parameters,
can also perform better in the transient regime.

## Acknowledgements

This work was supported by the project _Mathematical Modelling of Epidemics,
from PNCB CITMA, Republic of Cuba_

## References

  * [1] E. Aurell, G. Del Ferraro, E. Domínguez, and R. Mulet.  A cavity master equation for the continuous time dynamics of discrete spins models.  Physical Review E, 95:052119, 2017.
  * [2] NG van Kampen.  Stochastic Processes in Physics and Chemistry, volume 1.  Elsevier, Amsterdam, 1992.
  * [3] D. J. Daley and D. Vere-Jones.  An Introduction to the Theory of Point Processes Volume I: Elementary Theory and Methods.  Springer Science & Business Media, Berlin, 2002.
  * [4] B. Karrer and M. E. J. Newman.  Message passing approach for general epidemic models.  Physical Review E, 82:016101, 2010.
  * [5] Munik Shrestha, Samuel V Scarpino, and Cristopher Moore.  Message-passing approach for recurrent-state epidemic models on networks.  Physical Review E, 92(2):022821, 2015.
  * [6] E. Aurell, E. Domínguez, D. Machado, and R. Mulet.  Exploring the diluted ferromagnetic p-spin model with a cavity master equation.  Physical Review E, 97:05103(R), 2018.
  * [7] E. Aurell, E. Domínguez, D. Machado, and R. Mulet.  A theory of non-equilibrium local search on random satisfaction problems.  Physical Review Letters, 123:230602, 2019.
  * [8] Roy J Glauber.  Time-dependent statistics of the ising model.  Journal of mathematical physics, 4:294, 1963.
  * [9] Anthony CC Coolen, Reimer Kühn, and Peter Sollich.  Theory of neural information processing systems.  Oxford University Press, 2005.
  * [10] A. Crisanti and H. Sompolinsky.  Dynamics of spin systems with randomly asymmetric bonds: Ising spins and glanber dynamics.  Physical Review A, 37:4865, 1988.
  * [11] H. Rieger, M. Schreckenberg, and J. Zittartz.  Glauber dynamics of the asymmetric sk-model.  Z. Phys. B - Condensed Matter, 74:527–538, 1989.
  * [12] Bernard Derrida, Elizabeth Gardner, and Anne Zippelius.  An exactly solvable asymmetric neural network model.  EPL (Europhysics Letters), 4(2):167, 1987.
  * [13] E. Domínguez, D. Machado, and R. Mulet.  The cavity master equation: average and fixed point of the ferromagnetic model in random graphs.  Journal of Statistical Mechanics: Theory and Experiment, 2020:073304, 2020.
  * [14] ACC Coolen, SN Laughton, and D Sherrington.  Dynamical replica theory for disordered spin systems.  Physical Review B, 53:8184, 1996.
  * [15] SN Laughton, ACC Coolen, and D Sherrington.  Order-parameter flow in the sk spin-glass: Ii. inclusion of microscopic memory effects.  Journal of Physics A: Mathematical and General, 29(4):763, 1996\.
  * [16] Alexander Mozeika and ACC Coolen.  Dynamical replica analysis of processes on finitely connected random graphs: I. vertex covering.  Journal of Physics A: Mathematical and Theoretical, 41(11):115003, 2008.
  * [17] Alexander Mozeika and ACC Coolen.  Dynamical replica analysis of processes on finitely connected random graphs: Ii. dynamics in the griffiths phase of the diluted ising ferromagnet.  Journal of Physics A: Mathematical and Theoretical, 42(19):195006, 2009.
  * [18] Hidetoshi Nishimori and Michiko Yamana.  Dynamical probability distribution function of the sk model at high temperatures.  Journal of the Physical Society of Japan, 65(1):3–6, 1996.
  * [19] ACC Coolen, SN Laughton, and D Sherrington.  Perturbative large deviation analysis of non-equilibrium dynamics.  Journal of the Physical Society of Japan, 93:084001, 2014.
  * [20] Wolfgang Barthel, Alexander K. Hartmann, and Martin Weigt.  Solving satisfiability problems by fluctuations: The dynamics of stochastic local search algorithms.  Phys. Rev. E, 67:066104, 2003.
  * [21] Guilhem Semerjian and Rémi Monasson.  Relaxation and metastability in a local search procedure for the random satisfiability problem.  Phys. Rev. E, 67:066103, 2003.
  * [22] G. Semerjian and M. Weigt.  Approximation schemes for the dynamics of diluted spin models: the ising ferromagnet on a bethe lattice.  Journal of Physics A: Mathematical and General, 37(21):5525–5546, may 2004.
  * [23] C. Kipnis and C. Landim.  Scaling limits of interacting particle systems, volume 320.  Springer Verlag, 1999.
  * [24] Paul Erdős and Alfréd Rényi.  On random graphs.  Publicationes Mathematicae Debrecen, 6:290–297, 1959.
  * [25] Eric Cator and Piet Van Mieghem.  Second-order mean-field susceptible-infected-susceptible epidemic threshold.  Physical review E, 85(5):056111, 2012.
  * [26] Angélica S Mata and Silvio C Ferreira.  Pair quenched mean-field theory for the susceptible-infected-susceptible model on complex networks.  EPL (Europhysics Letters), 103(4):48003, 2013.
  * [27] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani.  Epidemic processes in complex networks.  Reviews of modern physics, 87(3):925, 2015.
  * [28] Romualdo Pastor-Satorras and Alessandro Vespignani.  Epidemic spreading in scale-free networks.  Physical review letters, 86(14):3200, 2001.
  * [29] Erik Cuevas.  An agent-based model to evaluate the covid-19 transmission risks in facilities.  Computers in biology and medicine, 121:103827, 2020.
  * [30] William Ogilvy Kermack and Anderson G McKendrick.  A contribution to the mathematical theory of epidemics.  Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character, 115(772):700–721, 1927.
  * [31] Piet Van Mieghem, Jasmina Omic, and Robert Kooij.  Virus spread in networks.  IEEE/ACM Transactions on Networking, 17(1):1–14, 2009.
  * [32] Claudio Castellano and Romualdo Pastor-Satorras.  Thresholds for epidemic spreading in networks.  Phys. Rev. Lett., 105:218701, Nov 2010.
  * [33] Diogo H. Silva, Silvio C. Ferreira, Wesley Cota, Romualdo Pastor-Satorras, and Claudio Castellano.  Spectral properties and the accuracy of mean-field approaches for epidemics on correlated power-law networks.  Phys. Rev. Research, 1:033024, Oct 2019.
  * [34] Diogo H. Silva, Francisco A. Rodrigues, and Silvio C. Ferreira.  High prevalence regimes in the pair-quenched mean-field theory for the susceptible-infected-susceptible model on networks.  Phys. Rev. E, 102:012313, Jul 2020.
  * [35] Ernesto Ortega, David Machado, and Alejandro Lage-Castellanos.  Dynamics of epidemic models from cavity master equations.  arXiv preprint arXiv:2006.15881v2 [physics.soc-ph], 2021.

##  Appendix A Connecting cavity messages with instantaneous magnitudes

For completeness, let’s reproduce the derivation of equation (46) that appears
in [1].

If we marginalize (30) on $X_{\partial i\setminus\\{i,j\\}}$ and combine the
result with (32) we get a cavity message passing equation:

$\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid X_{j})=\sum_{\\{X_{k}\\},k\in\partial
i\setminus j}\Phi_{i}^{t}(X_{i}|X_{\partial i})\prod_{k\in\partial i\setminus
j}\mu_{k\rightarrow(ki)}^{t}(X_{k}\mid X_{i})$ (97)

Where $X_{i}(t)$ is the history of spin $i$ to time $t$. To simplify the
notation we will sometimes write $\mu_{i\rightarrow(ij)}^{t}$ for the cavity
conditional message.

We have learned that probability densities in the random point processes
formulation can be written as (see equation (33)):

$\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid
X_{j})=\prod_{l_{i}=1}^{s_{i}}\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t_{l_{i}})\exp\\{-\int_{t_{0}}^{t}\lambda_{i\rightarrow(ij)}(X_{i},X_{j},\tau)d\tau\\}$
(98)

Changing indexes accordingly, we can use the same parametrization for the
other cavity messages in the update equation.

The interaction term $\Phi_{i}^{t}(X_{i}|X_{\partial i})$ can be interpreted
as the probability density of $X_{i}$ conditioned on the histories of spins in
$\partial i$:

$\Phi_{i}^{t}(X_{i}|X_{\partial
i})=\prod_{l_{i}=1}^{s_{i}}r_{i}(\sigma_{i}(t_{l_{i}}),\sigma_{\partial
i}(t_{l_{i}}))\exp\\{-\int_{t_{0}}^{t}r_{i}(\sigma_{i}(\tau),\sigma_{\partial
i}(\tau))d\tau\\}$ (99)

Here, $r_{i}$ is the real jumping rate of $i$. Under a Markov assumption this
is an instantaneous quantity, meaning that at time $\tau$ it depends only on
the values of $\sigma_{i}(\tau),\sigma_{\partial i}(\tau)$.

The trace on the right hand side of (97) can be written in more detail. Let
$F$ be the argument of the sum:

$\sum_{\\{X_{k}\\},k\in\partial i\setminus j}^{[t_{0},t]}F(X_{i},X_{\partial
i},t)=\sum_{\\{s_{k}\\},k\in\partial i\setminus
j}\left[\prod_{k=1}^{d}\int_{t_{0}}^{t}dt_{1}^{k}\int_{t_{1}^{k}}^{t}dt_{2}^{k}\ldots\int_{t_{s_{k}-1}^{k}}^{t}dt_{s_{k}}^{k}\right]F(X_{i},X_{\partial
i},t)$ (100)

In all equations, $s_{k}$ will be the number of jumps of the history of spin
$k$.

If we write (97) for every pair $(i,j)$ in the network we get a system of
coupled equations. Having these functions, we could describe all the dynamics
of the system. However (97) is a very involved expression and we need to
transform it to a more tractable one. With that objective, we will
differentiate both sides of (97).

Differentiation in this context should be handled carefully since increasing
$t$ means we are changing the sample space itself. Therefore, it is safer to
use the definition of differentiation for both sides of equation (97). We will
compute the limit:

$\lim_{\Delta t\rightarrow 0}\dfrac{\mu_{i\rightarrow(ij)}^{t+\Delta
t}(X_{i}\mid X_{j})-\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid X_{j})}{\Delta t}$
(101)

A very important question arises at this point. What is the relation of
$(X_{i}(t+\Delta t),X_{j}(t+\Delta t))$ and $(X_{i}(t),X_{j}(t))$? Or in other
words, what happens in the interval $(t,t+\Delta t)$? The answer is important
because expressions for $\mu_{i\rightarrow(ij)}^{t+\Delta t}$ are different
whether we consider there can be jumps in the small $\Delta t$ interval or
not. The first thing that makes sense to impose is that histories must agree
up to time $t$. In $(t,t+\Delta t)$ we can have several combinations.

An implicit assumption throughout all this theory is that on an infinitesimal
interval only two things can happen to a spin; it can stay on its current
state or make one and only one jump to the opposite orientation. Two or more
jumps are not allowed. Considering this we have four cases to analyze:

  * •

There are $s_{i},s_{j}$ jumps in $(t_{0},t)$ and neither $i$ nor $j$ jumps in
$(t,t+\Delta t)$. This occurs with a probability $(1-\lambda_{i}\Delta
t)(1-\lambda_{j}\Delta t)$.

  * •

There are $s_{i},s_{j}$ jumps in $(t_{0},t)$ and $i$ XOR $j$ jumps in
$(t,t+\Delta t)$. This occurs with a probability $(1-\lambda_{i}\Delta
t)(\lambda_{j}\Delta t)$ or $(1-\lambda_{j}\Delta t)(\lambda_{i}\Delta t)$.
These are two cases in one.

  * •

There are $s_{i},s_{j}$ jumps in $(t_{0},t)$ and both $i$ and $j$ jumps in
$(t,t+\Delta t)$. This has a probability of $\lambda_{j}\lambda_{i}\Delta
t^{2}$

When $\Delta t$ goes to zero, from the previous analysis we conclude that the
derivative should be computed, with probability 1, using the first option,
where histories for $i$ and $j$ have no jumps in the interval of length
$\Delta t$.

To differentiate the left hand side of (97) we can use the parametrization
(98):

$\displaystyle\mu_{i\rightarrow(ij)}^{t+\Delta t}$ $\displaystyle=$
$\displaystyle\prod_{l_{i}=1}^{s_{i}}\lambda_{i\rightarrow(ij)}(t_{l_{i}})\exp\\{-\int_{t_{0}}^{t+\Delta
t}\lambda_{i\rightarrow(ij)}(\tau)d\tau\\}$
$\displaystyle\mu_{i\rightarrow(ij)}^{t+\Delta t}$ $\displaystyle=$
$\displaystyle[1-\lambda_{i\rightarrow(ij)}(t)\Delta
t]\prod_{l_{i}=1}^{s_{i}}\lambda_{i\rightarrow(ij)}(t_{l_{i}})\exp\\{-\int_{t_{0}}^{t}\lambda_{i\rightarrow(ij)}(\tau)d\tau\\}+o(\Delta
t)$ $\displaystyle\mu_{i\rightarrow(ij)}^{t+\Delta t}$ $\displaystyle=$
$\displaystyle[1-\lambda_{i\rightarrow(ij)}(t)\Delta
t]\>\mu_{i\rightarrow(ij)}^{t}+o(\Delta t)$ (102)

Then, with probability 1, the time derivative of the left hand side of
equation (97) is equal to
$-\lambda_{i\rightarrow(ij)}(t)\>\mu_{i\rightarrow(ij)}^{t}$

To calculate the derivative of the right hand side of (97) we should compute:

$\lim_{\Delta t\rightarrow 0}\dfrac{\sum_{\\{X_{k}\\},k\in\partial i\setminus
j}^{[t_{0},t+\Delta t]}F(X_{i},X_{\partial i},t+\Delta
t)-\mu_{i\rightarrow(ij)}^{t}}{\Delta t}$ (103)

Let’s focus on the first term of the numerator in the previous expression. It
can be expanded to the first order in $\Delta t$. It is important to remember
that $\Delta t$ appears in the integration limits as well as in the integrand
$F$. In addition, we should keep in mind that all jumps for $X_{i}$ and
$X_{j}$ must occur before $t$. This restriction, however, does not apply to
the histories $X_{k}$ for $k$ in $\partial i\setminus j$.

The expansion of (103) can be explained as follows. First, let’s remember that
$F$ is the joint probability of $X_{i}$ and $\\{X_{k}\\}$ with $k\in\partial
i\setminus j$, conditioned on $X_{j}$. All histories in the term of interest
are in the interval $[t_{0},t+\Delta t]$. The expression:

$\displaystyle\sum_{\\{X_{k}\\},k\in\partial i\setminus j}^{[t_{0},t+\Delta
t]}F(X_{i},X_{\partial i},t+\Delta t)$ (104)

is the marginalization of the mentioned joint probability distribution. The
previous sum, to order $\Delta t$, has two main contributions. One comes from
summing over $\\{X_{k}\\}$ with all $X_{k}$ having no jumps in $[t,t+\Delta
t]$:

$A=\sum_{\\{X_{k}\\},k\in\partial i\setminus j}^{[t_{0},t]}F(X_{i},X_{\partial
i},t)\Big{\\{}1-\Big{[}\sum_{k}\lambda_{k\rightarrow(ki)}(t)+r_{i}(t)\Big{]}\Delta
t\Big{\\}}$ (105)

The other considers all the possibilities of having one of the $X_{k}$ with a
jump in the interval of length $\Delta t$:

$B=\sum_{k}\displaystyle\sum_{\\{X_{k}\\},k\in\partial i\setminus
j}^{[t_{0},t]}F(X_{i},X_{\partial i},t)\lambda_{k\rightarrow(ki)}(t)\Delta t$
(106)

Then:

$\displaystyle\sum_{\\{X_{k}\\},k\in\partial i\setminus j}^{[t_{0},t+\Delta
t]}F(X_{i},X_{\partial i},t+\Delta t)=A+B+o(\Delta t)$ (107)

Putting all together we see that B cancels out with the $\lambda$ part of A,
and the remaining term of order 1 is
$\mu_{i\rightarrow(ij)}(X_{i}(t)|X_{j}(t))$, which cancels when inserted in
the limit expression.

The final outcome of this differentiation process is:

$\lambda_{i\rightarrow(ij)}[X_{i},X_{j},t]\;\;\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid
X_{j})=\displaystyle\sum_{\\{X_{k}\\},k\in\partial i\setminus
j}^{[t_{0},t]}r_{i}[\sigma_{i}(t),\sigma_{j}(t),\sigma_{\partial i\setminus
j}(t)]F(X_{i},X_{\partial i},t)$ (108)

We can now marginalize the right-hand side of the above equation over all the
histories of the spins $k\in\partial i\setminus j$ by keeping the
configuration of these spins at the last time $t$ fixed. The result reads

$\lambda_{i\rightarrow(ij)}(X_{i},X_{j},t)\;\mu_{i\rightarrow(ij)}^{t}(X_{i}\mid
X_{j})=\displaystyle\sum_{\sigma_{\partial i\setminus
j}(t)}\\!\\!r_{i}(\sigma_{i}(t),\sigma_{j}(t),\sigma_{\partial i\setminus
j}(t))\,p^{t}(\sigma_{\partial i\setminus j},X_{i}|X_{j})$ (109)

where we introduced the function $p$ as the marginalization of the function
$F$ over all the spin histories of the neighbors of $i$ except $j$, with the
configuration at the final time fixed. Note that the notation
$\sigma_{\partial i\setminus j}(t)$ is equivalent to
$\\{\sigma_{k}(t)\\}k\in\partial i\setminus j$ and that in $p$ above appears
again explicitly the conditional nature of the probability distribution $F$.

Equation (109) represents the differential dynamic message-passing update
equation obtained by differentiating (97) in time. It connects the derivative
of the dynamic message $\mu_{i\rightarrow(ij)}^{t}$, and so the jumping rate
$\lambda_{i\rightarrow(ij)}(t)$ of spin $i$ in the cavity used to parametrize
the message in (98), with the transition rate of the same spin
$r_{i}(\sigma_{i}(t),\sigma_{j}(t),\sigma_{\partial i\setminus j}(t))$ at time
$t$.

##  Appendix B Average dynamic cavity equations in random regular graphs

Whenever we have initial conditions independent of the site, and homogeneous
node connectivity, the equations of our dynamic cavity method acquire a
spatial symmetry that will allow us to reduce them to a few average case
equations. Let’s show this through an example.

We will start with equation (III.4.1).

$\displaystyle\frac{dp^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})}{dt}$ $\displaystyle=$
$\displaystyle-r_{i}(\sigma_{i},\sigma_{\partial i})\,p^{t}(\sigma_{\partial
i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})+r_{i}(-\sigma_{i},\sigma_{\partial
i})\,p^{t}(\sigma_{\partial i\setminus
j},-\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ $\displaystyle-\sum_{l\in\partial
i\setminus j}\sum_{\sigma_{\partial l\setminus
i}}\Big{\\{}r_{l}(\sigma_{l},\sigma_{\partial l})\,p^{t}(\sigma_{\partial
l\setminus
i}\mid\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})\,p^{t}(\sigma_{\partial
i\setminus j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{l}(-\sigma_{l},\sigma_{\partial
l})\,p^{t}(\sigma_{\partial l\setminus
i}\mid-\sigma_{l}\>\big{|}\big{|}\>\sigma_{i})\,p^{t}(F_{l}[\sigma_{\partial
i\setminus j}],\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\Big{\\}}$

In a random regular graph with connectivity $K$, a spatial symmetry results
from choosing an initial condition for $p^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$ which does not depends on $i$ and
$j$. In that case, the equation governing the time evolution of
$p^{t}(\sigma_{\partial i\setminus j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$
wont depend on the values of $i,j$ for any time. We can write a single
equation for all sites:

$\displaystyle\frac{dp^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})}{dt}$
$\displaystyle=$
$\displaystyle-r(K,\hat{u}+\delta_{\sigma,-\sigma^{\prime}})\,p^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})+r(K,K-\hat{u}-\delta_{\sigma,-\sigma^{\prime}})\,p^{t}(K-\hat{u},-\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
(111)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-(K-1-\hat{u})\sum_{u^{\prime}=0}^{K-1}\binom{K-1}{u^{\prime}}r(K,u^{\prime})\,\Big{\\{}p^{t}(u^{\prime}\mid\sigma\>\big{|}\big{|}\sigma)\,p^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-p^{t}(u^{\prime}\mid-\sigma\>\big{|}\big{|}\sigma)\,p^{t}(\hat{u}+1,\sigma\>\big{|}\big{|}\>\sigma^{\prime})\Big{\\}}-$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\hat{u}\sum_{u^{\prime}=0}^{K-1}\binom{K-1}{u^{\prime}}r(K,u^{\prime})\,\Big{\\{}p^{t}(u^{\prime}\mid-\sigma\>\big{|}\big{|}\sigma)\,p^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-p^{t}(u^{\prime}\mid\sigma\>\big{|}\big{|}\sigma)\,p^{t}(\hat{u}-1,\sigma\>\big{|}\big{|}\>\sigma^{\prime})\Big{\\}}$

In (111), regardless of the site where they are defined, the probability
densities depend only on two connected spin variables, $\sigma$ and
$\sigma^{\prime}$, and on the number $\hat{u}=0,\ldots,K-1$ of unsatisfied
interactions that $\sigma$ has with its neighbors other than
$\sigma^{\prime}$. The sum in the fourth and fifth lines takes into account
that $\hat{u}$ unsatisfied interactions, all of which contribute equally to
the derivative. The same happens with second and third lines and the
contribution of the remaining $(K-1-\hat{u})$ satisfied interactions.

So far, we have not defined new probability densities, we just re-denoted
$p^{t}(\sigma_{\partial i\setminus j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})$
into $p^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})$ after realizing
that the values of sites $i$ and $j$ where irrelevant due to spatial symmetry,
and that it was not important to keep track of all the combinations of
$\sigma_{\partial i\setminus j}$, we needed only to record the number of
unsatisfied $\hat{u}$ interactions between $\sigma_{\partial i\setminus j}$
and $\sigma_{i}$. Now, we will introduce the total densities
$\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})=\binom{K-1}{\hat{u}}p^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})$,
noticing that these are analogous to the average probability densities in
(90). Actually, exactly as before, we have:

$\displaystyle\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})=\lim_{N\rightarrow\infty}\ll\frac{1}{NK}\sum_{i=1}^{N}\sum_{j\in\partial
i}\sum_{\sigma_{i}}\sum_{\sigma_{j}}\sum_{\sigma_{\partial i}\setminus
j}p_{\xi_{K}(N)}^{t}(\sigma_{\partial i\setminus
j},\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\times$
$\displaystyle\times\,\delta_{\sigma,\sigma_{i}}\delta_{\sigma^{\prime},\sigma_{j}}\>\>\delta_{(\sum_{k\in\partial
i\setminus j}\sigma_{k}),K-1-2\hat{u}}\gg_{\xi_{K}(N)}$ (112)

where $\xi_{K}(N)$ is the ensemble of random regular graphs with connectivity
$K$ and size $N$. The symbol $\ll\cdot\gg_{\xi_{K}(N)}$ represents an average
over this ensemble, with proper probabilistic weights.

The differential equation for this average probability densities is:

$\displaystyle\frac{d\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})}{dt}$
$\displaystyle=$
$\displaystyle-r(K,\hat{u}+\delta_{\sigma,-\sigma^{\prime}})\,\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})+r(K,K-\hat{u}-\delta_{\sigma,-\sigma^{\prime}})\,\hat{p}^{t}(K-\hat{u},-\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
(113)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-(K-1-\hat{u})\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})+(\hat{u}+1)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{p}^{t}(\hat{u}+1,\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\hat{u}\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})+(K-\hat{u})\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{p}^{t}(\hat{u}-1,\sigma\>\big{|}\big{|}\>\sigma^{\prime})$

where we have written the averages
$\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma^{\prime}}=\sum_{u^{\prime}=0}^{K-1}r(K,u^{\prime})\,\hat{p}^{t}(u^{\prime}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})$.
The conditional cavity probability densities in (113) are:

$\hat{p}^{t}(\hat{u}\mid\sigma\>\big{|}\big{|}\>\sigma^{\prime})=\frac{\hat{p}^{t}(\hat{u},\sigma\>\big{|}\big{|}\>\sigma^{\prime})}{\sum_{u^{\prime}}\hat{p}^{t}(u^{\prime},\sigma\>\big{|}\big{|}\>\sigma^{\prime})}$
(114)

Equations (113) and (114) form a closed system that can be numerically
integrated. This must be complemented with:

$\displaystyle\frac{d\hat{P}^{t}(\sigma,u)}{dt}$ $\displaystyle=$
$\displaystyle-r(K,u)\,\hat{P}^{t}(\sigma,u)+r(K,K-u)\,\hat{P}^{t}(-\sigma,K-u)-$
(115)
$\displaystyle-(K-u)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{P}^{t}(\sigma,u)+(u+1)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{P}^{t}(\sigma,u+1)-$
$\displaystyle-u\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{-\sigma,\sigma}\,\hat{P}^{t}(\sigma,u)+(K-u+1)\>\big{\langle}\,r(K,u^{\prime})\,\big{\rangle}_{\sigma,\sigma}\,\hat{P}^{t}(\sigma,\hat{u}-1)$

which can be derived similarly.

To finish writing the equations of the levels of approximation shown in main
text’s figure (7), we still need to address first and third levels: CME-1 and
CME-3. In the first case we start with equations (63) and (III.4.1)

$\displaystyle\frac{dp^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})}{dt}$
$\displaystyle=$ $\displaystyle-\sum_{\sigma_{\partial i\setminus
j}}\Big{\\{}r_{i}(\sigma_{i},\sigma_{\partial i})\Big{[}\prod_{k\in\partial
i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}\>\sigma_{i})\Big{]}p^{t}(\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})-$
(116)
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i}(-\sigma_{i},\sigma_{\partial
i})\Big{[}\prod_{k\in\partial i\setminus
j}p^{t}(\sigma_{k}\>\big{|}\big{|}-\sigma_{i})\Big{]}p^{t}(-\sigma_{i}\>\big{|}\big{|}\>\sigma_{j})\Big{\\}}$
$\displaystyle\frac{dP^{t}(\sigma_{i})}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{\sigma_{\partial
i}}\Big{\\{}r_{i}(\sigma_{i},\sigma_{\partial i})\,\Big{[}\prod_{k\in\partial
i}p^{t}(\sigma_{k}\>\big{|}\big{|}\sigma_{i})\Big{]}\,P^{t}(\sigma_{i})-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r_{i}(-\sigma_{i},\sigma_{\partial
i})\,\Big{[}\prod_{k\in\partial
i}p^{t}(\sigma_{k}\>\big{|}\big{|}-\sigma_{i})\Big{]}\,P^{t}(-\sigma_{i})\Big{\\}}$

It is easy so see that the average case versions of this equations at random
regular graphs with homogeneous initial conditions are:

$\displaystyle\frac{d\hat{p}^{t}(\sigma\>\big{|}\big{|}\>\sigma^{\prime})}{dt}$
$\displaystyle=$
$\displaystyle-\sum_{u^{\prime}=0}^{K-1}\Big{\\{}r(K,u^{\prime}+\delta_{\sigma,-\sigma^{\prime}})\big{[}\hat{p}(-\sigma\mid\sigma)\big{]}^{u^{\prime}}\,\big{[}\hat{p}(\sigma\mid\sigma)\big{]}^{K-1-u^{\prime}}\hat{p}^{t}(\sigma\>\big{|}\big{|}\>\sigma^{\prime})-$
(118)
$\displaystyle\>\>\>\>\>\>-r(K,u^{\prime}+\delta_{\sigma,\sigma^{\prime}})\big{[}\hat{p}(\sigma\mid-\sigma)\big{]}^{u^{\prime}}\,\big{[}\hat{p}(-\sigma\mid-\sigma)\big{]}^{K-1-u^{\prime}}\hat{p}^{t}(-\sigma\>\big{|}\big{|}\>\sigma^{\prime})\Big{\\}}$
$\displaystyle\frac{d\hat{P}^{t}(\sigma)}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{u^{\prime}=0}^{K}\Big{\\{}r(K,u^{\prime})\,\big{[}\hat{p}(-\sigma\mid\sigma)\big{]}^{u^{\prime}}\,\big{[}\hat{p}(\sigma\mid\sigma)\big{]}^{K-u^{\prime}}\,\hat{P}^{t}(\sigma)-$
$\displaystyle\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>-r(K,u^{\prime})\,\big{[}\hat{p}(\sigma\mid-\sigma)\big{]}^{u^{\prime}}\,\big{[}\hat{p}(-\sigma\mid-\sigma)\big{]}^{K-u^{\prime}}\,\hat{P}^{t}(-\sigma)\Big{\\}}$

Equations (118) and (B) also form a closed system and its numerical
integration is shown in figure (7) of main text. The equations in the third
level of approximation, CME-3, whose results appear also in that figure,
require some extra work.

Our probability densities are now
$p^{t}(\hat{u}_{2},\sigma_{2},\tilde{u}_{1},\sigma_{1}\>\big{|}\big{|}\>\sigma_{0})$,
which are defined over the connected sets illustrated in bottom-left panel of
(4). The spins $\sigma_{0}$, $\sigma_{1}$ and $\sigma_{2}$ are marked in
figure (10). The integer $\hat{u}_{2}$, which goes from zero to $K-1$,
represents the number of unsatisfied interactions between $\sigma_{2}$ and its
neighbors, without counting $\sigma_{1}$. Finally, the integer
$\tilde{u}_{1}$, which goes from zero to $K-2$, is the number of unsatisfied
interactions between $\sigma_{1}$ and its neighbors, without counting
$\sigma_{0}$ and $\sigma_{2}$.

Figure 10: The cavity probability density
$p^{t}(\hat{u}_{2},\sigma_{2},\tilde{u}_{1},\sigma_{1}\>\big{|}\big{|}\>\sigma_{0})$
is defined over a connected set in CS-3. This figure illustrates meaning of
each one of the variables $\hat{u}_{2}$, $\sigma_{2}$, $\tilde{u}_{1}$,
$\sigma_{1}$ and $\sigma_{0}$.

As we learned with (58), the exact equations for this densities are:

$\displaystyle\frac{dp^{t}(\hat{u}_{2},\sigma_{2},\tilde{u}_{1},\sigma_{1}\>\big{|}\big{|}\>\sigma_{0})}{dt}$
$\displaystyle=$