1 Algorithm

An ordinary Monte Carlo (MC) algorithm is widely used to compute thermodynamic properties by averaging over the Boltzmann distribution using the Metropolis algorithm [123]. The Metropolis algorithm works by generating trial moves at random and accepting or rejecting each move based on the Boltzman probability density $\exp(-\beta \cdot\Delta \textit{E})$, where $\beta = \frac{1}{\textit{k}_{B}\cdot T}$ ( $\textit{k}_{B}$ is the Boltzamn constant and T is the temperature). The ratio $(\frac{\Delta \textit{E}}{\textit{k}_{B}\cdot T})$ depends on the energy difference between two states, $\Delta\textit{E}=\mathcal{H}^{'}-\mathcal{H}$ ( $\mathcal{H}$ is the Hamiltonian of the system). In the usual Monte Carlo method for magnetic systems we generate our trial moves by drawing a vector $\mathbf{v}$ from an isotropic normal distribution, choosing a spin $\mathbf{S}_i$ at random, adding $\mathbf{S}_i\times\mathbf{v}$ to it and normalizing the result to obtain a trial spin $\mathbf{S}_i'$.

The probability density of the move depends only on the angle between $\mathbf{S}_i$ and $\mathbf{S}_i'$, which ensures reversibility. The variance of the $\mathbf{v}$ distribution controls the average step size and can be chosen at will to improve the ratio of accepted to rejected moves, similarly to the MC angle in Refs. [123,124]. For our Hamiltonian, the energy difference involves only spin $i$ and few neighboring spins to which it is coupled by exchange interaction, so that the decision to accept or reject the move can be made quickly. A sequence of $N$ moves, counting null moves, constitutes a sweep (or one Monte Carlo step); we compute quantities of interest once per sweep to average them.

The constrained Monte Carlo method has been suggested by P. Asselin from Seagate technology [125,126] in relation to the necessity to model high temperature spin dynamics in the heat-assisted magnetic recording. The innovation of the constrained Monte Carlo method is to modify the elementary steps of the random walk so as to conserve the average magnetization direction $\textbf{M}\equiv \bigl(\sum_{i}\textbf{S}_{i}\bigr)/\Vert\sum_{i}\textbf{S}_{i}\Vert$. In this way we sample the Boltzmann distribution over a submanifold of the full phase space. Thus we keep the system out of thermodynamic equilibrium in a controlled manner, while allowing its microscopic degrees of freedom to thermalize.

In the constrained Monte Carlo method the trial moves act on two spins at a time. The extra degrees of freedom allow us to fix $\mathbf{M}$ to any given unit vector, which we take here to be the positive $z$ axis since we can always reduce the problem to this case by means of a global rotation. In detail, the algorithm is the following:

1. Choose a primary spin $\mathbf{S}_i$ and a compensation spin $\mathbf{S}_j,$ not necessarily neighbors.
2. Rotate the primary spin normally similar to normal Metropolis MC algorithm, obtaining a new spin $\mathbf{S}_i'$.
3. Adjust the compensation spin's $x$ and $y$ components to preserve $M_x=M_y=0$,

   $$S_{jx}' = S_{jx} + S_{ix} - S_{ix}'$$ (57)

   $$S_{jy}' = S_{jy} + S_{iy} - S_{iy}'$$ (58)

   
   

4. Adjust the $z$ component,

   $$S_{jz}' = \sign \left(S_{jz}\right)\sqrt{1-S_{jx}'^2-S_{jy}'^2} \quad.$$ (59)

   
   
   If the argument of the square root is negative, stop and take a null move (reject the move).
5. Compute the new magnetization,

   $$M_z'= M_z + S_{iz}' + S_{jz}' - S_{iz} - S_{jz} \quad.$$ (60)

   
   
   If $M_z'\leq 0$, stop and take a null move.
6. Compute the energy difference $\Delta \mathcal{H}=\mathcal{H}'-\mathcal{H}$.
7. Compute the acceptance probability $P$,

   $$P=\min\Biggl[1, \biggl(\frac{M_z'}{M_z}\biggr)^2 \frac{\lvert S_{... ...{\lvert S_{jz}'\rvert } \exp\bigl(-\beta\Delta \mathcal{H}\bigr) \Biggr] \quad.$$ (61)

   
   

8. Accept the move with probability $P$ or take a null move with probability $1-P$.

In effect, we use the compensation spin to project the system back to its admissible manifold. The projection is not orthogonal and does not preserve measure. Consequently the Boltzmann ratio in step 7 is multiplied by a geometric correction, the ratio of two Jacobians (derived in details in Ref. [126]). The reversibility and the ergodicity of this new Monte Carlo method have been rigorously proved, see Ref. [126].