Add some doc warnings

docs/src/renormalization.md

@@ -42,22 +42,24 @@ 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
 
 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},
 ```

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

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

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.")