As far as I understand classical symmetry breaking occurs because of the inadequacy of the ensemble to account for the true thermodynamics of the system. Consider a simple classical Ising model: $$ H = -J \sum_{<i, j>} s_I s_j $$
the canonical partition function would count the configuration all spins up and all spin down as having the same probabilities. This is indeed true but when a system chooses to realize one of the configurations in the ensemble either all spins up or down, the thermal fluctuation in the ordered phase are not strong enough to flip all the spins and averaging over the canonical ensemble makes actually no sense. After applying the Hubbard-Stratonovich transformation one gets to Landau-Ginzburg free energy. The same problem arise here but one can perform a saddle point approximation (around the chosen minima) to account for the relevant fluctuation in the thermodynamical functions. The field-integral of the Landau-Ginzburg free energy still consider the two possibilities for the magnetization on equal footing, so if one wants to obtain physically sensible results (for the specific heat for example) one has to expand around one specific minima.
In the quantum Ising model we have the same Hamiltonian but now $s$ are quantum-mechanical spins. The ground state of the system is the superposition: $$ |\psi \rangle = \frac{|\uparrow \uparrow ... \uparrow \rangle + |\downarrow \downarrow ... \downarrow \rangle}{\sqrt{2}} $$
This state is not stable under little decoherence, thus realizing the spontaneous symmetry breaking. This discussions:
Spontaneous Quantym Symmetry Breaking
somehow support this point of view.
The other ways to express the presence of a phase with broken symmetry are indeed equivalent (I think): ODLRO or $\lim_{h \rightarrow 0} \lim_{N \rightarrow \infty} \langle s_I \rangle \neq 0$.
Of course if the ground state are really a degenerate manifold then one can choose a ground state which hasn't the symmetries of the Hamiltonian, but this isn't strictly necessary. Interactions in quantum systems typically break degeneracy producing a unique ground state which transforms as a singlet of the symmetry group $G$ of the Hamiltonian. This indeed happens in the above example.How can one talk about spontaneous symmetry breaking in those cases? Lesson from the Ising model is that one can identify a SSB ground state if it isn't stable under little decoherence.
Now I want to apply this all to the case of BEC, superfluids and superconductors. In those cases one typically works in the grand canonical ensemble. Let's start from BECs. One has the Hamiltonian
$$\mathcal{H} = \sum_{\vec{k}} \left( \epsilon_{\vec{k}} - \mu \right) a^{\dagger}_{\vec{k}} a_{\vec{k}}$$
For a fixed mean number of particles $N$ the coherent states $|\alpha \rangle = e^{-\frac{|\alpha|^2}{2}}\sum_{n} \frac{\alpha^n}{\sqrt{n!}} |n, 0, 0, ... \rangle$ have the same energy than the fixed number of particle state: $|n, 0, 0, ... \rangle$. By the way these are not eigenstate of the Hamiltonian but only satisfy the mean value condition and are only an ansatz, I probably fail to perceive the importance of them. $|\alpha|$ is fixed by the average number of particles. Now those coherent states are a $U(1)$ manifold of ''ground states''. Its notable how the ground state would have been unique in the canonical ensemble but in the grandcanonical it actually become degenerate. Now these are not orthogonal at all. I fail to see why a cooling non interacting gas of bosons should choose one of the coherent states as a ground state. Moreover as stated this degeneracy is absent in the canonical ensemble where the system will simply reach the state $|n, 0, 0, ... \rangle$. It seems to me that in the canonical ensemble there is no symmetry breaking and I don't completely understand how it is realized in the grandcanonical. It doesn't seem to me that the state $|n, 0, 0, ... \rangle$ is not stable under dechoerence.
In the context of interacting bose gas and superfluid I have the same problem in understanding what does the $U(1)$ symmetry breaking means. I understand the macroscopic description given in terms of Landau-Ginzburg functional, which have a manifold of degenerate ground state for the order parameter $\psi$ among which one and only one is realized but I fail to understand the microscopic origin of the so called $U(1)$ spontaneous symmetry breaking.
Looking at superconductivity we have the BCS Hamiltonian:
$$ \mathcal{H} = \sum_{\vec{k}} \left( \epsilon_{\vec{k}} - \mu \right) c^{\dagger}_{\vec{k}} c_{\vec{k}} + U_{0} \sum_{\vec{p}, \vec{q}} c^{\dagger}_{\vec{p}, \uparrow} c^{\dagger}_{-\vec{p}, \downarrow} c_{\vec{q}, \uparrow} c_{-\vec{q}, \downarrow} $$
One typically make the BCS ansatz finding the best value of $u_{\vec{k}}$ and $v_{\vec{k}}$ in $|BCS \rangle = \prod_{\vec{k}} \left( u_{\vec{k}} + v_{\vec{k}} c^{\dagger}_{\vec{k}, \uparrow} c^{\dagger}_{-\vec{k}, \downarrow} \right) |0 \rangle$ that minimize the energy as an atsatz for the ground state. $|BCS \rangle$ state is analog to coherent state for Cooper pairs. Once again working in the grandcanonical ensemble introduce degeneracy that probably wouldn't have been there in the canonical ensemble (is this true?). Than is the true ground state of the BCS Hamiltonian with fixed number of electrons not stable under decoherence? What is the microscopic mechanism that causes the spontaneous symmetry breaking analogous to the one found for the quantum Ising model in superfluid, BEC and superconductors?