From 3396ad8707d7865fe679173a2ed87b2bade6e10b Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Tue, 21 Apr 2026 13:55:39 +0200
Subject: [PATCH 01/13] windowhamiltonian poc

---
 src/MPSKit.jl                          |  3 ++-
 src/algorithms/propagator/corvector.jl |  6 ++---
 src/environments/finite_envs.jl        | 10 +++++++
 src/operators/windowhamiltonian.jl     | 36 ++++++++++++++++++++++++++
 test/states/windowmps.jl               |  7 ++---
 5 files changed, 55 insertions(+), 7 deletions(-)
 create mode 100644 src/operators/windowhamiltonian.jl

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index b29174216..694dd8024 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -19,7 +19,7 @@ export QP, LeftGaugedQP, RightGaugedQP
 export AbstractMPO
 export MPO, FiniteMPO, InfiniteMPO
 export JordanMPOTensor
-export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian
+export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian, WindowMPOHamiltonian
 export MultilineMPO
 export UntimedOperator, TimedOperator, MultipliedOperator, LazySum
 
@@ -123,6 +123,7 @@ include("operators/abstractmpo.jl")
 include("operators/mpo.jl")
 include("operators/jordanmpotensor.jl")
 include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
+include("operators/windowhamiltonian.jl")
 include("operators/ortho.jl")
 include("operators/multilinempo.jl")
 include("operators/projection.jl")
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index de04e8078..be457cf76 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -59,7 +59,7 @@ See also [`Jeckelmann`](@ref) for the original approach.
 struct NaiveInvert <: DDMRG_Flavour end
 
 function propagator(
-        A::AbstractFiniteMPS, z::Number, H::FiniteMPOHamiltonian,
+        A::AbstractFiniteMPS, z::Number, H,
         alg::DynamicalDMRG{NaiveInvert}; init = copy(A)
     )
     h_envs = environments(init, H, init) # environments for h
@@ -119,7 +119,7 @@ See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 struct Jeckelmann <: DDMRG_Flavour end
 
 function propagator(
-        A::AbstractFiniteMPS, z, H::FiniteMPOHamiltonian,
+        A::AbstractFiniteMPS, z, H,
         alg::DynamicalDMRG{Jeckelmann}; init = copy(A)
     )
     ω = real(z)
@@ -176,7 +176,7 @@ function propagator(
 end
 
 function squaredenvs(
-        state::AbstractFiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H, state)
+        state::AbstractFiniteMPS, H, envs = environments(state, H, state)
     )
     H² = conj(H) * H
     L = length(state)
diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index fee1555bb..1b4986a0b 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -56,6 +56,16 @@ function environments(
     )
     return initialize_environments(below, operator, above, leftstart, rightstart)
 end
+function environments(
+        below::WindowMPS, O::WindowMPOHamiltonian, above = nothing;
+        lenvs = environments(below.left_gs, O.left_ham, below.left_gs),
+        renvs = environments(below.right_gs, O.right_ham, below.right_gs)
+    )
+    leftstart = copy(lenvs.GLs[1])
+    rightstart = copy(renvs.GRs[end])
+
+    return environments(below, O.finite_ham, above; leftstart, rightstart)
+end
 
 environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} =
     environments(below, nothing, above)
diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
new file mode 100644
index 000000000..7eed2c458
--- /dev/null
+++ b/src/operators/windowhamiltonian.jl
@@ -0,0 +1,36 @@
+"""
+A WindowMPS is a finite MPS embedded between an infinite mps to the left, and an inifite mps to the right. 
+The associated hamiltonian has also an infinite part to the left, a finite hamiltonian in the middle, and an infinite part to the right.
+
+Acts simalar as just a finite hamiltonian, but we 'remember' the boundary hamiltonians.
+"""
+
+# todo - what is the required interface for abstractmpo?
+# support densempo windows?
+struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
+    left_ham :: InfiniteMPOHamiltonian{O}
+    finite_ham :: FiniteMPOHamiltonian{O}
+    right_ham :: InfiniteMPOHamiltonian{O}
+end
+
+#utility constructor
+function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
+    
+    # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
+    left_edge = (interval.start-1) % length(ham)
+    left_ham = InfiniteMPOHamiltonian([ham[i] for i in (left_edge-length(ham)+1):left_edge])
+    right_edge = (interval.stop+1)%length(ham)
+    right_ham = InfiniteMPOHamiltonian([ham[i] for i in right_edge:(right_edge+length(ham)-1)])
+
+    finite_ham = FiniteMPOHamiltonian([ham[i] for i in  interval])
+    WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
+end
+#
+Base.parent(h::WindowMPOHamiltonian) = h.finite_ham
+
+Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham))
+
+# some basic linalg
+for fun in (:(Base.:+), :(Base.:-), :(Base.:*))
+    @eval $fun(a::WindowMPOHamiltonian,b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham,b.left_ham),$fun(a.finite_ham,b.finite_ham),$fun(a.right_ham,b.right_ham))
+end
diff --git a/test/states/windowmps.jl b/test/states/windowmps.jl
index 3cbddf09d..2106dab0e 100644
--- a/test/states/windowmps.jl
+++ b/test/states/windowmps.jl
@@ -63,14 +63,15 @@ using TensorKit: ℙ
 
     e1 = expectation_value(window, (2, 3) => O)
 
-    window, envs, _ = find_groundstate(window, ham, DMRG(; verbosity = 0))
+    w_ham = WindowMPOHamiltonian(ham,1:length(window))
+    window, envs, _ = find_groundstate(window, w_ham, DMRG(; verbosity = 0))
 
     e2 = expectation_value(window, (2, 3) => O)
 
     @test real(e2) ≤ real(e1)
 
-    window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme = truncrank(20)), envs)
-    window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs)
+    window, envs = timestep(window, w_ham, 0.1, 0.0, TDVP2(; trscheme = truncrank(20)), envs)
+    window, envs = timestep(window, w_ham, 0.1, 0.0, TDVP(), envs)
 
     e3 = expectation_value(window, (2, 3) => O)
 

From ea455deee17c5d7bd959c20d7305de7eeb76b1ed Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Tue, 28 Apr 2026 20:29:27 +0200
Subject: [PATCH 02/13] test

---
 src/algorithms/expval.jl           |  7 +++++++
 src/operators/windowhamiltonian.jl |  5 +++++
 test/algorithms/dynamical_dmrg.jl  | 31 +++++++++++++++++++++++++++++-
 3 files changed, 42 insertions(+), 1 deletion(-)

diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index c147ba391..d60213af2 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -128,6 +128,13 @@ function expectation_value(
     return dot(ψ, H, ψ, envs) / dot(ψ, ψ)
 end
 
+function expectation_value(
+        ψ::WindowMPS, H::WindowMPOHamiltonian,
+        envs::AbstractMPSEnvironments = environments(ψ, H)
+    )
+    return dot(ψ, H, ψ, envs) / dot(ψ, ψ)
+end
+
 function expectation_value(
         ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
         envs::AbstractMPSEnvironments = environments(ψ, H, ψ)
diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index 7eed2c458..e9ee8c629 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -34,3 +34,8 @@ Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy
 for fun in (:(Base.:+), :(Base.:-), :(Base.:*))
     @eval $fun(a::WindowMPOHamiltonian,b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham,b.left_ham),$fun(a.finite_ham,b.finite_ham),$fun(a.right_ham,b.right_ham))
 end
+
+TensorKit.dot(
+        bra::WindowMPS, H::WindowMPOHamiltonian, ket::WindowMPS = bra,
+        envs = environments(bra, H, ket)
+        ) = dot(bra.window, H.finite_ham, ket.window,envs)
diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl
index 2b7a821b7..8fb1d712c 100644
--- a/test/algorithms/dynamical_dmrg.jl
+++ b/test/algorithms/dynamical_dmrg.jl
@@ -12,7 +12,7 @@ using TensorKit: ℙ
 
 verbosity_conv = 1
 
-@testset "Dynamical DMRG" verbose = true begin
+@testset "Dynamical DMRG FiniteMPS" verbose = true begin
     L = 10
     H = force_planar(-transverse_field_ising(; L, g = -4))
     gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^10), H; verbosity = verbosity_conv)
@@ -32,3 +32,32 @@ verbosity_conv = 1
         @test data ≈ predicted atol = 1.0e-8
     end
 end
+
+
+@testset "Dynamical DMRG WindowMPS" verbose = true begin
+    N = 20
+
+    H = transverse_field_ising(g = -4)
+    Ω = InfiniteMPS(ComplexSpace(2),ComplexSpace(20))
+
+    (Ω,_) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv))
+    XΩ = WindowMPS(Ω,N)
+    H_w = WindowMPOHamiltonian(H,1:N)
+    
+
+    gs_en =  expectation_value(XΩ, H_w)
+    vals = range(gs_en - 1.0,gs_en+1.0, length=5)
+    eta = 0.3im
+    predicted = [1 / (v + eta - gs_en) for v in vals]
+
+
+    @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert())
+        alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8)
+        data = map(vals) do v
+            result, = propagator(XΩ, v + eta, H_w, alg)
+            return result
+        end
+        @test data ≈ predicted atol = 1.0e-8
+    end
+end
+

From c8799d2852c318771b1f004b638d1c889b441b3b Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Wed, 29 Apr 2026 10:07:07 +0200
Subject: [PATCH 03/13] remove parent

---
 src/operators/windowhamiltonian.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index e9ee8c629..030d6c59a 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -25,8 +25,7 @@ function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
     finite_ham = FiniteMPOHamiltonian([ham[i] for i in  interval])
     WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
 end
-#
-Base.parent(h::WindowMPOHamiltonian) = h.finite_ham
+
 
 Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham))
 

From 92acdedfcaf2b45529f884ddf936481fcb762614 Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Wed, 29 Apr 2026 10:07:55 +0200
Subject: [PATCH 04/13] format

---
 src/operators/windowhamiltonian.jl | 28 ++++++++++++++--------------
 test/algorithms/dynamical_dmrg.jl  | 14 +++++++-------
 test/states/windowmps.jl           |  2 +-
 3 files changed, 22 insertions(+), 22 deletions(-)

diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index 030d6c59a..98f72390f 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -8,22 +8,22 @@ Acts simalar as just a finite hamiltonian, but we 'remember' the boundary hamilt
 # todo - what is the required interface for abstractmpo?
 # support densempo windows?
 struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
-    left_ham :: InfiniteMPOHamiltonian{O}
-    finite_ham :: FiniteMPOHamiltonian{O}
-    right_ham :: InfiniteMPOHamiltonian{O}
+    left_ham::InfiniteMPOHamiltonian{O}
+    finite_ham::FiniteMPOHamiltonian{O}
+    right_ham::InfiniteMPOHamiltonian{O}
 end
 
 #utility constructor
 function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
-    
+
     # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
-    left_edge = (interval.start-1) % length(ham)
-    left_ham = InfiniteMPOHamiltonian([ham[i] for i in (left_edge-length(ham)+1):left_edge])
-    right_edge = (interval.stop+1)%length(ham)
-    right_ham = InfiniteMPOHamiltonian([ham[i] for i in right_edge:(right_edge+length(ham)-1)])
+    left_edge = (interval.start - 1) % length(ham)
+    left_ham = InfiniteMPOHamiltonian([ham[i] for i in (left_edge - length(ham) + 1):left_edge])
+    right_edge = (interval.stop + 1) % length(ham)
+    right_ham = InfiniteMPOHamiltonian([ham[i] for i in right_edge:(right_edge + length(ham) - 1)])
 
-    finite_ham = FiniteMPOHamiltonian([ham[i] for i in  interval])
-    WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
+    finite_ham = FiniteMPOHamiltonian([ham[i] for i in interval])
+    return WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
 end
 
 
@@ -31,10 +31,10 @@ Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy
 
 # some basic linalg
 for fun in (:(Base.:+), :(Base.:-), :(Base.:*))
-    @eval $fun(a::WindowMPOHamiltonian,b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham,b.left_ham),$fun(a.finite_ham,b.finite_ham),$fun(a.right_ham,b.right_ham))
+    @eval $fun(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham, b.left_ham), $fun(a.finite_ham, b.finite_ham), $fun(a.right_ham, b.right_ham))
 end
 
 TensorKit.dot(
-        bra::WindowMPS, H::WindowMPOHamiltonian, ket::WindowMPS = bra,
-        envs = environments(bra, H, ket)
-        ) = dot(bra.window, H.finite_ham, ket.window,envs)
+    bra::WindowMPS, H::WindowMPOHamiltonian, ket::WindowMPS = bra,
+    envs = environments(bra, H, ket)
+) = dot(bra.window, H.finite_ham, ket.window, envs)
diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl
index 8fb1d712c..146c43af3 100644
--- a/test/algorithms/dynamical_dmrg.jl
+++ b/test/algorithms/dynamical_dmrg.jl
@@ -38,15 +38,15 @@ end
     N = 20
 
     H = transverse_field_ising(g = -4)
-    Ω = InfiniteMPS(ComplexSpace(2),ComplexSpace(20))
+    Ω = InfiniteMPS(ComplexSpace(2), ComplexSpace(20))
 
-    (Ω,_) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv))
-    XΩ = WindowMPS(Ω,N)
-    H_w = WindowMPOHamiltonian(H,1:N)
-    
+    (Ω, _) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv))
+    XΩ = WindowMPS(Ω, N)
+    H_w = WindowMPOHamiltonian(H, 1:N)
 
-    gs_en =  expectation_value(XΩ, H_w)
-    vals = range(gs_en - 1.0,gs_en+1.0, length=5)
+
+    gs_en = expectation_value(XΩ, H_w)
+    vals = range(gs_en - 1.0, gs_en + 1.0, length = 5)
     eta = 0.3im
     predicted = [1 / (v + eta - gs_en) for v in vals]
 
diff --git a/test/states/windowmps.jl b/test/states/windowmps.jl
index 2106dab0e..ab37bbab8 100644
--- a/test/states/windowmps.jl
+++ b/test/states/windowmps.jl
@@ -63,7 +63,7 @@ using TensorKit: ℙ
 
     e1 = expectation_value(window, (2, 3) => O)
 
-    w_ham = WindowMPOHamiltonian(ham,1:length(window))
+    w_ham = WindowMPOHamiltonian(ham, 1:length(window))
     window, envs, _ = find_groundstate(window, w_ham, DMRG(; verbosity = 0))
 
     e2 = expectation_value(window, (2, 3) => O)

From a1ed562cd721ea09abcdd22c71fbe18dc956cd7a Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Fri, 1 May 2026 09:35:06 +0200
Subject: [PATCH 05/13] ::Number

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

diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index be457cf76..86439dc0f 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -119,7 +119,7 @@ See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 struct Jeckelmann <: DDMRG_Flavour end
 
 function propagator(
-        A::AbstractFiniteMPS, z, H,
+        A::AbstractFiniteMPS, z::Number, H,
         alg::DynamicalDMRG{Jeckelmann}; init = copy(A)
     )
     ω = real(z)

From bdf1bce9da3a6fd4d712f25c193f412f21437579 Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Fri, 1 May 2026 09:52:33 +0200
Subject: [PATCH 06/13] constructor in line with windowhamiltonian

---
 src/states/windowmps.jl | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl
index 5c0bab9d0..fa23acdf1 100644
--- a/src/states/windowmps.jl
+++ b/src/states/windowmps.jl
@@ -106,18 +106,27 @@ function WindowMPS(N::Int, V::VectorSpace, args...; kwargs...)
     return WindowMPS(fill(V, N), args...; kwargs...)
 end
 
-function WindowMPS(ψ::InfiniteMPS{A, B}, L::Int) where {A, B}
+
+WindowMPS(ψ::InfiniteMPS, L::Int) = WindowMPS(ψ, 1:L)
+
+function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A,B}
+
+    # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
+    left_edge = (interval.start - 1) % length(ψ)
+    right_edge = (interval.stop + 1) % length(ψ)
+
+    L = length(interval)
     CLs = Vector{Union{Missing, B}}(missing, L + 1)
     ALs = Vector{Union{Missing, A}}(missing, L)
     ARs = Vector{Union{Missing, A}}(missing, L)
     ACs = Vector{Union{Missing, A}}(missing, L)
 
-    ALs .= ψ.AL[1:L]
-    ARs .= ψ.AR[1:L]
-    ACs .= ψ.AC[1:L]
-    CLs .= ψ.C[0:L]
+    ALs .= ψ.AL[interval]
+    ARs .= ψ.AR[interval]
+    ACs .= ψ.AC[interval]
+    CLs .= ψ.C[interval.start-1:interval.stop]
 
-    return WindowMPS(ψ, FiniteMPS(ALs, ARs, ACs, CLs), ψ)
+    return WindowMPS(circshift(ψ, -left_edge), FiniteMPS(ALs, ARs, ACs, CLs), circshift(ψ,-right_edge + 1))
 end
 
 #===========================================================================================

From 11cb2d5bba44ebcabd35c03f291376c748ac72d6 Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Fri, 1 May 2026 09:56:54 +0200
Subject: [PATCH 07/13] circshift

---
 src/operators/mpohamiltonian.jl    | 2 +-
 src/operators/windowhamiltonian.jl | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index 37c52b9fb..f89d46a9f 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -707,7 +707,7 @@ end
 function Base.similar(H::MPOHamiltonian, ::Type{TorA}) where {TorA <: Union{Number, DenseVector}}
     return MPOHamiltonian(similar.(parent(H), TorA))
 end
-
+Base.circshift(H::InfiniteMPOHamiltonian, shift::Integer) = InfiniteMPOHamiltonian(circshift(parent(copy(H)), shift))
 # Linear Algebra
 # --------------
 function Base.:+(
diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index 98f72390f..071e0eca3 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -18,15 +18,15 @@ function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
 
     # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
     left_edge = (interval.start - 1) % length(ham)
-    left_ham = InfiniteMPOHamiltonian([ham[i] for i in (left_edge - length(ham) + 1):left_edge])
+    left_ham = circshift(ham, -left_edge)
     right_edge = (interval.stop + 1) % length(ham)
-    right_ham = InfiniteMPOHamiltonian([ham[i] for i in right_edge:(right_edge + length(ham) - 1)])
+    right_ham = circshift(ham, -right_edge + 1)
 
     finite_ham = FiniteMPOHamiltonian([ham[i] for i in interval])
     return WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
 end
 
-
+Base.parent(h::WindowMPOHamiltonian) = h.finite_ham
 Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham))
 
 # some basic linalg

From 2dc3c8955e93b28b703b35872b053712adab0bd4 Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Fri, 1 May 2026 09:57:15 +0200
Subject: [PATCH 08/13] format

---
 src/states/windowmps.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl
index fa23acdf1..7c5d48327 100644
--- a/src/states/windowmps.jl
+++ b/src/states/windowmps.jl
@@ -109,7 +109,7 @@ end
 
 WindowMPS(ψ::InfiniteMPS, L::Int) = WindowMPS(ψ, 1:L)
 
-function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A,B}
+function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A, B}
 
     # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
     left_edge = (interval.start - 1) % length(ψ)
@@ -124,9 +124,9 @@ function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A,B}
     ALs .= ψ.AL[interval]
     ARs .= ψ.AR[interval]
     ACs .= ψ.AC[interval]
-    CLs .= ψ.C[interval.start-1:interval.stop]
+    CLs .= ψ.C[(interval.start - 1):interval.stop]
 
-    return WindowMPS(circshift(ψ, -left_edge), FiniteMPS(ALs, ARs, ACs, CLs), circshift(ψ,-right_edge + 1))
+    return WindowMPS(circshift(ψ, -left_edge), FiniteMPS(ALs, ARs, ACs, CLs), circshift(ψ, -right_edge + 1))
 end
 
 #===========================================================================================

From e75ae1fc9ae11b879c9cc267008b06ecdfbf0d08 Mon Sep 17 00:00:00 2001
From: maartenvd <maartenvd1994@gmail.com>
Date: Mon, 4 May 2026 14:43:42 +0200
Subject: [PATCH 09/13] attempt at better doc

---
 src/algorithms/propagator/corvector.jl | 23 +++++++++++++++--------
 1 file changed, 15 insertions(+), 8 deletions(-)

diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index 86439dc0f..02ea4cf20 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -35,7 +35,7 @@ end
 """
     propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init=copy(ψ₀))
 
-Calculate the propagator ``\\frac{1}{E₀ + z - H}|ψ₀⟩`` using the dynamical DMRG
+Calculate the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
 algorithm.
 """
 function propagator end
@@ -47,13 +47,11 @@ An alternative approach to the dynamical DMRG algorithm, without quadratic terms
 less controlled approximation.
 This algorithm minimizes the following cost function
 ```math
-⟨ψ|(H - E)|ψ⟩ - ⟨ψ|ψ₀⟩ - ⟨ψ₀|ψ⟩
-```
-which is equivalent to the original approach if
-```math
-|ψ₀⟩ = (H - E)|ψ⟩
+⟨ψ|(z - H)|ψ⟩ - ⟨ψ|ψ₀⟩ - ⟨ψ₀|ψ⟩
 ```
 
+Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩.
+
 See also [`Jeckelmann`](@ref) for the original approach.
 """
 struct NaiveInvert <: DDMRG_Flavour end
@@ -105,11 +103,20 @@ end
 """
 $(TYPEDEF)
 
-The original flavour of dynamical DMRG, which minimizes the following (quadratic) cost function:
+The original flavour of dynamical DMRG, which minimizes functional (14) from Jeckelmann2002.
+```math
+|| <ψ| (Re(z) - H)^2 + Im(z)^2 |ψ⟩ +Im(z) (<ψ₀|ψ⟩+<ψ|ψ₀⟩ ||
+```
+
+This would achieve a minima at
 ```math
-|| (H - E) |ψ₀⟩ - |ψ⟩ ||
+-Im(z) |ψ₀⟩ = ((Re(z) - H)^2 + Im(z)^2)|ψ⟩
 ```
 
+Together with equation (11) from that same paper we can determine the full propagator (z-H)^-1 |ψ₀>.
+
+Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩.
+
 See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 
 ## References

From a0b4e371d2757158c36c927edcba82753ed05772 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Tue, 30 Jun 2026 19:19:59 -0400
Subject: [PATCH 10/13] Fix WindowMPOHamiltonian environments to build finite
 envs directly

The WindowMPOHamiltonian environments method must call
initialize_environments on the finite Hamiltonian (there is no
environments(::WindowMPS, ::FiniteMPOHamiltonian, ...) method), and
expose leftstart/rightstart as keywords so squaredenvs (Jeckelmann
flavour) can pass squared boundary environments through.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/environments/finite_envs.jl | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index 1b4986a0b..785bec24d 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -59,12 +59,11 @@ end
 function environments(
         below::WindowMPS, O::WindowMPOHamiltonian, above = nothing;
         lenvs = environments(below.left_gs, O.left_ham, below.left_gs),
-        renvs = environments(below.right_gs, O.right_ham, below.right_gs)
+        renvs = environments(below.right_gs, O.right_ham, below.right_gs),
+        leftstart = copy(lenvs.GLs[1]),
+        rightstart = copy(renvs.GRs[end])
     )
-    leftstart = copy(lenvs.GLs[1])
-    rightstart = copy(renvs.GRs[end])
-
-    return environments(below, O.finite_ham, above; leftstart, rightstart)
+    return initialize_environments(below, O.finite_ham, above, leftstart, rightstart)
 end
 
 environments(below::S, above::S) where {S <: Union{FiniteMPS, WindowMPS}} =

From 75eda68f6d09d14b2d750950b00d04a582924308 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 1 Jul 2026 08:44:52 -0400
Subject: [PATCH 11/13] Fix WindowMPOHamiltonian linalg, squared-env warnings,
 and docs

- Fix `+`/`-` on WindowMPOHamiltonian: the finite window carries full Jordan
  virtual spaces on its boundaries, so summands must be block-diagonalized at
  every site (including edges), like InfiniteMPOHamiltonian addition. Add a
  dedicated `_add_finite_window` helper and a `VectorInterface.scale` method so
  unary `-`, scalar `*`/`/` work through the AbstractMPO fallbacks.
- Avoid spurious non-convergent fixed-point solves in `squaredenvs`: inline the
  `leftstart`/`rightstart` defaults in the WindowMPOHamiltonian environments
  method so the infinite environment of e.g. `conj(H) * H` is only computed when
  no explicit boundaries are passed.
- Attach and clean up the WindowMPOHamiltonian docstring (was silently detached),
  document the interval constructors of WindowMPOHamiltonian and WindowMPS, and
  fix the propagator/DynamicalDMRG docstrings (return value + LaTeX).
- Add test/operators/windowhamiltonian.jl covering non-trivial intervals
  (circshift alignment via offset invariance) and the linear algebra.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/algorithms/propagator/corvector.jl | 20 +++---
 src/environments/finite_envs.jl        | 10 +--
 src/operators/windowhamiltonian.jl     | 85 +++++++++++++++++++++++---
 src/states/windowmps.jl                |  7 +++
 test/operators/windowhamiltonian.jl    | 47 ++++++++++++++
 5 files changed, 148 insertions(+), 21 deletions(-)
 create mode 100644 test/operators/windowhamiltonian.jl

diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index 02ea4cf20..93bfd6637 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -35,8 +35,11 @@ end
 """
     propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init=copy(ψ₀))
 
-Calculate the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
+Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
 algorithm.
+
+Returns a tuple `(g, ψ)` where `g` is the approximation of the propagator matrix element
+``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and `ψ` is the MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``.
 """
 function propagator end
 
@@ -50,7 +53,7 @@ This algorithm minimizes the following cost function
 ⟨ψ|(z - H)|ψ⟩ - ⟨ψ|ψ₀⟩ - ⟨ψ₀|ψ⟩
 ```
 
-Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩.
+Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
 
 See also [`Jeckelmann`](@ref) for the original approach.
 """
@@ -104,18 +107,19 @@ end
 $(TYPEDEF)
 
 The original flavour of dynamical DMRG, which minimizes functional (14) from Jeckelmann2002.
+Writing ``ω = \\mathrm{Re}(z)`` and ``η = \\mathrm{Im}(z)``, this is
 ```math
-|| <ψ| (Re(z) - H)^2 + Im(z)^2 |ψ⟩ +Im(z) (<ψ₀|ψ⟩+<ψ|ψ₀⟩ ||
+W(ψ) = ⟨ψ|(ω - H)^2 + η^2|ψ⟩ + η(⟨ψ₀|ψ⟩ + ⟨ψ|ψ₀⟩)
 ```
-
-This would achieve a minima at
+which attains its minimum at
 ```math
--Im(z) |ψ₀⟩ = ((Re(z) - H)^2 + Im(z)^2)|ψ⟩
+((ω - H)^2 + η^2)|ψ⟩ = -η|ψ₀⟩
 ```
 
-Together with equation (11) from that same paper we can determine the full propagator (z-H)^-1 |ψ₀>.
+Together with equation (11) from that same paper we can determine the full propagator
+``\\frac{1}{z - H}|ψ₀⟩``.
 
-Returns the approximation of <ψ₀| (z-H)^-1 |ψ₀> and |ψ⟩.
+Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
 
 See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 
diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index 785bec24d..cb11cd53f 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -58,10 +58,12 @@ function environments(
 end
 function environments(
         below::WindowMPS, O::WindowMPOHamiltonian, above = nothing;
-        lenvs = environments(below.left_gs, O.left_ham, below.left_gs),
-        renvs = environments(below.right_gs, O.right_ham, below.right_gs),
-        leftstart = copy(lenvs.GLs[1]),
-        rightstart = copy(renvs.GRs[end])
+        # NOTE: the defaults are deliberately inlined rather than shared via `lenvs`/`renvs`
+        # kwargs: when explicit boundaries are passed (e.g. the squared environments in
+        # `squaredenvs`), we must not trigger the infinite fixed-point solve, which need not
+        # converge for e.g. `conj(H) * H`.
+        leftstart = copy(environments(below.left_gs, O.left_ham, below.left_gs).GLs[1]),
+        rightstart = copy(environments(below.right_gs, O.right_ham, below.right_gs).GRs[end])
     )
     return initialize_environments(below, O.finite_ham, above, leftstart, rightstart)
 end
diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index 071e0eca3..15ecd552d 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -1,19 +1,38 @@
+# todo - what is the required interface for abstractmpo?
+# support densempo windows?
 """
-A WindowMPS is a finite MPS embedded between an infinite mps to the left, and an inifite mps to the right. 
-The associated hamiltonian has also an infinite part to the left, a finite hamiltonian in the middle, and an infinite part to the right.
+$(TYPEDEF)
 
-Acts simalar as just a finite hamiltonian, but we 'remember' the boundary hamiltonians.
-"""
+The Hamiltonian counterpart of a [`WindowMPS`](@ref): a finite region embedded between an
+infinite environment to the left and to the right.
+It consists of an infinite Hamiltonian to the left, a finite Hamiltonian in the middle, and
+an infinite Hamiltonian to the right.
 
-# todo - what is the required interface for abstractmpo?
-# support densempo windows?
+Acts similar to just a finite Hamiltonian, but we "remember" the boundary Hamiltonians.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+## Constructors
+
+    WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
+
+Construct a `WindowMPOHamiltonian` by carving a finite `interval` out of an infinite
+Hamiltonian `ham`.
+The finite window consists of the sites in `interval`, while the left and right environments
+are copies of `ham` whose unit cells are circshifted so that they line up with the window
+boundaries.
+"""
 struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
+    "Hamiltonian acting on the infinite environment to the left of the window"
     left_ham::InfiniteMPOHamiltonian{O}
+    "Hamiltonian acting on the finite window"
     finite_ham::FiniteMPOHamiltonian{O}
+    "Hamiltonian acting on the infinite environment to the right of the window"
     right_ham::InfiniteMPOHamiltonian{O}
 end
 
-#utility constructor
 function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
 
     # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
@@ -30,8 +49,56 @@ Base.parent(h::WindowMPOHamiltonian) = h.finite_ham
 Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham))
 
 # some basic linalg
-for fun in (:(Base.:+), :(Base.:-), :(Base.:*))
-    @eval $fun(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian) = WindowMPOHamiltonian($fun(a.left_ham, b.left_ham), $fun(a.finite_ham, b.finite_ham), $fun(a.right_ham, b.right_ham))
+# NOTE: `+` cannot be delegated to the regular `FiniteMPOHamiltonian` addition: the finite
+# window carries the full (non-trivial) Jordan virtual spaces at its boundaries, so the two
+# summands have to be block-diagonalized at every site -- including the edges -- exactly like
+# for an `InfiniteMPOHamiltonian`. `-` reuses `+` through the `AbstractMPO` fallback
+# (`a - b == a + (-b)`), and scaling is space-preserving so it works out of the box.
+function Base.:+(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
+    return WindowMPOHamiltonian(
+        a.left_ham + b.left_ham,
+        _add_finite_window(a.finite_ham, b.finite_ham),
+        a.right_ham + b.right_ham
+    )
+end
+function Base.:*(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
+    return WindowMPOHamiltonian(
+        a.left_ham * b.left_ham, a.finite_ham * b.finite_ham, a.right_ham * b.right_ham
+    )
+end
+
+# Scaling a Jordan Hamiltonian scales every path exactly once; since each path starts in
+# exactly one of the three parts, scaling each part by `α` scales the whole operator by `α`.
+# This also powers unary `-`, `*` and `/` through the `AbstractMPO` fallbacks.
+function VectorInterface.scale(H::WindowMPOHamiltonian, α::Number)
+    return WindowMPOHamiltonian(
+        scale(H.left_ham, α), scale(H.finite_ham, α), scale(H.right_ham, α)
+    )
+end
+
+# Block-diagonal addition of two finite Jordan Hamiltonians that have non-trivial virtual
+# spaces on their boundaries (as produced by slicing an infinite Hamiltonian into a window).
+# Contrary to `Base.:+(::FiniteMPOHamiltonian, ::FiniteMPOHamiltonian)` the boundary spaces
+# are grown as well, mirroring the `InfiniteMPOHamiltonian` implementation.
+function _add_finite_window(
+        H₁::FiniteMPOHamiltonian{O}, H₂::FiniteMPOHamiltonian{O}
+    ) where {O <: JordanMPOTensor}
+    N = check_length(H₁, H₂)
+    H = similar(parent(H₁))
+    Vtriv = leftunitspace(first(physicalspace(H₁)))
+    for i in 1:N
+        A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
+        B = cat(H₁[i].B, H₂[i].B; dims = 1)
+        C = cat(H₁[i].C, H₂[i].C; dims = 3)
+        D = H₁[i].D + H₂[i].D
+
+        Vleft = ⊞(Vtriv, left_virtualspace(A), Vtriv)
+        Vright = ⊞(Vtriv, right_virtualspace(A), Vtriv)
+        V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright
+
+        H[i] = JordanMPOTensor(V, A, B, C, D)
+    end
+    return FiniteMPOHamiltonian(H)
 end
 
 TensorKit.dot(
diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl
index 7c5d48327..fb2396da6 100644
--- a/src/states/windowmps.jl
+++ b/src/states/windowmps.jl
@@ -106,7 +106,14 @@ function WindowMPS(N::Int, V::VectorSpace, args...; kwargs...)
     return WindowMPS(fill(V, N), args...; kwargs...)
 end
 
+"""
+    WindowMPS(ψ::InfiniteMPS, L::Int)
+    WindowMPS(ψ::InfiniteMPS, interval::UnitRange)
 
+Construct a [`WindowMPS`](@ref) from an infinite MPS `ψ` by promoting the sites in `interval`
+(or `1:L`) to the finite window while keeping `ψ` as the left and right infinite environments.
+The environment unit cells are circshifted so that they line up with the window boundaries.
+"""
 WindowMPS(ψ::InfiniteMPS, L::Int) = WindowMPS(ψ, 1:L)
 
 function WindowMPS(ψ::InfiniteMPS{A, B}, interval::UnitRange) where {A, B}
diff --git a/test/operators/windowhamiltonian.jl b/test/operators/windowhamiltonian.jl
new file mode 100644
index 000000000..8371ff7f0
--- /dev/null
+++ b/test/operators/windowhamiltonian.jl
@@ -0,0 +1,47 @@
+println("
+---------------------------------
+|   WindowMPOHamiltonian tests   |
+---------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+
+@testset "WindowMPOHamiltonian" begin
+    # a uniform state, represented with a two-site unit cell, so that carving out a window
+    # at different offsets exercises the (non-trivial) circshift alignment of the infinite
+    # environments while the physics stays translationally invariant.
+    H1 = force_planar(transverse_field_ising(; g = 1.5))
+    gs1, = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^12]), H1, VUMPS(; verbosity = 0))
+    gs = repeat(gs1, 2)
+    H = repeat(H1, 2)
+
+    L = 6
+
+    @testset "interval offset (circshift)" begin
+        # the window energy must be independent of the offset for a uniform state
+        energies = map((1:L, 2:(L + 1), 3:(L + 2))) do interval
+            ψ = WindowMPS(gs, interval)
+            Hw = WindowMPOHamiltonian(H, interval)
+            @test length(Hw) == length(interval)
+            return expectation_value(ψ, Hw)
+        end
+        @test energies[1] ≈ energies[2] atol = 1.0e-8
+        @test energies[1] ≈ energies[3] atol = 1.0e-8
+    end
+
+    @testset "linear algebra" begin
+        ψ = WindowMPS(gs, 1:L)
+        Hw = WindowMPOHamiltonian(H, 1:L)
+        e = expectation_value(ψ, Hw)
+
+        @test expectation_value(ψ, Hw + Hw) ≈ 2 * e atol = 1.0e-8
+        @test expectation_value(ψ, Hw - Hw) ≈ 0 atol = 1.0e-8
+        @test expectation_value(ψ, 3 * Hw) ≈ 3 * e atol = 1.0e-8
+        @test expectation_value(ψ, Hw * 3) ≈ 3 * e atol = 1.0e-8
+        @test expectation_value(ψ, -Hw) ≈ -e atol = 1.0e-8
+    end
+end

From 5f018ad64c3b689c425d53b780c297defbe3f665 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 1 Jul 2026 10:15:50 -0400
Subject: [PATCH 12/13] Test WindowMPOHamiltonian addition with distinct
 summands

Add two Hamiltonians with different Jordan bond dimensions (Ising ZZ string
vs XY hopping) so the boundary/bulk virtual spaces of the summands differ and
must be block-diagonalized consistently, checking expectation-value additivity
in both orders.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 test/operators/windowhamiltonian.jl | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/test/operators/windowhamiltonian.jl b/test/operators/windowhamiltonian.jl
index 8371ff7f0..46d504e2a 100644
--- a/test/operators/windowhamiltonian.jl
+++ b/test/operators/windowhamiltonian.jl
@@ -44,4 +44,20 @@ using TensorKit: ℙ
         @test expectation_value(ψ, Hw * 3) ≈ 3 * e atol = 1.0e-8
         @test expectation_value(ψ, -Hw) ≈ -e atol = 1.0e-8
     end
+
+    @testset "distinct summands" begin
+        # add two genuinely different Hamiltonians: the ZZ Ising string has a single bulk
+        # level while the XY hopping needs two, so the boundary/bulk virtual spaces of the
+        # two summands differ and have to be block-diagonalized consistently.
+        ψ = WindowMPS(gs1, 1:L)
+        Ha = WindowMPOHamiltonian(H1, 1:L)
+        Hb = WindowMPOHamiltonian(force_planar(XY_model(; g = 0.7)), 1:L)
+        @test right_virtualspace(Ha.finite_ham, 3) != right_virtualspace(Hb.finite_ham, 3)
+
+        ea = expectation_value(ψ, Ha)
+        eb = expectation_value(ψ, Hb)
+        @test expectation_value(ψ, Ha + Hb) ≈ ea + eb atol = 1.0e-8
+        @test expectation_value(ψ, Hb + Ha) ≈ ea + eb atol = 1.0e-8
+        @test expectation_value(ψ, Ha - Hb) ≈ ea - eb atol = 1.0e-8
+    end
 end

From ad7ba761e6d9a79f6e38a367161b3c5ea569e6fb Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 1 Jul 2026 10:29:26 -0400
Subject: [PATCH 13/13] Inline _add_finite_window into WindowMPOHamiltonian +
 and trim comments

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/operators/windowhamiltonian.jl | 61 ++++++++++--------------------
 1 file changed, 20 insertions(+), 41 deletions(-)

diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index 15ecd552d..f4ace7c45 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -34,13 +34,10 @@ struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
 end
 
 function WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
-
-    # to make sure the interval corresponds with finite_ham, it is important that the unitcell of the left/right hamiltonians is circshifted correctly
     left_edge = (interval.start - 1) % length(ham)
     left_ham = circshift(ham, -left_edge)
     right_edge = (interval.stop + 1) % length(ham)
     right_ham = circshift(ham, -right_edge + 1)
-
     finite_ham = FiniteMPOHamiltonian([ham[i] for i in interval])
     return WindowMPOHamiltonian(left_ham, finite_ham, right_ham)
 end
@@ -48,17 +45,28 @@ end
 Base.parent(h::WindowMPOHamiltonian) = h.finite_ham
 Base.copy(h::WindowMPOHamiltonian) = WindowMPOHamiltonian(copy(h.left_ham), copy(h.finite_ham), copy(h.right_ham))
 
-# some basic linalg
-# NOTE: `+` cannot be delegated to the regular `FiniteMPOHamiltonian` addition: the finite
-# window carries the full (non-trivial) Jordan virtual spaces at its boundaries, so the two
-# summands have to be block-diagonalized at every site -- including the edges -- exactly like
-# for an `InfiniteMPOHamiltonian`. `-` reuses `+` through the `AbstractMPO` fallback
-# (`a - b == a + (-b)`), and scaling is space-preserving so it works out of the box.
 function Base.:+(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
+    # the finite window carries full Jordan virtual spaces at its boundaries, so it has to be
+    # block-diagonalized at every site (like an InfiniteMPOHamiltonian) rather than through
+    # the regular FiniteMPOHamiltonian addition
+    Ha, Hb = a.finite_ham, b.finite_ham
+    N = check_length(Ha, Hb)
+    finite_ham = similar(parent(Ha))
+    Vtriv = leftunitspace(first(physicalspace(Ha)))
+    for i in 1:N
+        A = cat(Ha[i].A, Hb[i].A; dims = (1, 4))
+        B = cat(Ha[i].B, Hb[i].B; dims = 1)
+        C = cat(Ha[i].C, Hb[i].C; dims = 3)
+        D = Ha[i].D + Hb[i].D
+
+        Vleft = ⊞(Vtriv, left_virtualspace(A), Vtriv)
+        Vright = ⊞(Vtriv, right_virtualspace(A), Vtriv)
+        V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright
+
+        finite_ham[i] = JordanMPOTensor(V, A, B, C, D)
+    end
     return WindowMPOHamiltonian(
-        a.left_ham + b.left_ham,
-        _add_finite_window(a.finite_ham, b.finite_ham),
-        a.right_ham + b.right_ham
+        a.left_ham + b.left_ham, FiniteMPOHamiltonian(finite_ham), a.right_ham + b.right_ham
     )
 end
 function Base.:*(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
@@ -66,41 +74,12 @@ function Base.:*(a::WindowMPOHamiltonian, b::WindowMPOHamiltonian)
         a.left_ham * b.left_ham, a.finite_ham * b.finite_ham, a.right_ham * b.right_ham
     )
 end
-
-# Scaling a Jordan Hamiltonian scales every path exactly once; since each path starts in
-# exactly one of the three parts, scaling each part by `α` scales the whole operator by `α`.
-# This also powers unary `-`, `*` and `/` through the `AbstractMPO` fallbacks.
 function VectorInterface.scale(H::WindowMPOHamiltonian, α::Number)
     return WindowMPOHamiltonian(
         scale(H.left_ham, α), scale(H.finite_ham, α), scale(H.right_ham, α)
     )
 end
 
-# Block-diagonal addition of two finite Jordan Hamiltonians that have non-trivial virtual
-# spaces on their boundaries (as produced by slicing an infinite Hamiltonian into a window).
-# Contrary to `Base.:+(::FiniteMPOHamiltonian, ::FiniteMPOHamiltonian)` the boundary spaces
-# are grown as well, mirroring the `InfiniteMPOHamiltonian` implementation.
-function _add_finite_window(
-        H₁::FiniteMPOHamiltonian{O}, H₂::FiniteMPOHamiltonian{O}
-    ) where {O <: JordanMPOTensor}
-    N = check_length(H₁, H₂)
-    H = similar(parent(H₁))
-    Vtriv = leftunitspace(first(physicalspace(H₁)))
-    for i in 1:N
-        A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
-        B = cat(H₁[i].B, H₂[i].B; dims = 1)
-        C = cat(H₁[i].C, H₂[i].C; dims = 3)
-        D = H₁[i].D + H₂[i].D
-
-        Vleft = ⊞(Vtriv, left_virtualspace(A), Vtriv)
-        Vright = ⊞(Vtriv, right_virtualspace(A), Vtriv)
-        V = Vleft ⊗ physicalspace(A) ← physicalspace(A) ⊗ Vright
-
-        H[i] = JordanMPOTensor(V, A, B, C, D)
-    end
-    return FiniteMPOHamiltonian(H)
-end
-
 TensorKit.dot(
     bra::WindowMPS, H::WindowMPOHamiltonian, ket::WindowMPS = bra,
     envs = environments(bra, H, ket)