diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 0139d60d0..e04679871 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -28,10 +28,13 @@ using LoggingExtras
 import TupleTools
 
 using MPSKit
-using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, ProductTransferMatrix
+using MPSKit:
+    MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor,
+    ProductTransferMatrix, TransferMatrix
 using MPSKit: InfiniteEnvironments
 import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace
 import MPSKit: infinite_temperature_density_matrix
+import MPSKit: fuser
 
 import TensorKitTensors.SpinOperators as SO
 import TensorKitTensors.FermionOperators as FO
@@ -73,6 +76,8 @@ include("operators/infinitepepo.jl")
 include("operators/transfermatrix.jl")
 include("operators/localoperator.jl")
 include("operators/localcircuit.jl")
+include("operators/mpo_observable.jl")
+
 include("operators/lattices/squarelattice.jl")
 include("operators/models.jl")
 
@@ -107,6 +112,12 @@ include("algorithms/contractions/correlator/peps.jl")
 include("algorithms/contractions/correlator/pepo_purified.jl")
 include("algorithms/contractions/correlator/pepo_1layer.jl")
 
+include("algorithms/contractions/mpo_path/pepo_1layer.jl")
+include("algorithms/contractions/window/tools.jl")
+include("algorithms/contractions/window/pepo_1layer.jl")
+include("algorithms/contractions/window/twosite/caching.jl")
+include("algorithms/contractions/window/twosite/pepo_1layer.jl")
+
 include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/projectors.jl")
@@ -139,6 +150,8 @@ include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
 include("algorithms/correlator_adapters.jl")
 include("algorithms/correlators.jl")
+include("algorithms/expval_approx.jl")
+include("algorithms/correlator_approx.jl")
 
 include("algorithms/optimization/fixed_point_differentiation.jl")
 include("algorithms/optimization/peps_optimization.jl")
@@ -156,9 +169,10 @@ export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
 export initialize_random_c4v_env, initialize_singlet_c4v_env
-export LocalOperator, physicalspace
+export LocalOperator, MPOObservable, physicalspace
 export product_peps
-export reduced_densitymatrix, expectation_value, network_value, cost_function
+export reduced_densitymatrix, expectation_value_approx, correlator_approx
+export expectation_value, network_value, cost_function
 export correlator, correlation_length
 export leading_boundary
 export PEPSOptimize, FixedPointGradient, GeomSum, ManualIter
diff --git a/src/algorithms/contractions/mpo_path/pepo_1layer.jl b/src/algorithms/contractions/mpo_path/pepo_1layer.jl
new file mode 100644
index 000000000..7628b0413
--- /dev/null
+++ b/src/algorithms/contractions/mpo_path/pepo_1layer.jl
@@ -0,0 +1,238 @@
+"""
+Check that the physical legs of a first OBC-MPO tensor match the PEPO site's physical space.
+"""
+function _check_pepo_first_physicalspace(A, op)
+    physicalspace(A) == space(op, 1) == space(op, 2)' ||
+        throw(SpaceMismatch("first MPO tensor physical space does not match PEPO site"))
+    return nothing
+end
+
+"""
+Check that the physical legs of a last OBC-MPO tensor match the PEPO site's physical space.
+"""
+function _check_pepo_last_physicalspace(A, op)
+    physicalspace(A) == space(op, 2) == space(op, 3)' ||
+        throw(SpaceMismatch("last MPO tensor physical space does not match PEPO site"))
+    return nothing
+end
+
+"""
+Check that the physical legs of a middle MPO tensor match the PEPO site's physical space.
+"""
+function _check_pepo_middle_physicalspace(A, op)
+    physicalspace(A) == space(op, 2) == space(op, 3)' ||
+        throw(SpaceMismatch("middle MPO tensor physical space does not match PEPO site"))
+    return nothing
+end
+
+"""
+Convert a symbolic cardinal path direction to the corresponding PEPO virtual-leg index.
+"""
+function _mpo_path_direction(direction::Symbol)
+    direction === :north && return NORTH
+    direction === :east && return EAST
+    direction === :south && return SOUTH
+    direction === :west && return WEST
+    throw(ArgumentError("invalid MPO path direction: $direction"))
+end
+
+"""
+Return the tensor-expression label for the PEPO virtual leg in a cardinal direction.
+"""
+function _mpo_path_virtual_label(direction::Symbol)
+    return (:N, :E, :S, :W)[_mpo_path_direction(direction)]
+end
+
+"""
+Canonicalize an incoming MPO fuser so its fused PEPO leg has standard dualness.
+"""
+function _mpo_path_incoming_fuser(F, direction::Int)
+    direction in (NORTH, EAST) && return F
+    direction in (SOUTH, WEST) && return twist(flip(F, 1), 1)
+    throw(ArgumentError("invalid MPO path direction index: $direction"))
+end
+
+"""
+Canonicalize an outgoing MPO fuser so its fused PEPO leg has standard dualness and braiding.
+"""
+function _mpo_path_outgoing_fuser(F, direction::Int)
+    direction in (NORTH, EAST) && return twist(flip(F, 1), 1)
+    direction in (SOUTH, WEST) && return twist(F, 3)
+    throw(ArgumentError("invalid MPO path direction index: $direction"))
+end
+
+"""
+Build the `@tensor` labels from `(direction, suffix)` pairs
+used to fuse MPO virtual strings.
+
+- Each direction is `:north`, `:east`, `:south`, or `:west`.
+- Suffix `:l` marks an incoming MPO bond and `:r` an outgoing one.
+
+Examples:
+
+- `((:east, :r),)` produces `[W S; N Er]`.
+- `((:west, :l), (:north, :r))` produces `[Wl S; Nr E]`.
+"""
+function _mpo_path_result_expr(directions)
+    labels = [:N, :E, :S, :W]
+    for (direction, suffix) in directions
+        index = _mpo_path_direction(direction)
+        labels[index] = Symbol(labels[index], suffix)
+    end
+    return tensorexpr(:t, (labels[WEST], labels[SOUTH]), (labels[NORTH], labels[EAST]))
+end
+
+"""
+Act the first tensor `op` of an OBC-MPO on PEPO tensor `A` and fuse the
+outgoing MPO string with the virtual space of `A` along `direction`.
+"""
+@generated function mpo_path_first(A::PEPOTensor, op, ::Val{direction}) where {direction}
+    direction_index = _mpo_path_direction(direction)
+    virtual_label = _mpo_path_virtual_label(direction)
+    fused_label = Symbol(virtual_label, :r)
+
+    result_e = _mpo_path_result_expr(((direction, :r),))
+    op_e = tensorexpr(:op, :dout, (:din, :r))
+    A_e = tensorexpr(:A′, (:din, :dout), (:N, :E, :S, :W))
+    F_e = tensorexpr(:F, fused_label, (virtual_label, :r))
+    rhs = Expr(:call, :*, op_e, A_e, F_e)
+    contraction = macroexpand(
+        @__MODULE__, :(return @tensoropt $result_e := $rhs)
+    )
+
+    return quote
+        _check_pepo_first_physicalspace(A, op)
+        A′ = twistdual(A, 2)
+        F = _mpo_path_outgoing_fuser(
+            fuser(storagetype(A), domain(A, $direction_index)', space(op, 3)),
+            $direction_index,
+        )
+        $contraction
+    end
+end
+
+"""
+Act the last tensor `op` of an OBC-MPO on PEPO tensor `A` and fuse the
+incoming MPO string with the virtual space of `A` along `direction`.
+"""
+@generated function mpo_path_last(A::PEPOTensor, op, ::Val{direction}) where {direction}
+    direction_index = _mpo_path_direction(direction)
+    virtual_label = _mpo_path_virtual_label(direction)
+    fused_label = Symbol(virtual_label, :l)
+
+    result_e = _mpo_path_result_expr(((direction, :l),))
+    F_e = Expr(:call, :conj, tensorexpr(:F, fused_label, (virtual_label, :l)))
+    op_e = tensorexpr(:op, (:l, :dout), :din)
+    A_e = tensorexpr(:A′, (:din, :dout), (:N, :E, :S, :W))
+    rhs = Expr(:call, :*, F_e, op_e, A_e)
+    contraction = macroexpand(
+        @__MODULE__, :(return @tensoropt $result_e := $rhs)
+    )
+
+    return quote
+        _check_pepo_last_physicalspace(A, op)
+        A′ = twistdual(A, 2)
+        F = _mpo_path_incoming_fuser(
+            fuser(storagetype(A), domain(A, $direction_index), space(op, 1)'),
+            $direction_index,
+        )
+        $contraction
+    end
+end
+
+"""
+Act the middle tensor `op` of an MPO on PEPO tensor `A` and fuse the
+incoming and the outgoing MPO string with the virtual space of `A` along
+`directions = (incoming, outgoing)`.
+"""
+@generated function mpo_path_middle(
+        A::PEPOTensor, op, ::Val{directions}
+    ) where {directions}
+    incoming, outgoing = directions
+    incoming == outgoing &&
+        throw(ArgumentError("MPO path should enter and exit in different directions"))
+
+    incoming_index = _mpo_path_direction(incoming)
+    outgoing_index = _mpo_path_direction(outgoing)
+    incoming_label = _mpo_path_virtual_label(incoming)
+    outgoing_label = _mpo_path_virtual_label(outgoing)
+    fused_incoming_label = Symbol(incoming_label, :l)
+    fused_outgoing_label = Symbol(outgoing_label, :r)
+
+    result_e = _mpo_path_result_expr(((incoming, :l), (outgoing, :r)))
+    Fin_e = Expr(
+        :call, :conj,
+        tensorexpr(:Fin, fused_incoming_label, (incoming_label, :l)),
+    )
+    op_e = tensorexpr(:op, (:l, :dout), (:din, :r))
+    A_e = tensorexpr(:A′, (:din, :dout), (:N, :E, :S, :W))
+    Fout_e = tensorexpr(
+        :Fout, fused_outgoing_label, (outgoing_label, :r)
+    )
+    rhs = Expr(:call, :*, Fin_e, op_e, A_e, Fout_e)
+    contraction = macroexpand(
+        @__MODULE__, :(return @tensoropt $result_e := $rhs)
+    )
+
+    return quote
+        _check_pepo_middle_physicalspace(A, op)
+        A′ = twistdual(A, 2)
+        Fin = _mpo_path_incoming_fuser(
+            fuser(storagetype(A), domain(A, $incoming_index), space(op, 1)'),
+            $incoming_index,
+        )
+        Fout = _mpo_path_outgoing_fuser(
+            fuser(storagetype(A), domain(A, $outgoing_index)', space(op, 4)),
+            $outgoing_index,
+        )
+        $contraction
+    end
+end
+
+"""
+Route an MPO virtual string with `stringspace` through a PEPO tensor `A`
+along `directions = (incoming, outgoing)`.
+"""
+@generated function mpo_path_string(
+        A::PEPOTensor, stringspace::ElementarySpace, ::Val{directions}
+    ) where {directions}
+    incoming, outgoing = directions
+    incoming == outgoing &&
+        throw(ArgumentError("MPO path should enter and exit in different directions"))
+
+    incoming_index = _mpo_path_direction(incoming)
+    outgoing_index = _mpo_path_direction(outgoing)
+    incoming_label = _mpo_path_virtual_label(incoming)
+    outgoing_label = _mpo_path_virtual_label(outgoing)
+    fused_incoming_label = Symbol(incoming_label, :l)
+    fused_outgoing_label = Symbol(outgoing_label, :r)
+
+    result_e = _mpo_path_result_expr(((incoming, :l), (outgoing, :r)))
+    Fin_e = Expr(
+        :call, :conj,
+        tensorexpr(:Fin, fused_incoming_label, (incoming_label, :l)),
+    )
+    O_e = tensorexpr(:O, (:W, :S), (:N, :E))
+    I_e = tensorexpr(:I, :l, :r)
+    Fout_e = tensorexpr(
+        :Fout, fused_outgoing_label, (outgoing_label, :r)
+    )
+    rhs = Expr(:call, :*, Fin_e, O_e, I_e, Fout_e)
+    contraction = macroexpand(
+        @__MODULE__, :(return @tensoropt $result_e := $rhs)
+    )
+
+    return quote
+        O = trace_physicalspaces(A)
+        I = id(storagetype(A), stringspace)
+        Fin = _mpo_path_incoming_fuser(
+            fuser(storagetype(A), domain(A, $incoming_index), stringspace'),
+            $incoming_index,
+        )
+        Fout = _mpo_path_outgoing_fuser(
+            fuser(storagetype(A), domain(A, $outgoing_index)', stringspace'),
+            $outgoing_index,
+        )
+        $contraction
+    end
+end
diff --git a/src/algorithms/contractions/window/pepo_1layer.jl b/src/algorithms/contractions/window/pepo_1layer.jl
new file mode 100644
index 000000000..8531553c7
--- /dev/null
+++ b/src/algorithms/contractions/window/pepo_1layer.jl
@@ -0,0 +1,189 @@
+# Approximate finite-window contractions for single-layer PEPO networks.
+
+"""
+Validate that the network is a single-layer PEPO and that the sweep direction is supported.
+"""
+function _check_window_inputs(ρ::InfinitePEPO, direction::Symbol)
+    size(ρ, 3) == 1 || throw(DimensionMismatch("only single-layer PEPO contractions are supported"))
+    direction in (:auto, :rows, :columns) ||
+        throw(ArgumentError("invalid sweep direction: $direction"))
+    return nothing
+end
+
+"""
+Return a PEPO and CTMRG environment with standard virtual-space dualness without mutating the inputs.
+"""
+function standardize_dualness(ρ::InfinitePEPO, env::CTMRGEnv)
+    isdual_easts, isdual_norths = _check_virtual_dualness(ρ)
+    all(isdual_easts) && all(isdual_norths) && return ρ, env
+
+    nrows, ncols = size(ρ, 1), size(ρ, 2)
+    tensors = map(CartesianIndices(unitcell(ρ))) do site
+        row, col, layer = Tuple(site)
+        directions = Int[]
+        !isdual_norths[row, col, layer] && push!(directions, NORTH)
+        !isdual_easts[row, col, layer] && push!(directions, EAST)
+        !isdual_norths[_next(row, nrows), col, layer] && push!(directions, SOUTH)
+        !isdual_easts[row, _prev(col, ncols), layer] && push!(directions, WEST)
+        A = unitcell(ρ)[site]
+        return isempty(directions) ? A : flip_virtualspace(A, directions)
+    end
+    ρ′ = InfinitePEPO(tensors)
+
+    edges = map(CartesianIndices(env.edges)) do index
+        direction, row, col = Tuple(index)
+        should_flip = if direction == NORTH
+            !isdual_norths[_next(row, nrows), col, 1]
+        elseif direction == EAST
+            !isdual_easts[row, _prev(col, ncols), 1]
+        elseif direction == SOUTH
+            !isdual_norths[row, col, 1]
+        else
+            !isdual_easts[row, col, 1]
+        end
+        E = env.edges[index]
+        return should_flip ? flip(E, 2) : E
+    end
+    env′ = CTMRGEnv(copy(env.corners), edges)
+    return ρ′, env′
+end
+
+"""
+Contract an MPO observable in its enclosing window, rotating column sweeps into row sweeps.
+"""
+function _expectation_value_approx(
+        ρ::InfinitePEPO, observable::MPOObservable, env::CTMRGEnv,
+        alg::WindowApprox, direction::Symbol,
+    )
+    _check_window_inputs(ρ, direction)
+    rowrange, colrange = _window_ranges(observable)
+    sweep = direction === :auto ? (length(colrange) > length(rowrange) ? :rows : :columns) : direction
+    if sweep === :rows
+        return _expectation_value_approx_rows(
+            ρ, observable, env, rowrange, colrange, alg
+        )
+    else
+        unitcell = size(ρ)[1:2]
+        sites = siterotl90.(observable.sites, Ref(unitcell))
+        path = siterotl90.(observable.path, Ref(unitcell))
+        rotated_observable = MPOObservable(sites, path, observable.mpo)
+        rotated_rowrange, rotated_colrange = _window_ranges(rotated_observable)
+        return _expectation_value_approx_rows(
+            rotl90(ρ), rotated_observable, rotl90(env),
+            rotated_rowrange, rotated_colrange, alg
+        )
+    end
+end
+
+"""
+Build a local tensor for row MPOs without observables by tracing the PEPO physical legs.
+"""
+function _window_site_tensor(
+        ρ::InfinitePEPO, ::Nothing, row::Int, col::Int,
+    )
+    return trace_physicalspaces(ρ[row, col, 1])
+end
+
+"""
+Build the local row-MPO tensor at one window site, inserting the observable tensor or routed
+string when the site lies on the MPO path and tracing the PEPO physical legs otherwise.
+"""
+function _window_site_tensor(
+        ρ::InfinitePEPO, observable::MPOObservable, row::Int, col::Int,
+    )
+    A = ρ[row, col, 1]
+    site = CartesianIndex(row, col)
+    path_index = findfirst(==(site), observable.path)
+    isnothing(path_index) && return trace_physicalspaces(A)
+
+    mpo_index = findfirst(==(site), observable.sites)
+    # sites with string passing by (cannot be first/last site)
+    if isnothing(mpo_index)
+        incoming = _step_direction(observable.path[path_index], observable.path[path_index - 1])
+        outgoing = _step_direction(observable.path[path_index], observable.path[path_index + 1])
+        next_mpo_index = count(in(observable.sites), @view observable.path[1:path_index]) + 1
+        stringspace = space(observable.mpo[next_mpo_index], 1)
+        return mpo_path_string(A, stringspace, Val((incoming, outgoing)))
+    end
+
+    # sites acted on by the MPO
+    op = observable.mpo[mpo_index]
+    if mpo_index == 1
+        direction = _step_direction(observable.path[1], observable.path[2])
+        return mpo_path_first(A, op, Val(direction))
+    elseif mpo_index == length(observable.mpo)
+        direction = _step_direction(observable.path[end], observable.path[end - 1])
+        return mpo_path_last(A, op, Val(direction))
+    else
+        incoming = _step_direction(observable.path[path_index], observable.path[path_index - 1])
+        outgoing = _step_direction(observable.path[path_index], observable.path[path_index + 1])
+        return mpo_path_middle(A, op, Val((incoming, outgoing)))
+    end
+end
+
+"""
+Contract and normalize an MPO observable using row-oriented window boundary contractions.
+"""
+function _expectation_value_approx_rows(
+        ρ::InfinitePEPO, observable::MPOObservable, env::CTMRGEnv,
+        rowrange::UnitRange{Int}, colrange::UnitRange{Int}, alg::WindowApprox,
+    )
+    ρ, env = standardize_dualness(ρ, env)
+    numerator = _contract_window_rows(ρ, observable, env, rowrange, colrange, alg)
+    norm = _contract_window_rows(ρ, nothing, env, rowrange, colrange, alg)
+    return numerator / norm
+end
+
+"""
+Contract a complete PEPO window row by row from north to south,
+optionally inserting an MPO observable.
+"""
+function _contract_window_rows(
+        ρ::InfinitePEPO, observable::Union{Nothing, MPOObservable},
+        env::CTMRGEnv, rowrange::UnitRange{Int}, colrange::UnitRange{Int},
+        alg::WindowApprox,
+    )
+    ψ = _north_boundary_mps(env, first(rowrange), colrange)
+    for row in rowrange
+        W = _row_mpo(ρ, observable, env, row, colrange)
+        ψ = _approximate_window_step(W, ψ, alg)
+    end
+    south = _south_boundary_mps(env, last(rowrange), colrange)
+    return dot(south, ψ)
+end
+
+"""
+Build one finite row MPO from west/east CTMRG edges and the PEPO tensors inside the window.
+
+Convention of west, east CTM edges and the PF tensors:
+```
+    [1 2; 3]    [1 2; 3 4]     [1 2; 3]
+    3               3               1
+    ↓               ↓               ↑
+    E₄-←-2      1-←-O-←-4       2-←-C₂
+    ↓               ↓               ↑
+    1               2               3
+```
+Legs 1, 3 need to be flipped to match standard MPS convention
+"""
+function _row_mpo(
+        ρ::InfinitePEPO, observable::Union{Nothing, MPOObservable},
+        env::CTMRGEnv, row::Int, colrange::UnitRange{Int},
+    )
+    cmin, cmax = first(colrange), last(colrange)
+    W = repartition(edge(env, WEST, row, cmin - 1), 1, 2)
+    tensors = [insertleftunit(W, 1)]
+    append!(
+        tensors,
+        (
+            _window_site_tensor(ρ, observable, row, col)
+                for col in colrange
+        ),
+    )
+    E = permute(
+        flip(edge(env, EAST, row, cmax + 1), (1, 3)),
+        ((2, 3), (1,))
+    )
+    push!(tensors, insertrightunit(E, 3))
+    return FiniteMPO(tensors)
+end
diff --git a/src/algorithms/contractions/window/tools.jl b/src/algorithms/contractions/window/tools.jl
new file mode 100644
index 000000000..72ab61b8a
--- /dev/null
+++ b/src/algorithms/contractions/window/tools.jl
@@ -0,0 +1,96 @@
+"""
+Bundle the zip-up contraction and optional DMRG refinement
+algorithms used after each finite window MPO-MPS contraction.
+"""
+struct WindowApprox{Z, D}
+    zipup::Z
+    dmrg::D
+end
+
+# TODO: generalize the following to multi-layer networks
+
+"""
+Apply a finite MPO to a finite MPS with zip-up truncation and optional DMRG refinement.
+"""
+function _approximate_window_step(W::FiniteMPO, ψ::FiniteMPS, alg::WindowApprox)
+    ψ′, = approximate((W, ψ), alg.zipup)
+    isnothing(alg.dmrg) && return ψ′
+    ψ′, = approximate(ψ′, (W, ψ), alg.dmrg)
+    return ψ′
+end
+
+"""
+Convert a south-boundary MPS tensor into the stored bra representation expected by `dot` using a planar repartition.
+"""
+function _bra_mps_tensor(A::MPSTensor)
+    return repartition(A', 2, 1; copy = true)
+end
+
+"""
+Construct the planar adjoint of a finite MPO while restoring MPSKit's local MPO leg partition.
+"""
+function _adjoint_mpo(W::FiniteMPO)
+    return FiniteMPO(map(A -> transpose(A', ((3, 1), (4, 2)); copy = true), parent(W)))
+end
+
+"""
+Build the finite MPS representing the north CTMRG boundary of a window.
+
+Convention of CTM tensors on the north boundary is
+```
+    [1; 2]      [1 2; 3]        [1; 2]
+    C₁-←-2      1-←-E₁-←-3      1-←-C₂
+    ↓               ↓               ↓
+    1               2               2
+```
+Leg 2 of C₂ needs to be flipped to match standard MPS convention.
+"""
+function _north_boundary_mps(
+        env::CTMRGEnv, row::Int, colrange::UnitRange{Int},
+    )
+    r = row - 1
+    cmin, cmax = first(colrange), last(colrange)
+    Cwest = insertleftunit(corner(env, NORTHWEST, r, cmin - 1), 1)
+    tensors = [Cwest]
+    append!(tensors, (edge(env, NORTH, r, col) for col in colrange))
+    Ceast = repartition(
+        flip(corner(env, NORTHEAST, r, cmax + 1), 2), 2, 0
+    )
+    push!(tensors, insertleftunit(Ceast, 3))
+    return FiniteMPS(tensors)
+end
+
+"""
+Build the finite MPS representing the adjointed south CTMRG boundary of a window.
+
+Convention of CTM tensors on the south boundary is
+```
+    [1; 2]      [1 2; 3]        [1; 2]
+    2               2               1
+    ↓               ↓               ↑
+    C₄-→-1      3-→-E₃-→-1      2-→-C₃
+```
+Leg 1 of C₃ needs to be flipped to match standard MPS convention.
+Then, their adjoints are
+```
+    [1; 2]      [1; 2 3]        [1; 2]
+    C̄₄-←-2      1-←-Ē₃-←-2      1-←-C̄₃
+    ↓               ↓               ↓
+    1               3               2
+```
+The edge tensors then need a further repartition of indices.
+"""
+function _south_boundary_mps(
+        env::CTMRGEnv, row::Int, colrange::UnitRange{Int},
+    )
+    r = row + 1
+    cmin, cmax = first(colrange), last(colrange)
+    Cwest = insertleftunit(corner(env, SOUTHWEST, r, cmin - 1)', 1)
+    tensors = [Cwest]
+    append!(tensors, (_bra_mps_tensor(edge(env, SOUTH, r, col)) for col in colrange))
+    Ceast = repartition(
+        flip(corner(env, SOUTHEAST, r, cmax + 1), 1)', 2, 0
+    )
+    push!(tensors, insertleftunit(Ceast, 3))
+    return FiniteMPS(tensors)
+end
diff --git a/src/algorithms/contractions/window/twosite/caching.jl b/src/algorithms/contractions/window/twosite/caching.jl
new file mode 100644
index 000000000..28f1e40a5
--- /dev/null
+++ b/src/algorithms/contractions/window/twosite/caching.jl
@@ -0,0 +1,124 @@
+"""
+Put 2-site bonds in groups to reuse partial contractions in `correlator_approx`.
+
+For every bond `(first_site, second_site)` in `bonds`:
+
+- After ordering the two sites, `source`/`target` is the first/second site.
+- `swapped` is `false` when `source == first_site`, and `true` otherwise.
+- Bonds with the same `source` and `swapped` are grouped together.
+- In each group, the inner dict records each bond's position in `bonds`.
+
+For example, the ordered bonds
+
+```julia
+CI = CartesianIndex
+bonds = [
+    (CI(1, 1), CI(1, 3)),
+    (CI(1, 1), CI(2, 2)),
+    (CI(1, 2), CI(2, 2)), # another source site
+    (CI(2, 2), CI(1, 1)), # reversed site order
+]
+```
+
+are grouped as
+
+```julia
+(CI(1, 1), false) => Dict(
+    CI(1, 3) => 1,
+    CI(2, 2) => 2,
+)
+(CI(1, 1), true) => Dict(CI(2, 2) => 4)
+(CI(1, 2), false) => Dict(CI(2, 2) => 3)
+```
+"""
+function _twosite_source_groups(
+        bonds::Vector{NTuple{2, CartesianIndex{2}}},
+    )
+    groups = Dict{
+        Tuple{CartesianIndex{2}, Bool},
+        Dict{CartesianIndex{2}, Int},
+    }()
+    for (i, (first_site, second_site)) in enumerate(bonds)
+        swapped = !issorted((first_site, second_site); by = Tuple)
+        source, target = swapped ? (second_site, first_site) : (first_site, second_site)
+        targets = get!(Dict{CartesianIndex{2}, Int}, groups, (source, swapped))
+        targets[target] = i
+    end
+    return groups
+end
+
+"""
+Group the targets associated with one source by their row coordinate. The returned outer
+dictionary maps each target row to a dictionary whose entries retain the original
+`target => result_position` mapping.
+
+This lookup lets the source contraction close all targets in the current row together while
+propagating a single open MPO string between rows.
+"""
+function _twosite_targets_by_row(targets::Dict{CartesianIndex{2}, Int})
+    targets_by_row = Dict{Int, Dict{CartesianIndex{2}, Int}}()
+    for (target, position) in targets
+        row_targets = get!(Dict{CartesianIndex{2}, Int}, targets_by_row, target[1])
+        row_targets[target] = position
+    end
+    return targets_by_row
+end
+
+"""
+Cache the row MPOs without observables and boundary contractions shared by measurements in one window.
+
+The fields contain:
+
+- `rowrange` and `colrange`: the coordinate ranges defining the window.
+- `row_mpos`: the row MPOs without observables, one per row and including the west and east CTMRG edges.
+- `north_prefixes`: `nrows + 1` north boundary MPSs.
+  Entry `k` is above row `k` in the window.
+- `south_suffixes`: `nrows + 1` adjointed south boundary MPSs.
+  Entry `k + 1` is below row `k` in the window.
+- `norm`: the approximate contraction of the window with no observable inserted.
+"""
+struct WindowRowCache{M <: FiniteMPO, N <: FiniteMPS, S <: FiniteMPS, T <: Number}
+    rowrange::UnitRange{Int}
+    colrange::UnitRange{Int}
+    row_mpos::Dict{Int, M}
+    north_prefixes::Vector{N}
+    south_suffixes::Vector{S}
+    norm::T
+end
+
+"""
+Precompute the row MPOs without observables and north/south boundary contractions reused when measuring many two-site bonds in one window.
+"""
+function _window_row_cache(
+        ρ::InfinitePEPO, env::CTMRGEnv,
+        rowrange::UnitRange{Int}, colrange::UnitRange{Int}, alg::WindowApprox,
+    )::WindowRowCache
+    row_mpos = Dict(
+        row => _row_mpo(ρ, nothing, env, row, colrange)
+            for row in rowrange
+    )
+    nrows = length(rowrange)
+    north = _north_boundary_mps(env, first(rowrange), colrange)
+    south = _south_boundary_mps(env, last(rowrange), colrange)
+
+    north_prefixes = Vector{typeof(north)}(undef, nrows + 1)
+    north_prefixes[1] = north
+    for (k, row) in enumerate(rowrange)
+        north_prefixes[k + 1] = _approximate_window_step(
+            row_mpos[row], north_prefixes[k], alg
+        )
+    end
+
+    south_suffixes = Vector{typeof(south)}(undef, nrows + 1)
+    south_suffixes[end] = south
+    for (k, row) in Iterators.reverse(enumerate(rowrange))
+        W = _adjoint_mpo(row_mpos[row])
+        south_suffixes[k] = _approximate_window_step(
+            W, south_suffixes[k + 1], alg
+        )
+    end
+    norm = dot(south, north_prefixes[end])
+    return WindowRowCache(
+        rowrange, colrange, row_mpos, north_prefixes, south_suffixes, norm
+    )
+end
diff --git a/src/algorithms/contractions/window/twosite/pepo_1layer.jl b/src/algorithms/contractions/window/twosite/pepo_1layer.jl
new file mode 100644
index 000000000..c380e3090
--- /dev/null
+++ b/src/algorithms/contractions/window/twosite/pepo_1layer.jl
@@ -0,0 +1,251 @@
+# Source-cached dense two-site contractions for single-layer PEPO networks.
+
+"""
+Validate a dense two-site measurement, choose its sweep orientation, and dispatch to the
+row-oriented source-cached contraction.
+"""
+function _correlator_approx(
+        ρ::InfinitePEPO, op::AbstractTensorMap,
+        bonds::Vector{NTuple{2, CartesianIndex{2}}},
+        env::CTMRGEnv, alg::WindowApprox, direction::Symbol,
+    )
+    _check_window_inputs(ρ, direction)
+    numout(op) == numin(op) == 2 ||
+        throw(ArgumentError("correlator_approx requires a two-site operator"))
+    for bond in bonds, (i, site) in enumerate(bond)
+        V = physicalspace(ρ, Tuple(site)...)
+        V == codomain(op)[i] == domain(op)[i] ||
+            throw(SpaceMismatch("operator physical space does not match PEPO site $site"))
+    end
+
+    rowrange, colrange = _window_ranges(bonds)
+    sweep = direction === :auto ? (length(colrange) > length(rowrange) ? :rows : :columns) : direction
+    if sweep === :rows
+        return _correlator_approx_rows(
+            ρ, op, bonds, env, rowrange, colrange, alg
+        )
+    end
+    # rotate column-wise contraction to reuse row-wise code
+    unitcell = size(ρ)[1:2]
+    rotated_bonds = map(bonds) do bond
+        return (siterotl90(bond[1], unitcell), siterotl90(bond[2], unitcell))
+    end
+    rotated_rowrange, rotated_colrange = _window_ranges(rotated_bonds)
+    return _correlator_approx_rows(
+        rotl90(ρ), op, rotated_bonds, rotl90(env),
+        rotated_rowrange, rotated_colrange, alg
+    )
+end
+
+"""
+Measure all ordered bonds in one row-oriented window using shared row MPOs without observables, shared boundaries, and one exactly decomposed MPO for each operator-leg ordering.
+"""
+function _correlator_approx_rows(
+        ρ::InfinitePEPO, op::AbstractTensorMap,
+        bonds::Vector{NTuple{2, CartesianIndex{2}}}, env::CTMRGEnv,
+        rowrange::UnitRange{Int}, colrange::UnitRange{Int}, alg::WindowApprox,
+    )
+    ρ, env = standardize_dualness(ρ, env)
+    cache = _window_row_cache(ρ, env, rowrange, colrange, alg)
+    groups = _twosite_source_groups(bonds)
+    mpo = gate_to_mpo(op; trunc = notrunc())
+    swapped_mpo = if any(key[2] for key in keys(groups))
+        swapped_op = permute(op, ((2, 1), (4, 3)))
+        gate_to_mpo(swapped_op; trunc = notrunc())
+    end
+
+    T = promote_type(scalartype(op), typeof(cache.norm))
+    numerators = zeros(T, length(bonds))
+    for ((source, swapped), targets) in groups
+        source_mpo = swapped ? something(swapped_mpo) : mpo
+        _contract_twosite_source!(
+            numerators, ρ, source_mpo, source, targets, env, cache, alg,
+        )
+    end
+    return numerators ./ cache.norm
+end
+
+"""
+Contract the correlator numerator for all targets associated with the same source
+and one ordering of the dense operator, writing each result into `numerators`.
+"""
+function _contract_twosite_source!(
+        numerators::Vector{<:Number}, ρ::InfinitePEPO,
+        mpo::AbstractVector{<:AbstractTensorMap},
+        source::CartesianIndex{2}, targets::Dict{CartesianIndex{2}, Int},
+        env::CTMRGEnv, cache::WindowRowCache, alg::WindowApprox,
+    )
+    # grouping targets by which row they are in
+    targets_by_row = _twosite_targets_by_row(targets)
+    source_idx = source[1] - first(cache.rowrange) + 1
+    north = cache.north_prefixes[source_idx]
+
+    # Close targets in the same row as the source
+    if haskey(targets_by_row, source[1])
+        _contract_twosite_target_row!(
+            numerators, ρ, mpo, source, targets_by_row[source[1]], north, cache
+        )
+    end
+    last_target_row = maximum(keys(targets_by_row))
+    last_target_row == source[1] && return numerators
+
+    # Open the MPO string toward the south for targets in later rows.
+    A = ρ[source[1], source[2], 1]
+    source_tensor = mpo_path_first(A, mpo[1], Val(:south))
+    W = _row_mpo_with_site(ρ, source_tensor, env, source[1], source[2], cache.colrange)
+    north = _approximate_window_step(W, north, alg)
+
+    stringspace = space(mpo[2], 1)
+    for row in (source[1] + 1):last_target_row
+        # Close every target in this row
+        if haskey(targets_by_row, row)
+            _contract_twosite_target_row!(
+                numerators, ρ, mpo, source, targets_by_row[row], north, cache
+            )
+        end
+        row == last_target_row && break
+        # Carry the open string down to the next row
+        A = ρ[row, source[2], 1]
+        string_tensor = mpo_path_string(A, stringspace, Val((:north, :south)))
+        W = _row_mpo_with_site(ρ, string_tensor, env, row, source[2], cache.colrange)
+        north = _approximate_window_step(W, north, alg)
+    end
+    return numerators
+end
+
+"""
+Contract the correlator numerator for all targets in one row with a
+shared open-string north state, writing the results into `numerators`.
+"""
+function _contract_twosite_target_row!(
+        numerators::Vector{<:Number}, ρ::InfinitePEPO,
+        mpo::AbstractVector{<:AbstractTensorMap},
+        source::CartesianIndex{2}, targets::Dict{CartesianIndex{2}, Int},
+        north::FiniteMPS, cache::WindowRowCache,
+    )
+    row = first(keys(targets))[1]
+    row_idx = row - first(cache.rowrange) + 1
+    south = cache.south_suffixes[row_idx + 1]
+    envs = environments(south, cache.row_mpos[row], north)
+    source_site = _window_mps_site(source[2], cache.colrange)
+    stringspace = space(mpo[2], 1)
+
+    # close the target right at the column of the incoming string
+    same_col = get(targets, CartesianIndex(row, source[2]), nothing)
+    if !isnothing(same_col)
+        target_tensor = mpo_path_last(ρ[row, source[2], 1], mpo[2], Val(:north))
+        value = _contract_window_site(envs, north, south, source_site, target_tensor)
+        numerators[same_col] = value
+    end
+
+    # Close targets on the right of the incoming string from left to right
+    right_targets = [target for target in keys(targets) if target[2] > source[2]]
+    if !isempty(right_targets)
+        sort!(right_targets; by = x -> x[2])
+        A = ρ[row, source[2], 1]
+        source_tensor = if row == source[1]
+            mpo_path_first(A, mpo[1], Val(:east))
+        else
+            mpo_path_string(A, stringspace, Val((:north, :east)))
+        end
+        left = leftenv(envs, source_site, south) *
+            TransferMatrix(north.AC[source_site], source_tensor, south.AC[source_site])
+        previous_col = source[2]
+        for target in right_targets
+            target_col = target[2]
+            for col in (previous_col + 1):(target_col - 1)
+                site = _window_mps_site(col, cache.colrange)
+                string_tensor = mpo_path_string(ρ[row, col, 1], stringspace, Val((:west, :east)))
+                left = left * TransferMatrix(north.AR[site], string_tensor, south.AR[site])
+            end
+
+            target_site = _window_mps_site(target_col, cache.colrange)
+            target_tensor = mpo_path_last(ρ[row, target_col, 1], mpo[2], Val(:west))
+            target_left = left * TransferMatrix(north.AR[target_site], target_tensor, south.AR[target_site])
+            value = _contract_transfer_boundaries(target_left, rightenv(envs, target_site, south))
+            numerators[targets[target]] = value
+
+            string_tensor = mpo_path_string(ρ[row, target_col, 1], stringspace, Val((:west, :east)))
+            left = left * TransferMatrix(north.AR[target_site], string_tensor, south.AR[target_site])
+            previous_col = target_col
+        end
+    end
+
+    # Close targets on the left of the incoming string from right to left
+    left_targets = [target for target in keys(targets) if target[2] < source[2]]
+    if !isempty(left_targets)
+        sort!(left_targets; by = x -> x[2], rev = true)
+        A = ρ[row, source[2], 1]
+        source_tensor = mpo_path_string(A, stringspace, Val((:north, :west)))
+        right = TransferMatrix(
+            north.AC[source_site], source_tensor, south.AC[source_site]
+        ) * rightenv(envs, source_site, south)
+        previous_col = source[2]
+        for target in left_targets
+            target_col = target[2]
+            for col in (previous_col - 1):-1:(target_col + 1)
+                site = _window_mps_site(col, cache.colrange)
+                string_tensor = mpo_path_string(ρ[row, col, 1], stringspace, Val((:east, :west)))
+                right = TransferMatrix(north.AL[site], string_tensor, south.AL[site]) * right
+            end
+
+            target_site = _window_mps_site(target_col, cache.colrange)
+            target_tensor = mpo_path_last(ρ[row, target_col, 1], mpo[2], Val(:east))
+            target_right = TransferMatrix(north.AL[target_site], target_tensor, south.AL[target_site]) * right
+            value = _contract_transfer_boundaries(leftenv(envs, target_site, south), target_right)
+            numerators[targets[target]] = value
+
+            string_tensor = mpo_path_string(ρ[row, target_col, 1], stringspace, Val((:east, :west)))
+            right = TransferMatrix(north.AL[target_site], string_tensor, south.AL[target_site]) * right
+            previous_col = target_col
+        end
+    end
+    return numerators
+end
+
+"""
+Map a PEPO column coordinate to its finite-MPS site number,
+which includes an additional west edge CTM tensor.
+"""
+_window_mps_site(col::Int, colrange::UnitRange{Int}) = col - first(colrange) + 2
+
+"""
+Build a finite row MPO by replacing one site tensor in one of the row MPOs without observables.
+"""
+function _row_mpo_with_site(
+        ρ::InfinitePEPO, tensor::MPOTensor, env::CTMRGEnv,
+        row::Int, col::Int, colrange::UnitRange{Int},
+    )
+    W = _row_mpo(ρ, nothing, env, row, colrange)
+    parent(W)[_window_mps_site(col, colrange)] = tensor
+    return W
+end
+
+"""
+Contract one modified row site between precomputed left and right MPS environments.
+"""
+function _contract_window_site(
+        envs::MPSKit.FiniteEnvironments, north::FiniteMPS, south::FiniteMPS,
+        site::Int, tensor::MPOTensor,
+    )
+    left = leftenv(envs, site, south) *
+        TransferMatrix(north.AC[site], tensor, south.AC[site])
+    return _contract_transfer_boundaries(left, rightenv(envs, site, south))
+end
+
+"""
+Contract the left and right transfer-matrix environments to a scalar.
+```
+    (north)
+    ┌-←-- 3 --←-┐
+    |           |
+    L-←-- 2 --←-R
+    |           |
+    └-→-- 1 --→-┘
+    (south)
+```
+"""
+function _contract_transfer_boundaries(left::MPSTensor, right::MPSTensor)
+    # The three bonds close around the window without crossing
+    return @plansor left[1 2; 3] * right[3 2; 1]
+end
diff --git a/src/algorithms/correlator_approx.jl b/src/algorithms/correlator_approx.jl
new file mode 100644
index 000000000..36499ba28
--- /dev/null
+++ b/src/algorithms/correlator_approx.jl
@@ -0,0 +1,57 @@
+# Approximate finite-window two-site correlators
+# ----------------------------------------------
+
+"""
+$(SIGNATURES)
+
+Approximately measure a dense two-site operator on one or more ordered pairs of sites in a single-layer PEPO.
+The first and second operator legs act on the first and second sites of each pair, respectively.
+Multiple pairs are evaluated in one shared window and reuse open-string boundary contractions.
+Ordered pairs in a batched call must be unique.
+"""
+function correlator_approx(
+        ρ::InfinitePEPO, op::AbstractTensorMap, bond::Tuple, env::CTMRGEnv;
+        trunc = _approx_trunc(env), maxiter::Int = 1, direction::Symbol = :auto,
+    )
+    return only(
+        correlator_approx(ρ, op, [bond], env; trunc, maxiter, direction)
+    )
+end
+
+function correlator_approx(
+        ρ::InfinitePEPO, op::AbstractTensorMap, bonds::AbstractVector,
+        env::CTMRGEnv;
+        trunc = _approx_trunc(env), maxiter::Int = 1, direction::Symbol = :auto,
+    )
+    bonds′ = _approx_twosite_bonds(bonds)
+    return _correlator_approx(
+        ρ, op, bonds′, env,
+        WindowApprox(Zipup(; trunc), _approx_dmrg(maxiter)), direction
+    )
+end
+
+"""
+Validate and regularize a nonempty collection of ordered two-site bonds to Cartesian indices.
+"""
+function _approx_twosite_bonds(bonds)
+    isempty(bonds) && throw(ArgumentError("correlator_approx requires at least one bond"))
+    bonds′ = NTuple{2, CartesianIndex{2}}[]
+    sizehint!(bonds′, length(bonds))
+    for bond in bonds
+        length(bond) == 2 || throw(ArgumentError("each bond should contain two sites"))
+        first_site = _mpo_observable_site(bond[1])
+        second_site = _mpo_observable_site(bond[2])
+        first_site != second_site ||
+            throw(ArgumentError("the sites of a bond should be distinct"))
+        push!(bonds′, (first_site, second_site))
+    end
+    allunique(bonds′) || throw(ArgumentError("bonds should be unique"))
+    return bonds′
+end
+
+"""
+Return the row and column ranges enclosing every endpoint in a collection of bonds.
+"""
+function _window_ranges(bonds::Vector{NTuple{2, CartesianIndex{2}}})
+    return _window_ranges(Iterators.flatten(bonds))
+end
diff --git a/src/algorithms/expval_approx.jl b/src/algorithms/expval_approx.jl
new file mode 100644
index 000000000..9876365d6
--- /dev/null
+++ b/src/algorithms/expval_approx.jl
@@ -0,0 +1,60 @@
+# Approximate finite-window expectation values
+# --------------------------------------------
+
+"""
+$(SIGNATURES)
+
+Approximately measure the expectation value of an open-boundary `MPOObservable` in a single-layer PEPO using finite boundary MPS/MPO zipup sweeps.
+The zipup truncation is controlled by `trunc`, which defaults to `truncrank(χ)` with `χ` the largest CTMRG boundary dimension.
+After each zipup step, the result is refined by a single-site DMRG approximation step with `maxiter` sweeps.
+Set `maxiter = 0` to disable this refinement.
+"""
+function expectation_value_approx(
+        ρ::InfinitePEPO, observable::MPOObservable, env::CTMRGEnv;
+        trunc = _approx_trunc(env), maxiter::Int = 1, direction::Symbol = :auto,
+    )
+    return _expectation_value_approx(
+        ρ, observable, env,
+        WindowApprox(Zipup(; trunc), _approx_dmrg(maxiter)), direction
+    )
+end
+
+"""
+Return the row and column ranges enclosing an MPO observable's complete routed path.
+"""
+function _window_ranges(observable::MPOObservable)
+    return _window_ranges(observable.path)
+end
+
+"""
+Return the smallest row and column ranges containing a collection of lattice sites.
+"""
+function _window_ranges(sites)
+    rows = getindex.(sites, 1)
+    cols = getindex.(sites, 2)
+    return UnitRange(extrema(rows)...), UnitRange(extrema(cols)...)
+end
+
+"""
+Return the largest CTMRG boundary-space dimension appearing in the corner tensors.
+"""
+function _ctmrg_boundary_chi(env::CTMRGEnv)
+    χ = 0
+    for C in env.corners
+        χ = max(χ, dim(space(C, 1)), dim(space(C, 2)))
+    end
+    return χ
+end
+
+"""
+Construct the default rank truncation from the largest CTMRG boundary dimension.
+"""
+_approx_trunc(env::CTMRGEnv) = truncrank(_ctmrg_boundary_chi(env))
+
+"""
+Construct the optional one-site DMRG refinement, or disable refinement for zero iterations.
+"""
+function _approx_dmrg(maxiter::Int)
+    maxiter >= 0 || throw(ArgumentError("maxiter should be nonnegative"))
+    return iszero(maxiter) ? nothing : DMRG(; maxiter, verbosity = 0)
+end
diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index 83c1a8ff1..652f39f1b 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -50,6 +50,19 @@ end
 # Default to Any for eltype: needs to be abstract anyways so not that much to gain
 LocalOperator(lattice, terms) = LocalOperator{Any}(lattice, terms)
 LocalOperator(lattice, terms::Pair...) = LocalOperator(lattice, terms)
+
+"""
+Sort operator sites using the default `CartesianIndex` ordering. When the order changes,
+permute the corresponding output and input physical legs by the same ordering.
+"""
+function _sort_op_sites(sites::Vector{CartesianIndex{2}}, op::AbstractTensorMap)
+    issorted(sites) && return sites, op
+    order = sortperm(sites)
+    sites′ = sites[order]
+    op′ = permute(op, (Tuple(order), Tuple(order) .+ numout(op)))
+    return sites′, op′
+end
+
 # TODO: add terms beyond AbstractTensorMap
 # e.g. tensor product of 1-site operators, MPOs
 add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap) = add_term!(operator, collect(inds), term)
@@ -68,12 +81,7 @@ function add_term!(
     end
     norm(term) <= atol && return operator # skip adding negligible terms
 
-    # permute input
-    if !issorted(inds)
-        I = sortperm(inds)
-        inds = inds[I]
-        term = permute(term, (Tuple(I), Tuple(I) .+ numout(term)))
-    end
+    inds, term = _sort_op_sites(inds, term)
 
     # translate coordinates
     _shift_into_unitcell!(inds, size(operator))
diff --git a/src/operators/mpo_observable.jl b/src/operators/mpo_observable.jl
new file mode 100644
index 000000000..561200f3c
--- /dev/null
+++ b/src/operators/mpo_observable.jl
@@ -0,0 +1,185 @@
+"""
+$(TYPEDEF)
+
+An open-boundary MPO embedded along a non-self-intersecting nearest-neighbor path on the
+square lattice. The tensor `mpo[k]` acts on `sites[k]`; `path` also contains intermediate
+sites that only carry the MPO string.
+
+The first and last MPO tensors use the reduced endpoint partitions `(1, 2)` and `(2, 1)`;
+all intermediate tensors use the standard `(2, 2)` MPO partition.
+"""
+struct MPOObservable{M}
+    sites::Vector{CartesianIndex{2}}
+    path::Vector{CartesianIndex{2}}
+    mpo::Vector{M}
+
+    function MPOObservable(
+            sites::Vector{CartesianIndex{2}}, path::Vector{CartesianIndex{2}},
+            mpo::Vector{M},
+        ) where {M}
+        _validate_mpo_observable(sites, path, mpo)
+        return new{M}(sites, path, mpo)
+    end
+end
+
+"""
+Construct an MPO observable from ordered operator sites and MPO tensors. The tensor `mpo[k]`
+acts on `sites[k]`; consecutive sites are connected by horizontal-first shortest paths.
+"""
+function MPOObservable(sites, mpo)
+    sites′ = CartesianIndex{2}[_mpo_observable_site(site) for site in sites]
+    mpo′ = collect(mpo)
+    path = _route_mpo_observable(sites′)
+    return MPOObservable(sites′, path, mpo′)
+end
+
+"""
+Construct an MPO observable from a dense operator. Operator sites are ordered in the same
+way as `LocalOperator` terms, and consecutive sites are connected by horizontal-first
+shortest paths. The lattice is a periodically indexed matrix of elementary physical spaces,
+as for `LocalOperator`.
+"""
+function MPOObservable(
+        sites, op::AbstractTensorMap, lattice::Matrix{<:ElementarySpace};
+        trunc = trunctol(; atol = MPSKit.Defaults.tol),
+    )
+    sites′ = CartesianIndex{2}[_mpo_observable_site(site) for site in sites]
+    length(sites′) >= 2 || throw(ArgumentError("an MPO observable requires at least two sites"))
+    allunique(sites′) || throw(ArgumentError("operator sites should be unique"))
+    length(sites′) == numin(op) == numout(op) ||
+        throw(ArgumentError("number of operator legs should match the number of sites"))
+
+    Nr, Nc = size(lattice)
+    sites′, op′ = _sort_op_sites(sites′, op)
+    for (i, site) in enumerate(sites′)
+        V = lattice[mod1(site[1], Nr), mod1(site[2], Nc)]
+        V == domain(op′)[i] == codomain(op′)[i] ||
+            throw(SpaceMismatch("operator physical space does not match lattice site $site"))
+    end
+
+    mpo = gate_to_mpo(op′; trunc)
+    length(mpo) == length(sites′) ||
+        throw(ArgumentError("expected an MPO decomposition matching the number of sites"))
+    return MPOObservable(sites′, mpo)
+end
+
+"""
+Normalize a supported lattice-site representation to a two-dimensional Cartesian index.
+"""
+_mpo_observable_site(site::CartesianIndex{2}) = site
+_mpo_observable_site(site::Tuple{Int, Int}) = CartesianIndex(site)
+_mpo_observable_site(site) =
+    throw(ArgumentError("MPO sites should be CartesianIndex{2} or (row, col) tuples"))
+
+"""
+Validate the topology, tensor partitions, physical placement, and adjacent string spaces of
+an already routed open-boundary MPO observable.
+"""
+function _validate_mpo_observable(sites, path, mpo)
+    length(sites) >= 2 || throw(ArgumentError("an MPO observable requires at least two sites"))
+    length(sites) == length(mpo) ||
+        throw(ArgumentError("the MPO should contain one tensor for every operator site"))
+    allunique(sites) || throw(ArgumentError("operator sites should be unique"))
+    allunique(path) || throw(ArgumentError("the MPO path should not intersect itself"))
+    all(op -> op isa AbstractTensorMap, mpo) ||
+        throw(ArgumentError("all MPO entries should be tensor maps"))
+
+    for (from, to) in zip(path, Iterators.drop(path, 1))
+        _step_direction(from, to)
+    end
+    path_positions = indexin(sites, path)
+    all(!isnothing, path_positions) && issorted(path_positions) ||
+        throw(ArgumentError("the MPO path should contain operator sites in MPO order"))
+    first(path_positions) == 1 && last(path_positions) == length(path) ||
+        throw(ArgumentError("the MPO path should start and end at operator sites"))
+
+    numout(first(mpo)) == 1 && numin(first(mpo)) == 2 ||
+        throw(ArgumentError("the first MPO tensor should have partition (1, 2)"))
+    numout(last(mpo)) == 2 && numin(last(mpo)) == 1 ||
+        throw(ArgumentError("the last MPO tensor should have partition (2, 1)"))
+    for op in @view mpo[2:(end - 1)]
+        numout(op) == 2 && numin(op) == 2 ||
+            throw(ArgumentError("middle MPO tensors should have partition (2, 2)"))
+    end
+    for k in 1:(length(mpo) - 1)
+        _mpo_right_stringspace(mpo[k])' == space(mpo[k + 1], 1) ||
+            throw(SpaceMismatch("incompatible MPO string spaces between sites $k and $(k + 1)"))
+    end
+    return nothing
+end
+
+function VI.scalartype(observable::MPOObservable)
+    return promote_type((scalartype(op) for op in observable.mpo)...)
+end
+
+"""
+Connect ordered operator sites by horizontal-first shortest paths while rejecting crossings
+and routes that pass through later operator sites.
+"""
+function _route_mpo_observable(sites)
+    length(sites) >= 2 || throw(ArgumentError("an MPO observable requires at least two sites"))
+    allunique(sites) || throw(ArgumentError("operator sites should be unique"))
+
+    path = CartesianIndex{2}[first(sites)]
+    measured = Set(sites)
+    visited = Set(path)
+
+    for k in 1:(length(sites) - 1)
+        segment = _l_path(sites[k], sites[k + 1])
+
+        for site in @view segment[2:(end - 1)]
+            site in measured &&
+                throw(ArgumentError("MPO path passes through another operator site"))
+            site in visited &&
+                throw(ArgumentError("MPO path intersects itself at site $site"))
+            push!(path, site)
+            push!(visited, site)
+        end
+
+        site = last(segment)
+        site in visited && throw(ArgumentError("MPO path intersects itself at site $site"))
+        push!(path, site)
+        push!(visited, site)
+    end
+    return path
+end
+
+"""
+Return the right MPO string space in the local tensor's stored leg orientation.
+"""
+_mpo_right_stringspace(op) = space(op, numind(op))
+
+"""
+Build a deterministic shortest nearest-neighbor path between two sites. The path traverses
+the horizontal separation first and then the vertical separation.
+"""
+function _l_path(start::CartesianIndex{2}, stop::CartesianIndex{2})
+    start == stop && throw(ArgumentError("MPO path sites should be unique"))
+
+    path = CartesianIndex{2}[start]
+    row, col = Tuple(start)
+    stoprow, stopcol = Tuple(stop)
+
+    while col != stopcol
+        col += sign(stopcol - col)
+        push!(path, CartesianIndex(row, col))
+    end
+    while row != stoprow
+        row += sign(stoprow - row)
+        push!(path, CartesianIndex(row, col))
+    end
+    return path
+end
+
+"""
+Return the cardinal direction of a nearest-neighbor step from one site to another. Throw an
+error when the sites are not nearest neighbors.
+"""
+function _step_direction(from::CartesianIndex{2}, to::CartesianIndex{2})
+    delta = to - from
+    delta == CartesianIndex(0, 1) && return :east
+    delta == CartesianIndex(0, -1) && return :west
+    delta == CartesianIndex(1, 0) && return :south
+    delta == CartesianIndex(-1, 0) && return :north
+    throw(ArgumentError("MPO path should use nearest-neighbor steps"))
+end
diff --git a/test/toolbox/correlator_approx.jl b/test/toolbox/correlator_approx.jl
new file mode 100644
index 000000000..28429f086
--- /dev/null
+++ b/test/toolbox/correlator_approx.jl
@@ -0,0 +1,112 @@
+using TensorKit
+using PEPSKit
+using MPSKit
+using Test
+using Random
+
+const CI = CartesianIndex
+
+"""
+Contract `⟨op⟩` on each bond in `bonds` independently without caching
+"""
+function _shared_window_reference(
+        op::AbstractTensorMap, bonds::Vector{NTuple{2, CI{2}}}, ρ, env
+    )
+    lattice = physicalspace(ρ)
+    observables = [MPOObservable(collect(pair), op, lattice) for pair in bonds]
+    rowrange, colrange = PEPSKit._window_ranges(bonds)
+    alg = PEPSKit.WindowApprox(Zipup(; trunc = notrunc()), nothing)
+    norm = PEPSKit._contract_window_rows(
+        ρ, nothing, env, rowrange, colrange, alg
+    )
+    return map(observables) do observable
+        numerator = PEPSKit._contract_window_rows(
+            ρ, observable, env, rowrange, colrange, alg
+        )
+        return numerator / norm
+    end
+end
+
+bonds = [
+    # source at (1, 1), bond sites in order
+    ## target in the same row as source
+    (CI(1, 1), CI(1, 2)),
+    ## target in a later row, on both sides of or just below the source
+    (CI(1, 1), CI(2, 0)), (CI(1, 1), CI(2, 1)), (CI(1, 1), CI(2, 2)),
+    # source still at (1, 1), but bond sites need swapping
+    (CI(2, 2), CI(1, 1)),
+    # another source at (1, 2), bond sites in order
+    (CI(1, 2), CI(2, 0)), (CI(1, 2), CI(2, 2)),
+]
+
+spaces = Dict(
+    U1Irrep => (
+        U1Space(1 => 2, -1 => 1),
+        U1Space(1 => 1, 0 => 1, -1 => 2),
+        U1Space(1 => 1, 0 => 1, -1 => 2),
+    ),
+    FermionParity => (
+        Vect[FermionParity](0 => 1, 1 => 1),
+        Vect[FermionParity](0 => 1, 1 => 2),
+        Vect[FermionParity](0 => 2, 1 => 2),
+    )
+)
+
+@testset "Two-site source grouping" begin
+    groups = PEPSKit._twosite_source_groups(bonds)
+    @test length(groups) == 3 &&
+        groups[(CI(1, 1), false)][CI(2, 2)] == 4 &&
+        haskey(groups, (CI(1, 1), true))
+end
+
+@testset "Single-layer PEPO ($S)" for S in keys(spaces)
+    Random.seed!(1234)
+    d, D, χ = spaces[S]
+    ρ = InfinitePEPO(d, D; unitcell = (2, 2, 1))
+    lattice = physicalspace(ρ)
+    env = CTMRGEnv(InfinitePartitionFunction(ρ), χ)
+    trunc = notrunc()
+
+    O² = rand(ComplexF64, d^2, d^2)
+    id² = isomorphism(d, d) ⊗ isomorphism(d, d)
+
+    # bonds to be measured should be unique
+    @test_throws ArgumentError correlator_approx(ρ, O², fill(bonds[3], 2), env)
+
+    exact_same_row = expectation_value(ρ, LocalOperator(lattice, bonds[1] => O²), env)
+    @test correlator_approx(
+        ρ, O², bonds[1], env; trunc, maxiter = 0, direction = :rows
+    ) ≈ exact_same_row
+
+    exact_reversed = expectation_value(ρ, LocalOperator(lattice, bonds[5] => O²), env)
+    @test correlator_approx(
+        ρ, O², bonds[5], env; trunc, maxiter = 0, direction = :columns
+    ) ≈ exact_reversed
+
+    vals_ref = _shared_window_reference(O², bonds, ρ, env)
+    vals_rows = correlator_approx(
+        ρ, O², bonds, env; trunc, maxiter = 0, direction = :rows
+    )
+    @test vals_rows ≈ vals_ref
+    vals_columns = correlator_approx(
+        ρ, O², bonds, env; trunc, maxiter = 0, direction = :columns
+    )
+    @test vals_columns ≈ vals_rows
+    vals_auto = correlator_approx(
+        ρ, O², bonds, env; trunc, maxiter = 0, direction = :auto
+    )
+    @test vals_auto ≈ vals_columns
+
+    @test correlator_approx(ρ, id², bonds, env; trunc, direction = :rows) ≈
+        ones(length(bonds))
+
+    alg = PEPSKit.WindowApprox(Zipup(; trunc), nothing)
+    cache = PEPSKit._window_row_cache(ρ, env, 1:2, 1:2, alg)
+    @test all(eachindex(cache.north_prefixes)) do k
+        dot(cache.south_suffixes[k], cache.north_prefixes[k]) ≈ cache.norm
+    end
+
+    W = cache.row_mpos[first(cache.rowrange)]
+    W_adjoint = PEPSKit._adjoint_mpo(W)
+    @test convert(TensorMap, W_adjoint) ≈ convert(TensorMap, W)'
+end
diff --git a/test/toolbox/correlator_approx_phys.jl b/test/toolbox/correlator_approx_phys.jl
new file mode 100644
index 000000000..62aad5ea7
--- /dev/null
+++ b/test/toolbox/correlator_approx_phys.jl
@@ -0,0 +1,39 @@
+using TensorKit
+using PEPSKit
+using Test
+using TensorKitTensors.SpinOperators: S_exchange
+
+const CI = CartesianIndex
+
+@testset "correlator_approx for physical state" begin
+    Nr, Nc = 2, 2
+    lattice = InfiniteSquare(Nr, Nc)
+    sym = Trivial
+    ham = j1_j2_model(Float64, sym, lattice; J1 = 1.0, J2 = 0.0, sublattice = false)
+    op = S_exchange(Float64, sym)
+    lattice = physicalspace(ham)
+
+    ρ = PEPSKit.infinite_temperature_density_matrix(ham)
+    state_trunc = truncrank(4) & truncerror(; atol = 1.0e-12)
+    su_alg = SimpleUpdate(; trunc = state_trunc, purified = false)
+    ρ, = time_evolve(ρ, ham, 1.0e-2, 50, su_alg, SUWeight(ρ))
+
+    network = InfinitePartitionFunction(ρ)
+    env = initialize_ctmrg_environment(network, ProductStateInitialization())
+    env_trunc = truncrank(8) & truncerror(; atol = 1.0e-12)
+    env, = leading_boundary(env, network; alg = :SequentialCTMRG, trunc = env_trunc)
+
+    bonds = [
+        (CI(1, 1), CI(1, 2)),
+        (CI(1, 1), CI(2, 0)), (CI(1, 1), CI(2, 1)), (CI(1, 1), CI(2, 2)),
+        (CI(1, 1), CI(3, 2)),
+    ]
+    cor_exact = map(bonds) do bond
+        O = LocalOperator(lattice, bond => op)
+        return expectation_value(ρ, O, env)
+    end
+    cor_trunc = correlator_approx(ρ, op, bonds, env; trunc = env_trunc, maxiter = 1)
+    @info "Exact:" cor_exact
+    @info "Approx:" cor_trunc
+    @test cor_trunc ≈ cor_exact rtol = 1.0e-3
+end
diff --git a/test/toolbox/expval_approx.jl b/test/toolbox/expval_approx.jl
new file mode 100644
index 000000000..2f754c74b
--- /dev/null
+++ b/test/toolbox/expval_approx.jl
@@ -0,0 +1,52 @@
+using TensorKit
+using PEPSKit
+using MPSKit
+using Test
+using Random
+
+const CI = CartesianIndex
+
+spaces = Dict(
+    U1Irrep => (
+        U1Space(1 => 2, -1 => 1),
+        U1Space(1 => 1, 0 => 1, -1 => 2),
+        U1Space(1 => 1, 0 => 1, -1 => 2),
+    ),
+    FermionParity => (
+        Vect[FermionParity](0 => 1, 1 => 1),
+        Vect[FermionParity](0 => 1, 1 => 2),
+        Vect[FermionParity](0 => 2, 1 => 2),
+    ),
+)
+
+sites_list = (
+    [CI(1, 1), CI(1, 2)], # horizontal
+    [CI(1, 1), CI(2, 1)], # vertical
+    [CI(1, 1), CI(2, 2)], # turned
+    [CI(2, 2), CI(1, 1)], # reversed turned
+    [CI(2, 1), CI(1, 1), CI(1, 2), CI(2, 2)], # U-shaped
+)
+
+@testset "Single-layer PEPO ($S)" for S in keys(spaces)
+    Random.seed!(1234)
+
+    d, D, χ = spaces[S]
+    ρ = InfinitePEPO(d, D; unitcell = (2, 2, 1))
+    env = CTMRGEnv(InfinitePartitionFunction(ρ), χ)
+    trunc = notrunc()
+
+    for sites in sites_list
+        n = length(sites)
+        op = randn(ComplexF64, d^n → d^n)
+        mpo = PEPSKit.gate_to_mpo(op; trunc)
+        observable = MPOObservable(sites, mpo)
+        exact = expectation_value(
+            ρ, LocalOperator(physicalspace(ρ), sites => op), env
+        )
+        for direction in (:rows, :columns)
+            @test expectation_value_approx(
+                ρ, observable, env; trunc, maxiter = 0, direction
+            ) ≈ exact
+        end
+    end
+end
diff --git a/test/toolbox/mpo_routing.jl b/test/toolbox/mpo_routing.jl
new file mode 100644
index 000000000..2b9a4a1ca
--- /dev/null
+++ b/test/toolbox/mpo_routing.jl
@@ -0,0 +1,93 @@
+using PEPSKit
+using TensorKit
+using Test
+using Random
+
+const directions = (:north, :east, :south, :west)
+
+spaces = Dict(
+    U1Irrep => (
+        Rep[U₁](0 => 1, 1 => 1),
+        Rep[U₁](0 => 1, 1 => 1, -1 => 1),
+        Rep[U₁](1 => 1),
+    ),
+    FermionParity => (
+        Vect[FermionParity](0 => 1, 1 => 1),
+        Vect[FermionParity](0 => 1, 1 => 1),
+        Vect[FermionParity](1 => 1),
+    ),
+)
+
+@testset "Fuser flips and twists" begin
+    d = Vect[FermionParity](0 => 1, 1 => 1)
+    D = Vect[FermionParity](0 => 1, 1 => 1)
+    ρ = InfinitePEPO(d, D; unitcell = (2, 2, 1))
+    A = ρ[1, 1, 1]
+    Bh = ρ[1, 2, 1]
+    Bv = ρ[2, 1, 1]
+    A′ = PEPSKit.twistdual(A, 2)
+    Bh′ = PEPSKit.twistdual(Bh, 2)
+    Bv′ = PEPSKit.twistdual(Bv, 2)
+
+    op = randn(ComplexF64, d^2 → d^2)
+    mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc())
+
+    first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(:east))
+    last_tensor = PEPSKit.mpo_path_last(Bh, last(mpo), Val(:west))
+    @tensor exact[W1 S1 S2; N1 N2 E2] := op[po1 po2; pi1 pi2] *
+        A′[pi1 po1; N1 x S1 W1] * Bh′[pi2 po2; N2 E2 S2 x]
+    @tensor routed[W1 S1 S2; N1 N2 E2] :=
+        first_tensor[W1 S1; N1 x] * last_tensor[x S2; N2 E2]
+    @test routed ≈ exact
+
+    first_tensor = PEPSKit.mpo_path_first(Bh, first(mpo), Val(:west))
+    last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(:east))
+    @tensor exact[W2 S1 S2; N1 E1 N2] := op[po1 po2; pi1 pi2] *
+        Bh′[pi1 po1; N1 E1 S1 x] * A′[pi2 po2; N2 x S2 W2]
+    @tensor routed[W2 S1 S2; N1 E1 N2] :=
+        first_tensor[x S1; N1 E1] * last_tensor[W2 S2; N2 x]
+    @test routed ≈ exact
+
+    first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(:south))
+    last_tensor = PEPSKit.mpo_path_last(Bv, last(mpo), Val(:north))
+    @tensor exact[W1 W2 S2; N1 E1 E2] := op[po1 po2; pi1 pi2] *
+        A′[pi1 po1; N1 E1 x W1] * Bv′[pi2 po2; x E2 S2 W2]
+    @tensor routed[W1 W2 S2; N1 E1 E2] :=
+        first_tensor[W1 x; N1 E1] * last_tensor[W2 S2; x E2]
+    @test routed ≈ exact
+
+    first_tensor = PEPSKit.mpo_path_first(Bv, first(mpo), Val(:north))
+    last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(:south))
+    @tensor exact[W1 S1 W2; E1 N2 E2] := op[po1 po2; pi1 pi2] *
+        Bv′[pi1 po1; x E1 S1 W1] * A′[pi2 po2; N2 E2 x W2]
+    @tensor routed[W1 S1 W2; E1 N2 E2] :=
+        first_tensor[W1 S1; x E1] * last_tensor[W2 x; N2 E2]
+    @test routed ≈ exact
+end
+
+@testset "Routing identities ($S)" for S in keys(spaces)
+    Random.seed!(1234)
+    d, D, stringspace = spaces[S]
+    ρ = InfinitePEPO(d, D; unitcell = (1, 1, 1))
+    op = rand(ComplexF64, d^2 → d^2)
+    mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc())
+
+    A = ρ[1, 1, 1]
+    for direction in directions
+        first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(direction))
+        last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(direction))
+        @test (numout(first_tensor), numin(first_tensor)) == (2, 2)
+        @test (numout(last_tensor), numin(last_tensor)) == (2, 2)
+    end
+
+    middle = TensorMap(TensorKit.BraidingTensor{ComplexF64}(d, stringspace))
+    for incoming in directions, outgoing in directions
+        incoming == outgoing && continue
+        tensor = PEPSKit.mpo_path_middle(A, middle, Val((incoming, outgoing)))
+        string_tensor = PEPSKit.mpo_path_string(
+            A, stringspace, Val((incoming, outgoing))
+        )
+        @test (numout(tensor), numin(tensor)) == (2, 2)
+        @test string_tensor ≈ tensor
+    end
+end
diff --git a/test/types/mpo_observable.jl b/test/types/mpo_observable.jl
new file mode 100644
index 000000000..f6498a74e
--- /dev/null
+++ b/test/types/mpo_observable.jl
@@ -0,0 +1,98 @@
+using PEPSKit
+using TensorKit
+using Test
+using Random
+
+const CI = CartesianIndex
+
+function _check_observable_bonds(observable::MPOObservable)
+    @test length(observable.sites) == length(observable.mpo)
+    path_positions = indexin(observable.sites, observable.path)
+    @test all(!isnothing, path_positions)
+    @test issorted(path_positions)
+    @test first(path_positions) == 1
+    @test last(path_positions) == length(observable.path)
+    for k in 1:(length(observable.mpo) - 1)
+        @test PEPSKit._mpo_right_stringspace(observable.mpo[k])' ==
+            space(observable.mpo[k + 1], 1)
+    end
+    return nothing
+end
+
+Random.seed!(1234)
+d = ℂ^2
+lattice = fill(d, 2, 2)
+
+@testset "Manhattan paths and directions" begin
+    # first move horizontally
+    @test PEPSKit._l_path(CI(2, 1), CI(2, 4)) ==
+        CI.([(2, 1), (2, 2), (2, 3), (2, 4)])
+    @test PEPSKit._l_path(CI(1, 3), CI(4, 3)) ==
+        CI.([(1, 3), (2, 3), (3, 3), (4, 3)])
+    @test PEPSKit._l_path(CI(3, 4), CI(1, 1)) ==
+        CI.([(3, 4), (3, 3), (3, 2), (3, 1), (2, 1), (1, 1)])
+    # nearest neighbor directions
+    origin = CI(2, 2)
+    @test PEPSKit._step_direction(origin, CI(1, 2)) === :north
+    @test PEPSKit._step_direction(origin, CI(2, 3)) === :east
+    @test PEPSKit._step_direction(origin, CI(3, 2)) === :south
+    @test PEPSKit._step_direction(origin, CI(2, 1)) === :west
+end
+
+@testset "Routed MPO tensors" begin
+    op2 = rand(ComplexF64, d^2, d^2)
+    observable2 = MPOObservable([CI(1, 1), CI(3, 3)], op2, lattice)
+    @test observable2.path ==
+        CI.([(1, 1), (1, 2), (1, 3), (2, 3), (3, 3)])
+    @test observable2.sites == CI.([(1, 1), (3, 3)])
+    _check_observable_bonds(observable2)
+
+    original2 = PEPSKit.gate_to_mpo(op2)
+    @test first(observable2.mpo) ≈ first(original2)
+    @test last(observable2.mpo) ≈ last(original2)
+    @test length(observable2.mpo) == 2
+
+    sites4 = [CI(4, 2), CI(1, 3), CI(3, 4), CI(1, 1)]
+    op4 = rand(ComplexF64, d^4, d^4)
+    observable4 = MPOObservable(sites4, op4, lattice)
+    @test observable4.path == CI.(
+        [
+            (1, 1), (1, 2), (2, 2), (3, 2), (4, 2), (4, 3),
+            (3, 3), (2, 3), (1, 3), (1, 4), (2, 4), (3, 4),
+        ]
+    )
+    _check_observable_bonds(observable4)
+
+    plaquette_sites = CI.([(2, 1), (1, 1), (1, 2), (2, 2)])
+    plaquette_mpo = PEPSKit.gate_to_mpo(rand(ComplexF64, d^4, d^4))
+    plaquette = MPOObservable(plaquette_sites, plaquette_sites, plaquette_mpo)
+    @test plaquette.sites == plaquette_sites
+    @test plaquette.path == plaquette_sites
+    _check_observable_bonds(plaquette)
+end
+
+@testset "Explicit MPOs and validation" begin
+    op = rand(ComplexF64, d^2, d^2)
+    mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc())
+
+    observable = MPOObservable([(2, 2), (1, 1)], mpo)
+    @test observable.sites == CI.([(2, 2), (1, 1)])
+    @test observable.path == CI.([(2, 2), (2, 1), (1, 1)])
+    @test observable.mpo !== mpo
+    _check_observable_bonds(observable)
+
+    # Invalid MPO with virtual space mismatch
+    rank_one_mpo = PEPSKit.gate_to_mpo(op; trunc = truncrank(1))
+    @test_throws SpaceMismatch MPOObservable(
+        [CI(1, 1), CI(1, 2)], [first(mpo), last(rank_one_mpo)]
+    )
+
+    # the observable must act on two diffrent sites
+    @test_throws ArgumentError MPOObservable([CI(1, 1), CI(1, 1)], mpo)
+    @test_throws ArgumentError MPOObservable([CI(1, 1), CI(1, 1)], op, lattice)
+    # Routing between consecutive tensors must not cross a later operator site.
+    @test_throws ArgumentError MPOObservable(
+        [CI(1, 1), CI(1, 3), CI(1, 2)],
+        PEPSKit.gate_to_mpo(rand(ComplexF64, d^3, d^3)),
+    )
+end
diff --git a/test/types/standardize_dualness.jl b/test/types/standardize_dualness.jl
new file mode 100644
index 000000000..5b394e39f
--- /dev/null
+++ b/test/types/standardize_dualness.jl
@@ -0,0 +1,39 @@
+using TensorKit
+using PEPSKit
+using Test
+
+@testset "Standardize virtual-space dualness" begin
+    d = U1Space(0 => 1)
+    D = U1Space(0 => 1)
+    χ = U1Space(0 => 1)
+    Pspaces = fill(d, 2, 2, 1)
+
+    patterns = (
+        (fill(D', 2, 2, 1), fill(D', 2, 2, 1)),
+        (
+            reshape([D, D', D', D], 2, 2, 1),
+            reshape([D', D, D, D'], 2, 2, 1),
+        ),
+    )
+    for (Nspaces, Espaces) in patterns
+        ρ = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces)
+        env = CTMRGEnv(randn, ComplexF64, InfinitePartitionFunction(ρ), χ)
+        standardized_ρ, standardized_env = PEPSKit.standardize_dualness(ρ, env)
+
+        for row in axes(ρ, 1), col in axes(ρ, 2)
+            A = standardized_ρ[row, col, 1]
+            edge_coordinates = (
+                PEPSKit.NORTH => (row - 1, col),
+                PEPSKit.EAST => (row, col + 1),
+                PEPSKit.SOUTH => (row + 1, col),
+                PEPSKit.WEST => (row, col - 1),
+            )
+            for (direction, coordinates) in edge_coordinates
+                V = PEPSKit.virtualspace(A, direction)
+                E = PEPSKit.edge(standardized_env, direction, coordinates...)
+                @test isdual(V) == (direction in (PEPSKit.NORTH, PEPSKit.EAST))
+                @test V == space(E, 2)'
+            end
+        end
+    end
+end