text:
# The Blume-Emery-Griffiths neural network: dynamics for arbitrary temperature

## 1 Introduction

Recently, we have studied the Blume-Emery-Griffiths (BEG) neural network dynamics at zero temperature using the probabilistic signal-to-noise analysis (SNA) . The interest of such a study arises from the fact that the BEG-model optimises the mutual information content for three-state networks ,. In this way the best retrieval properties including, e.g., the largest retrieval overlap, loading capacity, basin of attraction and convergence time, are guaranteed in comparison with other three-state networks studied in the literature . As has been pointed out in  the system allows for new stable stationary states, i.e. quadrupolar states, when the temperature becomes high enough. Therefore, it is interesting to study the dynamics of the model at arbitrary temperatures. This can be done in the framework of the SNA by introducing auxiliary thermal fields  to express the stochastic dynamics within the gain function formulation of the deterministic dynamics. The calculations are rather involved because the feedback in these systems  causes the appearance of discrete noise, besides the Gaussian one, involving the neurons at all previous time steps. This prevents a closed-form solution of the dynamics although it allows for a recursive scheme to find the order parameters, starting from the time evolution of the distribution of the local fields. Hence, it seems worthwhile to also apply the generating functional approach (GFA) ,  to solve this feedback dynamics. This method enables us to find all relevant physical order parameters at any time step via the derivatives of the generating functional.

Both methods are then compared for the first few time steps and the results are analytically shown to be identical up to the third time step. Numerical simulations are performed for systems up to $`7000`$ neurons and they confirm the analytical results we have obtained. Beyond the third time step it is found numerically that the SNA approach certainly leads to different results than those obtained through simulations for values of the model parameters corresponding to spin-glass behaviour of the system. This is traced back to the treatment of the residual overlap representing the influence of the non-condensed patterns in the time evolution of the system. Finally, as in the $`T=0`$ case, we derive fixed-point equations in the stationary limit and compare them with existing results within a replica symmetric thermodynamic approach .

The paper is organized as follows. In Section 2 we introduce the BEG model and the order parameters of interest. Section 3 is devoted to the solution of the dynamics of the model with the SNA. The first few time steps are obtained for the main overlap, the neural activity, the activity overlap and the variance of the residual overlaps. Moreover, the stationary limit of the dynamics is discussed. Furthermore, numerical simulations of this dynamics are presented for the first three time steps showing agreement with the theoretical findings. Section 4 discusses the GFA approach to the BEG dynamics. In Section 5 we compare both methods analytically up to the third time step and numerically beyond that. Some concluding remarks are given in Section 6. All technical details are refered to Appendices.

## 2 The BEG model

Consider a neural network consisting of $`N`$ neurons which can take values $`\{\sigma _i|i=1,\mathrm{},N\}`$ from the discrete set $`𝒮=\{s_k|k=1,2,3\}\{1,0,+1\}`$. The $`p=\alpha N`$ patterns to be stored in this network are supposed to be a collection of independent and identically distributed random variables (i.i.d.r.v.), $`\{\xi _i^\mu |\mu =1,\mathrm{},p;i=1,\mathrm{},N\}`$ with a probability distribution

$$p(\xi _i^\mu )=\frac{a}{2}\delta (\xi _i^\mu 1)+\frac{a}{2}\delta (\xi _i^\mu +1)+(1a)\delta (\xi _i^\mu )$$
(1)
with $`a`$ the activity of the patterns so that

$$\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{i}{\overset{N}{}}(\xi _i^\mu )^2=a.$$
(2)
Given the network configuration at time $`t`$, $`𝝈_N(t)\{\sigma _j(t)|j=1,\mathrm{},N\}`$, the neurons are updated according to the stochastic parallel spin-flip dynamics defined by the transition probabilities

$$\text{Prob}\left(\sigma _i(t+1)=s_k𝒮|𝝈_N(t)\right)=\frac{\mathrm{exp}\{\beta ϵ_i[s_k|𝝈_N(t)]\}}{_{s𝒮}\mathrm{exp}\{\beta ϵ_i[s|𝝈_N(t)]\}}.$$
(3)
The configuration $`𝝈_N(0)`$ is chosen as input. The energy potential $`ϵ_i[s|𝝈_N(t)]`$ is defined by

$$ϵ_i[s|𝝈_N(t)]=sh_i(𝝈_N(t))s^2\theta _i(𝝈_N(t)),$$
(4)
where the following local fields in neuron $`i`$ carry all the information

$$h_i(𝝈_N(t))=\underset{ji}{}J_{ij}\sigma _j(t),\theta _i(𝝈_N(t))=\underset{ji}{}K_{ij}\sigma _j^2(t).$$
(5)
The synaptic couplings $`J_{ij}`$ and $`K_{ij}`$ are of the Hebb-type

$$J_{ij}=\frac{1}{a^2N}\underset{\mu =1}{\overset{p}{}}\xi _i^\mu \xi _j^\mu ,K_{ij}=\frac{1}{N}\underset{\mu =1}{\overset{p}{}}\eta _i^\mu \eta _j^\mu $$
(6)
with

$$\eta _i^\mu =\frac{1}{\stackrel{~}{a}}((\xi _i^\mu )^2a),\stackrel{~}{a}=a(1a).$$
(7)
The order parameters relevant for the present discussion of this model are

$$m_N^\mu (t)=\frac{1}{aN}\underset{i}{}\xi _i^\mu \sigma _i(t),l_N^\mu (t)=\frac{1}{N}\underset{i}{}\eta _i^\mu \sigma _i^2(t),q_N(t)=\frac{1}{N}\underset{i}{}\sigma _i^2(t)$$
(8)
where $`m_N^\mu `$ is the retrieval overlap, $`l_N^\mu (t)`$ related to the activity overlap, and $`q_N(t)`$ the neural activity , .

## 3 The SNA for arbitrary temperatures

### 3.1 Recursive dynamical scheme

In  we discussed the SNA for the BEG model at temperature $`T=0`$. In order to generalize this approach to finite temperatures we introduce auxiliary thermal fields  in order to express the stochastic dynamics within the gain function formulation of the deterministic dynamics. Consider for each $`i`$ and $`t`$ the auxiliary thermal fields, $`𝚿_i(t)=\{\psi _i^l(t)|l=1,2,3\}`$, which are i.i.d.r.v. with respect to $`i`$ and $`t`$, such that

$$\begin{array}{c}\text{Prob}[\sigma _i(t+1)=s_k𝒮|𝝈_N(t)]=\hfill \\ \hfill \left[\underset{lk}{}\mathrm{\Theta }(ϵ_i[s_l|𝝈_N(t)]\psi _i^l(t+1)ϵ_i[s_k|𝝈_N(t)]+\psi _i^k(t+1))\right]_\mathrm{\Psi }\end{array}$$
(9)
with $`[]_\mathrm{\Psi }`$ the expectation value with respect to $`𝚿_i(t)`$. The joint probability density of the $`𝚿_i(t)`$ with respect to $`l`$ follows from comparing (3)-(4) with (9) and the usual normalisation to $`1`$. The result is

$$f_\beta [𝚿_i(t)]=\beta ^23!\frac{\mathrm{exp}[\beta _{l=1}^3\psi _i^l(t)]}{\{_{l=1}^3\mathrm{exp}[\beta \psi _i^l(t)]\}^3}\delta \left(\underset{l=1}{\overset{3}{}}\psi _i^l(t)\right).$$
(10)
To make the link between the stochastic dynamics and the deterministic dynamics in a convenient way we observe that for each realization of the $`𝚿_i(t+1)`$ with density $`f_\beta ()`$ the network state evolves according to the deterministic rule

$`\sigma _i(t+1)`$ $`={\displaystyle \underset{k=1}{\overset{3}{}}}s_k{\displaystyle \underset{lk}{}}\mathrm{\Theta }(ϵ_i[s_l|𝝈_N(t)]\psi _i^l(t+1)ϵ_i[s_k|𝝈_N(t)]+\psi _i^k(t+1))`$
$`g_\mathrm{\Psi }[h_i(𝝈_N(t)),\theta _i(𝝈_N(t))].`$ (11)
In this way the problem of deriving recursion relations for $`T0`$ becomes tractable. We start from the following initial conditions

$`E[\sigma _i(0)]=0,Var[\sigma _i(0)]=q_0,`$ (12)
$`E[\xi _i^\mu \sigma _j(0)]=\delta _{i,j}\delta _{\mu ,1}m_0^1a,m_0^1>0,E[\eta _i^\mu \sigma _j^2(0)]=\delta _{i,j}\delta _{\mu ,1}l_0^1.`$ (13)
So we assume that the initial configuration is correlated with only one (i.e. the condensed) stored pattern. By the law of large numbers (LLN) Eqs.(8) and (12)-(13) determine the order parameters $`m_N^1(0),q_N^1(0)`$ and $`l_N^1(0)`$ at $`t=0`$ in the limit $`N\mathrm{}`$. In the following we leave out the pattern index $`1`$.

For a general time step we find from Eq.(11) and the LLN in the limit $`N\mathrm{}`$ for the order parameters (eqs.(8))

$`m^1(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`{\displaystyle \frac{1}{a}}\left[\xi _ig_𝚿(h_i(t),\theta _i(t))\right]_\mathrm{\Psi }`$ (14)
$`q(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`\left[g_𝚿^2(h_i(t),\theta _i(t))\right]_\mathrm{\Psi }`$ (15)
$`l^1(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`\left[\eta _ig_𝚿^2(h_i(t),\theta _i(t))\right]_\mathrm{\Psi }`$ (16)
where $`h_i(t)=lim_N\mathrm{}h_i(𝝈_N(t))`$ (with an analogous formula for $`\theta _i(t)`$), and where the convergence is in probability. In the above the average, denoted by $``$, has to be taken over both the distribution of the $`\{\xi _i^\mu \}`$ and the $`\{\sigma _i(0)\}`$. Note that the $`\{\sigma _i(0)\}`$ is hidden in the local fields through the updating rule (11). For the Eqs. (14)-(16) the average over the thermal fields can be done explicitly by using the relations (9)-(10) leading to

$`m^1(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`{\displaystyle \frac{1}{a}}\xi _iV_\beta (h_i(t),\theta _i(t))`$ (17)
$`q(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`W_\beta (h_i(t),\theta _i(t))`$ (18)
$`l^1(t+1)`$ $`\stackrel{𝒫r}{=}`$ $`\eta _iW_\beta (h_i(t),\theta _i(t))`$ (19)
with

$`V_\beta (h_i(t),\theta _i(t))`$ $`={\displaystyle \frac{\mathrm{sinh}(\beta h_i(t))}{\frac{1}{2}\mathrm{exp}\left(\beta \theta _i(t)\right)+\mathrm{cosh}(\beta h_i(t))}}`$ (20)
$`W_\beta (h_i(t),\theta _i(t))`$ $`={\displaystyle \frac{\mathrm{cosh}(\beta h_i(t))}{\frac{1}{2}\mathrm{exp}\left(\beta \theta _i(t)\right)+\mathrm{cosh}(\beta h_i(t))}}.`$ (21)
In order to obtain the recursion relations for the local fields we can then follow the derivation of the zero temperature case  (see Section 3) step by step. We do not repeat this calculation here but write down the final results

$`h_i(t+1)`$ $`={\displaystyle \frac{\xi _i}{a}}m(t+1)+\chi _h(t)\left\{h_i(t){\displaystyle \frac{1}{a}}\xi _im(t)+{\displaystyle \frac{\alpha }{a}}\sigma _i(t)\right\}+𝒩(0,{\displaystyle \frac{\alpha }{a^2}}q(t+1))`$ (22)
$`\theta _i(t+1)`$ $`=\eta _il(t+1)+\chi _\theta (t)\left\{\theta _i(t)\eta _il(t)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\sigma _i^2(t)\right\}+𝒩(0,{\displaystyle \frac{\alpha }{\stackrel{~}{a}^2}}q(t+1)).`$ (23)
From this it is clear that the local fields at time $`t+1`$ consist out of a discrete part and a normally distributed part

$`h_i(t)`$ $`=`$ $`M_i(t)+𝒩(0,\alpha aD(t))`$ (24)
$`\theta _i(t)`$ $`=`$ $`L_i(t)+𝒩(0,{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}E(t))`$ (25)
where

$`M_i(t+1)`$ $`=`$ $`\chi _h(t)\left[M_i(t){\displaystyle \frac{\xi _i}{a}}m(t)+{\displaystyle \frac{\alpha }{a}}\sigma _i(t)\right]+{\displaystyle \frac{\xi _i}{a}}m(t+1)`$ (26)
$`L_i(t+1)`$ $`=`$ $`\chi _\theta (t)\left[L_i(t)\eta _il(t)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\sigma _i^2(t)\right]+\eta _il(t+1)`$ (27)
with $`D(t)`$ and $`E(t)`$ the variances of the residual overlaps defined by

$$r^\mu (t)=\underset{N\mathrm{}}{lim}\frac{1}{a^2\sqrt{N}}\underset{j}{}\xi _j^\mu \sigma _j(t),s^\mu (t)=\underset{N\mathrm{}}{lim}\frac{1}{\sqrt{N}}\underset{j}{}\eta _j^\mu \sigma _j^2(t),\mu >1,$$
(28)
and taking into account the noise produced by the non-condensed patterns in the local fields. They satisfy the recursion relations

$`D(t+1)`$ $`=`$ $`{\displaystyle \frac{q(t+1)}{a^3}}+\chi _h^2(t)D(t)+2\chi _h(t)\text{Cov}[\stackrel{~}{r}^\mu (t),r^\mu (t)]`$ (29)
$`E(t+1)`$ $`=`$ $`{\displaystyle \frac{q(t+1)}{\stackrel{~}{a}}}+\chi _\theta ^2(t)E(t)+2\chi _\theta (t)\text{Cov}[\stackrel{~}{s}^\mu (t),s^\mu (t)].`$ (30)
In these expressions, the modified residual overlaps $`\stackrel{~}{r}^\mu ,\stackrel{~}{s}^\mu ,\mu >1`$ are introduced because in the residual overlaps the local fields strongly depend upon $`\xi _i^\mu `$ and $`\eta _i^\mu `$ respectively. Therefore, the local fields for finite $`N`$ have been modified as follows

$$\stackrel{~}{h}_{N,i}^\mu (t)=h_{N,i}(t)\frac{1}{\sqrt{N}}\xi _i^\mu r_N^\mu (t),\stackrel{~}{\theta }_{N,i}^\mu (t)=\theta _{N,i}(t)\frac{1}{\sqrt{N}}\eta _i^\mu s_N^\mu (t),$$
(31)
such that, employing the central limit theorem (CLT)

$`\stackrel{~}{r}^\mu (t)`$ $`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{a^2\sqrt{N}}}{\displaystyle \underset{j}{}}\xi _j^\mu g_\mathrm{\Psi }[\stackrel{~}{h}_{N,j}(t),\stackrel{~}{\theta }_{N,j}(t)]\stackrel{𝒫r}{=}𝒩(0,{\displaystyle \frac{q(t+1)}{a^3}})`$ (32)
$`\stackrel{~}{s}^\mu (t)`$ $`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{j}{}}\eta _j^\mu g_\mathrm{\Psi }^2[\stackrel{~}{h}_{N,j}(t),\stackrel{~}{\theta }_{N,j}(t)]\stackrel{𝒫r}{=}𝒩(0,{\displaystyle \frac{q(t+1)}{\stackrel{~}{a}}}),`$ (33)
where we have now assumed that the $`\stackrel{~}{r}^\mu ,\stackrel{~}{s}^\mu ,\mu >1`$ depend only weakly on $`\xi _i^\mu `$ and $`\eta _i^\mu `$ respectively. We remark that in the thermodynamic limit the density distributions of the modified local fields $`\stackrel{~}{h}_i(t)`$ and $`\stackrel{~}{\theta }_i(t)`$ equal the density distributions of the local fields $`h_i(t)`$ and $`\theta _i(t)`$ itself. The $`\chi _h(t)`$ and $`\chi _\theta (t)`$ are the susceptibilities corresponding to the fields $`h_i(t)`$ and $`\theta _i(t)`$

$`\chi _h(t)`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle 𝑑\widehat{h}𝑑\widehat{\theta }\rho _{h_i(t)}(\widehat{h})\rho _{\theta _i(t)}(\widehat{\theta })\left(\frac{V_\beta (h(t),\theta (t))}{h(t)}\right)}|_{(\widehat{h},\widehat{\theta })}`$ (34)
$`\chi _\theta (t)`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}{\displaystyle 𝑑\widehat{h}𝑑\widehat{\theta }\rho _{h_i(t)}(\widehat{h})\rho _{\theta _i(t)}(\widehat{\theta })\left(\frac{W_\beta (h(t),\theta (t))}{\theta (t)}\right)}|_{(\widehat{h},\widehat{\theta })}.`$ (35)
In these expressions, $`\rho _{h_i(t)}(x)`$ and $`\rho _{\theta _i(t)}(x)`$ are the probability densities of the local fields

$$\begin{array}{cc}\hfill \rho _{h_i(t)}(z)& =\left(\underset{s=0}{\overset{t2}{}}dx_sdy_s\right)𝑑x_t𝑑y_t\frac{1}{\sqrt{det(2\pi 𝐂_h)}}\frac{1}{\sqrt{det(2\pi 𝐂_\theta )}}\hfill \\ & \times \mathrm{exp}\left(\frac{1}{2}\mathrm{𝐱𝐂}_h^1𝐱^{}\frac{1}{2}\mathrm{𝐲𝐂}_\theta ^1𝐲^{}\right)\delta \left(zM_i(t)\sqrt{\alpha aD(t)}x_t\right)\hfill \end{array}$$
(36)
$$\begin{array}{cc}\hfill \rho _{\theta _i(t)}(z)& =\left(\underset{s=0}{\overset{t2}{}}dx_sdy_s\right)𝑑x_t𝑑y_t\frac{1}{\sqrt{det(2\pi 𝐂_h)}}\frac{1}{\sqrt{det(2\pi 𝐂_\theta )}}\hfill \\ & \times \mathrm{exp}\left(\frac{1}{2}\mathrm{𝐱𝐂}_h^1𝐱^{}\frac{1}{2}\mathrm{𝐲𝐂}_\theta ^1𝐲^{}\right)\delta \left(zL_i(t)\sqrt{\frac{\alpha }{\stackrel{~}{a}}E(t)}y_t\right)\hfill \end{array}$$
(37)
with $`𝐱=(x_0,\mathrm{},x_{t2},x_t)`$ and $`𝐲=(y_0,\mathrm{},y_{t2},y_t)`$ two sets of correlated normally distributed variables which we choose to normalize. The correlations matrices between the different time steps are then given by the elements $`𝐂_h(t,t^{})=E[x_tx_t^{}]`$ and $`𝐂_\theta (t,t^{})=E[y_ty_t^{}]`$.

This concludes the explanation of the recursive scheme from which the order parameters of the system can be obtained. A practical difficulty that remains is the calculation of the explicit correlations at different time steps. As an illustration we present the first few time steps in Appendix A.

### 3.2 Stationary equations

Stationary equations for non-zero temperature can, in principle, be obtained from the zero temperature expressions obtained in  by introducing auxiliary thermal fields and averaging over them. However, extreme care has to be taken concerning these thermal fluctuations. The point is that in the case of the stationary limit at zero temperature it is well known that there is no difference beween the neural activity and the Edwards-Anderson order parameter, but for non-zero temperatures they become distinct. Therefore, in order to obtain the stationary limit, we assume that for $`t\mathrm{}`$ the local fields can be replaced by their thermal averages

$$\underset{t\mathrm{}}{lim}[\sigma _i^n(t+1)]_\mathrm{\Psi }=\underset{t\mathrm{}}{lim}[g_\mathrm{\Psi }^n(h_i(t),\theta _i(t))]_\mathrm{\Psi }\underset{t\mathrm{}}{lim}[g_\mathrm{\Psi }^n([h_i(t)]_\mathrm{\Psi },[\theta _i(t)]_\mathrm{\Psi })]_\mathrm{\Psi }n=1,2.$$
(38)
The first step is to find the thermally averaged fields. Starting from the residual overlaps (28) we write (we focus on the $`h`$-quantities)

$$[r^\mu (t+1)]_\mathrm{\Psi }=[\stackrel{~}{r}^\mu (t)]_\mathrm{\Psi }+[\chi _h(t)r^\mu (t)]_\mathrm{\Psi }$$
(39)
where now

$$[\stackrel{~}{r}^\mu (t)]_\psi =\underset{N\mathrm{}}{lim}\frac{1}{a^2\sqrt{N}}\underset{i}{}\xi _i^\mu [g_\mathrm{\Psi }(\stackrel{~}{h}_{N,i}^\mu (t),\stackrel{~}{\theta }_{N,i}^\mu (t))]_\mathrm{\Psi }\stackrel{𝒫r}{=}𝒩(0,\frac{1}{a^3}q_1(t+1))$$
(40)
with

$$q_1(t)=[\sigma _i(t)]_\mathrm{\Psi }^2.$$
(41)
The recursion relation for the field then becomes

$$[h_i(t+1)]_\mathrm{\Psi }=\frac{\xi _i}{a}m(t+1)+\chi _h(t)\left[h_i(t)\frac{\xi _i}{a}m(t)+\frac{\alpha }{a}\sigma _i(t)\right]_\mathrm{\Psi }+𝒩(0,\frac{q_1(t+1)}{a})$$
(42)
where we have used the LLN for $`m(t)`$ and $`\chi _h(t)`$, viz.

$$\begin{array}{cc}\hfill m(t)& =\frac{1}{aN}\underset{i}{}\xi _i\sigma _i(t)\stackrel{𝒫r}{=}\frac{1}{a}\xi _i[\sigma _i(t)]_\mathrm{\Psi }\hfill \\ \hfill \chi _h(t)& =\frac{1}{a^2N}\underset{i}{}(\xi _i^\mu )^2\frac{g_\mathrm{\Psi }(t)}{h_i(t)}\stackrel{𝒫r}{=}\frac{1}{a}\left[\frac{g_\mathrm{\Psi }(\mathrm{})}{h}\right]_\mathrm{\Psi }.\hfill \end{array}$$
(43)
Completely analogous equations can be written down for $`\theta _i(t+1)_\mathrm{\Psi }`$ and $`\stackrel{~}{s}^\mu (t+1)_\mathrm{\Psi }`$ with

$$p_1(t)=[\sigma _i^2(t)]_\mathrm{\Psi }^2.$$
(44)
Next, in order to obtain the stationary limit we assume that the thermal averages of the local fields and of the order parameters do not change in time. This leads to

$$\begin{array}{cc}\hfill [h]_\mathrm{\Psi }& =\frac{\xi }{a}m+𝒩(0,\frac{\alpha q_1}{a^2(1\chi _h)^2})+\frac{\alpha }{a}\eta _h[\sigma ]_\mathrm{\Psi }\stackrel{~}{h}+\frac{\alpha }{a}\eta _h[\sigma ]_\mathrm{\Psi }\hfill \\ \hfill [\theta ]_\mathrm{\Psi }& =\eta l+𝒩(0,\frac{\alpha p_1}{\stackrel{~}{a}^2(1\chi _\theta )^2})+\frac{\alpha }{\stackrel{~}{a}}\eta _\theta [\sigma ^2]_\mathrm{\Psi }\stackrel{~}{\theta }+\frac{\alpha }{\stackrel{~}{a}}\eta _\theta [\sigma ^2]_\mathrm{\Psi }\hfill \end{array}$$
(45)
with $`\eta _x=\chi _x/(1\chi _x)`$, $`x=h,\theta `$. The solution for $`[\sigma ]_\mathrm{\Psi }`$ is a self-consistent equation of the form

$$[\sigma ]_\mathrm{\Psi }=V_\beta (\stackrel{~}{h}+\frac{\alpha }{a}\eta _h[\sigma ]_\mathrm{\Psi },\stackrel{~}{\theta }+\frac{\alpha }{\stackrel{~}{a}}\eta _\theta [\sigma ^2]_\mathrm{\Psi }).$$
(46)
In order to solve this equation we consider the case $`T=0`$. There the averages over the thermal noise disappear, yielding

$$\underset{\beta \mathrm{}}{lim}V_\beta (h,\theta )=\text{sign}(h)\mathrm{\Theta }(|h|+\theta )$$
(47)
and a Maxwell construction  leads to

$$\sigma =\text{sign}(h)\mathrm{\Theta }(|h|+\theta +\mathrm{\Delta })$$
(48)
where

$$\mathrm{\Delta }=\frac{1}{2}\left(\frac{\alpha }{a}\eta _h+\frac{\alpha }{\stackrel{~}{a}}\eta _\theta \right).$$
(49)
We introduce then a set of auxiliary thermal fields and average over them to arrive at the fixed point equations in the limit $`N\mathrm{}`$ (dropping the index i)

$`m`$ $`=`$ $`{\displaystyle \frac{\xi }{a}}{\displaystyle DzDyV_\beta (h(z),\theta (y)+\mathrm{\Delta })}`$ (50)
$`l`$ $`=`$ $`\eta {\displaystyle DzDyW_\beta (h(z),\theta (y)+\mathrm{\Delta })}`$ (51)
$`q`$ $`=`$ $`{\displaystyle DzDyW_\beta (h(z),\theta (y)+\mathrm{\Delta })}`$ (52)
$`\chi _h`$ $`=`$ $`{\displaystyle \frac{1\chi _h}{\sqrt{\alpha q_1}}}{\displaystyle DzzDyV_\beta (h(z),\theta (y)+\mathrm{\Delta })}`$ (53)
$`\chi _\theta `$ $`=`$ $`{\displaystyle \frac{1\chi _\theta }{\sqrt{\alpha p_1}}}{\displaystyle DzDyyW_\beta (h(z),\theta (y)+\mathrm{\Delta })}`$ (54)
with

$$h(z)=\frac{\xi }{a}m+\frac{\sqrt{\alpha q_1}}{a(1\chi _h)}z,\theta (y)=\eta l+\frac{\sqrt{\alpha p_1}}{\stackrel{~}{a}(1\chi _\theta )}y.$$
(55)
and $`q_1`$ and $`p_1`$ given by

$$\begin{array}{cc}\hfill q_1& =DzDyV_\beta ^2(h(z),\theta (y)+\mathrm{\Delta })\hfill \\ \hfill p_1& =DzDyW_\beta ^2(h(z),\theta (y)+\mathrm{\Delta }).\hfill \end{array}$$
(56)
Furthermore, in agreement with the fluctuation-dissipation theorem

$$\chi _h=\frac{\beta }{a}(qq_1)\chi _\theta =\frac{\beta }{\stackrel{~}{a}}(qp_1)$$
(57)
with $`\chi _h`$ and $`\chi _\theta `$ the susceptibilities defined before, in the stationary limit. It is interesting to remark that the equations (50)-(54) are the same as the fixed-point equations derived from a replica symmetric thermodynamic approach with a hamiltonian for sequential dynamics .

### 3.3 Numerical results and simulations

The results derived in Sections 2-3 and apppendix A have been studied numerically and compared with simulations for systems up to $`N=7000`$ neurons averaged over $`500`$ runs.

As in the $`T=0`$ case the first few time steps, given by the explicit formula in appendix A, and the simulations agree very well over the whole range of $`\alpha `$. As an illustration we refer to Fig. 1 presenting the order parameters as a function of $`\alpha `$ for uniform patterns, $`m_0=l_0=0.6,q_0=0.5`$ and $`T=0.5`$. We remark that the critical capacity for this system is $`\alpha _c0.06`$ (). We then learn that the first time steps do give us a reasonable estimate for the critical capacity, especially through the order parameter $`l`$.

In Fig. 2 we examine the order parameters $`m`$ and $`l`$ in the quadrupolar phase ($`m=0,l>0`$) versus the paramagnetic phase ($`m=0,l=0`$), for several values of $`\alpha `$. We see that a few time steps do give us already the characteristic behaviour. When time increases $`m`$ decreases while $`l`$ differentiates between the phases, as is seen in the theoretical results as well as in the simulations. For the quadrupolar phase ($`\alpha =0.001`$) $`l`$ increases, deep inside the paramagnetic phase ($`\alpha =0.1`$) $`l`$ decreases, while in the intermediate region ($`\alpha =0.01`$) the rate of increase of $`l`$ quickly diminishes and $`l`$ itself goes to zero.

At this point, we remark that there is a small but visible discrepancy between the theory and simulations especially in $`l(3)`$. It is of the order $`O(10^3)`$ and attributed to finite-size effects. In this respect we report that earlier calculations  of further time steps of the overlap order parameter in the Hopfield model also showed a discrepancy between simulations and numerical results, e.g., for $`\alpha =0.10`$ the difference $`m_{th}(3)m_{sim}(3)`$ is again of the order $`O(10^3)`$. However, for $`0.1<m_0<0.4`$, the difference $`m_{th}(4)m_{sim}(4)`$ is about $`0.018`$ and always positive both for $`N=6000`$ and $`N=7000`$. This, and the fact that the SNA does not give a closed form solution of the dynamics, has motivated us to look in Section 4 at a second, alternative approach for solving the dynamics.

## 4 The generating functional approach

The second approach we use to discuss the dynamics of the BEG-model is an extension of the generating functional approach (GFA)   (for a recent review see, e.g.,  and references therein). It is an exact procedure based on the idea of looking at the probability for finding a microscopic path in time

$$\text{Prob}[𝝈_N(0),\mathrm{},𝝈_N(t)]=\text{Prob}(𝝈(0))\underset{t^{}=1}{\overset{t}{}}\underset{i=1}{\overset{N}{}}\text{Prob}\left(\sigma _i(t^{})=s_k𝒮|𝝈_N(t^{}1)\right)$$
(58)
where the transition probabilities are given by (3)-(4) with the local fields $`h_i(𝝈_N(t))`$ and $`\theta _i(𝝈_N(t))`$ replaced by

$$h_i(𝝈_N(t))=\underset{ji}{}J_{ij}\sigma _j(t)+\gamma _{h,i}(t)\theta _i(t)=\underset{ji}{}K_{ij}\sigma _j^2(t)+\gamma _{\theta ,i}(t).$$
(59)
Here the $`\gamma _{h,i}(t)`$ and $`\gamma _{\theta ,i}(t)`$ are time-dependent external fields introduced in order to define a response function.

The basic tool is then the generating functional

$`Z_N[\mathrm{\Psi },\mathrm{\Phi }]`$ $`=\mathrm{exp}\left\{i{\displaystyle \underset{s=1}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\psi _i(s)\sigma _i(s)+\varphi _i(s)\sigma _i^2(s)\right)\right\}_{paths}`$
$`={\displaystyle \underset{𝝈_N(0)}{}}\mathrm{}{\displaystyle \underset{𝝈_N(t)}{}}\text{Prob}[𝝈_N(0),\mathrm{},𝝈_N(t)]`$
$`\times \mathrm{exp}\left\{i{\displaystyle \underset{s=1}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\psi _i(s)\sigma _i(s)+\varphi _i(s)\sigma _i^2(s)\right)\right\}`$ (60)
from which we can find all the relevant observables by calculating derivatives with respect to $`\{\psi _i(t)\}`$ and $`\{\varphi _i(t)\}`$ and taking these zero afterwards. In the sequel we omit the subscript indicating the path average. This path average of the spin, the correlation functions and the response functions are given by

$$\begin{array}{cc}\sigma _i(t)=i\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\psi _i(t)},& \sigma _i^2(t)=i\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\varphi _i(t)},\\ C_{h,ij}(t,t^{})=\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\psi _i(t)\psi _j(t^{})},& C_{\theta ,ij}(t,t^{})=\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\varphi _i(t)\varphi _j(t^{})},\\ G_{h,ij}(t,t^{})=i\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\psi _i(t)\gamma _{h,j}(t^{})},& G_{\theta ,ij}(t,t^{})=i\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2Z_N[\mathrm{\Psi },\mathrm{\Phi }]}{\varphi _i(t)\gamma _{\theta ,j}(t^{})}.\end{array}$$
(61)
The generating function is averaged over all pattern realizations, i.e. over the disorder, before the path average is taken. In the limit $`N\mathrm{}`$, this results in a theory with effective single spin local fields which, for the model at hand, are given by

$`h(t)`$ $`={\displaystyle \frac{\xi }{a}}m(t)+\gamma _h(t)+{\displaystyle \frac{\alpha }{a}}{\displaystyle \underset{s^{}=0}{\overset{t}{}}}R_h(t,s^{})\sigma (s^{})+\sqrt{\alpha }z_h(t)`$ (62)
$`\theta (t)`$ $`=\eta l(t)+\gamma _\theta (t)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}{\displaystyle \underset{s^{}=0}{\overset{t}{}}}R_\theta (t,s^{})\sigma ^2(s^{})+\sqrt{\alpha }z_\theta (t)`$ (63)
where the noises $`z_h(t),z_\theta (t)`$ are a set of correlated normally distributed variables with the measure

$`w(𝐳_h,𝐳_\theta )`$ $`={\displaystyle \frac{\mathrm{exp}\left\{\frac{1}{2}_{s,s^{}=0}^{t1}z_h(s)[𝐰_h^1](s,s^{})z_h(s^{})\frac{1}{2}_{s,s^{}=0}^{t1}z_\theta (s)[𝐰_\theta ^1](s,s^{})z_\theta (s^{})\right\}}{\sqrt{(2\pi )^tdet(𝐰_h)}\sqrt{(2\pi )^tdet(𝐰_\theta )}}}`$ (64)
$`𝐰_h`$ $`=\left(a𝐈𝐆_h\right)^1𝐂_h\left(a𝐈𝐆_h^{}\right)^1`$ (65)
$`𝐰_\theta `$ $`=\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1𝐂_\theta \left(\stackrel{~}{a}𝐈𝐆_\theta ^{}\right)^1.`$ (66)
The retarded self-interactions appearing in (62)-(63) are then given by

$$𝐑_h=\left(a𝐈𝐆_h\right)^1𝐆_h,𝐑_\theta =\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1𝐆_\theta .$$
(67)
Furthermore, the order parameters can be written as

$`m(t)`$ $`=`$ $`i\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{aN}}{\displaystyle \underset{i}{}}\xi _i{\displaystyle \frac{\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(t)}}={\displaystyle \frac{1}{a}}\xi \sigma (t)`$ (68)
$`l(t)`$ $`=`$ $`i\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\eta _i{\displaystyle \frac{\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(t)}}=\eta \sigma ^2(t)`$ (69)
$`C_h(t,t^{})`$ $`=`$ $`\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \frac{^2\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(t)\psi _i(t^{})}}=\sigma (t)\sigma (t^{})`$ (70)
$`C_\theta (t,t^{})`$ $`=`$ $`\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \frac{^2\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(t)\varphi _i(t^{})}}=\sigma ^2(t)\sigma ^2(t^{})`$ (71)
$`G_h(t,t^{})`$ $`=`$ $`i\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \frac{^2\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(t)\gamma _{h,i}(t^{})}}={\displaystyle \frac{}{\gamma _h}}\sigma (t^{})`$ (72)
$`G_\theta (t,t^{})`$ $`=`$ $`i\underset{N\mathrm{}}{lim}\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \frac{^2\overline{Z_N[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(t)\gamma _{\theta ,i}(t^{})}}={\displaystyle \frac{}{\gamma _\theta }}\sigma ^2(t^{})`$ (73)
where the overline denotes the disorder average and $``$ denotes the averages over the noises $`𝐳_h,𝐳_\theta `$, and the initial conditions and the condensed pattern, according to the measure

$$f[\sigma (0),\mathrm{},\sigma (t)]=𝑑𝐳_h𝑑𝐳_\theta w[𝐳_h,𝐳_\theta ]\text{Prob}[\sigma (0),\mathrm{},\sigma (t)|𝐳_h,𝐳_\theta ]f[\sigma (0),\mathrm{},\sigma (t)]_\xi $$
(74)
Here

$$\begin{array}{c}\text{Prob}[\sigma (0),\mathrm{},\sigma (t)|𝐳_h,𝐳_\theta ]\hfill \\ \hfill =\underset{\sigma (0),\mathrm{},\sigma (t)}{}\text{Prob}(\sigma (0))\underset{s=0}{\overset{t1}{}}\frac{\mathrm{exp}\left(\beta \sigma (s+1)h(s)+\beta \sigma ^2(s+1)\theta (s)\right)}{2\mathrm{exp}\left(\beta \theta (s)\right)\mathrm{cosh}(\beta h(s))+1}\end{array}$$
(75)
with $`h(s)`$ and $`\theta (s)`$ given by (62) and (63), $`\gamma _h,\gamma _\theta =0`$. For more details of this calculation we refer to Appendix B.

For later convenience we remark that $`G_h(t,t^{})=0`$ for $`tt^{}`$ due to causality (analogously for $`G_\theta (t,t^{})`$), and that for $`t>t^{}`$

$`G_h(t,t^{})`$ $`=`$ $`\beta \sigma (t)\left[\sigma (t^{}+1)V_\beta (h(t^{}),\theta (t^{}))\right]`$ (76)
$`G_\theta (t,t^{})`$ $`=`$ $`\beta \sigma ^2(t)\left[\sigma ^2(t^{}+1)W_\beta (h(t^{}),\theta (t^{}))\right].`$ (77)
This set of equations completely solves the dynamics for the BEG-model. We calculate the first few time steps in Appendix C, in order to compare the results with those of the SNA method.

## 5 Comparison between SNA and GFA

By comparing the explicit calculations described in the Appendices A and C for the first two time steps in the dynamics of the BEG model we find that they are identical and that in the GFA response functions, restricting ourselves for the moment to the $`h`$-field,

$$G_h(1,0)=a\chi _h(0)G_h(2,1)=a\chi _h(1).$$
(78)
The third time step requires a bit more work. Looking at the response function in general (recall (76)) and considering $`t^{}=t2`$, the sums over $`\sigma (t)`$ and $`\sigma (t1)`$ can be done explicitly because of the specific time behaviour of the local fields, leading to

$$G_h(t,t2)=0.$$
(79)
Using this, together with the causality requirement on $`G_h`$ and noting that a similar treatment is valid for $`G_\theta `$, we find that the result for the third time step in the GFA method coincides with the corresponding SNA result. Comparing the next time steps explicitly is more cumbersome. Also from the local fields as expressed by Eqs.(22)-(23) in the SNA approach and Eqs.(62)-(63) in the GFA method no immediate obvious analytic connection is apparent. Therefore, a separate detailed analytic study of these relationships for the simpler Hopfield model is appropriate .

In the sequel we present some strong numerical evidence that the SNA and the GFA approach will lead to slightly different results for certain values of the model parameters, specifying spin-glass behaviour, from the fourth time step onwards. This is consistent with our remark at the end of Section 3 that a difference, although a very small one, has been noticed between the SNA approach and simulations at zero temperature in the fourth time step of the parallel Hopfield dynamics . This result was not entirely conclusive because of the possible entanglement with finite-size effects.

Reviewing the SNA approach one knows that one of the basic ingredients is the introduction of a modified residual overlap $`\stackrel{~}{r}(t)`$ which converges to a normal distribution in the thermodynamic limit (recall eqs. (32)-(33)).

We have done some numerical experiments for different values of the model parameters comparing this limiting normal distribution in the SNA method with simulations for different time steps. Some results are shown in Fig. 3. We compare the normal distribution with simulations for time steps $`t=2`$ and $`t=9`$ for systems with $`N=2000`$ neurons averaged over $`250`$ runs for the initial conditions $`m_0=l_0=0.6,q_0=0.5,a=2/3`$ and temperature $`T=0.2`$ as a function of $`\alpha `$. We note that the critical $`\alpha `$ for these parameter values is $`\alpha _c=0.086`$ .

We conclude that in the retrieval region ($`\alpha _c<0.086`$) the simulations results coincide quite well with the limiting distribution, while in the spin-glass region, certainly from $`\alpha 0.11`$ onwards, the results for $`t=9`$ start diverting systematically. We remark that a set of analogous results can be obtained for $`\stackrel{~}{s}(t)`$. This may indicate that for $`t>3`$ the SNA works quite well in a large part of the retrieval phase since most of the time the system reaches the attractor already after a few time steps. However, in the spin-glass phase the SNA method does not estimate correctly the noise of the non-condensed patterns represented by $`\stackrel{~}{r}(t)`$, suggesting that the assumption on the residual noise discussed above is no longer fulfilled. Resolving these issues analytically is not immediately straightforward and beyond the scope of the present work. They will be worked out first for the Hopfield model .

## 6 Concluding remarks

The parallel dynamics of the fully connected BEG neural network model for finite temperatures has been solved using both the signal to noise analysis and the generating functional approach. It is shown that both methods give identical results up to the first three time steps.

Furthermore, by studying the distribution of the residual noise in the signal to noise analysis numerically and comparing it with further time steps up to $`t=9`$ it is found that the signal to noise method as it is used in the literature does not estimate correctly the modified residual overlaps in the spin-glass region. However, in a large part of the retrieval region it can be considered as accurate. Numerical simulations confirm these findings. The signal to noise approach can be made exact by carefully studying the long time correlations .

## Acknowledgments

The authors thank A.C.C. Coolen for useful discussions. This work has been supported in part by the Fund of Scientific Research, Flanders-Belgium.

## Appendix A The first three time steps in the SNA

Employing the general recursive scheme of Section 3 evolution equations are derived for the first three time steps of the dynamics in the SNA approach. The following initial conditions are taken

$`m^1(0)`$ $`=`$ $`m_0,l^1(0)=l_0,q(0)=q_0,`$ (80)
$`h_i(0)`$ $`=`$ $`{\displaystyle \frac{\xi _i^1}{a}}m_0+𝒩(0,{\displaystyle \frac{\alpha q_0}{a^2}}),\theta _i(0)=\eta _i^1l_0+𝒩(0,{\displaystyle \frac{\alpha q_0}{\stackrel{~}{a}^2}}),`$ (81)
$`D(0)`$ $`=`$ $`{\displaystyle \frac{q_0}{a^3}},E(0)={\displaystyle \frac{q_0}{\stackrel{~}{a}}}.`$ (82)
From now we forget about the superscript $`1`$ to indicate the condensed pattern. Furthermore, we do not write down all details but we refer to  for the analogous $`T=0`$ case.

### A.1 First time step

From (17)-(19) we immediately find

$`m(1)`$ $`=`$ $`{\displaystyle \frac{1}{a}}\xi {\displaystyle DzDyV_\beta (h_0^{}(z),\theta _0^{}(y))},`$ (83)
$`q(1)`$ $`=`$ $`{\displaystyle DzDyW_\beta (h_0^{}(z),\theta _0^{}(y))},`$ (84)
$`l(1)`$ $`=`$ $`\eta {\displaystyle DzDyW_\beta (h_0^{}(z),\theta _0^{}(y))},`$ (85)
where

$$h_0^{}(z)=\frac{1}{a}\xi m_0+\sqrt{\alpha aD(0)}z,\theta _0^{}(y)=\eta l_0+\sqrt{\frac{\alpha }{\stackrel{~}{a}}E(0)}y$$
(86)
and where $`Dx=dx\mathrm{exp}(x^2/2)/\sqrt{2\pi }`$ is the Gaussian measure. Because of the modification of the local fields we know that

$`\text{Cov}[\stackrel{~}{r}^\mu (0),r^\mu (0)]=E\left[\stackrel{~}{r}^\mu (0)r^\mu (0)\right]`$ $`=`$ $`E\left[{\displaystyle \frac{1}{a^3}}\sigma _i(0)g_\mathrm{\Psi }(\stackrel{~}{h}_i^\mu (0),\stackrel{~}{\theta }_i^\mu (0))\right]`$ (87)
$`\text{Cov}[\stackrel{~}{s}^\mu (0),s^\mu (0)]=E\left[\stackrel{~}{s}^\mu (0)s^\mu (0)\right]`$ $`=`$ $`E\left[{\displaystyle \frac{1}{\stackrel{~}{a}}}\sigma _i^2(0)g_\mathrm{\Psi }^2(\stackrel{~}{h}_i^\mu (0),\stackrel{~}{\theta }_i^\mu (0))\right]`$ (88)
Using the recursion relations (29)-(30) we obtain

$`D(1)`$ $`=`$ $`{\displaystyle \frac{q(1)}{a^3}}+\chi _h^2(0)D(0)+2\chi _h(0)R(1,0)`$ (89)
$`E(1)`$ $`=`$ $`{\displaystyle \frac{q(1)}{\stackrel{~}{a}}}+\chi _\theta ^2(0)E(0)+2\chi _\theta (0)S(1,0),`$ (90)
where the correlation parameters $`R(t,t^{})`$ and $`S(t,t^{})`$ are defined by

$`R(t,t^{})`$ $`=`$ $`{\displaystyle \frac{1}{a^3}}E\left[g_\mathrm{\Psi }(\stackrel{~}{h}(t1),\stackrel{~}{\theta }(t1))g_\mathrm{\Psi }(\stackrel{~}{h}(t^{}1),\stackrel{~}{\theta }(t^{}1))\right]`$ (91)
$`R(t,0)`$ $`=`$ $`{\displaystyle \frac{1}{a^3}}E\left[g_\mathrm{\Psi }(\stackrel{~}{h}(t1),\stackrel{~}{\theta }(t1))\sigma (0)\right]`$ (92)
$`S(t,t^{})`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}E\left[g_\mathrm{\Psi }^2(\stackrel{~}{h}(t1),\stackrel{~}{\theta }(t1))g_\mathrm{\Psi }^2(\stackrel{~}{h}(t^{}1),\stackrel{~}{\theta }(t^{}1))\right]`$ (93)
$`S(t,0)`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}E\left[g_\mathrm{\Psi }^2(\stackrel{~}{h}(t1),\stackrel{~}{\theta }(t1))\sigma ^2(0)\right]`$ (94)
Since the density distribution of $`h`$ and $`\stackrel{~}{h}`$ are the same in the limit $`N\mathrm{}`$ we have for the first time step

$$R(1,0)=\frac{1}{a^3}\sigma (0)DzDyV_\beta (h_0^{}(z),\theta _0^{}(y))$$
(95)
$$S(1,0)=\frac{1}{\stackrel{~}{a}}\sigma ^2(0)DzDyW_\beta (h_0^{}(z),\theta _0^{}(y))$$
(96)
Furthermore, since there are no feedback correlations at time zero

$`\chi _h(0)`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle DzDy\left(\frac{V_\beta (h(t),\theta (t))}{h(t)}\right)}|_{(h_0^{}(z),\theta _0^{}(z))}`$ (97)
$`\chi _\theta (0)`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}{\displaystyle DzDy\left(\frac{W_\beta (h(t),\theta (t))}{\theta (t)}\right)}|_{(h_0^{}(z),\theta _0^{}(z))}`$ (98)
### A.2 Second time step

First, we need the distribution of the local fields at time step $`1`$ using (24)-(27)

$`h_i(1)`$ $`=`$ $`{\displaystyle \frac{\xi _i}{a}}m(1)+{\displaystyle \frac{\alpha }{a}}\chi _h(0)\sigma _i(0)+𝒩(0,\alpha aD(1))`$ (99)
$`\theta _i(1)`$ $`=`$ $`\eta _il(1)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\chi _\theta (0)\sigma _i^2(0)+𝒩(0,{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}E(1))`$ (100)
These results enable us to write down the important order parameters at time $`2`$:

$`m(2)`$ $`=`$ $`{\displaystyle \frac{1}{a}}\xi {\displaystyle DzDyV_\beta (h_1^{}(z),\theta _1^{}(y))}`$ (101)
$`q(2)`$ $`=`$ $`{\displaystyle DzDyW_\beta (h_1^{}(z),\theta _1^{}(y))}`$ (102)
$`l(2)`$ $`=`$ $`\eta {\displaystyle DzDyW_\beta (h_1^{}(z),\theta _1^{}(y))}`$ (103)
defining

$`h_1^{}(z)`$ $`=`$ $`{\displaystyle \frac{1}{a}}\xi m(1)+{\displaystyle \frac{\alpha }{a}}\chi _h(0)\sigma (0)+\sqrt{\alpha aD(1)}z`$ (104)
$`\theta _1^{}(y)`$ $`=`$ $`\eta l(1)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\chi _\theta (0)\sigma ^2(0)+\sqrt{{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}E(1)}y`$ (105)
The calculation of the variance of the residual overlaps proceeds as follows. From the recursion relations (29)-(30) we get

$`D(2)`$ $`=`$ $`{\displaystyle \frac{q(2)}{a^3}}+\chi _h^2(1)D(1)+2\chi _h(1)(R(2,1)+\chi _h(0)R(2,0))`$ (106)
$`E(2)`$ $`=`$ $`{\displaystyle \frac{q(2)}{\stackrel{~}{a}}}+\chi _\theta ^2(1)E(1)+2\chi _\theta (1)(S(2,1)+\chi _\theta (0)S(2,0))`$ (107)
where

$`R(2,0)`$ $`=`$ $`{\displaystyle \frac{1}{a^3}}\sigma (0){\displaystyle DzDyV_\beta (h_1^{}(z),\theta _1^{}(y))}`$ (108)
$`S(2,0)`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}\sigma ^2(0){\displaystyle DzDyW_\beta (h_1^{}(z),\theta _1^{}(y))}`$ (109)
To obtain the further correlations introduced we remark that because $`\sigma (1)`$ appears in the expressions, we have to introduce the correlation between time steps $`0`$ and $`1`$. The relevant correlation coefficients are given by

$`\rho _h(t,t^{})`$ $`=`$ $`{\displaystyle \frac{E[(h(t)M(t))(h(t^{})M(t^{}))]}{\sqrt{\alpha aD(t)}\sqrt{\alpha aD(t^{})}}}`$ (110)
$`\rho _\theta (t,t^{})`$ $`=`$ $`{\displaystyle \frac{E[(\theta (t)L(t))(\theta (t^{})L(t^{}))]}{\sqrt{(\alpha /\stackrel{~}{a})E(t)}\sqrt{(\alpha /\stackrel{~}{a})E(t^{})}}}`$ (111)
and the joint distribution

$$D\omega _x^{a,b}(z,y)=\frac{dzdy}{2\pi \sqrt{1\rho _x^2(a,b)}}\mathrm{exp}\left(\frac{z^22zy\rho _x(a,b)+y^2}{2(1\rho _x^2(a,b))}\right)$$
(112)
Then we arrive at

$`R(2,1)`$ $`={\displaystyle \frac{1}{a^3}}{\displaystyle D\omega _h^{1,0}(z,s)D\omega _\theta ^{1,0}(y,t)V_\beta (h_1^{}(z),\theta _1^{}(y))V_\beta (h_0^{}(s),\theta _0^{}(t))}`$ (113)
$`S(2,1)`$ $`={\displaystyle \frac{1}{\stackrel{~}{a}}}{\displaystyle D\omega _h^{1,0}(z,s)D\omega _\theta ^{1,0}(y,t)W_\beta (h_1^{}(z),\theta _1^{}(y))W_\beta (h_0^{}(s),\theta _0^{}(t))}`$ (114)
Finally, the susceptibilities at time step $`1`$ are given by

$`\chi _h(1)`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle DzDy\left(\frac{V_\beta (h(t),\theta (t))}{h(t)}\right)}|_{(h_1^{}(z),\theta _1^{}(y))}`$ (115)
$`\chi _\theta (1)`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{a}}}{\displaystyle DzDy\left(\frac{W_\beta (h(t),\theta (t))}{\theta (t)}\right)}|_{(h_1^{}(z),\theta _1^{}(y))}`$ (116)
### A.3 Third time step

We have all the quantities needed to write down the local fields at time $`t=2`$

$`h_i(2)`$ $`=`$ $`{\displaystyle \frac{\xi _i}{a}}m(2)+{\displaystyle \frac{\alpha }{a}}\chi _h(1)(\sigma _i(1)+\chi _h(0)\sigma (0))+𝒩(0,\alpha aD(2))`$ (117)
$`\theta _i(2)`$ $`=`$ $`\eta _il(2)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\chi _\theta (1)(\sigma _i^2(1)+\chi _\theta (0)\sigma ^2(0))+𝒩(0,{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}E(2))`$ (118)
This gives for the order parameters

$`m(3)`$ $`={\displaystyle \frac{1}{a}}\xi {\displaystyle D\omega _h^{2,0}(z,s)D\omega _\theta ^{2,0}(y,t)\underset{\sigma (1)}{}V_\beta (h_2^{}(z),\theta _2^{}(y))U_\beta (h_0^{}(s),\theta _0^{}(t),\sigma (1))}`$ (119)
$`q(3)`$ $`={\displaystyle D\omega _h^{2,0}(z,s)D\omega _\theta ^{2,0}(y,t)\underset{\sigma (1)}{}W_\beta (h_2^{}(z),\theta _2^{}(y))U_\beta (h_0^{}(s),\theta _0^{}(t),\sigma (1))}`$ (120)
$`l(3)`$ $`=\eta {\displaystyle D\omega _h^{2,0}(z,s)D\omega _\theta ^{2,0}(y,t)\underset{\sigma (1)}{}W_\beta (h_2^{}(z),\theta _2^{}(y))U_\beta (h_0^{}(s),\theta _0^{}(t),\sigma (1))}`$ (121)
with

$$U_\beta (h_0^{}(s),\theta _0^{}(t),\sigma (1))=\frac{\mathrm{exp}[\beta \sigma (1)h_0^{}(s)+\beta \sigma ^2(1)\theta _0^{}(t)]}{2\mathrm{exp}(\beta \theta _0^{}(t))\mathrm{cosh}(\beta h_0^{}(s))+1}$$
(122)
and

$`h_2^{}(z)`$ $`=`$ $`{\displaystyle \frac{\xi }{a}}m(2)+{\displaystyle \frac{\alpha }{a}}\chi _h(1)(\sigma (1)+\chi _h(0)\sigma (0))+\sqrt{\alpha aD(2)}z`$ (123)
$`\theta _2^{}(y)`$ $`=`$ $`\eta l(2)+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\chi _\theta (1)(\sigma ^2(1)+\chi _\theta (0)\sigma ^2(0))+\sqrt{{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}E(2)}y`$ (124)
In an analogous way further time steps could be calculated.

## Appendix B The complete GFA calculation

Starting from the path (58) with the transition probabilities given by (9) the generating function (60) can be written as

$`Z[\mathrm{\Psi },\mathrm{\Phi }]`$ $`={\displaystyle \underset{𝝈(0)}{}}\mathrm{}{\displaystyle \underset{𝝈(t)}{}}\text{Prob}(𝝈(0)){\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{s=1}{\overset{t}{}}}\mathrm{exp}\left[i\sigma _i(s)\left(\psi _i(s)+i\beta h_i(𝝈(s1))\right)\right]`$
$`\times \mathrm{exp}\left[i\sigma _i^2(s)\left(\varphi _i(s)+i\beta \theta _i(𝝈(s1))\right)\right]`$
$`\times \mathrm{exp}\left[\mathrm{ln}\left(2\mathrm{exp}(\beta \theta _i(𝝈(s1)))\mathrm{cosh}(\beta h_i(𝝈(s1)))+1\right)\right]`$ (125)
where we have left out the subscript $`N`$. We then isolate the local fields by inserting the appropriate delta distributions

$`Z[\mathrm{\Psi },\mathrm{\Phi }]`$ $`={\displaystyle \underset{𝝈(0)}{}}\mathrm{}{\displaystyle \underset{𝝈(t)}{}}\text{Prob}(𝝈(0)){\displaystyle \frac{d𝐡d\widehat{𝐡}d\theta d\widehat{\theta }}{(2\pi )^{2N(t+1)}}}`$
$`\times {\displaystyle \underset{s=1}{\overset{t}{}}}\{\mathrm{exp}[i{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{h}_i(s)(h_i(s)\gamma _{h,i}(s))+i{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{\theta }_i(s)(\theta _i(s)\gamma _{\theta ,i}(s))]`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i(s)(\psi _i(s)+i\beta h_i(s1))i{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i^2(s)(\varphi _i(s)+i\beta \theta _i(s1))\right]`$
$`\times \mathrm{exp}[\mathrm{ln}\{2\mathrm{exp}\left(\beta \theta _i(s1)\right)\mathrm{cosh}(\beta h_i(s1))+1\}]\text{}\}`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{h}_i(0)h_i(0)+i{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{\theta }_i(0)\theta _i(0)\right]`$
$`\times \mathrm{exp}[]`$ (126)
where $`\mathrm{exp}[]`$ includes all $`J`$ and $`K`$ dependent terms

$`\mathrm{exp}[]`$ $``$ $`\mathrm{exp}\left[iN{\displaystyle \underset{s=0}{\overset{t}{}}}\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{h}_i(s)\right)\left({\displaystyle \underset{j=1}{\overset{N}{}}}J_{ij}\sigma _j(s)\right)\right]`$ (127)
$`\times `$ $`\mathrm{exp}\left[iN{\displaystyle \underset{s=0}{\overset{t}{}}}\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{\theta }_i(s)\right)\left({\displaystyle \underset{j=1}{\overset{N}{}}}K_{ij}\sigma _j^2(s)\right)\right]`$
$`\times `$ $`\mathrm{exp}\left[i{\displaystyle \underset{s=0}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\widehat{h}_i(s)J_{ii}\sigma _i(s)+\widehat{\theta }_i(s)K_{ii}\sigma _i^2(s)\right)\right]`$
We remark that in the expression above, the diagonal terms are excluded by the introduction of the last factor, which we denote by $`M`$.

We first concentrate on the terms containing the disorder, i.e., the terms containing $`J_{ij}`$ and $`K_{ij}`$, of which we need to know the limit $`N\mathrm{}`$. Using the learning rule (6) and taking for simplicity one condensed pattern again, say $`\mu =1`$, we can study the condensed and the non-condensed part of the noise (127) by splitting

$$\mathrm{exp}[]=\mathrm{exp}[_c+_{nc}].$$
(128)
With respect to the condensed pattern it is seen that the following single time variables are relevant

$`m(s)`$ $``$ $`{\displaystyle \frac{1}{aN}}{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _i^1\sigma _i(s)k_h(s){\displaystyle \frac{1}{aN}}{\displaystyle \underset{i=1}{\overset{N}{}}}\xi _i^1\widehat{h}_i(s)`$ (129)
$`l(s)`$ $``$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\eta _i^1\sigma _i^2(s)k_\theta (s){\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\eta _i^1\widehat{\theta }_i(s)`$ (130)
To achieve site factorisation we then introduce appropriate delta distributions for (129)-(130) and their conjugated variables denoted with a hat. Consequently,

$`\mathrm{exp}[_c]={\displaystyle 𝑑𝐦𝑑\widehat{𝐦}𝑑𝐤_h𝑑\widehat{𝐤}_h𝑑𝐥𝑑\widehat{𝐥}𝑑𝐤_\theta 𝑑\widehat{𝐤}_\theta \frac{N^{4(t+1)}}{(2\pi )^{4(t+1)}}}`$
$`\times \mathrm{exp}\{iN(𝐤_𝐡𝐦+𝐤_\theta 𝐥))\}`$
$`\times \mathrm{exp}\left\{iN\left(\widehat{𝐦}𝐦+\widehat{𝐤}_𝐡𝐤_𝐡+\widehat{𝐥}𝐥+\widehat{𝐤}_\theta 𝐤_\theta \right)\right\}`$
$`\times \mathrm{exp}\left\{i{\displaystyle \underset{s=0}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\widehat{m}(s){\displaystyle \frac{\xi _i}{a}}^1\sigma _i(s)+\widehat{k}_h(s){\displaystyle \frac{\xi _i}{a}}^1\widehat{h}_i(s)+\widehat{l}(s)\eta _i^1\sigma _i^2(s)+\widehat{k}_\theta (s)\eta _i^1\widehat{\theta }_i(s)\right)\right\}`$ (131)
With respect to the non-condensed patterns, we first introduce gaussian integration $`D𝐱_{s=0}^t\left((2\pi )^{1/2}dx(s)\mathrm{exp}[\frac{1}{2}x^2(s)]\right)`$ to arrive at (we forget about superindex 1 for the condensed pattern)

$`\mathrm{exp}[_{nc}]`$ $`={\displaystyle }D𝐱_1D𝐱_2D𝐲_1D𝐲_2`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{1}{\sqrt{2iN}}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\xi _i}{a}}{\displaystyle \underset{s=0}{\overset{t}{}}}\left(x_1(\widehat{h}_i(s)+\sigma _i(s))+iy_1(s)(\widehat{h}_i(s)\sigma _i(s))\right)\right\}`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{1}{\sqrt{2iN}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\eta _i{\displaystyle \underset{s=0}{\overset{t}{}}}(x_2(\widehat{\theta }_i(s)+\sigma _i^2(s))+iy_2(s)(\widehat{\theta }_i(s)\sigma _i^2(s)))\right\}_{\xi _i}^{p1}`$
$`\times \mathrm{exp}\left\{i\alpha {\displaystyle \underset{s,s^{}=0}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{1}{a}}\widehat{h}_i(s^{})\sigma _i(s)+{\displaystyle \frac{1}{\stackrel{~}{a}}}\widehat{\theta }_i(s^{})\sigma _i^2(s)\right)\delta _{ss^{}}\right\}`$ (132)
with clear notation. We have also rewritten $`M`$ in a more suitable way. Performing the average over the patterns and expanding the exponentials up to order $`O(N^1)`$ it is seen that the following two time variables turn up

$`q(s,s^{})`$ $``$ $`{\displaystyle \frac{1}{Na}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i(s)\sigma _i(s^{})Q_h(s,s^{}){\displaystyle \frac{1}{Na}}{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{h}_i(s)\widehat{h}_i(s^{})`$ (133)
$`p(s,s^{})`$ $``$ $`{\displaystyle \frac{1}{N\stackrel{~}{a}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i^2(s)\sigma _i^2(s^{})Q_\theta (s,s^{}){\displaystyle \frac{1}{N\stackrel{~}{a}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{\theta }_i(s)\widehat{\theta }_i(s^{})`$ (134)
$`K_h(s,s^{})`$ $``$ $`{\displaystyle \frac{1}{Na}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i(s)\widehat{h}_i(s^{})K_\theta (s,s^{}){\displaystyle \frac{1}{N\stackrel{~}{a}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma _i^2(s)\widehat{\theta }_i(s^{})`$ (135)
Introducing again appropriate delta distributions we obtain after some matrix algebra

$`\mathrm{exp}[_{nc}]=\left({\displaystyle \frac{N}{2\pi }}\right)^{12(t+1)}{\displaystyle 𝑑𝐪𝑑\widehat{𝐪}𝑑𝐐_h𝑑\widehat{𝐐}_h𝑑𝐊_h𝑑\widehat{𝐊}_h𝑑𝐩𝑑\widehat{𝐩}𝑑𝐐_\theta 𝑑\widehat{𝐐}_\theta 𝑑𝐊_\theta 𝑑\widehat{𝐊}_\theta }`$
$`\times \mathrm{exp}\left[iN\text{Tr}\left(\widehat{𝐪}𝐪^{}+\widehat{𝐐}_𝐡𝐐_{𝐡}^{}{}_{}{}^{}+\widehat{𝐊}_𝐡𝐊_{𝐡}^{}{}_{}{}^{}+\widehat{𝐩}𝐩^{}+\widehat{𝐐}_\theta 𝐐_{\theta }^{}{}_{}{}^{}+\widehat{𝐊}_\theta 𝐊_{\theta }^{}{}_{}{}^{}\right)\right]`$
$`\times \mathrm{exp}\left[i\alpha N\text{Tr}\left(𝐊_𝐡+𝐊_\theta \right)\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{i}{a}}{\displaystyle \underset{s,s^{}}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\widehat{q}(s,s^{})\sigma _i(s)\sigma _i(s^{})+\widehat{Q}_h(s,s^{})\widehat{h}_i(s)\widehat{h}_i(s^{})+\widehat{K}_h(s,s^{})\sigma _i(s)\widehat{h}_i(s^{})\right)\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{i}{\stackrel{~}{a}}}{\displaystyle \underset{s,s^{}}{\overset{t}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(\widehat{p}(s,s^{})\sigma _i^2(s)\sigma _i^2(s^{})+\widehat{Q}_\theta (s,s^{})\widehat{\theta }_i(s)\widehat{\theta }_i(s^{})+\widehat{K}_\theta (s,s^{})\sigma _i^2(s)\widehat{\theta }_i(s^{})\right)\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{p}{2}}\left[\mathrm{ln}\left(𝐐_h𝐪(i𝐈𝐊_h^{})(i𝐈𝐊_h)\right)+\mathrm{ln}\left(𝐐_\theta 𝐩(i𝐈𝐊_\theta ^{})(i𝐈𝐊_\theta )\right)\right]\right]`$
with $`𝐈`$ the unit matrix of dimension $`(t+1)`$.

This concludes the discussion of the averaging over the disorder with as result a generating functional of the form $`Z\mathrm{exp}N\mathrm{}`$ which can be evaluated in the limit $`N\mathrm{}`$ via the saddle-point method.

$$\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}\{\mathrm{}\}\mathrm{exp}\left[N\left(\alpha \mathrm{ln}\mathrm{\Omega }+\mathrm{\Xi }+\mathrm{\Lambda }\right)\right]$$
(137)
with

$$\mathrm{\Omega }=\left([𝐐_h𝐪(i𝐈𝐊_h^{})(i𝐈𝐊_h)][𝐐_\theta 𝐩(i𝐈𝐊_\theta ^{})(i𝐈𝐊_\theta )]\right)^{\frac{1}{2}}$$
(138)
$$\begin{array}{c}\mathrm{\Xi }=i\left(\widehat{𝐦}𝐦+\widehat{𝐤}_𝐡𝐤_𝐡+𝐤_𝐡𝐦+\widehat{𝐥}𝐥+\widehat{𝐤}_\theta 𝐤_\theta +𝐤_\theta 𝐥\right)\hfill \\ \hfill +i\text{Tr}\left(\widehat{𝐪}𝐪^{}+\widehat{𝐐}_𝐡𝐐_{𝐡}^{}{}_{}{}^{}+\widehat{𝐊}_𝐡𝐊_{𝐡}^{}{}_{}{}^{}+\widehat{𝐩}𝐩^{}+\widehat{𝐐}_\theta 𝐐_{\theta }^{}{}_{}{}^{}+\widehat{𝐊}_\theta 𝐊_{\theta }^{}{}_{}{}^{}+\alpha (𝐊_𝐡+𝐊_\theta )\right)\end{array}$$
(139)
$`\mathrm{\Lambda }={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}\{{\displaystyle }d𝐡d\widehat{𝐡}d𝜽d\widehat{𝜽}{\displaystyle \underset{𝝈(0)\mathrm{}𝝈(t)}{}}\text{Prob}(𝝈(0))`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{s=0}{\overset{t}{}}}\widehat{h}_i(s)(h_i(s)\gamma _{h,i}(s)\widehat{k}_h(s){\displaystyle \frac{\xi _i}{a}})\right]`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{s=0}{\overset{t}{}}}\widehat{\theta }_i(s)(\theta _i(s)\gamma _{\theta ,i}(s)\widehat{k}_\theta (s)\eta _i){\displaystyle \underset{s=1}{\overset{t}{}}}\mathrm{ln}\left(2\mathrm{exp}\left(\beta \theta (s1)\right)\mathrm{cosh}(\beta h(s1))+1\right)\right]`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{s=0}{\overset{t}{}}}\left(\sigma (s)(\psi _i(s)+\widehat{m}(s){\displaystyle \frac{\xi _i}{a}}+i\beta h(s1))+\sigma ^2(s)(\varphi _i(s)+\widehat{l}(s)\eta _i+i\beta \theta (s1))\right)\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{i}{a}}{\displaystyle \underset{s,s^{}=0}{\overset{t}{}}}\widehat{q}(s,s^{})\sigma (s)\sigma (s^{})+\widehat{K}_h(s,s^{})\sigma (s)\widehat{h}(s^{})+\widehat{Q}_h(s,s^{})\widehat{h}(s)\widehat{h}(s^{})\right]`$
$`\times \mathrm{exp}[{\displaystyle \frac{i}{\stackrel{~}{a}}}{\displaystyle \underset{s,s^{}=0}{\overset{t}{}}}\widehat{p}(s,s^{})\sigma ^2(s)\sigma ^2(s^{})+\widehat{K}_\theta (s,s^{})\sigma ^2(s)\widehat{\theta }(s^{})+\widehat{Q}_\theta (s,s^{})\widehat{\theta }(s)\widehat{\theta }(s^{})]\}`$ (140)
At this point, some remarks are in order. First, we have been able to factorize the generating functional with respect to the sites. Further, the external fields are only present in $`\mathrm{\Lambda }`$. Therefore it is useful to introduce the effective single spin measure

$$f_{,i}=\frac{𝑑𝐡𝑑\widehat{𝐡}𝑑𝜽𝑑\widehat{𝜽}_{\sigma (0)\mathrm{}\sigma (t)}\text{Prob}(\sigma (0))_i[𝝈,\xi _i,𝐡,\widehat{𝐡},𝜽,\widehat{𝜽}]f}{𝑑𝐡𝑑\widehat{𝐡}𝑑𝜽𝑑\widehat{𝜽}_{\sigma (0)\mathrm{}\sigma (t)}\text{Prob}(\sigma (0))_i[𝝈,\xi _i,𝐡,\widehat{𝐡},𝜽,\widehat{𝜽}]}$$
(141)
where the auxiliary fields $`\psi _i`$ and $`\varphi _i`$ have been taken to be zero and $`_i`$ is given by

$`_i[\mathrm{}]`$ $`=\mathrm{exp}\left[i{\displaystyle \underset{s=0}{\overset{t}{}}}\widehat{h}(s)(h(s)\gamma _h(s)m(s){\displaystyle \frac{\xi _i}{a}})+i{\displaystyle \underset{i=1}{\overset{N}{}}}\widehat{\theta }(s)(\theta (s)\gamma _\theta (s)l(s)\eta _i)\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \underset{s=1}{\overset{t}{}}}\mathrm{ln}\left(2\mathrm{exp}\left(\beta \theta (s1)\right)\mathrm{cosh}(\beta h(s1))+1\right)\right]`$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{s=1}{\overset{t}{}}}\sigma (s)\left(i\beta h(s1)+k_h(s){\displaystyle \frac{\xi _i}{a}}\right)+\sigma ^2(s)\left(i\beta \theta (s1)+k_\theta (s)\eta _i\right)\right]`$
$`\times \mathrm{exp}\left[i\sigma (0)k_h(0){\displaystyle \frac{\xi _i}{a}}i\sigma ^2(0)k_\theta (0)\eta _i\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{i}{a}}{\displaystyle \underset{s,s^{}}{\overset{t}{}}}\widehat{q}(s,s^{})\sigma (s)\sigma (s^{})+\widehat{Q}_h(s,s^{})\widehat{h}(s)\widehat{h}(s^{})+\widehat{K}_h(s,s^{})\sigma (s)\widehat{h}(s^{})\right]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{i}{\stackrel{~}{a}}}{\displaystyle \underset{s,s^{}}{\overset{t}{}}}\left(\widehat{p}(s,s^{})\sigma ^2(s)\sigma ^2(s^{})+\widehat{Q}_\theta (s,s^{})\widehat{\theta }(s)\widehat{\theta }(s^{})+\widehat{K}_\theta (s,s^{})\sigma ^2(s)\widehat{\theta }(s^{})\right)\right]`$ (142)
The parameters $`𝐦,𝐤_h,𝐥,𝐤_\theta ,\widehat{𝐪},\widehat{𝐐}_h,\widehat{𝐊}_h,\widehat{𝐩},\widehat{𝐐}_\theta ,\widehat{𝐊}_\theta `$ in $`_i[\mathrm{}]`$ are defined as the solutions of the saddle-point equation $`d(\mathrm{\Xi }+\mathrm{\Lambda })=0`$. Furthermore, we have dropped the subindices $`i`$ from the external fields since they are taken to be site-independent.

Before determining the saddle-point solutions we write down the following identities

$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(s)}=\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{𝑑𝐦\mathrm{}𝑑\widehat{𝐊}_\theta e^{N(\mathrm{\Xi }+\mathrm{\Lambda })+\alpha N\mathrm{ln}(\mathrm{\Omega })+O(\mathrm{ln}(N))}\left[N\frac{\mathrm{\Lambda }}{\psi _i(s)}\right]}{𝑑𝐦\mathrm{}𝑑\widehat{𝐊}_\theta e^{N(\mathrm{\Xi }+\mathrm{\Lambda })+\alpha N\mathrm{ln}(\mathrm{\Omega })+O(\mathrm{ln}(N))}}=i\sigma (s)_{,i}$$
(143)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(s)}=i\sigma ^2(s)_{,i}$$
(144)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(s)\psi _j(s^{})}=\delta _{ij}[\sigma (s)\sigma (s^{})_{,i}\sigma (s)_{,i}\sigma (s^{})_{,i}]\sigma (s)_{,i}\sigma (s^{})_{,j}$$
(145)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(s)\varphi _j(s^{})}=\delta _{ij}[\sigma ^2(s)\sigma ^2(s^{})_{,i}\sigma ^2(s)_{,i}\sigma ^2(s^{})_{,i}]\sigma ^2(s)_{,i}\sigma ^2(s^{})_{,j}$$
(146)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(s)\varphi _j(s^{})}=\delta _{ij}[\sigma (s)\sigma ^2(s^{})_{,i}\sigma (s)_{,i}\sigma ^2(s^{})_{,i}]\sigma (s)_{,i}\sigma ^2(s^{})_{,j}$$
(147)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\psi _i(s)\gamma _{h,j}(s^{})}=\delta _{ij}[\sigma (s)\widehat{h}(s^{})_{,i}\sigma (s)_{,i}\widehat{h}(s^{})_{,i}]\sigma (s)_{,i}\widehat{h}(s^{})_{,j}$$
(148)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{^2\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\varphi _i(s)\gamma _{\theta ,j}(s^{})}=\delta _{ij}[\sigma ^2(s)\widehat{\theta }(s^{})_{,i}\sigma ^2(s)_{,i}\widehat{\theta }(s^{})_{,i}]\sigma ^2(s)_{,i}\widehat{\theta }(s^{})_{,j}$$
(149)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\gamma _{h,i}(s)}=i\widehat{h}(s)_{,i}$$
(150)
$$\underset{\mathrm{\Psi },\mathrm{\Phi }0}{lim}\frac{\overline{Z[\mathrm{\Psi },\mathrm{\Phi }]}}{\gamma _{\theta ,i}(s)}=i\widehat{\theta }(s)_{,i}$$
(151)
where we recall that the overline denotes the disorder average and where we used the normalization $`\overline{Z[0,0]}=1`$, we find that the last two equations are equal to $`0`$, since the limit $`\mathrm{\Psi },\mathrm{\Phi }0`$ can be interchanged with the derivative, due to the fact that the local fields are completely unrelated with the external fields. This leads to $`\widehat{h}(s)_{,i}=0`$, $`\widehat{\theta }(s)_{,i}=0`$, $`\widehat{h}(s)\widehat{h}(s^{})_{,i}=0`$ and $`\widehat{\theta }(s)\widehat{\theta }(s^{})_{,i}=0`$.

Recalling (61) we then have obtained the following identities

$`\overline{\sigma _i(s)}`$ $`=`$ $`\sigma (s)_{,i}`$ (152)
$`\overline{\sigma _i^2(s)}`$ $`=`$ $`\sigma ^2(s)_{,i}`$ (153)
$`\overline{C_{h,ij}(s,s^{})}`$ $`=`$ $`\delta _{ij}\sigma (s)\sigma (s^{})_{,i}+(1\delta _{ij})\sigma (s)_{,i}\sigma (s^{})_{,j}`$ (154)
$`\overline{C_{\theta ,ij}(s,s^{})}`$ $`=`$ $`\delta _{ij}\sigma ^2(s)\sigma ^2(s^{})_{,i}+(1\delta _{ij})\sigma ^2(s)_{,i}\sigma ^2(s^{})_{,j}`$ (155)
$`\overline{G_{h,ij}(s,s^{})}`$ $`=`$ $`\delta _{ij}\sigma (s)\widehat{h}(s^{})_{,i}`$ (156)
$`\overline{G_{\theta ,ij}(s,s^{})}`$ $`=`$ $`\delta _{ij}\sigma ^2(s)\widehat{\theta }(s^{})_{,i}`$ (157)
from which we obtain by using the LLN

$`C_h(s,s^{})`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{C_{h,ii}(s,s^{})}=\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma \left(s\right)\sigma \left(s^{}\right)_{,i}=\sigma \left(s\right)\sigma \left(s^{}\right)_{}_\xi `$ (158)
$`C_\theta (s,s^{})`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{C_{\theta ,ii}(s,s^{})}=\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\sigma ^2\left(s\right)\sigma ^2\left(s^{}\right)_{,i}=\sigma ^2\left(s\right)\sigma ^2\left(s^{}\right)_{}_\xi `$ (159)
$`G_h(s,s^{})`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{G_{h,ii}(s,s^{})}=\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{}{\gamma _{h,i}\left(s^{}\right)}}\sigma \left(s\right)_{,i}=i\sigma \left(s\right)\widehat{h}\left(s^{}\right)_{}_\xi `$ (160)
$`G_\theta (s,s^{})`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{G_{\theta ,ii}(s,s^{})}=\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{}{\gamma _{\theta ,i}\left(s^{}\right)}}\sigma ^2\left(s\right)_{,i}=i\sigma ^2\left(s\right)\widehat{\theta }\left(s^{}\right)_{}_\xi `$ (161)
where $`_\xi `$ denotes the average over the condensed pattern. Finally, we write down the saddle-point equations. Variation with respect to $`𝐦,𝐤_h,𝐥`$ and $`𝐤_\theta `$ gives

$$\widehat{m}(s)=k_h(s)\widehat{k}_h(s)=m(s)\widehat{l}(s)=k_\theta (s)\widehat{k}_\theta (s)=l(s)$$
(162)
Next, we calculate the variation with respect to $`\widehat{𝐦},\widehat{𝐤}_h,\widehat{𝐥}`$ and $`\widehat{𝐤}_\theta `$, leading to

$`m(s)`$ $`=`$ $`\xi \sigma (s)_{}_\xi k_h(s)=\xi \widehat{h}(s)_{}_\xi =0`$ (163)
$`l(s)`$ $`=`$ $`\eta \sigma ^2(s)_{}_\xi k_\theta (s)=\eta \widehat{\theta }(s)_{}_\xi =0`$ (164)
Variation with respect to $`\widehat{𝐪}`$, $`\widehat{𝐐}_h`$ and $`\widehat{𝐊}_h`$ gives

$$q(s,s^{})=\frac{1}{a}C_h(s,s^{})Q_h(s,s^{})=0K_h(s,s^{})=\frac{i}{a}G_h(s,s^{})$$
(165)
and with respect to $`\widehat{𝐩}`$, $`\widehat{𝐐}_\theta `$ and $`\widehat{𝐊}_\theta `$

$$p(s,s^{})=\frac{1}{\stackrel{~}{a}}C_\theta (s,s^{})Q_\theta (s,s^{})=0K_\theta (s,s^{})=\frac{i}{\stackrel{~}{a}}G_\theta (s,s^{})$$
(166)
Finally, variation with respect to of $`𝐪`$, $`𝐐_h`$, $`𝐊_h`$, $`𝐩`$, $`𝐐_\theta `$ and $`𝐊_\theta `$ i leads to, after some algebra

$`\widehat{q}(s,s^{})`$ $`=`$ $`0\widehat{p}(s,s^{})=0`$ (167)
$`\widehat{Q}_h(s,s^{})`$ $`=`$ $`{\displaystyle \frac{i\alpha a}{2}}\left[\left(a𝐈𝐆_h\right)^1𝐂_h\left(a𝐈𝐆_h^{}\right)^1\right](s,s^{})`$ (168)
$`\widehat{K}_h(s,s^{})`$ $`=`$ $`\alpha a\left[\left(a𝐈𝐆_h\right)^1\right](s,s^{})\alpha 𝐈(s,s^{})`$ (169)
$`\widehat{Q}_\theta (s,s^{})`$ $`=`$ $`{\displaystyle \frac{i\alpha \stackrel{~}{a}}{2}}\left[\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1𝐂_\theta \left(\stackrel{~}{a}𝐈𝐆_\theta ^{}\right)^1\right](s,s^{})`$ (170)
$`\widehat{K}_\theta (s,s^{})`$ $`=`$ $`\alpha \stackrel{~}{a}\left[\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1\right](s,s^{})\alpha 𝐈(s,s^{})`$ (171)
After the introduction of the saddle-point solutions, the single spin measure (141) is expressed in terms of the physical order parameters, $`\{𝐦`$ , $`𝐥`$, $`𝐂_h`$, $`𝐂_\theta `$, $`𝐆_h`$, $`𝐆_\theta \}`$. Rewriting the average in the following way

$$\begin{array}{c}f[\sigma (0),\mathrm{},\sigma (t)]dh(0)\mathrm{}dh(t1)d\theta (0)\mathrm{}d\theta (t1)\hfill \\ \hfill \underset{\sigma (0)\mathrm{}\sigma (t)}{}\text{Prob}[h(0),\mathrm{},h(t1),\theta (0),\mathrm{},\theta (t1);\sigma (0),\mathrm{},\sigma (t)]f[\sigma (0),\mathrm{},\sigma (t)]_\xi \end{array}$$
(172)
with

$`\text{Prob}\left[\mathrm{}\right]\text{Prob}(\sigma (0)){\displaystyle \underset{s=0}{\overset{t1}{}}dz_h(s)dz_\theta (s)}`$
$`\times \delta \left[h(s)\gamma _h(s)m(s)\xi \alpha {\displaystyle \underset{s^{}=0}{\overset{t}{}}}\left[\left(a𝐈𝐆_h\right)^1\right](s,s^{})\sigma (s^{})+{\displaystyle \frac{\alpha }{a}}\sigma (s)\sqrt{\alpha }z_h(s)\right]`$
$`\times \delta \left[\theta (s)\gamma _\theta (s)l(s)\eta \alpha {\displaystyle \underset{s^{}=0}{\overset{t}{}}}\left[\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1\right](s,s^{})\sigma ^2(s^{})+{\displaystyle \frac{\alpha }{\stackrel{~}{a}}}\sigma ^2(s)\sqrt{\alpha }z_\theta (s)\right]`$
$`\times {\displaystyle \frac{\mathrm{exp}\left(\beta \sigma (s+1)h(s)+\beta \sigma ^2(s+1)\theta (s)\right)}{2\mathrm{exp}\left(\beta \theta (s)\right)\mathrm{cosh}(\beta h(s))+1}}`$
$`\times {\displaystyle }d\widehat{𝐡}d\widehat{𝜽}\mathrm{exp}[{\displaystyle \frac{\alpha }{2}}{\displaystyle \underset{s,s^{}=0}{\overset{t}{}}}\left[(a𝐈𝐆_h)^1𝐂_h(a𝐈𝐆_h^{})^1\right](s,s^{})\widehat{h}(s)\widehat{h}(s^{})]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{\alpha }{2}}{\displaystyle \underset{s,s^{}=0}{\overset{t}{}}}\left[\left(\stackrel{~}{a}𝐈𝐆_\theta \right)^1𝐂_\theta \left(\stackrel{~}{a}𝐈𝐆_\theta ^{}\right)^1\right](s,s^{})\widehat{\theta }(s)\widehat{\theta }(s^{})\right]`$
$`\times \mathrm{exp}\left[i\sqrt{\alpha }{\displaystyle \underset{s=0}{\overset{t}{}}}\left(\widehat{h}(s)z_h(s)+\widehat{\theta }(s)z_\theta (s)\right)\right]`$ (173)
the integration over $`\widehat{𝐡}`$ and $`\widehat{𝜽}`$ is easy since the exponents multiplying the terms $`\widehat{h}(s)\widehat{h}(s^{})`$ and $`\widehat{\theta }(s)\widehat{\theta }(s^{})`$ are symmetric matrices. Doing the integration and including the normalization factors, we find the results of Section 4.

## Appendix C First two time steps in the GFA

We use the GFA results to calculate the first two time steps of the dynamics explicitly. The initial conditions are chosen as follows: $`\gamma (0)=\beta (0)=0`$, i.e., no external field, $`m(0)=m_0`$, $`l(0)=l_0`$, $`C_h(0,0)=C_\theta (0,0)=q_0`$, $`G_h(0,0)=G_\theta (0,0)=G_h(0,1)=G_\theta (0,1)=0`$.

### C.1 First time step

For time 1, the relevant probabilities are

$`w[z_h(0),z_\theta (0)]`$ $`={\displaystyle \frac{a\stackrel{~}{a}}{2\pi q_0}}\mathrm{exp}\left[{\displaystyle \frac{a^2}{2q_0}}(z_h(0))^2{\displaystyle \frac{\stackrel{~}{a}^2}{2q_0}}(z_\theta (0))^2\right]`$ (174)
$`\text{Prob}[\sigma (0),\sigma (1)|z_h(0),z_\theta (0)]`$ $`=\text{Prob}(\sigma (0)){\displaystyle \frac{\mathrm{exp}\left(\beta \sigma (1)h(0)+\beta \sigma ^2(1)\theta (0)\right)}{2\mathrm{exp}\left(\beta \theta (0)\right)\mathrm{cosh}(\beta h(0))+1}}`$ (175)
with $`h(0)=\xi m_0/a+\sqrt{\alpha }z_h(0)`$ and $`\theta (0)=l_0\eta +\sqrt{\alpha }z_\theta (0)`$. Then, using (74) and (75),

$`m(1)={\displaystyle \frac{1}{a}}\xi \sigma (1)`$ $`={\displaystyle }{\displaystyle }Dz_h(0)Dz_\theta (0)\{{\displaystyle \frac{\mathrm{exp}\left(\beta \theta ^+(0)\right)\mathrm{sinh}(\beta h^+(0))}{2\mathrm{exp}\left(\beta \theta ^+(0)\right)\mathrm{cosh}(\beta h^+(0))+1}}`$
$`+{\displaystyle \frac{\mathrm{exp}\left(\beta \theta ^+(0)\right)\mathrm{sinh}(\beta h^{}(0))}{2\mathrm{exp}\left(\beta \theta ^+(0)\right)\mathrm{cosh}(\beta h^{}(0))+1}}\}`$ (176)
with, as usual, $`Dz`$ a gaussian measure defined in (174) and where $`h^\pm (0)=\pm m_0/a+\sqrt{\alpha }z_h(0)`$ and $`\theta ^+(0)=l_0/a+\sqrt{\alpha }z_\theta (0)`$. In the same way we find the values for $`l(1)`$ and $`q(1)`$, i.e. $`l(1)=\eta \sigma ^2(1)`$ and $`q(1)=\sigma ^2(1)`$.

The correlation and response functions are

$`C_h(1,0)`$ $`=`$ $`C_h(0,1)=\sigma (0)\sigma (1)`$ (177)
$`C_\theta (1,0)`$ $`=`$ $`C_\theta (0,1)=\sigma (0)\sigma (1)`$ (178)
$`G_h(1,0)`$ $`=`$ $`\beta \sigma (1)\left[\sigma (1)V_\beta (h(0),\theta (0))\right]`$ (179)
$`G_\theta (t,t^{})`$ $`=`$ $`\beta \sigma ^2(1)\left[\sigma ^2(1)W_\beta (h(0),\theta (0))\right]`$ (180)
### C.2 Second time step

In this time step correlations start to appear.

$$𝐂_{h,\theta }^{t=2}=\left(\begin{array}{cc}q_0& C_{h,\theta }(1,0)\\ C_{h,\theta }(1,0)& q(1)\end{array}\right)$$

$$\left(a𝐈𝐆_h^{t=2}\right)^1=\frac{1}{a^2}\left(\begin{array}{cc}a& 0\\ G_h(1,0)& a\end{array}\right)\left(\stackrel{~}{a}𝐈𝐆_\theta ^{t=2}\right)^1=\frac{1}{\stackrel{~}{a}^2}\left(\begin{array}{cc}\stackrel{~}{a}& 0\\ G_h(1,0)& \stackrel{~}{a}\end{array}\right)$$

The explicit single-site average for time $`2`$ is determined by

$`w[z_h(0),z_h(1),z_\theta (0)z_\theta (1)]`$ $`={\displaystyle \frac{1}{4\pi ^2\sqrt{(q_0A_hB_h^2)(q_0A_\theta B_\theta ^2)}}}`$
$`\times \mathrm{exp}[{\displaystyle \frac{1}{2C_h}}(z_h(0))^2A_h+z_h(1))^2q_02z_h(0)z_h(1)B_h]`$
$`\times \mathrm{exp}[{\displaystyle \frac{1}{2C_\theta }}(z_\theta (0))^2A_\theta +z_\theta (1))^2q_02z_\theta (0)z_\theta (1)B_\theta ]`$ (181)
$`\text{Prob}[\sigma (0),\sigma (1)|z_h(0),z_\theta (0)]`$ $`=\text{Prob}(\sigma (0)){\displaystyle \frac{\mathrm{exp}\left(\beta \sigma (1)h(0)+\beta \sigma ^2(1)\theta (0)\right)}{2\mathrm{exp}\left(\beta \theta (0)\right)\mathrm{cosh}(\beta h(0))+1}}`$
$`\times {\displaystyle \frac{\mathrm{exp}\left(\beta \sigma (2)h(1)+\beta \sigma ^2(2)\theta (1)\right)}{2\mathrm{exp}\left(\beta \theta (1)\right)\mathrm{cosh}(\beta h(1))+1}}`$ (182)
with the local fields at time $`t=1`$ given by

$`h(1)`$ $`=`$ $`{\displaystyle \frac{\xi }{a}}m(1)+{\displaystyle \frac{\alpha }{a^2}}G_h(1,0)\sigma (0)+\sqrt{\alpha }z_h(1)`$ (183)
$`\theta (1)`$ $`=`$ $`l(1)\eta +{\displaystyle \frac{\alpha }{\stackrel{~}{a}^2}}G_h(1,0)\sigma ^2(0)+\sqrt{\alpha }z_\theta (1)`$ (184)
and the coeficients defined by

$`A_h`$ $``$ $`a^2q(1)+2aC_h(1,0)G_h(1,0)+q_0(G_h(1,0))^2`$ (185)
$`B_h`$ $``$ $`a^2C_h(1,0)+aq_0G_h(1,0)`$ (186)
$`C_h`$ $``$ $`q_0q(1)C_h^2(1,0)`$ (187)
$`A_\theta `$ $``$ $`\stackrel{~}{a}^2q(1)+2\stackrel{~}{a}C_\theta (1,0)G_\theta (1,0)+q_0(G_\theta (1,0))^2`$ (188)
$`B_h`$ $``$ $`\stackrel{~}{a}^2C_\theta (1,0)+\stackrel{~}{a}q_0G_\theta (1,0)`$ (189)
$`C_h`$ $``$ $`q_0q(1)C_\theta ^2(1,0)`$ (190)
From these relations any quantity can be computed. We present the final expressions for $`m(2)`$, $`l(2)`$ and $`q(2)`$.

$`m(2)`$ $`={\displaystyle \frac{1}{a}}\xi \sigma (2)`$ (191)
$`={\displaystyle \frac{\stackrel{~}{a}}{2\pi \sqrt{A_hA_\theta }}}{\displaystyle 𝑑z_h(1)𝑑z_\theta (1)\mathrm{exp}\left[\frac{1}{2}\left(\frac{(z_h(1))^2}{A_h}+\frac{(z_\theta (1))^2}{A_\theta }\right)\right]}`$
$`\times \{AV_\beta (h^{++}(1),\theta ^{++}(1))CV_\beta (h^+(1),\theta ^{++}(1))+DV_\beta (h^{+0}(1),\theta ^{+0}(1))`$
$`DV_\beta (h^0(1),\theta ^{+0}(1))+CV_\beta (h^+(1),\theta ^{++}(1))AV_\beta (h^{}(1),\theta ^{++}(1))\}`$ (192)
where the coeficients $`A`$, $`C`$ and $`D`$ coming from the distribution

$$\begin{array}{c}\text{Prob}(\sigma (0),\xi )=A\delta (\sigma (0)1)\delta (\xi 1)+B\delta (\sigma (0)1)\delta (\xi )\hfill \\ \hfill +C\delta (\sigma (0)1)\delta (\xi +1)+D\delta (\sigma (0))\delta (\xi 1)+E\delta (\sigma (0))\delta (\xi )\\ \hfill +D\delta (\sigma (0))\delta (\xi +1)+C\delta (\sigma (0)+1)\delta (\xi 1)\\ \hfill +B\delta (\sigma (0)+1)\delta (\xi )+A\delta (\sigma (0)+1)\delta (\xi +1)\end{array}$$
(193)
defined by

$`A`$ $`=`$ $`{\displaystyle \frac{a}{4}}(m_0+(1a)l_0+q_0)`$ (194)
$`B`$ $`=`$ $`{\displaystyle \frac{1a}{2}}(q_0al_0)`$ (195)
$`C`$ $`=`$ $`{\displaystyle \frac{a}{4}}(m_0+(1a)l_0+q_0)`$ (196)
$`D`$ $`=`$ $`{\displaystyle \frac{a}{2}}(1(1a)l_0+q_0)`$ (197)
$`E`$ $`=`$ $`(1a)(1+al_0q_0)`$ (198)
The local field contributions appearing in these final equations have the following form, due to the average over the initial conditions

$`h^{uv}(1)`$ $`=`$ $`u{\displaystyle \frac{m(1)}{a}}+v{\displaystyle \frac{\alpha }{a^2}}G_h(1,0)+\sqrt{\alpha }z_h(1)`$ (199)
$`\theta ^{xy}(1)`$ $`=`$ $`{\displaystyle \frac{xa}{\stackrel{~}{a}}}l(1)+y{\displaystyle \frac{\alpha }{\stackrel{~}{a}^2}}G_\theta (1,0)+\sqrt{\alpha }z_\theta (1)`$ (200)
where $`u,v\{1,0,1\}`$ and $`x,y\{0,1\}`$.

In the same way we find $`q(2)`$ and $`l(2)`$.