directly with the covariance matrix is given by the relation $$\label{eq:eq2}
\frac{\partial G}{\partial J_{ij}} = - m_{ij}.$$ The Plefka expansion is used to expand the Gibbs free energy in power of the coupling $J_{ij}$ assumed to be small. Multiplying all coupling $J_{ij}$ by $\alpha$ yields the following cluster expansion: $$\begin{aligned}
G[{{\bf m}},\alpha{{\bf J}}] &= {{\bf h}}^{T}({{\bf m}},\alpha){{\bf m}}+ F[{{\bf h}}({{\bf m}},\alpha),\alpha{{\bf J}}]\label{eq:plefka1}\\[0.2cm]
&= G_0[{{\bf m}}] + \sum_{n=0}^\infty\frac{\alpha^n}{n!}G_n[{{\bf m}},{{\bf J}}]\label{eq:plefka2}\end{aligned}$$ where each term $G_n$ corresponds to cluster contributions of size $n$ in the number of links $J_{ij}$ involved, and ${{\bf h}}({{\bf m}},\alpha)$ depends implicitly on $\alpha$ in order to always fulfill (\[eq:constraints\]). This precisely is the Plefka expansion, and each term of the expansion (\[eq:plefka2\]) can be obtained by successive derivation of (\[eq:plefka1\]). We have $$G_0[{{\bf m}}] = \sum_i\frac{1+m_i}{2}\log\frac{1+m_i}{2}+\frac{1-m_i}{2}\log\frac{1-m_i}{2}.$$ Letting $$H_J {\stackrel{\mbox{\upshape\tiny def}}{=}}\sum_{i<j} J_{ij}s_is_j,$$ using (\[eq:constraints\]), the two first derivatives of (\[eq:plefka1\]) w.r.t $\alpha$ read $$\begin{aligned}
\label{eq:dG1}
\frac{dG[{{\bf m}},\alpha{{\bf J}}]}{d\alpha} &= -{\mathbb E}_\alpha\bigl(H_J\bigr), \\
\label{eq:dG2}
\frac{d^2G[{{\bf m}},\alpha{{\bf J}}]}{d\alpha^2} &= -\operatorname{Var}_\alpha\bigl(H_J\bigr)
-\sum_i\frac{dh_i({{\bf m}},\alpha)}{d\alpha}\operatorname{Cov}_\alpha\bigl(H_J,s_i\bigr),\end{aligned}$$ where subscript $\alpha$ indicates that expectations, variance and covariance are taken at given $\alpha$. To get successive derivatives of ${{\bf h}}({{\bf m}},\alpha)$ one can use (\[eq:dual1\]). Another possibility is to express the fact that ${{\bf m}}$ is fixed, $$\begin{aligned}
