diff --git a/Project.toml b/Project.toml
index 6c60958d..f0b72efe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "GradedArrays"
 uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
-version = "0.13.4"
+version = "0.13.5"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -42,6 +42,6 @@ SUNRepresentations = "0.3, 0.4"
 SparseArraysBase = "0.10"
 SplitApplyCombine = "1.2.3"
 StridedViews = "0.5"
-TensorAlgebra = "0.15"
+TensorAlgebra = "0.16"
 TensorKitSectors = "0.3"
 julia = "1.10"
diff --git a/src/fusion.jl b/src/fusion.jl
index 538e8698..a349f616 100644
--- a/src/fusion.jl
+++ b/src/fusion.jl
@@ -120,12 +120,14 @@ end
 
 # ========================  AbelianGradedArray unmatricize  ========================
 
+# `unmatricize` receives the domain axes codomain-facing (un-dualized); a graded array stores
+# them dualized, so `conj` re-dualizes them before they are placed.
 function TensorAlgebra.unmatricize(
         ::SectorFusion, m::AbstractSectorDelta,
         codomain_axes::Tuple{Vararg{SectorRange}},
         domain_axes::Tuple{Vararg{SectorRange}}
     )
-    return AbelianSectorDelta{eltype(m)}((codomain_axes..., domain_axes...))
+    return AbelianSectorDelta{eltype(m)}((codomain_axes..., conj.(domain_axes)...))
 end
 
 # Unmatricize a 2D sector array back to an N-D AbelianSectorArray. The
@@ -158,7 +160,9 @@ function TensorAlgebra.unmatricize(
     )
     K = length(codomain_axes)
     N = K + length(domain_axes)
-    a = similar(m, (codomain_axes..., domain_axes...))
+    # The domain axes arrive codomain-facing (un-dualized); dualize them into the stored
+    # convention with `conj` as they are flattened into the destination axes.
+    a = similar(m, (codomain_axes..., conj.(domain_axes)...))
     return TensorAlgebra.unmatricizeperm!(
         SectorFusion(), a, m, ntuple(identity, Val(K)), ntuple(i -> K + i, Val(N - K))
     )
diff --git a/test/Project.toml b/test/Project.toml
index a96e2b98..fd26809f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -32,7 +32,7 @@ SUNRepresentations = "0.3, 0.4"
 SafeTestsets = "0.1"
 SparseArraysBase = "0.10"
 Suppressor = "0.2.8"
-TensorAlgebra = "0.15"
+TensorAlgebra = "0.16"
 TensorKitSectors = "0.3"
 Test = "1.10"
 TestExtras = "0.3.1"
diff --git a/test/test_fermionic.jl b/test/test_fermionic.jl
index f6fd7d73..184ca005 100644
--- a/test/test_fermionic.jl
+++ b/test/test_fermionic.jl
@@ -426,7 +426,9 @@ end
         codomain_axes = ntuple(i -> axes(a_dest)[invperm_codomain[i]], K)
         domain_axes = ntuple(i -> axes(a_dest)[invperm_domain[i]], ndims(a_dest) - K)
         return permutedims(
-            unmatricize(m, codomain_axes, domain_axes),
+            # `axes(a_dest)` are stored (dualized) domain axes, but `unmatricize` takes them
+            # codomain-facing, so un-dualize with `conj` before the call.
+            unmatricize(m, codomain_axes, conj.(domain_axes)),
             invperm((invperm_codomain..., invperm_domain...))
         )
     end
diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index 52015940..c7ffdf07 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -240,8 +240,10 @@ end
     @test collect(keys(fsm.domain)) == [U1(0), U1(1)]
     @test collect(keys(fsm.blocks)) == [U1(0), U1(1)]
 
-    # Round-trip through `unmatricize` recovers the original blocks.
-    a_back = unmatricize(fsm, (r, r), (dual(r),))
+    # Round-trip through `unmatricize` recovers the original blocks. The domain axes are
+    # passed codomain-facing (un-dualized), so the original `dual(r)` domain axis is given
+    # as `r`.
+    a_back = unmatricize(fsm, (r, r), (r,))
     @test a_back isa AbelianGradedArray
     @test ndims(a_back) == 3
     for I in eachblockstoredindex(a)
@@ -368,8 +370,9 @@ end
     @test sm isa SectorMatrix
     @test ndims(sm) == 2
 
-    # Unmatricize back to 3D
-    s_back = unmatricize(sm, (codomain_ax,), (domain_ax1, domain_ax2))
+    # Unmatricize back to 3D. The domain axes are passed codomain-facing (un-dualized),
+    # so the stored `conj`-ed domain axes are given as their un-dualized counterparts.
+    s_back = unmatricize(sm, (codomain_ax,), (conj(domain_ax1), conj(domain_ax2)))
     @test s_back isa AbelianSectorArray
     @test ndims(s_back) == 3
     @test size(s_back) == size(s)