Difference between a Hamiltonian and its mean-field form? In most solid state physics text books, BCS theory is introduced in the following way:
The authors introduce a mean-field parameter (based on the superconduncting pair) and upon that, they reduce the four-operator interaction term into a two-operator one. Thus a mean-field Hamiltonian is obtained. Solving the (quadratic) mean-field Hamiltonian, we get the BCS results.
I don't understand the relationship between the mean-field Hamiltonian and its original form. If Hamiltonian has a symmetry that is absent in the ground-state, a spontaneous symmetry-breaking has occurred. But the BCS mean-field Hamiltonian does not obey particle number conservation anymore.
Is the ground-state of the mean-field Hamiltonian still similar to the ground-state of the original Hamiltonian?
I am totally confused whether the particle number is conserved or not.
 A: Mathematically speaking, the answer is yes. The reason is that the ground state must minimise the energy of the system. Mean field is just a technique to calculate (sometimes only approximately) the true ground state. If the mean-field approximation is correct, you expect to reach the true ground state of the system, or to be close to the true ground state.
In mean-field treatment, you first guess a form of an order parameter. This somehow means you guess which kind of ground state you expect. If the ground state of the mean-field Hamiltonian were not the one of the original Hamiltonian (or close to it), you would typically get diverging quantities somewhere during the calculation. 
This reflects the fact that fluctuations becomes too strong for the ground state you expected to be realised in the system under study. Indeed, when writing the mean-field approximation, you kill the fluctuations at first. This is just what you do when writing
$$ccc^{\dagger}c^{\dagger}\rightarrow\left\langle cc\right\rangle c^{\dagger}c^{\dagger}+cc\left\langle c^{\dagger}c^{\dagger}\right\rangle $$
the term $ccc^{\dagger}c^{\dagger}$ have more informations coded in it than the averaged form on the right-hand-side. 
In short, if you get divergencies while calculating things, it's because the fluctuations are too strong. This in turn means that the system would not remains in the ground state you expected at first when choosing the mean-field approximation. Thus it will transit to an other state, and this state should be of lower energy, hence would be a better approximation as a ground state.
What you can do then is to try an other ground state, for instance in our simplistic example one could have chosen an other configuration for the mean-field approximation
$$ccc^{\dagger}c^{\dagger}\rightarrow\left\langle cc^{\dagger}\right\rangle c^{\dagger}c+cc^{\dagger}\left\langle cc^{\dagger}\right\rangle $$
instead. If all mean field treatments with all possibilities of order-parameter are exhausted without any success, then it means for sure that your system does not accept mean-field treatment as an accurate approximation.
