Hi Yi, the questions and issues you raise here are good ones. First of all, regarding your three ideas, I would say there is no single best one, but it depends on various details like the kind of order you are studying, the sizes of systems you can reach, and the accuracy of your calculations. Let us get back to that in a moment.
One point you raised is about SU(2) conservation of DMRG and the idea that not having it could prevent seeing symmetry breaking. In fact, the situation is the exact opposite: conserving a symmetry means the code is not allowed to break that symmetry, even if the actual physical system would. For example, if one conserves Z2 spin parity in the transverse-field Ising model (something you can do in ITensor since it is an Abelian symmetry), then the ground states DMRG will return in the symmetry breaking ferromagnetic phase will actually be linear combinations of all spins up plus all spins down, thus not actually breaking the symmetry. Whereas if you studied the system with all symmetry conservation turned off, then the code could return a ground state which is just all spins up, say (or all spins down and which one depends on subtle details of your initial state).
This ties into your other point about finite systems not being able to break a symmetry. Technically this is absolutely correct, however in practice the energy difference between a state which breaks the symmetry and one which does not is extremely tiny, and goes to zero exponentially with system size. So in DMRG, unless you study rather small systems and push the calculation to very high accuracy, often the approximate ground state DMRG returns will in fact break the symmetry if that is what the system wants to do. But also in practice DMRG might not break the symmetry in a "clean" way, and could return a state which is some mixture of two or more of the symmetry-equivalent ground states. Which brings me to my next point.
The safest "best practice" for observing symmetry breaking with DMRG is usually your third idea, about applying a field. Except here the best idea is usually to apply the field only at the edges of the system, so that you minimally disrupt the physics in the bulk of the system. See for example this review article which discusses this idea in more detail: https://arxiv.org/abs/1105.1374
With the edge-pinning idea, you don't have to adjust the size of the field, but can just adjust the length of your system and do finite-size scaling, which is nice.
Measuring correlation functions can also be a good idea, but because you are measuring the square of two operators, this quantity can be affected more by approximations such as the truncation of an MPS that is done by DMRG. The number that you obtain from this method is smaller and so can be harder to converge. But in many cases it can still be a good method to use, because it can be an easier method for detecting complicated orders which might be more easily missed by just guessing an order parameter to measure.
So briefly, I would usually recommend a combination of (1) and (3), with the fields in (3) only applied at the edges. But (2) can also be a good option. Ideally trying both together and comparing the results is the best, especially when first exploring a new system and trying to understand it.
Hope that's helpful -
Miles