diff --git a/CLAUDE.md b/CLAUDE.md
index 85791461..60b0815a 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -5,21 +5,26 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 ## Development Commands
 
 ### Testing
+
 - `make test` - Run the full test suite
 - `julia --project=test -e 'include("test/runtests.jl")'` - Run tests directly
 
 ### Development Environment
+
 - `make devrepl` - Start interactive REPL with development environment (recommended for development)
 - `julia -i --banner=no devrepl.jl` - Alternative way to start development REPL
 
 ### Documentation
+
 - `make docs` - Build documentation
 
 ### Code Quality
+
 - `make codestyle` - Apply JuliaFormatter to entire project
 - `julia --project=test -e 'using JuliaFormatter; format(".", verbose=true)'` - Format code directly
 
 ### Cleanup
+
 - `make clean` - Clean build/doc/testing artifacts
 - `make distclean` - Restore to clean checkout state
 
@@ -28,6 +33,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 This is a Julia package for quantum system time propagation within the JuliaQuantumControl ecosystem.
 
 ### Core Structure
+
 - **src/QuantumPropagators.jl** - Main module definition with submodule includes
 - **src/interfaces/** - Interface definitions for operators, states, generators, propagators, etc.
 - **Submodules** organized by functionality:
@@ -43,27 +49,36 @@ This is a Julia package for quantum system time propagation within the JuliaQuan
   - `Generators` (generators.jl) - Hamiltonian/Liouvillian generators
 
 ### Key Propagation Methods
+
 - **Chebyshev propagation** (cheby_propagator.jl) - Polynomial expansion method
 - **Newton propagation** (newton_propagator.jl) - Newton interpolation method
 - **Matrix exponential** (exp_propagator.jl) - Direct exponentiation
 - **ODE integration** (ode_function.jl) - Interface to OrdinaryDiffEq.jl via extensions
 
 ### High-Level Interface
+
 - `propagate()` function for single time evolution
 - `propagate_sequence()` for sequential propagations with different parameters
 
 ### Extensions
+
 The package uses Julia's extension system for optional dependencies:
 - `QuantumPropagatorsODEExt` - OrdinaryDiffEq.jl integration
 - `QuantumPropagatorsRecursiveArrayToolsExt` - RecursiveArrayTools.jl support
 - `QuantumPropagatorsStaticArraysExt` - StaticArrays.jl optimization
 
 ### Testing Framework
+
 Uses SafeTestsets for isolated test execution. Tests are comprehensive and include:
 - Interface validation tests
 - Individual propagator method tests
 - Integration tests for complete propagation workflows
 
+Running `make test` prints out coverage data in a table.
+
+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
@@ -72,6 +87,7 @@ The project uses a sophisticated development setup:
 - Code formatting with JuliaFormatter
 
 ### Documentation
+
 Uses Documenter.jl with comprehensive API documentation and examples. Documentation includes detailed method explanations.
 
 ## Docstrings
diff --git a/src/QuantumPropagators.jl b/src/QuantumPropagators.jl
index 642d818b..4d5573ac 100644
--- a/src/QuantumPropagators.jl
+++ b/src/QuantumPropagators.jl
@@ -38,11 +38,13 @@ export init_prop, reinit_prop!, prop_step!
 
 #! format: off
 module Interfaces
-    export supports_inplace
+    export supports_inplace, supports_vector_interface, supports_matrix_interface
     export check_operator, check_state, check_tlist, check_amplitude
     export check_control, check_generator, check_propagator
     export check_parameterized_function, check_parameterized
     include("interfaces/supports_inplace.jl")
+    include("interfaces/supports_vector_interface.jl")
+    include("interfaces/supports_matrix_interface.jl")
     include("interfaces/utils.jl")
     include("interfaces/state.jl")
     include("interfaces/tlist.jl")
diff --git a/src/interfaces/operator.jl b/src/interfaces/operator.jl
index dd5845d4..81e0d5f7 100644
--- a/src/interfaces/operator.jl
+++ b/src/interfaces/operator.jl
@@ -17,7 +17,6 @@ verifies the given `op` relative to `state`. The `state` must pass
 An "operator" is any object that [`evaluate`](@ref) returns when evaluating a
 time-dependent dynamic generator. The specific requirements for `op` are:
 
-* `eltype(op)` must be defined and return a numeric type
 * `size(op)` must be defined and return a tuple of integers
 * `size(op, dim)` must be defined for each dimension and be consistent with
   `size(op)`
@@ -26,9 +25,11 @@ time-dependent dynamic generator. The specific requirements for `op` are:
 * `op * state` must be defined
 * The [`QuantumPropagators.Interfaces.supports_inplace`](@ref) method must be
   defined for `op`. If it returns `true`, it must be possible to evaluate a
-  generator in-place into the existing `op`. See [`check_generator`](@ref).
+  generator in-place into the existing `op`. See
+  [`QuantumPropagators.Interfaces.check_generator`](@ref).
 
-If [`QuantumPropagators.Interfaces.supports_inplace(state)`](@ref):
+If [`QuantumPropagators.Interfaces.supports_inplace(state)`](@ref
+QuantumPropagators.Interfaces.supports_inplace):
 
 * The 3-argument `LinearAlgebra.mul!` must apply `op` to the given `state`
 * The 5-argument `LinearAlgebra.mul!` must apply `op` to the given `state`
@@ -40,6 +41,25 @@ If `for_expval` (typically required for optimal control):
 * `LinearAlgebra.dot(state, op, state)` must return return a number
 * `dot(state, op, state)` must match `dot(state, op * state)`, if applicable
 
+If [`QuantumPropagators.Interfaces.supports_matrix_interface(op)`](@ref
+QuantumPropagators.Interfaces.supports_matrix_interface) is `true`, the
+operator must implement the
+[Abstract Array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+for two-dimensional arrays:
+
+* `eltype(op)` must be defined and return a numeric type
+* `getindex(op, i, j)` must be defined and return elements matching `eltype`
+* `length(op)` must equal `prod(size(op))`
+* `iterate(op)` must be defined
+* `similar(op)` must be defined and return a mutable array with the same shape
+  and element type
+* `similar(op, ::Type{S})` must return a mutable array with the same shape and
+  element type `S`
+* `similar(op, dims::Dims)` must return a mutable array with the same element
+  type and the given dimensions
+* `similar(op, ::Type{S}, dims::Dims)` must return a mutable array with the
+  given element type and dimensions
+
 The function returns `true` for a valid operator and `false` for an invalid
 operator. Unless `quiet=true`, it will log an error to indicate which of the
 conditions failed.
@@ -74,20 +94,6 @@ function check_operator(
         success = false
     end
 
-    try
-        T = eltype(op)
-        if !(T isa Type && T <: Number)
-            quiet || @error "$(px)`eltype(op)` must return a numeric type, not $T"
-            success = false
-        end
-    catch exc
-        quiet || @error(
-            "$(px)`eltype(op)` must be defined.",
-            exception = (exc, catch_abbreviated_backtrace())
-        )
-        success = false
-    end
-
     try
         s = size(op)
         if !(s isa Tuple)
@@ -171,7 +177,6 @@ function check_operator(
         success = false
     end
 
-
     if supports_inplace(state)
 
         try
@@ -254,6 +259,179 @@ function check_operator(
 
     end
 
+    if supports_matrix_interface(op)
+
+        try
+            T = eltype(op)
+            if !(T isa Type && T <: Number)
+                quiet || @error "$(px)`eltype(op)` must return a numeric type, not $T"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`eltype(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            s = size(op)
+            if length(s) == 2 && all(d -> d > 0, s)
+                val = op[1, 1]
+                T = eltype(op)
+                if T isa Type && T <: Number && !(val isa T)
+                    quiet ||
+                        @error "$(px)`op[1, 1]` must return a value of type `eltype(op)=$T`, not $(typeof(val))"
+                    success = false
+                end
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`getindex(op, i, j)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            l = length(op)
+            s = size(op)
+            if l != prod(s)
+                quiet ||
+                    @error "$(px)`length(op)` must equal `prod(size(op))`: $l ≠ $(prod(s))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`length(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            itr = iterate(op)
+            s = size(op)
+            if isnothing(itr) && prod(s) > 0
+                quiet ||
+                    @error "$(px)`iterate(op)` must not return `nothing` for a non-empty operator"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`iterate(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            op2 = similar(op)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != size(op)
+                quiet ||
+                    @error "$(px)`similar(op)` must return an array with the same shape: size $(size(op2)) ≠ $(size(op))"
+                success = false
+            end
+            if eltype(op2) != eltype(op)
+                quiet ||
+                    @error "$(px)`similar(op)` must return an array with the same element type: $(eltype(op2)) ≠ $(eltype(op))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(op) == ComplexF64) ? ComplexF32 : ComplexF64
+            op2 = similar(op, S)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != size(op)
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return an array with the same shape: size $(size(op2)) ≠ $(size(op))"
+                success = false
+            end
+            if eltype(op2) != S
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return an array with element type $S, got $(eltype(op2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, ::Type{S})` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            dims = size(op)
+            op2 = similar(op, dims)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != dims
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return an array with size $dims, got $(size(op2))"
+                success = false
+            end
+            if eltype(op2) != eltype(op)
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return an array with the same element type: $(eltype(op2)) ≠ $(eltype(op))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(op) == ComplexF64) ? ComplexF32 : ComplexF64
+            dims = size(op)
+            op2 = similar(op, S, dims)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != dims
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return an array with size $dims, got $(size(op2))"
+                success = false
+            end
+            if eltype(op2) != S
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return an array with element type $S, got $(eltype(op2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, ::Type{S}, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+    end
+
     return success
 
 end
diff --git a/src/interfaces/state.jl b/src/interfaces/state.jl
index e141aecb..f8e88f72 100644
--- a/src/interfaces/state.jl
+++ b/src/interfaces/state.jl
@@ -49,6 +49,34 @@ If `normalized` (not required by default):
 
 * `LinearAlgebra.norm(state)` must be 1
 
+If [`QuantumPropagators.Interfaces.supports_vector_interface(state)`](@ref
+QuantumPropagators.Interfaces.supports_vector_interface) is `true`, the state
+must implement the
+[Abstract Array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+for one-dimensional arrays:
+
+* `eltype(state)` must be defined and return a numeric type
+* `getindex(state, i)` must be defined and return elements matching `eltype`
+* `length(state)` must equal `prod(size(state))`
+* `iterate(state)` must be defined
+* `similar(state)` must be defined and return a mutable vector with the same
+  length and element type.
+* `similar(state, ::Type{S})` must return a mutable vector with the same length
+  and element type `S`
+* `similar(state, dims::Dims)` must return a mutable array with the same
+  element type and the given dimensions
+* `similar(state, ::Type{S}, dims::Dims)` must return a mutable array with the
+  given element type and dimensions
+
+If additionally [`QuantumPropagators.Interfaces.supports_inplace(state)`](@ref
+QuantumPropagators.Interfaces.supports_inplace) is `true` (read-write vector):
+
+* `setindex!(state, v, i)` must be defined
+
+Note: When `supports_inplace(state)` is `true`, `similar(state)` is already
+required to return the *same type* as `state` (see above). The vector interface
+checks are weaker and complementary.
+
 It is strongly recommended to always support immutable operations (also for
 mutable states)
 
@@ -357,6 +385,197 @@ function check_state(
         end
     end
 
+    if supports_vector_interface(state)
+
+        try
+            T = eltype(state)
+            if !(T isa Type && T <: Number)
+                quiet || @error "$(px)`eltype(state)` must return a numeric type, not $T"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`eltype(state)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            s = size(state)
+            if length(s) >= 1 && prod(s) > 0
+                val = state[1]
+                T = eltype(state)
+                if T isa Type && T <: Number && !(val isa T)
+                    quiet ||
+                        @error "$(px)`state[1]` must return a value of type `eltype(state)=$T`, not $(typeof(val))"
+                    success = false
+                end
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`getindex(state, i)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            l = length(state)
+            s = size(state)
+            if l != prod(s)
+                quiet ||
+                    @error "$(px)`length(state)` must equal `prod(size(state))`: $l ≠ $(prod(s))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`length(state)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            itr = iterate(state)
+            s = size(state)
+            if isnothing(itr) && prod(s) > 0
+                quiet ||
+                    @error "$(px)`iterate(state)` must not return `nothing` for a non-empty state"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`iterate(state)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            st2 = similar(state)
+            if !supports_inplace(st2)
+                quiet ||
+                    @error "$(px)`similar(state)` must return a mutable vector (`supports_inplace` must be `true`), got $(typeof(st2))"
+                success = false
+            end
+            if size(st2) != size(state)
+                quiet ||
+                    @error "$(px)`similar(state)` must return a vector with the same shape: size $(size(st2)) ≠ $(size(state))"
+                success = false
+            end
+            if eltype(st2) != eltype(state)
+                quiet ||
+                    @error "$(px)`similar(state)` must return a vector with the same element type: $(eltype(st2)) ≠ $(eltype(state))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(state)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(state) == ComplexF64) ? ComplexF32 : ComplexF64
+            st2 = similar(state, S)
+            if !supports_inplace(st2)
+                quiet ||
+                    @error "$(px)`similar(state, $S)` must return a mutable vector (`supports_inplace` must be `true`), got $(typeof(st2))"
+                success = false
+            end
+            if size(st2) != size(state)
+                quiet ||
+                    @error "$(px)`similar(state, $S)` must return a vector with the same shape: size $(size(st2)) ≠ $(size(state))"
+                success = false
+            end
+            if eltype(st2) != S
+                quiet ||
+                    @error "$(px)`similar(state, $S)` must return a vector with element type $S, got $(eltype(st2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(state, ::Type{S})` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            dims = size(state)
+            st2 = similar(state, dims)
+            if !supports_inplace(st2)
+                quiet ||
+                    @error "$(px)`similar(state, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(st2))"
+                success = false
+            end
+            if size(st2) != dims
+                quiet ||
+                    @error "$(px)`similar(state, dims)` must return an array with size $dims, got $(size(st2))"
+                success = false
+            end
+            if eltype(st2) != eltype(state)
+                quiet ||
+                    @error "$(px)`similar(state, dims)` must return an array with the same element type: $(eltype(st2)) ≠ $(eltype(state))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(state, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(state) == ComplexF64) ? ComplexF32 : ComplexF64
+            dims = size(state)
+            st2 = similar(state, S, dims)
+            if !supports_inplace(st2)
+                quiet ||
+                    @error "$(px)`similar(state, $S, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(st2))"
+                success = false
+            end
+            if size(st2) != dims
+                quiet ||
+                    @error "$(px)`similar(state, $S, dims)` must return an array with size $dims, got $(size(st2))"
+                success = false
+            end
+            if eltype(st2) != S
+                quiet ||
+                    @error "$(px)`similar(state, $S, dims)` must return an array with element type $S, got $(eltype(st2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(state, ::Type{S}, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        if supports_inplace(state)
+
+            try
+                s = size(state)
+                if prod(s) > 0
+                    v = state[1]
+                    state[1] = v  # write back the same value
+                end
+            catch exc
+                quiet || @error(
+                    "$(px)`setindex!(state, v, i)` must be defined for an in-place state.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
+
+        end
+
+    end
+
     return success
 
 end
diff --git a/src/interfaces/supports_inplace.jl b/src/interfaces/supports_inplace.jl
index 56e69d9b..10c4128a 100644
--- a/src/interfaces/supports_inplace.jl
+++ b/src/interfaces/supports_inplace.jl
@@ -1,6 +1,7 @@
 import ..Operator
 import ..ScaledOperator
 import LinearAlgebra
+import SparseArrays: SparseMatrixCSC
 
 """Indicate whether a type supports in-place operations.
 
@@ -39,6 +40,7 @@ supports_inplace(::Type{T}) where {T<:AbstractVector} = ismutabletype(T)  # fall
 supports_inplace(::Type{<:Matrix}) = true
 supports_inplace(::Type{<:Operator}) = true
 supports_inplace(::Type{<:LinearAlgebra.Diagonal}) = true
+supports_inplace(::Type{<:SparseMatrixCSC}) = true
 supports_inplace(::Type{<:ScaledOperator{<:Any,OT}}) where {OT} = supports_inplace(OT)
 supports_inplace(::Type{T}) where {T<:AbstractMatrix} = ismutabletype(T)  # fallback
 
diff --git a/src/interfaces/supports_matrix_interface.jl b/src/interfaces/supports_matrix_interface.jl
new file mode 100644
index 00000000..f09b464a
--- /dev/null
+++ b/src/interfaces/supports_matrix_interface.jl
@@ -0,0 +1,27 @@
+"""Indicate whether a type implements the 2D array interface.
+
+```julia
+supports_matrix_interface(::Type{T})
+```
+
+returns `true` if `T` implements the
+[Abstract Array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+for two-dimensional arrays. This is `true` for all subtypes of
+`AbstractMatrix`, but may also be `true` for types that implement an array
+interface (`size`, `getindex`, etc.) without declaring themselves subtypes
+of `AbstractMatrix`. Calling `supports_matrix_interface` on an instance
+`x` also works via a convenience fallback that forwards to
+`supports_matrix_interface(typeof(x))`.
+
+Depending the value of [`supports_inplace`](@ref), `T` may implement a
+read-write matrix (`setindex!` etc.) or a read-only matrix.
+
+The matrix interface is encouraged for [operators](@ref Operators), and the
+specific conditions of the required interface in that context are checked via
+[`QuantumPropagators.Interfaces.check_operator`](@ref).
+"""
+supports_matrix_interface(::Type{<:AbstractMatrix}) = true
+
+supports_matrix_interface(::Type{T}) where {T} = false
+
+supports_matrix_interface(x) = supports_matrix_interface(typeof(x))
diff --git a/src/interfaces/supports_vector_interface.jl b/src/interfaces/supports_vector_interface.jl
new file mode 100644
index 00000000..f026904c
--- /dev/null
+++ b/src/interfaces/supports_vector_interface.jl
@@ -0,0 +1,28 @@
+"""Indicate whether a type implements the 1D array interface.
+
+```julia
+supports_vector_interface(::Type{T})
+```
+
+returns `true` if `T` implements the
+[Abstract Array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+for one-dimensional arrays. This is `true` for all subtypes of
+`AbstractVector`, but may also be `true` for types that implement an array
+interface (`size`, `getindex`, etc.) without declaring themselves subtypes
+of `AbstractVector`. Calling `supports_vector_interface` on an instance
+`x` also works via a convenience fallback that forwards to
+`supports_vector_interface(typeof(x))`.
+
+Depending the value of [`supports_inplace`](@ref), `T` may implement a
+read-write vector (`setindex!` etc.) or a read-only vector.
+
+The vector interface is encouraged for quantum states, and the
+specific conditions of the required interface in that context are checked via
+[`QuantumPropagators.Interfaces.check_state`](@ref).
+
+"""
+supports_vector_interface(::Type{<:AbstractVector}) = true
+
+supports_vector_interface(::Type{T}) where {T} = false
+
+supports_vector_interface(x) = supports_vector_interface(typeof(x))
diff --git a/test/test_invalid_interfaces.jl b/test/test_invalid_interfaces.jl
index 97a522dc..05e792a7 100644
--- a/test/test_invalid_interfaces.jl
+++ b/test/test_invalid_interfaces.jl
@@ -7,6 +7,8 @@ using QuantumPropagators: QuantumPropagators, init_prop
 using LinearAlgebra
 import QuantumPropagators.Interfaces: supports_inplace
 using QuantumPropagators.Interfaces:
+    supports_matrix_interface,
+    supports_vector_interface,
     check_amplitude,
     check_control,
     check_operator,
@@ -14,6 +16,12 @@ using QuantumPropagators.Interfaces:
     check_state,
     check_propagator
 
+# Helper type: an immutable object with wrong size/eltype, used by "bad similar" tests
+struct ImmutableResult end
+QuantumPropagators.Interfaces.supports_inplace(::Type{ImmutableResult}) = false
+Base.size(::ImmutableResult) = (1,)
+Base.eltype(::Type{ImmutableResult}) = Float32
+
 
 @testset "Invalid amplitude" begin
 
@@ -102,6 +110,299 @@ end
 end
 
 
+@testset "Invalid operator with matrix interface" begin
+
+    struct InvalidMatrixOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{InvalidMatrixOp}) = true
+    QuantumPropagators.Interfaces.supports_matrix_interface(::Type{InvalidMatrixOp}) = true
+    Base.size(::InvalidMatrixOp) = (4, 4)
+    Base.size(::InvalidMatrixOp, d::Int) = 4
+    QuantumPropagators.Controls.get_controls(::InvalidMatrixOp) = ()
+    QuantumPropagators.Controls.evaluate(op::InvalidMatrixOp, args...; kwargs...) = op
+    Base.:(*)(::InvalidMatrixOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::InvalidMatrixOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::InvalidMatrixOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::InvalidMatrixOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+    operator = InvalidMatrixOp()
+
+    captured = IOCapture.capture() do
+        check_operator(operator; state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`eltype(op)` must return a numeric type")
+    @test contains(captured.output, "`getindex(op, i, j)` must be defined")
+    @test contains(captured.output, "`length(op)` must be defined")
+    @test contains(captured.output, "`iterate(op)` must be defined")
+    @test contains(captured.output, "`similar(op)` must be defined")
+    @test contains(captured.output, "`similar(op, ::Type{S})` must be defined")
+    @test contains(captured.output, "`similar(op, dims::Dims)` must be defined")
+    @test contains(captured.output, "`similar(op, ::Type{S}, dims::Dims)` must be defined")
+
+end
+
+
+@testset "Invalid operator with bad similar (matrix interface)" begin
+
+    # Operator that fully implements the matrix interface, but `similar` returns
+    # an immutable object with wrong size and eltype.
+    struct BadSimilarMatrixOp
+        data::Matrix{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSimilarMatrixOp}) = true
+    QuantumPropagators.Interfaces.supports_matrix_interface(::Type{BadSimilarMatrixOp}) =
+        true
+    Base.size(op::BadSimilarMatrixOp) = size(op.data)
+    Base.size(op::BadSimilarMatrixOp, d::Int) = size(op.data, d)
+    Base.eltype(::Type{BadSimilarMatrixOp}) = ComplexF64
+    Base.getindex(op::BadSimilarMatrixOp, i, j) = op.data[i, j]
+    Base.length(op::BadSimilarMatrixOp) = length(op.data)
+    Base.iterate(op::BadSimilarMatrixOp, args...) = iterate(op.data, args...)
+    Base.similar(::BadSimilarMatrixOp, args...) = ImmutableResult()
+    QuantumPropagators.Controls.get_controls(::BadSimilarMatrixOp) = ()
+    QuantumPropagators.Controls.evaluate(op::BadSimilarMatrixOp, args...; kwargs...) = op
+    Base.:(*)(op::BadSimilarMatrixOp, Ψ) = op.data * Ψ
+    LinearAlgebra.mul!(ϕ, op::BadSimilarMatrixOp, Ψ) = mul!(ϕ, op.data, Ψ)
+    LinearAlgebra.mul!(ϕ, op::BadSimilarMatrixOp, Ψ, α, β) = mul!(ϕ, op.data, Ψ, α, β)
+    LinearAlgebra.dot(a, op::BadSimilarMatrixOp, b) = dot(a, op.data, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+    H = zeros(ComplexF64, 4, 4)
+    H[1, 1] = 1.0
+    operator = BadSimilarMatrixOp(H)
+
+    captured = IOCapture.capture() do
+        check_operator(operator; state, tlist)
+    end
+    @test captured.value ≡ false
+    # similar(op): wrong mutability, shape, and eltype
+    @test contains(
+        captured.output,
+        "`similar(op)` must return a mutable array (`supports_inplace` must be `true`)"
+    )
+    @test contains(
+        captured.output,
+        "`similar(op)` must return an array with the same shape"
+    )
+    @test contains(
+        captured.output,
+        "`similar(op)` must return an array with the same element type"
+    )
+    # similar(op, S):
+    @test contains(captured.output, "`similar(op, ComplexF32)` must return a mutable array")
+    @test contains(
+        captured.output,
+        "`similar(op, ComplexF32)` must return an array with the same shape"
+    )
+    @test contains(
+        captured.output,
+        "`similar(op, ComplexF32)` must return an array with element type ComplexF32"
+    )
+    # similar(op, dims):
+    @test contains(captured.output, "`similar(op, dims)` must return a mutable array")
+    @test contains(captured.output, "`similar(op, dims)` must return an array with size")
+    @test contains(
+        captured.output,
+        "`similar(op, dims)` must return an array with the same element type"
+    )
+    # similar(op, S, dims):
+    @test contains(
+        captured.output,
+        "`similar(op, ComplexF32, dims)` must return a mutable array"
+    )
+    @test contains(
+        captured.output,
+        "`similar(op, ComplexF32, dims)` must return an array with size"
+    )
+    @test contains(
+        captured.output,
+        "`similar(op, ComplexF32, dims)` must return an array with element type ComplexF32"
+    )
+
+end
+
+
+@testset "Invalid operator with wrong returns" begin
+
+    # Operator where `*` returns wrong type, `mul!` returns wrong reference,
+    # and `dot` returns a mismatched value.
+    struct WrongReturnOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{WrongReturnOp}) = true
+    QuantumPropagators.Interfaces.supports_matrix_interface(::Type{WrongReturnOp}) = false
+    Base.size(::WrongReturnOp) = (4, 4)
+    Base.size(::WrongReturnOp, d::Int) = 4
+    QuantumPropagators.Controls.get_controls(::WrongReturnOp) = ()
+    QuantumPropagators.Controls.evaluate(op::WrongReturnOp, args...; kwargs...) = op
+    Base.:(*)(::WrongReturnOp, Ψ::Vector{ComplexF64}) = real.(Ψ)  # wrong type
+    function LinearAlgebra.mul!(ϕ, ::WrongReturnOp, Ψ)
+        ϕ .= 3 .* Ψ
+        return Ψ  # wrong: should return ϕ
+    end
+    function LinearAlgebra.mul!(ϕ, ::WrongReturnOp, Ψ, α, β)
+        ϕ .= 3 .* Ψ
+        return Ψ  # wrong: should return ϕ
+    end
+    LinearAlgebra.dot(a, ::WrongReturnOp, b) = 999.0 + 0im  # wrong value
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+    operator = WrongReturnOp()
+
+    captured = IOCapture.capture() do
+        check_operator(operator; state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(
+        captured.output,
+        "`op * state` must return an object of the same type as `state`"
+    )
+    @test contains(captured.output, "`mul!(ϕ, op, state)` must return the resulting ϕ")
+    @test contains(captured.output, "`mul!(ϕ, op, state)` must match `op * state`")
+    @test contains(
+        captured.output,
+        "`mul!(ϕ, op, state, α, β)` must return the resulting ϕ"
+    )
+    @test contains(
+        captured.output,
+        "`mul!(ϕ, op, state, α, β)` must match β*ϕ + α*op*state"
+    )
+    @test contains(
+        captured.output,
+        "`dot(state, op, state)` must match `dot(state, op * state)`"
+    )
+
+end
+
+
+@testset "Invalid operator with bad size" begin
+
+    struct BadSizeOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSizeOp}) = true
+    Base.size(::BadSizeOp) = 4  # not a Tuple!
+    Base.size(::BadSizeOp, d::Int) = 4
+    QuantumPropagators.Controls.get_controls(::BadSizeOp) = ()
+    QuantumPropagators.Controls.evaluate(op::BadSizeOp, args...; kwargs...) = op
+    Base.:(*)(::BadSizeOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::BadSizeOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::BadSizeOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::BadSizeOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(BadSizeOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op)` must return a tuple")
+
+    struct BadSizeOp2 end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSizeOp2}) = true
+    Base.size(::BadSizeOp2) = (1.5, 2.5)  # non-integer tuple
+    Base.size(::BadSizeOp2, d::Int) = [1.5, 2.5][d]
+    QuantumPropagators.Controls.get_controls(::BadSizeOp2) = ()
+    QuantumPropagators.Controls.evaluate(op::BadSizeOp2, args...; kwargs...) = op
+    Base.:(*)(::BadSizeOp2, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::BadSizeOp2, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::BadSizeOp2, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::BadSizeOp2, b) = dot(a, b)
+
+    captured = IOCapture.capture() do
+        check_operator(BadSizeOp2(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op)` must return a tuple of integers")
+
+end
+
+
+@testset "Invalid operator with bad size dimensions" begin
+
+    # size(op, dim) returns non-integer
+    struct BadSizeDimOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSizeDimOp}) = true
+    Base.size(::BadSizeDimOp) = (4, 4)
+    Base.size(::BadSizeDimOp, d::Int) = 4.0  # Float64, not Int
+    QuantumPropagators.Controls.get_controls(::BadSizeDimOp) = ()
+    QuantumPropagators.Controls.evaluate(op::BadSizeDimOp, args...; kwargs...) = op
+    Base.:(*)(::BadSizeDimOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::BadSizeDimOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::BadSizeDimOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::BadSizeDimOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(BadSizeDimOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op, 1)` must return an integer, not Float64")
+
+    # size(op, dim) inconsistent with size(op)
+    struct BadSizeDimOp2 end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSizeDimOp2}) = true
+    Base.size(::BadSizeDimOp2) = (4, 4)
+    Base.size(::BadSizeDimOp2, d::Int) = 5  # inconsistent: 5 ≠ 4
+    QuantumPropagators.Controls.get_controls(::BadSizeDimOp2) = ()
+    QuantumPropagators.Controls.evaluate(op::BadSizeDimOp2, args...; kwargs...) = op
+    Base.:(*)(::BadSizeDimOp2, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::BadSizeDimOp2, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::BadSizeDimOp2, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::BadSizeDimOp2, b) = dot(a, b)
+
+    captured = IOCapture.capture() do
+        check_operator(BadSizeDimOp2(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op, 1)` must be consistent with `size(op)`")
+
+end
+
+
+@testset "Invalid operator with partial matrix interface" begin
+
+    # Operator that has matrix interface but getindex returns wrong type,
+    # length is wrong, iterate returns nothing, and get_controls is non-empty.
+    struct PartialMatrixOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{PartialMatrixOp}) = false
+    QuantumPropagators.Interfaces.supports_matrix_interface(::Type{PartialMatrixOp}) = true
+    Base.size(::PartialMatrixOp) = (4, 4)
+    Base.size(::PartialMatrixOp, d::Int) = 4
+    Base.eltype(::Type{PartialMatrixOp}) = ComplexF64
+    Base.getindex(::PartialMatrixOp, i, j) = 1.0  # Float64, not ComplexF64!
+    Base.length(::PartialMatrixOp) = 99  # wrong: should be 16
+    Base.iterate(::PartialMatrixOp) = nothing  # wrong for non-empty
+    Base.similar(::PartialMatrixOp, args...) = zeros(ComplexF64, 4, 4)  # correct
+    QuantumPropagators.Controls.get_controls(::PartialMatrixOp) = (1.0,)  # non-empty!
+    QuantumPropagators.Controls.evaluate(op::PartialMatrixOp, args...; kwargs...) = op
+    Base.:(*)(::PartialMatrixOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.dot(a, ::PartialMatrixOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(PartialMatrixOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "get_controls(op) must return an empty tuple")
+    @test contains(
+        captured.output,
+        "`op[1, 1]` must return a value of type `eltype(op)=ComplexF64`"
+    )
+    @test contains(captured.output, "`length(op)` must equal `prod(size(op))`")
+    @test contains(
+        captured.output,
+        "`iterate(op)` must not return `nothing` for a non-empty operator"
+    )
+
+end
+
+
 @testset "Invalid generator" begin
 
     struct InvalidGenerator end
@@ -213,6 +514,524 @@ end
 end
 
 
+@testset "Invalid state with vector interface" begin
+
+    struct InvalidVectorState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{InvalidVectorState}) = true
+    QuantumPropagators.Interfaces.supports_vector_interface(::Type{InvalidVectorState}) =
+        true
+    Base.copy(s::InvalidVectorState) = InvalidVectorState(copy(s.data))
+    Base.similar(s::InvalidVectorState) = InvalidVectorState(similar(s.data))
+    Base.copyto!(a::InvalidVectorState, b::InvalidVectorState) =
+        (copyto!(a.data, b.data); a)
+    LinearAlgebra.dot(a::InvalidVectorState, b::InvalidVectorState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::InvalidVectorState) = norm(a.data)
+    Base.:+(a::InvalidVectorState, b::InvalidVectorState) =
+        InvalidVectorState(a.data + b.data)
+    Base.:-(a::InvalidVectorState, b::InvalidVectorState) =
+        InvalidVectorState(a.data - b.data)
+    Base.:*(α::Number, s::InvalidVectorState) = InvalidVectorState(α * s.data)
+    Base.zero(s::InvalidVectorState) = InvalidVectorState(zero(s.data))
+    Base.fill!(s::InvalidVectorState, v) = (fill!(s.data, v); s)
+    LinearAlgebra.lmul!(c, s::InvalidVectorState) = (lmul!(c, s.data); s)
+    LinearAlgebra.axpy!(c, a::InvalidVectorState, b::InvalidVectorState) =
+        (axpy!(c, a.data, b.data); b)
+    # Deliberately don't define: eltype, getindex, setindex!, length, iterate, size
+
+    state = InvalidVectorState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`eltype(state)` must return a numeric type")
+    @test contains(captured.output, "`getindex(state, i)` must be defined")
+    @test contains(captured.output, "`length(state)` must be defined")
+    @test contains(captured.output, "`iterate(state)` must be defined")
+    @test contains(captured.output, "`setindex!(state, v, i)` must be defined")
+
+end
+
+
+@testset "Invalid state with vector interface (no similar)" begin
+
+    struct InvalidVectorState2
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{InvalidVectorState2}) = false
+    QuantumPropagators.Interfaces.supports_vector_interface(::Type{InvalidVectorState2}) =
+        true
+    Base.copy(s::InvalidVectorState2) = InvalidVectorState2(copy(s.data))
+    LinearAlgebra.dot(a::InvalidVectorState2, b::InvalidVectorState2) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::InvalidVectorState2) = norm(a.data)
+    Base.:+(a::InvalidVectorState2, b::InvalidVectorState2) =
+        InvalidVectorState2(a.data + b.data)
+    Base.:-(a::InvalidVectorState2, b::InvalidVectorState2) =
+        InvalidVectorState2(a.data - b.data)
+    Base.:*(α::Number, s::InvalidVectorState2) = InvalidVectorState2(α * s.data)
+    Base.zero(s::InvalidVectorState2) = InvalidVectorState2(zero(s.data))
+    # Deliberately don't define: eltype, getindex, length, iterate, similar, size
+
+    state = InvalidVectorState2(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`eltype(state)` must return a numeric type")
+    @test contains(captured.output, "`getindex(state, i)` must be defined")
+    @test contains(captured.output, "`length(state)` must be defined")
+    @test contains(captured.output, "`iterate(state)` must be defined")
+    @test contains(captured.output, "`similar(state)` must be defined")
+    @test contains(captured.output, "`similar(state, ::Type{S})` must be defined")
+    @test contains(captured.output, "`similar(state, dims::Dims)` must be defined")
+    @test contains(
+        captured.output,
+        "`similar(state, ::Type{S}, dims::Dims)` must be defined"
+    )
+
+end
+
+
+@testset "Invalid state with bad similar (vector interface)" begin
+
+    # State that fully implements the Hilbert space and array interface, but
+    # `similar` returns an immutable object with wrong size and eltype.
+    struct BadSimilarVectorState3
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSimilarVectorState3}) = false
+    QuantumPropagators.Interfaces.supports_vector_interface(
+        ::Type{BadSimilarVectorState3}
+    ) = true
+    Base.copy(s::BadSimilarVectorState3) = BadSimilarVectorState3(copy(s.data))
+    LinearAlgebra.dot(a::BadSimilarVectorState3, b::BadSimilarVectorState3) =
+        dot(a.data, b.data)
+    LinearAlgebra.norm(a::BadSimilarVectorState3) = norm(a.data)
+    Base.:+(a::BadSimilarVectorState3, b::BadSimilarVectorState3) =
+        BadSimilarVectorState3(a.data + b.data)
+    Base.:-(a::BadSimilarVectorState3, b::BadSimilarVectorState3) =
+        BadSimilarVectorState3(a.data - b.data)
+    Base.:*(α::Number, s::BadSimilarVectorState3) = BadSimilarVectorState3(α * s.data)
+    Base.zero(s::BadSimilarVectorState3) = BadSimilarVectorState3(zero(s.data))
+    Base.eltype(::Type{BadSimilarVectorState3}) = ComplexF64
+    Base.size(s::BadSimilarVectorState3) = size(s.data)
+    Base.getindex(s::BadSimilarVectorState3, i) = s.data[i]
+    Base.length(s::BadSimilarVectorState3) = length(s.data)
+    Base.iterate(s::BadSimilarVectorState3, args...) = iterate(s.data, args...)
+    Base.similar(::BadSimilarVectorState3, args...) = ImmutableResult()
+
+    state = BadSimilarVectorState3(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    # similar(state): wrong mutability, shape, and eltype
+    @test contains(
+        captured.output,
+        "`similar(state)` must return a mutable vector (`supports_inplace` must be `true`)"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state)` must return a vector with the same shape"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state)` must return a vector with the same element type"
+    )
+    # similar(state, S):
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32)` must return a mutable vector"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32)` must return a vector with the same shape"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32)` must return a vector with element type ComplexF32"
+    )
+    # similar(state, dims):
+    @test contains(captured.output, "`similar(state, dims)` must return a mutable array")
+    @test contains(captured.output, "`similar(state, dims)` must return an array with size")
+    @test contains(
+        captured.output,
+        "`similar(state, dims)` must return an array with the same element type"
+    )
+    # similar(state, S, dims):
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32, dims)` must return a mutable array"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32, dims)` must return an array with size"
+    )
+    @test contains(
+        captured.output,
+        "`similar(state, ComplexF32, dims)` must return an array with element type"
+    )
+
+end
+
+
+@testset "Invalid state non-unit norm" begin
+
+    state = ComplexF64[2, 0, 0, 0]
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`norm(state)` must be 1, not 2.0")
+
+end
+
+
+@testset "Invalid state with partial vector interface" begin
+
+    # State that supports vector interface but getindex returns wrong type,
+    # length is wrong, and iterate returns nothing.
+    struct PartialVectorState4
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{PartialVectorState4}) = false
+    QuantumPropagators.Interfaces.supports_vector_interface(::Type{PartialVectorState4}) =
+        true
+    Base.copy(s::PartialVectorState4) = PartialVectorState4(copy(s.data))
+    LinearAlgebra.dot(a::PartialVectorState4, b::PartialVectorState4) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::PartialVectorState4) = norm(a.data)
+    Base.:+(a::PartialVectorState4, b::PartialVectorState4) =
+        PartialVectorState4(a.data + b.data)
+    Base.:-(a::PartialVectorState4, b::PartialVectorState4) =
+        PartialVectorState4(a.data - b.data)
+    Base.:*(α::Number, s::PartialVectorState4) = PartialVectorState4(α * s.data)
+    Base.zero(s::PartialVectorState4) = PartialVectorState4(zero(s.data))
+    Base.eltype(::Type{PartialVectorState4}) = ComplexF64
+    Base.size(s::PartialVectorState4) = size(s.data)
+    Base.getindex(s::PartialVectorState4, i) = real(s.data[i])  # Float64, not ComplexF64!
+    Base.length(::PartialVectorState4) = 99  # wrong: should be 4
+    Base.iterate(::PartialVectorState4) = nothing  # wrong for non-empty
+    Base.similar(::PartialVectorState4, args...) = zeros(ComplexF64, 4)
+
+    state = PartialVectorState4(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(
+        captured.output,
+        "`state[1]` must return a value of type `eltype(state)=ComplexF64`"
+    )
+    @test contains(captured.output, "`length(state)` must equal `prod(size(state))`")
+    @test contains(
+        captured.output,
+        "`iterate(state)` must not return `nothing` for a non-empty state"
+    )
+
+end
+
+
+@testset "Invalid operator with throwing size" begin
+
+    # Operator where size(op) itself throws
+    struct ThrowingSizeOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{ThrowingSizeOp}) = true
+    Base.size(::ThrowingSizeOp) = error("size not implemented")
+    QuantumPropagators.Controls.get_controls(::ThrowingSizeOp) = ()
+    QuantumPropagators.Controls.evaluate(op::ThrowingSizeOp, args...; kwargs...) = op
+    Base.:(*)(::ThrowingSizeOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::ThrowingSizeOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::ThrowingSizeOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::ThrowingSizeOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(ThrowingSizeOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op)` must be defined")
+
+end
+
+
+@testset "Invalid operator with throwing evaluate" begin
+
+    # Operator where evaluate throws
+    struct ThrowingEvalOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{ThrowingEvalOp}) = true
+    Base.size(::ThrowingEvalOp) = (4, 4)
+    Base.size(::ThrowingEvalOp, d::Int) = 4
+    QuantumPropagators.Controls.get_controls(::ThrowingEvalOp) = ()
+    QuantumPropagators.Controls.evaluate(::ThrowingEvalOp, args...; kwargs...) =
+        error("cannot evaluate")
+    Base.:(*)(::ThrowingEvalOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::ThrowingEvalOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::ThrowingEvalOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::ThrowingEvalOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(ThrowingEvalOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "op must not be time-dependent")
+
+end
+
+
+@testset "Invalid operator with throwing size(op, dim)" begin
+
+    # Operator where size(op) works but size(op, dim) throws
+    struct ThrowingSizeDimOp end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{ThrowingSizeDimOp}) = true
+    Base.size(::ThrowingSizeDimOp) = (4, 4)
+    Base.size(::ThrowingSizeDimOp, d::Int) = error("size not implemented for dim")
+    QuantumPropagators.Controls.get_controls(::ThrowingSizeDimOp) = ()
+    QuantumPropagators.Controls.evaluate(op::ThrowingSizeDimOp, args...; kwargs...) = op
+    Base.:(*)(::ThrowingSizeDimOp, Ψ::Vector{ComplexF64}) = Ψ
+    LinearAlgebra.mul!(ϕ, ::ThrowingSizeDimOp, Ψ) = copyto!(ϕ, Ψ)
+    LinearAlgebra.mul!(ϕ, ::ThrowingSizeDimOp, Ψ, α, β) = (lmul!(β, ϕ); axpy!(α, Ψ, ϕ))
+    LinearAlgebra.dot(a, ::ThrowingSizeDimOp, b) = dot(a, b)
+
+    state = ComplexF64[1, 0, 0, 0]
+    tlist = collect(range(0, 10, length = 101))
+
+    captured = IOCapture.capture() do
+        check_operator(ThrowingSizeDimOp(); state, tlist)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`size(op, 1)` must be defined")
+
+end
+
+
+@testset "Invalid state with wrong copy type" begin
+
+    # State where copy returns the wrong type (inner Vector instead of wrapper)
+    struct WrongCopyState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{WrongCopyState}) = true
+    Base.copy(s::WrongCopyState) = copy(s.data)  # wrong: returns Vector, not WrongCopyState
+    Base.similar(s::WrongCopyState) = WrongCopyState(similar(s.data))
+    Base.copyto!(a::WrongCopyState, b::WrongCopyState) = (copyto!(a.data, b.data); a)
+    LinearAlgebra.dot(a::WrongCopyState, b::WrongCopyState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::WrongCopyState) = norm(a.data)
+    Base.:+(a::WrongCopyState, b::WrongCopyState) = WrongCopyState(a.data + b.data)
+    Base.:-(a::WrongCopyState, b::WrongCopyState) = WrongCopyState(a.data - b.data)
+    Base.:*(α::Number, s::WrongCopyState) = WrongCopyState(α * s.data)
+    Base.zero(s::WrongCopyState) = WrongCopyState(zero(s.data))
+    Base.fill!(s::WrongCopyState, v) = (fill!(s.data, v); s)
+    LinearAlgebra.lmul!(c, s::WrongCopyState) = (lmul!(c, s.data); s)
+    LinearAlgebra.axpy!(c, a::WrongCopyState, b::WrongCopyState) =
+        (axpy!(c, a.data, b.data); b)
+
+    state = WrongCopyState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(
+        captured.output,
+        "`copy(state)::Vector{ComplexF64}` must have the same type as `state::"
+    )
+
+end
+
+
+@testset "Invalid state with constant norm" begin
+
+    # State where norm always returns a constant (inconsistent with dot)
+    struct ConstantNormState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{ConstantNormState}) = false
+    Base.copy(s::ConstantNormState) = ConstantNormState(copy(s.data))
+    LinearAlgebra.dot(a::ConstantNormState, b::ConstantNormState) = dot(a.data, b.data)
+    LinearAlgebra.norm(::ConstantNormState) = 999.0  # always 999
+    Base.:+(a::ConstantNormState, b::ConstantNormState) = ConstantNormState(a.data + b.data)
+    Base.:-(a::ConstantNormState, b::ConstantNormState) = ConstantNormState(a.data - b.data)
+    Base.:*(α::Number, s::ConstantNormState) = ConstantNormState(α * s.data)
+    Base.zero(s::ConstantNormState) = ConstantNormState(zero(s.data))
+
+    state = ConstantNormState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "must match `√(state⋅state)")
+    @test contains(captured.output, "`state - state` must have norm 0")
+
+end
+
+
+@testset "Invalid state with squared norm" begin
+
+    # State where norm = sum(abs2) which violates triangle inequality
+    struct SquaredNormState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{SquaredNormState}) = false
+    Base.copy(s::SquaredNormState) = SquaredNormState(copy(s.data))
+    LinearAlgebra.dot(a::SquaredNormState, b::SquaredNormState) = dot(a.data, b.data)
+    LinearAlgebra.norm(s::SquaredNormState) = sum(abs2, s.data)
+    Base.:+(a::SquaredNormState, b::SquaredNormState) = SquaredNormState(a.data + b.data)
+    Base.:-(a::SquaredNormState, b::SquaredNormState) = SquaredNormState(a.data - b.data)
+    Base.:*(α::Number, s::SquaredNormState) = SquaredNormState(α * s.data)
+    Base.zero(s::SquaredNormState) = SquaredNormState(zero(s.data))
+
+    state = SquaredNormState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "must fulfill the triangle inequality")
+
+end
+
+
+@testset "Invalid state with similar returning wrong type" begin
+
+    # State where similar returns a plain Vector (wrong type) in inplace section
+    struct BadSimilarTypeState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadSimilarTypeState}) = true
+    Base.copy(s::BadSimilarTypeState) = BadSimilarTypeState(copy(s.data))
+    Base.similar(s::BadSimilarTypeState) = similar(s.data)  # returns Vector, not wrapper
+    Base.copyto!(a::BadSimilarTypeState, b::BadSimilarTypeState) =
+        (copyto!(a.data, b.data); a)
+    LinearAlgebra.dot(a::BadSimilarTypeState, b::BadSimilarTypeState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::BadSimilarTypeState) = norm(a.data)
+    Base.:+(a::BadSimilarTypeState, b::BadSimilarTypeState) =
+        BadSimilarTypeState(a.data + b.data)
+    Base.:-(a::BadSimilarTypeState, b::BadSimilarTypeState) =
+        BadSimilarTypeState(a.data - b.data)
+    Base.:*(α::Number, s::BadSimilarTypeState) = BadSimilarTypeState(α * s.data)
+    Base.zero(s::BadSimilarTypeState) = BadSimilarTypeState(zero(s.data))
+    Base.fill!(s::BadSimilarTypeState, v) = (fill!(s.data, v); s)
+    LinearAlgebra.lmul!(c, s::BadSimilarTypeState) = (lmul!(c, s.data); s)
+    LinearAlgebra.axpy!(c, a::BadSimilarTypeState, b::BadSimilarTypeState) =
+        (axpy!(c, a.data, b.data); b)
+
+    state = BadSimilarTypeState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(
+        captured.output,
+        "`similar(state)::Vector{ComplexF64}` must have the same type as `state::"
+    )
+
+end
+
+
+@testset "Invalid state with unfaithful copyto!" begin
+
+    # State where copyto! doesn't actually copy the data
+    struct BadCopytoState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadCopytoState}) = true
+    Base.copy(s::BadCopytoState) = BadCopytoState(copy(s.data))
+    Base.similar(s::BadCopytoState) = BadCopytoState(similar(s.data))
+    Base.copyto!(a::BadCopytoState, b::BadCopytoState) = a  # doesn't actually copy!
+    LinearAlgebra.dot(a::BadCopytoState, b::BadCopytoState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::BadCopytoState) = norm(a.data)
+    Base.:+(a::BadCopytoState, b::BadCopytoState) = BadCopytoState(a.data + b.data)
+    Base.:-(a::BadCopytoState, b::BadCopytoState) = BadCopytoState(a.data - b.data)
+    Base.:*(α::Number, s::BadCopytoState) = BadCopytoState(α * s.data)
+    Base.zero(s::BadCopytoState) = BadCopytoState(zero(s.data))
+    Base.fill!(s::BadCopytoState, v) = (fill!(s.data, v); s)
+    LinearAlgebra.lmul!(c, s::BadCopytoState) = (lmul!(c, s.data); s)
+    LinearAlgebra.axpy!(c, a::BadCopytoState, b::BadCopytoState) =
+        (axpy!(c, a.data, b.data); b)
+
+    state = BadCopytoState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(
+        captured.output,
+        "must have norm 0, where `ϕ = similar(state); copyto!(ϕ, state)`"
+    )
+
+end
+
+
+@testset "Invalid state with fill! not zeroing" begin
+
+    # State where fill! returns the state (correct type) but doesn't zero the data
+    struct BadFillNormState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BadFillNormState}) = true
+    Base.copy(s::BadFillNormState) = BadFillNormState(copy(s.data))
+    Base.similar(s::BadFillNormState) = BadFillNormState(similar(s.data))
+    Base.copyto!(a::BadFillNormState, b::BadFillNormState) = (copyto!(a.data, b.data); a)
+    LinearAlgebra.dot(a::BadFillNormState, b::BadFillNormState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::BadFillNormState) = norm(a.data)
+    Base.:+(a::BadFillNormState, b::BadFillNormState) = BadFillNormState(a.data + b.data)
+    Base.:-(a::BadFillNormState, b::BadFillNormState) = BadFillNormState(a.data - b.data)
+    Base.:*(α::Number, s::BadFillNormState) = BadFillNormState(α * s.data)
+    Base.zero(s::BadFillNormState) = BadFillNormState(zero(s.data))
+    Base.fill!(s::BadFillNormState, v) = s  # returns state but doesn't fill!
+    LinearAlgebra.lmul!(c, s::BadFillNormState) = (lmul!(c, s.data); s)
+    LinearAlgebra.axpy!(c, a::BadFillNormState, b::BadFillNormState) =
+        (axpy!(c, a.data, b.data); b)
+
+    state = BadFillNormState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`fill!(state, 0.0)` must have norm 0")
+
+end
+
+
+@testset "Invalid state with buggy inplace operations" begin
+
+    # State where fill!/lmul!/axpy! are implemented but produce wrong results
+    struct BuggyInplaceState
+        data::Vector{ComplexF64}
+    end
+    QuantumPropagators.Interfaces.supports_inplace(::Type{BuggyInplaceState}) = true
+    Base.copy(s::BuggyInplaceState) = BuggyInplaceState(copy(s.data))
+    Base.similar(s::BuggyInplaceState) = BuggyInplaceState(similar(s.data))
+    Base.copyto!(a::BuggyInplaceState, b::BuggyInplaceState) = (copyto!(a.data, b.data); a)
+    LinearAlgebra.dot(a::BuggyInplaceState, b::BuggyInplaceState) = dot(a.data, b.data)
+    LinearAlgebra.norm(a::BuggyInplaceState) = norm(a.data)
+    Base.:+(a::BuggyInplaceState, b::BuggyInplaceState) = BuggyInplaceState(a.data + b.data)
+    Base.:-(a::BuggyInplaceState, b::BuggyInplaceState) = BuggyInplaceState(a.data - b.data)
+    Base.:*(α::Number, s::BuggyInplaceState) = BuggyInplaceState(α * s.data)
+    Base.zero(s::BuggyInplaceState) = BuggyInplaceState(zero(s.data))
+    # fill! returns wrong type (nothing instead of the state)
+    Base.fill!(s::BuggyInplaceState, v) = (fill!(s.data, v); nothing)
+    # lmul! doesn't actually scale (just returns the state unchanged)
+    LinearAlgebra.lmul!(c, s::BuggyInplaceState) = s
+    # axpy! doesn't actually add (just returns the destination unchanged)
+    LinearAlgebra.axpy!(c, a::BuggyInplaceState, b::BuggyInplaceState) = b
+
+    state = BuggyInplaceState(ComplexF64[1, 0, 0, 0])
+    captured = IOCapture.capture() do
+        check_state(state; normalized = true)
+    end
+    @test captured.value ≡ false
+    @test contains(captured.output, "`fill!(state, 0.0)` must return the filled state")
+    @test contains(captured.output, "`norm(state)` must have absolute homogeneity")
+    @test contains(captured.output, "`axpy!(a, state, ϕ)` must match `ϕ += a * state`")
+
+end
+
+
 @testset "Invalid propagator" begin
 
     struct InvalidPropagatorEmpty end