diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 15ff5402c..e86f71f25 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -126,15 +126,35 @@ function changebond!(
         # a sweep that already holds one should pass it: this is called once per site
         allocator = default_allocator(ψ, SerialScheduler()),
     )
-    bond = site
+    # two-site update from the projected effective Hamiltonian
+    AC2 = AC2_projection(site, ψ, H, ψ, envs; alg.backend, allocator)
+    return expand_from_projection!(site, Val(:right), ψ, AC2, alg; normalize)
+end
+function changebond!(
+        site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs;
+        normalize::Bool = true,
+        # a sweep that already holds one should pass it: this is called once per site
+        allocator = default_allocator(ψ, SerialScheduler()),
+    )
+    # two-site update from the projected effective Hamiltonian
+    AC2 = AC2_projection(site - 1, ψ, H, ψ, envs; alg.backend, allocator)
+    return expand_from_projection!(site, Val(:left), ψ, AC2, alg; normalize)
+end
+
+"""
+    expand_from_projection!(site, direction, ψ, AC2, alg::OptimalExpand; normalize = true)
+
+Expand the bond of `site` in `direction` of `ψ` from `AC2`, the projected two-site update of that
+bond, as [`changebond!`](@ref) does after forming it.
+"""
+function expand_from_projection!(
+        site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, AC2, alg::OptimalExpand; normalize::Bool = true
+    )
     left = ψ.AC[site]
     right = ψ.AR[site + 1]
     NL = left_null(left)
     NR = right_null!(_transpose_tail(right; copy = true))
 
-    # two-site update from the projected effective Hamiltonian
-    AC2 = AC2_projection(bond, ψ, H, ψ, envs; alg.backend, allocator)
-
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * AC2 * adjoint(NR)
     nrm = norm(g2)
@@ -155,21 +175,14 @@ function changebond!(
     ψ.AC[site + 1] = (nc, nar)
     return ψ
 end
-function changebond!(
-        site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::OptimalExpand, envs;
-        normalize::Bool = true,
-        # a sweep that already holds one should pass it: this is called once per site
-        allocator = default_allocator(ψ, SerialScheduler()),
+function expand_from_projection!(
+        site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, AC2, alg::OptimalExpand; normalize::Bool = true
     )
-    bond = site - 1
     left = ψ.AL[site - 1]
     right = ψ.AC[site]
     NL = left_null(left)
     NR = right_null!(_transpose_tail(right; copy = true))
 
-    # two-site update from the projected effective Hamiltonian
-    AC2 = AC2_projection(bond, ψ, H, ψ, envs; alg.backend, allocator)
-
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * AC2 * adjoint(NR)
     nrm = norm(g2)
diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl
index 51e680d04..9d760cce5 100644
--- a/src/algorithms/derivatives/hamiltonian_derivatives.jl
+++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl
@@ -11,6 +11,10 @@ const _HAM_MPS_TYPES = Union{
 
 Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
 In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
+
+`GL_I` and `GR_E` are the left environment in the not-started channel and the right
+environment in the finished channel. They are identities, stored as `missing`, if both MPSs
+are the same.
 """
 struct JordanMPO_AC_Hamiltonian{O1, O2, O3, Bk <: AbstractBackend, Al} <: DerivativeOperator
     D::Union{O1, Missing} # onsite
@@ -19,24 +23,28 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3, Bk <: AbstractBackend, Al} <: Deriva
     C::Union{O2, Missing} # starting
     B::Union{O2, Missing} # ending
     A::Union{O3, Missing} # continuing
+    GL_I::Union{O1, Missing} # not started, left environment
+    GR_E::Union{O1, Missing} # finished, right environment
     backend::Bk           # contraction backend used by the matvec
     allocator::Al         # scratch-buffer allocator used by the matvec
 
     function JordanMPO_AC_Hamiltonian{O1, O2, O3, Bk, Al}(
             D::Union{O1, Missing}, I::Union{O1, Missing}, E::Union{O1, Missing},
             C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing},
+            GL_I::Union{O1, Missing}, GR_E::Union{O1, Missing},
             backend::Bk, allocator::Al
         ) where {O1, O2, O3, Bk <: AbstractBackend, Al}
-        return new{O1, O2, O3, Bk, Al}(D, I, E, C, B, A, backend, allocator)
+        return new{O1, O2, O3, Bk, Al}(D, I, E, C, B, A, GL_I, GR_E, backend, allocator)
     end
 end
 function JordanMPO_AC_Hamiltonian{O1, O2, O3}(
-        D, I, E, C, B, A, backend = DefaultBackend(), allocator = DefaultAllocator()
+        D, I, E, C, B, A, GL_I, GR_E, backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {O1, O2, O3}
     return JordanMPO_AC_Hamiltonian{O1, O2, O3, typeof(backend), typeof(allocator)}(
         ismissing(D) ? D : convert(O1, D), ismissing(I) ? I : convert(O1, I),
         ismissing(E) ? E : convert(O1, E), ismissing(C) ? C : convert(O2, C),
         ismissing(B) ? B : convert(O2, B), ismissing(A) ? A : convert(O3, A),
+        ismissing(GL_I) ? GL_I : convert(O1, GL_I), ismissing(GR_E) ? GR_E : convert(O1, GR_E),
         backend, allocator
     )
 end
@@ -46,114 +54,40 @@ function AC_hamiltonian(
         prepare::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
-    @assert below === above "JordanMPO assumptions break"
+    prepare && has_halves(envs) &&
+        return hamiltonian_from_halves(Val(1), site, below, envs, backend, allocator)
     GL = leftenv(envs, site, below; backend, allocator)
     GR = rightenv(envs, site, below; backend, allocator)
     W = operator[site]
-    H_AC = JordanMPO_AC_Hamiltonian(GL, W, GR; backend, allocator)
+    H_AC = JordanMPO_AC_Hamiltonian(GL, W, GR; samestate = below === above, fold = prepare, backend, allocator)
     return prepare ? prepare_operator!!(H_AC) : H_AC
 end
 
+"""
+    JordanMPO_AC_Hamiltonian(GL, W, GR; samestate = true, fold = true, backend, allocator)
+
+Construct the single-site derivative from its environments. If `samestate`, the bra and ket
+are the same state, whose not-started left and finished right environments are identities.
+If `fold`, the smaller terms are folded into the larger ones, which only pays off for an operator
+that is applied repeatedly.
+"""
 function JordanMPO_AC_Hamiltonian(
-        GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor;
+        GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor; samestate::Bool = true, fold::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
-    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
-    # them once and reuse the locals throughout
-    WA, WB, WC, WD = W.A, W.B, W.C, W.D
-    GL2 = GL[2:(end - 1)]
-    GR2 = GR[2:(end - 1)]
-
-    # onsite
-    D = nonzero_length(WD) > 0 ? only(WD) : missing
-
-    # not started
-    I = size(W, 4) == 1 ? missing : removeunit(GR[1], 2)
-
-    # finished
-    E = size(W, 1) == 1 ? missing : removeunit(GL[end], 2)
-
-    # starting
-    C = if nonzero_length(WC) > 0
-        @plansor backend = backend allocator = allocator starting[-1 -2; -3 -4] ≔ WC[-1; -3 1] * GR2[-4 1; -2]
-        only(starting)
-    else
-        missing
-    end
-
-    # ending
-    B = if nonzero_length(WB) > 0
-        @plansor backend = backend allocator = allocator ending[-1 -2; -3 -4] ≔ GL2[-1 1; -3] * WB[1 -2; -4]
-        only(ending)
-    else
-        missing
-    end
-
-    # continuing
-    A = MPO_AC_Hamiltonian(GL2, WA, GR2, backend, allocator)
-
-    # obtaining storagetype of environments since these should have already mixed
-    # the types of the operator and state
-    S = spacetype(GL)
-    M = storagetype(GL)
-    O1 = tensormaptype(S, 1, 1, M)
-    O2 = tensormaptype(S, 2, 2, M)
-    O3 = typeof(A)
-
+    WA = W.A
+    A = MPO_AC_Hamiltonian(GL[2:(end - 1)], WA, GR[2:(end - 1)], backend, allocator)
+    L = jordan_hamiltonian_left(Val(1), GL, W, samestate, fold, backend, allocator)
+    R = jordan_hamiltonian_right(Val(1), GR, W, samestate, fold, backend, allocator)
     # specialization for nearest neighbours
-    nonzero_length(WA) == 0 && (A = missing)
-
-    return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A, backend, allocator)
+    return assemble_operator_halves(L, R, nonzero_length(WA) == 0 ? missing : A, typeof(A), backend, allocator)
 end
 
-function prepare_operator!!(
-        H::JordanMPO_AC_Hamiltonian{O1, O2, O3}
-    ) where {O1, O2, O3}
-    backend, allocator = H.backend, H.allocator
-    C::Union{Missing, O2} = H.C
-    B::Union{Missing, O2} = H.B
-
-    # onsite
-    D::Union{Missing, O1} = if ismissing(H.D)
-        missing
-    elseif !ismissing(C)
-        Id = TensorKit.id(storagetype(C), space(C, 2))
-        @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += H.D[-1; -3] * Id[-2; -4]
-        missing
-    elseif !ismissing(B)
-        Id = TensorKit.id(storagetype(B), space(B, 1))
-        @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += Id[-1; -3] * H.D[-2; -4]
-        missing
-    else
-        H.D
-    end
-
-    # not_started
-    I::Union{Missing, O1} = if ismissing(H.I)
-        missing
-    elseif !ismissing(C)
-        Id = id(storagetype(C), space(C, 1))
-        @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += Id[-1; -3] * H.I[-4; -2]
-        missing
-    else
-        H.I
-    end
-
-    # finished
-    E::Union{Missing, O1} = if ismissing(H.E)
-        missing
-    elseif !ismissing(B)
-        Id = id(storagetype(B), space(B, 2))
-        @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += H.E[-1; -3] * Id[-2; -4]
-        missing
-    else
-        H.E
-    end
-
-    O3′ = prepared_operator_type(O3)
-    A = ismissing(H.A) ? H.A : prepare_operator!!(H.A)
-
-    return JordanMPO_AC_Hamiltonian{O1, O2, O3′}(D, I, E, C, B, A, backend, allocator)::JordanMPO_AC_Hamiltonian{O1, O2, O3′}
+function prepare_operator!!(H::JordanMPO_AC_Hamiltonian{O1, O2, O3}) where {O1, O2, O3}
+    A = ismissing(H.A) ? missing : prepare_operator!!(H.A)
+    return JordanMPO_AC_Hamiltonian{O1, O2, prepared_operator_type(O3)}(
+        H.D, H.I, H.E, H.C, H.B, A, H.GL_I, H.GR_E, H.backend, H.allocator
+    )
 end
 
 
@@ -162,8 +96,10 @@ end
 """
     JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}
 
-Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
+Efficient operator for representing the two-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
 In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
+
+`GL_I` and `GR_E` are as for [`JordanMPO_AC_Hamiltonian`](@ref).
 """
 struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, Bk <: AbstractBackend, Al} <: DerivativeOperator
     II::Union{O1, Missing} # not_started
@@ -176,6 +112,8 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, Bk <: AbstractBackend, Al} <: D
     BE::Union{O2, Missing} # ending left
     DE::Union{O1, Missing} # onsite left
     EE::Union{O1, Missing} # finished
+    GL_I::Union{O1, Missing} # not started, left environment
+    GR_E::Union{O1, Missing} # finished, right environment
     backend::Bk            # contraction backend used by the matvec
     allocator::Al          # scratch-buffer allocator used by the matvec
 
@@ -184,13 +122,14 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, Bk <: AbstractBackend, Al} <: D
             CB::Union{O2, Missing}, CA::Union{O3, Missing},
             AB::Union{O3, Missing}, AA::Union{O4, Missing},
             BE::Union{O2, Missing}, DE::Union{O1, Missing}, EE::Union{O1, Missing},
+            GL_I::Union{O1, Missing}, GR_E::Union{O1, Missing},
             backend::Bk, allocator::Al
         ) where {O1, O2, O3, O4, Bk <: AbstractBackend, Al}
-        return new{O1, O2, O3, O4, Bk, Al}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE, backend, allocator)
+        return new{O1, O2, O3, O4, Bk, Al}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE, GL_I, GR_E, backend, allocator)
     end
 end
 function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(
-        II, IC, ID, CB, CA, AB, AA, BE, DE, EE,
+        II, IC, ID, CB, CA, AB, AA, BE, DE, EE, GL_I, GR_E,
         backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {O1, O2, O3, O4}
     return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, typeof(backend), typeof(allocator)}(
@@ -199,6 +138,7 @@ function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(
         ismissing(CA) ? CA : convert(O3, CA), ismissing(AB) ? AB : convert(O3, AB),
         ismissing(AA) ? AA : convert(O4, AA), ismissing(BE) ? BE : convert(O2, BE),
         ismissing(DE) ? DE : convert(O1, DE), ismissing(EE) ? EE : convert(O1, EE),
+        ismissing(GL_I) ? GL_I : convert(O1, GL_I), ismissing(GR_E) ? GR_E : convert(O1, GR_E),
         backend, allocator
     )
 end
@@ -208,11 +148,12 @@ function AC2_hamiltonian(
         prepare::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
-    @assert below === above "JordanMPO assumptions break"
+    prepare && has_halves(envs) &&
+        return hamiltonian_from_halves(Val(2), site, below, envs, backend, allocator)
     GL = leftenv(envs, site, below; backend, allocator)
     GR = rightenv(envs, site + 1, below; backend, allocator)
     W1, W2 = operator[site], operator[site + 1]
-    H_AC2 = JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR; backend, allocator)
+    H_AC2 = JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR; samestate = below === above, fold = prepare, backend, allocator)
     return prepare ? prepare_operator!!(H_AC2) : H_AC2
 end
 
@@ -226,14 +167,14 @@ for f in (:AC_hamiltonian, :AC2_hamiltonian)
 end
 
 """
-    _connected_channels(A1, A2) -> Union{Nothing, NTuple{3, Vector{Int}}}
+    _connected_channels(A1, A2) -> NTuple{3, Vector{Int}}
 
 The continuing-channel indices that can carry a contribution across both bonds, as `(rows, mids, cols)`.
 A middle channel counts only if it is reachable from some row of `A1` *and* reaches some column of `A2`,
 and a row or column counts only if it meets such a middle channel.
 Everything else contributes exactly zero.
 
-Returns `nothing` when nothing connects, which is the nearest-neighbour case.
+Returns three empty index vectors when nothing connects, as in the nearest-neighbour case.
 
 The channel indices are small dense integers, so reachability is tracked in flat bit-flag arrays and read
 out with `findall`, which is already ascending.
@@ -241,7 +182,7 @@ That ordering is load-bearing: the three index vectors are used to slice `A1`/`A
 which only line up if all of them keep the original channel order.
 """
 function _connected_channels(A1, A2)
-    (nonzero_length(A1) == 0 || nonzero_length(A2) == 0) && return nothing
+    (nonzero_length(A1) == 0 || nonzero_length(A2) == 0) && return (Int[], Int[], Int[])
     keys1, keys2 = nonzero_keys(A1), nonzero_keys(A2)
     @assert size(A1, 4) == size(A2, 1) "A-blocks do not share a bond"
 
@@ -255,7 +196,7 @@ function _connected_channels(A1, A2)
         @inbounds flags[I[1]] |= FEEDS
     end
     mids = findall(==(BOTH), flags)
-    isempty(mids) && return nothing
+    isempty(mids) && return (Int[], Int[], Int[])
 
     rows = falses(size(A1, 1))
     for I in keys1
@@ -270,231 +211,407 @@ function _connected_channels(A1, A2)
 end
 
 function JordanMPO_AC2_Hamiltonian(
-        GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor;
+        GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor; samestate::Bool = true, fold::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
-    # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind
-    # them once and reuse the locals throughout
-    A1, B1, C1, D1 = W1.A, W1.B, W1.C, W1.D
-    A2, B2, C2, D2 = W2.A, W2.B, W2.C, W2.D
-    GL2 = GL[2:(end - 1)]
-    GR2 = GR[2:(end - 1)]
-
-    # not started
-    II = size(W2, 4) == 1 ? missing : transpose(removeunit(GR[1], 2))
+    pair = JordanPairData(W1, W2, storagetype(GL))
+    # continuing - continuing, restricted to the channels that can actually contribute
+    rows, mids, cols = pair.channels
+    GL2, GR2 = GL[2:(end - 1)], GR[2:(end - 1)]
+    O4 = typeof(MPO_AC2_Hamiltonian(GL2, pair.A1, pair.A2, GR2, backend, allocator))
+    AA = isempty(mids) ? missing : MPO_AC2_Hamiltonian(GL2[rows], pair.A1, pair.A2, GR2[cols], backend, allocator)
+    L = jordan_hamiltonian_left(Val(2), GL, pair, samestate, fold, backend, allocator)
+    R = jordan_hamiltonian_right(Val(2), GR, pair, samestate, fold, backend, allocator)
+    return assemble_operator_halves(L, R, AA, O4, backend, allocator)
+end
 
-    # finished
-    EE = size(W1, 1) == 1 ? missing : removeunit(GL[end], 2)
+function prepare_operator!!(H::JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}) where {O1, O2, O3, O4}
+    AA = ismissing(H.AA) ? missing : prepare_operator!!(H.AA)
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, prepared_operator_type(O4)}(
+        H.II, H.IC, H.ID, H.CB, H.CA, H.AB, AA, H.BE, H.DE, H.EE, H.GL_I, H.GR_E, H.backend, H.allocator
+    )
+end
 
-    # starting right
-    IC = if nonzero_length(C2) > 0
-        @plansor backend = backend allocator = allocator IC_[-1 -2; -3 -4] ≔ C2[-1; -3 1] * GR2[-4 1; -2]
-        only(IC_)
+# Actions
+# -------
+# Terms that start before the window act on `xL`, the ket with the not-started left
+# environment applied, and terms that finish within it on `xR`, with the finished right one.
+function _apply_boundaries(GL_I, x::MPSTensor, GR_E, backend, allocator)
+    ismissing(GL_I) && ismissing(GR_E) && return x, x
+    xL = if ismissing(GL_I)
+        x
     else
-        missing
+        @plansor backend = backend allocator = allocator xL[-1 -2; -3] := GL_I[-1; 1] * x[1 -2; -3]
     end
-
-    # onsite left
-    DE = nonzero_length(D1) > 0 ? only(D1) : missing
-
-    # onsite right
-    ID = nonzero_length(D2) > 0 ? only(D2) : missing
-
-    # starting left - ending right
-    CB = if nonzero_length(C1) > 0 && nonzero_length(B2) > 0
-        @plansor backend = backend allocator = allocator CB_[-1 -2; -3 -4] ≔ C1[-1; -3 1] * B2[1 -2; -4]
-        # have to convert to complex if hamiltonian is real but states are complex
-        scalartype(GL) <: Complex ? complex(only(CB_)) : only(CB_)
+    xR = if ismissing(GR_E)
+        x
     else
-        missing
+        @plansor backend = backend allocator = allocator xR[-1 -2; -3] := x[-1 -2; 1] * GR_E[-3; 1]
     end
-
-    # starting left - continuing right
-    CA = if nonzero_length(C1) > 0 && nonzero_length(A2) > 0
-        @plansor backend = backend allocator = allocator CA_[-1 -2 -3; -4 -5 -6] ≔ C1[-1; -4 1] * A2[1 -2; -5 2] *
-            GR2[-6 2; -3]
-        only(CA_)
+    return xL, xR
+end
+function _apply_boundaries(GL_I, x::MPOTensor, GR_E, backend, allocator)
+    ismissing(GL_I) && ismissing(GR_E) && return x, x
+    xL = if ismissing(GL_I)
+        x
     else
-        missing
+        @plansor backend = backend allocator = allocator xL[-1 -2; -3 -4] := GL_I[-1; 1] * x[1 -2; -3 -4]
     end
-
-    # continuing left - ending right
-    AB = if nonzero_length(A1) > 0 && nonzero_length(B2) > 0
-        @plansor backend = backend allocator = allocator AB_[-1 -2 -3; -4 -5 -6] ≔ GL2[-1 2; -4] * A1[2 -2; -5 1] *
-            B2[1 -3; -6]
-        only(AB_)
+    xR = if ismissing(GR_E)
+        x
     else
-        missing
+        @plansor backend = backend allocator = allocator xR[-1 -2; -3 -4] := x[-1 -2; 1 -4] * GR_E[-3; 1]
     end
+    return xL, xR
+end
+# The output lives on the bra, whose virtual spaces are those of `xL` on the left and `xR` on the right.
+_zero_output(x, xL, xR) = xL === xR ? zerovector(x) : zerovector!(similar(x, codomain(xL) ← domain(xR)))
+# The identity in the onsite folds of `prepare_operator!!` is the boundary environment.
+_boundary(G, t, i) = ismissing(G) ? TensorKit.id(storagetype(t), space(t, i)) : G
 
-    # ending left
-    BE = if nonzero_length(B1) > 0
-        @plansor backend = backend allocator = allocator BE_[-1 -2; -3 -4] ≔ GL2[-1 2; -3] * B1[2 -2; -4]
-        only(BE_)
-    else
-        missing
+function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor)
+    backend, allocator = H.backend, H.allocator
+    xL, xR = _apply_boundaries(H.GL_I, x, H.GR_E, backend, allocator)
+    y = ismissing(H.A) ? _zero_output(x, xL, xR) : H.A(x)
+
+    if !ismissing(H.D)
+        xLR = last(_apply_boundaries(missing, xL, H.GR_E, backend, allocator))
+        @plansor backend = backend allocator = allocator y[-1 -2; -3] += xLR[-1 1; -3] * H.D[-2; 1]
     end
+    ismissing(H.E) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.E[-1; 1] * xR[1 -2; -3]
+    ismissing(H.I) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += xL[-1 -2; 1] * H.I[1; -3]
+    ismissing(H.C) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += xL[-1 2; 1] * H.C[-2 -3; 2 1]
+    ismissing(H.B) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.B[-1 -2; 1 2] * xR[1 2; -3]
 
-    S = spacetype(GL)
-    M = storagetype(GL)
-    O1 = tensormaptype(S, 1, 1, M)
-    O2 = tensormaptype(S, 2, 2, M)
-    O3 = tensormaptype(S, 3, 3, M)
-    # slicing preserves the block-tensor types, so `AA`'s type does not depend on whether the
-    # channels end up restricted - no need to build a throwaway operator just to read it off
-    O4 = MPO_AC2_Hamiltonian{
-        typeof(GL2), typeof(A1), typeof(A2), typeof(GR2), typeof(backend), typeof(allocator),
-    }
+    return y
+end
 
-    # continuing - continuing, restricted to the channels that can actually contribute
-    channels = _connected_channels(A1, A2)
-    AA = if isnothing(channels)
-        missing
+function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor)
+    backend, allocator = H.backend, H.allocator
+    xL, xR = _apply_boundaries(H.GL_I, x, H.GR_E, backend, allocator)
+    y = ismissing(H.AA) ? _zero_output(x, xL, xR) : H.AA(x)
+    xLR = if ismissing(H.ID) && ismissing(H.CB) && ismissing(H.DE)
+        x
     else
-        rows, mids, cols = channels
-        MPO_AC2_Hamiltonian(
-            GL2[rows], A1[rows, 1:1, 1:1, mids], A2[mids, 1:1, 1:1, cols], GR2[cols],
-            backend, allocator
-        )
+        last(_apply_boundaries(missing, xL, H.GR_E, backend, allocator))
     end
 
-    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(
-        II, IC, ID,
-        CB, CA,
-        AB, AA,
-        BE, DE, EE,
-        backend, allocator
-    )
+    ismissing(H.II) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xL[-1 -2; 1 -4] * H.II[-3; 1]
+    ismissing(H.IC) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xL[-1 -2; 1 2] * H.IC[-4 -3; 2 1]
+    ismissing(H.ID) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xLR[-1 -2; -3 1] * H.ID[-4; 1]
+    ismissing(H.CB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xLR[-1 1; -3 2] * H.CB[-2 -4; 1 2]
+    ismissing(H.CA) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xL[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3]
+    ismissing(H.AB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xR[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3]
+    ismissing(H.BE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xR[1 2; -3 -4] * H.BE[-1 -2; 1 2]
+    ismissing(H.DE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xLR[-1 1; -3 -4] * H.DE[-2; 1]
+    ismissing(H.EE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += xR[1 -2; -3 -4] * H.EE[-1; 1]
 
+    return y
 end
 
-function prepare_operator!!(
-        H::JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}
-    ) where {O1, O2, O3, O4}
-    backend, allocator = H.backend, H.allocator
 
-    CA::Union{Missing, O3} = H.CA
-    AB::Union{Missing, O3} = H.AB
+# Left and right halves
+# ---------------------
+"""
+    JordanPairData
 
-    CB::Union{Missing, O2} = if !ismissing(CA) && !ismissing(H.CB)
-        Id = TensorKit.id(storagetype(H.CB), space(CA, 3))
-        @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += H.CB[-1 -2; -4 -5] * Id[-3; -6]
-        missing
-    elseif !ismissing(AB) && !ismissing(H.CB)
-        Id = TensorKit.id(storagetype(H.CB), space(AB, 1))
-        @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.CB[-2 -3; -5 -6] * Id[-1; -4]
-        missing
-    else
-        H.CB
+Operator-only data for the adjacent sites `i` and `i + 1` with operators `W1` and `W2`: the
+products `C_i * B_{i+1}`, `A_i * B_{i+1}` and `C_i * A_{i+1}`, each `missing` if absent, and
+the continuing blocks restricted to the connected `channels = (rows, mids, cols)` of
+`_connected_channels`.
+"""
+struct JordanPairData{W <: JordanMPOTensor, T <: MPOTensor, U, V, A <: MPOTensor}
+    W1::W
+    W2::W
+    CB::Union{Missing, T}
+    AB::Union{Missing, U}
+    CA::Union{Missing, V}
+    channels::NTuple{3, Vector{Int}}
+    A1::A
+    A2::A
+    function JordanPairData{W, T, U, V, A}(
+            W1, W2, CB, AB, CA, channels, A1, A2
+        ) where {W <: JordanMPOTensor, T <: MPOTensor, U, V, A <: MPOTensor}
+        return new{W, T, U, V, A}(W1, W2, CB, AB, CA, channels, A1, A2)
     end
+end
 
-    # starting right
-    IC::Union{Missing, O2} = if !ismissing(CA) && !ismissing(H.IC)
-        Id = TensorKit.id(storagetype(H.IC), space(CA, 1))
-        @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * H.IC[ -2 -3; -5 -6]
-        missing
+Base.@assume_effects :foldable function jordan_pair_type(::Type{W}, ::Type{M}) where {W <: JordanMPOTensor, M}
+    S = spacetype(W)
+    T = tensormaptype(S, 2, 2, M)
+    O32, O23 = tensormaptype(S, 3, 2, M), tensormaptype(S, 2, 3, M)
+    U = SparseBlockTensorMap{O32, scalartype(O32), S, 3, 2, 5}
+    V = SparseBlockTensorMap{O23, scalartype(O23), S, 2, 3, 5}
+    A = SparseBlockTensorMap{tensormaptype(S, 2, 2, storagetype(W)), scalartype(W), S, 2, 2, 4}
+    return JordanPairData{W, T, U, V, A}
+end
+
+function JordanPairData(W1::JordanMPOTensor, W2::JordanMPOTensor, ::Type{M}) where {M}
+    # block accessors recompute a fresh `SparseBlockTensorMap` on every access
+    A1, B2, C1, A2 = W1.A, W2.B, W1.C, W2.A
+    CB = if nonzero_length(C1) > 0 && nonzero_length(B2) > 0
+        @plansor cb[-1 -2; -3 -4] := C1[-1; -3 1] * B2[1 -2; -4]
+        only(cb)
     else
-        H.IC
+        missing
     end
-
-    # ending left
-    BE::Union{Missing, O2} = if !ismissing(AB) && !ismissing(H.BE)
-        Id = TensorKit.id(storagetype(H.BE), space(AB, 3))
-        @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.BE[-1 -2; -4 -5] * Id[-3; -6]
+    AB = if nonzero_length(A1) > 0 && nonzero_length(B2) > 0
+        @plansor ab[-1 -2 -3; -4 -5] := A1[-1 -2; -4 1] * B2[1 -3; -5]
+    else
         missing
+    end
+    CA = if nonzero_length(C1) > 0 && nonzero_length(A2) > 0
+        @plansor ca[-1 -2; -3 -4 -5] := C1[-1; -3 1] * A2[1 -2; -4 -5]
     else
-        H.BE
+        missing
     end
+    channels = _connected_channels(A1, A2)
+    rows, mids, cols = channels
+    A1r, A2r = A1[rows, 1:1, 1:1, mids], A2[mids, 1:1, 1:1, cols]
+    return jordan_pair_type(typeof(W1), M)(W1, W2, CB, AB, CA, channels, A1r, A2r)
+end
 
-    # onsite left
-    DE::Union{Missing, O1} = if !ismissing(BE) && !ismissing(H.DE)
-        Id = TensorKit.id(storagetype(H.DE), space(BE, 1))
-        @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += Id[-1; -3] * H.DE[-2; -4]
-        missing
-    elseif !ismissing(AB) && !ismissing(H.DE)
-        Id1 = id(storagetype(H.DE), space(AB, 1))
-        Id2 = id(storagetype(H.DE), space(AB, 3))
-        @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.DE[-2; -5] * Id2[-3; -6]
-        missing
-        # TODO: could also try in CA?
-    else
-        H.DE
+_jordan_types(G) = _jordan_types(typeof(G))
+function _jordan_types(::Type{G}) where {G <: AbstractTensorMap}
+    S, M = spacetype(G), storagetype(G)
+    return tensormaptype(S, 1, 1, M), tensormaptype(S, 2, 2, M), tensormaptype(S, 3, 3, M)
+end
+
+_onsite(W) = nonzero_length(W.D) > 0 ? only(W.D) : missing
+function _ending(GL, W, backend, allocator)
+    WB = W.B
+    nonzero_length(WB) > 0 || return missing
+    GL2 = GL[2:(end - 1)]
+    @plansor backend = backend allocator = allocator b[-1 -2; -3 -4] := GL2[-1 1; -3] * WB[1 -2; -4]
+    return only(b)
+end
+function _starting(GR, W, backend, allocator)
+    WC = W.C
+    nonzero_length(WC) > 0 || return missing
+    GR2 = GR[2:(end - 1)]
+    @plansor backend = backend allocator = allocator c[-1 -2; -3 -4] := WC[-1; -3 1] * GR2[-4 1; -2]
+    return only(c)
+end
+
+# The terms of `JordanMPO_AC_Hamiltonian` and `JordanMPO_AC2_Hamiltonian` that depend on only
+# one environment, excluding the continuing `A`/`AA` term, with the default backend and
+# allocator. Each onsite, not-started, finished or operator-only term is assigned to one half,
+# and, if `fold`, folded into the largest term of that half that applies to the same boundaries.
+function jordan_hamiltonian_left(::Val{1}, GL, W::JordanMPOTensor, samestate::Bool, fold::Bool, backend, allocator)
+    O1, O2, _ = _jordan_types(GL)
+    GL_I = samestate ? missing : removeunit(GL[1], 2)
+    E = size(W, 1) > 1 ? removeunit(GL[end], 2) : missing
+    # without folding, or without `C` to fold it into, the onsite term stays in this half
+    D = fold && nonzero_length(W.C) > 0 ? missing : _onsite(W)
+    B = _ending(GL, W, backend, allocator)
+    if fold && !ismissing(B)
+        if !ismissing(D)
+            GL_I′ = _boundary(GL_I, B, 1)
+            @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += GL_I′[-1; -3] * D[-2; -4]
+        end
+        if !ismissing(E)
+            Id = id(storagetype(B), space(B, 2))
+            @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += E[-1; -3] * Id[-2; -4]
+        end
+        D = E = missing
+    end
+    return JordanMPO_AC_Hamiltonian{O1, O2, Missing}(D, missing, E, missing, B, missing, GL_I, missing)
+end
+function jordan_hamiltonian_right(::Val{1}, GR, W::JordanMPOTensor, samestate::Bool, fold::Bool, backend, allocator)
+    O1, O2, _ = _jordan_types(GR)
+    GR_E = samestate ? missing : transpose(removeunit(GR[end], 2))
+    I = size(W, 4) > 1 ? removeunit(GR[1], 2) : missing
+    C = _starting(GR, W, backend, allocator)
+    if fold && !ismissing(C)
+        D = _onsite(W)
+        if !ismissing(D)
+            GR_E′ = _boundary(GR_E, C, 2)
+            @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += D[-1; -3] * GR_E′[-2; -4]
+        end
+        if !ismissing(I)
+            Id = id(storagetype(C), space(C, 1))
+            @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += Id[-1; -3] * I[-4; -2]
+        end
+        I = missing
     end
+    return JordanMPO_AC_Hamiltonian{O1, O2, Missing}(missing, I, missing, C, missing, missing, missing, GR_E)
+end
 
-    # onsite right
-    ID::Union{Missing, O1} = if !ismissing(IC) && !ismissing(H.ID)
-        Id = TensorKit.id(storagetype(H.ID), space(IC, 2))
-        @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += H.ID[-1; -3] * Id[-2; -4]
-        missing
-    elseif !ismissing(CA) && !ismissing(H.ID)
-        Id1 = TensorKit.id(storagetype(H.ID), space(CA, 1))
-        Id2 = TensorKit.id(storagetype(H.ID), space(CA, 3))
-        @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.ID[-2; -5] * Id2[-3; -6]
+function jordan_hamiltonian_left(::Val{2}, GL, pair::JordanPairData, samestate::Bool, fold::Bool, backend, allocator)
+    O1, O2, O3 = _jordan_types(GL)
+    W = pair.W1
+    GL_I = samestate ? missing : removeunit(GL[1], 2)
+    EE = size(W, 1) > 1 ? removeunit(GL[end], 2) : missing
+    DE = _onsite(W)
+    BE = _ending(GL, W, backend, allocator)
+    # `CB` belongs to the right half if it has `CA` to fold it into.
+    CB = fold && !ismissing(pair.CA) ? missing : pair.CB
+    AB = if ismissing(pair.AB)
         missing
     else
-        H.ID
+        GL2 = GL[2:(end - 1)]
+        @plansor backend = backend allocator = allocator ab[-1 -2 -3; -4 -5 -6] := GL2[-1 1; -4] * pair.AB[1 -2 -3; -5 -6]
+        only(ab)
     end
-
-    # finished
-    II::Union{Missing, O1} = if !ismissing(IC) && !ismissing(H.II)
-        I = id(storagetype(H.II), space(IC, 1))
-        @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += I[-1; -3] * H.II[-2; -4]
-        II = missing
-    elseif !ismissing(CA) && !ismissing(H.II)
-        I = id(storagetype(H.II), space(CA, 1) ⊗ space(CA, 2))
-        @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * H.II[-3; -6]
-        II = missing
-    else
-        H.II
+    if fold && !ismissing(AB)
+        GL_I′ = _boundary(GL_I, AB, 1)
+        Id = id(storagetype(AB), space(AB, 3))
+        ismissing(CB) || @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += GL_I′[-1; -4] * CB[-2 -3; -5 -6]
+        ismissing(BE) || @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += BE[-1 -2; -4 -5] * Id[-3; -6]
+        ismissing(DE) || @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += GL_I′[-1; -4] * DE[-2; -5] * Id[-3; -6]
+        if !ismissing(EE)
+            Id2 = id(storagetype(AB), space(AB, 2) ⊗ space(AB, 3))
+            @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += EE[-1; -4] * Id2[-2 -3; -5 -6]
+        end
+        CB = BE = DE = EE = missing
+    elseif fold && !ismissing(BE)
+        if !ismissing(DE)
+            GL_I′ = _boundary(GL_I, BE, 1)
+            @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += GL_I′[-1; -3] * DE[-2; -4]
+        end
+        if !ismissing(EE)
+            Id = id(storagetype(BE), space(BE, 2))
+            @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += EE[-1; -3] * Id[-2; -4]
+        end
+        DE = EE = missing
     end
-
-    # unstarted
-    EE::Union{Missing, O1} = if !ismissing(BE) && !ismissing(H.EE)
-        I = id(storagetype(H.EE), space(BE, 2))
-        @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += H.EE[-1; -3] * I[-2; -4]
-        EE = missing
-    elseif !ismissing(AB) && !ismissing(H.EE)
-        I = id(storagetype(H.EE), space(AB, 2) ⊗ space(AB, 3))
-        @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.EE[-1; -4] * I[-2 -3; -5 -6]
-        EE = missing
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, Missing}(
+        missing, missing, missing, CB, missing, AB, missing, BE, DE, EE, GL_I, missing
+    )
+end
+function jordan_hamiltonian_right(::Val{2}, GR, pair::JordanPairData, samestate::Bool, fold::Bool, backend, allocator)
+    O1, O2, O3 = _jordan_types(GR)
+    W = pair.W2
+    GR_E = samestate ? missing : transpose(removeunit(GR[end], 2))
+    II = size(W, 4) > 1 ? transpose(removeunit(GR[1], 2)) : missing
+    ID = _onsite(W)
+    IC = _starting(GR, W, backend, allocator)
+    CA = if ismissing(pair.CA)
+        missing
     else
-        H.EE
+        GR2 = GR[2:(end - 1)]
+        @plansor backend = backend allocator = allocator ca[-1 -2 -3; -4 -5 -6] := pair.CA[-1 -2; -4 -5 1] * GR2[-6 1; -3]
+        only(ca)
+    end
+    if fold && !ismissing(CA)
+        GR_E′ = _boundary(GR_E, CA, 3)
+        Id = id(storagetype(CA), space(CA, 1))
+        CB = pair.CB
+        ismissing(CB) || @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += CB[-1 -2; -4 -5] * GR_E′[-3; -6]
+        ismissing(IC) || @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * IC[-2 -3; -5 -6]
+        ismissing(ID) || @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * ID[-2; -5] * GR_E′[-3; -6]
+        if !ismissing(II)
+            Id2 = id(storagetype(CA), space(CA, 1) ⊗ space(CA, 2))
+            @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id2[-1 -2; -4 -5] * II[-3; -6]
+        end
+        IC = ID = II = missing
+    elseif fold && !ismissing(IC)
+        if !ismissing(ID)
+            GR_E′ = _boundary(GR_E, IC, 2)
+            @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += ID[-1; -3] * GR_E′[-2; -4]
+        end
+        if !ismissing(II)
+            Id = id(storagetype(IC), space(IC, 1))
+            @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += Id[-1; -3] * II[-2; -4]
+        end
+        ID = II = missing
     end
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, Missing}(
+        II, IC, ID, missing, CA, missing, missing, missing, missing, missing, missing, GR_E
+    )
+end
 
-    O4′ = prepared_operator_type(O4)
-    AA = prepare_operator!!(H.AA)
+function assemble_operator_halves(L::JordanMPO_AC_Hamiltonian{O1, O2}, R, A, ::Type{O3}, backend, allocator) where {O1, O2, O3}
+    return JordanMPO_AC_Hamiltonian{O1, O2, O3}(L.D, R.I, L.E, R.C, L.B, A, L.GL_I, R.GR_E, backend, allocator)
+end
+function assemble_operator_halves(L::JordanMPO_AC2_Hamiltonian{O1, O2, O3}, R, AA, ::Type{O4}, backend, allocator) where {O1, O2, O3, O4}
+    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(
+        R.II, R.IC, R.ID, L.CB, R.CA, L.AB, AA, L.BE, L.DE, L.EE, L.GL_I, R.GR_E, backend, allocator,
+    )
+end
 
-    return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4′}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE, backend, allocator)
+# Prepared halves for `FiniteEnvironments`
+# -----------------------------------------
+"""
+    JordanHalf
+
+One environment's prepared half of an effective `AC` or `AC2` Hamiltonian: the terms of
+`jordan_hamiltonian_left` or `jordan_hamiltonian_right`, and `continuing`, that side's
+prepared environment of the continuing term, or `missing` if absent.
+"""
+struct JordanHalf{H, C <: AbstractTensorMap}
+    terms::H
+    continuing::Union{Missing, C}
+    function JordanHalf{H, C}(terms, continuing) where {H, C <: AbstractTensorMap}
+        return new{H, C}(terms, continuing)
+    end
 end
 
-# Actions
-# -------
-function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor)
-    backend, allocator = H.backend, H.allocator
-    y = ismissing(H.A) ? zerovector(x) : H.A(x)
+function jordan_halftypes(::Val{1}, ::Type{G}, ::Type{W}) where {G, W}
+    O1, O2, _ = _jordan_types(G)
+    H = JordanMPO_AC_Hamiltonian{O1, O2, Missing, DefaultBackend, DefaultAllocator}
+    L, R = halftypes(Val(1), G, W)
+    return JordanHalf{H, L}, JordanHalf{H, R}
+end
+function jordan_halftypes(::Val{2}, ::Type{G}, ::Type{W}) where {G, W}
+    O1, O2, O3 = _jordan_types(G)
+    H = JordanMPO_AC2_Hamiltonian{O1, O2, O3, Missing, DefaultBackend, DefaultAllocator}
+    L, R = halftypes(Val(2), G, W)
+    return JordanHalf{H, L}, JordanHalf{H, R}
+end
 
-    ismissing(H.D) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1]
-    ismissing(H.E) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3]
-    ismissing(H.I) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3]
-    ismissing(H.C) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1]
-    ismissing(H.B) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3]
+function environment_halves(O::MPOHamiltonian, ::Type{G}, N) where {G}
+    W = eltype(O)
+    (G <: MPSTensor && isconcretetype(W)) || return nothing
+    AC = PreparedHalves{jordan_halftypes(Val(1), G, W)..., G}(N)
+    AC2 = PreparedHalves{jordan_halftypes(Val(2), G, W)..., G}(N)
+    pairs = fill!(Vector{Union{Nothing, jordan_pair_type(W, storagetype(G))}}(undef, max(N - 1, 0)), nothing)
+    return (; AC, AC2, pairs)
+end
 
-    return y
+function pair_data!(envs::FiniteEnvironments, i::Int)
+    pairs = envs.halves.pairs
+    pair = pairs[i]
+    isnothing(pair) || return pair
+    O = envs.operator
+    return pairs[i] = JordanPairData(O[i], O[i + 1], storagetype(eltype(envs.GLs)))
 end
 
-function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor)
-    backend, allocator = H.backend, H.allocator
-    y = ismissing(H.AA) ? zerovector(x) : H.AA(x)
-
-    ismissing(H.II) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1]
-    ismissing(H.IC) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1]
-    ismissing(H.ID) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1]
-    ismissing(H.CB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2]
-    ismissing(H.CA) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3]
-    ismissing(H.AB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3]
-    ismissing(H.BE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2]
-    ismissing(H.DE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1]
-    ismissing(H.EE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1]
+function prepare_operator_left(n::Val{1}, O::MPOHamiltonian, envs, GL, i, backend, allocator)
+    W = O[i]
+    WA = W.A
+    terms = jordan_hamiltonian_left(n, GL, W, isnothing(envs.above), true, backend, allocator)
+    continuing = nonzero_length(WA) > 0 ? prepare_operator_left(n, GL[2:(end - 1)], WA, backend, allocator) : missing
+    L, _ = jordan_halftypes(n, eltype(envs.GLs), eltype(O))
+    return L(terms, continuing)
+end
+function prepare_operator_right(n::Val{1}, O::MPOHamiltonian, envs, GR, i, backend, allocator)
+    W = O[i]
+    WA = W.A
+    terms = jordan_hamiltonian_right(n, GR, W, isnothing(envs.above), true, backend, allocator)
+    continuing = nonzero_length(WA) > 0 ? prepare_operator_right(n, GR[2:(end - 1)], WA, backend, allocator) : missing
+    _, R = jordan_halftypes(n, eltype(envs.GRs), eltype(O))
+    return R(terms, continuing)
+end
+function prepare_operator_left(n::Val{2}, O::MPOHamiltonian, envs, GL, i, backend, allocator)
+    pair = pair_data!(envs, i)
+    terms = jordan_hamiltonian_left(n, GL, pair, isnothing(envs.above), true, backend, allocator)
+    rows, mids, _ = pair.channels
+    continuing = isempty(mids) ? missing : prepare_operator_left(n, GL[2:(end - 1)][rows], pair.A1, backend, allocator)
+    L, _ = jordan_halftypes(n, eltype(envs.GLs), eltype(O))
+    return L(terms, continuing)
+end
+function prepare_operator_right(n::Val{2}, O::MPOHamiltonian, envs, GR, i, backend, allocator)
+    pair = pair_data!(envs, i - 1)
+    terms = jordan_hamiltonian_right(n, GR, pair, isnothing(envs.above), true, backend, allocator)
+    _, mids, cols = pair.channels
+    continuing = isempty(mids) ? missing : prepare_operator_right(n, GR[2:(end - 1)][cols], pair.A2, backend, allocator)
+    _, R = jordan_halftypes(n, eltype(envs.GRs), eltype(O))
+    return R(terms, continuing)
+end
 
-    return y
+function assemble_operator_halves(L::JordanHalf{<:Any, CL}, R::JordanHalf{<:Any, CR}, backend, allocator) where {CL, CR}
+    P = precomputed_type(CL, CR, typeof(backend), typeof(allocator))
+    A = ismissing(L.continuing) ? missing : P(L.continuing, R.continuing, backend, allocator)
+    return assemble_operator_halves(L.terms, R.terms, A, P, backend, allocator)
 end
diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 1e7bde691..784035a65 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -47,6 +47,8 @@ function AC_hamiltonian(
         site::Int, below, operator, above, envs; prepare::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
+    prepare && operator isa Union{FiniteMPO, InfiniteMPO} && has_halves(envs) &&
+        return hamiltonian_from_halves(Val(1), site, below, envs, backend, allocator)
     O = isnothing(operator) ? nothing : operator[site]
     H_AC = MPO_AC_Hamiltonian(
         leftenv(envs, site, below), O, rightenv(envs, site, below), backend, allocator
@@ -57,6 +59,8 @@ function AC2_hamiltonian(
         site::Int, below, operator, above, envs; prepare::Bool = true,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
+    prepare && operator isa Union{FiniteMPO, InfiniteMPO} && has_halves(envs) &&
+        return hamiltonian_from_halves(Val(2), site, below, envs, backend, allocator)
     O1, O2 = isnothing(operator) ? (nothing, nothing) : (operator[site], operator[site + 1])
     H_AC2 = MPO_AC2_Hamiltonian(
         leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below),
@@ -177,10 +181,13 @@ struct PrecomputedDerivative{
     allocator::A
 end
 function PrecomputedDerivative(L::AbstractTensorMap, R::AbstractTensorMap, backend, allocator)
+    return precomputed_type(typeof(L), typeof(R), typeof(backend), typeof(allocator))(L, R, backend, allocator)
+end
+Base.@assume_effects :foldable function precomputed_type(::Type{L}, ::Type{R}, ::Type{B}, ::Type{A}) where {L, R, B, A}
     S = TensorKit.check_spacetype(L, R)
     T = TensorOperations.promote_contract(scalartype(L), scalartype(R))
     M = TensorKit.promote_storagetype(T, L, R)
-    return PrecomputedDerivative{T, S, M, numout(L), numin(L), numout(R), numin(R), typeof(backend), typeof(allocator)}(L, R, backend, allocator)
+    return PrecomputedDerivative{T, S, M, numout(L), numin(L), numout(R), numin(R), B, A}
 end
 
 const PrecomputedCDerivative{T, S, M, B, A} = PrecomputedDerivative{T, S, M, 1, 2, 2, 1, B, A}
@@ -243,7 +250,10 @@ end
     return env
 end
 
-function _prepare_GL_O(GL, O, backend, allocator)
+# The prepared left and right environments of a `PrecomputedDerivative` of an `N`-site window:
+# a single-site operator is absorbed into the left one, and each operator of a two-site window
+# into its own side.
+function prepare_operator_left(::Val, GL, O, backend, allocator)
     cp = allocator_checkpoint!(allocator)
 
     TC = TensorOperations.promote_contract(scalartype(GL), scalartype(O))
@@ -264,7 +274,8 @@ function _prepare_GL_O(GL, O, backend, allocator)
     return leftenv
 end
 
-function _prepare_O_GR(O, GR, backend, allocator)
+prepare_operator_right(::Val{1}, GR, O, backend, allocator) = GR isa TensorMap ? GR : TensorMap(GR)
+function prepare_operator_right(::Val{2}, GR, O, backend, allocator)
     cp = allocator_checkpoint!(allocator)
 
     TC = TensorOperations.promote_contract(scalartype(O), scalartype(GR))
@@ -287,25 +298,19 @@ end
 
 function prepare_operator!!(H::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})
     backend, allocator = H.backend, H.allocator
-
-    GL, O = H.leftenv, H.operators[1]
-    leftenv = _prepare_GL_O(GL, O, backend, allocator)
-
-    rightenv = H.rightenv isa TensorMap ? H.rightenv : TensorMap(H.rightenv)
-
-    return prepared_operator_type(typeof(H))(leftenv, rightenv, backend, allocator)
+    O = only(H.operators)
+    L = prepare_operator_left(Val(1), H.leftenv, O, backend, allocator)
+    R = prepare_operator_right(Val(1), H.rightenv, O, backend, allocator)
+    return prepared_operator_type(typeof(H))(L, R, backend, allocator)
 end
 
 function prepare_operator!!(
         H::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor}
     )
     backend, allocator = H.backend, H.allocator
-
-    GL, O₁, O₂, GR = H.leftenv, H.operators[1], H.operators[2], H.rightenv
-    leftenv = _prepare_GL_O(GL, O₁, backend, allocator)
-    rightenv = _prepare_O_GR(O₂, GR, backend, allocator)
-
-    return prepared_operator_type(typeof(H))(leftenv, rightenv, backend, allocator)
+    L = prepare_operator_left(Val(2), H.leftenv, H.operators[1], backend, allocator)
+    R = prepare_operator_right(Val(2), H.rightenv, H.operators[2], backend, allocator)
+    return prepared_operator_type(typeof(H))(L, R, backend, allocator)
 end
 
 function (H::PrecomputedCDerivative)(x::MPSBondTensor)
@@ -389,3 +394,40 @@ const _ToPrepare = Union{
 function prepare_operator!!(H::Multiline{<:_ToPrepare})
     return Multiline(map(prepare_operator!!, parent(H)))
 end
+
+# Prepared halves for `FiniteEnvironments`
+# -----------------------------------------
+halfstorage(envs::FiniteEnvironments, ::Val{1}) = envs.halves.AC
+halfstorage(envs::FiniteEnvironments, ::Val{2}) = envs.halves.AC2
+
+# The types of the halves are those of the environments of the prepared operator.
+function halftypes(::Val{1}, ::Type{G}, ::Type{W}) where {G, W}
+    P = prepared_operator_type(MPO_AC_Hamiltonian{G, W, G, DefaultBackend, DefaultAllocator})
+    return fieldtype(P, :leftenv), fieldtype(P, :rightenv)
+end
+function halftypes(::Val{2}, ::Type{G}, ::Type{W}) where {G, W}
+    P = prepared_operator_type(MPO_AC2_Hamiltonian{G, W, W, G, DefaultBackend, DefaultAllocator})
+    return fieldtype(P, :leftenv), fieldtype(P, :rightenv)
+end
+
+function environment_halves(O::Union{FiniteMPO, InfiniteMPO}, ::Type{G}, N) where {G}
+    W = eltype(O)
+    (G <: MPSTensor && W <: MPOTensor && isconcretetype(W)) || return nothing
+    AC = PreparedHalves{halftypes(Val(1), G, W)..., G}(N)
+    AC2 = PreparedHalves{halftypes(Val(2), G, W)..., G}(N)
+    return (; AC, AC2)
+end
+
+prepare_operator_left(n::Val, O::Union{FiniteMPO, InfiniteMPO}, envs, GL, site, backend, allocator) =
+    prepare_operator_left(n, GL, O[site], backend, allocator)
+prepare_operator_right(n::Val, O::Union{FiniteMPO, InfiniteMPO}, envs, GR, site, backend, allocator) =
+    prepare_operator_right(n, GR, O[site], backend, allocator)
+
+assemble_operator_halves(L::AbstractTensorMap, R::AbstractTensorMap, backend, allocator) =
+    PrecomputedDerivative(L, R, backend, allocator)
+
+function hamiltonian_from_halves(n::Val{N}, site, below, envs, backend, allocator) where {N}
+    L = lefthalf(n, envs, site, below; backend, allocator)
+    R = righthalf(n, envs, site + N - 1, below; backend, allocator)
+    return assemble_operator_halves(L, R, backend, allocator)
+end
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 64f67aac0..51e77d6bf 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -124,14 +124,10 @@ function local_update!(
         ϵ_global, ϵ_trunc, decay_rate,
         iter, timeroutput, allocator
     )
-    ϵ_local = calc_galerkin(site, ψ, O, ψ, envs; alg.backend, allocator)
-
     # 1. expand
-    if !isnothing(alg.alg_expand)
-        @timeit timeroutput "expand" changebond!(
-            site, direction, ψ, O, alg.alg_expand, envs; allocator
-        )
-    end
+    ϵ_local = @timeit timeroutput "expand" galerkin_expand!(
+        site, direction, ψ, O, alg.alg_expand, envs; alg.backend, allocator
+    )
 
     # 2. local update
     alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc)
@@ -157,6 +153,31 @@ function local_update!(
     return ψ, ϵ_local, ϵ_trunc, decay_rate
 end
 
+"""
+    galerkin_expand!(site, direction, ψ, O, alg_expand, envs; backend, allocator) -> ϵ
+
+Return the Galerkin error at `site`, then expand the bond of `site` in `direction` with
+`alg_expand` (if any).
+"""
+function galerkin_expand!(site, direction, ψ, O, alg_expand, envs; backend, allocator)
+    # an expansion leaves the halves at `site` stale, so preparing them only pays off without one
+    prepare = isnothing(alg_expand)
+    ϵ = calc_galerkin(site, ψ, O, ψ, envs; prepare, backend, allocator)
+    isnothing(alg_expand) || changebond!(site, direction, ψ, O, alg_expand, envs; allocator)
+    return ϵ
+end
+# the single-site derivative is the two-site projection that the expansion forms anyway,
+# contracted with the isometry on the other site of the bond
+function galerkin_expand!(site, direction, ψ, O, alg_expand::OptimalExpand, envs; backend, allocator)
+    bond = direction === Val(:right) ? site : site - 1
+    AC2 = AC2_projection(bond, ψ, O, ψ, envs; alg_expand.backend, allocator)
+    AC´ = direction === Val(:right) ? AC2 * adjoint(_transpose_tail(ψ.AR[site + 1])) :
+        _transpose_front(adjoint(ψ.AL[site - 1]) * AC2)
+    ϵ = norm(project_complement!(normalize!(AC´), ψ.AL[site]))
+    expand_from_projection!(site, direction, ψ, AC2, alg_expand)
+    return ϵ
+end
+
 """
 $(TYPEDEF)
 
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index f5f0d5af4..8366f2522 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -59,12 +59,19 @@ Concretely, this is the overlap of the current state with the single-site deriva
 - `backend = DefaultBackend()`: backend for the tensor contractions of the derivative.
 - `allocator = DefaultAllocator()`: allocator serving their scratch space. A sweep that already
   holds one should pass it, rather than leaving this contraction to the garbage collector.
+- `prepare = false`: build the derivative at `pos` from the prepared halves cached in the
+  environments, if they hold them. Preparing costs more than one unprepared application, but an
+  effective Hamiltonian at `pos` built afterwards reuses the halves that are still current.
 """
 function calc_galerkin(
-        pos::Int, below, operator, above, envs;
+        pos::Int, below, operator, above, envs; prepare::Bool = false,
         backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
     )
-    AC´ = AC_projection(pos, below, operator, above, envs; backend, allocator)
+    AC´ = if prepare && has_halves(envs)
+        AC_hamiltonian(pos, below, operator, above, envs; backend, allocator) * above.AC[pos]
+    else
+        AC_projection(pos, below, operator, above, envs; backend, allocator)
+    end
     normalize!(AC´)
     return norm(project_complement!(AC´, below.AL[pos]))
 end
diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index 46afc6828..6c5f25496 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -4,7 +4,7 @@
 Environment manager for `FiniteMPS` and `WindowMPS`. This structure is responsible for automatically checking
 if the queried environment is still correctly cached and if not recalculates.
 """
-struct FiniteEnvironments{A, B, C, D} <: AbstractMPSEnvironments
+struct FiniteEnvironments{A, B, C, D, H} <: AbstractMPSEnvironments
     above::A
 
     operator::B #the operator
@@ -14,8 +14,38 @@ struct FiniteEnvironments{A, B, C, D} <: AbstractMPSEnvironments
 
     GLs::Vector{D}
     GRs::Vector{D}
+
+    halves::H # prepared halves of the effective Hamiltonians, see `lefthalf`
 end
 
+"""
+    PreparedHalves{L, R, G}
+
+The prepared left and right halves of the effective Hamiltonian of one window size at every site, each
+stored with the environment it was prepared from. Entries are `nothing` until first queried.
+"""
+struct PreparedHalves{L, R, G}
+    left::Vector{Union{Nothing, L}}
+    right::Vector{Union{Nothing, R}}
+    lsources::Vector{Union{Nothing, G}}
+    rsources::Vector{Union{Nothing, G}}
+    function PreparedHalves{L, R, G}(left, right, lsources, rsources) where {L, R, G}
+        return new{L, R, G}(left, right, lsources, rsources)
+    end
+end
+function PreparedHalves{L, R, G}(N::Int) where {L, R, G}
+    unset(T) = fill!(Vector{Union{Nothing, T}}(undef, N), nothing)
+    return PreparedHalves{L, R, G}(unset(L), unset(R), unset(G), unset(G))
+end
+
+"""
+    environment_halves(operator, G, N)
+
+The storage for the prepared halves of the effective Hamiltonians of `operator` on `N` sites
+with environments of type `G`, or `nothing` if `operator` has no prepared halves.
+"""
+environment_halves(operator, G, N) = nothing
+
 function initialize_environments(below::AbstractFiniteMPS, operator, above, leftstart, rightstart)
     above = above === below ? nothing : above
     N = length(below)
@@ -23,8 +53,9 @@ function initialize_environments(below::AbstractFiniteMPS, operator, above, left
     rightenvs = [i == N ? rightstart : similar(rightstart) for i in 0:N]
 
     t = similar(below.AL[1])
+    halves = environment_halves(operator, eltype(leftenvs), N)
     return FiniteEnvironments(
-        above, operator, fill(t, N), fill(t, N), leftenvs, rightenvs
+        above, operator, fill(t, N), fill(t, N), leftenvs, rightenvs, halves
     )
 end
 
@@ -157,3 +188,42 @@ function leftenv(
 
     return ca.GLs[ind]
 end
+
+"""
+    lefthalf(Val(N), envs, site, state; backend, allocator)
+    righthalf(Val(N), envs, site, state; backend, allocator)
+
+Return the prepared half of the effective Hamiltonian of an `N`-site window (`AC_hamiltonian`
+for `N = 1`, `AC2_hamiltonian` for `N = 2`) that depends only on the left (right) environment
+of `site` for `state`, where `site` is the first (last) site of the window.
+A half is prepared when first queried, and again whenever its environment has been replaced.
+
+See also [`leftenv`](@ref) and [`rightenv`](@ref).
+"""
+function lefthalf(
+        n::Val, envs::FiniteEnvironments, site::Int, state;
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    GL = leftenv(envs, site, state; backend, allocator)
+    halves = halfstorage(envs, n)
+    if halves.lsources[site] !== GL
+        halves.left[site] = prepare_operator_left(n, envs.operator, envs, GL, site, backend, allocator)
+        halves.lsources[site] = GL
+    end
+    return something(halves.left[site])
+end
+function righthalf(
+        n::Val, envs::FiniteEnvironments, site::Int, state;
+        backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()
+    )
+    GR = rightenv(envs, site, state; backend, allocator)
+    halves = halfstorage(envs, n)
+    if halves.rsources[site] !== GR
+        halves.right[site] = prepare_operator_right(n, envs.operator, envs, GR, site, backend, allocator)
+        halves.rsources[site] = GR
+    end
+    return something(halves.right[site])
+end
+
+has_halves(envs) = false
+has_halves(envs::FiniteEnvironments) = !isnothing(envs.halves)
diff --git a/test/groundstate/approximate.jl b/test/groundstate/approximate.jl
index 14e8d6801..567d9633e 100644
--- a/test/groundstate/approximate.jl
+++ b/test/groundstate/approximate.jl
@@ -132,6 +132,16 @@ end
         @test norm(O * ψ₁ - ψ₂) ≈ 0 atol = 0.001
     end
 
+    @testset "MPOHamiltonian * finitemps1 ≈ finitemps2" for alg in finite_algs
+        L = 6
+        Random.seed!(5)
+        H = long_range_ising(Float64; g = 2.0, L)
+        ϕ = FiniteMPS(randn, ComplexF64, L, ℂ^2, ℂ^2)
+        ψ, = approximate(FiniteMPS(randn, ComplexF64, L, ℂ^2, ℂ^8), (H, ϕ), alg)
+        Hϕ = FiniteMPO(H) * ϕ
+        @test abs(dot(ψ, Hϕ)) ≈ norm(ψ) * norm(Hϕ) rtol = 1.0e-8
+    end
+
     @testset "Finite MPO-MPS zip-up $(spacetype(pspace))" for (pspace, Dspace, _) in zipup_spacelist
         L = 6
 
diff --git a/test/groundstate/groundstate.jl b/test/groundstate/groundstate.jl
index 58e349692..55f6f66b4 100644
--- a/test/groundstate/groundstate.jl
+++ b/test/groundstate/groundstate.jl
@@ -23,6 +23,7 @@ verbosity_conv = 1
     models = [
         "nearest-neighbour" => force_planar(transverse_field_ising(; g, L)),
         "long-range, real scalartype" => force_planar(long_range_ising(Float64; g, L)),
+        "three-site" => force_planar(three_site_ising(; g, L)),
     ]
 
     @testset "$name" for (name, H) in models
diff --git a/test/hamiltonian/derivatives.jl b/test/hamiltonian/derivatives.jl
index 246637444..5459bbb81 100644
--- a/test/hamiltonian/derivatives.jl
+++ b/test/hamiltonian/derivatives.jl
@@ -10,6 +10,7 @@ using MPSKit
 using MPSKit: C_hamiltonian, AC_hamiltonian, AC2_hamiltonian
 using MPSKit: _transpose_front, _transpose_tail
 using TensorKit
+using TensorKit: ℙ
 using BlockTensorKit: nonzero_length
 using Random
 
@@ -120,3 +121,82 @@ end
         end
     end
 end
+
+two_site_window(ψ, site) = _transpose_front(ψ.AC[site]) * _transpose_tail(ψ.AR[site + 1])
+
+@testset "FiniteMPS derivatives follow changes to the state" verbose = true begin
+    L = 6
+    H1, H2 = transverse_field_ising(; g = 2.0, L), long_range_ising(Float64; g = 1.0, L)
+    models = [
+        "nearest-neighbour" => (H1, ℂ^3),
+        "long-range" => (H2, ℂ^3),
+        "three-site" => (three_site_ising(; L), ℂ^3),
+        "planar" => (force_planar(three_site_ising(; L)), ℙ^3),
+        "SU2" => (long_range_heisenberg(; L), Rep[SU₂](0 => 2, 1 // 2 => 2, 1 => 1)),
+        "FiniteMPO" => (FiniteMPO(three_site_ising(; L)), ℂ^3),
+        "LazySum" => (LazySum([H1, MultipliedOperator(H2, 1.3)]), ℂ^3),
+    ]
+
+    @testset "$name" for (name, (H, V)) in models
+        P = physicalspace(H isa LazySum ? first(H) : H, 1)
+        ψ = FiniteMPS(randn, ComplexF64, L, P, V)
+        envs = environments(ψ, H, ψ)
+        for k in (nothing, 3, 1, L)
+            isnothing(k) || (ψ.AC[k] = randn(scalartype(ψ), space(ψ.AC[k])))
+            for site in (5, 1, 4, 2, 3)
+                AC = ψ.AC[site]
+                @test AC_hamiltonian(site, ψ, H, ψ, envs)(AC) ≈
+                    AC_hamiltonian(site, ψ, H, ψ, envs; prepare = false)(AC)
+                AC2 = two_site_window(ψ, site)
+                @test AC2_hamiltonian(site, ψ, H, ψ, envs)(AC2) ≈
+                    AC2_hamiltonian(site, ψ, H, ψ, envs; prepare = false)(AC2)
+            end
+        end
+    end
+
+    @testset "WindowMPS" begin
+        H = force_planar(transverse_field_ising(; g = 2.0))
+        ψ = WindowMPS(InfiniteMPS(randn, ComplexF64, [ℙ^2], [ℙ^3]), 4)
+        Hw = WindowMPOHamiltonian(H, 1:4)
+        envs = environments(ψ, Hw)
+        for site in 1:3
+            AC2 = two_site_window(ψ, site)
+            @test AC2_hamiltonian(site, ψ, Hw, ψ, envs)(AC2) ≈
+                AC2_hamiltonian(site, ψ, Hw, ψ, envs; prepare = false)(AC2)
+        end
+    end
+end
+
+@testset "FiniteMPS derivatives with different bra and ket" verbose = true begin
+    L = 6
+    models = [
+        "nearest-neighbour" => (transverse_field_ising(; g = 2.0, L), ℂ^2, ℂ^3),
+        "long-range" => (long_range_ising(Float64; g = 2.0, L), ℂ^2, ℂ^3),
+        "three-site" => (three_site_ising(; L), ℂ^2, ℂ^3),
+        "planar" => (force_planar(three_site_ising(; L)), ℙ^2, ℙ^3),
+        "real" => (transverse_field_ising(Float64; g = 2.0, L), ℂ^2, ℂ^3),
+        "fermionic" => (kitaev_model(ComplexF64, Trivial; L), Vect[FermionParity](0 => 2, 1 => 2), Vect[FermionParity](0 => 3, 1 => 3)),
+        "SU2" => (long_range_heisenberg(; L), Rep[SU₂](0 => 2, 1 // 2 => 2), Rep[SU₂](0 => 3, 1 // 2 => 3)),
+    ]
+
+    @testset "$name" for (name, (H, Vbra, Vket)) in models
+        T, P = scalartype(H), physicalspace(H, 1)
+        bra = FiniteMPS(randn, T, L, P, Vbra)
+        ket = FiniteMPS(randn, T, L, P, Vket)
+        envs = environments(bra, H, ket)
+        # dense MPOs make no assumptions about the bra and ket
+        reference = environments(bra, FiniteMPO(H), ket)
+        for site in 1:L
+            AC = ket.AC[site]
+            y = AC_hamiltonian(site, bra, FiniteMPO(H), ket, reference; prepare = false)(AC)
+            @test AC_hamiltonian(site, bra, H, ket, envs; prepare = false)(AC) ≈ y
+            @test AC_hamiltonian(site, bra, H, ket, envs)(AC) ≈ y
+        end
+        for site in 1:(L - 1)
+            AC2 = two_site_window(ket, site)
+            y = AC2_hamiltonian(site, bra, FiniteMPO(H), ket, reference; prepare = false)(AC2)
+            @test AC2_hamiltonian(site, bra, H, ket, envs; prepare = false)(AC2) ≈ y
+            @test AC2_hamiltonian(site, bra, H, ket, envs)(AC2) ≈ y
+        end
+    end
+end
diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl
index 7fbf4badf..19981515c 100644
--- a/test/setup/testsetup.jl
+++ b/test/setup/testsetup.jl
@@ -21,7 +21,7 @@ export f_plus_f_min, f_min_f_plus, f_num, f_hopping
 export force_planar
 export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model, kitaev_model
-export long_range_ising, long_range_ising_infinite
+export long_range_ising, long_range_ising_infinite, three_site_ising, long_range_heisenberg
 export classical_ising_tensors, classical_ising, sixvertex
 export bad_initial_state
 export SCHEDULERS, with_scheduler
@@ -223,6 +223,21 @@ function long_range_ising_infinite(
     return InfiniteMPOHamiltonian(lattice, ((i,) => -g * X for i in 1:L)..., (1, L) => -ZZ)
 end
 
+# three-site terms give non-trivial A·A, C·A and A·B channels in the two-site derivatives
+function three_site_ising(T::Type{<:Number} = ComplexF64; g = 2.0, L = 6)
+    X = S_x(T, Trivial; spin = 1 // 2)
+    return transverse_field_ising(T; g, L) + FiniteMPOHamiltonian(
+        fill(space(X, 1), L), (1, 2, 5) => X ⊗ X ⊗ X, (2, 3, 6) => X ⊗ X ⊗ X, (1, 4) => X ⊗ X, (3, 6) => X ⊗ X
+    )
+end
+
+# a Heisenberg chain with longer-range exchanges, which give the SU(2) MPO non-trivial continuing channels
+function long_range_heisenberg(T::Type{<:Number} = ComplexF64, sym::Type{<:Sector} = SU2Irrep; L = 6)
+    H = heisenberg_XXX(T, sym; spin = 1 // 2, L)
+    h = S_exchange(T, sym; spin = 1 // 2)
+    return H + FiniteMPOHamiltonian(fill(space(h, 1), L), (1, 3) => h, (2, 5) => h)
+end
+
 function ising_bond_tensor(β)
     J = 1.0
     K = β * J