diff --git a/test/test_binning.jl b/test/test_binning.jl
index 570aebd66..cb5a4d835 100644
--- a/test/test_binning.jl
+++ b/test/test_binning.jl
@@ -61,14 +61,14 @@
     # TODO: Test broadening
     formula = intensity_formula(sc, :perp)
     is, counts = intensities_binned(sc, params, formula)
-    is_golden = [0.2452071781; 0.0864959953; 0.1158561543; 0.0299947084;;;; 0.1226035891; 0.0432479977; 0.0579280772; 0.0149973542;;;; -0.0; -0.0; 0.0; -0.0]
+    is_golden = size(sc.samplebuf,6) * [0.2452071781; 0.0864959953; 0.1158561543; 0.0299947084;;;; 0.1226035891; 0.0432479977; 0.0579280772; 0.0149973542;;;; -0.0; -0.0; 0.0; -0.0]
     @test is ≈ is_golden
     @test all(counts .== 1.)
 
 
     # SQTODO: need to implement correct bin widths in powder_average_binned!
     is, counts = powder_average_binned(sc, (0,6π,6π/4), formula)
-    is_golden = [8.9511865548 4.4755932774 -0.0; 35.9054210445 17.9527105222 -0.0; 102.2777600371 51.1388800185 -0.0; 91.4466208136 45.7233104068 -0.0]
+    is_golden = size(sc.samplebuf,6) * [8.9511865548 4.4755932774 -0.0; 35.9054210445 17.9527105222 -0.0; 102.2777600371 51.1388800185 -0.0; 91.4466208136 45.7233104068 -0.0]
     counts_golden = [3.0 3.0 3.0; 15.0 15.0 15.0; 28.0 28.0 28.0; 39.0 39.0 39.0]
     @test is ≈ is_golden
     @test counts == counts_golden
@@ -79,12 +79,12 @@
     golden_correlations = [(sz,sz),(sy,sz),(sx,sz),(sy,sy),(sx,sy),(sx,sx)]
     formula = intensity_formula((k,ω,c) -> c,sc,golden_correlations; kT = 4.7, formfactors = [FormFactor("Fe2")], return_type = StaticArrays.SVector{6,ComplexF64})
     is, counts = intensities_binned(sc, params, formula)
-    is_golden = [0.0180419588 -0.0; -0.0 -0.0; -0.0225463643 0.0; 0.0 0.0; -0.0384472716 0.0; 0.0 -0.0; 0.0281753523 -0.0; -0.0 0.0; 0.0480461242 0.0; -0.0 -0.0; 0.0819308319 -0.0; -0.0 0.0; 0.0235568126 -0.0; 0.0 0.0; -0.00045821 0.0; -0.0150200679 -0.0; -0.0069056941 -0.0; 0.0143894403 -0.0; 0.0095858637 -0.0; -0.0 -0.0; -0.0090405317 -0.0; -0.0046830351 -0.0; 0.0108140522 -0.0; -0.0 0.0]
+    is_golden = size(sc.samplebuf,6) * [0.0180419588 -0.0; -0.0 -0.0; -0.0225463643 0.0; 0.0 0.0; -0.0384472716 0.0; 0.0 -0.0; 0.0281753523 -0.0; -0.0 0.0; 0.0480461242 0.0; -0.0 -0.0; 0.0819308319 -0.0; -0.0 0.0; 0.0235568126 -0.0; 0.0 0.0; -0.00045821 0.0; -0.0150200679 -0.0; -0.0069056941 -0.0; 0.0143894403 -0.0; 0.0095858637 -0.0; -0.0 -0.0; -0.0090405317 -0.0; -0.0046830351 -0.0; 0.0108140522 -0.0; -0.0 0.0]
     @test reinterpret(Float64, is[1:2, 1, 1, 2:3]) ≈ is_golden
     @test all(counts .== 1.)
 
     is, counts = intensities_binned(sc, params, formula; integrated_kernel=integrated_lorentzian(fwhm=1.0))
-    is_golden = [0.0112630945 0.005987424; -0.0 -0.0; -0.0140750699 -0.0074822609; 0.0 0.0; -0.0240015653 -0.0127591532; 0.0 0.0; 0.0175890909 0.0093503029; -0.0 -0.0; 0.0299938627 0.0159446388; -0.0 -0.0; 0.0511471459 0.0271896546; -0.0 -0.0; 0.0147058647 0.0078175893; 0.0 0.0; -0.0002860478 -0.0001520621; -0.0093766118 -0.004984576; -0.0043110333 -0.0022917311; 0.0089829284 0.0047752952; 0.0059841889 0.0031811751; -0.0 -0.0; -0.0056437532 -0.0030002006; -0.0029234889 -0.0015541171; 0.0067509129 0.0035887631; -0.0 -0.0]
+    is_golden = size(sc.samplebuf,6) * [0.0112630945 0.005987424; -0.0 -0.0; -0.0140750699 -0.0074822609; 0.0 0.0; -0.0240015653 -0.0127591532; 0.0 0.0; 0.0175890909 0.0093503029; -0.0 -0.0; 0.0299938627 0.0159446388; -0.0 -0.0; 0.0511471459 0.0271896546; -0.0 -0.0; 0.0147058647 0.0078175893; 0.0 0.0; -0.0002860478 -0.0001520621; -0.0093766118 -0.004984576; -0.0043110333 -0.0022917311; 0.0089829284 0.0047752952; 0.0059841889 0.0031811751; -0.0 -0.0; -0.0056437532 -0.0030002006; -0.0029234889 -0.0015541171; 0.0067509129 0.0035887631; -0.0 -0.0]
     counts_golden = [0.7164129516; 0.7164129516; 0.7164129516; 0.7164129516;;;; 0.6814242206; 0.6814242206; 0.6814242206; 0.6814242206;;;; 0.5437318669; 0.5437318669; 0.5437318669; 0.5437318669]
     @test reinterpret(Float64, is[1:2, 1, 1, 2:3]) ≈ is_golden
     @test counts ≈ counts_golden
@@ -92,7 +92,7 @@
     # Test all components using :full
     formula = intensity_formula(sc, :full; kT = 4.7, formfactors = [FormFactor("Fe2")])
     is, counts = intensities_binned(sc, params, formula)
-    is2_golden = [0.0206923266 + 0.0im -0.0172987545 - 0.0089608308im -0.0132138143 + 0.0275337118im; -0.0172987545 + 0.0089608308im 0.0183422291 + 0.0im -0.0008767696 - 0.0287403967im; -0.0132138143 - 0.0275337118im -0.0008767696 + 0.0287403967im 0.0450751717 + 0.0im]
+    is2_golden = size(sc.samplebuf,6) * [0.0206923266 + 0.0im -0.0172987545 - 0.0089608308im -0.0132138143 + 0.0275337118im; -0.0172987545 + 0.0089608308im 0.0183422291 + 0.0im -0.0008767696 - 0.0287403967im; -0.0132138143 - 0.0275337118im -0.0008767696 + 0.0287403967im 0.0450751717 + 0.0im]
     @test is[2] ≈ is2_golden
     @test all(counts .== 1.)
 
diff --git a/test/test_intensities_setup.jl b/test/test_intensities_setup.jl
index 7206db3fa..a2c1e62cb 100644
--- a/test/test_intensities_setup.jl
+++ b/test/test_intensities_setup.jl
@@ -95,11 +95,11 @@ end
     sc = dynamical_correlations(sys; dt=1, nω=3, ωmax=1)
     add_sample!(sc, sys)
 
-    # A simple sum (rather than energy-integral) is possible because a `Δω`
-    # factor is already included in the intensity data. This convention will
-    # change per https://github.com/SunnySuite/Sunny.jl/issues/264
     sum_rule_ixs = Sunny.Trace(sc.observables).indices
     sub_lat_sum_rules = sum(sc.data[sum_rule_ixs,:,:,:,:,:,:],dims = [1,4,5,6,7])[1,:,:,1,1,1,1]
+    # Since sc.data is in units of S²/BZ/fs, each entry is already an estimate of S².
+    # To get the frequency-averaged estimate, we need to divide by the number of estimators:
+    sub_lat_sum_rules ./= prod(size(sc.data)[[4,5,6,7]])
 
     # SU(N) sum rule for S = 1/2:
     # ⟨∑ᵢSᵢ²⟩ = 3/4 on every site, but because we're classical, we
@@ -108,7 +108,7 @@ end
     # need to include the g factor. This is the equal-space-and-time correlation value:
     gS_squared = (2 * 1/2)^2
 
-    expected_sum = gS_squared * prod(sys.latsize)
+    expected_sum = gS_squared
     # This sum rule should hold for each sublattice, independently, and only
     # need to be taken over a single BZ (which is what sc.data contains) to hold:
     @test [sub_lat_sum_rules[i,i] for i = 1:Sunny.natoms(sc.crystal)] ≈ expected_sum * ones(ComplexF64,Sunny.natoms(sc.crystal))
@@ -132,5 +132,5 @@ end
     is, counts = intensities_binned(sc,params_pasr,formula)
     expected_multi_BZ_sum = gS_squared * prod(nbzs) * nfs # ⟨S⋅S⟩
     expected_multi_BZ_sum_times_natoms = expected_multi_BZ_sum * Sunny.natoms(sc.crystal) # Nₐ×⟨S⋅S⟩
-    @test sum(is ./ counts) ≈ expected_multi_BZ_sum_times_natoms / size(sc.data,7)
+    @test sum(is) ≈ expected_multi_BZ_sum_times_natoms
 end