diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 64bb570b0..1e8ea3fa9 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -24,6 +24,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs
   for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion
   that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252))
+- In-place `timestep!` for an `InfiniteMPS` with `TDVP`, which also enables `time_evolve!` on
+  infinite states. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the
   almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and
   `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
@@ -33,6 +35,17 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain
   instead of a single site, and return the gauge transformations alongside the operator, matching
   the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+- Internal: the finite `DMRG`/`DMRG2` ground-state search and approximation, `DynamicalDMRG`,
+  multiline `IDMRG`/`IDMRG2` approximation and `leading_boundary`, and `time_evolve` now run as
+  iterators whose iterations are a forward and a backward half-sweep (or a single time step),
+  sharing one structure across solvers. Finite `TDVP`, `TDVP2` and `BUG` steps share a single sweep
+  over per-site local updates, and infinite `TDVP` steps are split into local update, gauge and
+  environment steps on the same step state. Results and logging are unchanged. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- `time_evolve` with finite `TDVP`, `TDVP2` or `BUG`, or infinite `TDVP`, copies (or promotes to
+  complex) the initial state once and then takes every step in place, reusing one step state and
+  its allocator across steps, rather than copying the state at every step. A `finalize` callback
+  that keeps the state of an earlier step therefore sees it updated by later steps, as with
+  `time_evolve!`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 
 ### Deprecated
 
@@ -44,6 +57,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538))
 - The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when
   interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+- `InfiniteMPS(C, AR)` no longer throws a `MethodError`, which also fixes infinite `TDVP` with
+  `leftorthflag = false`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- Real-time infinite `TDVP` from a real `InfiniteMPS` with a real Hamiltonian no longer fails with
+  an `InexactError`: the environments are now rebuilt for the complex state, as for finite MPS.
+  ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- `TDVP`, `TDVP2` and `BUG` promote a real state to complex for imaginary-time evolution with a
+  complex time step, rather than failing. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
 
 ### Performance
 
diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index 465b40f1b..ed17ce1a0 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -1,74 +1,82 @@
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...))
+# Internal state of the finite DMRG/DMRG2 approximation, where `ϵ` is the largest relative
+# local change of the last sweep
+struct ApproximateState{S, O, E, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+function approximate!(
+        ψ::AbstractFiniteMPS, Oϕ, alg::Union{DMRG, DMRG2},
+        envs = environments(ψ, _environment_args(Oϕ)...)
+    )
     allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG2")
+    log = IterLog(alg)
+    it = IterativeSolver(alg, ApproximateState(ψ, Oϕ, envs, 0, 2 * alg.tol, allocator))
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = AC2_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(alg))
-
-                AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
-                ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
-
-                ψ.AC[pos] = (al, complex(c))
-                ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-            end
-
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
-
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, _, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+    state = it.state
+    info = AlgorithmInfo(; converged = state.ϵ <= alg.tol, localchange = state.ϵ, numiter = state.iter)
+    return state.mps, state.envs, info
 end
 
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...))
-    allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                AC = ψ.AC[pos]
-                ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::ApproximateState)
+    iter = state.iter + 1
+    state = ApproximateState(
+        state.mps, state.operator, state.envs, state.iter, 0.0, state.allocator
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ψ, envs = it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = ApproximateState(ψ, state.operator, envs, iter, state.ϵ, state.allocator)
+    return (ψ, envs, state.ϵ), it.state
+end
 
-                ψ.AC[pos] = AC′
-            end
+function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC2′ = AC2_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(it.alg))
 
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
+        AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+        ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AC[pos] = (al, complex(c))
+        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
     end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC′ = AC_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        AC = ψ.AC[pos]
+        ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+
+        ψ.AC[pos] = AC′
+    end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
 end
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 8b8079379..9fe55fbf7 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -1,199 +1,182 @@
 function approximate!(
-        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS}, alg::IDMRG,
-        envs = environments(ψ, toapprox...)
+        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS},
+        alg::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox...)
     )
     allocator = default_allocator(ψ, SerialScheduler())
-    log = IterLog("IDMRG")
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
+    alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
+    log = IterLog(alg)
+    ϵ_truncs = alg isa IDMRG2 ?
+        PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing
+    state = IDMRGState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator)
+    it = IterativeSolver(alg, state)
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            C_current = ψ.C[:, 0]
-
-            # left to right sweep
-            for col in 1:width(ψ)
-                for row in 1:size(ψ, 1)
-                    ψ.AC[row + 1, col] = AC_projection(
-                        CartesianIndex(row, col), ψ, toapprox, envs;
-                        alg.backend, allocator
-                    )
-                    normalize!(ψ.AC[row + 1, col])
-                    ψ.AL[row + 1, col], ψ.C[row + 1, col] = left_orth!(ψ.AC[row + 1, col])
-                end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-            end
-
-            # right to left sweep
-            for col in reverse(1:width(ψ))
-                for row in 1:size(ψ, 1)
-                    ψ.AC[row + 1, col] = AC_projection(
-                        CartesianIndex(row, col), ψ, toapprox, envs;
-                        alg.backend, allocator
-                    )
-                    normalize!(ψ.AC[row + 1, col])
-                    ψ.C[row + 1, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row + 1, col]))
-                    ψ.AR[row + 1, col] = _transpose_front(temp)
-                end
-                transfer_rightenv!(envs, ψ, toapprox, col - 1)
-            end
-            normalize!(envs, ψ, toapprox)
-
-            ϵ = norm(C_current - ψ.C[:, 0])
-
+        @log_initialization loginit!(log, it.ϵ)
+        for _ in Iterators.take(it, alg.maxiter)
+            ϵ = it.ϵ
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
+    (; iter, ϵ) = it.state
+
     # TODO: immediately compute in-place
     alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
     recalculate!(envs, ψ, toapprox)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
+    info = AlgorithmInfo(;
+        converged = ϵ <= alg.tol, bondresidual = ϵ,
+        truncation_errors = isnothing(ϵ_truncs) ? nothing : parent(ϵ_truncs), numiter = iter,
+    )
+    return ψ, envs, info
 end
 
-function approximate!(
-        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS},
-        alg::IDMRG2, envs = environments(ψ, toapprox...)
+function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGState{<:MultilineMPS, <:Tuple})
+    iter = state.iter + 1
+    C_current = state.mps.C[:, 0]
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ϵ = bond_change(C_current, state.mps.C[:, 0])
+    it.state = IDMRGState(
+        state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.energy,
+        state.timeroutput, state.allocator,
     )
-    allocator = default_allocator(ψ, SerialScheduler())
-    width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    ϵ::Float64 = 2 * alg.tol
-    log = IterLog("IDMRG2")
-    O, ϕ = toapprox
-    iter = 0
-    ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ)))
+    return (it.state.mps, it.state.envs, it.state.ϵ), it.state
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            C_current = ψ.C[:, 0]
-
-            # sweep from left to right
-            for site in 1:(width(ψ) - 1)
-                for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(
-                        CartesianIndex(row, site), ψ, toapprox, envs;
-                        kind = :ACAR, alg.backend, allocator
-                    )
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
-                end
-
-                transfer_leftenv!(envs, ψ, toapprox, site + 1)
-                transfer_rightenv!(envs, ψ, toapprox, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:right}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    for col in 1:width(ψ)
+        for row in 1:size(ψ, 1)
+            ψ.AC[row + 1, col] = AC_projection(
+                CartesianIndex(row, col), ψ, toapprox, envs;
+                alg.backend, allocator
+            )
+            normalize!(ψ.AC[row + 1, col])
+            ψ.AL[row + 1, col], ψ.C[row + 1, col] = left_orth!(ψ.AC[row + 1, col])
+        end
+        transfer_leftenv!(envs, ψ, toapprox, col + 1)
+    end
+    return state
+end
 
-            # update the edge
-            ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            for row in 1:size(ψ, 1)
-                AC2′ = AC2_projection(
-                    CartesianIndex(row, width(ψ)), ψ, toapprox, envs;
-                    kind = :ALAC, alg.backend, allocator
-                )
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
-
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
-
-                ψ.AC[row + 1, end] = _mul_tail(al, c)
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-                ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:left}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    for col in reverse(1:width(ψ))
+        for row in 1:size(ψ, 1)
+            ψ.AC[row + 1, col] = AC_projection(
+                CartesianIndex(row, col), ψ, toapprox, envs;
+                alg.backend, allocator
+            )
+            normalize!(ψ.AC[row + 1, col])
+            ψ.C[row + 1, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row + 1, col]))
+            ψ.AR[row + 1, col] = _transpose_front(temp)
+        end
+        transfer_rightenv!(envs, ψ, toapprox, col - 1)
+    end
+    normalize!(envs, ψ, toapprox)
+    return state
+end
 
-            end
-            # update environments
-            transfer_leftenv!(envs, ψ, toapprox, 1)
-            transfer_rightenv!(envs, ψ, toapprox, 0)
-
-            normalize!(envs, ψ, toapprox)
-
-            # sweep from right to left
-            for site in reverse(1:(width(ψ) - 1))
-                for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(
-                        CartesianIndex(row, site), ψ, toapprox, envs;
-                        kind = :ALAC, alg.backend, allocator
-                    )
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                end
-
-                transfer_leftenv!(envs, ψ, toapprox, site + 1)
-                transfer_rightenv!(envs, ψ, toapprox, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:right}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    ϵ_truncs = state.truncation_errors
+    for site in 1:(width(ψ) - 1)
+        for row in 1:size(ψ, 1)
+            AC2′ = AC2_projection(
+                CartesianIndex(row, site), ψ, toapprox, envs;
+                kind = :ACAR, alg.backend, allocator
+            )
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+            ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
+        end
 
-            # update the edge
-            ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
-            ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
-            for row in 1:size(ψ, 1)
-                AC2′ = AC2_projection(
-                    CartesianIndex(row, 0), ψ, toapprox, envs;
-                    kind = :ACAR, alg.backend, allocator
-                )
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
-
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
-
-                ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \ _transpose_tail(al * c))
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-            end
-            transfer_leftenv!(envs, ψ, toapprox, 1)
-            transfer_rightenv!(envs, ψ, toapprox, 0)
+        transfer_leftenv!(envs, ψ, toapprox, site + 1)
+        transfer_rightenv!(envs, ψ, toapprox, site)
+    end
 
-            normalize!(envs, ψ, toapprox)
+    # update the edge
+    ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    for row in 1:size(ψ, 1)
+        AC2′ = AC2_projection(
+            CartesianIndex(row, width(ψ)), ψ, toapprox, envs;
+            kind = :ALAC, alg.backend, allocator
+        )
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
+
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
+
+        ψ.AC[row + 1, end] = _mul_tail(al, c)
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+        ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+    end
+    transfer_leftenv!(envs, ψ, toapprox, 1)
+    transfer_rightenv!(envs, ψ, toapprox, 0)
 
-            # update error
-            ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
-                smallest = infimum(_firstspace(c1), _firstspace(c2))
-                e1 = isometry(_firstspace(c1), smallest)
-                e2 = isometry(_firstspace(c2), smallest)
-                return norm(e2' * c2 * e2 - e1' * c1 * e1)
-            end
+    normalize!(envs, ψ, toapprox)
+    return state
+end
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:left}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    ϵ_truncs = state.truncation_errors
+    for site in reverse(1:(width(ψ) - 1))
+        for row in 1:size(ψ, 1)
+            AC2′ = AC2_projection(
+                CartesianIndex(row, site), ψ, toapprox, envs;
+                kind = :ALAC, alg.backend, allocator
+            )
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
         end
+
+        transfer_leftenv!(envs, ψ, toapprox, site + 1)
+        transfer_rightenv!(envs, ψ, toapprox, site)
     end
 
-    # TODO: immediately compute in-place
-    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
-    copy!(ψ, ψ′) # ensure output destination is unchanged
+    # update the edge
+    ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
+    ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
+    for row in 1:size(ψ, 1)
+        AC2′ = AC2_projection(
+            CartesianIndex(row, 0), ψ, toapprox, envs;
+            kind = :ACAR, alg.backend, allocator
+        )
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
+
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
+
+        ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \ _transpose_tail(al * c))
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+    end
+    transfer_leftenv!(envs, ψ, toapprox, 1)
+    transfer_rightenv!(envs, ψ, toapprox, 0)
 
-    recalculate!(envs, ψ, toapprox)
-    info = AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), numiter = iter)
-    return ψ, envs, info
+    normalize!(envs, ψ, toapprox)
+    return state
 end
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index 475913e08..361c846cc 100644
--- a/src/algorithms/approximate/vomps.jl
+++ b/src/algorithms/approximate/vomps.jl
@@ -21,7 +21,7 @@ function approximate(
 end
 
 function _approximate_vomps(mps, toapprox, alg::VOMPS, envs)
-    log = IterLog("VOMPS")
+    log = IterLog(alg)
     iter = 0
     ϵ = calc_galerkin(mps, toapprox..., envs; alg.backend)
     alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index b2400fdfa..1b9575e3a 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -295,78 +295,99 @@ function find_groundstate!(
     return find_groundstate_sweep!(ψ, H, alg, envs, allocator, timeroutput)
 end
 
-function find_groundstate_sweep!(
-        ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput
-    )
-    name = string(nameof(typeof(alg)))
-    log = IterLog(name)
+struct DMRGState{S, O, E, R, V, D, T, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::R
+    local_errors::V
+    truncation_errors::V
+    decay_rates::D
+    timeroutput::T
+    allocator::A
+end
 
+function DMRGState(ψ, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput)
     Tr = real(scalartype(ψ))
     n = _num_updates(alg, ψ)
-    ϵ_locals = ones(Tr, n)      # local Galerkin errors (drive both the eigensolver tol and the stop test)
-    ϵ_global = maximum(ϵ_locals) # sweep-wide Galerkin error, the convergence measure
-    ϵ_truncs = zeros(Tr, n)     # local truncation error
-    decay_rates = zeros(n)      # local observed decay rate of eigensolver
-    fwd, bwd = _sweep_ranges(alg, ψ)
-    iter = 0
+    local_errors = ones(Tr, n)
+    return DMRGState(
+        ψ, H, envs, 0, maximum(local_errors), local_errors, zeros(Tr, n),
+        zeros(n), timeroutput, allocator,
+    )
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ_global, expectation_value(ψ, H, envs))
-        for outer iter in 1:(alg.maxiter)
-            @timeit timeroutput "sweep" begin
-                # left-to-right
-                for pos in fwd
-                    ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] =
-                        local_update!(
-                        pos, Val(:right),
-                        ψ, H, alg, envs,
-                        ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        iter, timeroutput, allocator
-                    )
-                    ϵ_global = maximum(ϵ_locals)
-                end
-
-                # right-to-left
-                for pos in bwd
-                    ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] =
-                        local_update!(
-                        pos, Val(:left),
-                        ψ, H, alg, envs,
-                        ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        iter, timeroutput, allocator
-                    )
-                    ϵ_global = maximum(ϵ_locals)
-                end
-            end
-            ϵ_global = maximum(ϵ_locals)
+function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, ϵ = state.mps, state.ϵ
+    for pos in sites
+        ψ, state.local_errors[pos], state.truncation_errors[pos], state.decay_rates[pos] =
+            local_update!(
+            pos, direction, ψ, state.operator, it.alg, state.envs,
+            ϵ, state.truncation_errors[pos], state.decay_rates[pos],
+            iter, state.timeroutput, state.allocator,
+        )
+        # The next local solver uses the errors accumulated so far in this sweep.
+        ϵ = maximum(state.local_errors)
+    end
+    return DMRGState(
+        ψ, state.operator, state.envs, state.iter, ϵ, state.local_errors,
+        state.truncation_errors, state.decay_rates, state.timeroutput, state.allocator,
+    )
+end
+
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState)
+    iter = state.iter + 1
+    timeroutput = state.timeroutput
+    state = @timeit timeroutput "sweep" begin
+        state = sweep!(it, state, Val(:right), iter)
+        sweep!(it, state, Val(:left), iter)
+    end
+    ψ, envs = @timeit timeroutput "finalize" it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = DMRGState(
+        ψ, state.operator, envs, iter, state.ϵ, state.local_errors,
+        state.truncation_errors, state.decay_rates, timeroutput, state.allocator,
+    )
+    return (ψ, envs, state.ϵ), it.state
+end
 
-            ψ, envs = @timeit timeroutput "finalize" alg.finalize(
-                iter, ψ, H, envs
-            )::Tuple{typeof(ψ), typeof(envs)}
+# Truncation sets the attainable floor of the Galerkin error.
+sweep_converged(alg, state::DMRGState) =
+    state.ϵ <= max(alg.tol, maximum(state.truncation_errors))
 
-            # Truncation-aware convergence: the Galerkin gradient cannot drop below the level set by
-            # the discarded weight, so a truncating scheme converges once `ϵ_global` reaches the
-            # truncation error rather than the (unreachable) bare `tol`. With no truncation
-            # (`ϵ_truncs .= 0`, e.g. single-site/QR gauge) this reduces to the plain `ϵ_global ≤ tol`.
-            if ϵ_global <= max(alg.tol, maximum(ϵ_truncs))
+function find_groundstate_sweep!(
+        ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput
+    )
+    log = IterLog(alg)
+    it = IterativeSolver(alg, DMRGState(ψ, H, alg, envs, allocator, timeroutput))
+
+    with_verbosity(; alg.verbosity) do
+        @log_initialization loginit!(log, it.ϵ, expectation_value(ψ, H, envs))
+        for (ψ, envs, ϵ) in Iterators.take(it, alg.maxiter)
+            if sweep_converged(alg, it.state)
                 @info TimerReport(timeroutput) _group = :mpskit_timing
-                @log_convergence logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_convergence logfinish!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
                 break
-            end
-            if iter == alg.maxiter
+            elseif it.iter == alg.maxiter
                 @info TimerReport(timeroutput) _group = :mpskit_timing
-                @log_nonconvergence logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_nonconvergence logcancel!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
             else
-                @log_iteration logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_iteration logiter!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
             end
         end
     end
 
+    state = it.state
     info = AlgorithmInfo(;
-        converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), galerkin = ϵ_global,
-        truncation_errors = _bond_truncation_errors(alg, ϵ_truncs), numiter = iter
+        converged = sweep_converged(alg, state), galerkin = state.ϵ,
+        truncation_errors = _bond_truncation_errors(alg, state.truncation_errors),
+        numiter = state.iter,
     )
-    return ψ, envs, info
+    return state.mps, state.envs, info
 end
 
 function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...)
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index a0c314148..afe939fc4 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -70,14 +70,14 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_bounda
     backend::B = Defaults.backend()
 end
 
-# Internal state of the IDMRG algorithm
+# Internal state of the IDMRG/IDMRG2 algorithms
 struct IDMRGState{S, O, E, V, T, TO, A}
     mps::S
     operator::O
     envs::E
     iter::Int
     ϵ::Float64 # TODO: Could be any <:Real
-    truncation_errors::V # per bond, of the most recent sweep only
+    truncation_errors::V
     energy::T
     timeroutput::TO
     allocator::A
@@ -99,7 +99,7 @@ end
 function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_type <: Union{<:IDMRG, <:IDMRG2}}
     (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= 2"))
     name = alg isa IDMRG ? "IDMRG" : "IDMRG2"
-    log = IterLog(name)
+    log = IterLog(alg)
     timeroutput = alg.verbosity > 3 ? TimerOutput(name) : NoTimerOutput()
     mps = copy(mps)
     iter = 0
@@ -147,24 +147,14 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_
 end
 
 function Base.iterate(
-        it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, <:Any, T} = it.state
+        it::IterativeSolver{alg_type}, state::IDMRGState{<:InfiniteMPS, <:Any, <:Any, <:Any, T}
     ) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, T}
     timeroutput = state.timeroutput
     ϵ_truncs = zero(state.truncation_errors) # fresh each sweep, filled by the sweep itself
     mps, envs, C_old, E_new = @timeit timeroutput "localupdate" localupdate_step!(it, state, ϵ_truncs)
 
     # error criterion
-    C = mps.C[0]
-    space_C_old = _firstspace(C_old)
-    space_C = _firstspace(C)
-    if space_C != space_C_old
-        smallest = infimum(space_C_old, space_C)
-        e1 = isometry(space_C_old, smallest)
-        e2 = isometry(space_C, smallest)
-        ϵ = norm(e2' * C * e2 - e1' * C_old * e1)
-    else
-        ϵ = norm(C - C_old)
-    end
+    ϵ = bond_change(C_old, mps.C[0])
 
     # New energy
     ΔE = (E_new - state.energy) / 2
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 9e7b144ef..2d1d8fd5e 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -67,7 +67,7 @@ function dominant_eigsolve(
         operator, mps, alg::VUMPS, envs = environments(mps, operator, mps, alg.alg_environments);
         which
     )
-    log = IterLog("VUMPS")
+    log = IterLog(alg)
     timeroutput = alg.verbosity > 3 ? TimerOutput("VUMPS") : NoTimerOutput()
     iter = 0
 
@@ -106,7 +106,7 @@ function dominant_eigsolve(
     return result
 end
 
-function Base.iterate(it::IterativeSolver{<:VUMPS}, state = it.state)
+function Base.iterate(it::IterativeSolver{<:VUMPS}, state::VUMPSState)
     timeroutput = state.timeroutput
     ACs = @timeit timeroutput "localupdate (parallel)" localupdate_step!(it, state)
     mps = @timeit timeroutput "gauge" gauge_step!(it, state, ACs)
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index bd72f892c..0704743e6 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -38,6 +38,8 @@ Used as the `algorithm` argument of [`propagator`](@ref).
     backend::B = Defaults.backend()
 end
 
+IterativeLoggers.IterLog(::DynamicalDMRG) = IterLog("DDMRG")
+
 """
     propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ)
 
@@ -69,6 +71,57 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac
 """
 struct NaiveInvert <: DDMRG_Flavour end
 
+# Internal state of the dynamical DMRG sweeps
+struct DDMRGState{S, T, Z, O, E, A}
+    mps::S
+    target::T
+    z::Z
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+# dynamical DMRG sweeps over the sites like single-site DMRG
+_sweep_ranges(::DynamicalDMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
+
+function Base.iterate(it::IterativeSolver{<:DynamicalDMRG}, state::DDMRGState)
+    iter = state.iter + 1
+    state = DDMRGState(
+        state.mps, state.target, state.z, state.operator, state.envs, state.iter, 0.0,
+        state.allocator,
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    it.state = DDMRGState(
+        state.mps, state.target, state.z, state.operator, state.envs, iter, state.ϵ,
+        state.allocator,
+    )
+    return (state.mps, state.ϵ), it.state
+end
+
+function _propagator_sweeps!(alg::DynamicalDMRG, state::DDMRGState)
+    log = IterLog(alg)
+    it = IterativeSolver(alg, state)
+
+    with_verbosity(; alg.verbosity) do
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, ϵ) in Iterators.take(it, alg.maxiter)
+            if ϵ <= alg.tol
+                @log_convergence logfinish!(log, it.iter, ϵ)
+                break
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
+            else
+                @log_iteration logiter!(log, it.iter, ϵ)
+            end
+        end
+    end
+
+    return it.state
+end
+
 function propagator(
         A::AbstractFiniteMPS, z::Number, H,
         alg::DynamicalDMRG{NaiveInvert}; init = copy(A)
@@ -77,41 +130,39 @@ function propagator(
     h_envs = environments(init, H, init) # environments for h
     mixedenvs = environments(init, A) # environments for <init | A>
 
-    ϵ = 2 * alg.tol
-    log = IterLog("DDMRG")
+    state = DDMRGState(init, A, z, H, (h_envs, mixedenvs), 0, 2 * alg.tol, allocator)
+    _propagator_sweeps!(alg, state)
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for iter in 1:(alg.maxiter)
-            ϵ = 0.0
+    return dot(A, init), init
+end
 
-            for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
+function sweep!(
+        it::IterativeSolver{<:DynamicalDMRG{NaiveInvert}}, state::DDMRGState, direction, iter
+    )
+    alg = it.alg
+    (; mps, target, z, operator, allocator) = state
+    init, A, H = mps, target, operator
+    h_envs, mixedenvs = state.envs
+    fwd, bwd = _sweep_ranges(alg, A)
+    ϵ = state.ϵ
 
-                H_AC = AC_hamiltonian(i, init, H, init, h_envs; alg.backend, allocator)
-                AC = init.AC[i]
-                AC′, convhist = linsolve(H_AC, -tos, AC, alg.solver, -z, one(z))
+    for i in (direction === Val(:right) ? fwd : bwd)
+        tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
 
-                ϵ = max(ϵ, norm(AC′ - AC))
-                init.AC[i] = AC′
+        H_AC = AC_hamiltonian(i, init, H, init, h_envs; alg.backend, allocator)
+        AC = init.AC[i]
+        AC′, convhist = linsolve(H_AC, -tos, AC, alg.solver, -z, one(z))
 
-                convhist.converged == 0 &&
-                    @warn "propagator ($i) failed to converge: normres = $(convhist.normres)"
-            end
+        ϵ = max(ϵ, norm(AC′ - AC))
+        init.AC[i] = AC′
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        convhist.converged == 0 &&
+            @warn "propagator ($i) failed to converge: normres = $(convhist.normres)"
     end
 
-    return dot(A, init), init
+    return DDMRGState(
+        init, A, z, operator, state.envs, state.iter, ϵ, allocator,
+    )
 end
 
 """
@@ -154,39 +205,8 @@ function propagator(
     H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2
     mixedenvs = environments(init, A) # environments for <init | A>
 
-    ϵ = 2 * alg.tol
-    log = IterLog("DDMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for iter in 1:(alg.maxiter)
-            ϵ = 0.0
-
-            for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
-                H1_AC = AC_hamiltonian(i, init, H, init, envs1; alg.backend, allocator)
-                H2_AC = AC_hamiltonian(i, init, H2, init, envs2; alg.backend, allocator)
-                H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1))
-                AC′, convhist = linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1)
-
-                ϵ = max(ϵ, norm(AC′ - init.AC[i]))
-                init.AC[i] = AC′
-
-                convhist.converged == 0 &&
-                    @warn "propagator ($i) failed to converge: normres $(convhist.normres)"
-            end
-
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
-    end
+    state = DDMRGState(init, A, z, (H, H2), (envs1, envs2, mixedenvs), 0, 2 * alg.tol, allocator)
+    _propagator_sweeps!(alg, state)
 
     a = dot(AC_projection(1, init, A, mixedenvs; alg.backend, allocator), init.AC[1])
     cb = leftenv(envs1, 1, A) * TransferMatrix(init.AL, H[1:length(A.AL)], A.AL)
@@ -200,6 +220,38 @@ function propagator(
     return v, init
 end
 
+function sweep!(
+        it::IterativeSolver{<:DynamicalDMRG{Jeckelmann}}, state::DDMRGState, direction, iter
+    )
+    alg = it.alg
+    (; mps, target, z, allocator) = state
+    init, A = mps, target
+    H, H2 = state.operator
+    envs1, envs2, mixedenvs = state.envs
+    ω = real(z)
+    η = imag(z)
+    fwd, bwd = _sweep_ranges(alg, A)
+    ϵ = state.ϵ
+
+    for i in (direction === Val(:right) ? fwd : bwd)
+        tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
+        H1_AC = AC_hamiltonian(i, init, H, init, envs1; alg.backend, allocator)
+        H2_AC = AC_hamiltonian(i, init, H2, init, envs2; alg.backend, allocator)
+        H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1))
+        AC′, convhist = linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1)
+
+        ϵ = max(ϵ, norm(AC′ - init.AC[i]))
+        init.AC[i] = AC′
+
+        convhist.converged == 0 &&
+            @warn "propagator ($i) failed to converge: normres $(convhist.normres)"
+    end
+
+    return DDMRGState(
+        init, A, z, state.operator, state.envs, state.iter, ϵ, allocator,
+    )
+end
+
 function squaredenvs(
         state::AbstractFiniteMPS, H, envs = environments(state, H, state)
     )
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
index 8e5a9b0f0..fbfcea38c 100644
--- a/src/algorithms/statmech/idmrg.jl
+++ b/src/algorithms/statmech/idmrg.jl
@@ -1,201 +1,193 @@
 function leading_boundary(
-        ψ::InfiniteMultilineMPS, operator, alg::IDMRG, envs = environments(ψ, operator, ψ)
+        ψ::InfiniteMultilineMPS, operator, alg::Union{IDMRG, IDMRG2},
+        envs = environments(ψ, operator, ψ)
     )
     allocator = default_allocator(ψ, SerialScheduler())
-    log = IterLog("IDMRG")
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
+    alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
+    log = IterLog(alg)
+    ϵ_truncs = alg isa IDMRG2 ?
+        PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing
+    state = IDMRGState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator)
+    it = IterativeSolver(alg, state)
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ, leading_eigenvalue(ψ, operator, envs))
-        for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
-            C_current = ψ.C[:, 0]
-
-            # left to right sweep
-            for col in 1:width(ψ)
-                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    ac = ψ.AC[row, col]
-                    (col == width(ψ)) && (ac = copy(ac)) # needed in next sweep
-                    ψ.AL[row, col], ψ.C[row, col] = left_orth!(ac)
-                end
-
-                transfer_leftenv!(envs, ψ, operator, ψ, col + 1)
-            end
-
-            # right to left sweep
-            for col in width(ψ):-1:1
-                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    ψ.C[row, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row, col]; copy = true))
-                    ψ.AR[row, col] = _transpose_front(temp)
-                end
-
-                transfer_rightenv!(envs, ψ, operator, ψ, col - 1)
-            end
-
-            normalize!(envs, ψ, operator, ψ)
-
-            ϵ = norm(C_current - ψ.C[:, 0])
-
+        @log_initialization loginit!(log, it.ϵ, _boundary_objective(alg, ψ, operator, envs))
+        for (ψ, envs, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+                @log_convergence logfinish!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
             else
-                @log_iteration logiter!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+                @log_iteration logiter!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
             end
         end
     end
 
+    (; iter, ϵ) = it.state
     alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
 
     recalculate!(envs, ψ, operator, ψ)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
+    info = AlgorithmInfo(;
+        converged = ϵ <= alg.tol, bondresidual = ϵ,
+        truncation_errors = isnothing(ϵ_truncs) ? nothing : parent(ϵ_truncs), numiter = iter,
+    )
+    return ψ, envs, info
 end
 
-function leading_boundary(
-        ψ::InfiniteMultilineMPS, operator, alg::IDMRG2, envs = environments(ψ, operator, ψ)
+# only single-site IDMRG reports the leading eigenvalue in its log
+# this bypasses issues with intermediate IDMRG2 states not being correct InfiniteMPS
+_boundary_objective(::IDMRG, ψ, operator, envs) = leading_eigenvalue(ψ, operator, envs)
+_boundary_objective(::IDMRG2, ψ, operator, envs) = nothing
+
+function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGState{<:MultilineMPS})
+    iter = state.iter + 1
+    C_current = state.mps.C[:, 0]
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ϵ = bond_change(C_current, state.mps.C[:, 0])
+    it.state = IDMRGState(
+        state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.energy,
+        state.timeroutput, state.allocator,
     )
-    allocator = default_allocator(ψ, SerialScheduler())
-    width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    ϵ::Float64 = 2 * alg.tol
-    log = IterLog("IDMRG2")
-    iter = 0
-    ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ)))
+    return (it.state.mps, it.state.envs, it.state.ϵ), it.state
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
-            C_current = ψ.C[:, 0]
-
-            # sweep from left to right
-            for site in 1:(width(ψ) - 1)
-                ac2 = AC2(ψ, site; kind = :ACAR)
-                h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
-                end
-
-                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
-                transfer_rightenv!(envs, ψ, operator, ψ, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS}, ::Val{:right}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    for col in 1:width(ψ)
+        Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            ac = ψ.AC[row, col]
+            (col == width(ψ)) && (ac = copy(ac)) # needed in next sweep
+            ψ.AL[row, col], ψ.C[row, col] = left_orth!(ac)
+        end
 
-            normalize!(envs, ψ, operator, ψ)
+        transfer_leftenv!(envs, ψ, operator, ψ, col + 1)
+    end
+    return state
+end
 
-            # update the edge
-            site = width(ψ)
-            ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            ac2 = AC2(ψ, site; kind = :ALAC)
-            h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS}, ::Val{:left}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    for col in width(ψ):-1:1
+        Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
 
-            for row in 1:size(ψ, 1)
-                al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
+        for row in 1:size(ψ, 1)
+            ψ.C[row, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row, col]; copy = true))
+            ψ.AR[row, col] = _transpose_front(temp)
+        end
 
-                ψ.AL[row + 1, site] = al
-                ψ.C[row + 1, site] = complex(c)
-                ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+        transfer_rightenv!(envs, ψ, operator, ψ, col - 1)
+    end
+    normalize!(envs, ψ, operator, ψ)
+    return state
+end
 
-                ψ.AC[row + 1, site] = _mul_tail(al, c)
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-                ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS}, ::Val{:right}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    ϵ_truncs = state.truncation_errors
+    for site in 1:(width(ψ) - 1)
+        ac2 = AC2(ψ, site; kind = :ACAR)
+        h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+            ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
+        end
 
-            # TODO: decide if we should compare at the half-sweep level?
-            # C_current = ψ.C[:, site]
+        transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+        transfer_rightenv!(envs, ψ, operator, ψ, site)
+    end
 
-            transfer_leftenv!(envs, ψ, operator, ψ, 1)
-            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    normalize!(envs, ψ, operator, ψ)
 
-            # sweep from right to left
-            for site in reverse(1:(width(ψ) - 1))
-                ac2 = AC2(ψ, site; kind = :ALAC)
-                h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+    # update the edge
+    site = width(ψ)
+    ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    ac2 = AC2(ψ, site; kind = :ALAC)
+    h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+    _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
 
-                for row in 1:size(ψ, 1)
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
+    for row in 1:size(ψ, 1)
+        al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
 
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                end
+        ψ.AL[row + 1, site] = al
+        ψ.C[row + 1, site] = complex(c)
+        ψ.AR[row + 1, site + 1] = _transpose_front(ar)
 
-                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
-                transfer_rightenv!(envs, ψ, operator, ψ, site)
-            end
+        ψ.AC[row + 1, site] = _mul_tail(al, c)
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+        ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+    end
 
-            normalize!(envs, ψ, operator, ψ)
+    transfer_leftenv!(envs, ψ, operator, ψ, 1)
+    transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    return state
+end
 
-            # update the edge
-            ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
-            ac2 = AC2(ψ, 0; kind = :ACAR)
-            h = AC2_hamiltonian(0, ψ, operator, ψ, envs; alg.backend, allocator)
-            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS}, ::Val{:left}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    ϵ_truncs = state.truncation_errors
+    for site in reverse(1:(width(ψ) - 1))
+        ac2 = AC2(ψ, site; kind = :ALAC)
+        h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+        end
 
-            for row in 1:size(ψ, 1)
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
+        transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+        transfer_rightenv!(envs, ψ, operator, ψ, site)
+    end
 
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
+    normalize!(envs, ψ, operator, ψ)
 
-                ψ.AR[row + 1, end] = _transpose_front(
-                    ψ.C[row + 1, end - 1] \ _transpose_tail(al * c)
-                )
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-            end
+    # update the edge
+    ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
+    ac2 = AC2(ψ, 0; kind = :ACAR)
+    h = AC2_hamiltonian(0, ψ, operator, ψ, envs; alg.backend, allocator)
+    _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
 
-            transfer_leftenv!(envs, ψ, operator, ψ, 1)
-            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    for row in 1:size(ψ, 1)
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
 
-            # update error
-            ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
-                smallest = infimum(_firstspace(c1), _firstspace(c2))
-                e1 = isometry(_firstspace(c1), smallest)
-                e2 = isometry(_firstspace(c2), smallest)
-                return norm(e2' * c2 * e2 - e1' * c1 * e1)
-            end
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AR[row + 1, end] = _transpose_front(
+            ψ.C[row + 1, end - 1] \ _transpose_tail(al * c)
+        )
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
     end
 
-    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
-
-    recalculate!(envs, ψ, operator, ψ)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), numiter = iter)
+    transfer_leftenv!(envs, ψ, operator, ψ, 1)
+    transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    return state
 end
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 55c006cef..7ea87bb6b 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -61,7 +61,7 @@ function dominant_eigsolve(
         which
     )
     @assert which === :LM "VOMPS only supports the LM eigenvalue problem"
-    log = IterLog("VOMPS")
+    log = IterLog(alg)
     iter = 0
     ϵ = calc_galerkin(mps, operator, mps, envs; alg.backend)
     alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
@@ -90,7 +90,7 @@ function dominant_eigsolve(
     end
 end
 
-function Base.iterate(it::IterativeSolver{<:VOMPS}, state)
+function Base.iterate(it::IterativeSolver{<:VOMPS}, state::VOMPSState)
     ACs = localupdate_step!(it, state)
     mps = gauge_step!(it, state, ACs)
     envs = envs_step!(it, state, mps)
diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl
index 24df2c962..f9e79e199 100644
--- a/src/algorithms/timestep/bug.jl
+++ b/src/algorithms/timestep/bug.jl
@@ -133,10 +133,12 @@ function _evolve_center(site, ψ, H, alg::BUG, envs, t, h, allocator; imaginary_
     return integrate(Heff, ψ.AC[site], t, h, alg.integrator; imaginary_evolution)
 end
 
-function local_update!(
-        site, direction::Val, ψ, H, alg::BUG, envs, t, h, allocator;
-        imaginary_evolution, normalize, timeroutput
-    )
+# a half-sweep evolves over `h = dt / 2`
+function local_update!(site, direction::Val, ψ, alg::BUG, state::TimeStepState)
+    (; operator, envs, t, allocator, timeroutput, imaginary_evolution, normalize) = state
+    H = operator
+    h = state.dt / 2
+
     # at the far end of the sweep there is no bond ahead to cut: evolve and finalize
     if site == _sweep_end(ψ, direction)
         AC = @timeit timeroutput "AC_integrate" _evolve_center(
@@ -165,57 +167,8 @@ function local_update!(
     return ψ, ϵ
 end
 
-function timestep!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG,
-        envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
-        imaginary_evolution::Bool = false, normalize::Bool = false,
-        timeroutput = NoTimerOutput()
-    )
-    L = length(ψ)
-    h = dt / 2
-
-    # the sweep is serial, so a single allocator serves all local updates
-    allocator = default_allocator(ψ, SerialScheduler())
-
-    ϵ_truncs = zeros(real(scalartype(ψ)), L - 1)
-
-    # left→right half-sweep (root = last site): `t → t + dt / 2`
-    @timeit timeroutput "half-sweep" for site in 1:L
-        ψ, ϵ = local_update!(
-            site, Val(:right), ψ, H, alg, envs, t, h, allocator;
-            imaginary_evolution, normalize, timeroutput
-        )
-        site < L && (ϵ_truncs[site] = ϵ)
-    end
-
-    # right→left half-sweep (root = first site): `t + dt / 2 → t + dt`
-    @timeit timeroutput "half-sweep" for site in L:-1:1
-        ψ, ϵ = local_update!(
-            site, Val(:left), ψ, H, alg, envs, t + h, h, allocator;
-            imaginary_evolution, normalize, timeroutput
-        )
-        site > 1 && (ϵ_truncs[site - 1] = ϵ)
-    end
-
-    return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs)
-end
+# like single-site TDVP, BUG visits every site and only evolves the center at the far end
+_sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1)
 
-# copying version
-function timestep(
-        ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,
-        alg::BUG, envs::AbstractMPSEnvironments...;
-        imaginary_evolution::Bool = false, normalize::Bool = false, kwargs...
-    )
-    isreal = (scalartype(ψ) <: Real && !(imaginary_evolution && scalartype(H) <: Real))
-    ψ′ = isreal ? complex(ψ) : copy(ψ)
-    if length(envs) != 0 && isreal
-        @warn "Currently cannot reuse real environments for complex evolution"
-        envs′ = environments(ψ′, H, ψ′)
-    elseif length(envs) == 1
-        envs′ = only(envs)
-    else
-        @assert length(envs) == 0 "Invalid signature"
-        envs′ = environments(ψ′, H, ψ′)
-    end
-    return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...)
-end
+_promotes_to_complex(ψ, H, dt, ::BUG; imaginary_evolution) =
+    scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt) && scalartype(H) <: Real)
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index af4707f2d..1c0c126a3 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -72,7 +72,56 @@ function TDVP(;
     )
 end
 
-function timestep(
+# Internal state of a single time step `t → t + dt`
+struct TimeStepState{S, O, E, T, D, V, SC, A, TO}
+    mps::S
+    operator::O
+    envs::E
+    t::T
+    dt::D
+    truncation_errors::V
+    scheduler::SC
+    allocator::A
+    timeroutput::TO
+    imaginary_evolution::Bool
+    normalize::Bool
+end
+
+# a single allocator serves all local updates of a step, and it can be carried over to later
+# steps together with the scheduler that selected it
+function TimeStepState(
+        ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput(),
+        scheduler = _step_scheduler(ψ), allocator = default_allocator(ψ, scheduler)
+    )
+    truncation_errors = _truncation_errors(ψ)
+    return TimeStepState(
+        ψ, H, envs, t, dt, truncation_errors, scheduler, allocator, timeroutput,
+        imaginary_evolution, normalize,
+    )
+end
+
+# a finite sweep is serial, whereas an infinite step updates all sites concurrently
+_step_scheduler(::AbstractFiniteMPS) = SerialScheduler()
+_step_scheduler(::InfiniteMPS) = Defaults.scheduler[]
+
+# per-bond truncation errors of a finite sweep; an infinite step never truncates
+_truncation_errors(ψ::AbstractFiniteMPS) = zeros(real(scalartype(ψ)), length(ψ) - 1)
+_truncation_errors(ψ::InfiniteMPS) = nothing
+
+# infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't report
+# `truncation_errors` (rather than report zeros that would look like a measurement)
+_step_info(state::TimeStepState) = isnothing(state.truncation_errors) ? AlgorithmInfo() :
+    AlgorithmInfo(; truncation_errors = state.truncation_errors)
+
+function _with_mps(state::TimeStepState, ψ; t = state.t)
+    return TimeStepState(
+        ψ, state.operator, state.envs, t, state.dt, state.truncation_errors,
+        state.scheduler, state.allocator, state.timeroutput, state.imaginary_evolution,
+        state.normalize,
+    )
+end
+
+function timestep!(
         ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP,
         envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
         leftorthflag = true, imaginary_evolution::Bool = false, normalize::Bool = false
@@ -80,26 +129,38 @@ function timestep(
     # `normalize` is accepted for signature uniformity with the finite integrators, but an
     # `InfiniteMPS` is always normalized to norm-1-per-site by the gauge/reconstruction below
     # (a structural gauge requirement, not information erasure), so the flag has no effect here.
-    # convert state to complex if necessary
-    if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt))
-        return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution, normalize)
-    end
-
     # the scheduler is read here rather than below, so that the allocator it selects is inferable
-    return _timestep_infinite(
+    return _timestep_infinite!(
         ψ, H, t, dt, alg, envs, Defaults.scheduler[]; leftorthflag, imaginary_evolution
     )
 end
 
-function _timestep_infinite(
+function _timestep_infinite!(
         ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, scheduler::Scheduler;
         leftorthflag, imaginary_evolution
     )
+    state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize = false, scheduler)
+    state = _timestep!(alg, state; leftorthflag)
+    # the gauge-fixing residual is controlled by `tolgauge`, not reported here
+    return state.mps, state.envs, _step_info(state)
+end
+
+# an infinite step evolves every center tensor at once, then restores the uniform gauge and
+# recomputes the environments
+function _timestep!(alg::TDVP, state::TimeStepState{<:InfiniteMPS}; leftorthflag = true)
+    ACs, Cs = localupdate_step!(alg, state)
+    ψ = gauge_step!(alg, state, ACs, Cs; leftorthflag)
+    recalculate!(state.envs, ψ, state.operator)
+    return _with_mps(state, ψ; t = state.t + state.dt)
+end
+
+function localupdate_step!(alg::TDVP, state::TimeStepState)
+    (; mps, operator, envs, t, dt, scheduler, allocator, imaginary_evolution) = state
+    ψ, H = mps, operator
     temp_ACs = similar(ψ.AC)
     temp_Cs = similar(ψ.C)
 
     # both sweeps together are a single unit of concurrent work, and share one allocator
-    allocator = default_allocator(ψ, scheduler)
     ac_sweep!() = tforeach(1:length(ψ); scheduler) do loc
         Hac = AC_hamiltonian(loc, ψ, H, ψ, envs; alg.backend, allocator)
         temp_ACs[loc] = integrate(Hac, ψ.AC[loc], t, dt, alg.integrator; imaginary_evolution)
@@ -121,101 +182,90 @@ function _timestep_infinite(
             Threads.@spawn c_sweep!()
         end
     end
+    return temp_ACs, temp_Cs
+end
 
+# the evolved center tensors are split into the new left (or right) isometries, from which the
+# uniform gauge is restored in place
+function gauge_step!(alg::TDVP, state::TimeStepState, ACs, Cs; leftorthflag)
+    ψ = state.mps
     if leftorthflag
-        regauge!.(temp_ACs, temp_Cs)
-        ψ′ = InfiniteMPS(temp_ACs, ψ.C[end]; tol = alg.tolgauge, maxiter = alg.gaugemaxiter)
+        regauge!.(ACs, Cs)
+        copyto!(ψ.AL, ACs)
+        gaugefix!(ψ, ψ.AL, ψ.C[end]; order = :R, tol = alg.tolgauge, maxiter = alg.gaugemaxiter)
     else
-        circshift!(temp_Cs, 1)
-        regauge!.(temp_Cs, temp_ACs)
-        ψ′ = InfiniteMPS(ψ.C[0], temp_ACs; tol = alg.tolgauge, maxiter = alg.gaugemaxiter)
+        circshift!(Cs, 1)
+        regauge!.(Cs, ACs)
+        copyto!(ψ.AR, ACs)
+        gaugefix!(ψ, ψ.AR, ψ.C[0]; order = :L, tol = alg.tolgauge, maxiter = alg.gaugemaxiter)
     end
+    mul!.(ψ.AC, ψ.AL, ψ.C)
+    return ψ
+end
 
-    recalculate!(envs, ψ′, H)
-    # infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't
-    # report `truncation_errors` (rather than report zeros that would look like a measurement)
-    # the gauge-fixing residual is controlled by `tolgauge`, not reported here
-    return ψ′, envs, AlgorithmInfo()
+# Single-site integrators visit every site, finishing a half-sweep by evolving only the center
+# tensor at its far end.
+_sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1)
+
+# bond cut by the update at `site`, or `nothing` at the far end of a single-site half-sweep
+_swept_bond(alg, site, ψ, ::Val{:right}) = site < length(ψ) ? site : nothing
+_swept_bond(alg, site, ψ, ::Val{:left}) = site > 1 ? site - 1 : nothing
+
+function sweep!(alg, state::TimeStepState, direction)
+    fwd, bwd = _sweep_ranges(alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ = state.mps
+    @timeit state.timeroutput "half-sweep" for site in sites
+        ψ, ϵ = local_update!(site, direction, ψ, alg, state)
+        bond = _swept_bond(alg, site, ψ, direction)
+        isnothing(bond) || (state.truncation_errors[bond] = ϵ)
+    end
+    # each half-sweep evolves the state over `dt / 2`
+    return _with_mps(state, ψ; t = state.t + state.dt / 2)
 end
 
-function timestep!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP,
-        envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
-        imaginary_evolution::Bool = false, normalize::Bool = false
-    )
-    # the sweep is serial, so a single allocator serves all local updates
-    allocator = default_allocator(ψ, SerialScheduler())
-    return _timestep_finite!(
-        ψ, H, t, dt, alg, envs, allocator; imaginary_evolution, normalize
-    )
+# a step is a left→right half-sweep `t → t + dt / 2` followed by a right→left half-sweep
+# `t + dt / 2 → t + dt`
+function _timestep!(alg, state::TimeStepState{<:AbstractFiniteMPS})
+    state = sweep!(alg, state, Val(:right))
+    return sweep!(alg, state, Val(:left))
 end
 
-function _timestep_finite!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator;
-        imaginary_evolution::Bool, normalize::Bool
-    )
-    ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1)
-
-    # sweep left to right
-    for i in 1:(length(ψ) - 1)
-        # 1. optionally expand the bond ahead of the local update (CBE)
-        isnothing(alg.alg_expand) ||
-            changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize, allocator)
-
-        # 2. evolve the (possibly expanded) center tensor forward
-        Hac = AC_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution)
-
-        # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the
-        #    enlarged bond back down) and move the center to i+1. By default the norm is
-        #    preserved; `normalize` renormalizes.
-        _, ϵ_truncs[i] = left_gauge!(ψ, i, AC, alg.alg_gauge; normalize)
-
-        # 4. evolve the bond tensor backward
-        Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        ψ.C[i] = integrate(
-            Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator;
-            imaginary_evolution
-        )
+function local_update!(site, direction::Val, ψ, alg::TDVP, state::TimeStepState)
+    (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state
+    H = operator
+    t_AC, t_C = t, t + dt / 2
+
+    # at the far end of the sweep there is no bond ahead: only evolve the center tensor
+    if site == _sweep_end(ψ, direction)
+        Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+        ψ.AC[site] = integrate(Hac, ψ.AC[site], t_AC, dt / 2, alg.integrator; imaginary_evolution)
+        return ψ, zero(real(scalartype(ψ)))
     end
 
-    # edge case
-    Hac = AC_hamiltonian(length(ψ), ψ, H, ψ, envs; alg.backend, allocator)
-    ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator; imaginary_evolution)
-
-    # sweep right to left
-    for i in length(ψ):-1:2
-        # 1. optionally expand the bond ahead of the local update (CBE)
-        isnothing(alg.alg_expand) ||
-            changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize, allocator)
-
-        # 2. evolve the (possibly expanded) center tensor forward
-        Hac = AC_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        AC = integrate(
-            Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator;
-            imaginary_evolution
-        )
-
-        # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (norm preserved by
-        #    default; `normalize` renormalizes)
-        _, ϵ_truncs[i - 1] = right_gauge!(ψ, i, AC, alg.alg_gauge; normalize)
-
-        # 4. evolve the bond tensor backward
-        Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-        ψ.C[i - 1] = integrate(
-            Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator;
-            imaginary_evolution
-        )
+    # 1. optionally expand the bond ahead of the local update (CBE)
+    isnothing(alg.alg_expand) ||
+        changebond!(site, direction, ψ, H, alg.alg_expand, envs; normalize, allocator)
+
+    # 2. evolve the (possibly expanded) center tensor forward
+    Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+    AC = integrate(Hac, ψ.AC[site], t_AC, dt / 2, alg.integrator; imaginary_evolution)
+
+    # 3. gauge: split AC onto the bond ahead and move the center.
+    #    By default the norm is preserved; `normalize` renormalizes.
+    if direction === Val(:right)
+        _, ϵ = left_gauge!(ψ, site, AC, alg.alg_gauge; normalize)
+        bond = site
+    else
+        _, ϵ = right_gauge!(ψ, site, AC, alg.alg_gauge; normalize)
+        bond = site - 1
     end
 
-    # edge case
-    Hac = AC_hamiltonian(1, ψ, H, ψ, envs; alg.backend, allocator)
-    ψ.AC[1] = integrate(
-        Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator;
-        imaginary_evolution
-    )
+    # 4. evolve the bond tensor backward
+    Hc = C_hamiltonian(bond, ψ, H, ψ, envs; alg.backend, allocator)
+    ψ.C[bond] = integrate(Hc, ψ.C[bond], t_C, -dt / 2, alg.integrator; imaginary_evolution)
 
-    return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs)
+    return ψ, ϵ
 end
 
 """
@@ -259,81 +309,39 @@ Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and
     backend::B = Defaults.backend()
 end
 
-function timestep!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2,
-        envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
-        imaginary_evolution::Bool = false, normalize::Bool = false
-    )
-    # the sweep is serial, so a single allocator serves all local updates
-    allocator = default_allocator(ψ, SerialScheduler())
-    return _timestep2_finite!(
-        ψ, H, t, dt, alg, envs, allocator; imaginary_evolution, normalize
-    )
-end
-
-function _timestep2_finite!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs, allocator;
-        imaginary_evolution::Bool, normalize::Bool
-    )
-    # the two-site center always has to be split back up, so the gauge is always a truncated SVD
-    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc)
-
-    ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1)
-
-    # sweep left to right
-    for i in 1:(length(ψ) - 1)
-        ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])
-        Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution)
+# two-site TDVP visits every bond, and each update cuts the bond it acts on
+_sweep_ranges(::TDVP2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 1):-1:1)
+_swept_bond(::TDVP2, pos, ψ, ::Val{:right}) = pos
+_swept_bond(::TDVP2, pos, ψ, ::Val{:left}) = pos
 
-        # the norm of the discarded singular values is the truncation error
-        _, ϵ_truncs[i] = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize)
+function local_update!(pos, direction::Val, ψ, alg::TDVP2, state::TimeStepState)
+    (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state
+    H = operator
+    t_AC2, t_AC = t, t + dt / 2
 
-        if i != (length(ψ) - 1)
-            Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator)
-            ψ.AC[i + 1] = integrate(
-                Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator;
-                imaginary_evolution
-            )
-        end
+    # 1. evolve the two-site center tensor at `(pos, pos + 1)` forward
+    ac2 = if direction === Val(:right)
+        _transpose_front(ψ.AC[pos]) * _transpose_tail(ψ.AR[pos + 1])
+    else
+        _transpose_front(ψ.AL[pos]) * _transpose_tail(ψ.AC[pos + 1])
     end
+    Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs; alg.backend, allocator)
+    ac2′ = integrate(Hac2, ac2, t_AC2, dt / 2, alg.integrator; imaginary_evolution)
 
-    # sweep right to left
-    for i in length(ψ):-1:2
-        ac2 = _transpose_front(ψ.AL[i - 1]) * _transpose_tail(ψ.AC[i])
-        Hac2 = AC2_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-        ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution)
-
-        _, ϵ_truncs[i - 1] = gauge2!(ψ, i - 1, Val(:left), ac2′, alg_gauge; normalize)
+    # 2. gauge: the two-site center always has to be split back up, so this is always a
+    #    truncated SVD, and the norm of the discarded singular values is the truncation error
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc)
+    _, ϵ = gauge2!(ψ, pos, direction, ac2′, alg_gauge; normalize)
 
-        if i != 2
-            Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-            ψ.AC[i - 1] = integrate(
-                Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator;
-                imaginary_evolution
-            )
-        end
+    # 3. evolve the new single-site center backward, except at the far end of the sweep
+    if direction === Val(:right) ? pos != length(ψ) - 1 : pos != 1
+        site = direction === Val(:right) ? pos + 1 : pos
+        Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+        ψ.AC[site] = integrate(Hac, ψ.AC[site], t_AC, -dt / 2, alg.integrator; imaginary_evolution)
     end
 
-    return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs)
+    return ψ, ϵ
 end
 
-# copying version
-function timestep(
-        ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,
-        alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...;
-        imaginary_evolution::Bool = false, normalize::Bool = false, kwargs...
-    )
-    isreal = (scalartype(ψ) <: Real && !imaginary_evolution)
-    ψ′ = isreal ? complex(ψ) : copy(ψ)
-    if length(envs) != 0 && isreal
-        @warn "Currently cannot reuse real environments for complex evolution"
-        envs′ = environments(ψ′, H, ψ′)
-    elseif length(envs) == 1
-        envs′ = only(envs)
-    else
-        @assert length(envs) == 0 "Invalid signature"
-        envs′ = environments(ψ′, H, ψ′)
-    end
-    return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...)
-end
+_promotes_to_complex(ψ, H, dt, ::Union{TDVP, TDVP2}; imaginary_evolution) =
+    scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt))
diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl
index 9d9e87f40..866bb8230 100644
--- a/src/algorithms/timestep/time_evolve.jl
+++ b/src/algorithms/timestep/time_evolve.jl
@@ -36,43 +36,94 @@ evolution at `verbosity ≥ 2`.
 """
 function time_evolve end, function time_evolve! end
 
+# Internal state of `time_evolve`
+struct TimeEvolveState{S, O, E, T, F, K, P}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    t_span::T
+    stepper::F
+    kwargs::K
+    step::P
+end
+
+# one iteration is a single time step followed by `finalize`, which receives the step's start time
+function Base.iterate(it::IterativeSolver{<:Any, <:TimeEvolveState}, state::TimeEvolveState)
+    state.iter >= length(state.t_span) - 1 && return nothing
+    iter = state.iter + 1
+    t = state.t_span[iter]
+    dt = state.t_span[iter + 1] - t
+
+    ψ, envs, info, step = _evolve_step(it.alg, state, t, dt)
+    ψ, envs = it.finalize(t, ψ, state.operator, envs)::Tuple{typeof(ψ), typeof(envs)}
+
+    it.state = TimeEvolveState(
+        ψ, state.operator, envs, iter, state.t_span, state.stepper, state.kwargs, step
+    )
+    return (ψ, envs, info), it.state
+end
+
+# only the first step copies the state, if `stepper` is `timestep`
+function _evolve_step(alg, state::TimeEvolveState, t, dt)
+    (; imaginary_evolution, normalize) = state.kwargs
+    H = state.operator
+    step = if isnothing(state.step)
+        ψ, envs = state.stepper === timestep ?
+            _timestep_copy(
+                state.mps, H, (state.envs,),
+                _promotes_to_complex(state.mps, H, dt, alg; imaginary_evolution)
+            ) : (state.mps, state.envs)
+        TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize)
+    else
+        TimeStepState(
+            state.mps, H, state.envs, t, dt; imaginary_evolution, normalize,
+            state.step.scheduler, state.step.allocator,
+        )
+    end
+    step = _timestep!(alg, step)
+    return step.mps, step.envs, _step_info(step), step
+end
+
 for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!))
     @eval function $time_evolve(
             ψ, H, t_span::AbstractVector{<:Number}, alg,
             envs = environments(ψ, H, ψ);
             verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false
         )
-        log = IterLog(string(nameof(typeof(alg))))
-        truncation_errors = []
-        ϵ_max = 0.0
-        with_verbosity(; verbosity) do
-            @log_initialization loginit!(log, 0.0, first(t_span))
-            for iter in 1:(length(t_span) - 1)
-                t = t_span[iter]
-                dt = t_span[iter + 1] - t
-
-                ψ, envs, info_step = $timestep(
-                    ψ, H, t, dt, alg, envs; imaginary_evolution, normalize
-                )
-                ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)}
-
-                # the log shows the largest per-bond error, or zero for a non-truncating algorithm
-                ϵ_step = 0.0
-                if haskey(info_step, :truncation_errors)
-                    push!(truncation_errors, info_step.truncation_errors)
-                    ϵ_step = Float64(maximum(info_step.truncation_errors; init = 0.0))
-                end
-                ϵ_max = max(ϵ_max, ϵ_step)
-                @log_iteration logiter!(log, iter, ϵ_step, t)
+        state = TimeEvolveState(
+            ψ, H, envs, 0, t_span, $timestep, (; imaginary_evolution, normalize), nothing
+        )
+        return _time_evolve(alg, state; verbosity)
+    end
+end
+
+function _time_evolve(alg, state::TimeEvolveState; verbosity::Int)
+    log = IterLog(alg)
+    # the state type is left abstract, since the first step may promote a real state to complex
+    it = IterativeSolver{typeof(alg), TimeEvolveState}(alg, state)
+    t_span = state.t_span
+    truncation_errors = []
+    ϵ_max = 0.0
+    with_verbosity(; verbosity) do
+        @log_initialization loginit!(log, 0.0, first(t_span))
+        for (_, _, info_step) in it
+            # the log shows the largest per-bond error, or zero for a non-truncating algorithm
+            ϵ_step = 0.0
+            if haskey(info_step, :truncation_errors)
+                push!(truncation_errors, info_step.truncation_errors)
+                ϵ_step = Float64(maximum(info_step.truncation_errors; init = 0.0))
             end
-            @log_convergence logfinish!(log, length(t_span), ϵ_max, t_span[end])
+            ϵ_max = max(ϵ_max, ϵ_step)
+            @log_iteration logiter!(log, it.iter, ϵ_step, t_span[it.iter])
         end
-        info = AlgorithmInfo(;
-            numiter = length(t_span) - 1,
-            truncation_errors = isempty(truncation_errors) ? nothing : copy(truncation_errors)
-        )
-        return ψ, envs, info
+        @log_convergence logfinish!(log, length(t_span), ϵ_max, t_span[end])
     end
+    info = AlgorithmInfo(;
+        numiter = length(t_span) - 1,
+        truncation_errors = isempty(truncation_errors) ? nothing : copy(truncation_errors)
+    )
+    return it.state.mps, it.state.envs, info
 end
 
 """
@@ -142,6 +193,43 @@ julia> round(real(expectation_value(ψ, 2 => X)); digits = 6)
 """
 function timestep end, function timestep! end
 
+function timestep!(
+        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg,
+        envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
+        imaginary_evolution::Bool = false, normalize::Bool = false,
+        timeroutput = NoTimerOutput()
+    )
+    state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput)
+    state = _timestep!(alg, state)
+    return state.mps, state.envs, _step_info(state)
+end
+
+function timestep(
+        ψ::Union{AbstractFiniteMPS, InfiniteMPS}, H, time::Number, timestep::Number,
+        alg, envs::AbstractMPSEnvironments...;
+        imaginary_evolution::Bool = false, normalize::Bool = false, kwargs...
+    )
+    promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution)
+    ψ′, envs′ = _timestep_copy(ψ, H, envs, promote)
+    return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...)
+end
+
+# Copy of `ψ` to step out of place, promoted to complex when the evolution leaves the reals,
+# together with environments for that copy.
+function _timestep_copy(ψ, H, envs, promote::Bool)
+    ψ′ = promote ? complex(ψ) : copy(ψ)
+    if length(envs) != 0 && promote
+        @warn "Currently cannot reuse real environments for complex evolution"
+        envs′ = environments(ψ′, H, ψ′)
+    elseif length(envs) == 1
+        envs′ = only(envs)
+    else
+        @assert length(envs) == 0 "Invalid signature"
+        envs′ = environments(ψ′, H, ψ′)
+    end
+    return ψ′, envs′
+end
+
 @doc """
     make_time_mpo(H::MPOHamiltonian, dt::Number, alg; kwargs...) -> O::MPO
 
diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl
index 353b544f4..f91d06835 100644
--- a/src/states/infinitemps.jl
+++ b/src/states/infinitemps.jl
@@ -239,7 +239,7 @@ function InfiniteMPS(C₀::MPSBondTensor, AR::AbstractVector{<:GenericMPSTensor}
     # initialize tensor storage
     AC = similar.(AR)
     AL = similar.(AR)
-    T = TensorOperations.promote_contract(eltype(AR), eltype(C₀))
+    T = TensorOperations.promote_contract(scalartype(AR), scalartype(C₀))
     TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false)
     C = similar(AR, TC)
     ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC)
diff --git a/src/states/ortho.jl b/src/states/ortho.jl
index d7a53a82d..966428cf4 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -61,6 +61,9 @@ Used as the `alg` argument of [`gaugefix!`](@ref).
     eig_miniter::Int = 10
 end
 
+IterativeLoggers.IterLog(::LeftCanonical) = IterLog("LC")
+IterativeLoggers.IterLog(::RightCanonical) = IterLog("RC")
+
 """
 $(TYPEDEF)
 
@@ -228,7 +231,7 @@ function uniform_leftorth!(
     C[end] = normalize!(C₀)
     return with_verbosity(; alg.verbosity) do
         # initialize algorithm and temporary variables
-        log = IterLog("LC")
+        log = IterLog(alg)
         A_tail = _transpose_tail.(A) # pre-transpose A
         CA_tail = similar.(A_tail)  # pre-allocate workspace
         state = (; AL, C, A, A_tail, CA_tail, iter = 0, ϵ = Inf, timeroutput, backend, allocator)
@@ -250,7 +253,7 @@ function uniform_leftorth!(
     end
 end
 
-function Base.iterate(it::IterativeSolver{LeftCanonical}, state = it.state)
+function Base.iterate(it::IterativeSolver{LeftCanonical}, state::NamedTuple)
     timeroutput = state.timeroutput
     C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state)
     C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state)
@@ -297,7 +300,7 @@ function uniform_rightorth!(
     C[end] = normalize!(C₀)
     return with_verbosity(; alg.verbosity) do
         # initialize algorithm and temporary variables
-        log = IterLog("RC")
+        log = IterLog(alg)
         AC_tail = _similar_tail.(A) # pre-allocate workspace
         state = (; AR, C, A, AC_tail, iter = 0, ϵ = Inf, timeroutput, backend, allocator)
         it = IterativeSolver(alg, state)
@@ -318,7 +321,7 @@ function uniform_rightorth!(
     end
 end
 
-function Base.iterate(it::IterativeSolver{RightCanonical}, state = it.state)
+function Base.iterate(it::IterativeSolver{RightCanonical}, state::NamedTuple)
     timeroutput = state.timeroutput
     C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state)
     C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state)
diff --git a/src/utility/iterativesolvers.jl b/src/utility/iterativesolvers.jl
index 8eb863d12..ce1aea330 100644
--- a/src/utility/iterativesolvers.jl
+++ b/src/utility/iterativesolvers.jl
@@ -19,3 +19,6 @@ function Base.getproperty(it::IterativeSolver{A, B}, name::Symbol) where {A, B}
 end
 
 Base.iterate(it::IterativeSolver) = iterate(it, it.state)
+
+# iteration logs are named after their algorithm unless the algorithm specifies otherwise
+IterativeLoggers.IterLog(alg::Algorithm) = IterLog(string(nameof(typeof(alg))))
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index d527e06c3..88dce8c20 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -40,6 +40,30 @@ end
 _firstspace(t::AbstractTensorMap) = space(t, 1)
 _lastspace(t::AbstractTensorMap) = space(t, numind(t))
 
+"""
+    bond_change(C_old, C_new)
+
+Norm of the change of a bond tensor. If the bond space changed, both tensors are compared on
+their common subspace, i.e. after restricting them with `isometry(V, infimum(V_old, V_new))`.
+For vectors of bond tensors, such as a column of a multiline MPS, the changes are summed.
+"""
+function bond_change(C_old::AbstractTensorMap, C_new::AbstractTensorMap)
+    V_old, V_new = _firstspace(C_old), _firstspace(C_new)
+    V_old == V_new && return norm(C_new - C_old)
+    # the isometry onto a subspace selects the leading rows and columns of every block
+    V = infimum(V_old, V_new)
+    ϵ² = zero(real(scalartype(C_new)))
+    for c in sectors(V)
+        d = dim(V, c)
+        b_old = view(block(C_old, c), 1:d, 1:d)
+        b_new = view(block(C_new, c), 1:d, 1:d)
+        ϵ² += dim(c) * norm(b_new - b_old)^2
+    end
+    return sqrt(ϵ²)
+end
+bond_change(C_old::AbstractVector, C_new::AbstractVector) =
+    sum(splat(bond_change), zip(C_old, C_new))
+
 """
     similar_scalartype(T::Type{<:AbstractTensorMap}, S::Type{<:Number})
 
diff --git a/test/states/infinitemps.jl b/test/states/infinitemps.jl
index 1226c395d..90f03c468 100644
--- a/test/states/infinitemps.jl
+++ b/test/states/infinitemps.jl
@@ -44,6 +44,12 @@ using Adapt
         @test TransferMatrix(ψ.AR[i], ψ.AL[i]) * r_RL(ψ, i) ≈ r_RL(ψ, i + 1)
         @test TransferMatrix(ψ.AR[i], ψ.AR[i]) * r_RR(ψ, i) ≈ r_RR(ψ, i + 1)
     end
+
+    # rebuilding from the left or right isometries and a bond tensor gives back the same state
+    for ψ′ in (InfiniteMPS(ψ.AL, ψ.C[end]; tol), InfiniteMPS(ψ.C[0], ψ.AR; tol))
+        @test scalartype(ψ′) == elt
+        @test abs(dot(ψ, ψ′)) ≈ 1
+    end
 end
 
 @testset "InfiniteMPS copying" begin
diff --git a/test/timeevolution/timestep.jl b/test/timeevolution/timestep.jl
index 3baebcf91..75d282957 100644
--- a/test/timeevolution/timestep.jl
+++ b/test/timeevolution/timestep.jl
@@ -82,11 +82,34 @@ maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 1:length(ψ))
     # the AC and C sweeps of the infinite integrator spawn over the unit cell, and additionally run
     # concurrently with each other, so they share one allocator: check both schedulers agree
     @testset "Infinite TDVP ($schedname)" for (schedname, scheduler) in SCHEDULERS
-        ψ, envs = with_scheduler(scheduler) do
-            return timestep(ψ₀, H, 0.0, dt, TDVP())
+        with_scheduler(scheduler) do
+            ψ, envs = timestep(ψ₀, H, 0.0, dt, TDVP())
+            E = expectation_value(ψ, H, envs)
+            @test E₀ ≈ E atol = 1.0e-2
+
+            ψ′ = copy(ψ₀)
+            ψ″, = timestep!(ψ′, H, 0.0, dt, TDVP())
+            @test ψ″ === ψ′
+            @test ψ″.AC == ψ.AC
+
+            # evolving in the right-orthogonal gauge gives the same state
+            ψr, = timestep(ψ₀, H, 0.0, dt, TDVP(); leftorthflag = false)
+            @test abs(dot(ψ, ψr)) ≈ 1 atol = 1.0e-8
         end
-        E = expectation_value(ψ, H, envs)
-        @test E₀ ≈ E atol = 1.0e-2
+    end
+
+    @testset "Infinite real TDVP" begin
+        Hr = force_planar(transverse_field_ising(Float64; g = 0.7))
+        ψr₀ = InfiniteMPS(rand, Float64, [ℙ^2], [ℙ^8])
+        # real-time evolution promotes the state, and with it the environments, to complex
+        ψ, envs = timestep(ψr₀, Hr, 0.0, dt, TDVP())
+        @test scalartype(ψ) <: Complex
+        ψc, = timestep(complex(ψr₀), Hr, 0.0, dt, TDVP())
+        @test abs(dot(ψ, ψc)) ≈ 1 atol = 1.0e-8
+        @test expectation_value(ψ, Hr, envs) ≈ expectation_value(ψc, Hr)
+        # imaginary-time evolution with a real step keeps it real
+        ψi, = timestep(ψr₀, Hr, 0.0, dt, TDVP(); imaginary_evolution = true)
+        @test scalartype(ψi) == Float64
     end
 
     Hlazy = LazySum([3 * deepcopy(H), 1.55 * deepcopy(H), -0.1 * deepcopy(H)])
@@ -329,5 +352,10 @@ end
         ψ, envs = time_evolve(ψ₀, H, t_span, TDVP(); verbosity = 2)
         E = expectation_value(ψ, H, envs)
         @test E₀ ≈ E atol = 1.0e-2
+
+        ψ′ = copy(ψ₀)
+        ψ″, = MPSKit.time_evolve!(ψ′, H, t_span, TDVP())
+        @test ψ″ === ψ′
+        @test ψ″.AC == ψ.AC
     end
 end