How can I simulate a ground state degenerate system numerically? I'm using numerical method like DMRG to simulate ground state of correlated systems. But the degeneracy of the ground state has long bothered me:
When degeneracy exists the ground state isn't unique. So the ground state found by DMRG can be different from that in experiments: they could be in different state but both in ground state! So I cannot faithfully compare this system with my results in DMRG. 
I understand that general diagonalization algorithms find a basis for vector space. If you want the most general state you should do the superposition. But I don't know how to do superposition with existed DMRG code and I think it's not so much practical to search the whole vector space numerically. So since different ground state would be found, how can I possibly simulate any ground state degenerate system numerically by DMRG? 
 A: I'll go ahead and say that one can roughly classify systems into


*

*Systems with a unique ground state.

*Systems that have multiple (possibly infinitely many) degenerate ground states, but that will tend to select a unique one by some mechanism.

*Systems with topological ground state degeneracy.


(I won't rule out there existing  additional classes, and I'd welcome learning about them in the comments, but this classification seems satisfactory for the current answer.)

Clearly, class 1 is easy to deal with. There might be some system where it's difficult to numerically resolve the gap to the first excited state, but in practice, with finite size effects and all, the gap will usually be larger than machine resolution. An example of such a system would be spins coupled to a polarizing magnetic field.
Class 2 represents includes some of our most well studied systems, such as the antiferromagnetic Ising model $H_I=\sum_{\langle i,j\rangle} S_i^zS_j^z$ on the square lattice, or the Heisenberg model $H_H=\sum_{\langle i,j\rangle} \mathbf{S}_i \cdot \mathbf{S}_j$ on a cubic lattice. Clearly both models allow for a set of degenerate ground states. You're correct that for small systems, and in principle also larger ones, the system could be in a superposition of different several such states. However, in condensed matter we're usually interested in phases of matter, a concept which is only really well-defined in the thermodynamic limit - as the system size tends to infinity.
In that limit we tend to have spontaneous symmetry breaking - that is, for whatever reason (e.g. a small perturbation from the environment) the system picks out a ground state of lower symmetry than the Hamiltonian. If this wasn't the case, a permanent magnet would be a superposition of a state with all spins along $+\hat{z}$ and a state with all spins along $-\hat{z}$, and averaging the magnetization would yield zero - clearly an unphysical result. In other words, the physically observable quantities in macroscopic systems tend to arise from a symmetry broken ground state, not a superposition of degenerate states. Finally, for completeness, I will also mention that in some systems there is also a phenomenon called order by disorder that helps reduce the ground state degeneracy by entropic mechanisms.
For a DMRG calculation (and some other methods), assuming you want to say something about large systems, you will often want to include a small symmetry breaking term to guide the system to the ground state you're interested in. Often one can pick a tiny magnetic field (uniform for a ferromagnetic state, staggered for a Néel state). This picture looks  very clean for the kind of magnetic states I've discussed above, one's related by a symmetry of the Hamiltonian, since the symmetry allows us to pick our coordinate system as we like. (That's essentially what we're doing by adding a magnetic field - picking out the natural spin quantization axis.) Now, the picture would be less clear in a system with accidental (or maybe pathological) degeneracies between qualitatively very different states. In that case picking the wrong state likely leads to predictions that disagree with the experiment, so one could just try again.
For smaller systems, where superposition might be a reality, I'd probably try to stick with exact diagonalization techniques myself. People studying systems from class 3 tend to utilize very sophisticated approaches, see e.g. this paper. But based on your comment I think class 2 interests you more.
A: In complement to Anyon's answer, let me add two things:


*

*You can re-run your DMRG code with different initial conditions; you would assume that this gives different ground states in every run.  If you have a system with symmetry breaking, you will typically get the symmetry broken states in DMRG, but even if not, you can compute their overlaps and thus determine the dimension and basis of the ground space.  (The only case which might be tricky is if the entanglement in the different ground states is vastly different -- in that case, the DMRG method might tend to always find the less entangled state.

*You can first find one ground state, and then run DMRG again with the additional constraint that the new state is orthogonal to the state you already found, which is just another quadratic constraint and can be easily incorporated. This is e.g. done to find excited states, but obviously it will also return a basis for the space of ground states.
