From aa4535d1927e514d122070e437576ec440cc93db Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 5 Aug 2026 10:34:39 +0200
Subject: [PATCH 01/11] resolve duplicated level

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

diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 8745146ab..9bf4abdfe 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -313,7 +313,7 @@ function (H::PrecomputedAC2Derivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 3})
     @plansor backend = backend allocator = allocator begin
         y_braided[-1 -2; -4 -6 -5 -3] := L[-1 -2; 1 2] * xR_braided[1 2; -4 -6 -5 -3]
     end
-    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 5, 5, 3))
+    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 6, 5, 3))
 end
 
 const _ToPrepare = Union{

From 90ed1e7c2b92d457d2d5d8884026bded6205b88a Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 5 Aug 2026 10:50:35 +0200
Subject: [PATCH 02/11] invert some braiding tensors + densify to braid with
 sumspace

---
 src/algorithms/expval.jl       | 3 ++-
 src/operators/mpo.jl           | 3 ++-
 src/transfermatrix/transfer.jl | 7 ++++---
 3 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index 191d2f93c..90674f712 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -123,6 +123,7 @@ function contract_mpo_expval2(
     return @plansor conj(A1bar[4 5; 6]) * conj(A2bar[6 7; 8]) * O[5 7; 2 3] * A1[4 2; 1] *
         A2[1 3; 8]
 end
+# the ancilla braid need to be inverted due to handedness of mpo.jl#L90
 function contract_mpo_expval2(
         A1::GenericMPSTensor{S, 3}, A2::GenericMPSTensor{S, 3},
         O::AbstractTensorMap{<:Any, S},
@@ -130,7 +131,7 @@ function contract_mpo_expval2(
     ) where {S}
     numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator"))
     return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ[9 6; 1 2] *
-        τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7]
+        τ'[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7]
 end
 
 function expectation_value(
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index 225c83756..36a90a0d3 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -87,7 +87,8 @@ function Base.convert(::Type{<:InfiniteMPS}, mpo::InfiniteMPO)
     return InfiniteMPS(map(_mpo_to_mps, parent(mpo)))
 end
 function _mpo_to_mps(O::MPOTensor)
-    @plansor A[-1 -2 -3; -4] := O[-1 -2; 1 2] * τ[1 2; -4 -3]
+    O′ = O isa AbstractBlockTensorMap ? TensorMap(O) : O
+    @plansor A[-1 -2 -3; -4] := O′[-1 -2; 1 2] * τ[1 2; -4 -3]
     return A isa AbstractBlockTensorMap ? TensorMap(A) : A
 end
 
diff --git a/src/transfermatrix/transfer.jl b/src/transfermatrix/transfer.jl
index aa3d10ed7..09b98d861 100644
--- a/src/transfermatrix/transfer.jl
+++ b/src/transfermatrix/transfer.jl
@@ -146,16 +146,17 @@ function transfer_right(v::MPOTensor, O::MPOTensor, A::MPSTensor, Ab::MPSTensor)
         conj(Ab[-4 2; 1]) * v[5 3; -3 1]
 end
 
-# desity matrix transfer
+# density matrix transfer
+# these braids need to be inverted due to handedness of mpo.jl#L90
 function transfer_left(
         x::MPSTensor, O::MPOTensor, A::GenericMPSTensor{<:Any, 3}, Ab::GenericMPSTensor{<:Any, 3}
     )
-    return @plansor y[-1 -2; -3] ≔ x[1 2; 6] * A[6 7 8; -3] * O[2 3; 7 5] * τ[5 4; 8 -2] *
+    return @plansor y[-1 -2; -3] ≔ x[1 2; 6] * A[6 7 8; -3] * O[2 3; 7 5] * τ'[5 4; 8 -2] *
         conj(Ab[1 3 4; -1])
 end
 function transfer_right(
         x::MPSTensor, O::MPOTensor, A::GenericMPSTensor{<:Any, 3}, Ab::GenericMPSTensor{<:Any, 3}
     )
-    return @plansor y[-1 -2; -3] ≔ A[-1 4 2; 1] * O[-2 6; 4 5] * τ[5 7; 2 3] *
+    return @plansor y[-1 -2; -3] ≔ A[-1 4 2; 1] * O[-2 6; 4 5] * τ'[5 7; 2 3] *
         conj(Ab[-3 6 7; 8]) * x[1 3; 8]
 end

From e4e40018339793c3a55a17adac916326dcb33b28 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 18:29:51 +0200
Subject: [PATCH 03/11] unprepared and prepared AC(2)ham braids

---
 src/algorithms/derivatives/mpo_derivatives.jl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 9bf4abdfe..12b7c845e 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -121,13 +121,14 @@ function (h::MPO_AC_Hamiltonian{<:MPSBondTensor, Nothing, <:MPSBondTensor})(x::M
     end
     return y
 end
+# braid needs to be inverted due to handedness of mpo.jl#L90
 function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})(
         x::GenericMPSTensor{<:Any, 3}
     )
     backend, allocator = h.backend, h.allocator
     @plansor backend = backend allocator = allocator begin
         y[-1 -2 -3; -4] ≔ h.leftenv[-1 7; 6] * x[6 4 2; 1] *
-            h.operators[1][7 -2; 4 5] * τ[5 -3; 2 3] * h.rightenv[1 3; -4]
+            h.operators[1][7 -2; 4 5] * τ'[5 -3; 2 3] * h.rightenv[1 3; -4]
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
@@ -151,14 +152,15 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
+# these braids need to be inverted due to handedness of mpo.jl#L90
 function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})(
         x::AbstractTensorMap{<:Any, <:Any, 3, 3}
     )
     backend, allocator = h.backend, h.allocator
     @plansor backend = backend allocator = allocator begin
         y[-1 -2 -3; -4 -5 -6] ≔ h.leftenv[-1 11; 10] * x[10 8 6; 1 2 4] *
-            h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ[9 -3; 6 7] *
-            h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3]
+            h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ'[9 -3; 6 7] *
+            h.operators[2][7 -6; 4 5] * τ'[5 -5; 2 3]
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
@@ -290,9 +292,9 @@ end
 function (H::PrecomputedACDerivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 1})
     backend, allocator = H.backend, H.allocator
     L, R = H.leftenv, H.rightenv
-
+    # braid needs to be inverted due to handedness of mpo.jl#L90
     @plansor backend = backend allocator = allocator begin
-        xR[-1 -2; -4 -5 -3] := x[-1 -2 3; 1] * R[1 2; -4] * τ[2 3; -5 -3]
+        xR[-1 -2; -4 -5 -3] := x[-1 -2 3; 1] * R[1 2; -4] * τ'[2 3; -5 -3]
     end
     xR_fused = fuse_legs(xR, 2, 1)
     @plansor backend = backend allocator = allocator begin
@@ -304,16 +306,16 @@ function (H::PrecomputedAC2Derivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 3})
     backend, allocator = H.backend, H.allocator
     L, R = H.leftenv, H.rightenv
 
-    x_braided = fuse_legs(braid(x, ((1, 2, 3, 5), (4, 6)), (1, 2, 3, 4, 5, 6)), 1, 2)
+    x_braided = fuse_legs(braid(x, ((1, 2, 3, 5), (4, 6)), (1, 2, 3, 4, 6, 5)), 1, 2)
 
     @plansor backend = backend allocator = allocator begin
         xR[-1 -2; -4 -6 -7 -5 -3] := x_braided[-1 -2 -3 -5; 1] * R[1 -7; -4 -6]
     end
-    @notensor xR_braided = braid(fuse_legs(xR, 2, 1), ((1, 4), (2, 3, 5, 6)), (1, 4, 6, 7, 5, 3))
+    @notensor xR_braided = braid(fuse_legs(xR, 2, 1), ((1, 4), (2, 3, 5, 6)), (1, 2, 3, 5, 6, 7))
     @plansor backend = backend allocator = allocator begin
         y_braided[-1 -2; -4 -6 -5 -3] := L[-1 -2; 1 2] * xR_braided[1 2; -4 -6 -5 -3]
     end
-    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 6, 5, 3))
+    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 5, 6, 3))
 end
 
 const _ToPrepare = Union{

From 7f141f576fad9eb7bd5753e57d5a3198f44e8693 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 18:41:02 +0200
Subject: [PATCH 04/11] tests that catch these fixes

---
 test/algorithms/anyonic_braiding.jl | 180 ++++++++++++++++++++++++++++
 1 file changed, 180 insertions(+)
 create mode 100644 test/algorithms/anyonic_braiding.jl

diff --git a/test/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl
new file mode 100644
index 000000000..ef80a9db6
--- /dev/null
+++ b/test/algorithms/anyonic_braiding.jl
@@ -0,0 +1,180 @@
+println("
+------------------------------------
+|   Anyonic braiding consistency    |
+------------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using BlockTensorKit
+using MPSKit: AC_hamiltonian, AC2_hamiltonian, AC2
+using TensorKit
+using KrylovKit: eigsolve
+using LinearAlgebra
+using Random
+
+# when reinterpreting an MPO as an MPS ρ, there are contractions which
+# braid the ancilla leg past the operator's virtual leg
+# this turned out to be incorrect in several places
+# the tests below check that the braiding is consistent with the dense trace for finite systems,
+# and with the MPO-level overlap transfer for infinite systems
+
+@testset "density-matrix sandwich (Fibonacci)" begin
+    V = Vect[FibonacciAnyon](:I => 1, :τ => 1)
+    ρ = randn(ComplexF64, V ⊗ V, V ⊗ V)
+    O = randn(ComplexF64, V ⊗ V, V ⊗ V)
+    ρ_mps = convert(FiniteMPS, FiniteMPO(ρ))
+
+    ref = tr(ρ' * O * ρ)
+    @test dot(FiniteMPO(ρ), FiniteMPO(O) * FiniteMPO(ρ)) ≈ ref
+    @test dot(ρ_mps, FiniteMPO(O), ρ_mps) ≈ ref
+end
+
+@testset "expectation values with density matrices (Fibonacci)" begin
+    V = Vect[FibonacciAnyon](:I => 1, :τ => 1)
+    L = 4
+    h = randn(ComplexF64, V ⊗ V, V ⊗ V)
+    h = h + h'
+    H = FiniteMPOHamiltonian(fill(V, L), [(i, i + 1) => h for i in 1:(L - 1)])
+    ρ = make_time_mpo(H, 0.1, TaylorCluster(; N = 2); imaginary_evolution = true)
+    ρd, Hd = convert(TensorMap, ρ), convert(TensorMap, H)
+    ref = tr(ρd' * Hd * ρd) / tr(ρd' * ρd)
+
+    # expectation_value here attempts to braid ancilla leg with mpo virtual leg,
+    # which is a sumspace/plain pair and thus fails
+    @test_broken expectation_value(ρ, H) ≈ ref
+
+    # free workaround: flatten jordan tensors into plain mpo tensors
+    # don't convert to tensormap, this is an exponential cost in L
+    H_dense = FiniteMPO([W isa AbstractBlockTensorMap ? TensorMap(W) : W for W in parent(H)])
+    @test expectation_value(ρ, H_dense) ≈ ref
+
+    # these require the densification in _mpo_to_mps (mpo.jl#90)
+    @test expectation_value(ρ, FiniteMPO(Hd)) ≈ ref
+    @test expectation_value(ρ, (1, 2) => h) + expectation_value(ρ, (2, 3) => h) +
+        expectation_value(ρ, (3, 4) => h) ≈ ref
+end
+
+@testset "braiding consistency in derivative operators (Fibonacci)" begin
+    V = Vect[FibonacciAnyon](:I => 1, :τ => 1)
+    L = 3
+    T = ComplexF64
+    ρ, O = FiniteMPO(randn(T, V^L, V^L)), FiniteMPO(randn(T, V^L, V^L))
+    ρ_mps = convert(FiniteMPS, ρ)
+    exact = convert(FiniteMPS, O * ρ)
+    ψ0 = FiniteMPS(randn, T, fill(V ⊗ V', L), Vect[FibonacciAnyon](:I => 24, :τ => 24))
+    fid(ψ) = abs(dot(ψ, exact)) / (norm(ψ) * norm(exact))
+
+    # unprepared AC reached by one-site approximate
+    ψ1, = approximate(ψ0, (O, ρ_mps), DMRG(; tol = 1.0e-10, maxiter = 50)) # simply doesn't converge without the fix
+    @test fid(ψ1) ≈ 1 atol = 1.0e-6
+
+    # unprepared AC2 reached by two-site approximate
+    ψ2, = approximate(ψ0, (O, ρ_mps), DMRG2(; tol = 1.0e-10, maxiter = 50, trunc = truncrank(64))) # simply doesn't converge without the fix
+    @test fid(ψ2) ≈ 1 atol = 1.0e-6
+
+    # prepared AC, reached by any eigensolver
+    # τ' fix on mpo_derivatives.jl#252 allows the eigensolver to converge at all,
+    # even though the eigenproblem is degenerate over the ancilla leg
+    A = randn(T, V^L, V^L)
+    Oh = FiniteMPO(A + A')
+    ψ3, _, ϵ3 = find_groundstate(ρ_mps, Oh, DMRG(; tol = 1.0e-10, maxiter = 50))
+    @test ϵ3 ≤ 1.0e-10
+
+    # prepared AC2, reached by two-site TDVP
+    Ohd = convert(TensorMap, Oh)
+    ρd = randn(T, V^L, V^L)
+    ρ_mps2 = convert(FiniteMPS, FiniteMPO(ρd))
+    dt = 0.02
+    ρ_exact_d = exp(-im * dt * Ohd) * ρd
+    target = convert(FiniteMPS, FiniteMPO(ρ_exact_d))
+    ψt, = timestep(ρ_mps2, Oh, 0.0, dt, TDVP2(; trunc = truncrank(64)))
+    tdvp2_fid = abs(dot(ψt, target)) / (norm(ψt) * norm(target))
+    @test tdvp2_fid ≈ 1 atol = 1.0e-6 # fidelity was previously high, but not ≈ 1
+end
+
+const braid_spaces = (
+    ℂ^2,
+    Vect[Z2Irrep](0 => 1, 1 => 1),
+    Vect[FermionParity](0 => 1, 1 => 1),
+    Vect[FibonacciAnyon](:I => 1, :τ => 1),
+)
+
+# build L-site operator that's the identity everywhere except at site i
+# where it's op with physical space V
+_embed(L, op, i) = foldl(⊗, (ntuple(_ -> id(space(op, 1)), i - 1)..., op, ntuple(_ -> id(space(op, 1)), L - i - numout(op) + 1)...))
+
+@testset "density-matrix braiding, finite: $(sectortype(V))" for V in braid_spaces
+    L = 3
+    Random.seed!(4321)
+    ρd = randn(ComplexF64, V^L, V^L)
+    Od = randn(ComplexF64, V^L, V^L)
+    ρ, O = FiniteMPO(ρd), FiniteMPO(Od)
+    ψ = convert(FiniteMPS, ρ)
+    ref = tr(ρd' * Od * ρd)
+    nrm = tr(ρd' * ρd)
+
+    # both routes to <ρ|O|ρ> must reproduce the dense trace
+    @test dot(ρ, O * ρ) ≈ ref
+    @test dot(ψ, O, ψ) ≈ ref
+
+    # each local derivative operator must reproduce the same scalar
+    envs = environments(ψ, O, ψ)
+    for i in 1:L, prepare in (false, true)
+        h = AC_hamiltonian(i, ψ, O, ψ, envs; prepare)
+        @test dot(ψ.AC[i], h * ψ.AC[i]) ≈ ref
+    end
+    for i in 1:(L - 1)
+        x = AC2(ψ, i)
+        h = AC2_hamiltonian(i, ψ, O, ψ, envs; prepare = false)
+        @test dot(x, h * x) ≈ ref
+    end
+
+    # The prepared two-site derivative is a separate contraction whose crossings
+    # are encoded in `braid` levels rather than τ, along with fused legs
+    for i in 1:(L - 1)
+        x = AC2(ψ, i)
+        h = AC2_hamiltonian(i, ψ, O, ψ, envs; prepare = true)
+        @test dot(x, h * x) ≈ ref
+    end
+
+    # local observables measured on a density matrix
+    o1 = randn(ComplexF64, V, V)
+    o1 = o1 + o1'
+    @test expectation_value(ρ, 2 => o1) ≈ tr(ρd' * _embed(L, o1, 2) * ρd) / nrm
+    o2 = randn(ComplexF64, V ⊗ V, V ⊗ V)
+    o2 = o2 + o2'
+    for i in 1:(L - 1)
+        @test expectation_value(ρ, (i, i + 1) => o2) ≈ tr(ρd' * _embed(L, o2, i) * ρd) / nrm
+    end
+end
+
+# overlap per unit cell of two infinite MPOs as leading eigenvalue of
+# TM that contracts both physical legs directly
+# contains no braiding tensor -> independent reference for the MPS-MPO-MPS sandwich
+function mpo_overlap(mpo1, mpo2)
+    N = length(mpo1)
+    function step(v)
+        for i in 1:N
+            @plansor w[-1; -2] := v[1; 2] * conj(mpo1[i][1 3; 4 -1]) * mpo2[i][2 3; 4 -2]
+            v = w
+        end
+        return v
+    end
+    v0 = randn(ComplexF64, space(mpo1[N], 4)' ← space(mpo2[N], 4)')
+    vals, = eigsolve(step, v0, 1, :LM; krylovdim = 30)
+    return vals[1]
+end
+
+@testset "density-matrix braiding, infinite: $(sectortype(V))" for V in braid_spaces
+    Vv = V isa ComplexSpace ? ℂ^3 : V ⊕ V
+    for N in (1, 2)
+        Random.seed!(2024)
+        ρ = InfiniteMPO([randn(ComplexF64, Vv ⊗ V ← V ⊗ Vv) for _ in 1:N])
+        O = InfiniteMPO([randn(ComplexF64, V ⊗ V ← V ⊗ V) for _ in 1:N])
+        ref = mpo_overlap(ρ, O * ρ) / mpo_overlap(ρ, ρ)
+        ψ = convert(InfiniteMPS, ρ)
+        @test dot(ψ, O, ψ) ≈ ref
+    end
+end

From e6d874a1ba65ecaa13ecdc8763e59cf706e7ef30 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 19:58:01 +0200
Subject: [PATCH 05/11] add to changelog

---
 docs/src/changelog.md | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 893a670e7..953e85610 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -114,6 +114,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Fix hardcoding of number of physical spaces in the `changebonds` implementations for
   `FiniteMPS`, enabling its use for systems with composite physical spaces
   ([#514](https://github.com/QuantumKitHub/MPSKit.jl/pull/514))
+- A density matrix represented as a 3-leg MPO-as-MPS (`convert(FiniteMPS, ::FiniteMPO)`) braids its
+  ancilla leg past the operator's virtual leg in several places, and for anyonic sectors that crossing has a real handedness. For symmetric braiding, getting it backwards was silently harmless there, but several call sites had it
+  backwards, silently giving incorrect anyonic results: the density-matrix MPO transfer (`transfer_left`/
+  `transfer_right`), the two-site density-matrix observable (`contract_mpo_expval2`), and both the
+  unprepared and prepared `AC`/`AC2` derivative operators. This affected `approximate`, `DMRG2`, and
+  `timestep` on density matrices for anyonic sector types.
+- `_mpo_to_mps` threw a `SpaceMismatch` when converting a block/Jordan-structured MPO tensor (e.g.
+  from `make_time_mpo`) into an MPS, since it braided the physical leg directly against a
+  `SumSpace` virtual leg. It now densifies the tensor first, which the function already did to its
+  return value regardless, so this adds no extra cost.
 
 ### Performance
 

From a6471f7b8b491579705d5726385761217b061326 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 12 Aug 2026 09:53:38 +0200
Subject: [PATCH 06/11] import stuff

---
 test/algorithms/anyonic_braiding.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl
index ef80a9db6..7743cf1ca 100644
--- a/test/algorithms/anyonic_braiding.jl
+++ b/test/algorithms/anyonic_braiding.jl
@@ -8,9 +8,8 @@ using .TestSetup
 using Test, TestExtras
 using MPSKit
 using BlockTensorKit
-using MPSKit: AC_hamiltonian, AC2_hamiltonian, AC2
+using MPSKit: AC_hamiltonian, AC2_hamiltonian, AC2, eigsolve
 using TensorKit
-using KrylovKit: eigsolve
 using LinearAlgebra
 using Random
 

From 7338d894c7c70bdd170aa9677a1bd720f2e8066f Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 14 Aug 2026 16:45:23 +0200
Subject: [PATCH 07/11] invert some braids and clarify choice in convert

---
 src/algorithms/derivatives/mpo_derivatives.jl | 22 ++++++++++---------
 src/algorithms/expval.jl                      |  6 ++---
 src/operators/mpo.jl                          |  7 +++++-
 src/transfermatrix/transfer.jl                |  6 ++---
 4 files changed, 24 insertions(+), 17 deletions(-)

diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 12b7c845e..88b837d91 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -121,14 +121,14 @@ function (h::MPO_AC_Hamiltonian{<:MPSBondTensor, Nothing, <:MPSBondTensor})(x::M
     end
     return y
 end
-# braid needs to be inverted due to handedness of mpo.jl#L90
+# braid handedness dependent on MPO-to-MPS conversion
 function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})(
         x::GenericMPSTensor{<:Any, 3}
     )
     backend, allocator = h.backend, h.allocator
     @plansor backend = backend allocator = allocator begin
         y[-1 -2 -3; -4] ≔ h.leftenv[-1 7; 6] * x[6 4 2; 1] *
-            h.operators[1][7 -2; 4 5] * τ'[5 -3; 2 3] * h.rightenv[1 3; -4]
+            h.operators[1][7 -2; 4 5] * τ[5 -3; 2 3] * h.rightenv[1 3; -4]
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
@@ -152,15 +152,15 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
-# these braids need to be inverted due to handedness of mpo.jl#L90
+# braid handedness dependent on MPO-to-MPS conversion
 function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})(
         x::AbstractTensorMap{<:Any, <:Any, 3, 3}
     )
     backend, allocator = h.backend, h.allocator
     @plansor backend = backend allocator = allocator begin
         y[-1 -2 -3; -4 -5 -6] ≔ h.leftenv[-1 11; 10] * x[10 8 6; 1 2 4] *
-            h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ'[9 -3; 6 7] *
-            h.operators[2][7 -6; 4 5] * τ'[5 -5; 2 3]
+            h.rightenv[1 3; -4] * h.operators[1][11 -2; 8 9] * τ[9 -3; 6 7] *
+            h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3]
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
@@ -292,9 +292,9 @@ end
 function (H::PrecomputedACDerivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 1})
     backend, allocator = H.backend, H.allocator
     L, R = H.leftenv, H.rightenv
-    # braid needs to be inverted due to handedness of mpo.jl#L90
+    # braid handedness dependent on MPO-to-MPS conversion
     @plansor backend = backend allocator = allocator begin
-        xR[-1 -2; -4 -5 -3] := x[-1 -2 3; 1] * R[1 2; -4] * τ'[2 3; -5 -3]
+        xR[-1 -2; -4 -5 -3] := x[-1 -2 3; 1] * R[1 2; -4] * τ[2 3; -5 -3]
     end
     xR_fused = fuse_legs(xR, 2, 1)
     @plansor backend = backend allocator = allocator begin
@@ -302,20 +302,22 @@ function (H::PrecomputedACDerivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 1})
     end
     return y
 end
+
+# braid levels dependent on MPO-to-MPS conversion, see mpo.jl
 function (H::PrecomputedAC2Derivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 3})
     backend, allocator = H.backend, H.allocator
     L, R = H.leftenv, H.rightenv
 
-    x_braided = fuse_legs(braid(x, ((1, 2, 3, 5), (4, 6)), (1, 2, 3, 4, 6, 5)), 1, 2)
+    x_braided = fuse_legs(braid(x, ((1, 2, 3, 5), (4, 6)), (1, 2, 3, 4, 5, 6)), 1, 2)
 
     @plansor backend = backend allocator = allocator begin
         xR[-1 -2; -4 -6 -7 -5 -3] := x_braided[-1 -2 -3 -5; 1] * R[1 -7; -4 -6]
     end
-    @notensor xR_braided = braid(fuse_legs(xR, 2, 1), ((1, 4), (2, 3, 5, 6)), (1, 2, 3, 5, 6, 7))
+    @notensor xR_braided = braid(fuse_legs(xR, 2, 1), ((1, 4), (2, 3, 5, 6)), (1, 4, 6, 7, 5, 3))
     @plansor backend = backend allocator = allocator begin
         y_braided[-1 -2; -4 -6 -5 -3] := L[-1 -2; 1 2] * xR_braided[1 2; -4 -6 -5 -3]
     end
-    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 5, 6, 3))
+    return braid(y_braided, ((1, 2, 6), (3, 5, 4)), (1, 2, 4, 6, 5, 3))
 end
 
 const _ToPrepare = Union{
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index 90674f712..2b38e75dc 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -123,15 +123,15 @@ function contract_mpo_expval2(
     return @plansor conj(A1bar[4 5; 6]) * conj(A2bar[6 7; 8]) * O[5 7; 2 3] * A1[4 2; 1] *
         A2[1 3; 8]
 end
-# the ancilla braid need to be inverted due to handedness of mpo.jl#L90
+# the ancilla braid needs to be inverted due to handedness of converting MPO to MPS
 function contract_mpo_expval2(
         A1::GenericMPSTensor{S, 3}, A2::GenericMPSTensor{S, 3},
         O::AbstractTensorMap{<:Any, S},
         A1bar::GenericMPSTensor{S, 3} = A1, A2bar::GenericMPSTensor{S, 3} = A2
     ) where {S}
     numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator"))
-    return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ[9 6; 1 2] *
-        τ'[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7]
+    return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ'[9 6; 1 2] *
+        τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7]
 end
 
 function expectation_value(
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index 36a90a0d3..e3935c756 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -86,9 +86,14 @@ end
 function Base.convert(::Type{<:InfiniteMPS}, mpo::InfiniteMPO)
     return InfiniteMPS(map(_mpo_to_mps, parent(mpo)))
 end
+# bending O's own physical leg into an ancilla leg is a real crossing for anyonic sectors
+# overbraiding or underbraiding give different R-symbols, not just a relabeling
+# so it's a choice, and every other density-matrix-facing braid elsewhere
+# (transfer.jl, mpo_derivatives.jl, expval.jl) is defined relative to it
+# the current choice is to braid the ancilla leg over the MPO virtual leg
 function _mpo_to_mps(O::MPOTensor)
     O′ = O isa AbstractBlockTensorMap ? TensorMap(O) : O
-    @plansor A[-1 -2 -3; -4] := O′[-1 -2; 1 2] * τ[1 2; -4 -3]
+    @plansor A[-1 -2 -3; -4] := O′[-1 -2; 1 2] * τ'[1 2; -4 -3]
     return A isa AbstractBlockTensorMap ? TensorMap(A) : A
 end
 
diff --git a/src/transfermatrix/transfer.jl b/src/transfermatrix/transfer.jl
index 09b98d861..9cf61e363 100644
--- a/src/transfermatrix/transfer.jl
+++ b/src/transfermatrix/transfer.jl
@@ -147,16 +147,16 @@ function transfer_right(v::MPOTensor, O::MPOTensor, A::MPSTensor, Ab::MPSTensor)
 end
 
 # density matrix transfer
-# these braids need to be inverted due to handedness of mpo.jl#L90
+# braid handedness dependent on MPO-to-MPS conversion
 function transfer_left(
         x::MPSTensor, O::MPOTensor, A::GenericMPSTensor{<:Any, 3}, Ab::GenericMPSTensor{<:Any, 3}
     )
-    return @plansor y[-1 -2; -3] ≔ x[1 2; 6] * A[6 7 8; -3] * O[2 3; 7 5] * τ'[5 4; 8 -2] *
+    return @plansor y[-1 -2; -3] ≔ x[1 2; 6] * A[6 7 8; -3] * O[2 3; 7 5] * τ[5 4; 8 -2] *
         conj(Ab[1 3 4; -1])
 end
 function transfer_right(
         x::MPSTensor, O::MPOTensor, A::GenericMPSTensor{<:Any, 3}, Ab::GenericMPSTensor{<:Any, 3}
     )
-    return @plansor y[-1 -2; -3] ≔ A[-1 4 2; 1] * O[-2 6; 4 5] * τ'[5 7; 2 3] *
+    return @plansor y[-1 -2; -3] ≔ A[-1 4 2; 1] * O[-2 6; 4 5] * τ[5 7; 2 3] *
         conj(Ab[-3 6 7; 8]) * x[1 3; 8]
 end

From dbd8a0633e8ec3b0bf9ad63f7925c060b8a974e9 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 14 Aug 2026 16:46:14 +0200
Subject: [PATCH 08/11] make tests cheaper

---
 test/algorithms/anyonic_braiding.jl | 28 +++++++++++++---------------
 1 file changed, 13 insertions(+), 15 deletions(-)

diff --git a/test/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl
index 7743cf1ca..03e809a56 100644
--- a/test/algorithms/anyonic_braiding.jl
+++ b/test/algorithms/anyonic_braiding.jl
@@ -32,27 +32,27 @@ end
 
 @testset "expectation values with density matrices (Fibonacci)" begin
     V = Vect[FibonacciAnyon](:I => 1, :τ => 1)
-    L = 4
+    L = 3
     h = randn(ComplexF64, V ⊗ V, V ⊗ V)
     h = h + h'
     H = FiniteMPOHamiltonian(fill(V, L), [(i, i + 1) => h for i in 1:(L - 1)])
+
+    # free workaround: flatten jordan tensors into plain mpo tensors
+    # don't convert to tensormap, this is an exponential cost in L
+    H_dense = FiniteMPO([W isa AbstractBlockTensorMap ? TensorMap(W) : W for W in parent(H)])
     ρ = make_time_mpo(H, 0.1, TaylorCluster(; N = 2); imaginary_evolution = true)
-    ρd, Hd = convert(TensorMap, ρ), convert(TensorMap, H)
-    ref = tr(ρd' * Hd * ρd) / tr(ρd' * ρd)
+    ref = dot(ρ, H_dense * ρ) / norm(ρ)^2
 
     # expectation_value here attempts to braid ancilla leg with mpo virtual leg,
     # which is a sumspace/plain pair and thus fails
     @test_broken expectation_value(ρ, H) ≈ ref
 
-    # free workaround: flatten jordan tensors into plain mpo tensors
-    # don't convert to tensormap, this is an exponential cost in L
-    H_dense = FiniteMPO([W isa AbstractBlockTensorMap ? TensorMap(W) : W for W in parent(H)])
     @test expectation_value(ρ, H_dense) ≈ ref
 
     # these require the densification in _mpo_to_mps (mpo.jl#90)
+    Hd = convert(TensorMap, H)
     @test expectation_value(ρ, FiniteMPO(Hd)) ≈ ref
-    @test expectation_value(ρ, (1, 2) => h) + expectation_value(ρ, (2, 3) => h) +
-        expectation_value(ρ, (3, 4) => h) ≈ ref
+    @test expectation_value(ρ, (1, 2) => h) + expectation_value(ρ, (2, 3) => h) ≈ ref
 end
 
 @testset "braiding consistency in derivative operators (Fibonacci)" begin
@@ -94,10 +94,9 @@ end
 end
 
 const braid_spaces = (
-    ℂ^2,
-    Vect[Z2Irrep](0 => 1, 1 => 1),
-    Vect[FermionParity](0 => 1, 1 => 1),
-    Vect[FibonacciAnyon](:I => 1, :τ => 1),
+    ℂ^2, # bosonic
+    Vect[FermionParity](0 => 1, 1 => 1), # fermionic
+    Vect[FibonacciAnyon](:I => 1, :τ => 1), # anyonic
 )
 
 # build L-site operator that's the identity everywhere except at site i
@@ -144,9 +143,8 @@ _embed(L, op, i) = foldl(⊗, (ntuple(_ -> id(space(op, 1)), i - 1)..., op, ntup
     @test expectation_value(ρ, 2 => o1) ≈ tr(ρd' * _embed(L, o1, 2) * ρd) / nrm
     o2 = randn(ComplexF64, V ⊗ V, V ⊗ V)
     o2 = o2 + o2'
-    for i in 1:(L - 1)
-        @test expectation_value(ρ, (i, i + 1) => o2) ≈ tr(ρd' * _embed(L, o2, i) * ρd) / nrm
-    end
+    i = rand(1:(L - 1))
+    @test expectation_value(ρ, (i, i + 1) => o2) ≈ tr(ρd' * _embed(L, o2, i) * ρd) / nrm
 end
 
 # overlap per unit cell of two infinite MPOs as leading eigenvalue of

From d8bceaf1e54588780f93201e278b0b3c5efc27a9 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 18 Sep 2026 16:07:50 +0200
Subject: [PATCH 09/11] reduce yap

---
 docs/src/changelog.md | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 953e85610..87a5278dd 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -114,16 +114,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Fix hardcoding of number of physical spaces in the `changebonds` implementations for
   `FiniteMPS`, enabling its use for systems with composite physical spaces
   ([#514](https://github.com/QuantumKitHub/MPSKit.jl/pull/514))
-- A density matrix represented as a 3-leg MPO-as-MPS (`convert(FiniteMPS, ::FiniteMPO)`) braids its
-  ancilla leg past the operator's virtual leg in several places, and for anyonic sectors that crossing has a real handedness. For symmetric braiding, getting it backwards was silently harmless there, but several call sites had it
-  backwards, silently giving incorrect anyonic results: the density-matrix MPO transfer (`transfer_left`/
-  `transfer_right`), the two-site density-matrix observable (`contract_mpo_expval2`), and both the
-  unprepared and prepared `AC`/`AC2` derivative operators. This affected `approximate`, `DMRG2`, and
-  `timestep` on density matrices for anyonic sector types.
+- Fix density matrix methods with anyonic symmetries through consistency of braid orientation.
+  This affected `approximate`, `DMRG2`, and `timestep` on density matrices for anyonic sector types.
+  ([#509](https://github.com/QuantumKitHub/MPSKit.jl/pull/509))
 - `_mpo_to_mps` threw a `SpaceMismatch` when converting a block/Jordan-structured MPO tensor (e.g.
   from `make_time_mpo`) into an MPS, since it braided the physical leg directly against a
   `SumSpace` virtual leg. It now densifies the tensor first, which the function already did to its
-  return value regardless, so this adds no extra cost.
+  return value regardless, so this adds no extra cost. ([#509](https://github.com/QuantumKitHub/MPSKit.jl/pull/509))
 
 ### Performance
 

From 7eba218fb19d903e69f28914df3449cf45a2f60e Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 18 Sep 2026 20:08:08 +0200
Subject: [PATCH 10/11] fix broken test and use fast_tests

---
 test/algorithms/anyonic_braiding.jl | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/test/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl
index 03e809a56..22f59773b 100644
--- a/test/algorithms/anyonic_braiding.jl
+++ b/test/algorithms/anyonic_braiding.jl
@@ -43,10 +43,7 @@ end
     ρ = make_time_mpo(H, 0.1, TaylorCluster(; N = 2); imaginary_evolution = true)
     ref = dot(ρ, H_dense * ρ) / norm(ρ)^2
 
-    # expectation_value here attempts to braid ancilla leg with mpo virtual leg,
-    # which is a sumspace/plain pair and thus fails
-    @test_broken expectation_value(ρ, H) ≈ ref
-
+    @test expectation_value(ρ, H) ≈ ref
     @test expectation_value(ρ, H_dense) ≈ ref
 
     # these require the densification in _mpo_to_mps (mpo.jl#90)
@@ -93,11 +90,15 @@ end
     @test tdvp2_fid ≈ 1 atol = 1.0e-6 # fidelity was previously high, but not ≈ 1
 end
 
-const braid_spaces = (
-    ℂ^2, # bosonic
-    Vect[FermionParity](0 => 1, 1 => 1), # fermionic
-    Vect[FibonacciAnyon](:I => 1, :τ => 1), # anyonic
-)
+braid_spaces = if fast_tests
+    (Vect[FibonacciAnyon](:I => 1, :τ => 1),)
+else
+    (
+        ℂ^2, # bosonic
+        Vect[FermionParity](0 => 1, 1 => 1), # fermionic
+        Vect[FibonacciAnyon](:I => 1, :τ => 1), # anyonic
+    )
+end
 
 # build L-site operator that's the identity everywhere except at site i
 # where it's op with physical space V

From bfd030e8c18579089ef23909c32d86fbc2212bfa Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 18 Sep 2026 20:08:48 +0200
Subject: [PATCH 11/11] precompile fib

---
 LocalPreferences.toml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/LocalPreferences.toml b/LocalPreferences.toml
index cdb275e21..dcd3133f6 100644
--- a/LocalPreferences.toml
+++ b/LocalPreferences.toml
@@ -1,2 +1,5 @@
 [TensorOperations]
 precompile_workload = true
+
+[TensorKit]
+precompile_sectors = ["Trivial", "Z2Irrep", "SU2Irrep", "FermionParity", "FibonacciAnyon"]