\frac{dm_i}{d\alpha} = 0 & = -\frac{d}{d\alpha}\frac{\partial F[{{\bf h}}(\alpha),\alpha{{\bf J}}]}{\partial h_i}\nonumber\\[0.2cm]
& = \sum_{j}h'_j(\alpha)\operatorname{Cov}_\alpha(s_i,s_j)+\operatorname{Cov}_\alpha(H_J,s_i), \end{aligned}$$ giving $$\label{eq:dha}
h'_i(\alpha) = - \sum_j[\chi_\alpha^{-1}]_{ij}\operatorname{Cov}_\alpha(H_J,s_j),$$ where $\chi_\alpha$ is the susceptibility delivered by the model when $\alpha\ne 0$. To get the first two terms in the Plefka expansion, we need to compute these quantities at $\alpha=0$: $$\begin{aligned}
\operatorname{Var}\bigl(H_J\bigr) &= \sum_{i<k,j} J_{ij}J_{jk}m_im_k(1-m_j^2)+\sum_{i<j}J_{ij}^2(1-m_i^2m_j^2),\\[0.2cm]
\operatorname{Cov}\bigl(H_J,s_i\bigr) &= \sum_{j}J_{ij}m_j(1-m_i^2),\\[0.2cm]
h'_i(0) &= - \sum_{j} J_{ij}m_j,\end{aligned}$$ (by convention $J_{ii}=0$ in these sums). The first and second orders then finally read: $$G_1[{{\bf m}},{{\bf J}}] = -\sum_{i<j}J_{ij}m_im_j,\qquad\qquad G_2[{{\bf m}},{{\bf J}}] = -\sum_{i<j}J_{ij}^2(1-m_i^2)(1-m_j^2),$$ and correspond respectively to the mean field and to the TAP approximation. Higher order terms have been computed in [@GeYe].

At this point we are in position to find an approximate solution to the inverse Ising problem, either by inverting equation (\[eq:dual2\]) or (\[eq:eq2\]). To get a solution at a given order $n$ in the coupling, solving (\[eq:eq2\]) requires $G$ at order $n+1$, while it is needed at order $n$ in (\[eq:dual2\]).

Taking the expression of $G$ up to second order gives $$\frac{\partial G}{\partial J_{ij}} = -m_im_j-J_{ij}(1-m_i^2)(1-m_j^2),$$ and (\[eq:eq2\]) leads directly to the basic mean-field solution: $$\label{eq:JMF}
J_{ij}^{MF} = \frac{{{\hat \chi}}_{ij}}{(1-{\hat m}_i^2)(1-{\hat m}_j^2)}.$$ At this level of approximation for $G$, using (\[eq:dual1\]) we 