diff --git a/Project.toml b/Project.toml
index 9e7d1ca..8cc9f16 100644
--- a/Project.toml
+++ b/Project.toml
@@ -41,7 +41,7 @@ SimpleGraphConverter = "0.1.0"
 SplitApplyCombine = "1.2.3"
 Statistics = "1.11.1"
 StatsBase = "0.34.4"
-TensorAlgebra = "0.17.5"
+TensorAlgebra = "0.17.6"
 TypeParameterAccessors = "0.3.10, 0.4"
 VectorInterface = "0.5"
 julia = "1.10"
diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl
index 05f15cd..b9ce458 100644
--- a/src/Apply/simple_update.jl
+++ b/src/Apply/simple_update.jl
@@ -37,12 +37,17 @@ function simple_update(
 
         # The environments are hermitian only up to numerical noise, so project before
         # the square roots (which require hermitian input).
-        sqrt_invsqrt = env -> sqrth_invsqrth_safe(
-            project_hermitian(env, (inds(env)[1],), (inds(env)[2],)),
-            (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0
-        )
-        sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1)
-        sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2)
+        sqrt_invsqrt = function (env, ψᵥ)
+            ket_ind, bra_ind = commonind(env, ψᵥ), uniqueind(env, ψᵥ)
+            # Factorize from bra to ket (codomain = bra), the bipartition in which the
+            # fermionic message is PSD.
+            return sqrth_invsqrth_safe(
+                project_hermitian(env, (ket_ind,), (bra_ind,)),
+                (bra_ind,), (ket_ind,); atol = sqrt_cutoff, rtol = 0
+            )
+        end
+        sqrt_inv_sqrt_envs_v1 = map(env -> sqrt_invsqrt(env, ψ⃗[1]), envs_v1)
+        sqrt_inv_sqrt_envs_v2 = map(env -> sqrt_invsqrt(env, ψ⃗[2]), envs_v2)
         sqrt_envs_v1, inv_sqrt_envs_v1 = first.(sqrt_inv_sqrt_envs_v1), last.(sqrt_inv_sqrt_envs_v1)
         sqrt_envs_v2, inv_sqrt_envs_v2 = first.(sqrt_inv_sqrt_envs_v2), last.(sqrt_inv_sqrt_envs_v2)
 
diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl
index dfb8f04..e16fcfc 100644
--- a/src/MessagePassing/abstractbeliefpropagationcache.jl
+++ b/src/MessagePassing/abstractbeliefpropagationcache.jl
@@ -180,7 +180,14 @@ function updated_message(
     updated_message = contract_network(contract_list; sequence)
 
     if alg.kwargs.normalize
-        message_norm = sum(updated_message)
+        # A doubled (ket/bra) network message is a bond operator: normalize by its trace,
+        # which is sign-correct for fermionic bonds (the entrywise `sum` can flip the sign).
+        # A single-layer network message is a vector with no bra/ket pairing, so use `sum`.
+        message_norm = if is_operator(updated_message)
+            tr(updated_message, operator_inds(updated_message)...)
+        else
+            sum(updated_message)
+        end
         if !iszero(message_norm)
             updated_message = updated_message / message_norm
         end
diff --git a/src/Ops.jl b/src/Ops.jl
index 1d2a67c..8f6e724 100644
--- a/src/Ops.jl
+++ b/src/Ops.jl
@@ -59,14 +59,17 @@ end
 
 # Embed a `d^n × d^n` operator matrix (computational basis, first site most
 # significant) into an `ITensor` with codomain `prime.(sites)` (outputs) and
-# domain `sites` (inputs). `project` is checked, so on graded sites an operator
-# that is not symmetric under the site index's grading throws an `InexactError`.
+# domain `sites` (inputs). Routed through `project_aux` (like `state`): a graded
+# operator that is odd under the site grading (e.g. `X`/`Y`, bare `c`/`c†`) has no
+# parity-block-diagonal `(prime.(sites), sites)` map, so it gets a trailing auxiliary
+# charge leg instead of throwing an `InexactError`. Even operators and dense (ungraded)
+# sites project cleanly, so they are unaffected (no aux leg).
 function _op_matrix_to_itensor(M::AbstractMatrix, sites::Tuple)
     ds = length.(sites)
     n = length(sites)
     A = reshape(Matrix{ComplexF64}(M), (reverse(ds)..., reverse(ds)...))
     A = permutedims(A, (reverse(1:n)..., reverse((n + 1):(2n))...))
-    return project(A, ITensorBase.prime.(sites), sites)
+    return project_aux(A, ITensorBase.prime.(sites), sites)
 end
 
 # Top-level `op(name, sites...; kwargs...)`. The identity is dimension-general (used
diff --git a/src/itensors.jl b/src/itensors.jl
index c31d5de..1b9c8da 100644
--- a/src/itensors.jl
+++ b/src/itensors.jl
@@ -6,16 +6,60 @@ using Adapt: Adapt
 using ITensorBase: ITensorBase, ITensor, Index, inds, noprime, plev, prime, space, unnamed
 using TensorAlgebra: TensorAlgebra, project, tryproject
 
-function project_aux(v::AbstractVector{<:Number}, i::Index)
-    length(v) == length(i) ||
-        error(
-        "state vector has dimension $(length(v)) but the site index has dimension $(length(i))"
+# Project `a` onto the symmetry-restricted space given by `codomain`/`domain`, mirroring
+# `TensorAlgebra.project`/`tryproject` (three-argument operator form; `a` is indexed positionally
+# as `(codomain..., domain...)`). When the plain projection is charge-forbidden — a parity-odd
+# state or operator — the residual charge is absorbed into a trailing auxiliary index instead of
+# throwing, so the returned ITensor carries that extra aux leg last.
+function project_aux(a::AbstractArray, codomain, domain)
+    return @something tryproject(a, codomain, domain) begin
+        projected_a = project(reshape(a, (size(a)..., 1)), space.(codomain), space.(domain))
+        aux = Index(TensorAlgebra.axes(projected_a, ndims(a) + 1))
+        ITensor(projected_a, (codomain..., domain..., aux))
+    end
+end
+project_aux(a::AbstractArray, codomain) = project_aux(a, codomain, ())
+project_aux(v::AbstractVector{<:Number}, i::Index) = project_aux(v, (i,))
+
+# Build two tensors that share one contractible auxiliary leg, so `t1 * t2` gives the fermion
+# string of e.g. `c†ᵢcⱼ` (a `c†` and a `c`) with no `flip`. The aux is minted by projecting `a1`
+# with a trailing dummy axis, which absorbs `a1`'s residual charge, then attaching it to `a2`.
+# This bare-axis core takes spaces (no names), so it could live beside `TensorAlgebra.project`.
+function project_pair(a1::AbstractArray, codomain1, domain1, a2::AbstractArray, codomain2, domain2)
+    p1 = project(reshape(a1, (size(a1)..., 1)), codomain1, domain1)
+    aux = TensorAlgebra.axes(p1, ndims(a1) + 1)
+    p2 = project(reshape(a2, (size(a2)..., 1)), codomain2, (domain2..., aux))
+    return p1, p2
+end
+
+# Named worker: strip to bare axes, call the core, and reattach names — the one extra step vs
+# `project`'s wrapper is minting a single fresh aux `Index` shared by both outputs, so `t1 * t2`
+# contracts them by name.
+function project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2)
+    p1, p2 = project_pair(
+        a1, space.(codomain1), space.(domain1),
+        a2, space.(codomain2), space.(domain2),
     )
-    ψ = tryproject(v, (i,))
-    isnothing(ψ) || return ψ
-    raw = project(reshape(v, (length(v), 1)), (space(i),), ())
-    aux = Index(TensorAlgebra.axes(raw, 2))
-    return ITensor(raw, (i, aux))
+    aux = Index(TensorAlgebra.axes(p1, ndims(a1) + 1))
+    t1 = ITensor(p1, (codomain1..., domain1..., aux))
+    t2 = ITensor(p2, (codomain2..., domain2..., aux))
+    return t1, t2
+end
+
+# Named layer (mirrors ITensorBase's named `project`): two entries read operand 1's flavor from
+# whichever of its sides is non-empty (an empty `codomain1` is the all-domain/state mirror of the
+# empty-domain form), so a named operand never falls through to the bare-axis core.
+function project_pair(
+        a1::AbstractArray, codomain1::Tuple{Index, Vararg{Index}}, domain1::Tuple{Vararg{Index}},
+        a2::AbstractArray, codomain2::Tuple{Vararg{Index}}, domain2::Tuple{Vararg{Index}}
+    )
+    return project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2)
+end
+function project_pair(
+        a1::AbstractArray, codomain1::Tuple{}, domain1::Tuple{Index, Vararg{Index}},
+        a2::AbstractArray, codomain2::Tuple{Vararg{Index}}, domain2::Tuple{Vararg{Index}}
+    )
+    return project_pair_itensor(a1, codomain1, domain1, a2, codomain2, domain2)
 end
 
 function onehot(eltype::Type, (i, p)::Pair{<:Index})
@@ -65,6 +109,14 @@ delta(is::Tuple) = delta(Float64, is)
 delta(is::Index...) = delta(Float64, is)
 delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 
+# Whether `a` is an operator tensor: its indices are exactly the plev-0 indices paired with
+# their primes. A plain state/vector has the plev-0 indices but not their primed partners.
+function is_operator(a::ITensor)
+    domain = filter(i -> plev(i) == 0, inds(a))
+    codomain = prime.(domain)
+    return issetequal(inds(a), [codomain; domain])
+end
+
 # The codomain/domain bipartition of an operator tensor: each plev-0 index paired with its
 # prime. Viewing the operator as this square map is what `tr` factors through.
 function operator_inds(a::ITensor)