# Spin Hamiltonian with large fourth order terms: Triple well potentials and
Bloch sphere visualization

D S Lohr-Robles1, M Grether2, E López-Moreno2, P O Hess3,4 1Independent
researcher, Mexico  2Facultad de Ciencias, Universidad Nacional Autónoma de
México, 04510 Mexico-City, Mexico  3Instituto de Ciencias Nucleares,
Universidad Nacional Autónoma de México, A.P. 70-543, 04510 Mexico-City,
Mexico  4Frankfurt Institute for Advanced Studies, J. W. von Goethe
University, Hessen, Germany<EMAIL_ADDRESS>
###### Abstract

We present a study of a general spin Hamiltonian with terms up to fourth
order. With the coherent states the semiclassical potential is obtained and
with catastrophe theory its parameter space is constructed. When the fourth
order parameters are large enough the parameter space has regions where the
semiclassical potential has three wells. By applying an oscillating magnetic
field a trajectory in parameter space crosses the Maxwell set multiple times
resulting in many ground state quantum phase transitions. Using the coherent
states we are able to visualize the localization of the ground state on the
Bloch sphere as the magnetic field is varied.

Conflict of Interest declaration: The authors declare that they have NO
affiliations with or involvement in any organization or entity with any
financial interest in the subject matter or materials discussed in this
manuscript.

Keywords: Quantum phase transitions, algebraic models, semiclassical
approximation

##  1 Introduction

The use of spin Hamiltonian to model many body systems has experienced a wide
variety of applications in different areas of physics. By an appropriate
fitting of the interacting parameters the phenomena of interest may be
accurately described, assigning physical meaning to each of the terms. For
this reason, theoretical work when the Hamiltonian has higher order
contributing terms can be of use, in particular when the value of the
corresponding parameters are large enough to induce a meaningful change, and
how this contribution affects the phase transitions of the system. The physics
of single molecule magnets (SMMs) can be accurately described using a spin
Hamiltonian [1]. The applications of SMMs in new technologies as qubits
remains a possibility [2, 3, 4, 5, 6, 7].

Spin Hamiltonians have also been used in the description of quantum phase
transitions (QPTs) and excited states quantum phase transitions (ESQPTs) as a
function of the Hamiltonian parameters [8], and they have also been applied to
model systems of very large spin in ferritin proteins [9]. As a result these
models are very useful as a way to test new phenomena. The motivation of the
present contribution is to study phase transitions in SMMs and the search of a
tractable manner to manipulate them.

The present paper is structured as follows: In section 2 the spin Hamiltonian
is introduced with terms up to fourth order, the coherent states are defined
and the parameter space of the semiclassical potential is obtained. In section
3 the master equation that defines the population of states as a function of
time is introduced, together with the concept of relaxation time of the
material. In section 4 the concept of fidelity and fidelity susceptibility and
their relation to QPTs is discussed. A series of results for three different
examples are presented in section 5. In section 6, hypothetical examples of
systems with large values of the fourth order parameters are considered, which
lead to systems with three stability wells in the semiclassical potentials;
two of these examples are presented and studied. The concept of the Bloch
sphere is introduced as a way to visualize and describe the localization of
the eigenstates of the Hamiltonian with the intention to construct and
manipulate qubits. Finally, in section 7 conclusions are drawn and possible
future work is discussed.

##  2 Spin Hamiltonian

We consider the following expression for a spin Hamiltonian, designed to
describe magnetic molecules, with operators up to fourth order:

$H$ $\displaystyle=$
$\displaystyle-\frac{g\mu_{B}}{S}\vec{{\mbox{\boldmath$B$}}}\cdot\vec{{\mbox{\boldmath$S$}}}+\frac{1}{S(2S-1)}\left(D\left({\mbox{\boldmath$S$}}^{2}_{z}-\frac{1}{3}S(S+1)\right)+E({\mbox{\boldmath$S$}}^{2}_{x}-{\mbox{\boldmath$S$}}^{2}_{y})\right)$
(2) $\displaystyle+\frac{{\mbox{\boldmath$H$}}_{4}}{S(2S-1)(2S-2)(2S-3)}$

where $g$ is the Landé factor, $\mu_{B}$ is the Bohr magneton, $D$ and $E$ are
the second order anisotropy constants of the molecule and
${\mbox{\boldmath$H$}}_{4}$ indicates fourth order anisotropy terms. The first
order interaction is given by the Zeeman term in (2), with
$\vec{{\mbox{\boldmath$B$}}}=(B_{x},B_{y},B_{z})$ an external magnetic field
and $\vec{{\mbox{\boldmath$S$}}}=(S_{x},S_{y},S_{z})$ the spin vector, with
$S_{i}$, $i=x,y,z$, the spin components operators. The explicit form of the
fourth order terms is given by:

$\displaystyle{\mbox{\boldmath$H$}}_{4}$ $\displaystyle=$ $\displaystyle
B_{4}^{0}{\mbox{\boldmath$O$}}_{4}^{0}+B_{4}^{2}{\mbox{\boldmath$O$}}_{4}^{2}+B_{4}^{3}{\mbox{\boldmath$O$}}_{4}^{3}+B_{4}^{4}{\mbox{\boldmath$O$}}_{4}^{4}$
(3) $\displaystyle=$ $\displaystyle
B_{4}^{0}\left(35{\mbox{\boldmath$S$}}_{z}^{4}+(25-30S(S+1)){\mbox{\boldmath$S$}}_{z}^{2}+3S^{2}(S+1)^{2}-6S(S+1)\right)$
(6)
$\displaystyle+\frac{B_{4}^{2}}{4}\left((7{\mbox{\boldmath$S$}}_{z}^{2}-S(S+1)-5)({\mbox{\boldmath$S$}}_{+}^{2}+{\mbox{\boldmath$S$}}_{-}^{2})+({\mbox{\boldmath$S$}}_{+}^{2}+{\mbox{\boldmath$S$}}_{-}^{2})(7{\mbox{\boldmath$S$}}_{z}^{2}-S(S+1)-5)\right)$
$\displaystyle+\frac{B_{4}^{3}}{4}\left({\mbox{\boldmath$S$}}_{z}({\mbox{\boldmath$S$}}_{+}^{3}+{\mbox{\boldmath$S$}}_{-}^{3})+({\mbox{\boldmath$S$}}_{+}^{3}+{\mbox{\boldmath$S$}}_{-}^{3}){\mbox{\boldmath$S$}}_{z}\right)+\frac{B_{4}^{4}}{2}({\mbox{\boldmath$S$}}_{+}^{4}+{\mbox{\boldmath$S$}}_{-}^{4}),$

where ${\mbox{\boldmath$O$}}_{q}^{k}$ are the Steven operators [1, 10, 11] and
$B_{q}^{k}$ the corresponding fourth order anisotropy constants. It is
important to notice that the values of these constants are smaller than that
of the second order ones, and the non-zero values are dependent on the
symmetric properties of the material: tetragonal symmetry ($k=0,4$),
orthorhombic symmetry ($k=0,2,4$), and trigonal symmetry ($k=0,3$) [10, 11].
To each of the terms in the Hamiltonian we add a factor of the form
$2\prod_{q}(2S-q+1)^{-1}$ for the $q$-th order interaction term. This is to
ensure that the semiclassical potential calculated in the next subsection is
independent of $S$.

Using the complete basis $|SM\rangle$ of the eigenvectors of the $S_{z}$
operator, satisfying the eigenvalue equation
$S_{z}|SM\rangle=M\hbar|SM\rangle$, the Hamiltonian matrix of (2) is
constructed and its eigenvalues $E_{k}$, which satisfy

${\mbox{\boldmath$H$}}|\psi_{k}\rangle=E_{k}|\psi_{k}\rangle,$ (7)

are obtained by diagonalization. The eigenvectors $|\psi_{k}\rangle$ are
expressed as linear combinations in the $|SM\rangle$ basis as:

$|\psi_{k}\rangle=\sum_{M=-S}^{S}c_{k,M}|SM\rangle,$ (8)

where the square of the absolute value of the coefficients $c_{k,M}$ represent
the probability of finding the state $|\psi_{k}\rangle$ with spin projection
$M$, with $M=-S,-S+1,\ldots,S-1,S$. In the following we will simplify the
notation defining $|SM\rangle\equiv|M\rangle$.

The Hamiltonian of the system written in (2) depends on nine parameters: Three
free parameters for the interaction with the magnetic field and six
anisotropic constants (two second order and four fourth order). The non-free
parameters are adjusted to experiment. We want to study the effect of these
parameters in the structure of the eigenstates of the Hamiltonian. Using the
$\mathrm{SU}(2)$ coherent states [12] we can obtain the semiclassical
potential of the system as the expectation value of the Hamiltonian in the
coherent state basis. We will find that by studying the critical points of the
semiclassical potential we can obtain information about the behaviour of the
eigenvalues and eigenvectors. Catastrophe theory is a very useful method [13,
14, 15] for the categorization of functions that depend on parameters, as it
presents a classification of the stable singularities of a potential functions
leading to the construction of separatrices in parameter space, denoting the
different phases of the system, identifiable as elementary catastrophes: fold,
cusp, swallowtail, butterfly, etc., and their respective regions of stability.
The catastrophe theory has been applied with great success to nuclear and
particle physics [16, 17].

###  2.1 Coherent states and semiclassical potential

The semiclassical potential is obtained by calculating the expectation value
of the Hamiltonian in the atomic coherent state basis:
$V(\theta,\phi)=\langle\zeta|{\mbox{\boldmath$H$}}|\zeta\rangle$, where the
coherent state $|\zeta\rangle$ is defined as a rotation of lowest weighted
state $|-S\rangle$ of angle $\theta$ about an axis
$\hat{n}=(\sin\phi,-\cos\phi,0)$ in angular momentum space [12]:

$\displaystyle|\zeta\rangle$ $\displaystyle=$ $\displaystyle
R(\theta,\phi)|-S\rangle$ (9) $\displaystyle=$
$\displaystyle(1+|\zeta|^{2})^{-S}e^{\zeta S_{+}}|-S\rangle$ (10)
$\displaystyle=$
$\displaystyle\sum_{M=-S}^{S}\left(\frac{(2S)!}{(S+M)!(S-M)!}\right)^{1/2}\left(\cos\frac{\theta}{2}\right)^{S-M}\left(e^{-\mathrm{i}\phi}\sin\frac{\theta}{2}\right)^{S+M}|M\rangle$
(11)

with $\zeta=e^{-i\phi}\tan(\theta/2)$.

The semiclassical potential is explicitly given by:

$\displaystyle V(\theta,\phi;D,E,B_{i})$ $\displaystyle=$
$\displaystyle\frac{D}{12}(1+3\cos 2\theta)+\frac{E}{2}\cos
2\phi\sin^{2}\theta$ (15)
$\displaystyle-g\mu_{B}\left(-B_{z}\cos\theta+B_{x}\cos\phi\sin\theta+B_{y}\sin\phi\sin\theta\right)$
$\displaystyle+\frac{1}{8}\left(\frac{B_{4}^{0}}{8}(35\cos 4\theta+20\cos
2\theta+9)+\frac{B_{4}^{2}}{2}(7\cos 2\theta+5)\cos
2\phi\sin^{2}\theta\right.$ $\displaystyle\left.-B_{4}^{3}\cos
3\phi\cos\theta\sin^{3}\theta+B_{4}^{4}\cos 4\phi\sin^{4}\theta\right),$

which is a function of two angular variables $(\theta,\phi)$ and independent
of the value of spin $S$ of the system. Therefore, the parameter spaces
obtained in these sections are valid for all values of $S$. Because our
interest lies on how the changes of parameters can lead the system to change
from one phase to another, we will find that very telling results, e.g.
indications of QPTs, can be obtained even for small values of $S$.

As a starting point we will restrict ourselves to the case when the magnetic
field is constrained in the $xz$-plane: $\vec{B}=(B_{x},0,B_{z})$, as
previously done in [18]. We can see in (15) that when $B_{y}=0$, the values
$\phi_{c}=0,\pi$, are critical points for all values of the parameters. Thus,
we are able to substitute these values in (15) and focus on the study of the
following one-dimensional potential function:

$\displaystyle V(\theta,\phi_{c};r_{i})$ $\displaystyle=$ $\displaystyle
r_{1}\cos\phi_{c}\sin\theta+r_{2}\cos\theta+r_{3}\cos 2\theta+r_{4}\cos
4\theta$ (17) $\displaystyle+r_{5}\cos\phi_{c}(2\sin 2\theta-\sin 4\theta),$

where we defined the new parameters $r_{i}$ as:

$\displaystyle r_{1}$ $\displaystyle=$ $\displaystyle-g\mu_{B}B_{x}$ (18)
$\displaystyle r_{2}$ $\displaystyle=$ $\displaystyle g\mu_{B}B_{z}$ (19)
$\displaystyle r_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(D-E)+\frac{1}{16}(5B_{4}^{0}+B_{4}^{2}-B_{4}^{4})$
(20) $\displaystyle r_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{64}(35B_{4}^{0}-7B_{4}^{2}+B_{4}^{4})$ (21)
$\displaystyle r_{5}$ $\displaystyle=$ $\displaystyle-\frac{1}{64}B_{4}^{3}.$
(22)

Inversely, once the parameters $r_{i}$ are known, the magnetic field
components and $D$, $E$ can be deduced.

###  2.2 Bifurcation and Maxwell sets

The bifurcation set is a subspace in parameter space where critical points of
the potential function begin to emerge. Thus, in its vicinity we can find two
phases (regions), in one there exist a stability point $\theta_{c}$, while in
the other there is no such point. The bifurcation set can be found considering
the critical manifold, which is the hypersurface of critical points
$\theta_{c}$ satisfying:

$\left.\frac{d}{d\theta}V(\theta,\phi_{c};r_{i})\right|_{\theta=\theta_{c}}=0$
(23)

spanned by a continuous variation of the parameters $r_{i}$. The set of points
where the mapping of the critical manifold to the parameter space is singular
is defined as the bifurcation set. A straightforward calculation of this
singular mapping result in the following set of parametric functions:

$\displaystyle r_{1}(x,\phi_{c};r_{3},r_{4},r_{5})$ $\displaystyle=$
$\displaystyle r_{3}\cos\phi_{c}(3\sin x-\sin 3x)+2r_{4}\cos\phi_{c}(5\sin
3x-3\sin 5x)$ (25) $\displaystyle-6r_{5}(\cos x-2\cos 3x+\cos 5x)$
$\displaystyle r_{2}(x,\phi_{c};r_{3},r_{4},r_{5})$ $\displaystyle=$
$\displaystyle-r_{3}(3\cos x+\cos 3x)-2r_{4}(5\cos 3x+3\cos 5x)$ (27)
$\displaystyle+4r_{5}\cos\phi_{c}\sin x(2+7\cos 2x+3\cos 4x),$

where $x=\theta_{c}$ is the critical point and satisfies (23). These
parametric functions allow us to draw the bifurcation set in $(r_{1},r_{2})$
parameter space for a given set of values $(r_{3},r_{4},r_{5})$. One can also
solve these equations to obtain the bifurcation set in any other
$(r_{i},r_{j})$ parameter space of interest.

The Maxwell set is the subspace in parameter space where two or more extrema
(minima or maxima) of the potential function have the same value, i.e.
$V(\theta_{1},\phi_{c};r_{i})=V(\theta_{2},\phi_{c};r_{i})$; thus in its
vicinity two phases exist, in one $V(\theta_{1})>V(\theta_{2})$, while in the
other $V(\theta_{1})<V(\theta_{2})$, i.e. there is a change of the dominant
stability point (phase). Similarly to the bifurcation set, the Maxwell set can
be found by considering the singular mapping of a hypersurface to the
parameter space. The particular hypersurface is defined as the set of roots
$\theta_{i}$, such that $V(\theta_{i},\phi_{c};r_{i})+V_{0}=0$, where $V_{0}$
is an arbitrary real number, spanned by a continuous variation of the
parameters $r_{i}$. Then one has to find the values of $r_{i}$ such that the
above description is true for two values $\theta_{1}$ and $\theta_{2}$.
Performing the calculation is very involved. Fortunately, using the routines
of the software MATHEMATICA we could resolve the problem and we found the
following set of functions for the Maxwell set:

$\displaystyle r_{1}(\theta_{1},\theta_{2};r_{4},r_{5})$
$\displaystyle=8\sec\phi_{c}\sin\theta_{1}\sin\theta_{2}(\sin\theta_{1}+\sin\theta_{2})\Bigg{(}2r_{4}(1-\cos(\theta_{1}+\theta_{2}))$
(28) $\displaystyle\quad+r_{5}\cos
3\phi_{c}\bigg{(}2\cot(\theta_{1}+\theta_{2})-\cos
2(\theta_{1}+\theta_{2})\csc(\theta_{1}+\theta_{2})\bigg{)}\Bigg{)}$ (29)
$\displaystyle r_{2}(\theta_{1},\theta_{2};r_{4},r_{5})$
$\displaystyle=2\cos\frac{\theta_{1}-\theta_{2}}{2}\Bigg{(}32r_{4}\cos\theta_{1}\cos\theta_{2}\cos^{3}\frac{\theta_{1}+\theta_{2}}{2}\Bigg{)}$
(30) $\displaystyle\quad-r_{5}\cos
3\phi_{c}\csc\frac{\theta_{1}+\theta_{2}}{2}\bigg{(}2-2\cos 2\theta_{1}-2\cos
2\theta_{2}-\cos(\theta_{1}+\theta_{2})$ (31) $\displaystyle\quad-2\cos
2(\theta_{1}+\theta_{2})-\cos
3(\theta_{1}+\theta_{2})-\cos(3\theta_{1}+\theta_{2})-\cos(\theta_{1}+3\theta_{2})\bigg{)}$
(32) $\displaystyle r_{3}(\theta_{1},\theta_{2};r_{4},r_{5})$
$\displaystyle=-2r_{4}\bigg{(}3\cos 2\theta_{1}+3\cos
2\theta_{2}+4\cos(\theta_{1}+\theta_{2})\bigg{)}$ (33)
$\displaystyle\quad+\frac{r_{5}}{2}\cos
3\phi_{c}\csc\frac{\theta_{1}+\theta_{2}}{2}\sec\frac{\theta_{1}+\theta_{2}}{2}\bigg{(}2\cos(\theta_{1}+\theta_{2})-4\cos(2\theta_{1}+2\theta_{2})$
(34)
$\displaystyle\quad-3\big{(}\cos(3\theta_{1}+\theta_{2})+\cos(\theta_{1}+3\theta_{2})\big{)}\bigg{)}$
(35)

and we are able to draw the Maxwell set as a parametric function in parameter
space $(r_{1},r_{2})$ given the critical points $(\theta_{1},\theta_{2})$ that
satisfy (28) for a given set of values $(r_{3},r_{4},r_{5})$.

In the following applications, examples of phases are shown, and the
corresponding separatrices of the bifurcation and Maxwell sets will be
illustrated.

##  3 Magnetization and master equation

The labelling of the eigenstates of the Hamiltonian (2) is an important issue
to address. In cases when the $D$ parameter is the most relevant the low
energy states have a dominant coefficient $c_{k,M}$ in the expression (8),
while the other coefficients are small. In this cases it is possible to
maintain the labelling of states of the spin projections, which is required
when studying transitions between states in the interaction of the system with
a magnetic field. In this framework the states labelled $M=\pm S$, are the
lowest states in their respective potential wells.

The population $p_{m}$ of the state labelled $m$ at a time $t$ is described by
the master equation [1]:

$\frac{dp_{m}(t)}{dt}=\sum_{m^{\prime}}\big{(}\gamma_{m^{\prime}m}p_{m^{\prime}}(t)-\gamma_{mm^{\prime}}p_{m}(t)\big{)}$
(36)

with $\gamma_{mm^{\prime}}$ the transition rates of going from the state
$m^{\prime}$ to $m$, and are given by

$\displaystyle\gamma_{mm^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{3}{\pi\hbar^{4}\rho
c_{s}^{5}}\frac{(E_{m^{\prime}}-E_{m})^{3}}{e^{(E_{m^{\prime}}-E_{m})/k_{B}T}-1}\Big{(}|D_{1}|^{2}\big{(}|\langle\psi_{m^{\prime}}|S_{+}^{2}|\psi_{m}\rangle|^{2}+|\langle\psi_{m^{\prime}}|S_{-}^{2}|\psi_{m}\rangle|^{2}\big{)}$
(38)
$\displaystyle+|D_{2}|^{2}\big{(}|\langle\psi_{m^{\prime}}|\\{S_{+},S_{z}\\}|\psi_{m}\rangle|^{2}+|\langle\psi_{m^{\prime}}|\\{S_{-},S_{z}\\}|\psi_{m}\rangle|^{2}\big{)}\Big{)}$

where $\rho$ is the density of the material, $c_{s}$ is the velocity of sound
in the material, $D_{1}$ and $D_{2}$ the spin-phonon couplings parameters, and
$\\{\cdot,\cdot\\}$ is the anti-commutator [1, 19].

The differential equation in (36) can be written in matrix form as:

$\frac{d\vec{p}(t)}{dt}=\vec{p}(t){\mbox{\boldmath$G$}}$ (39)

where $G$ is the transition rate matrix, and its respective matrix elements
are given by:

$G_{mm^{\prime}}=\gamma_{mm^{\prime}}-\delta_{mm^{\prime}}\sum_{k}\gamma_{mk}.$
(40)

In order to solve the master equation (36) one has to take into account that
the variation of the magnetic field in (2) is time dependent: $|\vec{B}|=\nu
t+B_{0}$, with $\nu=dB/dt$ the sweeping rate of magnetization and $B_{0}$ the
initial magnetization. The study of the master equation has been considered
for various SMMs in order to model the magnetization hysteresis loops [20, 21,
19, 22]. Monte Carlo methods have proven to be a strong computational
technique when dealing with time dependent transition rates [23, 24, 25, 26,
27]. In particular, in [22] they used a kinetic Montecarlo variable size step
method to solve the master equation, which we will follow in the present
paper.

The magnetization hysteresis loops are obtained using the solutions of the
master equation (36) as:

$M=\sum_{m}p_{m}(t)\langle\psi_{m}|S_{z}|\psi_{m}\rangle.$ (41)

###  3.1 Relaxation time and relaxation rate

Measurements of the relaxation time and rate as a function of the magnetic
field have been of considerable importance in the understanding of SMMs [28,
29, 30, 31, 32]. The relaxation time $\tau$ of the material as a function of
the magnetic field can be obtained with the eigenvalues of the transition rate
matrix in (40) as done in [1, 33, 34, 35]. It is proven in [1] that one of the
eigenvalues is equal to zero and corresponds to the equilibrium state, while
the others are negative. The relaxation time is then identified as the
reciprocal of the smallest non-zero eigenvalue $g_{min}$ [34]:

$\tau=-\frac{1}{g_{min}}.$ (42)

The relaxation time tells us how long does the magnetization due to the
alignment of the molecules persists in the presence of the magnetic field.

The relaxation rate $\Gamma$ is in turn the reciprocal of the relaxation time
[1, 34]:

$\Gamma=\frac{1}{\tau}$ (43)

In [34] the peaks of the relaxation rate as a function of time are fitted with
a Lorentzian function.

##  4 Fidelity and fidelity susceptibility

Quantum fidelity is a measure of how much does a quantum state, that depends
on a parameter, resembles itself with a small variation of the parameter [36],
e.g. the magnetic field amplitude. This quantity has proven to be very useful
for the identification of QPTs, characterized by a sudden drop in fidelity at
the critical point of the parameter, and has been used successfully to detect
QPTs in quantum many-body systems [36, 37, 38, 39, 40, 41, 42, 43]. In order
to calculate the quantum fidelity of a state as a function of the parameter
one needs to introduce an arbitrary small value $dB$.

The fidelity of a state $|\psi_{k}\rangle$ is defined as [36]:

$F_{k}=|\langle\psi_{k}(B_{i}-dB_{i})|\psi_{k}(B_{i}+dB_{i})\rangle|^{2}$ (44)

where $B_{i}$ is the magnetic field parameter and $dB_{i}$ represents a small
increment. Using this expression one is able to calculate the value of the
fidelity of the $k$-th state as a function of a varying magnetic field.

The quantum fidelity susceptibility $\chi_{F_{k}}$ is defined as the second
order coefficient of the Taylor series expansion of the quantum fidelity in
(44) about $dB=0$ [36]. It contains all the information of the quantum
fidelity and has the advantage that it is independent of the arbitrary value
of $dB$, i.e. it does not depend on any external value and it is therefore
more suitable for calculations.

The fidelity susceptibility can be explicitly written as [36]

$\chi_{F_{k}}=2\sum_{m\neq
k}\frac{|\langle\psi_{m}|S_{z}|\psi_{k}\rangle|^{2}}{(E_{m}-E_{k})^{2}}$ (45)

where $S_{z}$ is related to the parameter of the interaction, i.e. the
magnetic field.

##  5 Results

In this section we will present results of the physical phenomena describe in
the previous sections for three different test cases of magnetic molecules: a)
$\mathrm{Fe}_{8}$ SMM, b) $\mathrm{Fe}_{4}$ SMM, c) Arbitrary parameters, not
related to a specific magnetic molecule, with a large fourth order term.

###  5.1 $\mathrm{Fe}_{8}$ SMM

We start with the parameter values given in the review [10] and modify them a
little in order to better fit the hysteresis plots in [44]: $S=10$,
$D/k_{B}=-0.295$ K, $E/k_{B}=0.056$ K, $B_{4}^{0}/k_{B}=1.15\times 10^{-6}$ K,
$B_{4}^{2}/k_{B}=-1.15\times 10^{-6}$ K, $B_{4}^{4}/k_{B}=-2.18\times 10^{-5}$
K, $B_{x}=0.02$ T, $B_{y}=B_{4}^{3}=0$. In this case we considered the
magnetic field to be not completely aligned with the $z$-axis and the easy
magnetization axis of the molecules, but making a small $3^{\circ}$ angle with
respect to it; this was done to recreate the hysteresis steps in the presence
of a very small transverse field.

The parameters in the transition rates $\gamma_{mm^{\prime}}$ in (38) used
are: $3/(\pi\hbar^{4}\rho c_{s}^{5})=3.13\times
10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$, $D_{1}=D_{2}=0.26$ K, and $T=0.04$ K.

In figure 1(a) we show the bifurcation set, and the horizontal dashed line
corresponds to the parameter values of $\mathrm{Fe}_{8}$. When the bifurcation
set is crossed the avoided level crossings of higher energy levels start to
occur, as a consequence of the double well structure of the semiclassical
energy surface. At zero magnetic field both potential wells are at the same
depth and there is an avoided level crossing of the ground state, which can be
seen in figure 1(d). In figures 1(b) and 1(e) we plotted the relaxation time
and relaxation rate as functions of the magnetic field, respectively. The
drops of the relaxation time, which are peaks in the relaxation rate, occur at
the avoided level crossings of the eigenstate with label $m=-10$, as seen in
figure 1(d). A similar behaviour occur with the fidelity and fidelity
susceptibility seen in figures 1(c) and 1(f), respectively. In figure 1(g) we
plotted the hysteresis loop for two values of the sweeping rate of the
magnetic field in the $z$ direction. These results can be compared with the
experimental results obtained in figure 2 in [44]. The transition rates
defined in (38) for the master equation in this model does not adequately
reproduce the transitions at zero magnetic field for pure tunneling, i.e. at
low temperature, as described in [22], and to solve this issue and additional
constant correction term is added to the transition rates
$\gamma_{mm^{\prime}}$ for these cases. Following the methodology described in
[22], we added the constant correction term
$\gamma_{t}=0.002\;\mathrm{s^{-1}}$ in order to obtain the hysteresis in
figure 1(g). The location of the steps in the hysteresis loops also coincide
with the location of the drops of the relaxation time.

Figure 1: (a) Bifurcation set in $(r_{2}/k_{B},r_{3}/k_{B})$ parameter space
and the dotted line represents the parameters for the $\mathrm{Fe}_{8}$
example: $S=10$, $D/k_{B}=-0.295$ K, $E/k_{B}=0.056$ K,
$B_{4}^{0}/k_{B}=1.15\times 10^{-6}$ K, $B_{4}^{2}/k_{B}=-1.15\times 10^{-6}$
K, $B_{4}^{4}/k_{B}=-2.18\times 10^{-5}$ K, $B_{x}=0.01$ T,
$B_{y}=B_{4}^{3}=0$; and $3/(\pi\hbar^{4}\rho c_{s}^{5})=3.13\times
10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$, $D_{1}=D_{2}=0.26$ K, and $T=0.04$ K.
(b) Relaxation time, (c) fidelity, (d) energy levels, (e) relaxation rate, (f)
fidelity susceptibility, (g) hysteresis loop for sweeping rates: $0.007$ T/s
and $0.056$ T/s, as a functions of $B_{z}$. To calculate the fidelity we used
$dB_{z}=0.001$.

###  5.2 $\mathrm{Fe}_{4}$ SMM

We use the parameters obtained in [45], with the addition of a fourth order
term $B_{4}^{3}$ permitted by the tetragonal symmetry of the molecule to
better fit the experimental hysteresis loops. The parameters used for
$\mathrm{Fe}_{4}$ were: $S=5$, $D/k_{B}=-0.601$ K, $E/k_{B}=0.024$ K,
$B_{4}^{0}/k_{B}=2.88\times 10^{-5}$ K, $B_{4}^{3}/k_{B}=0.0004$ K,
$B_{x}=B_{y}=B_{4}^{2}=B_{4}^{4}=0$.

The parameters in the transition rates $\gamma_{mm^{\prime}}$ in (38) used
are: $3/(\pi\hbar^{4}\rho c_{s}^{5})=4.1\times
10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$, $D_{1}=D_{2}=0.28$ K, and $T=0.04$ K.
For this case we also add a constant correction term $\gamma_{t}$ to the
transition rate to enhance the tunnelling between states $m=\pm 5$ and $m=\pm
4$, as was done in [22]. Depending on the magnetic sweeping rate we have:
$\gamma_{t}=0.1\;\mathrm{s}^{-1}$ for the sweeping rate
$0.001\;\mathrm{T}\mathrm{s}^{-1}$ and $\gamma_{t}=0.6\;\mathrm{s}^{-1}$ for
the sweeping rate $0.017\;\mathrm{T}\mathrm{s}^{-1}$.

In figure 2 we show the same results as in figure 1 but for the
$\mathrm{Fe}_{4}$ SMM. In this case we can see in figures 2(b) and 2(e) the
consequences of the constant correction term for the transition rate between
states $m=-5$ and $m=5$ at zero magnetic field, were the value of the
relaxation time is very small in the vicinity of $B_{z}=0$, while the
relaxation rate has an artificial peak in that vicinity. In figure 2(g) we
plotted the hysteresis loop for two values of the sweeping rate of the
magnetic field in the $z$ direction. These results can be compared with the
experimental results obtained in figure 3 in [45].

Figure 2: (a) Bifurcation set in $(r_{2}/k_{B},r_{3}/k_{B})$ parameter space
and the dotted line represents the parameters for the $\mathrm{Fe}_{4}$
example: $S=5$, $D/k_{B}=-0.601$ K, $E/k_{B}=0.024$ K,
$B_{4}^{0}/k_{B}=2.88\times 10^{-5}$ K, $B_{4}^{3}/k_{B}=0.0004$ K,
$B_{x}=B_{y}=B_{4}^{2}=B_{4}^{4}=0$; and $3/(\pi\hbar^{4}\rho
c_{s}^{5})=4.1\times 10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$,
$D_{1}=D_{2}=0.28$ K, and $T=0.04$ K. (b) Relaxation time, (c) fidelity, (d)
energy levels, (e) relaxation rate, (f) fidelity susceptibility, (g)
hysteresis loop for sweeping rates: $0.001$ T/s and $0.017$ T/s, as a
functions of $B_{z}$. To calculate the fidelity we used $dB_{z}=0.001$.

###  5.3 A case of an arbitrary set of parameters with large fourth order term

A system with arbitrary parameters is considered in order to investigate the
effects of large values of fourth order terms. Here we consider the
parameters: $S=5$, $D/k_{B}=-0.5$ K, $E/k_{B}=0$ K,
$B_{4}^{4}/k_{B}=-1.5\times 10^{-2}$ K, $B_{x}=0.001$ T,
$B_{y}=B_{4}^{0}=B_{4}^{2}=B_{4}^{3}=0$.

The parameters in the transition rates $\gamma_{mm^{\prime}}$ in (38) used
are: $3/(\pi\hbar^{4}\rho c_{s}^{5})=1.831\times
10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$, $D_{1}=D_{2}=0.5$ K, and $T=0.1$ K.

In figure 3 we show the parameter space and we can see that in this case the
horizontal dashed line corresponding to the parameter values crosses a
butterfly catastrophe structure of the bifurcation set. The energy levels are
shown in figure 3(d) were the ground state energy has strong avoided level
crossings for non-zero values of the magnetic field. This is reflected in the
relaxation time and relaxation rate in figures 3(b) and 3(e), respectively,
were in the insets we can see the respective drops and peaks about
$B_{z}\approx 5$ T and $B_{z}\approx 5.6$ T. This structure can also be seen
in the fidelity and fidelity susceptibility in figures 3(c) and 3(f),
respectively, were the drops in the fidelity and the peaks of the fidelity
susceptibility appear more clearly.

Figure 3: (a) Bifurcation set in $(r_{2}/k_{B},r_{3}/k_{B})$ parameter space
and the dotted line represents the parameters in the test case: $S=5$,
$D/k_{B}=-0.5$ K, $E/k_{B}=0$ K, $B_{4}^{4}/k_{B}=-1.5\times 10^{-2}$ K,
$B_{x}=0.001$ T, $B_{y}=B_{4}^{0}=B_{4}^{2}=B_{4}^{3}=0$; and
$3/(\pi\hbar^{4}\rho c_{s}^{5})=1.831\times
10^{3}\;\mathrm{K}^{-5}\mathrm{s}^{-1}$, $D_{1}=D_{2}=0.5$ K, and $T=0.1$ K.
(b) Relaxation time, (c) fidelity, (d) energy levels, (e) relaxation rate, (f)
fidelity susceptibility as a functions of $B_{z}$. To calculate the fidelity
we used $B_{z}=0.001$.

##  6 Triple well potentials

The spin Hamiltonian in consideration has two free parameters, namely the $x$
and $z$ components of magnetic field, while the parameters of higher order
terms are fixed by the system of study. However it should be noted that the
methods described above allow for the study of QPTs in parameter spaces
involving free parameters of higher order terms in a similar fashion.

In this section we will perform a theoretical study and explore the complete
parameter space for regions of interesting structural stability. Instead of
freely varying the Hamiltonian parameters, by constructing first the
semiclassical potential, and the separatrices in parameter space, we can have
an idea of the behaviour and structure of the energy levels of the
Hamiltonian, and it is this aspect where the usefulness and advantage of the
methodology resides.

Two example cases are presented next where triple well potentials are found:

  * Case I:

Parameters: $S=10$, $D=-0.5$ K, $E=0.04$ K, $B_{4}^{4}=-0.007$ K,
$B_{4}^{0}=B_{4}^{2}=B_{4}^{3}=0$. In figure 4(a) we show the parameter space
$(r_{1},r_{2})$ with the bifurcation set in green and the Maxwell set in red
(minima) and dark red (maxima). This particular shape of the bifurcation set
results from the overlap of two butterfly catastrophes one for each of the
critical points $\phi_{c}=0,\pi$. This overlap is a consequence of the large
value of the fourth order parameters $B_{4}^{4}$, which result in a large
value of $r_{4}$.

  * Case II:

Parameters: $S=10$, $D=0.4$ K, $E=-0.03$ K, $B_{4}^{0}=-0.00012$ K,
$B_{4}^{2}=B_{4}^{3}=B_{4}^{4}=0$. In this case the parameter space consists
of two cusps along the $r_{1}$-axis and two butterfly catastrophes along the
$r_{2}$-axis. In figure 4(b) we show one of the butterfly separatrix in
parameter space $(r_{1},r_{2})$. At the intersection of two Maxwell set lines
we find that the semiclassical potential has three equally deep potential
wells.

Figure 4: (a) and (b) Parameter space $(r_{1},r_{2})$ for the cases I and II,
respectively. The bifurcation set is depicted in green and Maxwell set in red
(minima) and dark red (maxima). (a) The varying magnetic field
$\vec{B}=(15.63\cos\omega t,0,15.63\sin\omega t)$ is shown as a black circle
with an arrow. The white dots where the Maxwell set for the global minimum is
crossed are labelled $1$ to $4$. (b) The varying magnetic field
$\vec{B}=(4.47\cos\omega t,0,4.47\sin\omega t-32.9)$ is shown as a black
circle with an arrow. The white dots where the Maxwell set for the global
minimum is crossed are labelled $1$ to $3$, and the white dot labelled $4$
correspond to the point where the semiclassical potential has three equally
deep potential wells. (c) and (d) Energy levels as a function of $\omega t$,
and in black lines the minima and maxima of the semiclassical potential are
also plotted. In (c) and (d) the avoided level crossings of the ground state
are indicated with a dotted vertical line, and the crossing with the Maxwell
set with a dashed vertical line. At the bottom, schematic pictures of the
semiclassical potential for each region are shown. (e) Ground state energy as
a function of $S$ for the parameter values used in (a) for $\omega t=0.4$. The
red horizontal line represents the semiclassical minimum.

###  6.1 Magnetic field trajectories

The magnetic field components $x$ and $z$ are free and we can use them to
create a trajectory in the parameter space, i.e., we can manipulate the system
such that it passes from one phase to another one. We will do it by a time
varying magnetic field.

A quantum phase transition occurs when the parameters $r_{1}$ and $r_{2}$ are
varied from one region to another, crossing a Maxwell set for the global
minimum. We consider the variation of the magnetic field to have the following
form: $\vec{B}=(|B|\cos\omega t+B_{x,0},0,|B|\sin\omega t+B_{z,0})$. In
figures 4(a) and 4(b) this can be seen as the black circles, with the arrow
indicating its direction: For case I we used $\vec{B}=(15.63\cos\omega
t,0,15.63\sin\omega t)$, and for case II we used $\vec{B}=(4.47\cos\omega
t,0,4.47\sin\omega t-32.9)$. The white dots indicate the points where a
Maxwell set for the global minimum is crossed. The effect of crossing the
Maxwell set can be seen in figures 4(c) and 4(d), where the energy levels as a
function of $\omega t$ are plotted. Here we can see how an avoided level
crossing of the ground state (dotted line) occurs near the crossing of the
Maxwell set (dashed line). In black lines we plotted the maxima and minima of
the semiclassical potential as a function of $\omega t$, and at the bottom an
schematic picture of the semiclassical potential for each region is shown.

An additional trajectory is considered in figure 4(a), where $B_{x}$ is fixed
at zero or almost zero, and the component $B_{z}$ is varied passing through
the point $4$, which correspond to the value of $B_{z}$ for which the
semiclassical potential of the system has three equally deep stable wells.

In figure 4(e) we have plotted the value of the ground state for one value of
the magnetic field at $\omega t=0.4$ for case I as a function of $S$. As $S$
increases the ground state approaches the value of the global minimum of the
semiclassical potential shown as an horizontal red line.

Figure 5: Squared amplitude $|\langle\psi_{1}|\zeta\rangle|^{2}$ of the ground
state eigenvector plotted on the Bloch sphere for various values of $\omega t$
for the trajectory shown in figure 4(a). Below each sphere the corresponding
semiclassical potential as a function of $\theta$ is shown.

###  6.2 Bloch sphere visualization

As shown in the schematic representation of the semiclassical potential in
each of the phases,as the Maxwell set is crossed, the deepest potential well
changes from one to another. In this subsection we present a way to visualize
the localization of the eigenstates of the Hamiltonian across these quantum
phase transitions. Using the expression of the eigenvectors of the Hamiltonian
in the basis of $|M\rangle$ in (8) and the definition of the coherent states
in (9) we can define the complex function $\langle\psi_{k}|\zeta\rangle$ as
the projection of the eigenvector on the coherent states:

$\langle\psi_{k}|\zeta\rangle=\sum_{M=-S}^{S}c_{k,M}^{*}\left(\frac{(2S)!}{(S+M)!(S-M)!}\right)^{1/2}\left(\cos\frac{\theta}{2}\right)^{S-M}\left(e^{-\mathrm{i}\phi}\sin\frac{\theta}{2}\right)^{S+M},$
(46)

where according to the definition of the coherent states $(\theta,\phi)$ are
the angular coordinates on the Bloch sphere. The latitudes of the Bloch sphere
are interpreted as the different spin projections $M$, with the south pole
corresponding to $M=S$ and the north pole to $M=-S$, and for the case when $S$
is an integer we have that $M=0$ corresponds to the equator. The projection of
the coherent states onto the eigenstates has been treated in [46].

To show the information provided by (46) we will apply it to the trajectory
shown in figure 4(a). In figure 5 we show $|\langle\psi_{1}|\zeta\rangle|^{2}$
for the ground state for various values of $\omega t$, to depict points
before, during and after a transition. Note, that for some values of $\omega
t$ the coherent state is localized to a definite area on the Bloch sphere.
Thus, controlling $\omega t$, stopping when a certain value is reached, we can
shift the coherent state at a certain area. In other words, we have a switch.

To complement this vision, in figure 5, we also show the semiclassical
potential as a function of $(\theta,\phi_{c})$. At the start, $\omega t=0.3$,
the ground state is localized at the equator, with $\theta\approx\pi/2$ and
$\phi=0$, or in cartesian coordinates about $(x=1,y=0,z=0)$. As the magnetic
field moves to $\omega t=1.1$ we can see the localization of the ground state
starts to shift to the south pole. We can interpret this as a competition
between the two wells close in energy, until, at $\omega t=1.6$ the well at
about $\theta\approx\pi$ is the deepest. In turn, in figure 5 for that value,
now the eigenstate of the ground state is localized in the south pole. As
$\omega t$ increases the ground state continues to travel along the Bloch
sphere: It started in the equator at about $(x=1,y=0,z=0)$, then it goes to
the south pole, then to the equator at about $(x=-1,y=0,z=0)$, then to the
north pole, and finally back to its starting point.

Figure 6: (a), (b) and (c) Energy levels, for $S=20$, as a function of $r_{2}$
for the trajectory that crosses point $4$, which corresponds to the point
where the potential has three equally deep wells, shown in figure 1(b). The
vertical dashed line corresponds to the semiclassical triple point at
$r_{2}=-44.3$. The vertical dotted line corresponds to the avoided level
crossing of the ground state at $r_{2}\approx-44.80887$. (d), (e) and (f):
Squared amplitude $|\langle\psi_{k}|\zeta\rangle|^{2}$ plotted on the Bloch
sphere, viewed from the top of the north pole, for the first three states at
the avoided level crossing of the ground state.

Another interesting feature of these example cases is the presence of the
triple point $4$ in figure 4(b), where the semiclassical potential has three
equally deep wells. In figure 6 we show the structure of the eigenstates as a
trajectory in parameters passes through that point. In the top row we plotted
the energy levels as a function of $r_{2}$, and show a zoom view of the
avoided level crossing at the triple point. At the bottom of the middle plot
schematic pictures of the potential at the different regions are shown. The
crossing of the semiclassical triple point is depicted as a dashed vertical
line at $r_{2}=-44.3$, while the avoided level crossing of the ground state is
depicted as a dotted vertical line at $r_{2}\approx-44.80887$. The minima and
maxima of the semiclassical potential is plotted as a function of $r_{2}$ and
shown as black lines. In the bottom row the
$|\langle\psi_{k}|\zeta\rangle|^{2}$ functions are shown on the Bloch sphere,
viewed from the top of the north pole, for the first three states. Here we can
see how at the point of the avoided level crossing of the ground state, the
ground state distribution has components at the north pole, and at opposite
sides of the equator at $(x=1,y=0,z=0)$ and $(x=-1,y=0,z=0)$, consequence of
the three wells. The second excited states also shares this behaviour, while
the distribution of the first excited state is localized at both sides of the
equator at $(x=1,y=0,z=0)$ and $(x=-1,y=0,z=0)$.

##  7 Conclusions

In the present work we have discussed and studied many properties of spin
Hamiltonians with terms up to fourth order and found how, by the variation of
the external magnetic field, the system can go through the different phases in
parameter space. We applied the powerful method of catastrophe theory.

To provide a real physical system, we studied two examples of single molecule
magnets and showed that the sudden changes in fidelity and fidelity
susceptibility are related with the inverse peaks in the relaxation rate and
the peaks in the relaxation rate. For each of these examples the parameter
space was constructed and it can be seen that the separatrices determine the
structure of the eigenvalues of the Hamiltonian as a function of the magnetic
field. We obtained a hysteresis for two magnetic molecules and demonstrated
that the model applied is able to describe observation.

When considering large values of the fourth order parameters the separatrices
in parameter space obtain a more complicated structure, which permits the
inclusion of more phases, and in particular there exists regions where the
semiclassical potential has three stability wells.

Next, we studied the manner on how to manipulate the position of the coherent
state on the Bloch sphere, using a varying magnetic field. An oscillatory
external magnetic field allows the system to travel through many of these
different phases and the localization of the eigenstates can be visualized on
the Bloch sphere, as well as how it changes along the sphere with the varying
magnetic field, by performing a projection of the coherent states onto the
eigenstates. This result can be viewed as the manipulation of a qubit state by
means of a varying magnetic field, choosing the value of $\omega t$ in order
to reach the desired orientation on the Bloch sphere. In one of the examples
we found a point in parameter space which corresponds to the presence of three
equally deep semiclassical potential wells. Using the visualization of the
distribution of the ground state on the Bloch sphere, we found that at the
avoided level crossing point of the ground state with the first excited state,
three separated regions on the sphere (north pole and $(\pm 1,0,0)$) have non
zero contribution which is a consequence of the three competing stability
points.

We acknowledge the financial support from PAPIIT-DGAPA IN117923, and IN116824.

## References

## References

  * [1] Gatteschi D, Sessoli R and Villain J 2006 Molecular Nanomagnets (Oxford: Oxford University Press)
  * [2] Jiang S, Goss K, Cervetti C and Bogani L 2012 Sci. China Chem. 55 867–882
  * [3] Jenkins M D, Duan Y, Diosdado B, García-Ripoll J J, Gaita-Ariño A, Giménez-Saiz C, Alonso P J, Coronado E and Luis F 2017 Phys. Rev. B 95 064423
  * [4] Moreno-Pineda E, Godfrin C, Balestro F, Wernsdorfer W and Ruben M 2018 Chem. Soc. Rev. 47 501–513
  * [5] Gaita-Ariño A, Luis F, Hill S and Coronado E 2019 Nat. Chem. 11 301–309
  * [6] Moreno-Pineda E and Wernsdorfer W 2021 Nat. Rev. Phys. 3 645–659
  * [7] Zhang Z, Wang Y, Wang H, Liu H and Dong L 2021 Nanoscale Res. Lett. 16 77
  * [8] Cejnar P, Stránský P, Macek M and Kloc M 2021 J. Phys. A: Math. Theor. 54 133001
  * [9] Hagen W R 2024 Molecules 29 2254
  * [10] Gatteschi D and Sessoli R 2003 Angew. Chem. Int. Ed. 42 268–297
  * [11] Cornia A, Gatteschi D and Sessoli R 2001 Coord. Chem. Rev. 219-221 573–604
  * [12] Arecchi F T, Courtens E, Gilmore R and Thomas H 1972 Phys. Rev. A 6(6) 2211–2237
  * [13] Thom R 1975 Structural Stability and Morphogenesis (Reading: W. A. Benjamin)
  * [14] Gilmore R 1981 Catastrophe Theory for Scientists and Engineers (New York: Wiley)
  * [15] Arnold V I 1986 Catastrophe Theory (Berlin: Springer-Verlag)
  * [16] López-Moreno E and Castaños O 1996 Phys. Rev. C 54(5) 2374–2384
  * [17] Lohr-Robles D S, López-Moreno E and Hess P O 2021 Nucl. Phys. A 1016 122335
  * [18] Lohr-Robles D S, López-Moreno E and Hess P O 2023 J. Phys. A: Theor. Math. 56 505301
  * [19] Mannini M, Pineider F, Danieli C, Totti F, Sorace L, Sainctavit Ph, Arrio M A, Otero E, Joly L, Cezar J C, Cornia A and Sessoli R 2010 Nature 468 417–421
  * [20] Fernández J F, Bartolomé J and Luis F 1998 J. Appl. Phys. 83 6940–6942
  * [21] Luis F, Bartolomé J and Fernández J F 1998 Phys. Rev. B 57 505
  * [22] Serrano G, et al. 2020 Nat. Mater. 19 546–551
  * [23] Fichthorn K A and Weinberg W H 1991 J. Chem. Phys. 95 1090–1096
  * [24] Jansen A P J 1995 Comput. Phys. Commun. 86 1–12
  * [25] Liu G B and Liu B G 2009 Appl. Phys. Lett. 95 183110
  * [26] Liu G B and Liu B G 2010 Phys. Rev. B 82 134410
  * [27] Cuppen H M, Karssemeijer L J and Lamberts T 2013 Chem. Rev.  113 8840–8871
  * [28] Thomas L, Lionti F, Ballou R, Gatteschi D, Sessoli R and Barbara B 1996 Phys. Rev. B 66 073309
  * [29] Wernsdorfer W, Caneschi A, Sessoli R, Gatteschi D, Cornia A, Villar V and Paulsen C 2000 Phys. Rev. Lett. 84 2965
  * [30] Ueda M, Maegawa S and Kitagawa S 2002 Phys. Rev. B 66 073309
  * [31] Dressel M, Gorshunov B, Rajagopal K, Vongtragool S and Mukhin A A 2002 Phys. Rev. B 67 060405(R)
  * [32] Adams S T, da Silva Neto E H, Datta S, Ware J F, Lampropoulos C, Christou G, Myaesoedov Y, Zeldov E, and Friedman J R 2013 Phys. Rev. Lett. 110 087205
  * [33] Villain J, Hartman-Boutron F, Sessoli R and Rettori A 1994 Europhys. Lett. 27 159
  * [34] Leuenberger M N and Loss D 1999 Europhys. Lett. 46 692
  * [35] Leuenberger M N and Loss D 2000 Phys. Rev. B 61 1286
  * [36] Gu S J 2010 Int. J. Mod. Phys. B 24 4371–4458
  * [37] Zanardi P and Paunković N 2006 Phys. Rev. E 74(3) 031123
  * [38] Buonsante P and Vezzani A 2007 Phys. Rev. Lett. 98(11) 110601
  * [39] You W L, Li Y W and Gu S J 2007 Phys. Rev. E 76(2) 022101
  * [40] Tzeng Y C, Hung H H, Chen Y C and Yang M F 2008 Phys. Rev. A 77(6) 062321
  * [41] Tian L J, Zhu C Q, Zhang H B and Qin L G 2011 Chin. Phys. B 20 040302
  * [42] Plötz P, Lubasch M and Wimberger S 2011 Physica A 390 1363–1369
  * [43] Rams M M and Damski B 2011 Phys. Rev. Lett. 106(5) 055701
  * [44] Wernsdorfer W, Sessoli R, Caneschi A, Gatteschi D, Cornia A and Mailly D 2000 J. Appl. Phys. 87 5481–5486
  * [45] Vergnani L, Barra A L, Neugebauer P, Rodriguez-Douton M J, Sessoli R, Sorace L, Wernsdorfer W and Cornia A 2012 Chem. Eur. J. 18 3390–3398
  * [46] López-Moreno E and Grether M 2014 Quantum Stud.: Math. Found. 1 203–211