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"]
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 893a670e7..87a5278dd 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -114,6 +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))
+- 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. ([#509](https://github.com/QuantumKitHub/MPSKit.jl/pull/509))
 
 ### Performance
 
diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 8745146ab..88b837d91 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -121,6 +121,7 @@ function (h::MPO_AC_Hamiltonian{<:MPSBondTensor, Nothing, <:MPSBondTensor})(x::M
     end
     return y
 end
+# braid handedness dependent on MPO-to-MPS conversion
 function (h::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor})(
         x::GenericMPSTensor{<:Any, 3}
     )
@@ -151,6 +152,7 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end
+# braid handedness dependent on MPO-to-MPS conversion
 function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor})(
         x::AbstractTensorMap{<:Any, <:Any, 3, 3}
     )
@@ -290,7 +292,7 @@ end
 function (H::PrecomputedACDerivative)(x::AbstractTensorMap{<:Any, <:Any, 3, 1})
     backend, allocator = H.backend, H.allocator
     L, R = H.leftenv, H.rightenv
-
+    # 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]
     end
@@ -300,6 +302,8 @@ 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
@@ -313,7 +317,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{
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index 191d2f93c..2b38e75dc 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -123,13 +123,14 @@ 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 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] *
+    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
 
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index 225c83756..e3935c756 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -86,8 +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)
-    @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..9cf61e363 100644
--- a/src/transfermatrix/transfer.jl
+++ b/src/transfermatrix/transfer.jl
@@ -146,7 +146,8 @@ 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
+# braid handedness dependent on MPO-to-MPS conversion
 function transfer_left(
         x::MPSTensor, O::MPOTensor, A::GenericMPSTensor{<:Any, 3}, Ab::GenericMPSTensor{<:Any, 3}
     )
diff --git a/test/algorithms/anyonic_braiding.jl b/test/algorithms/anyonic_braiding.jl
new file mode 100644
index 000000000..22f59773b
--- /dev/null
+++ b/test/algorithms/anyonic_braiding.jl
@@ -0,0 +1,178 @@
+println("
+------------------------------------
+|   Anyonic braiding consistency    |
+------------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using BlockTensorKit
+using MPSKit: AC_hamiltonian, AC2_hamiltonian, AC2, eigsolve
+using TensorKit
+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 = 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)
+    ref = dot(ρ, H_dense * ρ) / norm(ρ)^2
+
+    @test expectation_value(ρ, H) ≈ ref
+    @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) ≈ 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
+
+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
+_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'
+    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
+# 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