Chapter 3 General Optimisation Algorithm to Find Classical Magnetic Ground States

In order to make predictions of physical properties such as magnetic torque and linear magnetisation, calculations of magnetic structure using a gradient descent algorithm are performed in this thesis. This allowed development of a classical understanding of the pseudo-spin model of the materials which are discussed later. This approach treats spins as classical vectors and allows calculations to be performed much faster than quantum spin models such as exact-diagonalisation, which is intractable for many three-dimensional systems. A classical, three-dimensional pseudospin model is particularly useful when the system orders, and classical models are good approximations [81, 59].

The Hamiltonian for a collection of interacting spins in some applied magnetic field $𝑩$ can be written as

	$\mathscr{H}$	$={\displaystyle \sum _{i\alpha \beta ij}}{J}_{ij}^{\alpha \beta }{S}_i^{\alpha }{S}_j^{\beta }$	$-$	${\displaystyle \sum _{i\alpha \beta }}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_i^{\beta }$		(3.1)
		$={\displaystyle \sum _{ij}}{𝑺}_i^T{𝕁}_{ij}{𝑺}_j$	$-$	${\displaystyle \sum _i}{𝑩}^T{𝔾}_i{𝑺}_i$		(3.2)

where ${𝕁}_{ij}$ is the exchange tensor between sites $i$ and $j$ and ${𝑺}_i^T$ is the transpose of fixed-length vector spin at site $i$. ${𝔾}_i$ with elements $g_i^{\alpha \beta }$ is the g-factor at site $i$ in units of ${\mu }_B$. The classical ground state is the set of spins ${\psi }_0=\{{𝑺}_i\}$ that minimises $\mathscr{H}$, under the condition that all spins have a fixed length $s_i$:

$$|{𝑺}_i|=s_i\forall i$$

(3.3)

Finding this ground state is fundamentally an optimisation problem in which the energy of the system is minimised over the state space. Numerical solutions to these problems are well studied, in physics as well as computer science. A number of algorithms exist to solve these, such as branch and bound, simulated annealing, and gradient descent.

An approach to finding the exact solution in zero field is the Luttinger-Tisza method [63], which works by replacing the strong condition (3.3) with the so-called weak condition

$$\sum _i{|{𝑺}_i|}^2-s_i^2=0$$

(3.4)

which is necessary but insufficient for the validity of (3.3). The Hamiltonian can be exactly solved under this condition using a Lagrange multiplier method [64]. If the solution also satisfies the strong condition then it is the true ground state.

A commonly used numerical method for solving these problems in physics is simulated annealing [54]. This iterative solution involves considering the current state ${𝒙}_n$ with energy $E$ and some ‘neighbouring’ state ${𝒙}_n^{\ast }$ with energy $E^{\ast }$. The state is updated to this new value with probability $P(E,E^{\ast },T)$ which is commonly

(3.5)

where $T$ is a ‘temperature’ parameter which decreases over time. This algorithm provably converges if the cooling schedule is sufficiently slow, however it is challenging to determine the ideal cooling rate a priori.

A different method, used in this thesis, is based on gradient descent (GD). GD works by adjusting the state in the direction which locally minimises the energy. This method reliably and efficiently finds local minima in energy, but some care is needed to find global minima. A more in-depth explanation of the algorithm is in a later section.

Numerical mean field calculations of magnetic structures are particularly useful when considering systems with many exchange terms, or with applied magnetic fields along non-symmetry directions. In order to demonstrate their utility, I will use the example of a square lattice antiferromagnet with nearest-neighbour anisotropic XXZ exchange. A schematic diagram of such a lattice can be seen in Figure 3.1, which shows the primitive lattice vectors $𝒂$ and $𝒃$.

In this system, the total energy can be written as

$$\mathscr{H}=\frac{1}{2}\sum _{i,j}J\left(S_i^x{S}_j^x+S_i^y{S}_j^y+\mathrm{\Delta }{S}_i^z{S}_j^z\right)-\sum _{i}g{\mu }_B𝑩\cdot {𝑺}_i$$

(3.6)

where the half prefactor is to correct for double counting of bonds, $J>0$ parametrises the antiferromagnetic exchange, and $\mathrm{\Delta }$ allows for exchange anisotropy. In this example $\mathrm{\Delta }>1$ indicates easy-axis Ising-type anisotropy.

The energy on a bond can be minimized by having the adjacent spins be antiparallel and collinear with the Ising $𝒛$ axis, such that the bond energy is $-\mathrm{\Delta }J{S}^2$. On the square lattice, this can be satisfied for every bond by adopting a Néel order, such that adjacent spins are all always antiparallel, shown as alternating $+$ and $-$ symbols on Figure 3.1. Because this structure minimises the energy on every bond, it must be the ground state when $𝑩=0$.

The magnetic unit cell vectors are

(3.7)

with

(3.8)

showing that there are two primitive lattice points per magnetic unit cell. As such, the magnetic unit cell volume (in this case, area) is twice as large as the primitive unit cell volume. These magnetic unit cell vectors can be seen in Figure 3.1b, the outline of the unit cell can be seen in Figure 3.1c, with the two magnetic sublattices labelled $A$ and $B$.

In order to understand how the spin orientations vary with applied magnetic field in this system, it is assumed that the magnetic unit cell does not change with applied magnetic field. The mean field energy per magnetic unit cell (Eqn. 3.6) can then be rewitten in terms of the spin orientations of the two sublattices ${𝑺}_A$ and ${𝑺}_B$

$$\frac{\mathscr{H}}{N}=nJ\left(S_A^x{S}_B^x+S_A^y{S}_B^y+\mathrm{\Delta }{S}_A^z{S}_B^z\right)-g{\mu }_B𝑩\cdot \left({𝑺}_A+{𝑺}_B\right)$$

(3.9)

where $n=4$ is the number of nearest-neighbours on each site, and $N$ is the total number of magnetic unit cells in the crystal. In this simple case, constructing the Hamiltonian as a function of sublattice spins is straightforward by inspection. In the most general case with longer range interactions and arbitrary numbers of sublattices the formalism is non-trivial and this is discussed in detail in Section 3.3.

One interesting feature of this system is that it undergoes a spin-flop transition when $\mathrm{\Delta }>1$, this is a metamagnetic phase transition which occurs when a sufficiently strong field is applied along the Ising $𝒛$ axis. The zero-field ground state has spins along the $𝒛$ axis, so when a field is applied the Zeeman energy

$$E_{\text{Z}}=-\frac{1}{2}g{\mu }_B𝑩\cdot \left({𝑺}_A+{𝑺}_B\right)$$

(3.10)

does not decrease with a magnetic field along $𝒛$. The spins are also not able to cant towards the field as they would if the field was applied in the $xy$ plane. Determining the field at which this occurs, and the field dependence of the magnetisation can be done analytically and is demonstrated in many textbooks on the subject [13]. In this model, the spin-flop field $B_{\text{SF}}$ is found to be

$${\mu }_{B}g{B}_{\text{SF}}=SnJ\sqrt{{\mathrm{\Delta }}^2-1}$$

(3.11)

where a spin-flop occurs for all values of $\mathrm{\Delta }>1$. At this field, the magnetisation changes discontinuously from $0$ at to a finite value at $B>B_{\text{SF}}$.

However, the question of how the system behaves for non-axial fields, some arbitrary angle away from the $𝒛$ axis, cannot be answered so clearly analytically. For this reason, it is useful to develop numerical methods for solving the problem. A natural approach is an iterative algorithm called gradient descent. This algorithm works by starting from a state and finding the direction in which the spins can be adjusted that decreases the energy the most. This is repeated until a minimum in energy is found with satisfactory precision.

In this example, Equation 3.9 can be differentiated to give

$$\frac{1}{N}{\nabla }_{{𝑺}_A}\mathscr{H}=nJ\left(\begin{array}{c}{S}_B^x\\ S_B^y\\ \mathrm{\Delta }{S}_B^z\end{array}\right)-g{\mu }_B𝑩$$

(3.12)

which can be computed, then the sublattice spin ${𝑺}_A$ adjusted along the $-{\nabla }_{{𝑺}_A}\mathscr{H}$ direction, such that the energy decrease is maximised. Unfortunately, this runs into the problem that the magnitude of the spins $\left|{𝑺}_A\right|$ are unconstrained. For instance, in the Néel order ground state with no field

	${𝑺}_A$	$=S\widehat{z}$
	${𝑺}_B$	$=-S\widehat{z}$
	${\displaystyle \frac{1}{N}}{\nabla }_{{𝑺}_A}\mathscr{H}$	$=-\mathrm{\Delta }nJS\widehat{z}$

so adjusting ${𝑺}_A$ by $-\gamma {\nabla }_{{𝑺}_A}\mathscr{H}$ with small $\gamma >0$ will increase the magnitude of $\left|{𝑺}_A\right|$

$$\left|{𝑺}_A-\gamma {\nabla }_{{𝑺}_A}\mathscr{H}\right|=\left|{𝑺}_A\right|-\gamma {𝑺}_A\cdot {\nabla }_{{𝑺}_A}\mathscr{H}+𝒪({\gamma }^2)$$

(3.13)

To solve this problem, I note that in the limit of small step size $\gamma \ll 1$, $\left|{𝑺}_A\right|$ will not change if ${𝑺}_A\cdot {\nabla }_{{𝑺}_A}\mathscr{H}=0$. This means that the sublattice spins will be rotated in the plane containing the spin vector and the gradient ${\nabla }_{{𝑺}_A}\mathscr{H}$. I approach the formalism for this in a later section, but in summary this treats the condition $\left|{𝑺}_A\right|=S$ as a spherical manifold $ℳ$, and the adjustments to the state occur in the tangent space of this manifold $T_{𝒑}ℳ$.

In order to keep the step direction in this plane, I project the energy minimising vector $-{\nabla }_{{𝑺}_A}\mathscr{H}$ into the plane and use that as the step direction, and then normalise the new vector again to avoid accumulation of $𝒪({\gamma }^2)$ errors.

$${𝑺}_{A,n+1}={\mathrm{proj}}_S\left[{𝑺}_{A,n}-\gamma {\mathrm{proj}}_T{\nabla }_{{𝑺}_A}\mathscr{H}\right]$$

(3.14)

where ${\mathrm{proj}}_S$ indicates the projection onto a sphere with radius $S$ and ${\mathrm{proj}}_T$ indicates projection into the plane normal to ${𝑺}_A$. This can then be iterated for each spin until an energy minimum is found. Section 3.4 formalises this approach and derives a generalisation of the above expression.

This calculation was performed with the parameters in Table 3.1. From these, one expects a spin-flop transition to occur at ${\mu }_B{B}_{\text{SF}}=1\mathrm{meV}$ or $B_{\text{SF}}\approx 17.28\mathrm{T}$ and polarisation to occur at ${\mu }_B{B}_{\text{P}}=2.414\mathrm{meV}$ or $B_{\text{P}}\approx 41.71\mathrm{T}$. The results of these calculations can be seen in Figure 3.3. Fig. 3.3a) shows the evolution of magnetisation against magnetic field for field orientations parallel and perpendicular to the $𝒛$ axis. For fields perpendicular to $𝒛$, the magnetisation increases linearly in field from $0{\mu }_B$ at $0\mathrm{T}$ to $1{\mu }_B$ at $41.7(1)\mathrm{T}$. For fields parallel to $𝒛$, the magnetisation remains at $0{\mu }_B$ between $0\to 17.3(1)\mathrm{T}$, and then discontinuously joins the $𝑩⟂𝒛$ curve at higher fields.

Figure 3.2 shows how the projection of magnetisation onto the magnetic field direction

$$𝑴\cdot \frac{𝑩}{|𝑩|}$$

(3.15)

varies as as the magnetic field vector $𝑩$ varies in the $xz$ plane. This shows a spin-flop transition at the isolated point $𝑩=B_{\text{SF}}𝒛$, such that if the magnetic field is applied directly along the $𝒛$ axis, there is a phase transition, and discontinuous change in magnetisation at $|𝑩|=B_{\text{SF}}$. This also shows that for off-axis fields, there is no phase transition, but instead a smooth crossover between the two regimes. It show that the spin orientations which minimise the free energy take different angles with the magnetic field $\alpha \ne \beta$ except in the special cases of axial fields $\theta =0$, $\theta =\frac{\pi }{2}$, highlighting why this calculation is challenging to do analytically. For all field orientations, the system polarises at the same field $|𝑩|=B_{\text{P}}\approx 41.71\mathrm{T}$, which is in all cases a continuous phase transition.

Fig. 3.3b-d) show the effect of off-axis fields on this system. Fig. 3.3b) shows how the angle $\theta$ is defined in this setup: measuring an angle from the $𝒛$ axis, towards the $xy$ plane. At zero-field the spins are aligned parallel to the $𝒛$ axis. Fig. 3.3c) shows how magnetisation against field curves for $\theta$ between $0\mathrm{°}$ and $15\mathrm{°}$. At $\theta =0\mathrm{°}$ the field is exactly along axis and the $𝑩\parallel 𝒛$ curve discussed above is recovered. As the field gets further off axis, the discontinuity at $|𝑩|=B_{\text{SF}}$ is immediately replaced by a smooth crossover between the low and high-field regimes. At higher values of $\theta$, the crossover becomes less sharp. Fig. 3.3d) shows how magnetic torque, taken around the $𝒚=𝒛\times 𝒙$ axis shown in Fig. 3.3b), evolves with applied magnetic field at different values of $\theta$. At $\theta =0\mathrm{°}$ the field is exactly along the $𝒛$ axis, and the torque is constrained by symmetry to be exactly $0$ at all fields. For small values of $\theta$, there is a sharp peak in torque around the spin-flop field. As $\theta$ increases, the peak becomes broader, and increases slightly in field. At the polarisation field $B_{\text{P}}$, the magnetisation changes slope at all angles $\theta$, indicating a continuous phase transition. There is no anomaly in torque at $B_{\text{P}}$ at any angle $\theta$. These results reproduce what is expected theoretically [84, 14] has been seen before experimentally [113].

The previous section outlines the ground state determination method for the case of two sublattices in the magnetic unit cell and nearest-neighbour exchange. This is appropriate for many simple systems. This section demonstrates the appropriate formalism for multi-sublattice systems with exchanges beyond nearest-neighbours.

The general form of the a spin Hamiltonian with exchange and Zeeman terms is

	$\mathscr{H}$	$={\displaystyle \sum _{\alpha \beta ij}}{J}_{ij}^{\alpha \beta }{S}_i^{\alpha }{S}_j^{\beta }$	$-$	${\displaystyle \sum _{i\alpha \beta }}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_i^{\beta }$		(3.16)
		$={\displaystyle \sum _{ij}}{𝑺}_i^T{𝕁}_{ij}{𝑺}_j$	$-$	${\displaystyle \sum _i}{𝑩}^T{𝔾}_i{𝑺}_i$		(3.17)

allowing for all bilinear exchange terms $J_{ij}^{\alpha \beta }$ and anisotropic g-factors $g_i^{\alpha \beta }$. This form of the Hamiltonian requires indexing over all sites in the crystal. These sites however exist on a Bravis lattice with symmetry constraints which are hard to evaluate in this indexed form. An equivalent formalism for the Hamiltonian is to make this explicit, by summing over the sites in each primitive unit cell, and the lattice points in the crystal. This gives

$$\mathscr{H}=\sum _{\alpha \beta ij𝒖𝒗}{J}_{ij𝒖𝒗}^{\alpha \beta }{S}_{i𝒖}^{\alpha }{S}_{j𝒗}^{\beta }-\sum _{\alpha \beta i𝒖}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{𝒖}^{\beta }$$

(3.18)

where $i,j$ index over the sites in the unit cell and $𝒖,𝒗$ index over the lattice vectors in the crystal. Additionally, the value of the Hamiltonian must be invariant under the symmetry operations of the crystal, including translations by lattice vectors. In other words, it must be invariant under the transformation:

$${𝑺}_{i𝒖}\to {𝑺}_{i𝒖+𝒕}$$

(3.19)

where $𝒕$ is a lattice vector. This constrains the exchange parameters, allowing them to be rewritten:

$$J_{ij𝒖,𝒖+𝒕}\to J_{ij𝒕}$$

(3.20)

So including only the unique terms allowed by translational symmetry of the lattice:

$$\mathscr{H}=\sum _{\alpha \beta ij𝒖𝒕}{J}_{ij𝒕}^{\alpha \beta }{S}_{i𝒖}^{\alpha }{S}_{j𝒖+𝒕}^{\beta }-\sum _{\alpha \beta i𝒖}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i𝒖}^{\beta }$$

(3.21)

The exchange part of the Hamiltonian can then be fully specified by a set of exchange parameters (nearest-neighbour, next-nearest neighbour, and so on).

In order to convert the Hamiltonian into a sum over sublattices, a magnetic unit cell needs to be assumed. This means that the spins $S_{i𝒖}$ and $S_{i𝒗}$ at the same site index $i$, for two different lattice vectors $𝒖$, $𝒗$, are equivalent if $𝒖$ and $𝒗$ differ by a magnetic lattice vector. Formally:

$${𝑺}_{i𝒖}\equiv {𝑺}_{i𝒗}\text{if }\exists 𝒎\in 𝕄:𝒖=𝒗+𝒎$$

(3.22)

where $𝕄$ is the vector space of all magnetic lattice vectors.

This means that the lattice vectors $𝒖$, $𝒗$ can be treated as equivalent if they differ by a magnetic lattice vector. This is the same as performing the algebra in the quotient space $V/𝕄$, where $V$ is the vector space of primitive lattice vectors. It is worth noting that $𝕄$ is a vector subspace of $V$. A useful operator to define is the equivalence class of a vector in $V$

$$\left[𝒖\right]≔\{𝒖+𝒎\mid \forall 𝒎\in 𝕄\}$$

(3.23)

which is the set of all primitive lattice vectors which are equivalent under the magnetic unit cell.

To put it simply, by imposing a magnetic unit cell $𝕄$, two spins ${𝑺}_{i𝒖}$ and ${𝑺}_{i𝒗}$ at the same site index $i$, for two different primitive lattice vectors $𝒖$, $𝒗$, are equivalent if $\left[𝒖\right]=\left[𝒗\right]$. The set of equivalence classes here can be understood as the set of primitive lattice points in each magnetic unit cell. If there is only one site per primitive cell, then the equivalence classes are just the sublattices. If there is more than one site per primitive cell, the sublattices can be indexed by $i,\left[𝒖\right]$, and the number of sublattices is equal to the number of sites per primitive unit cell $\times$ the number of primitive lattice points per magnetic cell.

Because of this imposed symmetry, the spin at site $i,𝒖$ can be rewritten

$${𝑺}_{i𝒖}\to {𝑺}_{i\left[𝒖\right]}$$

(3.24)

and the Hamiltonian can be rewritten as

$$\mathscr{H}=\sum _{\alpha \beta ij𝒖𝒕}{J}_{ij𝒕}^{\alpha \beta }{S}_{i\left[𝒖\right]}^{\alpha }{S}_{j\left[𝒖+𝒕\right]}^{\beta }-\sum _{\alpha \beta i𝒖}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i\left[𝒖\right]}^{\beta }$$

(3.25)

This still sums over all of the sites in the crystal, but the sum can be simplified by noting that for a finite size magnetic unit cell, there are a finite number of equivalence classes. The sum is taken over the magnetic unit cells, and the sites per magnetic unit cell

	$\mathscr{H}$	$={\displaystyle \sum _{𝒎\in 𝕄}}{\displaystyle \sum _{\alpha \beta ij𝒑𝒕}}{J}_{ij𝒕}^{\alpha \beta }{S}_{i\left[𝒑+𝒎\right]}^{\alpha }{S}_{j\left[𝒑+𝒎+𝒕\right]}^{\beta }-{\displaystyle \sum _{𝒎\in 𝕄}}{\displaystyle \sum _{\alpha \beta i𝒑}}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i\left[𝒑+𝒎\right]}^{\beta }$		(3.26)
		$=N\left({\displaystyle \sum _{\alpha \beta ij𝒑𝒕}}{J}_{ij𝒕}^{\alpha \beta }{S}_{i\left[𝒑\right]}^{\alpha }{S}_{j\left[𝒑+𝒕\right]}^{\beta }-{\displaystyle \sum _{\alpha \beta i𝒑}}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i\left[𝒑\right]}^{\beta }\right)$		(3.27)

where the first sum over is taken over all magnetic unit cells in the crystal $𝒎\in 𝕄$, $N$ is the number of magnetic unit cells in the crystal and the sum over $𝒑$ is over the primitive lattice points per magnetic unit cell (the equivalence classes of $V/𝕄$). The above uses the property of the equivalence classes that $\left[𝒑+𝒎\right]=\left[𝒑\right]\forall 𝒎\in 𝕄$.

Here, the factor of $N$ can be safely removed by transforming $N\mathscr{H}\to \mathscr{H}$, such that

$$\mathscr{H}=\sum _{\alpha \beta ij𝒑𝒕}{J}_{ij𝒕}^{\alpha \beta }{S}_{i\left[𝒑\right]}^{\alpha }{S}_{j\left[𝒑+𝒕\right]}^{\beta }-\sum _{\alpha \beta i𝒑}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i\left[𝒑\right]}^{\beta }$$

(3.28)

The final hurdle to solve is the sum over $𝒕$, which still sums over all primitive lattice points. This is easily resolved as in the calculation, only a small number of interaction terms $J_{ij𝒕}^{\alpha \beta }$ are non-zero. This allows the sum to run over all non-zero interactions which have been defined (nearest-neighbour, next-nearest-neighbour, etc.), which is likely to be small in any real-world application, except for dipolar coupled systems which are not the focus of this thesis.

For the sake of completeness, the classical Hamiltonian can be rewritten in terms of the exchange interaction between the sublattices by defining

$${\tilde{J}}_{i\left[𝒑\right],j\left[𝒓\right]}^{\alpha \beta }≔\sum _{𝒎\in 𝕄}{J}_{ij\left(𝒓-𝒑+𝒎\right)}^{\alpha \beta }$$

(3.29)

such that

$$\mathscr{H}=\sum _{\alpha \beta i[𝒑]j[𝒓]}{\tilde{J}}_{i\left[𝒑\right],j\left[𝒓\right]}^{\alpha \beta }{S}_{i\left[𝒑\right]}^{\alpha }{S}_{j\left[𝒓\right]}^{\beta }-\sum _{\alpha \beta i[𝒑]}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_{i\left[𝒑\right]}^{\beta }$$

(3.30)

where $i[𝒑],j[𝒓]$ sum pairwise over the sublattices. This is the same form as Equation 3.9 in the example, and provides a good starting point for numerically solving these problems. The following section formalises a method for numerically minimising the mean field energy for Hamiltonians of this form.

The effect of applying a finite size magnetic unit cell is that only a subset of magnetic propagation vectors $𝒒$ are allowed, such that

$$𝒒\cdot 𝒎\in \mathbb{Z}\forall 𝒎\in 𝕄$$

(3.31)

This ‘rationalisation’ of k-space is the standard consequence of applying periodic boundary conditions to a crystal. It should be noted that incommensurate magnetic structures have irrational propagation vectors, so can never be solved with this method.

As the choice of magnetic unit cell $𝕄$ limits the solutions that can be found, it must be chosen from either prior theoretical work or empirical evidence of the propagation vector such as neutron diffraction experiments.

This section formalises a method for solving optimisation problems where the state space is not a vector space. This is needed for finding the ground state magnetic Hamiltonians of the form

	$\mathscr{H}$	$={\displaystyle \sum _{\alpha \beta ij}}{\tilde{J}}_{ij}^{\alpha \beta }{S}_i^{\alpha }{S}_j^{\beta }$	$-$	${\displaystyle \sum _{\alpha \beta i}}{g}_i^{\alpha \beta }{B}^{\alpha }{S}_i^{\beta }$	$+$	${\displaystyle \sum _i}{A}_i({𝑺}_i)$		(3.32)
		$={\displaystyle \sum _{ij}}{𝑺}_i^T{𝕁}_{ij}{𝑺}_j$	$-$	${\displaystyle \sum _i}{𝑩}^T{𝔾}_i{𝑺}_i$	$+$	${\displaystyle \sum _i}{A}_i({𝑺}_i)$		(3.33)

where $A_i({𝑺}_i)$ parametrises local anisotropy on site $i$, some complex anisotropy terms are discussed in Chapter 6.

These have the constraint that the sublattice spin vectors have a fixed magnitude i.e.

$$\left|{𝑺}_i\right|=s_i\forall i$$

(3.34)

As such, the space of allowed ${𝑺}_i$ do not form a vector space. Fortunately, as each spin is constrained to be a vector with a fixed length, it is equivalent to saying that each ${𝑺}_i$ is a point on a sphere of radius $s_i$. This forms a differentiable Riemannian manifold over which some optimisation algorithms such as gradient descent can be generalised as shown below.

Traditionally, gradient descent is an algorithm for finding the local minimum of a scalar, differentiable, function $F$ over a vector space

$$F:{ℝ}^d\to ℝ$$

(3.35)

The principle of the algorithm comes from defining the gradient operator $\nabla F$ such that its inner product with some vector direction $𝒗$

(3.36)

the derivative of $F$ along $𝒗$. This is equivalent to the more commonly used definition

$$\nabla F≔\sum _i\frac{\partial F}{\partial x_i}{\widehat{𝒆}}_i$$

(3.37)

where $\widehat{{𝒆}_i}$ is the $i$^thunit basis vector of the vector space, and $\frac{\partial F}{\partial x_i}$ is the partial derivative of $F$ with respect to the coordinate $x_i$.

For a small step $𝜹$, the change in $F$ is maximised when the step is parallel to the gradient $\nabla F$. This prompts the iterative solution

$${𝒙}_{n+1}={𝒙}_n-\gamma \nabla F({𝒙}_n)$$

(3.38)

where $\gamma$ is a learning rate, and is a free parameter.

In order to include the constant length spins constraint (Eqn. 3.3), modifications have to be made to the standard gradient descent algorithm. One option is to use a penalty method, which involves adding some smooth term to the Hamiltonian to penalise moving away from the constraint. In this case one could add the penalty term:

$$P=\sum _i\eta {(|{𝑺}_i|-1)}^2$$

(3.39)

where $\eta >0$ is a free parameter. This would add a cost to the optimiser to violate the constraint, so it will stay approximately satisfied.

Another solution, with has been implemented for this project, is to consider the constrained state space $\mathrm{\Psi }$ as a Riemannian manifold, described below. The gradient descent algorithm can be generalised to apply to any Riemannian manifold, not just the Euclidean geometry ${ℝ}^d$ shown above [2, 88, 87]. In this form the algorithm is sometimes known as Riemannian Optimisation. However, some of the steps which are trivial in the Euclidean case are not trivial in the general case.

A Riemannian manifold is a smooth, differentiable manifold whose tangent spaces have an inner product, which varies smoothly across the manifold. The tangent space $T_pℳ$ of a differentiable manifold $ℳ$at $p$ is the vector space spanned by the tangent vectors of all of the curves at $p$ on $ℳ$. If $ℳ$can be embedded in an ambient Euclidean space, this can be easily visualised as the (hyper-)plane spanned by the tangent vectors as in Figure 3.4a).

In this case, $F$ is instead a differentiable, scalar function over the manifold $ℳ$

$$F:ℳ\to ℝ$$

(3.40)

This allows the gradient $\nabla F$ to be defined again such that

(3.41)

is the derivative of $F$ along $𝒗$. It should be noted that the vectors $𝒗$ and $\nabla F$ are now elements of the tangent space $T_pℳ$. Clearly, a step along the vector $\nabla F$ will maximise the change in $F$, so the approach of gradient descent will work. However, unlike in the Euclidean case, the tangent space and the manifold have very different geometry. A tangent vector on $T_pℳ$ can be mapped to $ℳ$by the exponential map:

$${\exp }_p:T_pℳ\to ℳ$$

(3.42)

In some cases, this map may be computationally expensive to evaluate. Overall this gives the iteration step for gradient descent on an arbitrary Riemannian manifold:

$$p_{n+1}={\exp }_{p_n}(-\gamma \nabla F)$$

(3.43)

Applying this to the geometry of the problem in magnetism is straightforward. A point $\psi$ on the overall space $\mathrm{\Psi }$ is a collection of $N$ spins ${𝑺}_i$, which are constrained to be unit vectors. So each spin is a point on a sphere, that is a point on the ${\text{S}}^2$ spherical manifold. The overall space $\mathrm{\Psi }$ is therefore the product manifold

$$\mathrm{\Psi }=\underset{N\text{times}}{\underbrace{{\text{S}}^2\times {\text{S}}^2\times \mathrm{\dots }\times {\text{S}}^2}}$$

(3.44)

and is easily visualised as $N$ independent coordinates on a sphere. Formally:

	$\psi$	$=(s^1,\mathrm{\dots },s^N)\in \mathrm{\Psi }$		(3.45)
	$s^i$	$\in S^2$		(3.46)

The tangent space of a product manifold $X\times Y$ is the product of the tangent spaces of the individual manifolds [2]:

$$T_{(x,y)}(X\times Y)=T_{x}X\times T_{y}Y$$

(3.47)

The exponential map has a similar relationship [2]:

$${\exp }_{(x,y)}^{X\times Y}=({\exp }_x^X,{\exp }_y^Y)$$

(3.48)

This means that the update step with $N$ spins can be divided into $N$ update steps each on a sphere:

	${\psi }_{n+1}=(\mathrm{\dots },s_{n+1}^i,\mathrm{\dots })$	$={\exp }_{p_n}^{\mathrm{\Psi }}(-\gamma {\nabla }_{\mathrm{\Psi }}F)$		(3.49)
		$=(\mathrm{\dots },{\exp }_{s_n^i}^{S_i}(-\gamma {\nabla }_{S_i}F),\mathrm{\dots })$		(3.50)
	$\therefore s_{n+1}^i$	$={\exp }_{s_n^i}^{S_i}(-\gamma {\nabla }_{S_i}F)$		(3.51)

Thankfully, the spherical geometry allows this to be computed without much trouble. Figure 3.4a,b) show how the exponential map of the $d$-sphere ${\text{S}}^d$ can be cheaply approximated by projecting the first order approximation of the map, the point on the tangent (hyper-)plane, back onto the sphere:

	${\exp }_{𝒑}(𝒖)$	$=𝒑\cos |𝒖|+{\displaystyle \frac{𝒖}{|𝒖|}}\sin |𝒖|$		(3.52)
		$\approx {\mathrm{proj}}_S(𝒑+𝒖)={\displaystyle \frac{𝒑+𝒖}{|𝒑+𝒖|}}={\exp }_{𝒑}\left({\displaystyle \frac{\arctan (|𝒖|)}{|𝒖|}}𝒖\right)$		(3.53)

The effect of the approximation is to slightly lower the learning rate while the path along which the state is adjusted is unchanged.

By embedding the sphere in ${ℝ}^3$, the coordinate of the point on the sphere can be stored as $x$, $y$, $z$ components. The tangent space $T_p{S}^2$ is the plane $𝒗\cdot 𝒑=0$, and is a 2-dimensional vector subspace of the ambient ${ℝ}^3$. This allows for a very simple method for computing the gradient from the unconstrained Hamiltonian $\tilde{\mathscr{H}}:{ℝ}^3\to ℝ$. The unconstrained gradient $\nabla \tilde{\mathscr{H}}\in {ℝ}^3$ can be computed cheaply and projected into the tangent plane to give the gradient of the constrained Hamiltonian:

$${\nabla \mathscr{H}|}_{𝒑}={\mathrm{proj}}_{T_{𝒑}S}\nabla \tilde{\mathscr{H}}=\nabla \tilde{\mathscr{H}}-(\nabla \tilde{\mathscr{H}}\cdot 𝒑)𝒑$$

(3.54)

This gives the final expression for the iteration step in the gradient descent algorithm on the sphere as

	${𝒙}_{n+1}$	$={\mathrm{proj}}_S\left({𝒙}_n-\gamma \nabla \mathscr{H}\right)$		(3.55)
		$={\mathrm{proj}}_S\left({𝒙}_n-\gamma {\mathrm{proj}}_T\nabla \tilde{\mathscr{H}}\right)$		(3.56)
		$={\displaystyle \frac{{𝒙}_n-\gamma (\nabla \tilde{\mathscr{H}}-(\nabla \tilde{\mathscr{H}}\cdot {𝒙}_n){𝒙}_n)}{\left|{𝒙}_n-\gamma (\nabla \tilde{\mathscr{H}}-(\nabla \tilde{\mathscr{H}}\cdot {𝒙}_n){𝒙}_n)\right|}}$		(3.57)

where $S$ represents the sphere, and $T=T_{{𝒙}_n}{S}^2$ is the tangent space of the sphere at ${𝒙}_n$. This is the generalised form of the expression found in the example (Equation 3.14). In the implementation of the algorithm used in this thesis, $\gamma$ is chosen to keep the maximum angle a spin rotates below a threshold, which starts at $0.05\mathrm{rad}$, and is decreased until a specified precision is reached.

Numerically this can be computed with double precision floating point arithmetic very quickly, as it contains only linear operations, with the exception of the square-root in the magnitude. The final ingredient is the gradient of the unconstrained Hamiltonian with respect to an individual spin. Analytically this is

$${\nabla }_{{𝑺}_i}\tilde{\mathscr{H}}=\sum _{j\ne i}{𝕁}_{ij}{𝑺}_j+({𝕁}_{ii}+{𝕁}_{ii}^T){𝑺}_i-{𝔾}_i𝑩+{\nabla }_{{𝑺}_i}{A}_i({𝑺}_i)$$

(3.58)

where the term ${𝕁}_{ii}+{𝕁}_{ii}^T$ captures exchange interactions between two different sites in the same sublattice. For instance, on the square lattice example, a next-nearest-neighbour interactions would couple the $A$ and $B$ sublattices to themselves. The term ${\nabla }_{{𝑺}_i}{A}_i({𝑺}_i)$ captures single-ion anisotropy such as a $-D{S}_z^2$ term.

The time complexity of an algorithm is defined as the leading order scaling of the execution time as the size of the input grows, in the limit of large input size. For finding local minima this algorithm has a time complexity of $𝒪(N^2)$, quadratic in the number of spins and equal to the time complexity of computing the gradient of the Hamiltonian. To find global minima in energy, the initial point ${𝒙}_0$ can be guessed at many points across $ℳ$. However, in the case of magnetism, the number of points needed to cover the manifold scales exponentially in the number of spins in the calculation. With an overall time complexity of $𝒪(N^2{b}^N)$, where the branching factor $b$ depends on the exact form of Hamiltonian. In examples with only one local minima $b=1$ the scaling is very efficient, but examples with many local minima in energy are pathological for this algorithm.

In this case, Riemannian optimisation is a significant improvement over penalty methods in that is ensures that the constraint is satisfied exactly, without sacrificing speed. Compared to simulated annealing, it has the advantage of robustly determining local minima in energy, which is useful for tracking the evolution and stability of domains, and for visualising the evolution of the state when there are degenerate ground states. On the other hand, simulated annealing is able to efficiently find approximate global minima in energy, while gradient descent is exceptionally poor at the task for large numbers of spins and many local minima in energy. In such a case, a hybrid approach could work by using simulated annealing to produce an anstaz to seed the gradient descent algorithm, however in all cases in this thesis this was not required.

This algorithm is similar to variations of gradient descent used in many popular mean-field calculation programs. For example, SPINW’s [99] optmagsteep method rotates spins towards the local mean-field, similar to gradient descent. However, the precise formulation used in this thesis means that it provably converges on a local minimum in energy [2].

I used the principles outlined in this chapter to develop the Python library qafm, able to numerically find local and global minima in the mean-field energy with and without magnetic field applied. This can be performed on arbitrary crystal structures and commensurate magnetic unit cells, and is able to verify if given exchange tensors are symmetry allowed. This is used in later chapters to understand theoretically the response of different systems to magnetic fields at the mean-field level.