Interaction Renormalization

A unique feature of Sunny is its support for building classical models where quantum spin is represented as an $N$-level system, rather than just an expected dipole. This generalization can be important when modeling quantum spin Hamiltonians that include, e.g., a single-ion anisotropy, or a biquadratic coupling between sites. Sunny also supports constraining quantum spin to the space of pure dipoles; in this case, Sunny will automatically perform an interaction strength renormalization that maximizes accuracy.

Local operators

A quantum spin-$s$ state has $N = 2s + 1$ levels. Each local spin operator $\hat{S}^{\{x,y,z\}}$ is faithfully represented as an $N×N$ matrix. These matrices can be accessed using spin_matrices for a given label $s$. For example, the Pauli matrices are associated with $s = 1/2$.

When $s > 1/2$, it is possible to construct multipole moments beyond the spin-dipole. For example,

S = spin_matrices(3/2)
@assert S[3] ≈ diagm([3/2, 1/2, -1/2, -3/2])
@assert S[3]^2 ≈ diagm([9/4, 1/4, 1/4, 9/4])

If the operator -S[3]^2 is passed to set_onsite_coupling!, it would set an easy-axis anisotropy in the $\hat{z}$ direction.

Any Hermitian operator can be expanded in the basis of Stevens operators $\hat{\mathcal{O}}_{k,q}$ up to a constant shift. To see this expansion, use print_stevens_expansion:

print_stevens_expansion((S[1]^2 + S[2]^2)) # Prints -(1/3)𝒪₂₀ + 5/2

Alternatively, the same operator could have been constructed directly from stevens_matrices:

O = stevens_matrices(3/2)
@assert S[1]^2 + S[2]^2 ≈ -O[2, 0]/3 + (5/2)*I

See below for an explicit definition of Stevens operators as polynomials of the spin operators.

Renormalization procedure for :dipole mode

Sunny will typically operate in one of two modes: :SUN or :dipole. The former faithfully represents quantum spin as an SU(N) coherent-state which, for our purposes, is an $N$-component complex vector. In contrast, :dipole mode constrains the coherent-state to the space of pure dipoles. Here, Sunny will automatically renormalize the magnitude of each Stevens operator to achieve maximal consistency with :SUN mode. This procedure was derived in D. Dahlbom et al., [arXiv:2304.03874].

By way of illustration, consider a quantum operator $\hat{\mathcal{H}}_{\mathrm{local}}$ giving a single-ion anisotropy for one site. In Stevens operators,

\[\hat{\mathcal H}_{\mathrm{local}} = \sum_{k, q} A_{k,q} \hat{\mathcal{O}}_{k,q},\]

for some coefficients $A_{k,q}$.

In :SUN mode, Sunny will faithfully represent $\hat{\mathcal H}_{\mathrm{local}}$ as an $N×N$ matrix. In :dipole mode, the expected energy $\langle \hat{\mathcal H}_{\mathrm{local}} \rangle$ must somehow be approximated using the expected dipole data.

One approach is to formally take $s \to \infty$, and this yields the traditional classical limit of a spin system. In this limit spin operators commute and expectation values of polynomials become polynomials of expectation values. For example, $\langle \hat{S}^\alpha \hat{S}^\beta\rangle \to \langle \hat{S}^\alpha \rangle \langle \hat{S}^\beta\rangle$, because any corrections are damped by the factor $s^{-1} \to 0$. The expectation of a Stevens operator $\langle \hat{\mathcal{O}}_{k,q} \rangle$ would then become a classical Stevens function $\mathcal{O}_{k,q}(\langle\hat{\mathbf{S}}\rangle)$, i.e., a polynomial of the same form, but now applied to the expected dipole. Classical Stevens functions are constructed as homogeneous polynomials of order $k$, because lower-order terms would vanish in the limit $s \to \infty$.

In a real magnetic compound, however, the spin magnitude $s$ is not necessarily large. To obtain a better approximation, one should avoid the formal limit $s \to \infty$. Our approach is to start with the full dynamics of SU(N) coherent states, and then constrain it to the space of pure dipole states $|\boldsymbol{\Omega}\rangle$. The latter are defined as any states where the expected dipole 3-vector,

\[\boldsymbol{\Omega} ≡ \langle \boldsymbol{\Omega}| \hat{\mathbf{S}} | \boldsymbol{\Omega}\rangle,\]

has maximal magnitude $|\boldsymbol{\Omega}| = s$ and arbitrary direction.

For a pure dipole state, group theory dictates that expectations of the Stevens operators can be expressed as a renormalization of the classical Stevens functions,

\[\langle \boldsymbol{\Omega}| \hat{\mathcal{O}}_{k,q} | \boldsymbol{\Omega}\rangle = c_k \mathcal{O}_{k,q}(\boldsymbol{\Omega}).\]

At fixed $k$, the two sides must be proportional because they are both spin-$k$ irreducible representations of SO(3). The renormalization factors can be calculated explicitly:

\[\begin{align*} c_1 &= 1 \\ c_2 &= 1-\frac{1}{2}s^{-1} \\ c_3 &= 1-\frac{3}{2}s^{-1}+\frac{1}{2}s^{-2} \\ c_4 &= 1-3s^{-1}+\frac{11}{4}s^{-2}-\frac{3}{4}s^{-3} \\ c_5 &= 1-5s^{-1}+\frac{35}{4}s^{-2}-\frac{25}{4}s^{-3}+\frac{3}{2}s^{-4} \\ c_6 &= 1-\frac{15}{2}s^{-1}+\frac{85}{4}s^{-2}-\frac{225}{8}s^{-3}+\frac{137}{8}s^{-4}-\frac{15}{4}s^{-5} \\ &\vdots \end{align*}\]

Constrained to the space of dipoles, the expected local energy becomes

\[E_{\mathrm{local}}(\boldsymbol{\Omega}) = \langle \boldsymbol{\Omega}| \hat{\mathcal H}_{\mathrm{local}} | \boldsymbol{\Omega}\rangle = \sum_{k, q} c_k A_{k,q} \mathcal{O}_{k,q}(\boldsymbol{\Omega}).\]

It can be shown that SU(N) dynamics reduces to the usual Landau-Lifshitz dynamics of dipoles, but involving $E_{\mathrm{local}}(\boldsymbol{\Omega})$ as the classical Hamiltonian. The renormalization factors $c_k$ can therefore be interpreted as a correction to the traditional large-$s$ classical limit.

Renormalization also applies to the coupling between different sites. In Sunny, couplings will often be expressed as a polynomial of spin operators using set_pair_coupling!, but any such coupling can be decomposed as sum of tensor products of Stevens operators. Without loss of generality, consider a single coupling between two Stevens operators $\hat{\mathcal{H}}_\mathrm{coupling} = \hat{\mathcal{O}}_{k,q} \otimes \hat{\mathcal{O}}_{k',q'}$ along a bond connecting sites $i$ and $j$. Upon constraining to pure dipole states $|\boldsymbol{\Omega}_i\rangle$ and $|\boldsymbol{\Omega}_j\rangle$, the expected energy takes the form $E_\mathrm{coupling} = c_k c_k' \mathcal{O}_{k,q}(\boldsymbol{\Omega}_i) \mathcal{O}_{k',q'}(\boldsymbol{\Omega}_j)$, which now involves a product of renormalized Stevens functions.

Use :dipole_large_s mode to disable renormalization

Although we generally recommend the above renormalization procedure, there are circumstances where it is not desirable. Examples include reproducing a model-system study, or describing a micromagnetic system for which the $s\to\infty$ limit is a good approximation. To simulate dipoles without interaction strength renormalization, construct a System using the mode :dipole_large_s instead of :dipole. Symbolic operators in the large-$s$ limit can be constructed by passing Inf to either spin_matrices or stevens_matrices.

Definition of Stevens operators

The Stevens operators $\hat{\mathcal{O}}_{k,q}$ are defined as polynomials of angular momentum operators $\hat{S}_{\{x,y,z\}}$ in some spin-$s$ representation.

Using

\[\begin{align*} X &= \mathbf{\hat{S}} \cdot \mathbf{\hat{S}} = s (s+1) \\ \hat{S}_\pm &= \hat{S}_x \pm i \hat{S}_y \\ \phi_+ &= \frac{1}{4},\quad \phi_- = \frac{1}{4 i}, \end{align*}\]

the relevant Stevens operators are defined as,

\[\begin{align*} \hat{\mathcal{O}}_{0,0} & =1\\ \\ \hat{\mathcal{O}}_{1,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{1,0} & =\hat{S}_{z}\\ \\ \hat{\mathcal{O}}_{2,\pm2} & =\phi_{\pm}(\hat{S}_{+}^{2}\pm\hat{S}_{-}^{2})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{2,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})\hat{S}_{z}+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{2,0} & =3\hat{S}_{z}^{2}-X\\ \\ \hat{\mathcal{O}}_{3,\pm3} & =\phi_{\pm}(\hat{S}_{+}^{3}\pm\hat{S}_{-}^{3})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{3,\pm2} & =\phi_{\pm}(\hat{S}_{+}^{2}\pm\hat{S}_{-}^{2})\hat{S}_{z}+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{3,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})(5\hat{S}_{z}^{2}-X-1/2)+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{3,0} & =5\hat{S}_{z}^{3}-(3X-1)\hat{S}_{z}\\ \\ \hat{\mathcal{O}}_{4,\pm4} & =\phi_{\pm}(\hat{S}_{+}^{4}\pm\hat{S}_{-}^{4})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{4,\pm3} & =\phi_{\pm}(\hat{S}_{+}^{3}\pm\hat{S}_{-}^{3})\hat{S}_{z}+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{4,\pm2} & =\phi_{\pm}(\hat{S}_{+}^{2}\pm\hat{S}_{-}^{2})(7\hat{S}_{z}^{2}-(X+5))+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{4,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})(7\hat{S}_{z}^{3}-(3X+1)\hat{S}_{z})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{4,0} & =35\hat{S}_{z}^{4}-(30X-25)\hat{S}_{z}^{2}+(3X^{2}-6X)\\ \\ \hat{\mathcal{O}}_{5,\pm5} & =\phi_{\pm}(\hat{S}_{+}^{5}\pm\hat{S}_{-}^{5})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{5,\pm4} & =\phi_{\pm}(\hat{S}_{+}^{4}\pm\hat{S}_{-}^{4})\hat{S}_{z}+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{5,\pm3} & =\phi_{\pm}(\hat{S}_{+}^{3}\pm\hat{S}_{-}^{3})(9\hat{S}_{z}^{2}-(X+33/2))+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{5,\pm2} & =\phi_{\pm}(\hat{S}_{+}^{2}\pm\hat{S}_{-}^{2})(3\hat{S}_{z}^{3}-(X+6)\hat{S}_{z})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{5,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})(21\hat{S}_{z}^{4}-14X\hat{S}_{z}^{2}+(X^{2}-X+3/2))+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{5,0} & =63\hat{S}_{z}^{5}-(70X-105)\hat{S}_{z}^{3}+(15X^{2}-50X+12)\hat{S}_{z}\\ \\ \hat{\mathcal{O}}_{6,\pm6} & =\phi_{\pm}(\hat{S}_{+}^{6}\pm\hat{S}_{-}^{6})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,\pm5} & =\phi_{\pm}(\hat{S}_{+}^{5}\pm\hat{S}_{-}^{5})\hat{S}_{z}+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,\pm4} & =\phi_{\pm}(\hat{S}_{+}^{4}\pm\hat{S}_{-}^{4})(11\hat{S}_{z}^{2}-X-38)+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,\pm3} & =\phi_{\pm}(\hat{S}_{+}^{3}\pm\hat{S}_{-}^{3})(11\hat{S}_{z}^{3}-(3X+59)\hat{S}_{z})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,\pm2} & =\phi_{\pm}(\hat{S}_{+}^{2}\pm\hat{S}_{-}^{2})(33\hat{S}_{z}^{4}-(18X+123)\hat{S}_{z}^{2}+X^{2}+10X+102)+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,\pm1} & =\phi_{\pm}(\hat{S}_{+}\pm\hat{S}_{-})(33\hat{S}_{z}^{5}-(30X-15)\hat{S}_{z}^{3}+(5X^{2}-10X+12)\hat{S}_{z})+\mathrm{h.c.}\\ \hat{\mathcal{O}}_{6,0} & =231\hat{S}_{z}^{6}-(315X-735)\hat{S}_{z}^{4}+(105X^{2}-525X+294)\hat{S}_{z}^{2}-5X^{3}+40X^{2}-60X \end{align*}\]

Computer-generated tables of Stevens operators with $k > 6$ are available from C. Rudowicz and C. Y. Chung, J. Phys.: Condens. Matter 16, 5825 (2004), but these typically do not appear in magnetic simulations.

The case $k=1$ gives the dipole operators,

\[(\hat{\mathcal{O}}_{1,1}, \hat{\mathcal{O}}_{1,0}, \hat{\mathcal{O}}_{1,-1}) = (\hat{S}_{x}, \hat{S}_{z}, \hat{S}_{y}).\]

The case $k=2$ gives the quadrupole operators,

\[(\hat{\mathcal{O}}_{2,2}, \dots, \hat{\mathcal{O}}_{2,-2}) = \left(\hat{S}_x^2 - \hat{S}_y^2, \frac{\hat{S}_x \hat{S}_z + \hat{S}_z \hat{S}_x}{2}, 2\hat{S}_z^2-\hat{S}_x^2-\hat{S}_y^2, \frac{\hat{S}_y \hat{S}_z + \hat{S}_z \hat{S}_y}{2}, \hat{S}_x \hat{S}_y + \hat{S}_y \hat{S}_x\right).\]

For each $k$ value, the set of operators $\hat{\mathcal{O}}_{k,q'}$ for $q' = -k, \dots, k$ form an irreducible representation of the group of rotations O(3). That is, rotation will transform $\hat{\mathcal{O}}_{k,q}$ into a linear combination of $\hat{\mathcal{O}}_{k,q'}$ where $q'$ varies but $k$ remains fixed.

In taking the large-$s$ limit, each dipole operator is replaced by its expectation value $\boldsymbol{\Omega} = \langle \hat{\mathbf{S}} \rangle$, and only leading-order terms are retained. The operator $\hat{\mathcal{O}}_{k,q}$ becomes a homogeneous polynomial $O_{k,q}(\boldsymbol{\Omega})$ of order $k$ in the spin components $\Omega^\alpha$. One can see these polynomials by constructing stevens_matrices with the argument s = Inf. Due to the normalization constraint, each dipole can be expressed in polar angles, $(\theta, \phi)$. Then the Stevens functions $O_{k,q}(\boldsymbol{\Omega})$ correspond to the spherical harmonic functions $Y_l^m(\theta, \phi)$ where $l=k$ and $m=q$; this correspondence is valid up to $k$ and $q$-dependent rescaling factors.