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$,

    $\displaystyle S_{jx}'$ $\displaystyle = S_{jx} + S_{ix} - S_{ix}'$ (57)
    $\displaystyle S_{jy}'$ $\displaystyle = S_{jy} + S_{iy} - S_{iy}'$ (58)

  4. Adjust the $ z$ component,

    $\displaystyle 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,

    $\displaystyle 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$,

    $\displaystyle 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].

Rocio Yanes