As the title of the paper suggests, this is a spontaneous symmetry breaking. That is, it is not at the level of the Hamiltonian (which remains symmetric) but rather at the level of the state of the system. Just like in the 2d Ising model, or in many magnets - the Hamiltonian describing the system is symmetric, and therefore there are two degenerate ground states for example (in the 2d Ising model - all the spins pointing "up" or all of them pointing "down"). However, below some temperature, the system "chooses" one of the states over the other, and breaks the symmetry spontaneously. The fact that such phase-transitions (which are non-analytic behavior of the system) can happen is something that concentrated a lot of theoretical work in the previous century. It is dependent on taking the thermodynamic limit of system size to infinity, which introduces this non-analytic behavior.

One formal way to examine if a symmetry is broken spontaneously is to add a symmetry breaking term and then send it to zero while taking the size of the system to infinity, in the proper order. In your example, it will entail $$ H(\eta) = -\sum_{i,j} \left[\sigma^x_i \sigma^x_j + \sigma^y_i \sigma^y_j + \Delta \sigma^z_i \sigma^z_j\right] + \eta \sum_i \sigma^z_i$$ and then examine the magnetization per spin when the proper limits are taken $$m_z = \lim_{\eta\to 0^+}\lim_{N\to\infty}\frac{1}{N}\sum_i\langle \sigma^z_i\rangle$$ this is done in the paper you are referring to, as can be seen in the lower part of Fig.2 for example (I looked on the arXiv version of it)

As the title of the paper suggests, this is a spontaneous symmetry breaking. That is, it is not at the level of the Hamiltonian (which remains symmetric) but rather at the level of the state of the system. Just like in the 2d Ising model, or in many magnets - the Hamiltonian describing the system is symmetric, and therefore there are two degenerate ground states for example (in the 2d Ising model - all the spins pointing "up" or all of them pointing "down"). However, below some temperature, the system "chooses" one of the states over the other, and breaks the symmetry spontaneously. The fact that such phase-transitions (which are non-analytic behavior of the system) can happen is something that concentrated a lot of theoretical work in the previous century. It is dependent on taking the thermodynamic limit of system size to infinity, which introduces this non-analytic behavior.

One formal way to examine if a symmetry is broken spontaneously is to add a symmetry breaking term and then send it to zero while taking the size of the system to infinity, in the proper order. In your example, it will entail $$ H(\eta) = -\sum_{i,j} \left[\sigma^x_i \sigma^x_j + \sigma^y_i \sigma^y_j + \Delta \sigma^z_i \sigma^z_j\right] + \eta \sum_i \sigma^z_i$$ and then examine the magnetization per spin when the proper limits are taken $$m_z = \lim_{\eta\to 0^+}\lim_{N\to\infty}\frac{1}{N}\sum_i\langle \sigma^z_i\rangle$$ this is done in the paper you are referring to, as can be seen in the lower part of Fig.2 for example (I looked at the arXiv version of it)