I would like to offer a contrasting answer to that of @tparker. I want to emphasize the fact that it is actually not necessary to introduce any symmetry-breaking field, provided you use the proper setup. (Of course, you can use symmetry-breaking fields, or suitable boundary conditions to achieve that, but I believe that it is conceptually interesting that you don't have to.)
I'll only discuss the case of classical systems, as I am much more familiar with the relevant mathematical framework.
The goal is to define a Gibbs measure for an infinite system. This is necessary, as only truly infinite systems undergo genuine phase transitions (of course, large finite systems can display approximate, smoothed-out "phase transitions"). The usual definition relying on the Boltzmann weight $e^{-\beta H}$ is useless in this case, as the energy of an infinite system is usually undefined.
The most efficient framework to solve this problem, the so-called Dobrushin-Lanford-Ruelle theory, works as follows. I discuss the case of the Ising model on $\mathbb{Z}^2$ for simplicity, but the approach is completely general. One says that a probability measure $\mu$ on the set of infinite configurations $\{-1,1\}^{\mathbb{Z}^2}$ is an infinite-volume Gibbs measure if, for any finite set $\Lambda\subset\mathbb{Z}^2$ and any configuration $\eta\in\{-1,1\}^{\mathbb{Z}^2\setminus\Lambda}$, the conditional probability of seeing a configuration $\sigma$ inside $\Lambda$, given that the configuration outside $\Lambda$ is given by $\eta$, is given by the finite-volume Gibbs measure in $\Lambda$ with boundary condition $\eta$. The latter is well-defined (through the usual Boltzmann weight), since $\Lambda$ is finite.
For models with compact spins (such as the Ising model or the Heisenberg model), existence of infinite-volume Gibbs measures is guaranteed. Uniqueness, however, does not hold in general. For the Ising model on $\mathbb{Z}^2$, for example, there exists a critical value $\beta_{\rm c}\in(0,\infty)$ of the inverse temperature such that there is a unique infinite-volume measure when $\beta\leq \beta_{\rm c}$, while there are infinitely-many infinite-volume Gibbs measures when $\beta>\beta_{\rm c}$. It turns out that any infinite-volume Gibbs measure $\mu$ can be expressed as a convex combination of two of them: $\mu = \alpha \mu^+_\beta + (1-\alpha)\mu^-_\beta$, for some $0\leq\alpha\leq 1$. These two measures $\mu^+_\beta$ and $\mu^-_\beta$ thus contain all the relevant physics, and there are good reasons to consider them as the physically truly relevant ones. It turns out that the average magnetization under $\mu^+_\beta$ is equal to $m^*(\beta)>0$, the spontaneous magnetization (that is the value of the magnetization that you would get, had you first added a magnetic field $h>0$ and then let $h$ decrease to $0$). Under $\mu^-_\beta$, on the other hand, the average magnetization is $-m^*(\beta)<0$. In this precise sense, there is spontaneous symmetry breaking, even though the procedure described above does not explicitly break the symmetry at any step.
Let me now briefly link the above approach (still for the Ising model on $\mathbb{Z}^2$) with the one using an explicit symmetry breaking. A standard way to proceed in this case is to consider an increasing sequence of finite subsets $\Lambda_n\subset\mathbb{Z}^2$. We then consider the Gibbs measure $\mu^+_{\Lambda_n;\beta}$ associated to the Ising model in $\Lambda_n$, with $+$ boundary condition. One can then consider the (weak) limit of the probability measures $\mu^+_{\Lambda_N;\beta}$ as $\Lambda_n$ grows to cover $\mathbb{Z}^2$. It turns out that the limiting measure coincides with the infinite-volume Gibbs measure $\mu_\beta^+$ derived above. Of course, one recovers the measure $\mu^-_\beta$ by using a sequence with $-$ boundary condition.
So the two approached yield the same result, but I insist that the former does not require explicit symmetry breaking. (Moreover, it provides a much more powerful framework.)
