From 55d1d75beaec0c5d948301067848c51314c5e7d8 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 16:17:22 -0400
Subject: [PATCH 01/10] Test out other message factorizations

Tries an alternative factorization of the norm-network messages and their storage convention.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 Project.toml                               |  2 +-
 docs/Project.toml                          |  2 +-
 examples/Project.toml                      |  2 +-
 src/apply/apply_operators.jl               | 64 +++++++++++++---------
 src/beliefpropagation/beliefpropagation.jl | 11 ++--
 src/beliefpropagation/messagecache.jl      | 15 ++---
 test/Project.toml                          |  2 +-
 test/test_apply_operator.jl                | 29 +++++++++-
 test/test_beliefpropagation.jl             | 39 +++++++++++--
 9 files changed, 115 insertions(+), 51 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9cc7de74..8879f95b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
-version = "0.10.2"
+version = "0.11.0"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/docs/Project.toml b/docs/Project.toml
index 6adb17f4..eee8bc4e 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,5 +10,5 @@ path = ".."
 [compat]
 Documenter = "1"
 ITensorFormatter = "0.2.27"
-ITensorNetworksNext = "0.10"
+ITensorNetworksNext = "0.11"
 Literate = "2"
diff --git a/examples/Project.toml b/examples/Project.toml
index 406118b8..0fa22c56 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,4 +5,4 @@ ITensorNetworksNext = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
 path = ".."
 
 [compat]
-ITensorNetworksNext = "0.10"
+ITensorNetworksNext = "0.11"
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index f824b9d7..f9e6849e 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -7,7 +7,7 @@ using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, op
 using LinearAlgebra: norm
 using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
 using NamedGraphs: boundary_edges
-using TensorAlgebra.MatrixAlgebra: invsqrth_safe, sqrth_safe
+using TensorAlgebra.MatrixAlgebra: pow_diag_safe
 
 # === Top-level user entry point ===
 
@@ -205,34 +205,43 @@ end
 
 # === BP simple-update implementation ===
 
-# The odd-parity sign leaves a fermionic message positive semidefinite in only one
-# bipartition, so diagonalize it in the transposed (bra, ket) one, placing the
-# eigenvectors on the bra and ket legs ready to sandwich a function of the eigenvalues.
-# Shared by `message_root` and `message_gauge`.
+# A power of a diagonal (eigenvalue or singular value) tensor acts entrywise on the diagonal,
+# with no bipartition to choose, unlike a matrix function of a general fermionic tensor.
+function pow_diag(d, p; kwargs...)
+    return ITB.nameddims(pow_diag_safe(ITB.unnamed(d), p; kwargs...), dimnames(d))
+end
+
+# A norm-network message is stored as the operator `bra ← ket`, the bipartition in which it
+# is positive semidefinite (see `similar_message_environment`). Diagonalize it there,
+# `m = v d v'`; `v` carries the bond leg under the ket leg's name with the bra leg's arrow, so
+# `conj(v)` is the factor whose bond leg contracts into the vertex the message flows into.
 function message_eigen(message)
-    ket, bra = outputnames(message), inputnames(message)
-    hermitian_message = project_hermitian(ITB.state(message), ket, bra)
-    d, v = eigh_full(hermitian_message, bra, ket)
-    name_d′, name_d = dimnames(d)
-    v_ket = conj(v)
-    v_bra = replacedimnames(v, only(ket) => only(bra), name_d => name_d′)
-    return d, v_bra, v_ket
+    bra, ket = outputnames(message), inputnames(message)
+    hermitian_message = project_hermitian(ITB.state(message), bra, ket)
+    return eigh_full(hermitian_message, bra, ket)
 end
 
-# The balanced Hermitian root `v * √d * v'` of the message. This gauge works with
-# fermions; the asymmetric gauge `√d * v'` has an issue that is under investigation.
+# The asymmetric square root `F = √d v'` of the message, `m = F' F`, absorbed into the vertex
+# the message flows into.
 function message_root(message)
-    d, v_bra, v_ket = message_eigen(message)
-    name_d′, name_d = dimnames(d)
-    return v_bra * sqrth_safe(d, (name_d′,), (name_d,)) * v_ket
+    d, v = message_eigen(message)
+    return message_root(d, v)
 end
+message_root(d, v) = pow_diag(d, 1 // 2) * conj(v)
 
-# The message root paired with its inverse `v * √d⁻¹ * v'`, to gauge a bond and undo it.
+# The root `F` paired with its inverse under contraction, `G = (F' F)⁻¹ F' = √m⁻¹ v`, which
+# undoes the gauge on the updated tensor. Under composition `G = v √d⁻¹`, but that product runs
+# over the eigenvalue leg and misses the odd-parity sign a graded contraction over the bond leg
+# carries when the receiving vertex holds the dual bond leg; `√m⁻¹ v` is formed over the bond
+# leg, so `G F` is the identity on the bond in both bond directions without a twist.
 function message_gauge(message)
-    d, v_bra, v_ket = message_eigen(message)
+    d, v = message_eigen(message)
+    bra, ket = only(outputnames(message)), only(inputnames(message))
     name_d′, name_d = dimnames(d)
-    return v_bra * sqrth_safe(d, (name_d′,), (name_d,)) * v_ket,
-        v_bra * invsqrth_safe(d, (name_d′,), (name_d,)) * v_ket
+    v_bra = replacedimnames(v, ket => bra, name_d => name_d′)
+    inv_root = v_bra * pow_diag(d, -1 // 2) * conj(v)
+    return message_root(d, v),
+        replacedimnames(inv_root * replacedimnames(v, name_d => name_d′), bra => ket)
 end
 
 function apply_gate_bp!(
@@ -281,8 +290,8 @@ function apply_gate_bp_nsite!(
     edges_in = boundary_edges(state, vs; dir = :in)
     siv_v1 = [message_gauge(env[e]) for e in edges_in if dst(e) == v1]
     siv_v2 = [message_gauge(env[e]) for e in edges_in if dst(e) == v2]
-    gauges_v1, inv_gauges_v1 = first.(siv_v1), conj.(last.(siv_v1))
-    gauges_v2, inv_gauges_v2 = first.(siv_v2), conj.(last.(siv_v2))
+    gauges_v1, inv_gauges_v1 = first.(siv_v1), last.(siv_v1)
+    gauges_v2, inv_gauges_v2 = first.(siv_v2), last.(siv_v2)
 
     ψ_v1 = prod([[state[v1]]; gauges_v1])
     ψ_v2 = prod([[state[v2]]; gauges_v2])
@@ -295,19 +304,20 @@ function apply_gate_bp_nsite!(
         S = S / norm(S)
     end
     name_v1, name_v2 = dimnames(S)
-    sqrt_S = sqrth_safe(S, (name_v1,), (name_v2,); atol = 0, rtol = 0)
+    sqrt_S = pow_diag(S, 1 // 2; atol = 0, rtol = 0)
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
     dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    # The graded contraction of a factor with its conjugate carries the odd-parity sign.
+    # Stored as `bra ← ket`, the bipartition in which the graded contraction of a factor with
+    # its conjugate is positive semidefinite (the same convention as `message_update!`).
     env[v1 => v2] = operator(
-        replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v1,), (name_v2,)
+        replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v2,), (name_v1,)
     )
     env[v2 => v1] = operator(
-        replacedimnames(conj(R_v2), name_v1 => name_v2) * R_v2, (name_v1,), (name_v2,)
+        replacedimnames(conj(R_v2), name_v1 => name_v2) * R_v2, (name_v2,), (name_v1,)
     )
     return dest
 end
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index e6a86bfd..5ce1fbce 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -270,14 +270,15 @@ function message_update!(algorithm::SimpleMessageUpdate, cache, factors, edge)
 end
 
 # `NormNetwork`: the message is a doubled (ket/bra) bond operator. Contracting a plain vertex factor
-# with the incoming messages leaves the surviving bond legs dangling, so assign the ket/bra pairing
-# the norm network gives this edge (the same convention as `similar_message_environment`). Normalize
-# by the trace, which is sign-correct on fermionic bonds where the entrywise `sum` can flip the
-# odd-parity block's sign.
+# with the incoming messages leaves the surviving bond legs dangling, so assign the bra/ket pairing
+# the norm network gives this edge (the same convention as `similar_message_environment`). In that
+# bipartition the message is positive semidefinite, so its trace is a positive normalization; the
+# entrywise `sum`, or the trace in the other bipartition, is the fermionic supertrace and can vanish
+# or flip the sign.
 function message_update!(algorithm::SimpleMessageUpdate, cache, factors::NormNetwork, edge)
     new_tensor = updated_message(algorithm, cache, factors, edge)
     new_message = operator(
-        new_tensor, linknames(KetView(factors), edge), linknames(BraView(factors), edge)
+        new_tensor, linknames(BraView(factors), edge), linknames(KetView(factors), edge)
     )
     if algorithm.normalize
         message_norm = tr(new_message)
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index fc6bc611..e954ac03 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -203,14 +203,15 @@ function similar_message_environment(nn::NormNetwork)
             ketview = KetView(nn)
 
             ketnames = linknames(ketview, edge)
-            ketaxis = unnamed.(linkaxes(ketview, edge))
-
             branames = linknames(braview, edge)
-
-            # Bond leg (ket) = operator output, bra-layer leg = input. Built on the src-side ket
-            # axis, whose arrow is opposite the dst endpoint's bond, so the gauge contracts back
-            # into the destination state.
-            message = similar_operator(ketview[vertex], ketaxis, ketnames, branames)
+            braaxis = unnamed.(linkaxes(braview, edge))
+
+            # Bra-layer leg = operator output, bond (ket) leg = input: the bipartition in which a
+            # norm-network message is positive semidefinite, so `one` is the trivial environment
+            # and `tr` is the trace rather than the fermionic supertrace. Built on the src-side
+            # bra axis; `similar_operator` flips the input side, which gives the ket leg the
+            # src-side ket arrow so the gauge contracts into the destination state.
+            message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
 
             return edge => message
         end
diff --git a/test/Project.toml b/test/Project.toml
index fdf4adab..b4a867f7 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -29,7 +29,7 @@ Dictionaries = "0.4.5"
 GradedArrays = "0.16"
 Graphs = "1.13.1"
 ITensorBase = "0.14"
-ITensorNetworksNext = "0.10"
+ITensorNetworksNext = "0.11"
 ITensorPkgSkeleton = "0.3.42"
 MatrixAlgebraKit = "0.6"
 NamedGraphs = "0.14"
diff --git a/test/test_apply_operator.jl b/test/test_apply_operator.jl
index 9bee98eb..830e340e 100644
--- a/test/test_apply_operator.jl
+++ b/test/test_apply_operator.jl
@@ -1,8 +1,9 @@
 using GradedArrays: U1, gradedrange
 using Graphs: dst, edges, src, vertices
-using ITensorBase: ITensorBase as ITB, Index, name, operator, setname, uniquename
+using ITensorBase: ITensorBase as ITB, Index, inputnames, name, operator, outputnames,
+    replacedimnames, setname, uniquename
 using ITensorNetworksNext: NormNetwork, apply_operator, apply_operators, insertlink!,
-    message_environment, tensornetwork
+    message_environment, message_gauge, tensornetwork
 using MatrixAlgebraKit: svd_trunc, truncrank
 using NamedGraphs: named_cycle_graph, named_path_graph
 using Random: AbstractRNG
@@ -91,4 +92,28 @@ end
         @test prod(gated) ≈ ITB.apply(g2, ITB.apply(g1, prod(network))) rtol =
             eps(real(T))^(1 / 3)
     end
+
+    @testset "message gauge is a gauge transformation" begin
+        rng = StableRNG(123)
+        g = named_path_graph(N)
+        site_axes = Dict(v => Index(site_range) for v in vertices(g))
+        network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
+        ψ = prod(network)
+
+        # Both directions of one bond, so the vertex receiving the message holds the dual bond
+        # leg in one case and the nondual one in the other.
+        for (edge, recv, send) in ((2 => 3, 3, 2), (3 => 2, 2, 3))
+            message = env[edge]
+            bra, ket = only(outputnames(message)), only(inputnames(message))
+            F, G = message_gauge(message)
+            # `F' F` reproduces the message, and `G F` is the identity on the bond: absorbing
+            # `F` into the receiver and `G` into the sender leaves the state unchanged.
+            @test replacedimnames(conj(F), ket => bra) * F ≈ ITB.state(message) rtol =
+                eps(real(T))^(1 / 3)
+            gauged = copy(network)
+            gauged[recv] = network[recv] * F
+            gauged[send] = network[send] * G
+            @test prod(gauged) ≈ ψ rtol = eps(real(T))^(1 / 3)
+        end
+    end
 end
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index f1c634b7..1f3a08ba 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -2,16 +2,18 @@ import AlgorithmsInterface as AI
 using Base.Broadcast: materialize
 using DataGraphs: DataGraphs, DataGraph, edge_data, edge_data_type
 using Dictionaries: Dictionary, dictionary, set!
-using GradedArrays: U1, gradedrange
+using GradedArrays: U1, gradedrange, isdual
 using Graphs: AbstractGraph, add_vertex!, dst, edges, has_edge, has_vertex, nv, rem_edge!,
     src, vertices
-using ITensorBase: Greedy, ITensor, Index, inds, name, noprime, outputnames, prime
+using ITensorBase: Greedy, ITensor, Index, inds, inputnames, name, noprime, outputnames,
+    prime, replacedimnames, state
 using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, MessageCache,
     NormNetwork, SimpleMessageUpdate, StopWhenConverged, beliefpropagation,
-    bethe_free_energy, contract_network, contraction_order, edge_scalar, factor_tensors,
-    incoming_messages, insertlink!, linkinds, message_environment, messagecache,
-    region_scalar, subgraph, tensornetwork, updated_message, vertex_scalar, vertex_scalars
-using LinearAlgebra: LinearAlgebra
+    bethe_free_energy, bratensor, contract_network, contraction_order, edge_scalar,
+    factor_tensors, incoming_messages, insertlink!, kettensor, linkaxes, linkinds,
+    message_environment, message_root, messagecache, region_scalar, subgraph, tensornetwork,
+    updated_message, vertex_scalar, vertex_scalars
+using LinearAlgebra: LinearAlgebra, norm, tr
 using NamedGraphs: NamedEdge, all_edges, incident_edges, named_comb_tree, named_grid,
     named_path_graph, vertextype
 using StableRNGs: StableRNG
@@ -286,6 +288,31 @@ end
             z_exact = (ket * conj(ket))[]
             z_bp = exp(bethe_free_energy(nn, cache))
             @test z_bp ≈ z_exact rtol = eps(real(T))^(1 / 3)
+
+            # Messages are stored as `bra ← ket`, the bipartition in which they are positive
+            # semidefinite: the square root `F` exists, reproduces the message as `F' F`, and the
+            # trace that normalizes them is positive.
+            for msg in edge_data(cache)
+                bra, ket = only(outputnames(msg)), only(inputnames(msg))
+                F = message_root(msg)
+                @test replacedimnames(conj(F), ket => bra) * F ≈ state(msg) rtol =
+                    eps(real(T))^(1 / 3)
+                @test real(tr(msg)) > 0
+            end
+
+            # In that bipartition `one` is the trivial environment: paired with the ket and bra
+            # layers of the rest of the network it gives that part's plain squared norm. The two
+            # edges have the receiving vertex on opposite ends of its bond's arrow.
+            @test isdual(only(linkaxes(network, 1 => 2))) !=
+                isdual(only(linkaxes(network, 4 => 3)))
+            ones = message_environment(one, nn)
+            for (edge, rest) in ((1 => 2, 2:4), (4 => 3, 1:3))
+                layers =
+                    [[kettensor(nn, v) for v in rest]; [bratensor(nn, v) for v in rest]]
+                z_rest = contract_network([state(ones[edge]); layers])[]
+                @test z_rest ≈ norm(prod([network[v] for v in rest]))^2 rtol =
+                    eps(real(T))^(1 / 3)
+            end
         end
     end
 

From 799aeabbb5cd6251c47427f2c65cb439a89e51f2 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 16:20:28 -0400
Subject: [PATCH 02/10] Trim comments

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/apply/apply_operators.jl               | 17 ++---------------
 src/beliefpropagation/beliefpropagation.jl |  6 ++----
 src/beliefpropagation/messagecache.jl      |  7 ++-----
 test/test_apply_operator.jl                |  4 ----
 test/test_beliefpropagation.jl             |  6 ------
 5 files changed, 6 insertions(+), 34 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index f9e6849e..a5a1ced2 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -205,35 +205,24 @@ end
 
 # === BP simple-update implementation ===
 
-# A power of a diagonal (eigenvalue or singular value) tensor acts entrywise on the diagonal,
-# with no bipartition to choose, unlike a matrix function of a general fermionic tensor.
 function pow_diag(d, p; kwargs...)
     return ITB.nameddims(pow_diag_safe(ITB.unnamed(d), p; kwargs...), dimnames(d))
 end
 
-# A norm-network message is stored as the operator `bra ← ket`, the bipartition in which it
-# is positive semidefinite (see `similar_message_environment`). Diagonalize it there,
-# `m = v d v'`; `v` carries the bond leg under the ket leg's name with the bra leg's arrow, so
-# `conj(v)` is the factor whose bond leg contracts into the vertex the message flows into.
 function message_eigen(message)
     bra, ket = outputnames(message), inputnames(message)
     hermitian_message = project_hermitian(ITB.state(message), bra, ket)
     return eigh_full(hermitian_message, bra, ket)
 end
 
-# The asymmetric square root `F = √d v'` of the message, `m = F' F`, absorbed into the vertex
-# the message flows into.
 function message_root(message)
     d, v = message_eigen(message)
     return message_root(d, v)
 end
 message_root(d, v) = pow_diag(d, 1 // 2) * conj(v)
 
-# The root `F` paired with its inverse under contraction, `G = (F' F)⁻¹ F' = √m⁻¹ v`, which
-# undoes the gauge on the updated tensor. Under composition `G = v √d⁻¹`, but that product runs
-# over the eigenvalue leg and misses the odd-parity sign a graded contraction over the bond leg
-# carries when the receiving vertex holds the dual bond leg; `√m⁻¹ v` is formed over the bond
-# leg, so `G F` is the identity on the bond in both bond directions without a twist.
+# The root `F = √d v'` and its inverse under contraction `G = √m⁻¹ v`, formed over the bond leg
+# so that `G F` is the identity on the bond for either bond direction.
 function message_gauge(message)
     d, v = message_eigen(message)
     bra, ket = only(outputnames(message)), only(inputnames(message))
@@ -311,8 +300,6 @@ function apply_gate_bp_nsite!(
     dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
     dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
 
-    # Stored as `bra ← ket`, the bipartition in which the graded contraction of a factor with
-    # its conjugate is positive semidefinite (the same convention as `message_update!`).
     env[v1 => v2] = operator(
         replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v2,), (name_v1,)
     )
diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index 5ce1fbce..3741fb3a 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -271,10 +271,8 @@ end
 
 # `NormNetwork`: the message is a doubled (ket/bra) bond operator. Contracting a plain vertex factor
 # with the incoming messages leaves the surviving bond legs dangling, so assign the bra/ket pairing
-# the norm network gives this edge (the same convention as `similar_message_environment`). In that
-# bipartition the message is positive semidefinite, so its trace is a positive normalization; the
-# entrywise `sum`, or the trace in the other bipartition, is the fermionic supertrace and can vanish
-# or flip the sign.
+# the norm network gives this edge (the same convention as `similar_message_environment`), in which
+# the message is positive semidefinite and its trace is a positive normalization.
 function message_update!(algorithm::SimpleMessageUpdate, cache, factors::NormNetwork, edge)
     new_tensor = updated_message(algorithm, cache, factors, edge)
     new_message = operator(
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index e954ac03..2e486bb2 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -206,11 +206,8 @@ function similar_message_environment(nn::NormNetwork)
             branames = linknames(braview, edge)
             braaxis = unnamed.(linkaxes(braview, edge))
 
-            # Bra-layer leg = operator output, bond (ket) leg = input: the bipartition in which a
-            # norm-network message is positive semidefinite, so `one` is the trivial environment
-            # and `tr` is the trace rather than the fermionic supertrace. Built on the src-side
-            # bra axis; `similar_operator` flips the input side, which gives the ket leg the
-            # src-side ket arrow so the gauge contracts into the destination state.
+            # Bra leg = operator output, ket leg = input, the bipartition in which the message
+            # is positive semidefinite.
             message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
 
             return edge => message
diff --git a/test/test_apply_operator.jl b/test/test_apply_operator.jl
index 830e340e..7cac6ea4 100644
--- a/test/test_apply_operator.jl
+++ b/test/test_apply_operator.jl
@@ -100,14 +100,10 @@ end
         network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
         ψ = prod(network)
 
-        # Both directions of one bond, so the vertex receiving the message holds the dual bond
-        # leg in one case and the nondual one in the other.
         for (edge, recv, send) in ((2 => 3, 3, 2), (3 => 2, 2, 3))
             message = env[edge]
             bra, ket = only(outputnames(message)), only(inputnames(message))
             F, G = message_gauge(message)
-            # `F' F` reproduces the message, and `G F` is the identity on the bond: absorbing
-            # `F` into the receiver and `G` into the sender leaves the state unchanged.
             @test replacedimnames(conj(F), ket => bra) * F ≈ ITB.state(message) rtol =
                 eps(real(T))^(1 / 3)
             gauged = copy(network)
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index 1f3a08ba..ccaec9a0 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -289,9 +289,6 @@ end
             z_bp = exp(bethe_free_energy(nn, cache))
             @test z_bp ≈ z_exact rtol = eps(real(T))^(1 / 3)
 
-            # Messages are stored as `bra ← ket`, the bipartition in which they are positive
-            # semidefinite: the square root `F` exists, reproduces the message as `F' F`, and the
-            # trace that normalizes them is positive.
             for msg in edge_data(cache)
                 bra, ket = only(outputnames(msg)), only(inputnames(msg))
                 F = message_root(msg)
@@ -300,9 +297,6 @@ end
                 @test real(tr(msg)) > 0
             end
 
-            # In that bipartition `one` is the trivial environment: paired with the ket and bra
-            # layers of the rest of the network it gives that part's plain squared norm. The two
-            # edges have the receiving vertex on opposite ends of its bond's arrow.
             @test isdual(only(linkaxes(network, 1 => 2))) !=
                 isdual(only(linkaxes(network, 4 => 3)))
             ones = message_environment(one, nn)

From 247a61157001c0f3697049c4580ba68660e45d17 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 16:27:53 -0400
Subject: [PATCH 03/10] Form the inverse gauge without the inverse root of the
 message

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index a5a1ced2..51cb0555 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -221,16 +221,13 @@ function message_root(message)
 end
 message_root(d, v) = pow_diag(d, 1 // 2) * conj(v)
 
-# The root `F = √d v'` and its inverse under contraction `G = √m⁻¹ v`, formed over the bond leg
-# so that `G F` is the identity on the bond for either bond direction.
+# The root `F = √d v'` and its inverse under contraction `G = v √d⁻¹ (v' v)`, where `v' v` is
+# contracted over the bond leg so that `G F` is the identity on the bond for either bond direction.
 function message_gauge(message)
     d, v = message_eigen(message)
-    bra, ket = only(outputnames(message)), only(inputnames(message))
     name_d′, name_d = dimnames(d)
-    v_bra = replacedimnames(v, ket => bra, name_d => name_d′)
-    inv_root = v_bra * pow_diag(d, -1 // 2) * conj(v)
-    return message_root(d, v),
-        replacedimnames(inv_root * replacedimnames(v, name_d => name_d′), bra => ket)
+    v′ = replacedimnames(v, name_d => name_d′)
+    return message_root(d, v), v′ * pow_diag(d, -1 // 2) * (conj(v) * v′)
 end
 
 function apply_gate_bp!(

From 90598bc8f2adcf4e96290e4417e3139b9d86acba Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 16:29:48 -0400
Subject: [PATCH 04/10] Reorder the inverse gauge

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 51cb0555..327da4d5 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -221,13 +221,14 @@ function message_root(message)
 end
 message_root(d, v) = pow_diag(d, 1 // 2) * conj(v)
 
-# The root `F = √d v'` and its inverse under contraction `G = v √d⁻¹ (v' v)`, where `v' v` is
+# The root `F = √d v'` and its inverse under contraction `G = v (v' v) √d⁻¹`, where `v' v` is
 # contracted over the bond leg so that `G F` is the identity on the bond for either bond direction.
 function message_gauge(message)
     d, v = message_eigen(message)
     name_d′, name_d = dimnames(d)
     v′ = replacedimnames(v, name_d => name_d′)
-    return message_root(d, v), v′ * pow_diag(d, -1 // 2) * (conj(v) * v′)
+    inv_root = v * (conj(v) * v′) * pow_diag(d, -1 // 2)
+    return message_root(d, v), replacedimnames(inv_root, name_d => name_d′)
 end
 
 function apply_gate_bp!(

From 1162afd2540d32251fe5755725efbdfb3e2ad8e0 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 17:36:38 -0400
Subject: [PATCH 05/10] Gauge with the Hermitian message roots via apply

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/apply/apply_operators.jl   | 65 +++++++++++-----------------------
 test/Project.toml              |  2 ++
 test/test_apply_operator.jl    | 32 +++++++----------
 test/test_beliefpropagation.jl | 18 +++++-----
 4 files changed, 46 insertions(+), 71 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 327da4d5..b94acf52 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -2,12 +2,11 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, src, vertices
-using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
-    outputnames, replacedimnames
+using ITensorBase: AbstractITensor, apply, dimnames, inputnames, operator, replacedimnames
 using LinearAlgebra: norm
-using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
+using MatrixAlgebraKit: project_hermitian, qr_compact, svd_trunc
 using NamedGraphs: boundary_edges
-using TensorAlgebra.MatrixAlgebra: pow_diag_safe
+using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe
 
 # === Top-level user entry point ===
 
@@ -205,31 +204,7 @@ end
 
 # === BP simple-update implementation ===
 
-function pow_diag(d, p; kwargs...)
-    return ITB.nameddims(pow_diag_safe(ITB.unnamed(d), p; kwargs...), dimnames(d))
-end
-
-function message_eigen(message)
-    bra, ket = outputnames(message), inputnames(message)
-    hermitian_message = project_hermitian(ITB.state(message), bra, ket)
-    return eigh_full(hermitian_message, bra, ket)
-end
-
-function message_root(message)
-    d, v = message_eigen(message)
-    return message_root(d, v)
-end
-message_root(d, v) = pow_diag(d, 1 // 2) * conj(v)
-
-# The root `F = √d v'` and its inverse under contraction `G = v (v' v) √d⁻¹`, where `v' v` is
-# contracted over the bond leg so that `G F` is the identity on the bond for either bond direction.
-function message_gauge(message)
-    d, v = message_eigen(message)
-    name_d′, name_d = dimnames(d)
-    v′ = replacedimnames(v, name_d => name_d′)
-    inv_root = v * (conj(v) * v′) * pow_diag(d, -1 // 2)
-    return message_root(d, v), replacedimnames(inv_root, name_d => name_d′)
-end
+apply_gauges(gauges, ψ) = foldl((ψ, gauge) -> apply(gauge, ψ), gauges; init = ψ)
 
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
@@ -256,13 +231,13 @@ function apply_gate_bp_nsite!(
         normalize, kwargs...
     )
     v = only(vs)
-    ψv = ITB.apply(op, state[v])
+    ψv = apply(op, state[v])
     if normalize
-        gauges = [
-            message_root(env[e])
-                for e in boundary_edges(state, vs; dir = :in)
+        sqrt_messages = [
+            sqrth_safe(project_hermitian(env[e])) for
+                e in boundary_edges(state, vs; dir = :in)
         ]
-        ψv /= norm(prod([[ψv]; gauges]))
+        ψv /= norm(apply_gauges(sqrt_messages, ψv))
     end
     dest[v] = ψv
     return dest
@@ -275,28 +250,30 @@ function apply_gate_bp_nsite!(
     )
     v1, v2 = vs
     edges_in = boundary_edges(state, vs; dir = :in)
-    siv_v1 = [message_gauge(env[e]) for e in edges_in if dst(e) == v1]
-    siv_v2 = [message_gauge(env[e]) for e in edges_in if dst(e) == v2]
-    gauges_v1, inv_gauges_v1 = first.(siv_v1), last.(siv_v1)
-    gauges_v2, inv_gauges_v2 = first.(siv_v2), last.(siv_v2)
+    roots_v1 =
+        [sqrth_invsqrth_safe(project_hermitian(env[e])) for e in edges_in if dst(e) == v1]
+    roots_v2 =
+        [sqrth_invsqrth_safe(project_hermitian(env[e])) for e in edges_in if dst(e) == v2]
+    sqrt_messages_v1, invsqrt_messages_v1 = first.(roots_v1), last.(roots_v1)
+    sqrt_messages_v2, invsqrt_messages_v2 = first.(roots_v2), last.(roots_v2)
 
-    ψ_v1 = prod([[state[v1]]; gauges_v1])
-    ψ_v2 = prod([[state[v2]]; gauges_v2])
+    ψ_v1 = apply_gauges(sqrt_messages_v1, state[v1])
+    ψ_v2 = apply_gauges(sqrt_messages_v2, state[v2])
 
     Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(dimnames(ψ_v1), dimnames(ψ_v2), dimnames(op)))
     Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(dimnames(ψ_v2), dimnames(ψ_v1), dimnames(op)))
-    op_R_v1v2 = ITB.apply(op, R_v1 * R_v2)
+    op_R_v1v2 = apply(op, R_v1 * R_v2)
     U_v1, S, U_v2 = svd_trunc(op_R_v1v2, setdiff(dimnames(R_v1), dimnames(R_v2)); trunc)
     if normalize
         S = S / norm(S)
     end
     name_v1, name_v2 = dimnames(S)
-    sqrt_S = pow_diag(S, 1 // 2; atol = 0, rtol = 0)
+    sqrt_S = sqrth_safe(S, (name_v1,), (name_v2,); atol = 0, rtol = 0)
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
-    dest[v1] = prod([[Q_v1 * R_v1]; inv_gauges_v1])
-    dest[v2] = prod([[Q_v2 * R_v2]; inv_gauges_v2])
+    dest[v1] = apply_gauges(invsqrt_messages_v1, Q_v1 * R_v1)
+    dest[v2] = apply_gauges(invsqrt_messages_v2, Q_v2 * R_v2)
 
     env[v1 => v2] = operator(
         replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v2,), (name_v1,)
diff --git a/test/Project.toml b/test/Project.toml
index b4a867f7..709cc250 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,6 +15,7 @@ OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
@@ -37,5 +38,6 @@ OMEinsumContractionOrders = "1"
 QuadGK = "2.11.2"
 Random = "1.10"
 StableRNGs = "1"
+TensorAlgebra = "0.20"
 TensorKitSectors = "0.3"
 Test = "1.10"
diff --git a/test/test_apply_operator.jl b/test/test_apply_operator.jl
index 7cac6ea4..8c2048a9 100644
--- a/test/test_apply_operator.jl
+++ b/test/test_apply_operator.jl
@@ -1,13 +1,13 @@
 using GradedArrays: U1, gradedrange
 using Graphs: dst, edges, src, vertices
-using ITensorBase: ITensorBase as ITB, Index, inputnames, name, operator, outputnames,
-    replacedimnames, setname, uniquename
+using ITensorBase: Index, apply, name, operator, setname, uniquename
 using ITensorNetworksNext: NormNetwork, apply_operator, apply_operators, insertlink!,
-    message_environment, message_gauge, tensornetwork
+    message_environment, tensornetwork
 using MatrixAlgebraKit: svd_trunc, truncrank
 using NamedGraphs: named_cycle_graph, named_path_graph
 using Random: AbstractRNG
 using StableRNGs: StableRNG
+using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe
 using TensorKitSectors: FermionParity
 using Test: @test, @testset
 
@@ -61,7 +61,7 @@ end
                 randn_operator(rng, T, (site_axes[2], site_axes[3])),
             )
             gated, _ = apply_operator(gate, network, env)
-            @test prod(gated) ≈ ITB.apply(gate, prod(network)) rtol = eps(real(T))^(1 / 3)
+            @test prod(gated) ≈ apply(gate, prod(network)) rtol = eps(real(T))^(1 / 3)
         end
     end
 
@@ -72,7 +72,7 @@ end
         network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
 
         gate = randn_operator(rng, T, (site_axes[2], site_axes[3]))
-        gated_full = ITB.apply(gate, prod(network))
+        gated_full = apply(gate, prod(network))
         left = [name(site_axes[v]) for v in 1:2]
         U, S, Vt = svd_trunc(gated_full, left; trunc = truncrank(k))
         gated, _ = apply_operator(gate, network, env; trunc = truncrank(k))
@@ -89,27 +89,21 @@ end
         g1 = randn_operator(rng, T, (site_axes[2], site_axes[3]))
         g2 = randn_operator(rng, T, (site_axes[3], site_axes[4]))
         gated, _ = apply_operators([g1, g2], network, env)
-        @test prod(gated) ≈ ITB.apply(g2, ITB.apply(g1, prod(network))) rtol =
+        @test prod(gated) ≈ apply(g2, apply(g1, prod(network))) rtol =
             eps(real(T))^(1 / 3)
     end
 
-    @testset "message gauge is a gauge transformation" begin
+    @testset "message roots gauge a vertex and undo it" begin
         rng = StableRNG(123)
         g = named_path_graph(N)
         site_axes = Dict(v => Index(site_range) for v in vertices(g))
         network, env = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4))
-        ψ = prod(network)
-
-        for (edge, recv, send) in ((2 => 3, 3, 2), (3 => 2, 2, 3))
-            message = env[edge]
-            bra, ket = only(outputnames(message)), only(inputnames(message))
-            F, G = message_gauge(message)
-            @test replacedimnames(conj(F), ket => bra) * F ≈ ITB.state(message) rtol =
-                eps(real(T))^(1 / 3)
-            gauged = copy(network)
-            gauged[recv] = network[recv] * F
-            gauged[send] = network[send] * G
-            @test prod(gauged) ≈ ψ rtol = eps(real(T))^(1 / 3)
+        rtol = eps(real(T))^(1 / 3)
+
+        for (edge, v) in ((2 => 3, 3), (3 => 2, 2))
+            sqrt_message, invsqrt_message = sqrth_invsqrth_safe(env[edge])
+            @test apply(invsqrt_message, apply(sqrt_message, network[v])) ≈ network[v] rtol =
+                rtol
         end
     end
 end
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index ccaec9a0..d759b51f 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -5,18 +5,19 @@ using Dictionaries: Dictionary, dictionary, set!
 using GradedArrays: U1, gradedrange, isdual
 using Graphs: AbstractGraph, add_vertex!, dst, edges, has_edge, has_vertex, nv, rem_edge!,
     src, vertices
-using ITensorBase: Greedy, ITensor, Index, inds, inputnames, name, noprime, outputnames,
-    prime, replacedimnames, state
+using ITensorBase:
+    Greedy, ITensor, Index, apply, inds, name, noprime, outputnames, prime, state
 using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, MessageCache,
     NormNetwork, SimpleMessageUpdate, StopWhenConverged, beliefpropagation,
     bethe_free_energy, bratensor, contract_network, contraction_order, edge_scalar,
     factor_tensors, incoming_messages, insertlink!, kettensor, linkaxes, linkinds,
-    message_environment, message_root, messagecache, region_scalar, subgraph, tensornetwork,
+    message_environment, messagecache, region_scalar, subgraph, tensornetwork,
     updated_message, vertex_scalar, vertex_scalars
 using LinearAlgebra: LinearAlgebra, norm, tr
 using NamedGraphs: NamedEdge, all_edges, incident_edges, named_comb_tree, named_grid,
     named_path_graph, vertextype
 using StableRNGs: StableRNG
+using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe
 using TensorKitSectors: FermionParity
 using Test: @test, @testset
 
@@ -289,12 +290,13 @@ end
             z_bp = exp(bethe_free_energy(nn, cache))
             @test z_bp ≈ z_exact rtol = eps(real(T))^(1 / 3)
 
-            for msg in edge_data(cache)
-                bra, ket = only(outputnames(msg)), only(inputnames(msg))
-                F = message_root(msg)
-                @test replacedimnames(conj(F), ket => bra) * F ≈ state(msg) rtol =
-                    eps(real(T))^(1 / 3)
+            for edge in edges(cache)
+                msg = cache[edge]
                 @test real(tr(msg)) > 0
+                sqrt_msg, invsqrt_msg = sqrth_invsqrth_safe(msg)
+                v = dst(edge)
+                @test apply(invsqrt_msg, apply(sqrt_msg, network[v])) ≈ network[v] rtol =
+                    eps(real(T))^(1 / 3)
             end
 
             @test isdual(only(linkaxes(network, 1 => 2))) !=

From 3f09fc9d3c2ebad39cd084e6b28938d25848339c Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 17:42:28 -0400
Subject: [PATCH 06/10] Inline the gauge folds

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/apply/apply_operators.jl | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index b94acf52..3d87b10e 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -204,8 +204,6 @@ end
 
 # === BP simple-update implementation ===
 
-apply_gauges(gauges, ψ) = foldl((ψ, gauge) -> apply(gauge, ψ), gauges; init = ψ)
-
 function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
         state::AbstractITensorNetwork, env; kwargs...
@@ -237,7 +235,7 @@ function apply_gate_bp_nsite!(
             sqrth_safe(project_hermitian(env[e])) for
                 e in boundary_edges(state, vs; dir = :in)
         ]
-        ψv /= norm(apply_gauges(sqrt_messages, ψv))
+        ψv /= norm(foldl((ψ, m) -> apply(m, ψ), sqrt_messages; init = ψv))
     end
     dest[v] = ψv
     return dest
@@ -257,8 +255,8 @@ function apply_gate_bp_nsite!(
     sqrt_messages_v1, invsqrt_messages_v1 = first.(roots_v1), last.(roots_v1)
     sqrt_messages_v2, invsqrt_messages_v2 = first.(roots_v2), last.(roots_v2)
 
-    ψ_v1 = apply_gauges(sqrt_messages_v1, state[v1])
-    ψ_v2 = apply_gauges(sqrt_messages_v2, state[v2])
+    ψ_v1 = foldl((ψ, m) -> apply(m, ψ), sqrt_messages_v1; init = state[v1])
+    ψ_v2 = foldl((ψ, m) -> apply(m, ψ), sqrt_messages_v2; init = state[v2])
 
     Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(dimnames(ψ_v1), dimnames(ψ_v2), dimnames(op)))
     Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(dimnames(ψ_v2), dimnames(ψ_v1), dimnames(op)))
@@ -272,8 +270,8 @@ function apply_gate_bp_nsite!(
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)
     R_v2 = sqrt_S * U_v2
 
-    dest[v1] = apply_gauges(invsqrt_messages_v1, Q_v1 * R_v1)
-    dest[v2] = apply_gauges(invsqrt_messages_v2, Q_v2 * R_v2)
+    dest[v1] = foldl((ψ, m) -> apply(m, ψ), invsqrt_messages_v1; init = Q_v1 * R_v1)
+    dest[v2] = foldl((ψ, m) -> apply(m, ψ), invsqrt_messages_v2; init = Q_v2 * R_v2)
 
     env[v1 => v2] = operator(
         replacedimnames(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v2,), (name_v1,)

From ca9c1f17c845309e6992350f186925d17ac1d21f Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 18:05:12 -0400
Subject: [PATCH 07/10] Add a free-fermion test on a tree and a cycle

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 test/test_free_fermions.jl | 123 +++++++++++++++++++++++++++++++++++++
 1 file changed, 123 insertions(+)
 create mode 100644 test/test_free_fermions.jl

diff --git a/test/test_free_fermions.jl b/test/test_free_fermions.jl
new file mode 100644
index 00000000..660f4e22
--- /dev/null
+++ b/test/test_free_fermions.jl
@@ -0,0 +1,123 @@
+using GradedArrays: gradedrange
+using Graphs: dst, edges, src, vertices
+using ITensorBase: Index, apply, name, operator, prime
+using ITensorNetworksNext:
+    NormNetwork, apply_operators, insertlink!, message_environment, tensornetwork
+using LinearAlgebra: I
+using NamedGraphs: named_comb_tree, named_cycle_graph
+using TensorAlgebra: TensorAlgebra as TA
+using TensorKitSectors: FermionParity, Z2Irrep
+using Test: @test, @testset
+
+# Free fermions hopping on a graph, compared against the single-particle correlation matrix.
+# Fermion signs change the densities on a graph with a branching vertex or a loop, so both
+# geometries are covered. The circuit is not truncated, so the comparison is exact.
+
+# Two-site matrices in the basis `|n₁ n₂⟩` with the first site fastest.
+basis(n1, n2) = 1 + n1 + 2n2
+const hopping = let h = zeros(4, 4)
+    h[basis(1, 0), basis(0, 1)] = h[basis(0, 1), basis(1, 0)] = 1
+    h
+end
+const pair_creation = let p = zeros(4, 4)
+    p[basis(1, 1), basis(0, 0)] = 1
+    p
+end
+# The current `i(c₁† c₂ - c₂† c₁)`. The symmetric hopping `c₁† c₂ + c₂† c₁` averages to zero on a
+# bipartite graph for this circuit, so the current is the edge observable that is checked.
+const current = let j = zeros(ComplexF64, 4, 4)
+    j[basis(1, 0), basis(0, 1)] = im
+    j[basis(0, 1), basis(1, 0)] = -im
+    j
+end
+const density = [0 0; 0 1]
+
+function one_site_operator(matrix, s::Index)
+    return operator(TA.project(matrix, (prime(s),), (s,)), (name(prime(s)),), (name(s),))
+end
+function two_site_operator(matrix, s1::Index, s2::Index)
+    codomain, domain = (prime(s1), prime(s2)), (s1, s2)
+    return operator(
+        TA.project(reshape(matrix, 2, 2, 2, 2), codomain, domain), name.(codomain),
+        name.(domain)
+    )
+end
+
+hopping_gate(θ) = exp(-im * θ * hopping)
+
+expectation(op, ψ) = (conj(ψ) * apply(op, ψ))[] / (conj(ψ) * ψ)[]
+
+# `⟨cᵢ† cⱼ⟩` after the hopping gates `(v1, v2, θ)`, starting from the sites in `occupied` filled.
+function reference_correlations(vs, occupied, gates)
+    index = Dict(v => i for (i, v) in enumerate(vs))
+    C = zeros(ComplexF64, length(vs), length(vs))
+    for v in occupied
+        C[index[v], index[v]] = 1
+    end
+    for (v1, v2, θ) in gates
+        u = Matrix{ComplexF64}(I, length(vs), length(vs))
+        ij = [index[v1], index[v2]]
+        u[ij, ij] = exp(-im * θ * [0 1; 1 0])
+        C = conj(u) * C * transpose(u)
+    end
+    return C
+end
+
+# Fill the pairs of neighboring sites in `pairs` from the vacuum, apply the hopping gates `hops`,
+# and return the densities and the currents of the resulting state, with site sectors of type `S`.
+function circuit_observables(S, g, pairs, hops)
+    site_axes = Dict(v => Index(gradedrange([S(0) => 1, S(1) => 1])) for v in vertices(g))
+    network = tensornetwork(vertices(g)) do v
+        return TA.project([1.0 + 0im, 0], (site_axes[v],))
+    end
+    for edge in edges(g)
+        insertlink!(network, edge)
+    end
+    env = message_environment(one, NormNetwork(network))
+    fill_gates = [
+        two_site_operator(pair_creation, site_axes[v1], site_axes[v2]) for
+            (v1, v2) in pairs
+    ]
+    hop_gates = [
+        two_site_operator(hopping_gate(θ), site_axes[v1], site_axes[v2]) for
+            (v1, v2, θ) in hops
+    ]
+    network, env = apply_operators([fill_gates; hop_gates], network, env)
+    ψ = prod(network)
+    densities = Dict(
+        v => expectation(one_site_operator(density, site_axes[v]), ψ) for v in vertices(g)
+    )
+    currents = Dict(
+        e => expectation(
+                two_site_operator(current, site_axes[src(e)], site_axes[dst(e)]),
+                ψ
+            )
+            for e in edges(g)
+    )
+    return densities, currents
+end
+
+@testset "free fermions ($label)" for (label, g, pairs) in (
+        ("comb tree", named_comb_tree((3, 2)), [((1, 1), (1, 2)), ((3, 1), (3, 2))]),
+        ("cycle", named_cycle_graph(6), [(1, 2), (4, 5)]),
+    )
+    vs = collect(vertices(g))
+    index = Dict(v => i for (i, v) in enumerate(vs))
+    hops = [(src(e), dst(e), 0.2 + 0.1 * k) for k in 1:3 for e in edges(g)]
+    C = reference_correlations(vs, [v for pair in pairs for v in pair], hops)
+
+    densities, currents = circuit_observables(FermionParity, g, pairs, hops)
+    for v in vs
+        @test densities[v] ≈ C[index[v], index[v]] atol = 1.0e-8
+    end
+    for e in edges(g)
+        i, j = index[src(e)], index[dst(e)]
+        @test currents[e] ≈ im * (C[i, j] - C[j, i]) atol = 1.0e-8
+    end
+
+    # Hardcore bosons (a bosonic ℤ₂ grading) go through the same circuit without the fermion
+    # signs and come out with different densities on these graphs, so the checks above are
+    # sensitive to the signs.
+    boson_densities, _ = circuit_observables(Z2Irrep, g, pairs, hops)
+    @test maximum(abs(boson_densities[v] - C[index[v], index[v]]) for v in vs) > 0.05
+end

From e7026e0217e3ce876f35422f275860eb8f4dce00 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 18:36:28 -0400
Subject: [PATCH 08/10] Bump the patch version instead of the minor

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 Project.toml          | 2 +-
 docs/Project.toml     | 2 +-
 examples/Project.toml | 2 +-
 test/Project.toml     | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8879f95b..61958c46 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
-version = "0.11.0"
+version = "0.10.3"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/docs/Project.toml b/docs/Project.toml
index eee8bc4e..6adb17f4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,5 +10,5 @@ path = ".."
 [compat]
 Documenter = "1"
 ITensorFormatter = "0.2.27"
-ITensorNetworksNext = "0.11"
+ITensorNetworksNext = "0.10"
 Literate = "2"
diff --git a/examples/Project.toml b/examples/Project.toml
index 0fa22c56..406118b8 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,4 +5,4 @@ ITensorNetworksNext = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
 path = ".."
 
 [compat]
-ITensorNetworksNext = "0.11"
+ITensorNetworksNext = "0.10"
diff --git a/test/Project.toml b/test/Project.toml
index 709cc250..51471fcf 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -30,7 +30,7 @@ Dictionaries = "0.4.5"
 GradedArrays = "0.16"
 Graphs = "1.13.1"
 ITensorBase = "0.14"
-ITensorNetworksNext = "0.11"
+ITensorNetworksNext = "0.10"
 ITensorPkgSkeleton = "0.3.42"
 MatrixAlgebraKit = "0.6"
 NamedGraphs = "0.14"

From 37a5e13fe1b3d14c952074b69391b9f4bc1a157f Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Sep 2026 19:08:46 -0400
Subject: [PATCH 09/10] Build the free-fermion test from fermion operators

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 test/test_free_fermions.jl | 178 ++++++++++++++++---------------------
 1 file changed, 78 insertions(+), 100 deletions(-)

diff --git a/test/test_free_fermions.jl b/test/test_free_fermions.jl
index 660f4e22..6a07aa5b 100644
--- a/test/test_free_fermions.jl
+++ b/test/test_free_fermions.jl
@@ -1,123 +1,101 @@
-using GradedArrays: gradedrange
+using GradedArrays: U1
 using Graphs: dst, edges, src, vertices
-using ITensorBase: Index, apply, name, operator, prime
-using ITensorNetworksNext:
-    NormNetwork, apply_operators, insertlink!, message_environment, tensornetwork
-using LinearAlgebra: I
+using ITensorBase: Index, apply, inds, operator, prime
+using ITensorNetworksNext: NormNetwork, apply_operator, apply_operators, insertlink!,
+    message_environment, tensornetwork
 using NamedGraphs: named_comb_tree, named_cycle_graph
-using TensorAlgebra: TensorAlgebra as TA
-using TensorKitSectors: FermionParity, Z2Irrep
+using TensorAlgebra: project, project_aux
+using TensorKitSectors: FermionNumber
 using Test: @test, @testset
 
-# Free fermions hopping on a graph, compared against the single-particle correlation matrix.
-# Fermion signs change the densities on a graph with a branching vertex or a loop, so both
-# geometries are covered. The circuit is not truncated, so the comparison is exact.
+# Free fermions on a tree with a branching vertex and on a loop, where the fermion signs change
+# the observables, evolved in imaginary time bond by bond and compared with the exact Slater
+# determinant. The circuit is not truncated, so the comparison is exact.
 
-# Two-site matrices in the basis `|n₁ n₂⟩` with the first site fastest.
-basis(n1, n2) = 1 + n1 + 2n2
-const hopping = let h = zeros(4, 4)
-    h[basis(1, 0), basis(0, 1)] = h[basis(0, 1), basis(1, 0)] = 1
-    h
-end
-const pair_creation = let p = zeros(4, 4)
-    p[basis(1, 1), basis(0, 0)] = 1
-    p
-end
-# The current `i(c₁† c₂ - c₂† c₁)`. The symmetric hopping `c₁† c₂ + c₂† c₁` averages to zero on a
-# bipartite graph for this circuit, so the current is the edge observable that is checked.
-const current = let j = zeros(ComplexF64, 4, 4)
-    j[basis(1, 0), basis(0, 1)] = im
-    j[basis(0, 1), basis(1, 0)] = -im
-    j
-end
-const density = [0 0; 0 1]
+const cdag_matrix = Float64[0 0; 1 0]
+const c_matrix = Float64[0 1; 0 0]
 
-function one_site_operator(matrix, s::Index)
-    return operator(TA.project(matrix, (prime(s),), (s,)), (name(prime(s)),), (name(s),))
+# `project_aux` derives an auxiliary leg carrying the charge of a charge-shifting operator, so
+# `c†` conserves charge with that leg left dangling, and `c†ᵢ cⱼ` pairs a `c†` and a `c` over one
+# shared such leg.
+function cdag(s::Index)
+    return operator(project_aux(cdag_matrix, (prime(s),), (s,)), [prime(s)], [s])
 end
-function two_site_operator(matrix, s1::Index, s2::Index)
-    codomain, domain = (prime(s1), prime(s2)), (s1, s2)
-    return operator(
-        TA.project(reshape(matrix, 2, 2, 2, 2), codomain, domain), name.(codomain),
-        name.(domain)
-    )
+function cdag_c(si::Index, sj::Index)
+    a = project_aux(cdag_matrix, (prime(si),), (si,))
+    aux = last(inds(a))
+    b = project(reshape(c_matrix, 2, 2, 1), (prime(sj),), (sj, aux))
+    return operator(a, [prime(si)], [si]) * operator(b, [prime(sj)], [sj])
 end
+hopping(si::Index, sj::Index) = cdag_c(si, sj) + cdag_c(sj, si)
+number(s::Index) = operator(project([0 0; 0 1], (prime(s),), (s,)), [prime(s)], [s])
 
-hopping_gate(θ) = exp(-im * θ * hopping)
+expectation(o, ψ) = (conj(ψ) * apply(o, ψ))[] / (conj(ψ) * ψ)[]
 
-expectation(op, ψ) = (conj(ψ) * apply(op, ψ))[] / (conj(ψ) * ψ)[]
+# The densities on the vertices followed by the hoppings on the edges, of a state or of a
+# correlation matrix `⟨cᵢ† cⱼ⟩`.
+function observables(g, sites, ψ)
+    return [
+        [expectation(number(sites[v]), ψ) for v in vertices(g)];
+        [expectation(hopping(sites[src(e)], sites[dst(e)]), ψ) for e in edges(g)]
+    ]
+end
+function observables(g, C::AbstractMatrix)
+    index = Dict(v => i for (i, v) in enumerate(vertices(g)))
+    return [
+        [C[index[v], index[v]] for v in vertices(g)];
+        [
+            C[index[src(e)], index[dst(e)]] + C[index[dst(e)], index[src(e)]] for
+                e in edges(g)
+        ]
+    ]
+end
 
-# `⟨cᵢ† cⱼ⟩` after the hopping gates `(v1, v2, θ)`, starting from the sites in `occupied` filled.
-function reference_correlations(vs, occupied, gates)
-    index = Dict(v => i for (i, v) in enumerate(vs))
-    C = zeros(ComplexF64, length(vs), length(vs))
-    for v in occupied
-        C[index[v], index[v]] = 1
+# Create a fermion on each site in `occupied`, then apply `exp(-τ h)` on the bonds `(v1, v2, τ)`.
+function evolve(sector, g, occupied, bonds)
+    sites = Dict(v => Index([sector(0) => 1, sector(1) => 1]) for v in vertices(g))
+    ψ = tensornetwork(v -> ones(sites[v]), vertices(g))
+    for e in edges(g)
+        insertlink!(ψ, e)
     end
-    for (v1, v2, θ) in gates
-        u = Matrix{ComplexF64}(I, length(vs), length(vs))
-        ij = [index[v1], index[v2]]
-        u[ij, ij] = exp(-im * θ * [0 1; 1 0])
-        C = conj(u) * C * transpose(u)
+    env = message_environment(one, NormNetwork(ψ))
+    for v in occupied
+        ψ, env = apply_operator(cdag(sites[v]), ψ, env)
     end
-    return C
+    gates = [exp(-τ * hopping(sites[v1], sites[v2])) for (v1, v2, τ) in bonds]
+    ψ, env = apply_operators(gates, ψ, env)
+    return prod(ψ), sites
 end
 
-# Fill the pairs of neighboring sites in `pairs` from the vacuum, apply the hopping gates `hops`,
-# and return the densities and the currents of the resulting state, with site sectors of type `S`.
-function circuit_observables(S, g, pairs, hops)
-    site_axes = Dict(v => Index(gradedrange([S(0) => 1, S(1) => 1])) for v in vertices(g))
-    network = tensornetwork(vertices(g)) do v
-        return TA.project([1.0 + 0im, 0], (site_axes[v],))
+# The state stays a Slater determinant with orbitals `Φ`. Each gate multiplies `Φ` by `exp(-τ h)`
+# on its two sites, and `⟨cᵢ† cⱼ⟩ = (Φ (Φᵀ Φ)⁻¹ Φᵀ)ᵢⱼ`.
+function correlations(g, occupied, bonds)
+    index = Dict(v => i for (i, v) in enumerate(vertices(g)))
+    Φ = zeros(length(index), length(occupied))
+    for (k, v) in enumerate(occupied)
+        Φ[index[v], k] = 1
     end
-    for edge in edges(g)
-        insertlink!(network, edge)
+    for (v1, v2, τ) in bonds
+        ij = [index[v1], index[v2]]
+        Φ[ij, :] = exp(-τ * [0 1; 1 0]) * Φ[ij, :]
     end
-    env = message_environment(one, NormNetwork(network))
-    fill_gates = [
-        two_site_operator(pair_creation, site_axes[v1], site_axes[v2]) for
-            (v1, v2) in pairs
-    ]
-    hop_gates = [
-        two_site_operator(hopping_gate(θ), site_axes[v1], site_axes[v2]) for
-            (v1, v2, θ) in hops
-    ]
-    network, env = apply_operators([fill_gates; hop_gates], network, env)
-    ψ = prod(network)
-    densities = Dict(
-        v => expectation(one_site_operator(density, site_axes[v]), ψ) for v in vertices(g)
-    )
-    currents = Dict(
-        e => expectation(
-                two_site_operator(current, site_axes[src(e)], site_axes[dst(e)]),
-                ψ
-            )
-            for e in edges(g)
-    )
-    return densities, currents
+    return Φ * ((Φ' * Φ) \ Φ')
 end
 
-@testset "free fermions ($label)" for (label, g, pairs) in (
-        ("comb tree", named_comb_tree((3, 2)), [((1, 1), (1, 2)), ((3, 1), (3, 2))]),
-        ("cycle", named_cycle_graph(6), [(1, 2), (4, 5)]),
+# An even number of fermions on the cycle, where an odd number would be indistinguishable from
+# hardcore bosons.
+@testset "free fermions ($label)" for (label, g, occupied) in (
+        ("comb tree", named_comb_tree((3, 2)), [(2, 1), (2, 2)]),
+        ("cycle", named_cycle_graph(6), [1, 2, 4, 5]),
     )
-    vs = collect(vertices(g))
-    index = Dict(v => i for (i, v) in enumerate(vs))
-    hops = [(src(e), dst(e), 0.2 + 0.1 * k) for k in 1:3 for e in edges(g)]
-    C = reference_correlations(vs, [v for pair in pairs for v in pair], hops)
+    bonds = [(src(e), dst(e), 0.5) for _ in 1:3 for e in edges(g)]
+    reference = observables(g, correlations(g, occupied, bonds))
 
-    densities, currents = circuit_observables(FermionParity, g, pairs, hops)
-    for v in vs
-        @test densities[v] ≈ C[index[v], index[v]] atol = 1.0e-8
-    end
-    for e in edges(g)
-        i, j = index[src(e)], index[dst(e)]
-        @test currents[e] ≈ im * (C[i, j] - C[j, i]) atol = 1.0e-8
-    end
+    ψ, sites = evolve(FermionNumber, g, occupied, bonds)
+    @test observables(g, sites, ψ) ≈ reference atol = 1.0e-8
 
-    # Hardcore bosons (a bosonic ℤ₂ grading) go through the same circuit without the fermion
-    # signs and come out with different densities on these graphs, so the checks above are
-    # sensitive to the signs.
-    boson_densities, _ = circuit_observables(Z2Irrep, g, pairs, hops)
-    @test maximum(abs(boson_densities[v] - C[index[v], index[v]]) for v in vs) > 0.05
+    # Hardcore bosons go through the same circuit without the fermion signs and come out
+    # different on these graphs, so the check above is sensitive to the signs.
+    ψ_boson, sites_boson = evolve(U1, g, occupied, bonds)
+    @test maximum(abs.(observables(g, sites_boson, ψ_boson) .- reference)) > 0.05
 end

From 3a64dfcc6f77ca7addd549ed919215cadae82425 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 23 Sep 2026 09:02:32 -0400
Subject: [PATCH 10/10] Name the correlator helpers in the free-fermion test

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 test/test_free_fermions.jl | 86 +++++++++++++++++++++-----------------
 1 file changed, 47 insertions(+), 39 deletions(-)

diff --git a/test/test_free_fermions.jl b/test/test_free_fermions.jl
index 6a07aa5b..c4450b67 100644
--- a/test/test_free_fermions.jl
+++ b/test/test_free_fermions.jl
@@ -32,27 +32,40 @@ number(s::Index) = operator(project([0 0; 0 1], (prime(s),), (s,)), [prime(s)],
 
 expectation(o, ψ) = (conj(ψ) * apply(o, ψ))[] / (conj(ψ) * ψ)[]
 
-# The densities on the vertices followed by the hoppings on the edges, of a state or of a
-# correlation matrix `⟨cᵢ† cⱼ⟩`.
-function observables(g, sites, ψ)
-    return [
-        [expectation(number(sites[v]), ψ) for v in vertices(g)];
-        [expectation(hopping(sites[src(e)], sites[dst(e)]), ψ) for e in edges(g)]
-    ]
+# Densities on the vertices and hoppings on the edges from the full state `ψ`.
+function full_state_correlators(g, sites, ψ)
+    densities = Dict(v => expectation(number(sites[v]), ψ) for v in vertices(g))
+    hoppings = Dict(
+        e => expectation(hopping(sites[src(e)], sites[dst(e)]), ψ) for e in edges(g)
+    )
+    return densities, hoppings
 end
-function observables(g, C::AbstractMatrix)
+
+# Densities and hoppings of free fermions created on the sites in `occupied` and evolved by
+# `exp(-τ h)` on the bonds `(v1, v2, τ)`. The state stays a Slater determinant with orbitals `Φ`,
+# each gate multiplies `Φ` by `exp(-τ h)` on its two sites, and `⟨cᵢ† cⱼ⟩ = (Φ (Φᵀ Φ)⁻¹ Φᵀ)ᵢⱼ`.
+function free_fermion_correlators(g, occupied, bonds)
     index = Dict(v => i for (i, v) in enumerate(vertices(g)))
-    return [
-        [C[index[v], index[v]] for v in vertices(g)];
-        [
-            C[index[src(e)], index[dst(e)]] + C[index[dst(e)], index[src(e)]] for
-                e in edges(g)
-        ]
-    ]
+    Φ = zeros(length(index), length(occupied))
+    for (k, v) in enumerate(occupied)
+        Φ[index[v], k] = 1
+    end
+    for (v1, v2, τ) in bonds
+        ij = [index[v1], index[v2]]
+        Φ[ij, :] = exp(-τ * [0 1; 1 0]) * Φ[ij, :]
+    end
+    C = Φ * ((Φ' * Φ) \ Φ')
+    densities = Dict(v => C[index[v], index[v]] for v in vertices(g))
+    hoppings = Dict(
+        e => C[index[src(e)], index[dst(e)]] + C[index[dst(e)], index[src(e)]] for
+            e in edges(g)
+    )
+    return densities, hoppings
 end
 
-# Create a fermion on each site in `occupied`, then apply `exp(-τ h)` on the bonds `(v1, v2, τ)`.
-function evolve(sector, g, occupied, bonds)
+# Create a fermion on each site in `occupied`, then apply `exp(-τ h)` on the bonds `(v1, v2, τ)`
+# by belief-propagation simple update. Returns the network and its site indices.
+function bp_evolve(sector, g, occupied, bonds)
     sites = Dict(v => Index([sector(0) => 1, sector(1) => 1]) for v in vertices(g))
     ψ = tensornetwork(v -> ones(sites[v]), vertices(g))
     for e in edges(g)
@@ -64,22 +77,7 @@ function evolve(sector, g, occupied, bonds)
     end
     gates = [exp(-τ * hopping(sites[v1], sites[v2])) for (v1, v2, τ) in bonds]
     ψ, env = apply_operators(gates, ψ, env)
-    return prod(ψ), sites
-end
-
-# The state stays a Slater determinant with orbitals `Φ`. Each gate multiplies `Φ` by `exp(-τ h)`
-# on its two sites, and `⟨cᵢ† cⱼ⟩ = (Φ (Φᵀ Φ)⁻¹ Φᵀ)ᵢⱼ`.
-function correlations(g, occupied, bonds)
-    index = Dict(v => i for (i, v) in enumerate(vertices(g)))
-    Φ = zeros(length(index), length(occupied))
-    for (k, v) in enumerate(occupied)
-        Φ[index[v], k] = 1
-    end
-    for (v1, v2, τ) in bonds
-        ij = [index[v1], index[v2]]
-        Φ[ij, :] = exp(-τ * [0 1; 1 0]) * Φ[ij, :]
-    end
-    return Φ * ((Φ' * Φ) \ Φ')
+    return ψ, sites
 end
 
 # An even number of fermions on the cycle, where an odd number would be indistinguishable from
@@ -89,13 +87,23 @@ end
         ("cycle", named_cycle_graph(6), [1, 2, 4, 5]),
     )
     bonds = [(src(e), dst(e), 0.5) for _ in 1:3 for e in edges(g)]
-    reference = observables(g, correlations(g, occupied, bonds))
+    densities, hoppings = free_fermion_correlators(g, occupied, bonds)
 
-    ψ, sites = evolve(FermionNumber, g, occupied, bonds)
-    @test observables(g, sites, ψ) ≈ reference atol = 1.0e-8
+    ψ, sites = bp_evolve(FermionNumber, g, occupied, bonds)
+    bp_densities, bp_hoppings = full_state_correlators(g, sites, prod(ψ))
+    for v in vertices(g)
+        @test bp_densities[v] ≈ densities[v] atol = 1.0e-8
+    end
+    for e in edges(g)
+        @test bp_hoppings[e] ≈ hoppings[e] atol = 1.0e-8
+    end
 
     # Hardcore bosons go through the same circuit without the fermion signs and come out
-    # different on these graphs, so the check above is sensitive to the signs.
-    ψ_boson, sites_boson = evolve(U1, g, occupied, bonds)
-    @test maximum(abs.(observables(g, sites_boson, ψ_boson) .- reference)) > 0.05
+    # different on these graphs, so the checks above are sensitive to the signs.
+    ψ_boson, sites_boson = bp_evolve(U1, g, occupied, bonds)
+    boson_densities, boson_hoppings = full_state_correlators(g, sites_boson, prod(ψ_boson))
+    @test max(
+        maximum(abs(boson_densities[v] - densities[v]) for v in vertices(g)),
+        maximum(abs(boson_hoppings[e] - hoppings[e]) for e in edges(g))
+    ) > 0.05
 end