4 Numerical method

Since we are dealing with a multi-spin particle, the energy potential is multidimensional. Accordingly, in Ref. [74] the technique of Lagrange multiplier was introduced to represent the energy potential in terms of the coordinates of the particle's net magnetization $\vec{\nu} \equiv \sum_i \mathbf{S}_i/\vert\sum_i\mathbf{S}_i\vert$. This technique consists of adding the term $- {\cal N}\vec{\lambda}(\vec{\nu} - \vec{\nu}_0)$, to the total energy Eq. (2.1). This term produces an additional torque that forces the net magnetization to lie along the prescribed direction $\vec{\nu}_0$. The equilibrium state of the spin system is determined by solving the Landau-Lifshitz equation (without the precession term):

$$\frac{d\mathbf{S}_i}{dt}=-\alpha \mathbf{S}_i\times\left(\mathbf{S}_i\times\mathbf{H}_i\right).$$ (41)

The effective field $\mathbf{H}_i=-\delta \mathcal{H}/\delta \mathbf{S}_i$ now contains an additional term due to the Lagrange constraint. For the three components of the Lagrange parameter $\vec{\lambda}$ the set of equations (2.24) is augmented by the equations $d\vec{\lambda}/dt=\delta \mathcal{H}/\delta \vec{\lambda}$. The stationary points found with this method are also stationary points of the Hamiltonian (2.1), since for these points $\vec{\lambda}=0$. However, if the system has many metastable states, only part of these points, compatible with the behavior assumed by the direction of the net magnetization is determined. More precisely, in individual small particles where the exchange interaction is dominating, the deviations from the collinear state are small and thereby the individual spins adiabatically adjust to the net direction when the latter is rotated. In this case, it is possible to define a net magnetization and parameterize it, e.g., in the spherical system of coordinates as $\mathbf{\nu_0}(\theta,\varphi)$. The advantage of this technique is that it can produce highly non-collinear multi-dimensional stationary points [74,87,99]. Accordingly, in the present work and unlike the studies presented in Ref. [14], spin non-collinearities are taken into account. Moreover, in order to check the correct loci of the saddle points, we computed the eigenvalues and gradient of the Hessian matrix associated with the Hamiltonian (2.1).