diff --git a/CLAUDE.md b/CLAUDE.md
index 60b0815a..c9a25276 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -76,11 +76,20 @@ Uses SafeTestsets for isolated test execution. Tests are comprehensive and inclu
 
 Running `make test` prints out coverage data in a table.
 
+Tips for writing tests:
+
+- Use functions like `random_state_vector`, `random_dynamic_generator`, `random_matrix` from `QuantumControlTestUtils.RandomObjects` to generate random objects. These should always be given an explicit `StableRNG` as `rng`, with a unique seed
+
+- When a new seed is required, obtain one with `julia --project=test -e 'using QuantumControlTestUtils.RandomObjects: randseed; print(randseed())'`
+
+
 If necessary, detailed line-by-line coverage information can be obtained by running julia --project=test -e 'include("devrepl.jl"); generate_coverage_html()' after `make test`.
 This will produce html files inside the `coverage` subfolder, with `coverage/src` mirroring the structure of the `src` folder of `.jl` files. Lines with `<span class="tlaUNC">` are not covered. Ignore the raw tracefiles in the `.coverage` subfolder.
 
 ### Development Environment
+
 The project uses a sophisticated development setup:
+
 - Development REPL (devrepl.jl) with Revise.jl for hot reloading
 - Automatic dependency management via installorg.jl script
 - Integrated documentation building and serving
@@ -100,4 +109,6 @@ Each Julia function that is not explicitly private or has a name starting with a
 
 * When adding a new dependency to any `Project.toml` file, run `make distclean`, and then `make test/Manifest.toml`, `make docs/Manifest.toml`, etc. to recreate manifest files as necessary.
 
+* In order to get the documentation (docstring) of any function, e.g., `some_func`, in any package in the test environment, e.g., `SomeLib`, run `julia --project=. -e 'using REPL; using SomeLib; print(Base.doc(SomeLib.some_func))'`
+
 * Never commit any changes or ask to commit. I will always create git commits manually.
diff --git a/Project.toml b/Project.toml
index d3c0501c..bf3f866b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,17 +16,20 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extensions]
+QuantumPropagatorsExponentialUtilitiesExt = "ExponentialUtilities"
 QuantumPropagatorsODEExt = "OrdinaryDiffEq"
 QuantumPropagatorsRecursiveArrayToolsExt = "RecursiveArrayTools"
 QuantumPropagatorsStaticArraysExt = "StaticArrays"
 
 [compat]
 ArrayInterface = "7.0"
+ExponentialUtilities = "1.11"
 OffsetArrays = "1"
 OrdinaryDiffEq = "6.59"
 ProgressMeter = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index c62945d7..905dd99b 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,7 @@ DocInventories = "43dc2714-ed3b-44b5-b226-857eda1aa7de"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/make.jl b/docs/make.jl
index b0533025..9d809358 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,6 +3,7 @@ using Pkg
 using Documenter
 using QuantumPropagators
 using QuantumPropagators: AbstractPropagator, set_t!, set_state!
+import ExponentialUtilities  # ensure ExponentialUtilities extension is loaded
 import OrdinaryDiffEq  # ensure ODE extension is loaded
 using DocumenterCitations
 using DocumenterInterLinks
@@ -46,6 +47,7 @@ links = InterLinks(
     "ComponentArrays" => "https://sciml.github.io/ComponentArrays.jl/stable/",
     "RecursiveArrayTools" => "https://docs.sciml.ai/RecursiveArrayTools/stable/",
     "ArrayInterface" => "https://docs.sciml.ai/ArrayInterface/stable/",
+    "ExponentialUtilities" => "https://docs.sciml.ai/ExponentialUtilities/stable/",
     "qutip" => "https://qutip.readthedocs.io/en/qutip-5.0.x/",
 )
 
diff --git a/docs/src/methods.md b/docs/src/methods.md
index 0e02b1e7..f1cf31a1 100644
--- a/docs/src/methods.md
+++ b/docs/src/methods.md
@@ -18,10 +18,14 @@ As discussed in the [Overview](@ref overview_approaches), time propagation can b
    * [`Cheby`](@ref method_cheby) — expansion of ``\op{U} |Ψ⟩`` into Chebychev polynomials, valid if ``\op{H}`` has real eigenvalues
    * [`Newton`](@ref method_newton) – expansion of ``\op{U} |Ψ⟩`` into Newton polynomials, valid if ``\op{H}`` has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian)
 
-   The `ExpProp` method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomials series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the is truncated as soon as some desired precision is reached, which is machine precision by default).
+   The `ExpProp` method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomial series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the series is truncated as soon as some desired precision is reached, which is machine precision by default).
 
    Note that this high precision is *within the piecewise-constant approximation*. The discretization itself may introduce a non-negligible error compared to the time-continuous dynamics. There is tradeoff: A smaller `dt` decreases the discretization error, but the polynomial expansions are more effective with larger time steps.
 
+   There is also support for _external_ packages to efficiently evaluate the application of the piecewise-constant time evolution operator:
+
+   * [`ExponentialUtilities`](@ref method_exponentialutilities) – applies ``\op{U} |Ψ⟩`` using Krylov methods without explicitly forming ``\op{U}``
+
 2. By solving the equation of motion explicitly with an ODE solver.
 
    We support the use of any of the ODE solvers in [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/):
@@ -97,6 +101,58 @@ See [GoerzPhd2015; Chapter 3.2.1](@cite) for a detailed description of the metho
 * Closed quantum systems with piecewise constant dynamics
 
 
+## [`ExponentialUtilities`](@id method_exponentialutilities)
+
+The method requires that the [ExponentialUtilities.jl](https://docs.sciml.ai/ExponentialUtilities/stable/)
+package is loaded
+
+```
+using ExponentialUtilities
+```
+
+and then passed as `method=ExponentialUtilities` to
+[`propagate`](@ref) or [`init_prop`](@ref):
+
+```@docs
+init_prop(state, generator, tlist, method::Val{:ExponentialUtilities}; kwargs...)
+```
+
+This method evaluates ``\exp(-i \op{H} dt) |Ψ⟩`` via a Krylov
+[`expv`](@extref ExponentialUtilities :jl:function:`ExponentialUtilities.expv`)
+algorithm without explicitly forming the matrix exponential. It builds a
+Krylov subspace (via Arnoldi or Lanczos iteration) and then computes the
+action of the matrix exponential on the state within that subspace.
+This is often a good fit for larger systems or matrix-free operators where
+direct matrix exponentiation is too costly.
+
+The propagator requires that states and operators satisfy the `AbstractArray`
+interface (specifically,
+[`supports_vector_interface`](@ref QuantumPropagators.Interfaces.supports_vector_interface)
+for states and
+[`supports_matrix_interface`](@ref QuantumPropagators.Interfaces.supports_matrix_interface)
+for operators).
+
+**Advantages**
+
+* Avoids explicit construction of ``\op{U}``
+* Works with matrix-free operators that support `mul!`
+* Good scaling for large sparse systems
+* Supports both Hermitian (Lanczos) and non-Hermitian (Arnoldi) generators
+
+**Disadvantages**
+
+* Requires ExponentialUtilities.jl
+* Performance depends on Krylov subspace dimension and operator structure
+* Requires operators and states to implement an array interface, see
+  [`QuantumPropagators.Interfaces.supports_matrix_interface`](@ref) and
+  [`QuantumPropagators.Interfaces.supports_vector_interface`](@ref), respectively
+
+**When to use**
+
+* Large, sparse, or matrix-free generators
+* Systems where ``\op{U}`` is too expensive to form explicitly
+
+
 ## [`Newton`](@id method_newton)
 
 The method should be loaded with
diff --git a/docs/src/overview.md b/docs/src/overview.md
index 25d9fe9f..c231f172 100644
--- a/docs/src/overview.md
+++ b/docs/src/overview.md
@@ -77,7 +77,7 @@ This form differs from most textbooks by a factor of ``i``, but has the benefit
 
 There are two fundamental approaches to solving the Schrödinger equation (or any equation of motion):
 
-1. We can analytically solve the Schrödinger equation and then numerically evaluate the solution. Mathematically, this is the application of the time evolution operator ``\op{H}`` as ``|Ψ(t+dt)⟩ = \op{U}(t) |Ψ(t)⟩``. For a piecewise-constant ``\op{H}(t)``where there is a time-independent ``\op{H}`` in the interval ``[t, t+dt]``, the time evolution operator is well-known to be
+1. We can analytically solve the Schrödinger equation and then numerically evaluate the solution. Mathematically, this is the application of the time evolution operator ``\op{U}`` as ``|Ψ(t+dt)⟩ = \op{U}(t) |Ψ(t)⟩``. For a piecewise-constant ``\op{H}(t)``where there is a time-independent ``\op{H}`` in the interval ``[t, t+dt]``, the time evolution operator is well-known to be
 
    ```math
    \op{U} = \exp[-i \op{H} dt]\,.
@@ -112,7 +112,7 @@ Furthermore, as a fallback for very small system or for debugging,
 
 * `method=ExpProp`: Explicitly construct the time evolution operator by matrix exponentiation and apply it to the state.
 
-If the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) or (equivalently) th [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package is loaded, `QuantumPropagators` can delegate to it:
+If the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) or (equivalently) the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package is loaded, `QuantumPropagators` can delegate to it:
 
 ```julia
 using OrdinaryDiffEq
@@ -122,7 +122,11 @@ allows to pass
 
 * `method=OrdinaryDiffEq`: Use [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) as a backend with any of the [algorithms available for ODEs in DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/), `alg=Tsit5()` by default.
 
-Unlike any of the built-in methods, `OrdinaryDiffEq` is able to propagate for time-continuous generators. This is the default for that propagator (`pwc=false`). By setting `pwc=true` or `piecewise=true`) the ODE solvers can also be used for piecewise-constant Hamiltonians or Liouvillians, providing an alternative to the built-in `method=Cheby` and `method=Newton`.
+Unlike any of the built-in methods, `OrdinaryDiffEq` is able to propagate for time-continuous generators. This is the default for that propagator (`pwc=false`). By setting `pwc=true` or `piecewise=true`, the ODE solvers can also be used for piecewise-constant Hamiltonians or Liouvillians, providing an alternative to the built-in `method=Cheby` and `method=Newton`.
+
+If the [ExponentialUtilities.jl](https://docs.sciml.ai/ExponentialUtilities/stable/expv/) package is loaded, this enables
+
+* `method=ExponentialUtilities`: Evaluate the application of the piecewise-constant time evolution operator via [`ExponentialUtilities.expv`](@extref).
 
 See the more extended discussion of [Propagation Methods](@ref) for more details.
 
@@ -182,7 +186,7 @@ Most propagators support both an in-place and a not-in-place mode. These modes c
 
 In-place operations can be dramatically more efficient for large Hilbert space dimensions. On the other hand, not-in-place operations can be more efficient for small Hilbert spaces, in particular when a [static vector](@extref StaticArrays `SVector`) can be used to represent the state. Moreover, frameworks for automatic differentiation such as [Zygote](https://fluxml.ai/Zygote.jl/stable/) do not support in-place operations.
 
-When using custom structs for states, operators, or generators, the struct itself not not need to be mutable (according to [`Base.ismutable`](@extref Julia)) in order to support `inplace=true`. It only must support the in-place operations defined in the [formal interface](@ref QuantumPropagatorsInterfacesAPI) and indicate that support by defining [`QuantumPropagators.Interfaces.supports_inplace`](@ref).
+When using custom structs for states, operators, or generators, the struct itself does not need to be mutable (according to [`Base.ismutable`](@extref Julia)) in order to support `inplace=true`. It only must support the in-place operations defined in the [formal interface](@ref QuantumPropagatorsInterfacesAPI) and indicate that support by defining [`QuantumPropagators.Interfaces.supports_inplace`](@ref).
 Typically, in-place operations on immutable custom structs involve mutating the mutable properties of that struct.
 
 
diff --git a/ext/QuantumPropagatorsExponentialUtilitiesExt.jl b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl
new file mode 100644
index 00000000..1a63ece4
--- /dev/null
+++ b/ext/QuantumPropagatorsExponentialUtilitiesExt.jl
@@ -0,0 +1,212 @@
+module QuantumPropagatorsExponentialUtilitiesExt
+
+using LinearAlgebra
+using TimerOutputs: @timeit_debug, TimerOutput
+using ExponentialUtilities
+
+using QuantumPropagators:
+    QuantumPropagators,
+    ExponentialUtilitiesPropagator,
+    _pwc_get_max_genop,
+    _pwc_get_genop,
+    _pwc_set_genop!,
+    _pwc_advance_time!,
+    _pwc_process_parameters
+using QuantumPropagators.Controls: get_controls
+using QuantumPropagators.Interfaces: supports_inplace
+import QuantumPropagators: init_prop, prop_step!
+
+
+"""
+```julia
+using ExponentialUtilities
+
+expv_propagator = init_prop(
+    state,
+    generator,
+    tlist;
+    method=ExponentialUtilities,
+    inplace=QuantumPropagators.Interfaces.supports_inplace(state),
+    backward=false,
+    verbose=false,
+    parameters=nothing,
+    expv_kwargs=(;),
+    _...
+)
+```
+
+initializes an [`ExponentialUtilitiesPropagator`](@ref).
+
+# Method-specific keyword arguments
+
+* `expv_kwargs`: NamedTuple of keyword arguments forwarded to the underlying
+  [`ExponentialUtilities`](https://docs.sciml.ai/ExponentialUtilities/stable/)
+  routines. For in-place propagation, these are passed to both `arnoldi!`
+  (Krylov subspace construction) and `expv!` (exponentiation). For
+  out-of-place propagation, they are passed to `expv`, which forwards them
+  internally to `arnoldi`. The most useful keyword arguments are:
+
+  * `m::Int`: Maximum Krylov subspace dimension. Defaults to
+    `min(30, size(A, 1))`. Larger values improve accuracy but increase memory
+    and computation cost per step.
+  * `ishermitian::Bool`: Whether the operator is Hermitian. Defaults to
+    `LinearAlgebra.ishermitian(A)`. When `true`, the cheaper Lanczos iteration
+    is used instead of the full Arnoldi process. Note that the propagator
+    passes `-𝕚 dt` as the time argument, so a Hermitian generator ``Ĥ`` still
+    results in a non-Hermitian product ``-𝕚 Ĥ dt``; however,
+    `ExponentialUtilities` handles the complex time scaling internally.
+  * `tol::Real`: Tolerance for happy breakdown. Defaults to `1e-7`. The
+    Arnoldi iteration stops early when the norm of the next Krylov vector
+    drops below `tol * opnorm`.
+  * `opnorm`: Operator norm of `A`. By default, this is computed
+    automatically. Supplying it avoids redundant norm computations.
+  * `iop::Int`: Incomplete orthogonalization procedure length. Defaults to `0`
+    (full orthogonalization). A positive value limits the number of previous
+    Krylov vectors used for orthogonalization, reducing cost at the expense of
+    numerical stability.
+  * `mode::Symbol`: Termination strategy for `expv`. Either
+    `:happy_breakdown` (default) or `:error_estimate`. In `:happy_breakdown`
+    mode, the iteration relies on early termination when the Krylov subspace
+    captures the relevant dynamics. In `:error_estimate` mode, an adaptive
+    step-size strategy is used with additional tolerance parameters `rtol`
+    (relative tolerance, defaults to `√tol`).
+"""
+function init_prop(
+    state,
+    generator,
+    tlist,
+    method::Val{:ExponentialUtilities};
+    inplace = supports_inplace(state),
+    backward = false,
+    verbose = false,
+    parameters = nothing,
+    expv_kwargs = (;),
+    _...
+)
+    tlist = convert(Vector{Float64}, tlist)
+    controls = get_controls(generator)
+    G = _pwc_get_max_genop(generator, controls, tlist)
+
+    parameters = _pwc_process_parameters(parameters, controls, tlist)
+    timing_data = TimerOutput()
+
+    n = 1
+    t = tlist[1]
+    if backward
+        n = length(tlist) - 1
+        t = float(tlist[n+1])
+    end
+
+    if inplace
+        A₀ = QuantumPropagators.Controls.evaluate(generator, tlist, n)
+        is_herm = get(expv_kwargs, :ishermitian, ishermitian(A₀))
+        if haskey(expv_kwargs, :mode) && expv_kwargs[:mode] == :error_estimate && !is_herm
+            throw(
+                ArgumentError(
+                    "`mode=:error_estimate` in `expv_kwargs` requires a " *
+                    "Hermitian generator. Use the default " *
+                    "`mode=:happy_breakdown` for non-Hermitian generators."
+                )
+            )
+        end
+        Ks = ExponentialUtilities.arnoldi(A₀, state; expv_kwargs...)
+        if haskey(expv_kwargs, :mode) && expv_kwargs[:mode] == :error_estimate
+            cache = ExponentialUtilities.get_subspace_cache(Ks)
+        else
+            T = promote_type(eltype(A₀), eltype(state))
+            cache_type = is_herm ? real(T) : T
+            cache = ExponentialUtilities.ExpvCache{cache_type}(Ks.m)
+        end
+    else
+        Ks = nothing
+        cache = nothing
+    end
+
+    GT = typeof(generator)
+    OT = typeof(G)
+    ST = typeof(state)
+    KST = typeof(Ks)
+    CT = typeof(cache)
+
+    return ExponentialUtilitiesPropagator{GT,OT,ST,KST,CT}(
+        generator,
+        inplace ? copy(state) : state,
+        t,
+        n,
+        tlist,
+        parameters,
+        controls,
+        G,
+        Ks,
+        cache,
+        backward,
+        inplace,
+        expv_kwargs,
+        timing_data,
+    )
+end
+
+
+function prop_step!(propagator::ExponentialUtilitiesPropagator)
+    @timeit_debug propagator.timing_data "prop_step!" begin
+        n = propagator.n
+        tlist = getfield(propagator, :tlist)
+        (0 < n < length(tlist)) || return nothing
+        dt = tlist[n+1] - tlist[n]
+        if propagator.backward
+            dt = -dt
+        end
+
+        if propagator.inplace
+            if supports_inplace(propagator.genop)
+                _pwc_set_genop!(propagator, n)
+                H = propagator.genop
+            else
+                H = _pwc_get_genop(propagator, n)
+            end
+            @timeit_debug propagator.timing_data "expv" begin
+                if haskey(propagator.expv_kwargs, :mode) &&
+                   propagator.expv_kwargs[:mode] == :error_estimate
+                    ExponentialUtilities.expv!(
+                        propagator.state,
+                        -1im * dt,
+                        H,
+                        propagator.state,
+                        propagator.Ks,
+                        propagator.cache;
+                        propagator.expv_kwargs...,
+                    )
+                else
+                    ExponentialUtilities.arnoldi!(
+                        propagator.Ks,
+                        H,
+                        propagator.state;
+                        propagator.expv_kwargs...,
+                    )
+                    ExponentialUtilities.expv!(
+                        propagator.state,
+                        -1im * dt,
+                        propagator.Ks;
+                        cache = propagator.cache,
+                    )
+                end
+            end
+        else
+            H = _pwc_get_genop(propagator, n)
+            @timeit_debug propagator.timing_data "expv" begin
+                Ψ = ExponentialUtilities.expv(
+                    -1im * dt,
+                    H,
+                    propagator.state;
+                    propagator.expv_kwargs...
+                )
+            end
+            setfield!(propagator, :state, convert(typeof(propagator.state), Ψ))
+        end
+
+        _pwc_advance_time!(propagator)
+        return propagator.state
+    end
+end
+
+end
diff --git a/src/QuantumPropagators.jl b/src/QuantumPropagators.jl
index 4d5573ac..67a9f565 100644
--- a/src/QuantumPropagators.jl
+++ b/src/QuantumPropagators.jl
@@ -58,6 +58,7 @@ end
 #! format: on
 
 include("pwc_utils.jl")
+include("exponential_utilities_propagator.jl")
 include("cheby_propagator.jl")
 include("newton_propagator.jl")
 include("exp_propagator.jl")
diff --git a/src/exponential_utilities_propagator.jl b/src/exponential_utilities_propagator.jl
new file mode 100644
index 00000000..d0d25fc7
--- /dev/null
+++ b/src/exponential_utilities_propagator.jl
@@ -0,0 +1,27 @@
+using TimerOutputs: TimerOutput
+
+
+"""Propagator for Krylov `expv` propagation (`method=ExponentialUtilities`).
+
+This is a [`PWCPropagator`](@ref). Methods that depend on
+`ExponentialUtilities` are loaded via Julia package extensions.
+"""
+mutable struct ExponentialUtilitiesPropagator{GT,OT,ST,KST,CT} <: PWCPropagator
+    const generator::GT
+    state::ST
+    t::Float64  # time at which current `state` is defined
+    n::Int64 # index of next interval to propagate
+    const tlist::Vector{Float64}
+    parameters::AbstractDict
+    controls
+    genop::OT
+    Ks::KST
+    cache::CT
+    backward::Bool
+    inplace::Bool
+    expv_kwargs::NamedTuple
+    const timing_data::TimerOutput
+end
+
+
+set_t!(propagator::ExponentialUtilitiesPropagator, t) = _pwc_set_t!(propagator, t)
diff --git a/src/generators.jl b/src/generators.jl
index 7f96e579..e2e70cfc 100644
--- a/src/generators.jl
+++ b/src/generators.jl
@@ -176,6 +176,7 @@ function Base.copyto!(tgt::Operator, src::Operator)
     copyto!(tgt.coeffs, src.coeffs)
 end
 
+# All ops in an Operator must have the same size (they are summed together)
 Base.size(O::Operator) = size(O.ops[1])
 Base.size(O::Operator, dim::Integer) = size(O.ops[1], dim)
 Base.eltype(::Type{Operator{OT,CT}}) where {OT,CT} = promote_type(eltype(OT), CT)
@@ -255,7 +256,7 @@ function Base.show(io::IO, O::ScaledOperator{CT,OT}) where {CT,OT}
     print(io, "ScaledOperator($(O.coeff), $(O.operator))")
 end
 
-function Base.show(io::IO, ::MIME"text/plain", O::ScaledOperator{CT,OT}) where {CT,OT}
+function Base.show(io::IO, ::MIME"text/plain", O::ScaledOperator{CT,<:Operator}) where {CT}
     Base.summary(io, O)
     println(io, "\n operator.ops::$(typeof(O.operator.ops)):")
     for op in O.operator.ops
diff --git a/test/Project.toml b/test/Project.toml
index 3df51a9c..fa8c5de6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -9,6 +9,7 @@ DocInventories = "43dc2714-ed3b-44b5-b226-857eda1aa7de"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 GRAPE = "6b52fcaf-80fe-489a-93e9-9f92080510be"
 GRAPELinesearchAnalysis = "290eba36-e2d8-4488-81b6-f66cc44f2186"
 IOCapture = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
diff --git a/test/runtests.jl b/test/runtests.jl
index cc3ac140..3268b7ea 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -78,6 +78,11 @@ using SafeTestsets
         include("test_expprop.jl")
     end
 
+    println("\n* ExponentialUtilities baseline (test_exponential_utilities.jl):")
+    @time @safetestset "ExponentialUtilities (baseline)" begin
+        include("test_exponential_utilities.jl")
+    end
+
     println("\n* Propagator Interfaces (test_prop_interfaces.jl):")
     @time @safetestset "Propagator Interfaces" begin
         include("test_prop_interfaces.jl")
@@ -93,6 +98,11 @@ using SafeTestsets
         include("test_propagate_sequence.jl")
     end
 
+    println("\n* ExponentialUtilities propagation (test_exputils.jl):")
+    @time @safetestset "ExponentialUtilities (propagation)" begin
+        include("test_exputils.jl")
+    end
+
     println("\n* Time-dependent observables (test_timedependent_observables.jl):")
     @time @safetestset "Time-dependent observables" begin
         include("test_timedependent_observables.jl")
diff --git a/test/test_exponential_utilities.jl b/test/test_exponential_utilities.jl
new file mode 100644
index 00000000..b1bcb818
--- /dev/null
+++ b/test/test_exponential_utilities.jl
@@ -0,0 +1,651 @@
+using Test
+
+using QuantumPropagators.Controls: evaluate
+using QuantumPropagators.Generators: ScaledOperator
+using StableRNGs: StableRNG
+using QuantumControlTestUtils.RandomObjects:
+    random_state_vector, random_dynamic_generator, random_matrix
+using LinearAlgebra: norm, ishermitian
+import ExponentialUtilities
+import StaticArrays: SMatrix, SVector
+import OffsetArrays: OffsetVector
+using BenchmarkTools: @btimed
+
+RUN_BENCHMARKS = false
+
+
+@testset "expv with dense complex non-Hermitian matrix" begin
+
+    N = 10
+    rng = StableRNG(677968056)
+    A = random_matrix(N; rng, complex = true, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with dense real non-Hermitian matrix" begin
+
+    N = 10
+    rng = StableRNG(2945201349)
+    A = random_matrix(N; rng, complex = false, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with dense complex Hermitian matrix and complex/real time" begin
+
+    # Hermitian A triggers Lanczos instead of Arnoldi
+    N = 10
+    rng = StableRNG(697691674)
+    A = random_matrix(N; rng, hermitian = true)
+    @test ishermitian(A)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b1 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀)
+    end
+
+    # We can also check against the code path where we explicitly disable
+    # `ishermitian` detection
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; ishermitian = false)
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b2 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀; ishermitian = false)
+        @test b2.time > b1.time
+    end
+
+    # Instead of using complex time, we can also use real time and absorb the
+    # factor `-1im` in `A`
+    O = ScaledOperator(-1im, A)
+    @test !ishermitian(O)
+    Ψ = ExponentialUtilities.expv(dt, O, Ψ₀)
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b3 = @btimed ExponentialUtilities.expv($dt, $O, $Ψ₀)
+        @test b3.time > b1.time
+    end
+
+end
+
+
+@testset "expv with dense real symmetric matrix" begin
+
+    # Symmetric (Hermitian) `A` triggers Lanczos instead of Arnoldi
+    N = 10
+    rng = StableRNG(1056344208)
+    A = random_matrix(N; rng, complex = false, hermitian = true)
+    @test ishermitian(A)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b1 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀)
+    end
+    # We can also check against the code path where we explicitly disable
+    # `ishermitian` detection
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; ishermitian = false)
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b2 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀; ishermitian = false)
+        @test b2.time > b1.time
+    end
+
+end
+
+
+@testset "expv with sparse complex non-Hermitian matrix" begin
+
+    N = 50
+    rng = StableRNG(1788065498)
+    A = random_matrix(N; rng, density = 0.3, hermitian = false)
+    @test !ishermitian(A)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * Matrix(A)) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with sparse complex Hermitian matrix" begin
+
+    # Sparse Hermitian A triggers Lanczos
+    N = 50
+    rng = StableRNG(2648750134)
+    A = random_matrix(N; rng, density = 0.3, hermitian = true)
+    @test ishermitian(A)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * Matrix(A)) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with Operator" begin
+
+    N = 10
+    rng = StableRNG(677968056)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng)
+    Ψ₀ = random_state_vector(N; rng)
+    A = evaluate(H, tlist, 1)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * Array(A)) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with StaticArrays SMatrix and SVector" begin
+
+    N = 4
+    rng = StableRNG(814242294)
+    A = SMatrix{N,N}(random_matrix(N; rng, hermitian = true))
+    Ψ₀ = SVector{N}(random_state_vector(N; rng))
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    Ψ_expected = exp(t * Matrix(A)) * Vector(Ψ₀)
+    @test Ψ isa SVector{N}
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with pre-built KrylovSubspace" begin
+
+    # arnoldi builds the Krylov subspace once; expv(t, Ks) reuses it,
+    # which is efficient when applying the same operator at many time steps
+    N = 10
+    rng = StableRNG(3750746111)
+    A = random_matrix(N; rng, hermitian = true)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ks = ExponentialUtilities.arnoldi(A, Ψ₀)
+    Ψ = ExponentialUtilities.expv(t, Ks)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! in-place with pre-built KrylovSubspace" begin
+
+    N = 10
+    rng = StableRNG(1187394246)
+    A = random_matrix(N; rng, hermitian = true)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ks = ExponentialUtilities.arnoldi(A, Ψ₀)
+    Ψ = similar(Ψ₀, promote_type(typeof(t), eltype(A), eltype(Ψ₀)))
+    ExponentialUtilities.expv!(Ψ, t, Ks)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! in-place with dense complex Hermitian matrix and complex/real time" begin
+
+    N = 10
+    rng = StableRNG(3503257612)
+    A = random_matrix(N; rng, hermitian = true)
+    @test ishermitian(A)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ_expected = exp(t * A) * Ψ₀
+
+    # Baseline: expv
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+    @test norm(Ψ - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        b0 = @btimed ExponentialUtilities.expv($t, $A, $Ψ₀)
+    end
+
+    # Lanczos path: arnoldi detects ishermitian(A) and builds a symmetric
+    # tridiagonal H; expv! then uses the fast SymTridiagonal eigen path.
+    # For Hermitian A the H matrix is real: ExpvCache element type = real(eltype(A)).
+    Ks1 = ExponentialUtilities.arnoldi(A, Ψ₀)
+    w = similar(Ψ₀)
+    ExponentialUtilities.expv!(w, t, Ks1)
+    @test norm(w - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        cache1 = ExponentialUtilities.ExpvCache{real(eltype(A))}(Ks1.m)
+        b1 = @btimed ExponentialUtilities.expv!($w, $t, $Ks1; cache = $cache1)
+        @test b1.time < b0.time
+    end
+
+    # Arnoldi path: force full Arnoldi; H is upper Hessenberg, so expv! falls
+    # back to the generic Higham matrix exponential, which is slower.
+    # Arnoldi (including when forced on a Hermitian A) builds complex H:
+    # ExpvCache element type = eltype(A).
+    Ks2 = ExponentialUtilities.arnoldi(A, Ψ₀; ishermitian = false)
+    ExponentialUtilities.expv!(w, t, Ks2)
+    @test norm(w - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        cache2 = ExponentialUtilities.ExpvCache{eltype(A)}(Ks2.m)
+        b2 = @btimed ExponentialUtilities.expv!($w, $t, $Ks2; cache = $cache2)
+        # The timings here are a lot flakier, so we don't
+        # @test b2.time > b1.time
+    end
+
+    # Real-time path: absorb `-1im` into the operator via ScaledOperator;
+    # the scaled operator is not Hermitian so arnoldi uses Arnoldi iteration,
+    # giving an upper Hessenberg H (ComplexF64), same as Ks2.
+    O = ScaledOperator(-1im, A)
+    @test !ishermitian(O)
+    Ks3 = ExponentialUtilities.arnoldi(O, Ψ₀)
+    ExponentialUtilities.expv!(w, dt, Ks3)
+    @test norm(w - Ψ_expected) < 1e-12
+    if RUN_BENCHMARKS
+        cache3 = ExponentialUtilities.ExpvCache{eltype(O)}(Ks3.m)
+        b3 = @btimed ExponentialUtilities.expv!($w, $dt, $Ks3; cache = $cache3)
+        # The timings here are a lot flakier, so we don't
+        # @test b3.time > b1.time
+    end
+
+end
+
+
+@testset "expv propagation loop over time steps" begin
+
+    # Propagate Ψ through the first few time steps of a random dynamic
+    # generator, reusing a single KrylovSubspace.
+    # At each step: arnoldi! updates Ks in-place for the current (A, Ψ),
+    # then expv(t, Ks) returns the new Ψ without mutating the input.
+    N = 10
+    rng = StableRNG(1373376837)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    # Loop with arnoldi! + expv(t, Ks), reusing Ks across all steps
+    A₀ = evaluate(H, tlist, 1)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)
+    Ψ = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        ExponentialUtilities.arnoldi!(Ks, A, Ψ)
+        Ψ = ExponentialUtilities.expv(t, Ks)
+    end
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! propagation loop over time steps" begin
+
+    # Propagate Ψ through the first few time steps of a random dynamic
+    # generator, reusing a single KrylovSubspace and ExpvCache.
+    # At each step n: arnoldi! updates Ks in-place for the current (A, Ψ),
+    # then expv! overwrites Ψ. arnoldi! reads b before expv! writes w, so
+    # passing the same vector for both is safe.
+    # Ks and cache must remain type-consistent across the loop: the H element
+    # type is fixed at construction (real for Lanczos, complex for Arnoldi),
+    # so all arnoldi! calls must use the same ishermitian setting as the
+    # initial arnoldi call.
+    N = 10
+    rng = StableRNG(2266238742)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    # Loop with arnoldi! + expv!, reusing Ks and cache across all steps.
+    # The H element type is real(T) for Lanczos and T for Arnoldi, where
+    # T = promote_type(eltype(A), eltype(Ψ)) is the Krylov basis vector type.
+    # Hermitian A₀ → Lanczos → real H → ExpvCache{real(T)}
+    A₀ = evaluate(H, tlist, 1)
+    @test ishermitian(A₀)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)
+    T = promote_type(eltype(A₀), eltype(Ψ₀))
+    cache = ExponentialUtilities.ExpvCache{real(T)}(Ks.m)
+    Ψ = copy(Ψ₀)
+    if RUN_BENCHMARKS
+        # One might think that these in-place versions are completely
+        # non-allocating. Unfortunately, that's not the case. The allocations
+        # don't disappear if this is wrapped into a function.
+        b1 = @btimed ExponentialUtilities.arnoldi!($Ks, $A₀, $Ψ)
+        @test_broken b1.alloc == 0
+        b2 = @btimed ExponentialUtilities.expv!($Ψ, $t, $Ks; cache = $cache)
+        @test_broken b2.alloc == 0
+    end
+    Ψ = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        ExponentialUtilities.arnoldi!(Ks, A, Ψ)
+        ExponentialUtilities.expv!(Ψ, t, Ks; cache = cache)
+    end
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! propagation loop over time steps with non-Hermitian generator" begin
+
+    # Non-Hermitian generators require Arnoldi iteration (not Lanczos), giving
+    # an upper Hessenberg H matrix (ComplexF64). The ExpvCache element type
+    # must be ComplexF64 accordingly. Otherwise the pattern is the same as for
+    # Hermitian generators.
+    N = 10
+    rng = StableRNG(1573710086)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    # Loop with arnoldi! + expv!, reusing Ks and cache across all steps.
+    # The H element type is T for Arnoldi, where
+    # T = promote_type(eltype(A), eltype(Ψ)) is the Krylov basis vector type.
+    # Non-Hermitian A₀ → Arnoldi → H element type T → ExpvCache{T}
+    A₀ = evaluate(H, tlist, 1)
+    @test !ishermitian(A₀)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)
+    T = promote_type(eltype(A₀), eltype(Ψ₀))
+    cache = ExponentialUtilities.ExpvCache{T}(Ks.m)
+    Ψ = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        @test !ishermitian(A)
+        ExponentialUtilities.arnoldi!(Ks, A, Ψ)
+        ExponentialUtilities.expv!(Ψ, t, Ks; cache = cache)
+    end
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! propagation loop using error-estimate interface" begin
+
+    # Variation of the propagation loop using expv!(w, t, A, b, Ks, cache):
+    # this interface builds the Krylov subspace internally via Lanczos and
+    # terminates adaptively when Saad's error estimate drops below the
+    # requested tolerance (default atol=1e-8, rtol=1e-4).
+    # Only Hermitian operators are supported: get_subspace_cache requires a
+    # Lanczos (real H) subspace, which arnoldi builds automatically when
+    # ishermitian(A) is true.
+    # Passing w == b (same Ψ vector) is safe: b is fully consumed into
+    # V[:,1] before w is written at the end of the call.
+    N = 10
+    rng = StableRNG(106614078)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    A₀ = evaluate(H, tlist, 1)
+    @test ishermitian(A₀)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)  # auto-detects Hermitian → Lanczos
+    cache = ExponentialUtilities.get_subspace_cache(Ks)
+    Ψ = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        @test ishermitian(A)
+        ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache)
+    end
+    # accuracy is governed by the per-step rtol=1e-4 default
+    @test norm(Ψ - Ψ_expected) < 1e-3
+
+end
+
+
+@testset "expv! propagation loop using error-estimate interface with tight tolerance" begin
+
+    # Same error-estimate interface, but with atol=rtol=1e-14. The termination
+    # condition ε = atol + rtol*‖b‖ is now so tight that the algorithm uses
+    # all m = min(Ks.maxiter, N) Lanczos steps, spanning the full N-dimensional
+    # Krylov space and giving near-machine-precision results.
+    N = 10
+    rng = StableRNG(1653744503)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    A₀ = evaluate(H, tlist, 1)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)  # auto-detects Hermitian → Lanczos
+    cache = ExponentialUtilities.get_subspace_cache(Ks)
+    Ψ = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache; atol = 1e-14, rtol = 1e-14)
+    end
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv! propagation loop using error-estimate interface with non-Hermitian generator" begin
+
+    # The error-estimate interface expv!(w, t, A, b, Ks, cache) only supports
+    # Hermitian operators. For non-Hermitian generators, get_subspace_cache
+    # throws an ErrorException because HermitianSubspaceCache requires a
+    # Lanczos subspace with real H. Additionally, the method signature
+    # requires KrylovSubspace{T,B,B} (second and third type params equal),
+    # which only holds for Lanczos (Float64, Float64) not Arnoldi
+    # (ComplexF64, Float64).
+    N = 10
+    rng = StableRNG(208995058)
+    tlist = collect(range(0, 10, length = 101))
+    H = random_dynamic_generator(N, tlist; rng, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = tlist[2] - tlist[1]
+    t = -1im * dt
+    n_steps = 5
+
+    # Reference: sequential dense matrix exponentials
+    Ψ_expected = copy(Ψ₀)
+    for n = 1:n_steps
+        A = evaluate(H, tlist, n)
+        Ψ_expected = exp(t * Array(A)) * Ψ_expected
+    end
+
+    # get_subspace_cache throws ErrorException for non-Hermitian Ks;
+    # a MethodError on expv! dispatch would also occur since the method
+    # requires KrylovSubspace{T,B,B} but non-Hermitian Ks has B=ComplexF64 ≠ Float64
+    A₀ = evaluate(H, tlist, 1)
+    @test !ishermitian(A₀)
+    Ks = ExponentialUtilities.arnoldi(A₀, Ψ₀)
+    @test_broken begin
+        cache = ExponentialUtilities.get_subspace_cache(Ks)
+        Ψ = copy(Ψ₀)
+        for n = 1:n_steps
+            A = evaluate(H, tlist, n)
+            ExponentialUtilities.expv!(Ψ, t, A, Ψ, Ks, cache)
+        end
+        norm(Ψ - Ψ_expected) < 1e-3
+    end
+
+end
+
+
+@testset "expv with mode=:error_estimate" begin
+
+    # mode=:error_estimate uses adaptive Krylov iteration with an on-the-fly
+    # error estimate; default rtol ≈ 3e-4 (= sqrt(tol) with tol=1e-7)
+    N = 10
+    rng = StableRNG(3050745481)
+    A = random_matrix(N; rng, hermitian = true)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; mode = :error_estimate)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-4
+
+end
+
+
+@testset "expv with custom Krylov subspace size m" begin
+
+    N = 20
+    rng = StableRNG(1375059820)
+    A = random_matrix(N; rng, hermitian = true)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; m = 8)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv with iop (incomplete orthogonalization)" begin
+
+    # iop=k uses only the last k vectors for re-orthogonalization,
+    # reducing cost for large operators at the price of some accuracy
+    N = 20
+    rng = StableRNG(1918914874)
+    A = random_matrix(N; rng, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; iop = 2)
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-10
+
+end
+
+
+@testset "expv with custom opnorm" begin
+
+    # opnorm can be provided to avoid recomputing it internally,
+    # useful when the norm is expensive or already known
+    N = 10
+    rng = StableRNG(1355318934)
+    A = random_matrix(N; rng, hermitian = false)
+    Ψ₀ = random_state_vector(N; rng)
+    dt = 0.1
+    t = -1im * dt
+    Ψ = ExponentialUtilities.expv(t, A, Ψ₀; opnorm = norm(A, Inf))
+    Ψ_expected = exp(t * A) * Ψ₀
+    @test norm(Ψ - Ψ_expected) < 1e-12
+
+end
+
+
+@testset "expv return type with real state and complex time" begin
+
+    # When the state vector has real eltype but t is complex, expv returns a
+    # complex vector due to promote_type(typeof(t), eltype(A), eltype(b)).
+    # This verifies the type promotion behavior that motivates the convert()
+    # call in the non-inplace propagator path.
+    N = 10
+    rng = StableRNG(2296342638)
+    A = random_matrix(N; rng, complex = false, hermitian = true)
+    b_real = real.(random_state_vector(N; rng))
+    @test eltype(b_real) == Float64
+    dt = 0.1
+    t = -1im * dt
+
+    # expv with complex t and real b returns Vector{ComplexF64}
+    w = ExponentialUtilities.expv(t, A, b_real)
+    @test eltype(w) == ComplexF64
+    @test typeof(w) != typeof(b_real)
+    w_expected = exp(t * A) * b_real
+    @test norm(w - w_expected) < 1e-12
+
+    # expv with complex t and complex b returns the same type as the input
+    b_complex = complex.(b_real)
+    w2 = ExponentialUtilities.expv(t, A, b_complex)
+    @test eltype(w2) == ComplexF64
+    @test typeof(w2) == typeof(b_complex)
+    @test norm(w2 - w_expected) < 1e-12
+
+end
+
+
+@testset "expv with OffsetArray state vector (broken)" begin
+
+    # expv does not support OffsetVector as state vector b: arnoldi builds
+    # V = similar(b, T, (n, m+1)) which yields an OffsetMatrix, and the
+    # subsequent broadcast in expv! throws a DimensionMismatch
+    N = 9
+    rng = StableRNG(2841259413)
+    A = random_matrix(N; rng, hermitian = true)
+    Ψ₀ = OffsetVector(random_state_vector(N; rng), (-(N÷2)):(N÷2))
+    dt = 0.1
+    t = -1im * dt
+    @test_broken begin
+        Ψ = ExponentialUtilities.expv(t, A, Ψ₀)
+        Ψ_expected = exp(t * A) * collect(Ψ₀)
+        norm(Ψ - Ψ_expected) < 1e-12
+    end
+
+end
diff --git a/test/test_exputils.jl b/test/test_exputils.jl
new file mode 100644
index 00000000..68cde686
--- /dev/null
+++ b/test/test_exputils.jl
@@ -0,0 +1,172 @@
+using Test
+using LinearAlgebra
+using SparseArrays
+using QuantumPropagators
+using QuantumPropagators.Generators: liouvillian
+using QuantumPropagators: Cheby
+using ExponentialUtilities
+using QuantumGradientGenerators: GradGenerator, GradVector
+using StableRNGs: StableRNG
+
+
+@testset "Forward propagation (Hermitian)" begin
+    # Simple TLS Rabi cycling: compare with exact matrix exponential
+    Ψ₀ = ComplexF64[1, 0]
+    Ĥ = ComplexF64[0 0.5; 0.5 0]
+    tlist = collect(range(0, 1.5π, length = 101))
+    generator = (Ĥ,)
+
+    Ψ_out = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities)
+    Ψ_expected = ComplexF64[-1/√2, -1im/√2]
+
+    @test norm(Ψ_out - Ψ_expected) < 1e-9
+end
+
+
+@testset "Forward propagation (not in-place)" begin
+    Ψ₀ = ComplexF64[1, 0]
+    Ĥ = ComplexF64[0 0.5; 0.5 0]
+    tlist = collect(range(0, 1.5π, length = 101))
+    generator = (Ĥ,)
+
+    Ψ_out = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities, inplace = false)
+    Ψ_expected = ComplexF64[-1/√2, -1im/√2]
+
+    @test norm(Ψ_out - Ψ_expected) < 1e-9
+end
+
+
+@testset "Backward propagation" begin
+    Ψ₀ = ComplexF64[1, 0]
+    Ĥ = ComplexF64[0 0.5; 0.5 0]
+    tlist = collect(range(0, 1.5π, length = 101))
+    generator = (Ĥ,)
+
+    Ψ_fw = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities)
+    Ψ_bw = propagate(Ψ_fw, generator, tlist; method = ExponentialUtilities, backward = true)
+
+    @test norm(Ψ_bw - Ψ₀) < 1e-9
+end
+
+
+@testset "Liouvillian propagation" begin
+    rng = StableRNG(20260216)
+    # 2-level system with decay (Liouvillian)
+    σ_x = ComplexF64[0 1; 1 0]
+    σ_z = ComplexF64[1 0; 0 -1]
+    γ = 0.1  # decay rate
+    c_op = √γ * ComplexF64[0 1; 0 0]  # lowering operator
+
+    H₀ = 0.5 * σ_z
+    ℒ = liouvillian(H₀, [c_op]; convention = :TDSE)
+
+    # Start from |0⟩⟨0|
+    ρ₀ = ComplexF64[1 0; 0 0]
+    ρ⃗₀ = reshape(ρ₀, :)
+
+    tlist = collect(range(0, 10.0, length = 201))
+
+    ρ⃗_out = propagate(
+        ρ⃗₀,
+        (ℒ,),
+        tlist;
+        method = ExponentialUtilities,
+        expv_kwargs = (; ishermitian = false)
+    )
+
+    ρ_out = reshape(ρ⃗_out, 2, 2)
+
+    # Density matrix must have trace 1 and be positive semidefinite
+    @test tr(ρ_out) ≈ 1.0 atol = 1e-9
+    @test all(eigvals(Hermitian(ρ_out)) .>= -1e-10)
+
+    # Compare with exact solution via matrix exponential
+    L_full = Array(ℒ)
+    U = exp(-1im * L_full * tlist[end])
+    ρ⃗_expected = U * ρ⃗₀
+    @test norm(ρ⃗_out - ρ⃗_expected) < 1e-9
+end
+
+
+@testset "Liouvillian backward propagation" begin
+    σ_z = ComplexF64[1 0; 0 -1]
+    γ = 0.01
+    c_op = √γ * ComplexF64[0 1; 0 0]
+
+    H₀ = 0.5 * σ_z
+    ℒ = liouvillian(H₀, [c_op]; convention = :TDSE)
+
+    ρ₀ = ComplexF64[1 0; 0 0]
+    ρ⃗₀ = reshape(ρ₀, :)
+
+    tlist = collect(range(0, 1.0, length = 51))
+
+    ρ⃗_fw = propagate(
+        ρ⃗₀,
+        (ℒ,),
+        tlist;
+        method = ExponentialUtilities,
+        expv_kwargs = (; ishermitian = false)
+    )
+
+    ρ⃗_bw = propagate(
+        ρ⃗_fw,
+        (ℒ,),
+        tlist;
+        method = ExponentialUtilities,
+        backward = true,
+        expv_kwargs = (; ishermitian = false)
+    )
+
+    # For non-unitary dynamics, backward propagation won't exactly recover ρ₀,
+    # but applying exp(-iLdt) backward should still be consistent with exp(+iLdt)
+    # Compare with exact matrix exponential backward propagation
+    L_full = Array(ℒ)
+    U_bw = exp(+1im * L_full * tlist[end])
+    ρ⃗_bw_expected = U_bw * ρ⃗_fw
+    @test norm(ρ⃗_bw - ρ⃗_bw_expected) < 1e-9
+end
+
+
+@testset "Time-dependent generator" begin
+    rng = StableRNG(20260216)
+    Ψ₀ = ComplexF64[1, 0]
+    H₀ = ComplexF64[0 0; 0 1]
+    H₁ = ComplexF64[0 1; 1 0]
+
+    tlist = collect(range(0, 10.0, length = 201))
+    ϵ = 0.2 * ones(length(tlist))
+
+    generator = hamiltonian(H₀, (H₁, ϵ))
+
+    Ψ_out_expv = propagate(Ψ₀, generator, tlist; method = ExponentialUtilities)
+    Ψ_out_exp = propagate(Ψ₀, generator, tlist; method = :expprop)
+
+    @test norm(Ψ_out_expv - Ψ_out_exp) < 1e-9
+end
+
+@testset "GradGen propagation" begin
+    Ψ₀ = ComplexF64[1, 0]
+    H₀ = ComplexF64[0 0; 0 1]
+    H₁ = ComplexF64[0 1; 1 0]
+
+    tlist = collect(range(0, 10.0, length = 201))
+    ϵ = 0.2 * ones(length(tlist))
+
+    generator = hamiltonian(H₀, (H₁, ϵ))
+    G̃ = GradGenerator(generator)
+    Ψ̃₀ = GradVector(Ψ₀, 1)
+
+    Ψ̃_out = propagate(
+        Ψ̃₀,
+        G̃,
+        tlist;
+        method = ExponentialUtilities,
+        inplace = false,
+        expv_kwargs = (; ishermitian = false, tol = 1e-12)
+    )
+    Ψ_cheby = propagate(Ψ₀, generator, tlist; method = Cheby)
+
+    @test norm(Ψ̃_out.state - Ψ_cheby) < 1e-7
+    @test norm(Ψ̃_out.grad_states[1]) > 0
+end
diff --git a/test/test_prop_interfaces.jl b/test/test_prop_interfaces.jl
index a6f9a733..bc36191f 100644
--- a/test/test_prop_interfaces.jl
+++ b/test/test_prop_interfaces.jl
@@ -1,9 +1,10 @@
 using Test
 
 using QuantumControlTestUtils.RandomObjects: random_state_vector, random_dynamic_generator
-using QuantumPropagators: QuantumPropagators, init_prop
+using QuantumPropagators: QuantumPropagators, ExponentialUtilitiesPropagator, init_prop
 using QuantumPropagators.Interfaces: check_propagator
 using StableRNGs: StableRNG
+using ExponentialUtilities
 using OrdinaryDiffEq: OrdinaryDiffEq
 using QuantumPropagators.Shapes: flattop
 
@@ -149,6 +150,54 @@ end
 end
 
 
+@testset "ExponentialUtilities Propagator Interface" begin
+
+    N = 10
+    rng = StableRNG(677918057)
+    tlist = collect(range(0, 10, length = 101))
+    Ψ = random_state_vector(N; rng)
+    Ĥ = random_dynamic_generator(N, tlist; rng)
+
+    propagator = init_prop(
+        Ψ,
+        Ĥ,
+        tlist;
+        method = ExponentialUtilities,
+        backward = false,
+        inplace = true,
+        verbose = false
+    )
+
+    @test propagator isa ExponentialUtilitiesPropagator
+    @test check_propagator(propagator)
+
+    propagator = init_prop(
+        Ψ,
+        Ĥ,
+        tlist;
+        method = ExponentialUtilities,
+        backward = false,
+        inplace = false,
+        verbose = false
+    )
+
+    @test check_propagator(propagator)
+
+    propagator = init_prop(
+        Ψ,
+        Ĥ,
+        tlist;
+        method = ExponentialUtilities,
+        backward = true,
+        inplace = true,
+        verbose = false
+    )
+
+    @test check_propagator(propagator)
+
+end
+
+
 @testset "ODE Propagator Interface (time-continuous)" begin
 
     N = 10
diff --git a/test/test_propagate.jl b/test/test_propagate.jl
index cb855112..f4a4bc42 100644
--- a/test/test_propagate.jl
+++ b/test/test_propagate.jl
@@ -2,6 +2,7 @@ using Test
 using QuantumPropagators
 using QuantumPropagators: Cheby
 using QuantumPropagators.Storage
+using ExponentialUtilities
 using LinearAlgebra
 using UnicodePlots
 using StaticArrays: @SMatrix, SVector, @SVector
@@ -27,12 +28,21 @@ using StaticArrays: @SMatrix, SVector, @SVector
     @test eltype(storage) == ComplexF64
 
     Ψ_out = propagate(Ψ0, generator, tlist; method = :expprop, storage = storage)
+    Ψ_out_expv = propagate(Ψ0, generator, tlist; method = ExponentialUtilities)
+    Ψ_out_expv_ni =
+        propagate(Ψ0, generator, tlist; method = ExponentialUtilities, inplace = false)
+    generator_herm = (Hermitian(Ĥ),)
+    Ψ_out_expv_herm =
+        propagate(Ψ0, generator_herm, tlist; method = ExponentialUtilities, inplace = true)
     Ψ_expected = ComplexF64[-1/√2, -1im/√2]  # note the phases
 
     pop0 = abs.(storage[1, :]) .^ 2
     SHOWPLOT && println(lineplot(tlist ./ π, pop0, ylim = [0, 1], title = "fw prop"))
 
     @test norm(Ψ_out - Ψ_expected) < 1e-12
+    @test norm(Ψ_out_expv - Ψ_expected) < 1e-9
+    @test norm(Ψ_out_expv_ni - Ψ_expected) < 1e-9
+    @test norm(Ψ_out_expv_herm - Ψ_expected) < 1e-9
     @test pop0[end] ≈ 0.5
 
     # Propagating backward in time should exactly reverse the dynamics (since