Mean field theory formulation of Ising model I am having a bit of trouble with basic combinatorics pertaining to the Ising model and mean field theory. Specifically, I get that the Hamiltonian can be written
$H=-\frac{J}{2}\Sigma_{i,j} s_is_{i+j} - \mu \Sigma_i s_i B$
Where the notation in the first sum denotes sum over all spins, i, and it nearest neigbours, indexed by j. The factor of 1/2 comes from double counting.
However then in my lectures notes, and other texts, when the mean spin is replaced by s, and the 2D nearest neighbours are explicitly put in, this somehow becomes
$H=-\Sigma_i s_i (\mu B + 2DJs)$
But this seems to be off by a factors of 2!
If s is the mean spin, then the sum over the nearest neighbours is
$\Sigma_{j} s_{i+j} = 2Ds$
i.e. number of nearest neighbours multiplied by the mean spin. The factor of 2 then cancels with the factor of $\frac{1}{2}$ in the first form of the Hamiltonian, so I get that the new Hamiltonian should be 
$H=-\Sigma_i s_i (\mu B + DJs)$ ?
Where has the factor of 2 gone?
 A: This can be a bit confusing. But actually the Wikipedia page explains the counting fairly clearly. I prefer their notation (summing over nearest-neighbour pairs, i.e. edges of the lattice) to the double-counting notation. So the interaction term is written
$$
-J\sum_{\langle i,j\rangle} s_i s_j \tag{1}
$$
summing over all distinct nearest-neighbour pairs $\langle i,j\rangle$, so no factor $1/2$.
The approximation is obtained by writing 
$$
s_i = s + (s_i-s)
$$
and similarly for $s_j$. Then every pair term becomes
$$
s_is_j = s^2 + s(s_i-s) + s(s_j-s) + (s_i-s)(s_j-s) .
$$
The approximation consists of dropping the last term.
Also, all the terms in $s^2$ are constant, and can be dropped.
(Actually, more precisely, I should say that they can be 
taken out of the problem, and put back in again if we wish,
once the value of $s$ has been determined).
So we are left with
$$
-J\sum_{\langle i,j\rangle} s\,s_i + s\,s_j = -2J\sum_{\langle i,j\rangle} s\,s_i .
$$
Notice that there are two terms here, not one, for every edge.
Compare with eqn (1). This is where the factor of $2$ comes in.
Physically we are counting the fluctuation of $i$ in the mean field provided by $j$, and the fluctuation of $j$ in the mean field provided by $i$.
The $i$ and $j$ terms are equivalent, of course,
hence the second form of the equation.
If the dimensionality is $D$, there are $2D$ edges per site $i$,
but we take care to avoid double counting, so
$$
\sum_{\langle i,j\rangle} \rightarrow \frac{1}{2} \, 2D \sum_i .
$$
So the result may also be written
$$
-2DJs\, \sum_{i} s_i
$$
as you have found in notes and text books.
