If $\hat{m_z}=\frac{1}{N}\sum_i \hat{\sigma^z_i}$ is an order parameter for finite quantum system (transverse Ising model, say), then it will never break the $\mathbb{Z}_2$ symmetry since $\langle\psi_0|\hat{m_z}|\psi_0\rangle=0$ for the ground state $|\psi_0\rangle$ that satisfies the symmetry of system.
This is in fact not a unique problem to quantum systems. But for the classical systems, a quick fix that I'm aware of is replacing the scaling $m_z\sim \tau^\beta$ by $m_z^2\sim \tau^{2\beta}$, where $\tau$ is the control parameter and $\beta$ a critical exponent.
(a further possible step is the Binder cumulant/bimodality coefficient defined as the fourth moment scaling as $\tau^{4\beta}$ normalized by this second moment squared, so that no net scaling with the order parameter remains),
A naive, straightforward generalization to quantum systems is thus working with $\langle \hat{m}_z^2\rangle\sim \tau^{2\beta}$. This seems to be indeed what is done, for example in https://doi.org/10.1103/PhysRevB.87.174302 (eq. 54). But now I saw in a response to a different question here https://physics.stackexchange.com/a/460031/25292 that the scaling of the power of a variable should not be so trivial?
So, my questions are
- Is this procedure of squaring the operator of the order parameter and replacing $\beta\rightarrow 2\beta$ valid?
- In case not, what would be the alternative, for finite-size scaling purposes?
I just saw in an original paper on 1D quantum Ising, https://doi.org/10.1016/0003-4916(70)90270-8 that they used the square root of only the infinite-range correlation instead of from the full squared magnetisation. But it seems unclear what would be the result of this in a finite system.