diff --git a/docs/src/parallelism.md b/docs/src/parallelism.md
index 1ff1c53e6..a236b10fb 100644
--- a/docs/src/parallelism.md
+++ b/docs/src/parallelism.md
@@ -90,10 +90,9 @@ scs = [SampledCorrelations(sys; dt=0.1, energies, measure) for _ in 1:npar]
 
 !!! warning "Dealing with memory constraints"
 
-    If you have many threads available and are working with a large system, you
-    may not have enough memory to store all these systems and correlations. In
-    that case, simply reduce `npar` to a small enough value that you can make the
-    necessary allocations.
+  If many threads are available and the system sizes are large, there may not be
+  enough computer memory available to store the dynamical correlation data. In
+  this case, consider reducing the number of parallel threads, `npar`.
 
 When the `Threads.@threads` macro is applied before a `for` loop, the
 iterations of the loop will execute in parallel using the available threads.
diff --git a/docs/src/renormalization.md b/docs/src/renormalization.md
index 84558a39d..caa91f8ea 100644
--- a/docs/src/renormalization.md
+++ b/docs/src/renormalization.md
@@ -42,8 +42,9 @@ 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.
+The Stevens operators ``\mathcal{O}_{k, q}`` are [defined below](@ref
+"Definition of Stevens operators") as ``k``th order polynomials of the spin
+operators.
 
 ## Renormalization procedure for `:dipole` mode
 
@@ -51,13 +52,14 @@ 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]](https://arxiv.org/abs/2304.03874).
+will automatically renormalize the magnitude of each Stevens expectation value
+to achieve maximal consistency with `:SUN` mode. This procedure was derived in
+[D. Dahlbom et al., [arXiv:2304.03874]](https://arxiv.org/abs/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,
+site. It can be expanded in [Stevens operators](@ref "Definition of Stevens
+operators"),
 ```math
 \hat{\mathcal H}_{\mathrm{local}} = \sum_{k, q} A_{k,q} \hat{\mathcal{O}}_{k,q},
 ```
diff --git a/src/KPM/SpinWaveTheoryKPM.jl b/src/KPM/SpinWaveTheoryKPM.jl
index 8d8f9ec8f..20f54fb3e 100644
--- a/src/KPM/SpinWaveTheoryKPM.jl
+++ b/src/KPM/SpinWaveTheoryKPM.jl
@@ -18,9 +18,9 @@ Reasonable starting points are `1e-1` (more speed) or `1e-2` (more accuracy).
 
 !!! warning "Missing intensity at small quasi-particle energy"
 
-    The KPM calculation may mask intensities at small energies ``ω``. In
-    particular, such artifacts may arise near the Goldstone modes of an
-    ordered state with continuous symmetry.
+  The KPM-calculated intensities are unreliable at small energies ``ω``. In
+  particular, KPM may mask intensities that arise near the Goldstone modes of an
+  ordered state with continuous symmetry.
 
 References:
 
diff --git a/src/MonteCarlo/Samplers.jl b/src/MonteCarlo/Samplers.jl
index e4758ceab..f3c70f9fe 100644
--- a/src/MonteCarlo/Samplers.jl
+++ b/src/MonteCarlo/Samplers.jl
@@ -118,13 +118,15 @@ Multiple proposals can be mixed with the macro [`@mix_proposals`](@ref).
 The returned object stores fields `ΔE` and `ΔS`, which represent the cumulative
 change to the net energy and dipole, respectively.
 
-!!! warning "Efficiency considerations
-
-    A [`Langevin`](@ref) sampler is frequently much more efficient than a
-    `LocalSampler` for simulating Heisenberg-like spins that vary continuously. A
-    `LocalSampler` is appropriate in the special case that the spin states are
-    effectively discrete. E.g., [`propose_flip`](@ref) is very helpful simulating
-    Ising-like spins that arise due to a strong easy-axis anisotropy.
+!!! warning "Efficiency considerations"
+
+  Prefer [`Langevin`](@ref) sampling in most cases. Langevin dynamics will
+  usually be much more efficient for sampling Heisenberg-like spins that vary
+  continuously. `LocalSampler` is most useful for sampling from discrete spin
+  states. In particular, [`propose_flip`](@ref) may be required for sampling
+  Ising-like spins that arise due to a strong easy-axis anisotropy. For strong
+  but finite single-ion anisotropy, consider alternating between `Langevin` and
+  `LocalSampler` update steps.
 """
 mutable struct LocalSampler{F}
     kT      :: Float64   # Temperature
diff --git a/src/Operators/Stevens.jl b/src/Operators/Stevens.jl
index 5cb9e0514..fb0bb0c6b 100644
--- a/src/Operators/Stevens.jl
+++ b/src/Operators/Stevens.jl
@@ -225,16 +225,16 @@ end
 """
     stevens_matrices(s)
 
-Returns a generator of Stevens operators in the spin-`s` representation. The
-return value `O` can be indexed as `O[k,q]`, where ``0 ≤ k ≤ 6`` labels an irrep
-of SO(3) and ``-k ≤ q ≤ k``. This will produce an ``N×N`` matrix where ``N = 2s
-+ 1``. Linear combinations of Stevens operators can be used as a "physical
-basis" for decomposing local observables. To see this decomposition, use
+Returns the Stevens operators in the spin-`s` representation. The return value
+`O` can be indexed as `O[k,q]`, where ``0 ≤ k ≤ 6`` labels an irrep of SO(3) and
+``-k ≤ q ≤ k``. This will produce an ``N×N`` matrix where ``N = 2s + 1``. Linear
+combinations of Stevens operators can be used as a "physical basis" for
+decomposing local observables. To see this decomposition, use
 [`print_stevens_expansion`](@ref).
 
 If `s == Inf`, then symbolic operators will be returned. In this infinite
-dimensional limit, the Stevens operators become homogeneous polynomials of
-commuting spin operators.
+dimensional representation, the Stevens operators become homogeneous polynomials
+of commuting spin operators.
 
 # Example
 ```julia
diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl
index 1af12c793..1d6dc7c58 100644
--- a/src/System/Ewald.jl
+++ b/src/System/Ewald.jl
@@ -253,11 +253,11 @@ than `enable_dipole_dipole!`.
 
 !!! warning "Mutation of existing couplings"
 
-    This function will modify existing bilinear couplings between spins by
-    adding dipole-dipole interactions. It must therefore be called _after_
-    all other pair couplings have been specified. Conversely, any calls to
-    `set_exchange!`, `set_pair_coupling!`, etc. will irreversibly delete the
-    dipole-dipole interactions that have been introduced by this function.
+  This function will modify existing bilinear couplings between spins by adding
+  dipole-dipole interactions. It must therefore be called _after_ all other pair
+  couplings have been specified. Conversely, any calls to `set_exchange!`,
+  `set_pair_coupling!`, etc. will irreversibly delete the dipole-dipole
+  interactions that have been introduced by this function.
 """
 function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff, μ0_μB²=nothing) where N
     if isnothing(μ0_μB²)
diff --git a/src/System/OnsiteCoupling.jl b/src/System/OnsiteCoupling.jl
index fa2bcb542..473a4047b 100644
--- a/src/System/OnsiteCoupling.jl
+++ b/src/System/OnsiteCoupling.jl
@@ -114,6 +114,18 @@ set_onsite_coupling!(sys, S -> 20*(S[1]^4 + S[2]^4 + S[3]^4), i)
 O = stevens_matrices(spin_label(sys, i))
 set_onsite_coupling!(sys, O[4,0] + 5*O[4,4], i)
 ```
+
+!!! warning "Limitations arising from quantum spin operators"
+
+  Single-ion anisotropy is physically impossible for local moments with quantum
+  spin ``s = 1/2``. Consider, for example, that any Pauli matrix squared gives
+  the identity. More generally, one can verify that the ``k``th order Stevens
+  operators `O[k, q]` are zero whenever ``s < k/2``. Consequently, an anisotropy
+  quartic in the spin operators requires ``s ≥ 2`` and an anisotropy of sixth
+  order requires ``s ≥ 3``. To circumvent this physical limitation, Sunny
+  provides a mode `:dipole_uncorrected` that naïvely replaces quantum spin
+  operators with classical moments. See the documentation page [Interaction
+  Renormalization](@ref) for more information.
 """
 function set_onsite_coupling!(sys::System, op, i::Int)
     is_homogeneous(sys) || error("Use `set_onsite_coupling_at!` for an inhomogeneous system.")