Mean-Field Theory in Second Quantization Formalism Consider the Ising model in statistical physics
$$H=-J\sum_{\left<i,j\right>}s_{i}s_{j}-\mu h\sum_{i}s_{i}$$
In this case mean-field approximation is done by replacing the surrounding spins with their averages $s_{i}\rightarrow\left<s_{i}\right>\equiv s$ in such a way to get an effective Hamiltonian of decoupled spins
$$H_{\rm MF}=-\left(Jzs+\mu h\right)\sum_{i}s_{i}$$
On the other hand, consider the Bose-Hubbard model
$$H-\mu N=-t\sum_{\left<i,j\right>}b_{i}^{\dagger}b_{j}+\frac{U}{2}\sum_{i}n_{i}\left(n_{i}-1\right)-\mu\sum_{i}n_{i}$$
People usually assume that in the mean-field regime the destruction operator acquires an expectation value $\psi\equiv \left<b_{i}\right>$. What justifies this assumption? Why is this mean-field - what is the analogy to the Ising model mean-field procedure? In addition, people continue to state that the operators decouple as
$$b^{\dagger}_{i}b_{j}=\left<b^{\dagger}_{i}\right>b_{j}+b^{\dagger}_{i}\left<b_{j}\right>-\left<b^{\dagger}_{i}\right>\left<b_{j}\right>$$
How is this equation derived? I find in many books claims such as $b_{i}=\left<b_{i}\right>+\delta b_{i}$ trying to explain the above equation, but without much justification. To me this seems like an assumption on the matrix elements of $b_{i}-\left<b_{i}\right>$ being small, but I fail to connect this to any mean-field argument.
 A: Let me respond to your points in a different order.

In addition, people continue to state that the operators decouple as
$$b^{\dagger}_{i}b_{j}=\left<b^{\dagger}_{i}\right>b_{j}+b^{\dagger}_{i}\left<b_{j}\right>-\left<b^{\dagger}_{i}\right>\left<b_{j}\right>$$
How is this equation derived?

Writing $b_i=\langle b_i\rangle + \delta b_i = \langle b_i \rangle + \left(b_i-\langle b_i\rangle\right),$ we have
$$b_i^\dagger b_j = \langle b_i^\dagger\rangle b_j + b_i^\dagger \langle b_j\rangle - \langle b_i^\dagger\rangle\langle b_j\rangle + \delta b_i^\dagger \delta b_j$$
Note that $\delta b_j$ is to be interpreted as a deviation away from the mean-field value. For mean-field theory to be accurate, this fluctuation should be relatively small, which leads us to neglect the $\delta b_i^\dagger \delta b_j \approx 0$ term. This really is the central assumption of mean-field theory -  that some parameter doesn't deviate too much away from its average value. 

Consider the Ising model in statistical physics
  $$H=-J\sum_{\langle i,j\rangle} s_is_j -\mu h\sum_is_i$$
  In this case mean-field approximation is done by replacing the surrounding spins with their averages si→⟨si⟩≡s in such a way to get an effective Hamiltonian of decoupled spins...

Actually, let's try the same mean-field procedure as before: $s_i = \langle s_i\rangle + \left( s_i - \langle s_i\rangle\right)$. Then you get
$$H=-J\sum_{\langle i,j\rangle} \left( \langle s_i\rangle s_j + s_i\langle s_j\rangle - \langle s_i\rangle\langle s_j\rangle + \delta s_i \delta s_j \right) -\mu h \sum_i \left( \langle s_i\rangle + \left( s_i - \langle s_i\rangle \right) \right)$$
Now we can first neglect the $\delta s_i \delta s_j$ term. Then we can recognize that $\sum \langle s_i\rangle \langle s_j\rangle$ only shifts the energy by a constant, and can neglect also this part. This yields your $H_{MF}$. [Note that neglecting the constant part is a standard simplification, but one that may fail sometimes, for example if we want to compare the energy of different phases in order to construct a phase diagram.]
Now, back to the Bose-Hubbard model:

People usually assume that in the mean-field regime the destruction operator acquires an expectation value $\psi \equiv \langle b_i\rangle$. What justifies this assumption? Why is this mean-field - what is the analogy to the Ising model mean-field procedure?

This is a rather subtle question. In a state with a fixed number of particles you'll have $\langle b_i\rangle \equiv 0$. That is, vanishing vacuum expectation values for operators changing the particle number. In a superfluid or Bose-Einstein-condensate (BEC) phase this changes, because the state can be (to good approximation) described by coherent states. This makes $\langle b\rangle$ a good order parameter - it's non-zero in the BEC phase, and zero in the normal phase. Similarly, the magnetization $\langle s\rangle$ is a good order parameter for the Ising model. The mean-field procedure then proceeds (in both cases) by a mean-field approximation around the presumed finite value of the order parameter.
In the case of a magnet you can literally interpret the mean-field as a background magnetic field coming from the other spins. More generally, the mean-field is the average effect from the surrounding sites/particles. In the BEC case, the effect is due to the surrounding condensate.
