diff --git a/docs/make.jl b/docs/make.jl
index 1cc4a92..e8e405c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,11 +1,5 @@
 using Pkg
 
-DOCUMENTER_VERSION =
-    [p for (uuid, p) in Pkg.dependencies() if p.name == "Documenter"][1].version
-if DOCUMENTER_VERSION <= v"1.3.0"
-    Pkg.develop("Documenter")
-end
-
 using Documenter
 using QuantumPropagators
 using QuantumPropagators: AbstractPropagator, set_t!, set_state!
@@ -35,10 +29,6 @@ links = InterLinks(
         "https://github.com/KristofferC/TimerOutputs.jl",
         joinpath(@__DIR__, "src", "inventories", "TimerOutputs.toml"),
     ),
-    # We'll use `@extref` for links from docstrings to sections so that the
-    # docstrings can also be rendered as part of the QuantumControl
-    # documentation.
-    "QuantumPropagators" => "https://juliaquantumcontrol.github.io/QuantumPropagators.jl/$DEV_OR_STABLE",
     "QuantumControlBase" => "https://juliaquantumcontrol.github.io/QuantumControlBase.jl/$DEV_OR_STABLE",
     "ComponentArrays" => (
         "https://jonniedie.github.io/ComponentArrays.jl/stable/",
@@ -89,7 +79,8 @@ makedocs(;
                 "https://juliaquantumcontrol.github.io/QuantumControl.jl/dev/assets/topbar/topbar.js"
             ),
         ],
-        footer="[$NAME.jl]($GITHUB) v$VERSION docs powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl)."
+        footer="[$NAME.jl]($GITHUB) v$VERSION docs powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl).",
+        size_threshold_ignore=["api/quantumpropagators.md",]
     ),
     pages=[
         "Home" => "index.md",
@@ -101,7 +92,7 @@ makedocs(;
         hide("Benchmarks" => "benchmarks.md", [joinpath("benchmarks", "profiling.md")]),
         "API" => "api/quantumpropagators.md",
         "References" => "references.md",
-    ]
+    ],
 )
 
 println("Finished makedocs")
diff --git a/ext/QuantumPropagatorsODEExt.jl b/ext/QuantumPropagatorsODEExt.jl
index 710a31a..0d1328d 100644
--- a/ext/QuantumPropagatorsODEExt.jl
+++ b/ext/QuantumPropagatorsODEExt.jl
@@ -249,6 +249,7 @@ function set_state!(propagator::ODEPropagator, state)
             ODE.set_u!(propagator.integrator, state)
         end
     end
+    return propagator.state
 end
 
 
diff --git a/src/cheby.jl b/src/cheby.jl
index 07e8150..d6f626a 100644
--- a/src/cheby.jl
+++ b/src/cheby.jl
@@ -154,7 +154,7 @@ function cheby!(Ψ, H, dt, wrk; kwargs...)
 
     Δ = wrk.Δ
     β::Float64 = (Δ / 2) + E_min  # "normfactor"
-    @assert abs(dt) ≈ abs(wrk.dt) "wrk was initialized for dt=$(wrk.dt), not dt=$dt"
+    @assert abs(dt) ≈ abs(wrk.dt) "wrk was initialized for dt=$(wrk.dt), not dt=abs($dt)"
     if dt > 0
         c = -2im / Δ
     else
@@ -228,7 +228,7 @@ function cheby(Ψ, H, dt, wrk; kwargs...)
 
     Δ = wrk.Δ
     β::Float64 = (Δ / 2) + E_min  # "normfactor"
-    @assert dt ≈ wrk.dt "wrk was initialized for dt=$(wrk.dt), not dt=$dt"
+    @assert abs(dt) ≈ wrk.dt "wrk was initialized for dt=$(wrk.dt), not dt=abs($dt)"
     if dt > 0
         c = -2im / Δ
     else
@@ -237,9 +237,6 @@ function cheby(Ψ, H, dt, wrk; kwargs...)
     a = wrk.coeffs
     ϵ = wrk.limit
     @assert length(a) > 1 "Need at least 2 Chebychev coefficients"
-    v0 = wrk.v0
-    v1 = wrk.v1
-    v2 = wrk.v2
 
     v0 = Ψ
     Ψ = a[1] * v0
@@ -256,16 +253,17 @@ function cheby(Ψ, H, dt, wrk; kwargs...)
         @timeit_debug wrk.timing_data "matrix-vector product" begin
             v2 = H * v1
         end
-        v2 += -v1 * β
-        v2 = c * v2
         if check_normalization
+            v2 = c * (v2 - v1 * β)
             map_norm = abs(dot(v1, v2)) / (2 * norm(v1)^2)
             @assert(
                 map_norm <= (1.0 + ϵ),
                 "Incorrect normalization (E_min=$(E_min), Δ=$(Δ))"
             )
+            v2 += v0
+        else
+            v2 = c * (v2 - β * v1) + v0
         end
-        v2 += v0
 
         Ψ += a[i] * v2
 
diff --git a/src/cheby_propagator.jl b/src/cheby_propagator.jl
index c60a4e7..159a707 100644
--- a/src/cheby_propagator.jl
+++ b/src/cheby_propagator.jl
@@ -245,7 +245,7 @@ function reinit_prop!(
     transform_control_ranges=_transform_control_ranges,
     _...
 )
-    set_state!(propagator, state)
+    state = set_state!(propagator, state)
 
     wrk = propagator.wrk
     need_to_recalculate_cheby_coeffs = false
@@ -355,8 +355,8 @@ function prop_step!(propagator::ChebyPropagator)
         end
         tlist = getfield(propagator, :tlist)
         (0 < n < length(tlist)) || return nothing
-        _pwc_set_genop!(propagator, n)
         if propagator.inplace
+            H = _pwc_set_genop!(propagator, n)
             Cheby.cheby!(
                 Ψ,
                 H,
@@ -365,6 +365,7 @@ function prop_step!(propagator::ChebyPropagator)
                 check_normalization=propagator.check_normalization
             )
         else
+            H = _pwc_get_genop(propagator, n)
             Ψ = Cheby.cheby(
                 Ψ,
                 H,
diff --git a/src/interfaces/generator.jl b/src/interfaces/generator.jl
index 1b9a7ea..0c7cfad 100644
--- a/src/interfaces/generator.jl
+++ b/src/interfaces/generator.jl
@@ -4,11 +4,13 @@ using ..Generators: Generator
 """Check the dynamical `generator` for propagating `state` over `tlist`.
 
 ```julia
-@test check_generator(generator; state, tlist,
-                     for_mutable_state=true, for_immutable_state=true,
-                     for_pwc=true, for_time_continuous=false,
-                     for_expval=true, for_parameterization=false,
-                     atol=1e-14, quiet=false)
+@test check_generator(
+    generator; state, tlist,
+    for_mutable_operator=true, for_immutable_operator=true,
+    for_mutable_state=true, for_immutable_state=true,
+    for_pwc=true, for_time_continuous=false,
+    for_expval=true, for_parameterization=false,
+    atol=1e-14, quiet=false)
 ```
 
 verifies the given `generator`:
@@ -27,14 +29,16 @@ If `for_pwc` (default):
 * [`evaluate(generator, tlist, n)`](@ref evaluate) must return a valid
   operator ([`check_operator`](@ref)), with forwarded keyword arguments
   (including `for_expval`)
-* [`evaluate!(op, generator, tlist, n)`](@ref evaluate!) must be defined
+* If `for_mutable_operator`,
+  [`evaluate!(op, generator, tlist, n)`](@ref evaluate!) must be defined
 
 If `for_time_continuous`:
 
 * [`evaluate(generator, t)`](@ref evaluate) must return a valid
   operator ([`check_operator`](@ref)), with forwarded keyword arguments
   (including `for_expval`)
-* [`evaluate!(op, generator, t)`](@ref evaluate!) must be defined
+* If `for_mutable_operator`, [`evaluate!(op, generator, t)`](@ref evaluate!)
+  must be defined
 
 If `for_parameterization` (may require the `RecursiveArrayTools` package to be
 loaded):
@@ -51,6 +55,8 @@ function check_generator(
     generator;
     state,
     tlist,
+    for_mutable_operator=true,
+    for_immutable_operator=true,
     for_mutable_state=true,
     for_immutable_state=true,
     for_expval=true,
@@ -208,26 +214,27 @@ function check_generator(
             success = false
         end
 
-        try
-            op = evaluate(generator, tlist, 1)
-            evaluate!(op, generator, tlist, length(tlist) - 1)
-        catch exc
-            quiet || @error(
-                "$(px)`evaluate!(op, generator, tlist, n)` must be defined.",
-                exception = (exc, catch_abbreviated_backtrace())
-            )
-            success = false
-        end
-
-        try
-            op = evaluate(generator, tlist, 1)
-            evaluate!(op, generator, tlist, length(tlist) - 1; vals_dict)
-        catch exc
-            quiet || @error(
-                "$(px)`evaluate!(op, generator, tlist, n; vals_dict)` must be defined.",
-                exception = (exc, catch_abbreviated_backtrace())
-            )
-            success = false
+        if for_mutable_operator
+            try
+                op = evaluate(generator, tlist, 1)
+                evaluate!(op, generator, tlist, length(tlist) - 1)
+            catch exc
+                quiet || @error(
+                    "$(px)`evaluate!(op, generator, tlist, n)` must be defined.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
+            try
+                op = evaluate(generator, tlist, 1)
+                evaluate!(op, generator, tlist, length(tlist) - 1; vals_dict)
+            catch exc
+                quiet || @error(
+                    "$(px)`evaluate!(op, generator, tlist, n; vals_dict)` must be defined.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
         end
 
     end
@@ -284,28 +291,28 @@ function check_generator(
             success = false
         end
 
-        try
-            op = evaluate(generator, tlist[begin])
-            evaluate!(op, generator, tlist[end])
-        catch exc
-            quiet || @error(
-                "$(px)`evaluate!(op, generator, t)` must be defined.",
-                exception = (exc, catch_abbreviated_backtrace())
-            )
-            success = false
-        end
-
-        try
-            op = evaluate(generator, tlist[begin])
-            evaluate!(op, generator, tlist[end]; vals_dict)
-        catch exc
-            quiet || @error(
-                "$(px)`evaluate!(op, generator, t; vals_dict)` must be defined.",
-                exception = (exc, catch_abbreviated_backtrace())
-            )
-            success = false
+        if for_mutable_operator
+            try
+                op = evaluate(generator, tlist[begin])
+                evaluate!(op, generator, tlist[end])
+            catch exc
+                quiet || @error(
+                    "$(px)`evaluate!(op, generator, t)` must be defined.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
+            try
+                op = evaluate(generator, tlist[begin])
+                evaluate!(op, generator, tlist[end]; vals_dict)
+            catch exc
+                quiet || @error(
+                    "$(px)`evaluate!(op, generator, t; vals_dict)` must be defined.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
         end
-
     end
 
 
diff --git a/src/interfaces/propagator.jl b/src/interfaces/propagator.jl
index 3f7b9b2..a643242 100644
--- a/src/interfaces/propagator.jl
+++ b/src/interfaces/propagator.jl
@@ -39,6 +39,7 @@ initialized with [`init_prop`](@ref).
   `propagator.state`.
 * [`set_state!(propagator, state)`](@ref set_state!) for an in-place propagator
     must overwrite `propagator.state` in-place.
+* [`set_state!`](@ref) must return the set `propagator.state`
 * In a [`PiecewisePropagator`](@ref), `propagator.parameters` must be a dict
   mapping controls to a vector of values, one for each interval on
   `propagator.tlist`
@@ -236,7 +237,7 @@ function check_propagator(
     end
 
     try
-        set_state!(propagator, Ψ₀)
+        Ψ = set_state!(propagator, Ψ₀)
         if norm(propagator.state - Ψ₀) > atol
             quiet ||
                 @error "$(px)`set_state!(propagator, state)` must set `propagator.state`"
@@ -249,6 +250,11 @@ function check_propagator(
                 success = false
             end
         end
+        if propagator.state ≢ Ψ
+            quiet ||
+                @error "$(px)`set_state!(propagator, state)` must return `propagator.state`."
+            success = false
+        end
     catch exc
         quiet || @error(
             "$(px)`set_state!(propagator, state)` must be defined.",
diff --git a/src/newton_propagator.jl b/src/newton_propagator.jl
index 355527e..43cf623 100644
--- a/src/newton_propagator.jl
+++ b/src/newton_propagator.jl
@@ -119,7 +119,6 @@ init_prop(state, generator, tlist, method::Val{:newton}; kwargs...) =
 function prop_step!(propagator::NewtonPropagator)
     @timeit_debug propagator.timing_data "prop_step!" begin
         Ψ = propagator.state
-        H = propagator.genop
         n = propagator.n  # index of interval we're going to propagate
         tlist = getfield(propagator, :tlist)
         (0 < n < length(tlist)) || return nothing
@@ -127,8 +126,8 @@ function prop_step!(propagator::NewtonPropagator)
         if propagator.backward
             dt = -dt
         end
-        _pwc_set_genop!(propagator, n)
         if propagator.inplace
+            H = _pwc_set_genop!(propagator, n)
             Newton.newton!(
                 Ψ,
                 H,
diff --git a/src/ode_function.jl b/src/ode_function.jl
index 5c1bd41..cf12a2b 100644
--- a/src/ode_function.jl
+++ b/src/ode_function.jl
@@ -45,7 +45,8 @@ H = evaluate(f.generator, t; vals_dict=vals_dict)
 where `vals_dict` may be a dictionary mapping controls to values (set as the
 parameters `p` of the underlying ODE solver).
 
-If [`QuantumPropagators.enable_timings()` has been called](@extref QuantumPropagators TimerOutputs),
+If [`QuantumPropagators.enable_timings()`](@ref
+QuantumPropagators.enable_timings) has been called,
 profiling data is collected in `f.timing_data`.
 """
 function ode_function(generator::GT, tlist; c=-1im, _timing_data=TimerOutput()) where {GT}
@@ -74,9 +75,9 @@ end
 
 function (f::QuantumODEFunction)(u, p, t)
     @timeit_debug f.timing_data "operator evaluation" begin
-        evaluate!(f.operator, f.generator, t; vals_dict=p)
+        H = evaluate(f.generator, t; vals_dict=p)
     end
     @timeit_debug f.timing_data "matrix-vector product" begin
-        return f.c * f.operator * u
+        return f.c * H * u
     end
 end
diff --git a/src/propagate.jl b/src/propagate.jl
index 37fb430..57d8db8 100644
--- a/src/propagate.jl
+++ b/src/propagate.jl
@@ -175,6 +175,7 @@ function propagate(
     inplace=true, # cf. default of init_prop
     for_expval=true,  # undocumented
     for_immutable_state=true,  # undocumented
+    for_mutable_operator=inplace,  # undocumented
     for_mutable_state=inplace,  # undocumented
     kwargs...
 )
@@ -208,6 +209,7 @@ function propagate(
             state=state,
             tlist=tlist,
             for_immutable_state,
+            for_mutable_operator,
             for_mutable_state,
             for_expval,
             atol,
diff --git a/src/propagator.jl b/src/propagator.jl
index b9bfb28..ec19759 100644
--- a/src/propagator.jl
+++ b/src/propagator.jl
@@ -17,7 +17,9 @@ implement the following interface.
 * `backward`: Boolean flag to indicate whether the propagation moves forward or
   backward in time
 * `inplace`: Boolean flag to indicate whether `propagator.state` is modified
-  in-place or is recreated by every call of `prop_step!` or `set_state!`.
+  in-place or is recreated by every call of `prop_step!` or `set_state!`. With
+  `inplace=false`, the propagator should generally avoid in-place operations,
+  such as calls to [`QuantumPropagators.Controls.evaluate!`](@ref).
 
 Concrete `Propagator` types may have additional properties or fields, but these
 should be considered private.
@@ -332,8 +334,9 @@ function prop_step! end
 set_state!(propagator, state)
 ```
 
-sets the `propagator.state` property. In order to mutate the current state
-after a call to [`prop_step!`](@ref), the following pattern is recommended:
+sets the `propagator.state` property and returns `propagator.state`. In order
+to mutate the current state after a call to [`prop_step!`](@ref), the following
+pattern is recommended:
 
 ```
 Ψ = propagator.state
@@ -366,9 +369,11 @@ function set_state!(propagator::AbstractPropagator, state)
         if propagator.inplace
             copyto!(propagator.state, state)
         else
-            setfield!(propagator, :state, state)
+            T = typeof(propagator.state)
+            setfield!(propagator, :state, convert(T, state))
         end
     end
+    return propagator.state
 end
 
 
diff --git a/src/pwc_utils.jl b/src/pwc_utils.jl
index 16dcd49..0f669f8 100644
--- a/src/pwc_utils.jl
+++ b/src/pwc_utils.jl
@@ -88,13 +88,14 @@ function _pwc_set_genop!(propagator::PiecewisePropagator, n)
     generator = getfield(propagator, :generator)
     tlist = propagator.tlist
     evaluate!(propagator.genop, generator, tlist, n; vals_dict=vals_dict)
+    return propagator.genop
 end
 
 function _pwc_get_genop(propagator::PiecewisePropagator, n)
     vals_dict = IdDict(c => propagator.parameters[c][n] for c in propagator.controls)
     generator = getfield(propagator, :generator)
     tlist = propagator.tlist
-    evaluate(generator, tlist, n; vals_dict=vals_dict)
+    return evaluate(generator, tlist, n; vals_dict=vals_dict)
 end
 
 
diff --git a/src/storage.jl b/src/storage.jl
index 8041348..0297ad2 100644
--- a/src/storage.jl
+++ b/src/storage.jl
@@ -3,7 +3,7 @@ module Storage
 using LinearAlgebra
 
 export init_storage, map_observables, map_observable, write_to_storage!
-export get_from_storage!
+export get_from_storage!, get_from_storage
 
 """Create a `storage` array for propagation.
 
@@ -162,8 +162,28 @@ extracts data from the `storage` for the i'th time slot. Inverse of
 developer's responsibility
 that [`init_storage`](@ref),  [`write_to_storage!`](@ref), and
 `get_from_storage!` are compatible.
+
+To extract immutable `data`, the non-in-place version
+
+```julia
+data = get_from_storage(storage, i)
+```
+
+can be used.
 """
 get_from_storage!(data, storage::AbstractVector, i) = copyto!(data, storage[i])
 get_from_storage!(data, storage::Matrix, i) = copyto!(data, storage[:, i])
 
+
+"""Obtain immutable data from storage.
+
+```julia
+data = get_from_storage(storage, i)
+```
+
+See [`get_from_storage!`](@ref).
+"""
+get_from_storage(storage::AbstractVector, i) = storage[i]
+get_from_storage(storage::Matrix, i) = storage[:, i]
+
 end
diff --git a/test/Project.toml b/test/Project.toml
index 492b191..3f6b600 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -40,5 +40,5 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
-Documenter = "1.1"
+Documenter = "1.4"
 julia = "1.9"
diff --git a/test/test_propagate.jl b/test/test_propagate.jl
index 6fabb90..8806893 100644
--- a/test/test_propagate.jl
+++ b/test/test_propagate.jl
@@ -1,8 +1,10 @@
 using Test
 using QuantumPropagators
+using QuantumPropagators: Cheby
 using QuantumPropagators.Storage
 using LinearAlgebra
 using UnicodePlots
+using StaticArrays: @SMatrix, SVector, @SVector
 
 @testset "TLS Rabi Cycling" begin
 
@@ -59,6 +61,78 @@ using UnicodePlots
 end
 
 
+@testset "Immutable TLS" begin
+
+    # We use Rabi cycling in a TLS as a test case, which allows to to compare
+    # the propagation with the known analytic solution.
+
+    SHOWPLOT = false
+
+    Ψ0 = @SVector ComplexF64[1, 0]
+    Ĥ = @SMatrix ComplexF64[
+         0   0.5
+        0.5   0
+    ]
+    tlist = collect(range(0, 1.5π, length=101)) # 3π/2 pulse
+
+    generator = (Ĥ,)
+
+    storage = init_storage(Ψ0, tlist)
+    @test isa(storage, Vector)
+    @test eltype(storage) == SVector{2,ComplexF64}
+
+    Ψ_out = propagate(Ψ0, generator, tlist; inplace=false, method=:expprop, storage=storage)
+    Ψ_expected = @SVector ComplexF64[-1/√2, -1im/√2]  # note the phases
+
+    pop0 = abs2.(hcat((Array.(storage))...))
+    SHOWPLOT && println(lineplot(tlist ./ π, pop0', ylim=[0, 1], title="fw prop"))
+
+    @test norm(Ψ_out - Ψ_expected) < 1e-12
+    @test pop0[end] ≈ 0.5
+
+    Ψ_out_cheby =
+        propagate(Ψ0, generator, tlist; inplace=false, method=Cheby, storage=storage,)
+    @test norm(Ψ_out_cheby - Ψ_expected) < 1e-12
+
+    # Propagating backward in time should exactly reverse the dynamics (since
+    # they are unitary). Thus, we should end up back at the initial state, and
+    # the stored states (being filled back-to-front) should exactly match the
+    # stored states from the forward propagation.
+
+    storage_bw = init_storage(Ψ0, tlist)
+    Ψ_out_bw = propagate(
+        Ψ_out,
+        generator,
+        tlist;
+        inplace=false,
+        method=:expprop,
+        backward=true,
+        storage=storage_bw
+    )
+
+    pop0_bw = abs2.(hcat((Array.(storage))...))
+    SHOWPLOT && println(lineplot(tlist ./ π, pop0_bw', ylim=[0, 1], title="bw prop"))
+
+    @test norm(Ψ_out_bw - Ψ0) < 1e-12
+    @test pop0_bw[1] ≈ 1.0
+
+    @test norm(storage - storage_bw) < 1e-12
+
+    Ψ_out_bw_cheby = propagate(
+        Ψ_out,
+        generator,
+        tlist;
+        inplace=false,
+        method=Cheby,
+        backward=true,
+        storage=storage_bw
+    )
+    @test norm(Ψ_out_bw_cheby - Ψ0) < 1e-12
+
+
+end
+
+
 @testset "Optomech propagation" begin
     include("optomech.jl")
     Ψ0 = ket(0, N_cav) ⊗ ket(2, N_mech)