From 5b049af00dab0ad5404c1446ff6113b56c829a90 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Fri, 9 Jan 2026 18:55:08 +0000
Subject: [PATCH 01/19] Working on propagator struct implementation

---
 src/time_evolution/propagator.jl | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 src/time_evolution/propagator.jl

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
new file mode 100644
index 000000000..e69de29bb

From 8023dbcda56f55f4be58c88867ae62a57a1be46b Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Wed, 11 Feb 2026 03:41:37 +0000
Subject: [PATCH 02/19] Added propagator.jl and associated tests to
 test/core-test/propagator.jl.

---
 src/QuantumToolbox.jl            |   1 +
 src/time_evolution/propagator.jl | 303 +++++++++++++++++++++++++++++++
 test/core-test/propagator.jl     |  21 +++
 3 files changed, 325 insertions(+)

diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index dc5d17c13..1dd8e5034 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -124,6 +124,7 @@ include("time_evolution/mcsolve.jl")
 include("time_evolution/ssesolve.jl")
 include("time_evolution/smesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
+include("time_evolution/propagator.jl")
 
 ## Other functionalities
 include("correlations.jl")
diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index e69de29bb..9a2f48682 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -0,0 +1,303 @@
+export Propagator, propagator
+@doc raw"""
+    Propagator{HT, PT, DT, KWT}
+
+A callable struct representing a time-evolution propagator for a quantum system.
+
+It lazily computes and caches propagators over requested time intervals. For time-independent Hamiltonians
+([`QuantumObject`](@ref)), the propagator is computed via matrix exponentiation:
+
+```math
+\hat{U}(t, t_0) = e^{-i \hat{H} (t - t_0)}
+```
+
+For time-dependent Hamiltonians ([`QuantumObjectEvolution`](@ref)), the propagator is obtained by solving the
+Schrödinger equation ([`sesolve`](@ref)) for [`Operator`](@ref) types, or the master equation ([`mesolve`](@ref))
+for [`SuperOperator`](@ref) types, using an identity matrix as the initial condition.
+
+Previously computed propagators for sub-intervals are reused automatically to avoid redundant computation.
+
+# Fields
+
+- `H`: The Hamiltonian or Liouvillian of the system.
+- `props`: A dictionary mapping time intervals `[t0, t]` to their computed propagators.
+- `dims`: The dimensions of the Hilbert space.
+- `solver_kwargs`: Keyword arguments forwarded to the underlying solver ([`sesolve`](@ref) or [`mesolve`](@ref)).
+- `max_saved`: Maximum number of propagators to cache.
+- `threshhold`: Numerical tolerance for matching stored time intervals.
+- `remember_by_default`: Whether to cache newly computed propagators by default.
+
+# Usage
+
+A `Propagator` is callable. Use `U(t; t0=0.0)` to obtain the propagator from `t0` to `t`, or `U([t0, t])` for
+an interval. See [`propagator`](@ref) for construction.
+"""
+struct Propagator{
+    HT<:Union{Operator, SuperOperator},
+    PT<:AbstractQuantumObject,
+    DT<:AbstractDimensions,
+    KWT
+}   
+    H::AbstractQuantumObject{HT}
+    props::Dict{Vector, PT}
+    dims::AbstractArray
+    dimensions::DT
+    solver_kwargs::KWT
+    max_saved::Union{Integer, Float64}
+    threshhold::Float64
+    remember_by_default::Bool
+end
+
+
+@doc raw"""
+    propagator(
+        H::AbstractQuantumObject{HOpType},
+        t::Union{Nothing, Real} = nothing;
+        t0 = 0.0,
+        threshhold::Float64 = 1e-9,
+        max_saved::Union{Integer, Float64} = typemax(Int),
+        remember_by_default::Bool = true,
+        params = NullParameters(),
+        progress_bar::Union{Val, Bool} = Val(true),
+        inplace::Union{Val, Bool} = Val(true),
+        kwargs...,
+    )
+
+Construct a [`Propagator`](@ref) object for the quantum system described by Hamiltonian (or Liouvillian) `H`.
+
+If `t` is provided, the propagator from `t0` to `t` is immediately computed and cached. Otherwise, an empty
+`Propagator` is returned, ready to be evaluated lazily at arbitrary times.
+
+# Arguments
+
+- `H`: Hamiltonian or Liouvillian of the system ``\hat{H}``. It can be a [`QuantumObject`](@ref) or a
+  [`QuantumObjectEvolution`](@ref), with type [`Operator`](@ref) or [`SuperOperator`](@ref).
+- `t`: Optional final time. If given, the propagator for the interval `[t0, t]` is computed immediately.
+- `t0`: Initial time. Default is `0.0`.
+- `threshhold`: Numerical tolerance for matching cached time intervals. Default is `1e-9`.
+- `max_saved`: Maximum number of propagators to store in the cache. Can be an `Integer` or `Inf`. Default is `typemax(Int)`.
+- `remember_by_default`: Whether to automatically cache computed propagators. Default is `true`.
+- `params`: Parameters to pass to the underlying solver.
+- `progress_bar`: Whether to show a progress bar during time evolution. Using non-`Val` types might lead to type instabilities.
+- `inplace`: Whether to use inplace operations for the ODE solver. Default is `Val(true)`.
+- `kwargs`: Additional keyword arguments forwarded to the solver ([`sesolve`](@ref) or [`mesolve`](@ref)).
+
+# Returns
+
+- `U::Propagator`: A callable propagator object. Use `U(t; t0=0.0)` or `U([t0, t])` to evaluate.
+
+# Example
+
+```julia
+H = (ϵ / 2) * sigmaz() + (Ω / 2) * sigmax()
+U = propagator(H)      # create lazy propagator
+ψt = U(π) * basis(2, 0) # propagate |0⟩ to time π
+```
+"""
+function propagator(
+    H::AbstractQuantumObject{HOpType},
+    t::Union{Nothing, Real} = nothing;
+    t0 = 0.0,
+    threshhold::Float64=1e-9,
+    max_saved::Union{Integer, Float64}=typemax(Int),
+    remember_by_default::Bool=true,
+    params = NullParameters(),
+    progress_bar::Union{Val, Bool}=Val(true),
+    inplace::Union{Val, Bool}=Val(true),
+    kwargs...,
+) where HOpType<:Union{Operator, SuperOperator}
+
+    full_kwargs = (; params, progress_bar, inplace, kwargs...)
+
+    if !(max_saved isa Integer) && max_saved != Inf
+        max_saved = ceil(Int, max_saved)
+        @warn "max_saved should be an Integer or Inf. Setting to $max_saved."
+    end
+
+    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshhold, remember_by_default)
+
+    if t != nothing
+        U(t; t0 = t0, remember = true)
+    end
+    return U
+end
+
+
+@doc raw"""
+    (U::Propagator)(t; t0 = 0.0, remember = nothing, return_result = true, save_steps = true)
+
+Evaluate the propagator `U` from time `t0` to time `t`.
+
+This computes the time-evolution propagator ``\hat{U}(t, t_0)`` by combining any previously cached
+sub-interval propagators with newly computed ones for uncovered gaps. The full interval `[t0, t]`
+propagator is also cached separately when `remember` is enabled.
+
+# Arguments
+
+- `t`: The final time.
+- `t0`: The initial time. Default is `0.0`.
+- `remember`: Whether to cache newly computed propagators. If `nothing` (default), caching follows the
+  `remember_by_default` setting of the [`Propagator`](@ref), subject to the `max_saved` limit.
+- `return_result`: Whether to return the computed propagator. Default is `true`.
+- `save_steps`: Whether to cache intermediate sub-interval propagators in addition to the full `[t0, t]`
+  interval. Default is `true`. Set to `false` to only cache the composite result.
+
+# Returns
+
+- The propagator as an `AbstractQuantumObject` (if `return_result` is `true`).
+"""
+function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, return_result = true, save_steps = true)
+    intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold=U.threshhold)
+
+    prop = qeye_like(U.H)
+    if prop isa QobjEvo
+        prop = prop(0.0)
+    end
+
+    all_intervals = vcat(intervals.usable, intervals.to_compute)
+    sort!(all_intervals, by = first)
+
+    for interval in all_intervals
+        temp_prop = _propagator_compute_or_look_up(U, interval)
+        if (remember === nothing ? (U.remember_by_default && length(U.props)<U.max_saved) : remember) && !(interval in keys(U.props)) && save_steps
+            if length(U.props) >= U.max_saved
+                @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
+            end
+            U.props[interval] = temp_prop
+        end
+        prop = temp_prop * prop
+    end
+
+    if (remember === nothing ? (U.remember_by_default && length(U.props)<U.max_saved) : remember) && !([t0, t] in keys(U.props))
+        if length(U.props) >= U.max_saved
+            @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
+        end
+        U.props[[t0, t]] = prop
+    end
+
+    if return_result
+        return prop
+    end
+end
+
+@doc raw"""
+    (U::Propagator)(interval::Vector; kwargs...)
+
+Evaluate the propagator `U` over the time interval `[interval[1], interval[2]]`.
+
+This is a convenience method equivalent to `U(interval[2]; t0 = interval[1], kwargs...)`.
+See the other callable method of [`Propagator`](@ref) for details on keyword arguments.
+"""
+function (U::Propagator)(interval::Vector; kwargs...)
+    U(interval[2]; t0 = interval[1], kwargs...)
+end
+
+
+"""
+    _propagator_compute_or_look_up(U::Propagator{HT}, interval) where HT
+
+Look up a cached propagator for the given `interval`, or compute it if not found.
+
+For time-independent Hamiltonians ([`QuantumObject`](@ref)), the interval is shifted to `[0, Δt]` since the
+propagator depends only on the duration. For time-dependent Hamiltonians ([`QuantumObjectEvolution`](@ref)),
+the propagator is computed via [`sesolve`](@ref) (for [`Operator`](@ref)) or [`mesolve`](@ref) (for
+[`SuperOperator`](@ref)) using an identity matrix as the initial state.
+"""
+function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where HT<:Union{Operator, SuperOperator}
+    if !(U.H isa QobjEvo)
+        interval = [0.0, interval[2]-interval[1]]
+    end
+
+    if interval in keys(U.props)
+        return U.props[interval]
+    else
+        if !(U.H isa QobjEvo)
+            println("Computing propagator for interval [$(interval[1]), $(interval[2])] via matrix exponentiation...")
+            if HT<:Operator
+                return exp(-1im * U.H*(interval[2]-interval[1]))
+            else
+                return exp(U.H*(interval[2]-interval[1]))
+            end
+        end
+
+        if HT <: Operator
+            return sesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]     
+        else
+            return mesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]     
+        end
+    end
+end
+
+
+"""
+    _get_intervals_for_range(stored_intervals, target_interval; threshold=1e-9)
+
+Decompose `target_interval = [a, b]` into sub-intervals by reusing `stored_intervals` where possible.
+
+Returns a `NamedTuple` with:
+- `usable`: Stored intervals fully contained within `[a, b]` (within `threshold` tolerance).
+- `to_compute`: Gap intervals not covered by any stored interval that still need to be computed.
+"""
+function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_interval::Vector; threshold=1e-9) where T<:Vector
+    a, b = target_interval
+    
+    # Find stored intervals that are fully contained within target range (with fuzzy boundaries)
+    usable = [[s, e] for (s, e) in stored_intervals if s >= a - threshold && e <= b + threshold]
+    sort!(usable, by = first)
+    
+    # Merge usable intervals to find coverage
+    merged = Vector[]
+    for (s, e) in usable
+        if isempty(merged) || s > merged[end][2] + threshold
+            push!(merged, [s, e])
+        else
+            merged[end] = (merged[end][1], max(merged[end][2], e))
+        end
+    end
+    
+    # Find gaps that need to be computed
+    to_compute = Vector[]
+    current = a
+    for (s, e) in merged
+        if current < s - threshold
+            push!(to_compute, [current, s])
+        end
+        current = max(current, e)
+    end
+    if current < b - threshold
+        push!(to_compute, [current, b])
+    end
+    return (usable = usable, to_compute = to_compute)
+end
+
+
+function Base.show(io::IO, U::Propagator)
+    saved_times = String[]
+    times = collect(keys(U.props))
+    for i in 1:length(times)
+        push!(saved_times, "\n  $(times[i][1]) -> $(times[i][2])")
+    end
+    if length(saved_times) == 0
+        saved_times = ["None"]
+    end
+    println(
+    io,
+    "\nPropagator: H Type=",
+    U.H.type,
+    "   dims=",
+    _get_dims_string(U.dimensions),
+    "   size=",
+    size(U),
+    "\nSaved Propagators: ",
+    saved_times...
+    )
+end
+
+function Base.size(U::Propagator)
+    return size(U.H)
+end
+
+function Base.length(U::Propagator)
+    return length(U.H)
+end
+
diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl
index 4774d509b..39c65496a 100644
--- a/test/core-test/propagator.jl
+++ b/test/core-test/propagator.jl
@@ -18,3 +18,24 @@
         @test isapprox(U_me^n * ρ0, ρt[n]; atol = 1.0e-5)
     end
 end
+
+@testitem "Propagator (by propagator)" begin
+    ϵ0 = 1.0 * 2π
+    Ω = 0.8 * 2π
+    H = (ϵ0 / 2) * sigmaz() + (Ω / 2) * sigmax()
+    L = liouvillian(H)
+    ψ0 = basis(2, 0)
+    ρ0 = mat2vec(ket2dm(ψ0))
+
+    Δt = π / 5
+    tlist = 0:Δt:(2π)
+    ψt = sesolve(H, ψ0, tlist; progress_bar = Val(false)).states[2:end] # ignore the initial state
+    ρt = mesolve(H, ρ0, tlist; progress_bar = Val(false)).states[2:end] # ignore the initial state
+    U_se = propagator(H)(Δt)
+    U_me = propagator(L)(Δt)
+
+    for n in 1:(length(tlist) - 1)
+        @test isapprox(U_se^n * ψ0, ψt[n]; atol = 1.0e-5)
+        @test isapprox(U_me^n * ρ0, ρt[n]; atol = 1.0e-5)
+    end
+end
\ No newline at end of file

From 2f2478a1eb076363aab2774ae9e990a9923a4475 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 06:23:18 +0000
Subject: [PATCH 03/19] added `period` to Propagator.

---
 src/time_evolution/propagator.jl | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 9a2f48682..002b3e29d 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -26,6 +26,7 @@ Previously computed propagators for sub-intervals are reused automatically to av
 - `max_saved`: Maximum number of propagators to cache.
 - `threshhold`: Numerical tolerance for matching stored time intervals.
 - `remember_by_default`: Whether to cache newly computed propagators by default.
+- `period`: If the system is periodic in time, the propagator will calculate intervals modulo this period.
 
 # Usage
 
@@ -46,6 +47,7 @@ struct Propagator{
     max_saved::Union{Integer, Float64}
     threshhold::Float64
     remember_by_default::Bool
+    period::Real
 end
 
 
@@ -75,6 +77,7 @@ If `t` is provided, the propagator from `t0` to `t` is immediately computed and
 - `t`: Optional final time. If given, the propagator for the interval `[t0, t]` is computed immediately.
 - `t0`: Initial time. Default is `0.0`.
 - `threshhold`: Numerical tolerance for matching cached time intervals. Default is `1e-9`.
+- `period`: If the system is periodic in time, specify the period to automatically calculate intervals modulo this period. Default is `Inf` (no periodicity).
 - `max_saved`: Maximum number of propagators to store in the cache. Can be an `Integer` or `Inf`. Default is `typemax(Int)`.
 - `remember_by_default`: Whether to automatically cache computed propagators. Default is `true`.
 - `params`: Parameters to pass to the underlying solver.
@@ -99,6 +102,7 @@ function propagator(
     t::Union{Nothing, Real} = nothing;
     t0 = 0.0,
     threshhold::Float64=1e-9,
+    period::Real = Inf,
     max_saved::Union{Integer, Float64}=typemax(Int),
     remember_by_default::Bool=true,
     params = NullParameters(),
@@ -114,7 +118,7 @@ function propagator(
         @warn "max_saved should be an Integer or Inf. Setting to $max_saved."
     end
 
-    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshhold, remember_by_default)
+    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshhold, remember_by_default, period)
 
     if t != nothing
         U(t; t0 = t0, remember = true)
@@ -147,6 +151,10 @@ propagator is also cached separately when `remember` is enabled.
 - The propagator as an `AbstractQuantumObject` (if `return_result` is `true`).
 """
 function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, return_result = true, save_steps = true)
+    if U.period != Inf
+        t = mod(t, U.period)
+        t0 = mod(t0, U.period)
+    end
     intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold=U.threshhold)
 
     prop = qeye_like(U.H)
@@ -180,14 +188,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, re
     end
 end
 
-@doc raw"""
-    (U::Propagator)(interval::Vector; kwargs...)
-
-Evaluate the propagator `U` over the time interval `[interval[1], interval[2]]`.
 
-This is a convenience method equivalent to `U(interval[2]; t0 = interval[1], kwargs...)`.
-See the other callable method of [`Propagator`](@ref) for details on keyword arguments.
-"""
 function (U::Propagator)(interval::Vector; kwargs...)
     U(interval[2]; t0 = interval[1], kwargs...)
 end

From b72ac0a1938cd4ebb38a63ac3bba46f250b87fd4 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 06:30:49 +0000
Subject: [PATCH 04/19] ran format

---
 src/time_evolution/propagator.jl | 87 ++++++++++++++++----------------
 test/core-test/propagator.jl     |  2 +-
 2 files changed, 44 insertions(+), 45 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 002b3e29d..0fde19d6d 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -34,11 +34,11 @@ A `Propagator` is callable. Use `U(t; t0=0.0)` to obtain the propagator from `t0
 an interval. See [`propagator`](@ref) for construction.
 """
 struct Propagator{
-    HT<:Union{Operator, SuperOperator},
-    PT<:AbstractQuantumObject,
-    DT<:AbstractDimensions,
-    KWT
-}   
+        HT <: Union{Operator, SuperOperator},
+        PT <: AbstractQuantumObject,
+        DT <: AbstractDimensions,
+        KWT,
+    }
     H::AbstractQuantumObject{HT}
     props::Dict{Vector, PT}
     dims::AbstractArray
@@ -98,18 +98,18 @@ U = propagator(H)      # create lazy propagator
 ```
 """
 function propagator(
-    H::AbstractQuantumObject{HOpType},
-    t::Union{Nothing, Real} = nothing;
-    t0 = 0.0,
-    threshhold::Float64=1e-9,
-    period::Real = Inf,
-    max_saved::Union{Integer, Float64}=typemax(Int),
-    remember_by_default::Bool=true,
-    params = NullParameters(),
-    progress_bar::Union{Val, Bool}=Val(true),
-    inplace::Union{Val, Bool}=Val(true),
-    kwargs...,
-) where HOpType<:Union{Operator, SuperOperator}
+        H::AbstractQuantumObject{HOpType},
+        t::Union{Nothing, Real} = nothing;
+        t0 = 0.0,
+        threshhold::Float64 = 1.0e-9,
+        period::Real = Inf,
+        max_saved::Union{Integer, Float64} = typemax(Int),
+        remember_by_default::Bool = true,
+        params = NullParameters(),
+        progress_bar::Union{Val, Bool} = Val(true),
+        inplace::Union{Val, Bool} = Val(true),
+        kwargs...,
+    ) where {HOpType <: Union{Operator, SuperOperator}}
 
     full_kwargs = (; params, progress_bar, inplace, kwargs...)
 
@@ -150,12 +150,12 @@ propagator is also cached separately when `remember` is enabled.
 
 - The propagator as an `AbstractQuantumObject` (if `return_result` is `true`).
 """
-function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, return_result = true, save_steps = true)
+function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing, return_result = true, save_steps = true)
     if U.period != Inf
         t = mod(t, U.period)
         t0 = mod(t0, U.period)
     end
-    intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold=U.threshhold)
+    intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold = U.threshhold)
 
     prop = qeye_like(U.H)
     if prop isa QobjEvo
@@ -167,7 +167,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, re
 
     for interval in all_intervals
         temp_prop = _propagator_compute_or_look_up(U, interval)
-        if (remember === nothing ? (U.remember_by_default && length(U.props)<U.max_saved) : remember) && !(interval in keys(U.props)) && save_steps
+        if (remember === nothing ? (U.remember_by_default && length(U.props) < U.max_saved) : remember) && !(interval in keys(U.props)) && save_steps
             if length(U.props) >= U.max_saved
                 @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
             end
@@ -176,7 +176,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool}=nothing, re
         prop = temp_prop * prop
     end
 
-    if (remember === nothing ? (U.remember_by_default && length(U.props)<U.max_saved) : remember) && !([t0, t] in keys(U.props))
+    if (remember === nothing ? (U.remember_by_default && length(U.props) < U.max_saved) : remember) && !([t0, t] in keys(U.props))
         if length(U.props) >= U.max_saved
             @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
         end
@@ -190,7 +190,7 @@ end
 
 
 function (U::Propagator)(interval::Vector; kwargs...)
-    U(interval[2]; t0 = interval[1], kwargs...)
+    return U(interval[2]; t0 = interval[1], kwargs...)
 end
 
 
@@ -204,9 +204,9 @@ propagator depends only on the duration. For time-dependent Hamiltonians ([`Quan
 the propagator is computed via [`sesolve`](@ref) (for [`Operator`](@ref)) or [`mesolve`](@ref) (for
 [`SuperOperator`](@ref)) using an identity matrix as the initial state.
 """
-function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where HT<:Union{Operator, SuperOperator}
+function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <: Union{Operator, SuperOperator}}
     if !(U.H isa QobjEvo)
-        interval = [0.0, interval[2]-interval[1]]
+        interval = [0.0, interval[2] - interval[1]]
     end
 
     if interval in keys(U.props)
@@ -214,17 +214,17 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where HT<:U
     else
         if !(U.H isa QobjEvo)
             println("Computing propagator for interval [$(interval[1]), $(interval[2])] via matrix exponentiation...")
-            if HT<:Operator
-                return exp(-1im * U.H*(interval[2]-interval[1]))
+            if HT <: Operator
+                return exp(-1im * U.H * (interval[2] - interval[1]))
             else
-                return exp(U.H*(interval[2]-interval[1]))
+                return exp(U.H * (interval[2] - interval[1]))
             end
         end
 
         if HT <: Operator
-            return sesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]     
+            return sesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]
         else
-            return mesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]     
+            return mesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]
         end
     end
 end
@@ -239,13 +239,13 @@ Returns a `NamedTuple` with:
 - `usable`: Stored intervals fully contained within `[a, b]` (within `threshold` tolerance).
 - `to_compute`: Gap intervals not covered by any stored interval that still need to be computed.
 """
-function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_interval::Vector; threshold=1e-9) where T<:Vector
+function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_interval::Vector; threshold = 1.0e-9) where {T <: Vector}
     a, b = target_interval
-    
+
     # Find stored intervals that are fully contained within target range (with fuzzy boundaries)
     usable = [[s, e] for (s, e) in stored_intervals if s >= a - threshold && e <= b + threshold]
     sort!(usable, by = first)
-    
+
     # Merge usable intervals to find coverage
     merged = Vector[]
     for (s, e) in usable
@@ -255,7 +255,7 @@ function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_in
             merged[end] = (merged[end][1], max(merged[end][2], e))
         end
     end
-    
+
     # Find gaps that need to be computed
     to_compute = Vector[]
     current = a
@@ -281,16 +281,16 @@ function Base.show(io::IO, U::Propagator)
     if length(saved_times) == 0
         saved_times = ["None"]
     end
-    println(
-    io,
-    "\nPropagator: H Type=",
-    U.H.type,
-    "   dims=",
-    _get_dims_string(U.dimensions),
-    "   size=",
-    size(U),
-    "\nSaved Propagators: ",
-    saved_times...
+    return println(
+        io,
+        "\nPropagator: H Type=",
+        U.H.type,
+        "   dims=",
+        _get_dims_string(U.dimensions),
+        "   size=",
+        size(U),
+        "\nSaved Propagators: ",
+        saved_times...
     )
 end
 
@@ -301,4 +301,3 @@ end
 function Base.length(U::Propagator)
     return length(U.H)
 end
-
diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl
index 39c65496a..75aa2deaf 100644
--- a/test/core-test/propagator.jl
+++ b/test/core-test/propagator.jl
@@ -38,4 +38,4 @@ end
         @test isapprox(U_se^n * ψ0, ψt[n]; atol = 1.0e-5)
         @test isapprox(U_me^n * ρ0, ρt[n]; atol = 1.0e-5)
     end
-end
\ No newline at end of file
+end

From 2b3a8e2b6a4ee3786fa3ebb95adce17b07290684 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 06:40:28 +0000
Subject: [PATCH 05/19] fixed spelling issue

---
 src/time_evolution/propagator.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 0fde19d6d..4c0f58810 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -24,7 +24,7 @@ Previously computed propagators for sub-intervals are reused automatically to av
 - `dims`: The dimensions of the Hilbert space.
 - `solver_kwargs`: Keyword arguments forwarded to the underlying solver ([`sesolve`](@ref) or [`mesolve`](@ref)).
 - `max_saved`: Maximum number of propagators to cache.
-- `threshhold`: Numerical tolerance for matching stored time intervals.
+- `threshold`: Numerical tolerance for matching stored time intervals.
 - `remember_by_default`: Whether to cache newly computed propagators by default.
 - `period`: If the system is periodic in time, the propagator will calculate intervals modulo this period.
 
@@ -45,7 +45,7 @@ struct Propagator{
     dimensions::DT
     solver_kwargs::KWT
     max_saved::Union{Integer, Float64}
-    threshhold::Float64
+    threshold::Float64
     remember_by_default::Bool
     period::Real
 end
@@ -56,7 +56,7 @@ end
         H::AbstractQuantumObject{HOpType},
         t::Union{Nothing, Real} = nothing;
         t0 = 0.0,
-        threshhold::Float64 = 1e-9,
+        threshold::Float64 = 1e-9,
         max_saved::Union{Integer, Float64} = typemax(Int),
         remember_by_default::Bool = true,
         params = NullParameters(),
@@ -76,7 +76,7 @@ If `t` is provided, the propagator from `t0` to `t` is immediately computed and
   [`QuantumObjectEvolution`](@ref), with type [`Operator`](@ref) or [`SuperOperator`](@ref).
 - `t`: Optional final time. If given, the propagator for the interval `[t0, t]` is computed immediately.
 - `t0`: Initial time. Default is `0.0`.
-- `threshhold`: Numerical tolerance for matching cached time intervals. Default is `1e-9`.
+- `threshold`: Numerical tolerance for matching cached time intervals. Default is `1e-9`.
 - `period`: If the system is periodic in time, specify the period to automatically calculate intervals modulo this period. Default is `Inf` (no periodicity).
 - `max_saved`: Maximum number of propagators to store in the cache. Can be an `Integer` or `Inf`. Default is `typemax(Int)`.
 - `remember_by_default`: Whether to automatically cache computed propagators. Default is `true`.
@@ -101,7 +101,7 @@ function propagator(
         H::AbstractQuantumObject{HOpType},
         t::Union{Nothing, Real} = nothing;
         t0 = 0.0,
-        threshhold::Float64 = 1.0e-9,
+        threshold::Float64 = 1.0e-9,
         period::Real = Inf,
         max_saved::Union{Integer, Float64} = typemax(Int),
         remember_by_default::Bool = true,
@@ -118,7 +118,7 @@ function propagator(
         @warn "max_saved should be an Integer or Inf. Setting to $max_saved."
     end
 
-    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshhold, remember_by_default, period)
+    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, period)
 
     if t != nothing
         U(t; t0 = t0, remember = true)
@@ -155,7 +155,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
         t = mod(t, U.period)
         t0 = mod(t0, U.period)
     end
-    intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold = U.threshhold)
+    intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold = U.threshold)
 
     prop = qeye_like(U.H)
     if prop isa QobjEvo

From 45b017bfe49968293e3f642ca1e33f5bd9126090 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 06:53:58 +0000
Subject: [PATCH 06/19] removed print statement from debuging.

---
 src/time_evolution/propagator.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 4c0f58810..05059dcf7 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -213,7 +213,6 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <
         return U.props[interval]
     else
         if !(U.H isa QobjEvo)
-            println("Computing propagator for interval [$(interval[1]), $(interval[2])] via matrix exponentiation...")
             if HT <: Operator
                 return exp(-1im * U.H * (interval[2] - interval[1]))
             else

From d1970bb5e4ab6269d01e0e2662ef9b7c59a39cc2 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 07:15:55 +0000
Subject: [PATCH 07/19] updated `_propagator_compute_or_look_up` to fix the
 code quality issue.

---
 src/time_evolution/propagator.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 05059dcf7..ab8270efb 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -221,9 +221,9 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <
         end
 
         if HT <: Operator
-            return sesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]
+            return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, interval; U.solver_kwargs...).states[end]
         else
-            return mesolve(U.H, qeye_like(U.H)(0.0), interval; U.solver_kwargs...).states[end]
+            return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, interval; U.solver_kwargs...).states[end]
         end
     end
 end

From fa97360aabd628347f581b92693050eeb0ecdde0 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 17:01:31 +0000
Subject: [PATCH 08/19] updated saveat to only save the final state.

---
 src/time_evolution/propagator.jl | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index ab8270efb..cbf6c007e 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -48,6 +48,7 @@ struct Propagator{
     threshold::Float64
     remember_by_default::Bool
     period::Real
+    isconstant::Bool
 end
 
 
@@ -103,6 +104,7 @@ function propagator(
         t0 = 0.0,
         threshold::Float64 = 1.0e-9,
         period::Real = Inf,
+        isconstant::Bool = false,
         max_saved::Union{Integer, Float64} = typemax(Int),
         remember_by_default::Bool = true,
         params = NullParameters(),
@@ -118,7 +120,10 @@ function propagator(
         @warn "max_saved should be an Integer or Inf. Setting to $max_saved."
     end
 
-    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, period)
+    if !(H isa QobjEvo)
+        isconstant = true
+    end
+    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, period, isconstant)
 
     if t != nothing
         U(t; t0 = t0, remember = true)
@@ -205,14 +210,14 @@ the propagator is computed via [`sesolve`](@ref) (for [`Operator`](@ref)) or [`m
 [`SuperOperator`](@ref)) using an identity matrix as the initial state.
 """
 function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <: Union{Operator, SuperOperator}}
-    if !(U.H isa QobjEvo)
+    if U.isconstant
         interval = [0.0, interval[2] - interval[1]]
     end
 
     if interval in keys(U.props)
         return U.props[interval]
     else
-        if !(U.H isa QobjEvo)
+        if U.isconstant
             if HT <: Operator
                 return exp(-1im * U.H * (interval[2] - interval[1]))
             else
@@ -221,9 +226,9 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <
         end
 
         if HT <: Operator
-            return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, interval; U.solver_kwargs...).states[end]
+            return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, interval; saveat = [interval[2]], U.solver_kwargs...).states[end]
         else
-            return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, interval; U.solver_kwargs...).states[end]
+            return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, interval; saveat = [interval[2]], U.solver_kwargs...).states[end]
         end
     end
 end

From f31d8618194f9eda94ebe853e8420f1c83f10b01 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 17:41:20 +0000
Subject: [PATCH 09/19] updated show to include displaying the memory usage.

---
 src/time_evolution/propagator.jl | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index cbf6c007e..a313b4b1a 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -293,8 +293,12 @@ function Base.show(io::IO, U::Propagator)
         _get_dims_string(U.dimensions),
         "   size=",
         size(U),
+        "\nperiod: ",
+        U.period,
         "\nSaved Propagators: ",
-        saved_times...
+        saved_times...,
+        "\nMemory Usage: ",
+        Base.format_bytes(Base.summarysize(U.props))
     )
 end
 
@@ -305,3 +309,15 @@ end
 function Base.length(U::Propagator)
     return length(U.H)
 end
+
+function Base.pop!(U::Propagator, interval::Vector)
+    if U.period != Inf
+        interval = mod.(interval, U.period)
+    end
+    if interval in keys(U.props)
+        return pop!(U.props, interval)
+    else
+        @warn "Interval $interval not found in cache. No propagator removed."
+        return nothing
+    end
+end
\ No newline at end of file

From 6cbe4147333201cf456affa61308e5d41150d6a5 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 17:53:04 +0000
Subject: [PATCH 10/19] fixed a bug for when the period is finite.

---
 src/time_evolution/propagator.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index a313b4b1a..eeffad779 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -156,6 +156,7 @@ propagator is also cached separately when `remember` is enabled.
 - The propagator as an `AbstractQuantumObject` (if `return_result` is `true`).
 """
 function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing, return_result = true, save_steps = true)
+    ΔT = abs(t - t0)
     if U.period != Inf
         t = mod(t, U.period)
         t0 = mod(t0, U.period)
@@ -181,6 +182,10 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
         prop = temp_prop * prop
     end
 
+    if U.period != Inf
+        prop = prop^(ΔT / U.period)
+    end
+
     if (remember === nothing ? (U.remember_by_default && length(U.props) < U.max_saved) : remember) && !([t0, t] in keys(U.props))
         if length(U.props) >= U.max_saved
             @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."

From fb11c94292a3f39425bf905af0ea526a5d8ca7bd Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 17:53:06 +0000
Subject: [PATCH 11/19] minor changes.

---
 src/time_evolution/propagator.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index eeffad779..6b9ed3863 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -182,7 +182,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
         prop = temp_prop * prop
     end
 
-    if U.period != Inf
+    if U.period != Inf && (t - t0) > U.period
         prop = prop^(ΔT / U.period)
     end
 

From ca02c77b564576b7c8c7792a3b2ea2fa459259f2 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 18:02:39 +0000
Subject: [PATCH 12/19] removed the `period` functionality. Needed more work.

---
 src/time_evolution/propagator.jl | 20 +-------------------
 1 file changed, 1 insertion(+), 19 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 6b9ed3863..fc9efc746 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -26,7 +26,6 @@ Previously computed propagators for sub-intervals are reused automatically to av
 - `max_saved`: Maximum number of propagators to cache.
 - `threshold`: Numerical tolerance for matching stored time intervals.
 - `remember_by_default`: Whether to cache newly computed propagators by default.
-- `period`: If the system is periodic in time, the propagator will calculate intervals modulo this period.
 
 # Usage
 
@@ -47,7 +46,6 @@ struct Propagator{
     max_saved::Union{Integer, Float64}
     threshold::Float64
     remember_by_default::Bool
-    period::Real
     isconstant::Bool
 end
 
@@ -78,7 +76,6 @@ If `t` is provided, the propagator from `t0` to `t` is immediately computed and
 - `t`: Optional final time. If given, the propagator for the interval `[t0, t]` is computed immediately.
 - `t0`: Initial time. Default is `0.0`.
 - `threshold`: Numerical tolerance for matching cached time intervals. Default is `1e-9`.
-- `period`: If the system is periodic in time, specify the period to automatically calculate intervals modulo this period. Default is `Inf` (no periodicity).
 - `max_saved`: Maximum number of propagators to store in the cache. Can be an `Integer` or `Inf`. Default is `typemax(Int)`.
 - `remember_by_default`: Whether to automatically cache computed propagators. Default is `true`.
 - `params`: Parameters to pass to the underlying solver.
@@ -103,7 +100,6 @@ function propagator(
         t::Union{Nothing, Real} = nothing;
         t0 = 0.0,
         threshold::Float64 = 1.0e-9,
-        period::Real = Inf,
         isconstant::Bool = false,
         max_saved::Union{Integer, Float64} = typemax(Int),
         remember_by_default::Bool = true,
@@ -123,7 +119,7 @@ function propagator(
     if !(H isa QobjEvo)
         isconstant = true
     end
-    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, period, isconstant)
+    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, isconstant)
 
     if t != nothing
         U(t; t0 = t0, remember = true)
@@ -156,11 +152,6 @@ propagator is also cached separately when `remember` is enabled.
 - The propagator as an `AbstractQuantumObject` (if `return_result` is `true`).
 """
 function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing, return_result = true, save_steps = true)
-    ΔT = abs(t - t0)
-    if U.period != Inf
-        t = mod(t, U.period)
-        t0 = mod(t0, U.period)
-    end
     intervals = _get_intervals_for_range(collect(keys(U.props)), [t0, t]; threshold = U.threshold)
 
     prop = qeye_like(U.H)
@@ -182,10 +173,6 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
         prop = temp_prop * prop
     end
 
-    if U.period != Inf && (t - t0) > U.period
-        prop = prop^(ΔT / U.period)
-    end
-
     if (remember === nothing ? (U.remember_by_default && length(U.props) < U.max_saved) : remember) && !([t0, t] in keys(U.props))
         if length(U.props) >= U.max_saved
             @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
@@ -298,8 +285,6 @@ function Base.show(io::IO, U::Propagator)
         _get_dims_string(U.dimensions),
         "   size=",
         size(U),
-        "\nperiod: ",
-        U.period,
         "\nSaved Propagators: ",
         saved_times...,
         "\nMemory Usage: ",
@@ -316,9 +301,6 @@ function Base.length(U::Propagator)
 end
 
 function Base.pop!(U::Propagator, interval::Vector)
-    if U.period != Inf
-        interval = mod.(interval, U.period)
-    end
     if interval in keys(U.props)
         return pop!(U.props, interval)
     else

From 867e0a69e83f442259e2f4545401bbad6bc7e50c Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 18:05:00 +0000
Subject: [PATCH 13/19] fixed a bug in _get_intervals_for_range. Appended tuple
 instead of vector.

---
 src/time_evolution/propagator.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index fc9efc746..2eb7fdef8 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -248,7 +248,7 @@ function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_in
         if isempty(merged) || s > merged[end][2] + threshold
             push!(merged, [s, e])
         else
-            merged[end] = (merged[end][1], max(merged[end][2], e))
+            merged[end] = [merged[end][1], max(merged[end][2], e)]
         end
     end
 

From c89e3dd8c932271c937f9a550a5b6916d8e9c5bb Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Thu, 12 Feb 2026 20:35:46 +0000
Subject: [PATCH 14/19] changed a docstring.

---
 src/time_evolution/propagator.jl | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 2eb7fdef8..0d0afa0a6 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -86,14 +86,6 @@ If `t` is provided, the propagator from `t0` to `t` is immediately computed and
 # Returns
 
 - `U::Propagator`: A callable propagator object. Use `U(t; t0=0.0)` or `U([t0, t])` to evaluate.
-
-# Example
-
-```julia
-H = (ϵ / 2) * sigmaz() + (Ω / 2) * sigmax()
-U = propagator(H)      # create lazy propagator
-ψt = U(π) * basis(2, 0) # propagate |0⟩ to time π
-```
 """
 function propagator(
         H::AbstractQuantumObject{HOpType},

From 5cce96504ba48fb15d74f45774b32b04758248ff Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Mon, 27 Jul 2026 16:04:20 +0000
Subject: [PATCH 15/19] some small updates and fixes to keep it working with
 the newer versions.

---
 src/time_evolution/propagator.jl | 56 ++++++++++++++++++--------------
 1 file changed, 32 insertions(+), 24 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 0d0afa0a6..609f6b1c2 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -34,13 +34,14 @@ an interval. See [`propagator`](@ref) for construction.
 """
 struct Propagator{
         HT <: Union{Operator, SuperOperator},
-        PT <: AbstractQuantumObject,
-        DT <: AbstractDimensions,
+        HType <: AbstractQuantumObject{HT},
+        PT <: QuantumObject{HT},
+        DT <: Dimensions,
         KWT,
     }
-    H::AbstractQuantumObject{HT}
-    props::Dict{Vector, PT}
-    dims::AbstractArray
+    H::HType
+    props::Dict{NTuple{2, Float64}, PT}
+    dims::Tuple
     dimensions::DT
     solver_kwargs::KWT
     max_saved::Union{Integer, Float64}
@@ -88,7 +89,7 @@ If `t` is provided, the propagator from `t0` to `t` is immediately computed and
 - `U::Propagator`: A callable propagator object. Use `U(t; t0=0.0)` or `U([t0, t])` to evaluate.
 """
 function propagator(
-        H::AbstractQuantumObject{HOpType},
+        H::HType,
         t::Union{Nothing, Real} = nothing;
         t0 = 0.0,
         threshold::Float64 = 1.0e-9,
@@ -99,7 +100,7 @@ function propagator(
         progress_bar::Union{Val, Bool} = Val(true),
         inplace::Union{Val, Bool} = Val(true),
         kwargs...,
-    ) where {HOpType <: Union{Operator, SuperOperator}}
+    ) where {HOpType <: Union{Operator, SuperOperator}, HType <: AbstractQuantumObject{HOpType}}
 
     full_kwargs = (; params, progress_bar, inplace, kwargs...)
 
@@ -111,7 +112,8 @@ function propagator(
     if !(H isa QobjEvo)
         isconstant = true
     end
-    U = Propagator(H, Dict{Vector, AbstractQuantumObject}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, isconstant)
+    
+    U = Propagator(H, Dict{NTuple{2, Float64}, QuantumObject{HOpType}}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, isconstant)
 
     if t != nothing
         U(t; t0 = t0, remember = true)
@@ -169,11 +171,11 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
         if length(U.props) >= U.max_saved
             @warn "Maximum number of stored propagators reached, save is being forced because 'remember' is set to true."
         end
-        U.props[[t0, t]] = prop
+        U.props[(t0, t)] = prop
     end
 
     if return_result
-        return prop
+       return prop
     end
 end
 
@@ -209,14 +211,18 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <
             end
         end
 
-        if HT <: Operator
-            return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, interval; saveat = [interval[2]], U.solver_kwargs...).states[end]
-        else
-            return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, interval; saveat = [interval[2]], U.solver_kwargs...).states[end]
-        end
+        _get_new_propagator(U, interval)
+        
     end
 end
 
+function _get_new_propagator(U::Propagator{Operator}, interval)
+        return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, collect(interval); saveat = [interval[2]], U.solver_kwargs...).states[end]
+end
+function _get_new_propagator(U::Propagator{SuperOperator}, interval)
+    return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, collect(interval); saveat = [interval[2]], U.solver_kwargs...).states[end]
+end
+        
 
 """
     _get_intervals_for_range(stored_intervals, target_interval; threshold=1e-9)
@@ -227,34 +233,34 @@ Returns a `NamedTuple` with:
 - `usable`: Stored intervals fully contained within `[a, b]` (within `threshold` tolerance).
 - `to_compute`: Gap intervals not covered by any stored interval that still need to be computed.
 """
-function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_interval::Vector; threshold = 1.0e-9) where {T <: Vector}
+function _get_intervals_for_range(stored_intervals::AbstractVector{T}, target_interval::Vector; threshold = 1.0e-9) where {T <: Tuple}
     a, b = target_interval
 
     # Find stored intervals that are fully contained within target range (with fuzzy boundaries)
-    usable = [[s, e] for (s, e) in stored_intervals if s >= a - threshold && e <= b + threshold]
+    usable = [(s, e) for (s, e) in stored_intervals if s >= a - threshold && e <= b + threshold]
     sort!(usable, by = first)
 
     # Merge usable intervals to find coverage
-    merged = Vector[]
+    merged = Tuple[]
     for (s, e) in usable
         if isempty(merged) || s > merged[end][2] + threshold
-            push!(merged, [s, e])
+            push!(merged, (s, e))
         else
-            merged[end] = [merged[end][1], max(merged[end][2], e)]
+            merged[end] = (merged[end][1], max(merged[end][2], e))
         end
     end
 
     # Find gaps that need to be computed
-    to_compute = Vector[]
+    to_compute = Tuple[]
     current = a
     for (s, e) in merged
         if current < s - threshold
-            push!(to_compute, [current, s])
+            push!(to_compute, (current, s))
         end
         current = max(current, e)
     end
     if current < b - threshold
-        push!(to_compute, [current, b])
+        push!(to_compute, (current, b))
     end
     return (usable = usable, to_compute = to_compute)
 end
@@ -271,12 +277,14 @@ function Base.show(io::IO, U::Propagator)
     end
     return println(
         io,
-        "\nPropagator: H Type=",
+        "\nPropagator: type=",
         U.H.type,
         "   dims=",
         _get_dims_string(U.dimensions),
         "   size=",
         size(U),
+        "\nU Is ObjEvo: ",
+        (U.H isa QobjEvo),
         "\nSaved Propagators: ",
         saved_times...,
         "\nMemory Usage: ",

From 2c83a536969b51b53d18cb857e8f00d0edd77246 Mon Sep 17 00:00:00 2001
From: Gavin <gavin.rockwood@gmail.com>
Date: Tue, 28 Jul 2026 13:27:48 +0000
Subject: [PATCH 16/19] made changes, ran the "make test" and "make format".
 "make docs" is currently crashing.

---
 src/time_evolution/propagator.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
index 609f6b1c2..11b9916cc 100644
--- a/src/time_evolution/propagator.jl
+++ b/src/time_evolution/propagator.jl
@@ -112,7 +112,7 @@ function propagator(
     if !(H isa QobjEvo)
         isconstant = true
     end
-    
+
     U = Propagator(H, Dict{NTuple{2, Float64}, QuantumObject{HOpType}}(), H.dims, H.dimensions, full_kwargs, max_saved, threshold, remember_by_default, isconstant)
 
     if t != nothing
@@ -175,7 +175,7 @@ function (U::Propagator)(t; t0 = 0.0, remember::Union{Nothing, Bool} = nothing,
     end
 
     if return_result
-       return prop
+        return prop
     end
 end
 
@@ -212,17 +212,17 @@ function _propagator_compute_or_look_up(U::Propagator{HT}, interval) where {HT <
         end
 
         _get_new_propagator(U, interval)
-        
+
     end
 end
 
 function _get_new_propagator(U::Propagator{Operator}, interval)
-        return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, collect(interval); saveat = [interval[2]], U.solver_kwargs...).states[end]
+    return sesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{Operator}, collect(interval); saveat = [interval[2]], U.solver_kwargs...).states[end]
 end
 function _get_new_propagator(U::Propagator{SuperOperator}, interval)
     return mesolve(U.H, qeye_like(U.H)(0.0)::QuantumObject{SuperOperator}, collect(interval); saveat = [interval[2]], U.solver_kwargs...).states[end]
 end
-        
+
 
 """
     _get_intervals_for_range(stored_intervals, target_interval; threshold=1e-9)
@@ -283,7 +283,7 @@ function Base.show(io::IO, U::Propagator)
         _get_dims_string(U.dimensions),
         "   size=",
         size(U),
-        "\nU Is ObjEvo: ",
+        "\nH Is ObjEvo: ",
         (U.H isa QobjEvo),
         "\nSaved Propagators: ",
         saved_times...,
@@ -307,4 +307,4 @@ function Base.pop!(U::Propagator, interval::Vector)
         @warn "Interval $interval not found in cache. No propagator removed."
         return nothing
     end
-end
\ No newline at end of file
+end

From 2d69818089e4fad195b26f80a69f4981fb089963 Mon Sep 17 00:00:00 2001
From: Gavin Rockwood <gavin.rockwood@gmail.com>
Date: Tue, 28 Jul 2026 18:24:33 +0000
Subject: [PATCH 17/19] updated propagator.md and reran format and test.

---
 CHANGELOG.md                                      | 2 +-
 docs/src/users_guide/time_evolution/propagator.md | 7 +++++++
 2 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 259943418..7404d0114 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Add documentation about arbitrary precision computations. ([#745])
 - Fix `e_ops` expectation values being always stored as `ComplexF64` in `sesolve`, `mesolve`, `mcsolve`, `ssesolve`, and `smesolve`. They now follow the element type of the problem, so that arbitrary precision solutions return `sol.expect` with the requested precision instead of silently narrowing it to double precision. ([#745])
 - Fix `sesolve_map` and `mesolve_map` sharing mutable e_ops-saving callback state across trajectories when a custom `prob_func` is supplied, which could throw a `BoundsError` or produce incorrect results. Both now accept a `safetycopy` keyword (smart-defaulting to `true` for a custom `prob_func`, `false` otherwise) mirroring `SciMLBase.EnsembleProblem`. ([#645], [#747])
+- Added a propagator structure. 
 
 ## [v0.47.2]
 Release date: 2026-06-21
@@ -553,4 +554,3 @@ Release date: 2024-11-13
 [#736]: https://github.com/qutip/QuantumToolbox.jl/issues/736
 [#745]: https://github.com/qutip/QuantumToolbox.jl/issues/745
 [#747]: https://github.com/qutip/QuantumToolbox.jl/issues/747
-
diff --git a/docs/src/users_guide/time_evolution/propagator.md b/docs/src/users_guide/time_evolution/propagator.md
index c062c2711..9556c986a 100644
--- a/docs/src/users_guide/time_evolution/propagator.md
+++ b/docs/src/users_guide/time_evolution/propagator.md
@@ -59,3 +59,10 @@ n = 10
 
 ρt_vec = U_me^n * ρ0_vec
 ```
+## The Propagator Structure
+For problems that require many propagator evaluations, there is the [`Propagator`][@ref] object. Taking either a Hamiltonian operator or a Liouvillian superoperator, the propagator structure built as
+```@example propagator
+U = propagator(H)
+```
+and can be evaluated via ``U(tf, t0 = t0)``. By default, this call returns the calculated [`Qobject`](@ref) as well as stores it. Further evaluations will check whether or not the current evalutaion window overlaps with already calculated propagators and the time interval will be split up to make use of those that have already been computed.
+

From 814fac6ef2144907b6ec673e011fad44cfca4efd Mon Sep 17 00:00:00 2001
From: Gavin Rockwood <gavin.rockwood@gmail.com>
Date: Wed, 29 Jul 2026 14:39:05 +0000
Subject: [PATCH 18/19] fixed typo in docs.

---
 docs/src/users_guide/time_evolution/propagator.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/users_guide/time_evolution/propagator.md b/docs/src/users_guide/time_evolution/propagator.md
index 9556c986a..f76ea5cc6 100644
--- a/docs/src/users_guide/time_evolution/propagator.md
+++ b/docs/src/users_guide/time_evolution/propagator.md
@@ -64,5 +64,5 @@ For problems that require many propagator evaluations, there is the [`Propagator
 ```@example propagator
 U = propagator(H)
 ```
-and can be evaluated via ``U(tf, t0 = t0)``. By default, this call returns the calculated [`Qobject`](@ref) as well as stores it. Further evaluations will check whether or not the current evalutaion window overlaps with already calculated propagators and the time interval will be split up to make use of those that have already been computed.
+and can be evaluated via ``U(tf, t0 = t0)``. By default, this call returns the calculated [`Qobject`](@ref) as well as stores it. Further evaluations will check whether or not the current evaluation window overlaps with already calculated propagators and the time interval will be split up to make use of those that have already been computed.
 

From 54dfbafc69c77750e97eaf3d03b82ec68a3c5968 Mon Sep 17 00:00:00 2001
From: Gavin Rockwood <gavin.rockwood@gmail.com>
Date: Wed, 29 Jul 2026 15:06:12 +0000
Subject: [PATCH 19/19] fixed a doc error.

---
 docs/src/users_guide/time_evolution/propagator.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/users_guide/time_evolution/propagator.md b/docs/src/users_guide/time_evolution/propagator.md
index f76ea5cc6..3a5c34af7 100644
--- a/docs/src/users_guide/time_evolution/propagator.md
+++ b/docs/src/users_guide/time_evolution/propagator.md
@@ -64,5 +64,5 @@ For problems that require many propagator evaluations, there is the [`Propagator
 ```@example propagator
 U = propagator(H)
 ```
-and can be evaluated via ``U(tf, t0 = t0)``. By default, this call returns the calculated [`Qobject`](@ref) as well as stores it. Further evaluations will check whether or not the current evaluation window overlaps with already calculated propagators and the time interval will be split up to make use of those that have already been computed.
+and can be evaluated via ``U(tf, t0 = t0)``. By default, this call returns the calculated [`QuantumObject`](@ref) as well as stores it. Further evaluations will check whether or not the current evaluation window overlaps with already calculated propagators and the time interval will be split up to make use of those that have already been computed.