diff --git a/Project.toml b/Project.toml
index 3d334cd40..8f316ac4a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.4.4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CrystalInfoFramework = "6007d9b0-c6b2-11e8-0510-1d10e825f3f1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FilePathsBase = "48062228-2e41-5def-b9a4-89aafe57970f"
 Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
@@ -29,7 +28,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Colors = "0.12.8"
 CrystalInfoFramework = "0.5.0"
 DataStructures = "0.18.13"
-DynamicPolynomials = "0.4.6"
 FFTW = "1.4.5"
 FilePathsBase = "0.9.18"
 Inflate = "0.1.2"
diff --git a/docs/Project.toml b/docs/Project.toml
index d73436c04..f8cf1f1fb 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
diff --git a/docs/make.jl b/docs/make.jl
index 21f968781..e342df9b4 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,6 +2,9 @@
 
 using Literate, Documenter, Sunny
 
+import DynamicPolynomials # get symbolic functions
+import GLMakie # get plotting functions
+
 execute = true # set `false` to disable cell evaluation
 
 example_names = ["fei2_tutorial", "powder_averaging", "ising2d"]
diff --git a/docs/src/quick-start.md b/docs/src/quick-start.md
index be281a07b..9f3867edc 100644
--- a/docs/src/quick-start.md
+++ b/docs/src/quick-start.md
@@ -105,7 +105,7 @@ The next steps are typically the following
    crystal unit cells.
 2. Add interactions to the system using functions like
    [`set_external_field!`](@ref), [`set_exchange!`](@ref), and
-   [`set_anisotropy!`](@ref).
+   [`set_onsite!`](@ref).
 3. Perform Monte Carlo simulation to equilibrate the spin configuration. Options
    include the continuous [`Langevin`](@ref) dynamics, or single-spin flip
    updates with [`LocalSampler`](@ref). The former can efficiently handle
diff --git a/docs/src/versions.md b/docs/src/versions.md
index cd3bafc7e..b0453e6e9 100644
--- a/docs/src/versions.md
+++ b/docs/src/versions.md
@@ -1,3 +1,22 @@
+# Version 0.5.0
+
+This version includes many **breaking changes**.
+
+Replace `set_anisotropy!` with a new function [`set_onsite!`](@ref). The latter
+expects an explicit matrix representation for the local Hamiltonian. This can be
+constructed, e.g., as a linear combination of [`stevens_operators`](@ref), or as
+a polynomial of [`spin_operators`](@ref). To understand the mapping between
+these two, the new function [`print_stevens_expansion`](@ref) acts on an
+arbitrary local operator.
+
+Symbolic representations of operators are now hidden unless the package
+`DynamicPolynomials` is explicitly loaded by the user. The functionality of
+`print_anisotropy_as_stevens` has been replaced with
+[`print_classical_stevens_expansion`](@ref), while
+`print_anisotropy_as_classical_spins` has become
+[`print_classical_spin_polynomial`](@ref).
+
+
 # Version 0.4.3
 
 **Experimental** support for linear [`SpinWaveTheory`](@ref), implemented in
@@ -50,7 +69,7 @@ This update includes many breaking changes, and is missing some features of
 
 ### Creating a spin `System`
 
-`SpinSystem` has been renamed [`System`](@ref). Its constructor now has the form,
+Rename `SpinSystem` to [`System`](@ref). Its constructor now has the form,
 
 ```julia
 System(crystal, latsize, infos, mode)
@@ -66,13 +85,13 @@ The parameter `mode` is one of `:SUN` or `:dipole`.
 
 Interactions are now added mutably to an existing `System` using the following
 functions: [`set_external_field!`](@ref), [`set_exchange!`](@ref),
-[`set_anisotropy!`](@ref), [`enable_dipole_dipole!`](@ref).
+[`set_onsite!`](@ref), [`enable_dipole_dipole!`](@ref).
 
 As a convenience, one can use [`dmvec(D)`](@ref) to convert a DM vector to a
 $3×3$ antisymmetric exchange matrix.
 
 Fully general single-ion anisotropy is now possible. The function
-[`set_anisotropy!`](@ref) expects the single ion anisotropy to be expressed as a
+[`set_onsite!`](@ref) expects the single ion anisotropy to be expressed as a
 polynomial in symbolic spin operators [`𝒮`](@ref), or as a linear combination
 of symbolic Stevens operators [`𝒪`](@ref). For example, an easy axis anisotropy
 in the direction `n` may be written `D*(𝒮⋅n)^2`.
@@ -80,7 +99,7 @@ in the direction `n` may be written `D*(𝒮⋅n)^2`.
 Stevens operators `𝒪[k,q]` admit polynomial expression in spin operators
 `𝒮[α]`. Conversely, a polynomial of spin operators can be expressed as a linear
 combination of Stevens operators. To see this expansion use
-[`print_anisotropy_as_stevens`](@ref).
+`print_anisotropy_as_stevens`.
 
 
 ### Inhomogeneous field
diff --git a/examples/fei2_tutorial.jl b/examples/fei2_tutorial.jl
index f8dc48f0f..e025d0636 100644
--- a/examples/fei2_tutorial.jl
+++ b/examples/fei2_tutorial.jl
@@ -25,7 +25,6 @@
 # Begin by importing `Sunny` and `GLMakie`, a plotting package.
 
 using Sunny, GLMakie
-#md Makie.inline!(true); #hide
 #nb Sunny.offline_viewers()  # Inject Javascript code for additional plotting capabilities 
 
 # If you see an error `Package <X> not found in current path`, add the package
@@ -145,18 +144,17 @@ set_exchange!(sys, [J′2apm 0.0    0.0;
                     0.0    J′2apm 0.0;
                     0.0    0.0    J′2azz], Bond(1,1,[1,2,1]))
 
-# The function [`set_anisotropy!`](@ref) assigns a single-ion anisotropy. It
-# takes an abstract operator and an atom index. The operator may be a polynomial
-# of spin operators or a linear combination of Stevens operators. Sunny provides
-# special symbols for their construction: [`𝒮`](@ref) is a vector of the three
-# spin operators and [`𝒪`](@ref) are the symbolic Stevens operators. Here we
-# construct an easy-axis anisotropy.
+# The function [`set_onsite!`](@ref) assigns a single-ion anisotropy
+# operator. It can be constructed, e.g., from the matrices given by
+# [`spin_operators`](@ref) or [`stevens_operators`](@ref). Here we construct an
+# easy-axis anisotropy along the direction $\hat{z}$.
 
 D = 2.165
-set_anisotropy!(sys, -D*𝒮[3]^2, 1)
+S = spin_operators(sys, 1)
+set_onsite!(sys, -D*S[3]^2, 1)
 
 # Any anisotropy operator can be converted to a linear combination of Stevens
-# operators with [`print_anisotropy_as_stevens`](@ref).
+# operators with [`print_stevens_expansion`](@ref).
 
 # # Calculating structure factor intensities
 # In the remainder of this tutorial, we will examine Sunny's tools for
diff --git a/examples/powder_averaging.jl b/examples/powder_averaging.jl
index a00389e63..451817abf 100644
--- a/examples/powder_averaging.jl
+++ b/examples/powder_averaging.jl
@@ -7,7 +7,6 @@
 
 using Sunny, GLMakie
 using Statistics: mean
-#md Makie.inline!(true); #hide
 
 dims = (8,8,8)               # Lattice dimensions
 seed = 1                     # RNG seed for repeatable behavior 
diff --git a/src/Operators/Spin.jl b/src/Operators/Spin.jl
index 79a1969ba..8e8167ea5 100644
--- a/src/Operators/Spin.jl
+++ b/src/Operators/Spin.jl
@@ -1,4 +1,10 @@
-# Construct spin operators, i.e. generators of su(2), of dimension N
+"""
+    spin_matrices(N)
+
+Constructs the three spin operators, i.e. the generators of SU(2), in the
+`N`-dimensional irrep. See also [`spin_operators`](@ref), which determines the
+appropriate value of `N` for a given site index.
+"""
 function spin_matrices(N::Int)
     if N == 0
         return fill(zeros(ComplexF64,0,0), 3)
diff --git a/src/Operators/Stevens.jl b/src/Operators/Stevens.jl
index 4ba7c130b..7bf93ff85 100644
--- a/src/Operators/Stevens.jl
+++ b/src/Operators/Stevens.jl
@@ -167,3 +167,66 @@ function rotate_stevens_coefficients(c, R::Mat3)
     @assert norm(imag(c′)) < 1e-12
     return real(c′)
 end
+
+
+# Builds matrix representation of Stevens operators on demand
+struct StevensMatrices
+    N::Int
+end
+function Base.getindex(this::StevensMatrices, k::Int, q::Int)
+    k < 0  && error("Stevens operators 𝒪[k,q] require k >= 0.")
+    k > 6  && error("Stevens operators 𝒪[k,q] currently require k <= 6.")
+    !(-k <= q <= k) && error("Stevens operators 𝒪[k,q] require -k <= q <= k.")
+    if k == 0
+        return Matrix{ComplexF64}(I, this.N, this.N)
+    else
+        # Stevens operators are stored in descending order: k, k-1, ... -k.
+        # Therefore q=k should be the first index, q_idx=1.
+        q_idx = k - q + 1
+        # Build all (2k+1) matrices and extract the desired one. This seems a
+        # bit wasteful, but presumably this function won't be called in
+        # performance sensitive contexts.
+        return stevens_matrices(k; this.N)[q_idx]
+    end
+end
+
+
+"""
+    function print_stevens_expansion(op)
+
+Prints a local Hermitian operator as a linear combination of Stevens operators.
+This function works on explicit matrix representations. The analogous function
+in the large-``S`` classical limit is
+[`print_classical_stevens_expansion`](@ref).
+
+# Examples
+
+```julia
+S = spin_matrices(5)
+print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
+# Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + 102/5
+```
+"""
+function print_stevens_expansion(op::Matrix{ComplexF64})
+    op ≈ op' || error("Requires Hermitian operator")
+    terms = String[]
+
+    # Decompose op into Stevens coefficients
+    c = matrix_to_stevens_coefficients(op)
+    for k in 1:6
+        for (c_km, m) in zip(reverse(c[k]), -k:k)
+            abs(c_km) < 1e-12 && continue
+            push!(terms, *(coefficient_to_math_string(c_km), "𝒪", int_to_underscore_string.((k,m))...))
+        end
+    end
+
+    # Handle linear shift specially
+    c_00 = real(tr(op))/size(op, 1)
+    abs(c_00) > 1e-12 && push!(terms, number_to_math_string(c_00))
+
+    # Concatenate with plus signs
+    str = join(terms, " + ")
+    # Remove redundant plus signs and print
+    str = replace(str, "+ -" => "- ")
+    println(str)
+end
diff --git a/src/Operators/Symbolic.jl b/src/Operators/Symbolic.jl
index f30d4d5e5..838f94263 100644
--- a/src/Operators/Symbolic.jl
+++ b/src/Operators/Symbolic.jl
@@ -1,3 +1,5 @@
+const DP = DynamicPolynomials
+
 # It is convenient to present Stevens operators to the user in ascending order
 # for the index q = -k...k. Internally, however, the symbols must be stored in
 # descending order q = k...-k for consistency with the basis used for spin
@@ -5,17 +7,17 @@
 # to generate rotations of the Stevens operators via the Wigner D matrices.
 const stevens_operator_symbols = let
     # 𝒪₀ = identity
-    𝒪₁ = collect(reverse(DP.@ncpolyvar                          𝒪₁₋₁ 𝒪₁₀ 𝒪₁₁))
-    𝒪₂ = collect(reverse(DP.@ncpolyvar                     𝒪₂₋₂ 𝒪₂₋₁ 𝒪₂₀ 𝒪₂₁ 𝒪₂₂))
-    𝒪₃ = collect(reverse(DP.@ncpolyvar                𝒪₃₋₃ 𝒪₃₋₂ 𝒪₃₋₁ 𝒪₃₀ 𝒪₃₁ 𝒪₃₂ 𝒪₃₃))
-    𝒪₄ = collect(reverse(DP.@ncpolyvar           𝒪₄₋₄ 𝒪₄₋₃ 𝒪₄₋₂ 𝒪₄₋₁ 𝒪₄₀ 𝒪₄₁ 𝒪₄₂ 𝒪₄₃ 𝒪₄₄))
-    𝒪₅ = collect(reverse(DP.@ncpolyvar      𝒪₅₋₅ 𝒪₅₋₄ 𝒪₅₋₃ 𝒪₅₋₂ 𝒪₅₋₁ 𝒪₅₀ 𝒪₅₁ 𝒪₅₂ 𝒪₅₃ 𝒪₅₄ 𝒪₅₅))
-    𝒪₆ = collect(reverse(DP.@ncpolyvar 𝒪₆₋₆ 𝒪₆₋₅ 𝒪₆₋₄ 𝒪₆₋₃ 𝒪₆₋₂ 𝒪₆₋₁ 𝒪₆₀ 𝒪₆₁ 𝒪₆₂ 𝒪₆₃ 𝒪₆₄ 𝒪₆₅ 𝒪₆₆))
+    𝒪₁ = collect(DP.@ncpolyvar                     𝒪₁₁ 𝒪₁₀ 𝒪₁₋₁)
+    𝒪₂ = collect(DP.@ncpolyvar                 𝒪₂₂ 𝒪₂₁ 𝒪₂₀ 𝒪₂₋₁ 𝒪₂₋₂)
+    𝒪₃ = collect(DP.@ncpolyvar             𝒪₃₃ 𝒪₃₂ 𝒪₃₁ 𝒪₃₀ 𝒪₃₋₁ 𝒪₃₋₂ 𝒪₃₋₃)
+    𝒪₄ = collect(DP.@ncpolyvar         𝒪₄₄ 𝒪₄₃ 𝒪₄₂ 𝒪₄₁ 𝒪₄₀ 𝒪₄₋₁ 𝒪₄₋₂ 𝒪₄₋₃ 𝒪₄₋₄)
+    𝒪₅ = collect(DP.@ncpolyvar     𝒪₅₅ 𝒪₅₄ 𝒪₅₃ 𝒪₅₂ 𝒪₅₁ 𝒪₅₀ 𝒪₅₋₁ 𝒪₅₋₂ 𝒪₅₋₃ 𝒪₅₋₄ 𝒪₅₋₅)
+    𝒪₆ = collect(DP.@ncpolyvar 𝒪₆₆ 𝒪₆₅ 𝒪₆₄ 𝒪₆₃ 𝒪₆₂ 𝒪₆₁ 𝒪₆₀ 𝒪₆₋₁ 𝒪₆₋₂ 𝒪₆₋₃ 𝒪₆₋₄ 𝒪₆₋₅ 𝒪₆₋₆)
     [𝒪₁, 𝒪₂, 𝒪₃, 𝒪₄, 𝒪₅, 𝒪₆]
 end
 
 const spin_operator_symbols = let
-    SVector{3}(DP.@ncpolyvar 𝒮₁ 𝒮₂ 𝒮₃)
+    SVector{3}(reverse(DP.@ncpolyvar 𝒮₃ 𝒮₂ 𝒮₁))
 end
 
 const spin_squared_symbol = let
@@ -40,6 +42,7 @@ function Base.getindex(::StevensOpsAbstract, k::Int, q::Int)
     end
 end
 
+
 """
     𝒪[k,q]
 
@@ -66,7 +69,7 @@ function operator_to_matrix(p::DP.AbstractPolynomialLike; N)
     if !(rep ≈ rep')
         println("Warning: Symmetrizing non-Hermitian operator '$p'.")
     end
-    # Symmetrize in any case for more accuracy
+    # Symmetrize in any case for slightly more accuracy
     return (rep+rep')/2
 end
 function operator_to_matrix(p::Number; N)
@@ -74,7 +77,7 @@ function operator_to_matrix(p::Number; N)
 end
 
 
-##### Conversion of spin polynomial to linear combination of 'Stevens functions' #####
+##### Convert classical spin polynomial to linear combination of 'Stevens functions' #####
 
 # Construct Stevens operators in the classical limit, represented as polynomials
 # of spin expectation values
@@ -85,7 +88,7 @@ function stevens_classical(k::Int)
         # homogeneous polynomial of degree k
         𝒪 = sum(t for t in 𝒪 if DP.degree(t) == k)
         # Remaining coefficients must be real integers; make this explicit
-        𝒪 = DP.mapcoefficients(x -> Int(x), 𝒪)
+        𝒪 = DP.map_coefficients(x -> Int(x), 𝒪)
         return 𝒪
     end
 end
@@ -166,49 +169,16 @@ function operator_to_classical_stevens_coefficients(p, S)
     end
 end
 
-function rotate_operator(P::DP.AbstractPolynomialLike, R)
-    R = convert(Mat3, R)
-
-    # The spin operator vector rotates two equivalent ways:
-    #  1. S_α -> U' S_α U
-    #  2. S_α -> R_αβ S_β
-    #
-    # where U = exp(-i θ n⋅S), for the axis-angle rotation (n, θ). Apply the
-    # latter to transform symbolic spin operators.
-    𝒮′ = R * 𝒮
-
-    # Spherical tensors T_q rotate two equivalent ways:
-    #  1. T_q -> U' T_q U       (with U = exp(-i θ n⋅S) in dimension N irrep)
-    #  2. T_q -> D*_{q,q′} T_q′ (with D = exp(-i θ n⋅S) in dimension 2k+1 irrep)
-    #
-    # The Stevens operators 𝒪_q are linearly related to T_q via 𝒪 = α T.
-    # Therefore rotation on Stevens operators is 𝒪 -> α conj(D) α⁻¹ 𝒪.
-    local 𝒪 = stevens_operator_symbols
-    𝒪′ = map(𝒪) do 𝒪ₖ
-        k = Int((length(𝒪ₖ)-1)/2)
-        D = unitary_for_rotation(R; N=2k+1)
-        R_stevens = stevens_α[k] * conj(D) * stevens_αinv[k]
-        @assert norm(imag(R_stevens)) < 1e-12
-        real(R_stevens) * 𝒪ₖ
-    end
-
-    # Spin squared as a scalar may be introduced through
-    # operator_to_classical_stevens()
-    X = spin_squared_symbol
-
-    # Perform substitutions
-    P′ = P(𝒮 => 𝒮′, [𝒪[k] => 𝒪′[k] for k=1:6]..., X => X)
-
-    # Remove terms very near zero
-    return DP.mapcoefficients(P′) do c
-        abs(c) < 1e-12 ? zero(c) : c
-    end
-end
 
 ##### Printing of operators #####
 
 function pretty_print_operator(p::DP.AbstractPolynomialLike)
-    terms = map(zip(DP.coefficients(p), DP.monomials(p))) do (c, m)
+    # Iterator over coefficients and monomials
+    terms = zip(DP.coefficients(p), DP.monomials(p))
+    # Keep only terms with non-vanishing coefficients
+    terms = Iterators.filter(x -> abs(x[1]) > 1e-12, terms)
+    # Pretty-print each term
+    terms = map(terms) do (c, m)
         isone(m) ? number_to_math_string(c) : coefficient_to_math_string(c) * repr(m)
     end
     # Concatenate with plus signs
@@ -223,65 +193,64 @@ end
 
 
 """
-    function print_anisotropy_as_classical_spins(p)
-
-Prints a quantum operator (e.g. linear combination of Stevens operators) as a
-polynomial of spin expectation values in the classical limit.
-
-See also [`print_anisotropy_as_stevens`](@ref).
+    function print_classical_spin_polynomial(op)
+
+Prints an operator (e.g., a linear combination of Stevens operators `𝒪`) as a
+polynomial in the classical dipole components.
+            
+This function works in the "large-``S``" classical limit, which corresponds to
+replacing each spin operator with its expected dipole. There are ambiguities in
+defining this limit at sub-leading order in ``S``. The procedure in Sunny is as
+follows: First, uniquely decompose `op` as a linear combination of Stevens
+operators, each of which is defined as a polynomial in the spin operators. To
+take the large-``S`` limit, Sunny replaces each Stevens operator with a
+corresponding polynomial in the expected dipole components, keeping only leading
+order terms in ``S``. The resulting "classical Stevens functions" are
+essentially the spherical harmonics ``Yᵐₗ``, up to ``m``- and ``l``- dependent
+scaling factors.
+
+# Example
+
+```julia
+using DynamicPolynomials, Sunny
+
+print_classical_spin_polynomial((1/4)𝒪[4,4] + (1/20)𝒪[4,0] + (3/5)*(𝒮'*𝒮)^2)
+# Prints: 𝒮₁⁴ + 𝒮₂⁴ + 𝒮₃⁴
+```
+
+See also [`print_classical_stevens_expansion`](@ref) for the inverse mapping.
 """
-function print_anisotropy_as_classical_spins(p)
-    p = operator_to_classical_polynomial(p)
-    p = p(spin_classical_symbols => 𝒮)
-    pretty_print_operator(p)
+function print_classical_spin_polynomial(op)
+    op = operator_to_classical_polynomial(op)
+    op = op(spin_classical_symbols => 𝒮)
+    pretty_print_operator(op)
 end
 
 """
-    function print_anisotropy_as_stevens(p; N)
+    function print_classical_stevens_expansion(op)
 
-Prints a quantum operator (e.g. a polynomial of the spin operators `𝒮`) as a
-linear combination of Stevens operators. The parameter `N` specifies the
-dimension of the SU(_N_) representation, corresponding to quantum spin magnitude
-``S = (N-1)/2``. The special value `N = 0` indicates the large-``S`` classical
-limit.
+Prints an operator (e.g. a polynomial of the spin operators `𝒮`) as a linear
+combination of Stevens operators. This function works in the large-``S``
+classical limit, as described in the documentation for
+[`print_classical_spin_polynomial`](@ref).
 
-In the output, the symbol `X` denotes the spin operator magnitude squared.
-Quantum spin operators ``𝒮`` of any finite dimension satisfy ``X = |𝒮|^2 = S
-(S+1)``. To take the large-``S`` limit, however, we keep only leading order
-powers of ``S``, such that ``X = S^2``.
+In the output, the symbol `X` denotes the spin magnitude squared, which can be
+entered symbolically as `𝒮'*𝒮`.
 
-This function can be useful for understanding the conversions performed
-internally by [`set_anisotropy!`](@ref).
+# Examples
 
-For the inverse mapping, see [`print_anisotropy_as_classical_spins`](@ref).
-"""
-function print_anisotropy_as_stevens(p; N)
-    if N == 0
-        p′ = operator_to_classical_stevens(p)
-    else
-        Λ = operator_to_matrix(p; N)
-
-        # Stevens operators are orthogonal but not normalized. Pull out
-        # coefficients c one-by-one and accumulate into p′. These must be real
-        # because both Λ and 𝒪 are Hermitian.
-
-        # k = 0 term, for which 𝒪₀₀ = I.
-        p′ = real(tr(Λ)/N)
-
-        # Stevens operators are zero when k >= N
-        for k = 1:min(6, N-1)
-            for (𝒪mat, 𝒪sym) = zip(stevens_matrices(k; N), stevens_operator_symbols[k])
-                # See also: `matrix_to_stevens_coefficients`
-                c = real(tr(𝒪mat'*Λ) / tr(𝒪mat'*𝒪mat))
-                if abs(c) > 1e-12
-                    p′ += c*𝒪sym
-                end
-            end
-        end
+```julia
+using DynamicPolynomials, Sunny
 
-        # p′ should be faithful to p and its matrix representation Λ. This will
-        # fail if the spin polynomial order in p exceeds 6.
-        @assert operator_to_matrix(p′; N) ≈ Λ
-    end
-    pretty_print_operator(p′)
+print_classical_stevens_expansion(𝒮[1]^4 + 𝒮[2]^4 + 𝒮[3]^4)
+# Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + (3/5)X²
+```
+
+See also [`print_classical_spin_polynomial`](@ref) for the inverse mapping.
+
+The function [`print_stevens_expansion`](@ref) is analogous to this one, but
+expects a quantum operator in a finite-``S`` representation.
+"""
+function print_classical_stevens_expansion(op)
+    pretty_print_operator(operator_to_classical_stevens(op))
 end
diff --git a/src/Operators/TensorOperators.jl b/src/Operators/TensorOperators.jl
new file mode 100644
index 000000000..fa7ac7bac
--- /dev/null
+++ b/src/Operators/TensorOperators.jl
@@ -0,0 +1,126 @@
+# Produces a matrix representation of a tensor product of operators, C = A⊗B.
+# Like built-in `kron` but with permutation. Returns C_{acbd} = A_{ab} B_{cd}.
+function kron_operator(A::AbstractMatrix, B::AbstractMatrix)
+    TS = promote_type(eltype(A), eltype(B))
+    C = zeros(TS, size(A,1), size(B,1), size(A,2), size(B,2))
+    for ci in CartesianIndices(C)
+        a, c, b, d = Tuple(ci)
+        C[ci] = A[a,b] * B[c,d]
+    end
+    return reshape(C, size(A,1)*size(B,1), size(A,2)*size(B,2))
+end
+
+
+function degeneracy_groups(S, tol)
+    acc = UnitRange{Int}[]
+    isempty(S) && return acc
+
+    j0 = 1
+    for j = 2:lastindex(S)
+        if abs(S[j0] - S[j]) > tol
+            push!(acc, j0:(j-1))
+            j0 = j
+        end
+    end
+
+    push!(acc, j0:length(S))
+    return acc
+end
+
+# Use SVD to find the decomposition D = ∑ₖ Aₖ ⊗ Bₖ, where Aₖ is N₁×N₁ and Bₖ is
+# N₂×N₂. Returns the list of matrices [(A₁, B₁), (A₂, B₂), ...].
+function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
+    tol = 1e-12
+
+    @assert size(D, 1) == size(D, 2) == N1*N2
+    D̃ = permutedims(reshape(D, N1, N2, N1, N2), (1,3,2,4))
+    D̃ = reshape(D̃, N1*N1, N2*N2)
+    (; S, U, V) = svd(D̃)
+    ret = []
+
+    # Rotate columns of U and V within each degenerate subspace so that all
+    # columns (when reshaped) are Hermitian matrices
+    for range in degeneracy_groups(S, tol)
+        abs(S[first(range)]) < tol && break
+        
+        U_sub = view(U, :, range)
+        n = length(range)
+        Q = zeros(n, n)
+        for k in 1:n, k′ in 1:n
+            uk  = reshape(view(U_sub, :, k), N1, N1)
+            uk′ = reshape(view(U_sub, :, k′), N1, N1)
+            Q[k, k′] = conj(tr(uk * uk′))
+        end
+        
+        R = sqrt(Q)
+        @assert norm(R*R' - I) < 1e-12
+
+        @views U[:, range] = U[:, range] * R
+        @views V[:, range] = V[:, range] * R
+    end
+
+    # Check that rotation was valid
+    @assert U*diagm(S)*V' ≈ D̃
+
+    for (k, σ) in enumerate(S)
+        if abs(σ) > tol
+            u = reshape(U[:, k], N1, N1)
+            v = reshape(V[:, k], N2, N2)
+            # Check factors are really Hermitian
+            @assert norm(u - u') < tol
+            @assert norm(v - v') < tol
+            u = Hermitian(u+u')/2
+            v = Hermitian(v+v')/2
+            push!(ret, (σ*u, conj(v)))
+        end
+    end
+    return ret
+end
+
+function local_quantum_operators(A, B)
+    (isempty(A) || isempty(B)) && error("Nonempty lists required")
+
+    @assert allequal(size.(A))
+    @assert allequal(size.(B))
+
+    N1 = size(first(A), 1)
+    N2 = size(first(B), 1)
+
+    I1 = Ref(Matrix(I, N1, N1))
+    I2 = Ref(Matrix(I, N2, N2))
+    return (kron_operator.(A, I2), kron_operator.(I1, B))
+end
+
+
+#=
+
+# Returns the spin operators for two sites, with Hilbert space dimensions N₁ and
+# N₂, respectively.
+function spin_pair(N1, N2)
+    S1 = Sunny.spin_matrices(N1)
+    S2 = Sunny.spin_matrices(N2)
+    return local_quantum_operators(S1, S2)
+end
+
+N1 = 3
+N2 = 3
+
+# Check Kronecker product of operators
+A1 = randn(N1,N1)
+A2 = randn(N1,N1)
+B1 = randn(N2,N2)
+B2 = randn(N2,N2)
+@assert kron_operator(A1, B1) * kron_operator(A2, B2) ≈ kron_operator(A1*A2, B1*B2)
+
+# Check SVD decomposition
+S1, S2 = spin_pair(N1, N2)
+B = (S1' * S2)^2                                # biquadratic interaction
+D = svd_tensor_expansion(B, N1, N2)             # a sum of 9 tensor products
+@assert sum(kron_operator(d...) for d in D) ≈ B # consistency check
+
+
+B = S1' * S2
+D = svd_tensor_expansion(B, N1, N2)             # a sum of 9 tensor products
+@assert sum(kron_operator(d...) for d in D) ≈ B
+
+=#
\ No newline at end of file
diff --git a/src/Sunny.jl b/src/Sunny.jl
index 11253168e..4f43e1b9c 100644
--- a/src/Sunny.jl
+++ b/src/Sunny.jl
@@ -9,7 +9,6 @@ import FFTW
 import ProgressMeter: Progress, next!
 import Printf: @printf, @sprintf
 import Random: Random, randn!
-import DynamicPolynomials as DP
 import DataStructures: SortedDict
 import Optim
 
@@ -35,8 +34,8 @@ include("Util/CartesianIndicesShifted.jl")
 include("Operators/Spin.jl")
 include("Operators/Rotation.jl")
 include("Operators/Stevens.jl")
-include("Operators/Symbolic.jl")
-export 𝒪, 𝒮, rotate_operator, print_anisotropy_as_stevens, print_anisotropy_as_classical_spins
+include("Operators/TensorOperators.jl")
+export spin_matrices, rotate_operator, print_stevens_expansion
 
 include("Symmetry/LatticeUtils.jl")
 include("Symmetry/SymOp.jl")
@@ -50,9 +49,6 @@ include("Symmetry/Printing.jl")
 export Crystal, subcrystal, lattice_vectors, lattice_params, Bond, 
     reference_bonds, print_site, print_bond, print_symmetry_table,
     print_suggested_frame
-    # natoms, cell_volume, cell_type, coordination_number, displacement, distance,
-    # all_symmetry_related_bonds, all_symmetry_related_bonds_for_atom,
-    # all_symmetry_related_couplings, all_symmetry_related_couplings_for_atom
 
 include("Units.jl")
 export meV_per_K, Units
@@ -66,8 +62,8 @@ include("System/Ewald.jl")
 include("System/Interactions.jl")
 export SpinInfo, System, Site, all_sites, position_to_site,
     global_position, magnetic_moment, polarize_spin!, polarize_spins!, randomize_spins!, energy, forces,
-    set_external_field!, set_anisotropy!, set_exchange!, set_biquadratic!, dmvec, enable_dipole_dipole!,
-    to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_anisotropy_at!,
+    spin_operators, stevens_operators, set_external_field!, set_onsite!, set_exchange!, set_biquadratic!,
+    dmvec, enable_dipole_dipole!, to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_anisotropy_at!,
     symmetry_equivalent_bonds, set_exchange_at!, set_biquadratic_at!, remove_periodicity!
 
 include("Reshaping.jl")
@@ -110,12 +106,17 @@ include("MonteCarlo/WangLandau.jl")
 include("MonteCarlo/ParallelWangLandau.jl")
 export propose_uniform, propose_flip, propose_delta, @mix_proposals, LocalSampler
 
-# Makie (e.g., WGLMakie or GLMakie) is an optional dependency
 function __init__()
+    # Importing Makie (e.g., WGLMakie or GLMakie) will enable plotting
     @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
         include("Plotting.jl")
         export plot_spins
-            # plot_lattice, plot_bonds, plot_all_bonds, anim_integration, live_integration, live_langevin_integration
+    end
+
+    # Importing DynamicPolynomials will enable certain symbolic analysis
+    @require DynamicPolynomials="7c1d4256-1411-5781-91ec-d7bc3513ac07" begin
+        include("Operators/Symbolic.jl")
+        export 𝒪, 𝒮, print_classical_stevens_expansion, print_classical_spin_polynomial
     end
 end
 
diff --git a/src/Symmetry/Printing.jl b/src/Symmetry/Printing.jl
index 95c00b5c1..f04df3ff0 100644
--- a/src/Symmetry/Printing.jl
+++ b/src/Symmetry/Printing.jl
@@ -310,7 +310,7 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits,
             push!(lines, prefix * join(terms, " + "))
         end
     end
-    println("Allowed anisotropy in Stevens operators 𝒪[k,q]:")
+    println("Allowed anisotropy in Stevens operators:")
     println(join(lines, " +\n"))
 
     if R != I
diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl
index a0271f050..a544b9c6a 100644
--- a/src/System/PairExchange.jl
+++ b/src/System/PairExchange.jl
@@ -260,7 +260,7 @@ due to possible periodic wrapping. Resolve this ambiguity by passing an explicit
 
 See also [`set_exchange!`](@ref).
 """
-function set_exchange_at!(sys::System{N}, J, site1, site2; offset=nothing) where N
+function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; offset=nothing) where N
     site1 = to_cartesian(site1)
     site2 = to_cartesian(site2)
     bond = sites_to_internal_bond(sys, site1, site2, offset)
diff --git a/src/System/SingleIonAnisotropy.jl b/src/System/SingleIonAnisotropy.jl
index 197e426ee..08c465a3a 100644
--- a/src/System/SingleIonAnisotropy.jl
+++ b/src/System/SingleIonAnisotropy.jl
@@ -28,6 +28,27 @@ function SingleIonAnisotropy(sys::System, op, N)
     return SingleIonAnisotropy(matrep, stvexp)
 end
 
+function SingleIonAnisotropy(sys::System, matrep::Matrix{ComplexF64}, N)
+    # Remove trace
+    matrep -= (tr(matrep)/N)*I
+    if norm(matrep) < 1e-12
+        matrep = zero(matrep)
+    end
+    c = matrix_to_stevens_coefficients(matrep)
+
+    iszero(c[[1,3,5]]) || error("Single-ion anisotropy must be time-reversal invariant.")
+
+    if sys.mode == :dipole
+        λ = anisotropy_renormalization(N)
+        stvexp = StevensExpansion(λ[2]*c[2], λ[4]*c[4], λ[6]*c[6])
+    else
+        stvexp = StevensExpansion(c[2], c[4], c[6])
+    end
+    
+    return SingleIonAnisotropy(matrep, stvexp)
+end
+
+
 function anisotropy_renormalization(N)
     S = (N-1)/2
     return ((1), # k=1
@@ -78,15 +99,11 @@ end
 
 
 """
-    set_anisotropy!(sys::System, op, i::Int)
+    set_onsite!(sys::System, op::Matrix{ComplexF64}, i::Int)
 
 Set the single-ion anisotropy for the `i`th atom of every unit cell, as well as
-all symmetry-equivalent atoms. The parameter `op` may be a polynomial in
-symbolic spin operators `𝒮[α]`, or a linear combination of symbolic Stevens
-operators `𝒪[k,q]`.
-
-The characters `𝒮` and `𝒪` can be copy-pasted from this help message, or typed
-at a Julia terminal using `\\scrS` or `\\scrO` followed by tab-autocomplete.
+all symmetry-equivalent atoms. The local operator `op` may be constructed using
+[`spin_operators`](@ref) or [`stevens_operators`](@ref).
 
 For systems restricted to dipoles, the anisotropy operators interactions will
 automatically be renormalized to achieve maximum consistency with the more
@@ -95,25 +112,27 @@ variationally accurate SU(_N_) mode.
 # Examples
 ```julia
 # An easy axis anisotropy in the z-direction
-set_anisotropy!(sys, -D*𝒮[3]^3, i)
+S = spin_operators(sys, i)
+set_onsite!(sys, -D*S[3]^3, i)
 
 # The unique quartic single-ion anisotropy for a site with cubic point group
 # symmetry
-set_anisotropy!(sys, 𝒪[4,0] + 5𝒪[4,4], i)
+O = stevens_operators(sys, i)
+set_onsite!(sys, O[4,0] + 5O[4,4], i)
 
 # An equivalent expression of this quartic anisotropy, up to a constant shift
-set_anisotropy!(sys, 20*(𝒮[1]^4 + 𝒮[2]^4 + 𝒮[3]^4), i)
+set_onsite!(sys, 20*(S[1]^4 + S[2]^4 + S[3]^4), i)
 ```
 
-See also [`print_anisotropy_as_stevens`](@ref).
+See also [`spin_operators`](@ref).
 """
-function set_anisotropy!(sys::System{N}, op::DP.AbstractPolynomialLike, i::Int) where N
+function set_onsite!(sys::System{N}, op::Matrix{ComplexF64}, i::Int) where N
     is_homogeneous(sys) || error("Use `set_anisotropy_at!` for an inhomogeneous system.")
 
     # If `sys` has been reshaped, then operate first on `sys.origin`, which
     # contains full symmetry information.
     if !isnothing(sys.origin)
-        set_anisotropy!(sys.origin, op, i)
+        set_onsite!(sys.origin, op, i)
         set_interactions_from_origin!(sys)
         return
     end
@@ -159,15 +178,15 @@ end
 
 
 """
-    set_anisotropy_at!(sys::System, op, site::Site)
+    set_anisotropy_at!(sys::System, op::Matrix{ComplexF64}, site::Site)
 
 Sets the single-ion anisotropy operator `op` for a single [`Site`](@ref),
 ignoring crystal symmetry.  The system must support inhomogeneous interactions
 via [`to_inhomogeneous`](@ref).
 
-See also [`set_anisotropy!`](@ref).
+See also [`set_onsite!`](@ref).
 """
-function set_anisotropy_at!(sys::System{N}, op::DP.AbstractPolynomialLike, site) where N
+function set_anisotropy_at!(sys::System{N}, op::Matrix{ComplexF64}, site::Site) where N
     is_homogeneous(sys) && error("Use `to_inhomogeneous` first.")
     ints = interactions_inhomog(sys)
     site = to_cartesian(site)
@@ -175,7 +194,12 @@ function set_anisotropy_at!(sys::System{N}, op::DP.AbstractPolynomialLike, site)
 end
 
 
-# Evaluate a given linear combination of Stevens operators for classical spin s
+# Evaluate a given linear combination of Stevens operators in the large-S limit,
+# where each spin operator is replaced by its dipole expectation value. In this
+# limit, each Stevens operator O[ℓ,m](s) becomes a homogeneous polynomial in the
+# spin components sᵅ, and is equal to the spherical Harmonic Yₗᵐ(s) up to an
+# overall (l- and m-dependent) scaling factor. Also return the gradient of the
+# scalar output.
 function energy_and_gradient_for_classical_anisotropy(s::Vec3, stvexp::StevensExpansion)
     (; kmax, c2, c4, c6) = stvexp
 
diff --git a/src/System/System.jl b/src/System/System.jl
index 384e01d08..06f8ad777 100644
--- a/src/System/System.jl
+++ b/src/System/System.jl
@@ -138,7 +138,7 @@ case, it is convenient to construct the `Site` using [`position_to_site`](@ref),
 which always takes a position in fractional coordinates of the original lattice
 vectors.
 """
-const Site = NTuple{4, Int}
+const Site = Union{NTuple{4, Int}, CartesianIndex{4}}
 
 @inline to_cartesian(i::CartesianIndex{N}) where N = i
 @inline to_cartesian(i::NTuple{N, Int})    where N = CartesianIndex(i)
@@ -190,7 +190,7 @@ magnetic_moment(sys::System, site) = sys.units.μB * sys.gs[site] * sys.dipoles[
 volume(sys::System) = cell_volume(sys.crystal) * prod(sys.latsize)
 
 # The original crystal for a system, invariant under reshaping
-orig_crystal(sys) = isnothing(sys.origin) ? sys.crystal : sys.origin.crystal
+orig_crystal(sys) = something(sys.origin, sys).crystal
 
 # Position of a site in fractional coordinates of the original crystal
 function position(sys::System, site)
@@ -226,26 +226,8 @@ function position_to_site(sys::System, r)
 end
 
 
-# Given two [`Site`](@ref)s for a possibly reshaped system, return the
-# corresponding [`Bond`](@ref) for the original (unreshaped) crystal. This
-# `bond` can be used for symmetry analysis, e.g., as input to
-# [`print_bond`](@ref). See also [`site_to_atom`](@ref).
-#
-# Warning: This function will not be useful until site2 is wrapped properly.
-function sites_to_bond(sys::System{N}, site1, site2) where N
-    site1, site2 = to_cartesian.((site1, site2))
-
-    # Position of sites in multiples of original latvecs.
-    r1 = position(sys, site1)
-    r2 = position(sys, site2)
-
-    # Map to Bond for original crystal
-    return Bond(orig_crystal(sys), BondPos(r1, r2))
-end
-
 # Given a [`Site`](@ref)s for a possibly reshaped system, return the
-# corresponding atom index for the original (unreshaped) crystal. See also
-# [`sites_to_bond`](@ref).
+# corresponding atom index for the original (unreshaped) crystal.
 function site_to_atom(sys::System{N}, site) where N
     site = to_cartesian(site)
     r = position(sys, site)
@@ -437,3 +419,34 @@ function get_coherent_buffers(sys::System{N}, numrequested) where N
     end
     return view(sys.coherent_buffers, 1:numrequested)
 end
+
+"""
+    spin_operators(sys, i::Int)
+    spin_operators(sys, site::Int)
+
+Returns the three spin operators appropriate to an atom or [`Site`](@ref) index.
+Each is an ``N×N`` matrix of appropriate dimension ``N``.
+
+See also [`print_stevens_expansion`](@ref).
+"""
+spin_operators(sys::System{N}, i::Int) where N = spin_matrices(sys.Ns[i])
+spin_operators(sys::System{N}, site::Site) where N = spin_matrices(sys.Ns[to_atom(site)])
+
+"""
+    stevens_operators(sys, i::Int)
+    stevens_operators(sys, site::Int)
+
+Returns a generator of Stevens operators appropriate to an atom or
+[`Site`](@ref) index. The return value `O` can be indexed as `O[q,m]`, where ``0
+≤ q ≤ 6`` labels an irrep and ``m = -q, …, q``. This will produce an ``N×N``
+matrix of appropriate dimension ``N``.
+"""
+stevens_operators(sys::System{N}, i::Int) where N = StevensMatrices(sys.Ns[i])
+stevens_operators(sys::System{N}, site::Site) where N = StevensMatrices(sys.Ns[to_atom(site)])
+
+
+function spin_operators(sys::System{N}, b::Bond) where N
+    Si = spin_matrices(sys.Ns[b.i])
+    Sj = spin_matrices(sys.Ns[b.j])
+    return local_quantum_operators(Si, Sj)
+end
diff --git a/test/Project.toml b/test/Project.toml
index c835bed07..1b41f9cd3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 Ewalder = "2efd4204-53b6-4d30-8cf4-11a055f037eb"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 IOCapture = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/test/shared.jl b/test/shared.jl
index a6256d75f..af051776d 100644
--- a/test/shared.jl
+++ b/test/shared.jl
@@ -11,8 +11,9 @@ using Random, LinearAlgebra, IOCapture
 function add_linear_interactions!(sys, mode)
     set_external_field!(sys, (0.0, 1.0, 1.0))
     if mode == :SUN
-        # In SUN mode, anisotropy scales as ⟨Λ⟩ → κ ⟨Λ⟩.
-        set_anisotropy!(sys, 0.2*(𝒮[1]^4+𝒮[2]^4+𝒮[3]^4), 1)
+        # Kets scale as z → √κ z, so ⟨Λ⟩ → κ ⟨Λ⟩ is linear in κ
+        S = spin_operators(sys, 1)
+        set_onsite!(sys, 0.2*(S[1]^4+S[2]^4+S[3]^4), 1)
     end
 end
 
@@ -35,8 +36,9 @@ end
 
 function add_quartic_interactions!(sys, mode)
     if mode ∈ (:dipole, :large_S)
-        # In dipole mode, spins scale individually, S⁴ → κ⁴ S⁴
-        set_anisotropy!(sys, 0.2*(𝒮[1]^4+𝒮[2]^4+𝒮[3]^4), 1)
+        # Dipoles scale as ⟨S⟩ → κ ⟨S⟩, so ⟨S⟩⁴ → κ⁴ ⟨S⟩⁴ is quartic
+        S = spin_operators(sys, 1)
+        set_onsite!(sys, 0.2*(S[1]^4+S[2]^4+S[3]^4), 1)
     end
     if mode == :large_S
         set_biquadratic!(sys, 0.2, Bond(1, 3, [0, 0, 0]))
diff --git a/test/test_energy_consistency.jl b/test/test_energy_consistency.jl
index 770b9fff8..4e7286212 100644
--- a/test/test_energy_consistency.jl
+++ b/test/test_energy_consistency.jl
@@ -30,7 +30,9 @@
             if mode != :SUN
                set_biquadratic_at!(sys2, 0.7, (3,2,1,2), (3,1,1,3); offset=(0,-1,0))
             end
-            set_anisotropy_at!(sys2, 0.4*(𝒮[1]^4+𝒮[2]^4+𝒮[3]^4), (2,2,2,4))
+
+            S = spin_operators(sys2, 4)
+            set_anisotropy_at!(sys2, 0.4*(S[1]^4+S[2]^4+S[3]^4), (2,2,2,4))
             return sys2
         end
     end
diff --git a/test/test_lswt.jl b/test/test_lswt.jl
index 43a96ab97..eafe04bf3 100644
--- a/test/test_lswt.jl
+++ b/test/test_lswt.jl
@@ -28,7 +28,8 @@ end
 
     set_exchange!(sys, J,  Bond(1, 1, [1, 0, 0]))
     set_exchange!(sys, J′, Bond(1, 1, [0, 0, 1]))
-    set_anisotropy!(sys, D * 𝒮[3]^2, 1)
+    S = spin_operators(sys, 1)
+    set_onsite!(sys, D * S[3]^2, 1)
 
     Δt  = abs(0.05 / D)
     λ = 0.1
@@ -96,8 +97,9 @@ end
     D = D_ST / cov_factor
 
     set_exchange!(sys, J₁, Bond(1, 2, [0, 0, 0]))
-    Λ = D * (𝒮[1]^4 + 𝒮[2]^4 + 𝒮[3]^4)
-    set_anisotropy!(sys, Λ, 1)
+    S = spin_operators(sys, 1)
+    Λ = D * (S[1]^4 + S[2]^4 + S[3]^4)
+    set_onsite!(sys, Λ, 1)
 
     polarize_spin!(sys, (1, 1, 1), position_to_site(sys, (0, 0, 0)))
     polarize_spin!(sys, (1, -1, -1), position_to_site(sys, (1/2, 1/2, 0)))
diff --git a/test/test_operators.jl b/test/test_operators.jl
index d391465ab..7b3090d23 100644
--- a/test/test_operators.jl
+++ b/test/test_operators.jl
@@ -130,17 +130,17 @@ end
     # Check mapping between spherical tensors and Stevens operators
     for N in 2:7
         for k in 1:N-1
-            𝒪 = Sunny.stevens_matrices(k; N)
+            O = Sunny.stevens_matrices(k; N)
             T = spherical_tensors(k; N)
 
             # Check that Stevens operators are proper linear combination of
             # spherical tensors
-            @test 𝒪 ≈ Sunny.stevens_α[k] * T
+            @test O ≈ Sunny.stevens_α[k] * T
     
             # Check conversion of coefficients
             c = randn(2k+1)
             b = Sunny.transform_spherical_to_stevens_coefficients(k, c)
-            @test transpose(c)*T ≈ transpose(b)*𝒪
+            @test transpose(c)*T ≈ transpose(b)*O
         end
     end
 
@@ -167,7 +167,7 @@ end
 
     rng = Random.Xoshiro(0)
     R = Sunny.Mat3(Sunny.random_orthogonal(rng, 3; special=true))
-    N = 5
+    N = 7
 
     # Test axis-angle decomposition
     let
@@ -206,62 +206,60 @@ end
         @test c′1 ≈ c′2
     end
 
-    # Test that a symbolic operators rotate properly
-    let
-        p = randn(3)'*𝒮 + randn(5)'*Sunny.stevens_operator_symbols[2]
-        @test Sunny.operator_to_matrix(rotate_operator(p, R); N) ≈ rotate_operator(Sunny.operator_to_matrix(p; N), R)
-    end        
+    # Test evaluation of the classical Stevens functions (i.e. spherical
+    # harmonics) and their gradients
+    let 
+        using LinearAlgebra, FiniteDifferences
+
+        # Random dipole and Stevens coefficients
+        s = normalize(randn(Sunny.Vec3))
+        c = [randn(5), randn(9), randn(13)]
+        stvexp = Sunny.StevensExpansion(c[1], c[2], c[3])
+
+        # Rotate dipole and Stevens coefficients
+        s′ = R*s
+        c′ = Sunny.rotate_stevens_coefficients.(c, Ref(R))
+        stvexp′ = Sunny.StevensExpansion(c′[1], c′[2], c′[3])
+
+        # Verify that the energy is the same regardless of which is rotated
+        E1, _ = Sunny.energy_and_gradient_for_classical_anisotropy(s′, stvexp)
+        E2, _ = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp′)
+        @test E1 ≈ E2
+
+        # Verify that gradient agrees with finite differences
+        _, gradE1 = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp)
+        f(s) = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp)[1]
+        gradE2 = grad(central_fdm(5, 1), f, s)[1]
+
+        # When calculating gradE2, the value X = |S|^2 is treated as varying
+        # with S, such that dX/dS = 2S. Conversely, when calculating gradE1, the
+        # value X is treated as a constant, such that dX/dS = 0. In practice,
+        # gradE will be used to drive spin dynamics, for which |S| is constant,
+        # and the component of gradE parallel to S will be projected out anyway.
+        # Therefore we only need agreement in the components perpendicular to S.
+        gradE1 -= (gradE1⋅s)*s
+        gradE2 -= (gradE2⋅s)*s
+        @test gradE1 ≈ gradE2
+    end
 end
 
 
-@testitem "Symbolic operators" begin
-    using LinearAlgebra
+@testitem "Symbolics" begin
+    import IOCapture, DynamicPolynomials
 
-    # Internal conversion between spin and Stevens operators
-    let
-        J = randn(3,3)
-        J = J+J'
-        p = randn(3)'*𝒮 + 𝒮'*J*𝒮 +
-            randn(11)' * Sunny.stevens_operator_symbols[5] +
-            randn(13)' * Sunny.stevens_operator_symbols[6]
-        cp1 = p |> Sunny.operator_to_classical_polynomial
-        cp2 = p |> Sunny.operator_to_classical_stevens |> Sunny.operator_to_classical_polynomial
-        @test cp1 ≈ cp2
+    capt = IOCapture.capture() do
+        print_classical_spin_polynomial((1/4)𝒪[4,4] + (1/20)𝒪[4,0] + (3/5)*(𝒮'*𝒮)^2)
     end
+    @test capt.output == "𝒮₁⁴ + 𝒮₂⁴ + 𝒮₃⁴\n"
 
+    capt = IOCapture.capture() do
+        print_classical_stevens_expansion(𝒮[1]^4 + 𝒮[2]^4 + 𝒮[3]^4)
+    end
+    @test capt.output == "(1/20)𝒪₄₀ + (1/4)𝒪₄₄ + (3/5)X²\n"
 
-    # Test fast evaluation of Stevens functions
-    let
-        import DynamicPolynomials as DP
-
-        s = randn(Sunny.Vec3)
-        p = randn(5)' * Sunny.stevens_operator_symbols[2] + 
-            randn(9)' * Sunny.stevens_operator_symbols[4] +
-            randn(13)' * Sunny.stevens_operator_symbols[6]
-        (_, c2, _, c4, _, c6) = Sunny.operator_to_classical_stevens_coefficients(p, 1.0)
-
-        p_classical = Sunny.operator_to_classical_polynomial(p)
-        grad_p_classical = DP.differentiate(p_classical, Sunny.spin_classical_symbols)
-
-        E_ref = p_classical(Sunny.spin_classical_symbols => s)
-        gradE_ref = [g(Sunny.spin_classical_symbols => s) for g = grad_p_classical]
-
-        stvexp = Sunny.StevensExpansion(c2, c4, c6)
-        E, gradE = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp)
-
-        @test E ≈ E_ref
-
-        # Above, when calculating gradE_ref, the value X = |S|^2 is treated
-        # as varying with S, such that dX/dS = 2S. Conversely, when calculating
-        # gradE, the value X is treated as a constant, such that dX/dS = 0. In
-        # practice, gradE will be used to drive spin dynamics, for which |S| is
-        # constant, and the component of gradE parallel to S will be projected
-        # out anyway. Therefore we only need agreement between the parts of
-        # gradE and gradE_ref that are perpendicular to S.
-        gradE_ref -= (gradE_ref⋅s)*s / (s⋅s) # Orthogonalize to s
-        gradE -= (gradE⋅s)*s / (s⋅s)         # Orthogonalize to s
-        @test gradE_ref ≈ gradE
+    capt = IOCapture.capture() do
+        S = spin_matrices(5)
+        print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
     end
+    @test capt.output == "(1/20)𝒪₄₀ + (1/4)𝒪₄₄ + 102/5\n"
 end
-
-
diff --git a/test/test_resize.jl b/test/test_resize.jl
index 231a398e7..0adb479ef 100644
--- a/test/test_resize.jl
+++ b/test/test_resize.jl
@@ -185,8 +185,9 @@ end
     set_exchange!(sys, diagm([J′1pm, J′1pm, J′1zz]), Bond(1,1,[1,0,1]))
     set_exchange!(sys, diagm([J′2apm, J′2apm, J′2azz]), Bond(1,1,[1,2,1]))
 
+    S = spin_operators(sys, 1)
     D = 2.165
-    set_anisotropy!(sys, -D*𝒮[3]^2, 1)
+    set_onsite!(sys, -D*S[3]^2, 1)
 
     # periodic ground state for FeI2
     s = Sunny.Vec3(1, 0, 0)
diff --git a/test/test_samplers.jl b/test/test_samplers.jl
index 35f03fe28..0f73d1e7b 100644
--- a/test/test_samplers.jl
+++ b/test/test_samplers.jl
@@ -25,11 +25,10 @@
 
 
     function su3_anisotropy_model(; L=20, D=1.0, seed)
-        Λ = D*𝒮[3]^2
         cryst = asymmetric_crystal()
-
         sys = System(cryst, (L,1,1), [SpinInfo(1, S=1)], :SUN; seed)
-        set_anisotropy!(sys, Λ, 1)
+        S = spin_operators(sys, 1)
+        set_onsite!(sys, D*S[3]^2, 1)
         randomize_spins!(sys)
 
         return sys
@@ -40,9 +39,10 @@
         sys = System(cryst, (L,1,1), [SpinInfo(1, S=2)], :SUN; seed)
         randomize_spins!(sys)
 
+        S = spin_operators(sys, 1)
         R = Sunny.random_orthogonal(sys.rng, 3; special=true)
-        Λ = Sunny.rotate_operator(D*(𝒮[3]^2-(1/5)*𝒮[3]^4), R)
-        set_anisotropy!(sys, Λ, 1)
+        Λ = Sunny.rotate_operator(D*(S[3]^2-(1/5)*S[3]^4), R)
+        set_onsite!(sys, Λ, 1)
 
         return sys
     end
diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl
index a097305f3..7241ddb95 100644
--- a/test/test_symmetry.jl
+++ b/test/test_symmetry.jl
@@ -154,21 +154,22 @@ end
         k = 6
         i = 1
         cryst = Sunny.diamond_crystal()
+        O = Sunny.StevensMatrices(N)
 
         # print_site(cryst, i)
-        Λ = 𝒪[6,0]-21𝒪[6,4]
+        Λ = O[6,0]-21O[6,4]
         @test Sunny.is_anisotropy_valid(cryst, i, Λ)
 
         R = [normalize([1 1 -2]); normalize([-1 1 0]); normalize([1 1 1])]
         # print_site(cryst, i; R)
-        Λ = 𝒪[6,0]-(35/√8)*𝒪[6,3]+(77/8)*𝒪[6,6]
+        Λ = O[6,0]-(35/√8)*O[6,3]+(77/8)*O[6,6]
         Λ′ = rotate_operator(Λ, R)
         @test Sunny.is_anisotropy_valid(cryst, i, Λ′)
 
         latvecs = lattice_vectors(1.0, 1.1, 1.0, 90, 90, 90)
         cryst = Crystal(latvecs, [[0., 0., 0.]])
         # print_site(cryst, i)
-        Λ = randn()*(𝒪[6,0]-21𝒪[6,4]) + randn()*(𝒪[6,2]+(16/5)*𝒪[6,4]+(11/5)*𝒪[6,6])
+        Λ = randn()*(O[6,0]-21O[6,4]) + randn()*(O[6,2]+(16/5)*O[6,4]+(11/5)*O[6,6])
         @test Sunny.is_anisotropy_valid(cryst, i, Λ)
     end
 
@@ -180,38 +181,36 @@ end
         positions = [[0.1, 0, 0], [0, 0.1, 0], [-0.1, 0, 0], [0, -0.1, 0]]
         cryst = Crystal(latvecs, positions)
 
-        # Most general allowed anisotropy for this crystal
-        Λ = randn(9)'*[𝒪[2,0], 𝒪[2,2], 𝒪[4,0], 𝒪[4,2], 𝒪[4,4], 𝒪[6,0], 𝒪[6,2], 𝒪[6,4], 𝒪[6,6]]
-
-        # Test anisotropy invariance in dipole mode
-        sys = System(cryst, (1,1,1), [SpinInfo(1, S=1)], :dipole)
-        set_anisotropy!(sys, Λ, 1)
-        randomize_spins!(sys)
-        E1 = energy(sys)
-        # By circularly shifting the atom index, we are effectively rotating
-        # site positions by π/2 clockwise (see comment above `positions`)
-        sys.dipoles .= circshift(sys.dipoles, (0,0,0,1))
-        # Rotate spin vectors correspondingly
-        R = Sunny.Mat3([0 1 0; -1 0 0; 0 0 1])
-        sys.dipoles .= [R*d for d in sys.dipoles]
-        E2 = energy(sys)
-        @test E1 ≈ E2
-
-        # Test anisotropy invariance in SU(N) mode
-        sys = System(cryst, (1,1,1), [SpinInfo(1, S=2)], :SUN)
-        N = only(unique(sys.Ns))
-        set_anisotropy!(sys, Λ, 1)
-        randomize_spins!(sys)
-        E1 = energy(sys)
-        # By circularly shifting the atom index, we are effectively rotating
-        # site positions by π/2 clockwise (see comment above `positions`)
-        sys.coherents .= circshift(sys.coherents, (0,0,0,1))
-        # Rotate kets correspondingly
-        U = Sunny.unitary_for_rotation(R; N)
-        sys.coherents .= [U*z for z in sys.coherents]
-        sys.dipoles .= Sunny.expected_spin.(sys.coherents)
-        E2 = energy(sys)
-        @test E1 ≈ E2
+        for mode in (:dipole, :SUN)
+            sys = System(cryst, (1,1,1), [SpinInfo(1, S=2)], mode)
+            randomize_spins!(sys)
+
+            # Most general allowed anisotropy for this crystal
+            c = randn(9)
+            O = stevens_operators(sys, 1)
+            Λ = sum(c .* [O[2,0], O[2,2], O[4,0], O[4,2], O[4,4], O[6,0], O[6,2], O[6,4], O[6,6]])
+            set_onsite!(sys, Λ, 1)
+
+            E1 = energy(sys)
+
+            # By circularly shifting the atom index, we are effectively rotating
+            # site positions by π/2 clockwise (see comment above `positions`).
+            # Rotate spin state correspondingly
+            R = Sunny.Mat3([0 1 0; -1 0 0; 0 0 1])
+            sys.dipoles .= circshift(sys.dipoles, (0,0,0,1))
+            sys.dipoles .= [R*d for d in sys.dipoles]
+
+            # If coherents are present, perform same operation
+            if mode == :SUN
+                U = Sunny.unitary_for_rotation(R; N=sys.Ns[1])
+                sys.coherents .= circshift(sys.coherents, (0,0,0,1))
+                sys.coherents .= [U*z for z in sys.coherents]
+            end
+
+            # Verify energy is invariant
+            E2 = energy(sys)
+            @test E1 ≈ E2
+        end
     end
 end
 
@@ -220,15 +219,17 @@ end
     latvecs = lattice_vectors(1.0, 1.1, 1.0, 90, 90, 90)
     cryst = Crystal(latvecs, [[0., 0., 0.]])
     
-    i = 1
-    Λ = randn()*(𝒪[2,0]+3𝒪[2,2]) +
-        randn()*(𝒪[4,0]-5𝒪[4,2]) + randn()*(𝒪[4,0]+5𝒪[4,4]) +
-        randn()*(𝒪[6,0]-21𝒪[6,4]) + randn()*(𝒪[6,0]+(105/16)𝒪[6,2]+(231/16)𝒪[6,6])
     
     # Dipole system with renormalized anisotropy
     sys0 = System(cryst, (1,1,1), [SpinInfo(1, S=3)], :dipole)
     randomize_spins!(sys0)
-    set_anisotropy!(sys0, Λ, i)
+
+    i = 1
+    O = stevens_operators(sys0, i)
+    Λ = randn()*(O[2,0]+3O[2,2]) +
+        randn()*(O[4,0]-5O[4,2]) + randn()*(O[4,0]+5O[4,4]) +
+        randn()*(O[6,0]-21O[6,4]) + randn()*(O[6,0]+(105/16)O[6,2]+(231/16)O[6,6])
+    set_onsite!(sys0, Λ, i)
     E0 = energy(sys0)
     
     # Corresponding SU(N) system
@@ -236,7 +237,7 @@ end
     for site in all_sites(sys)
         polarize_spin!(sys, sys0.dipoles[site], site)
     end
-    set_anisotropy!(sys, Λ, i)
+    set_onsite!(sys, Λ, i)
     E = energy(sys)
     
     @test E ≈ E0