diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index b048a941..64f8c07d 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -55,6 +55,7 @@ include("inner.jl")
 include("normalize.jl")
 include("expect.jl")
 include("environment.jl")
+include("initialize_cache.jl")
 include("exports.jl")
 
 end
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index faddc4c9..4cee74b1 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -517,7 +517,7 @@ function inner_network(
     return BilinearFormNetwork(A, x, y; kwargs...)
 end
 
-norm_sqr_network(ψ::AbstractITensorNetwork) = inner_network(ψ, ψ)
+norm_sqr_network(ψ::AbstractITensorNetwork) = QuadraticFormNetwork(ψ)
 
 #
 # Printing
diff --git a/src/caches/abstractbeliefpropagationcache.jl b/src/caches/abstractbeliefpropagationcache.jl
index d682b2a7..29eb6542 100644
--- a/src/caches/abstractbeliefpropagationcache.jl
+++ b/src/caches/abstractbeliefpropagationcache.jl
@@ -1,7 +1,8 @@
 using Adapt: Adapt, adapt, adapt_structure
 using DataGraphs: DataGraphs, underlying_graph, vertex_data
+using Dictionaries: Dictionary
 using Graphs: Graphs, IsDirected, dst, src
-using ITensors: commoninds, delta, dir
+using ITensors: dir
 using LinearAlgebra: diag, dot
 using NDTensors: NDTensors
 using NamedGraphs.GraphsExtensions: subgraph
@@ -34,14 +35,6 @@ function message_diff(message_a::Vector{ITensor}, message_b::Vector{ITensor})
     return 1 - f
 end
 
-function default_message(datatype::Type{<:AbstractArray}, inds_e)
-    return [adapt(datatype, denseblocks(delta(i))) for i in inds_e]
-end
-
-function default_message(elt::Type{<:Number}, inds_e)
-    return default_message(Vector{elt}, inds_e)
-end
-default_messages(ptn::PartitionedGraph) = Dictionary()
 @traitfn default_bp_maxiter(g::::(!IsDirected)) = is_tree(g) ? 1 : nothing
 @traitfn function default_bp_maxiter(g::::IsDirected)
     return default_bp_maxiter(undirected_graph(underlying_graph(g)))
@@ -50,11 +43,6 @@ default_partitioned_vertices(ψ::AbstractITensorNetwork) = group(v -> v, vertice
 
 partitioned_tensornetwork(bpc::AbstractBeliefPropagationCache) = not_implemented()
 messages(bpc::AbstractBeliefPropagationCache) = not_implemented()
-function default_message(
-        bpc::AbstractBeliefPropagationCache, edge::QuotientEdge; kwargs...
-    )
-    return not_implemented()
-end
 default_update_alg(bpc::AbstractBeliefPropagationCache) = not_implemented()
 default_message_update_alg(bpc::AbstractBeliefPropagationCache) = not_implemented()
 Base.copy(bpc::AbstractBeliefPropagationCache) = not_implemented()
@@ -162,11 +150,6 @@ function PartitionedGraphs.quotientedge(
     return PartitionedGraphs.quotientedge(partitioned_tensornetwork(bpc), edge)
 end
 
-function linkinds(bpc::AbstractBeliefPropagationCache, pe::QuotientEdge)
-    pitn = partitioned_tensornetwork(bpc)
-    return commoninds(subgraph(pitn, src(pe)), subgraph(pitn, dst(pe)))
-end
-
 NDTensors.scalartype(bpc::AbstractBeliefPropagationCache) = scalartype(tensornetwork(bpc))
 
 """
@@ -187,12 +170,11 @@ function update_factor(bpc, vertex, factor)
     return bpc
 end
 
-function message(bpc::AbstractBeliefPropagationCache, edge::QuotientEdge; kwargs...)
-    mts = messages(bpc)
-    return get(() -> default_message(bpc, edge; kwargs...), mts, edge)
+function message(bpc::AbstractBeliefPropagationCache, edge::QuotientEdge)
+    return messages(bpc)[edge]
 end
-function messages(bpc::AbstractBeliefPropagationCache, edges; kwargs...)
-    return map(edge -> message(bpc, edge; kwargs...), edges)
+function messages(bpc::AbstractBeliefPropagationCache, edges)
+    return map(edge -> message(bpc, edge), edges)
 end
 function set_messages!(bpc::AbstractBeliefPropagationCache, quotientedges_messages)
     ms = messages(bpc)
diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl
index 0bb9ff93..6ee4eb8c 100644
--- a/src/caches/beliefpropagationcache.jl
+++ b/src/caches/beliefpropagationcache.jl
@@ -1,4 +1,5 @@
 using DataGraphs: DataGraphs, set_vertex_data!
+using Dictionaries: Dictionary
 using Graphs: IsDirected
 using ITensors: dir
 using LinearAlgebra: diag, dot
@@ -9,14 +10,6 @@ using NamedGraphs.PartitionedGraphs: AbstractPartitionedGraph, PartitionedGraph,
 using SimpleTraits: SimpleTraits, @traitfn, Not
 using SplitApplyCombine: group
 
-function default_cache_construction_kwargs(alg::Algorithm"bp", ψ::AbstractITensorNetwork)
-    return (; partitioned_vertices = default_partitioned_vertices(ψ))
-end
-
-function default_cache_construction_kwargs(alg::Algorithm"bp", pg::PartitionedGraph)
-    return (;)
-end
-
 struct BeliefPropagationCache{V, PV, PTN <: AbstractPartitionedGraph{V, PV}, MTS} <:
     AbstractBeliefPropagationCache{V, PV}
     partitioned_tensornetwork::PTN
@@ -24,7 +17,7 @@ struct BeliefPropagationCache{V, PV, PTN <: AbstractPartitionedGraph{V, PV}, MTS
 end
 
 #Constructors...
-function BeliefPropagationCache(ptn::PartitionedGraph; messages = default_messages(ptn))
+function BeliefPropagationCache(ptn::PartitionedGraph; messages = Dictionary())
     return BeliefPropagationCache(ptn, messages)
 end
 
@@ -41,20 +34,12 @@ function BeliefPropagationCache(
     return BeliefPropagationCache(tn, partitioned_vertices; kwargs...)
 end
 
-function cache(alg::Algorithm"bp", tn; kwargs...)
-    return BeliefPropagationCache(tn; kwargs...)
-end
-
 function partitioned_tensornetwork(bp_cache::BeliefPropagationCache)
     return bp_cache.partitioned_tensornetwork
 end
 
 messages(bp_cache::BeliefPropagationCache) = bp_cache.messages
 
-function default_message(bp_cache::BeliefPropagationCache, edge::QuotientEdge)
-    return default_message(datatype(bp_cache), linkinds(bp_cache, edge))
-end
-
 function Base.copy(bp_cache::BeliefPropagationCache)
     return BeliefPropagationCache(
         copy(partitioned_tensornetwork(bp_cache)), copy(messages(bp_cache))
diff --git a/src/contract.jl b/src/contract.jl
index 8c6f4f57..3592700e 100644
--- a/src/contract.jl
+++ b/src/contract.jl
@@ -42,12 +42,12 @@ function logscalar(
         alg::Algorithm,
         tn::AbstractITensorNetwork;
         (cache!) = nothing,
-        cache_construction_kwargs = default_cache_construction_kwargs(alg, tn),
+        cache_construction_kwargs = (;),
         update_cache = isnothing(cache!),
         cache_update_kwargs = (;)
     )
     if isnothing(cache!)
-        cache! = Ref(cache(alg, tn; cache_construction_kwargs...))
+        cache! = Ref(initialize_cache(alg, tn; cache_construction_kwargs...))
     end
 
     if update_cache
diff --git a/src/environment.jl b/src/environment.jl
index b7a0f6d4..e9a92bd8 100644
--- a/src/environment.jl
+++ b/src/environment.jl
@@ -24,11 +24,11 @@ function environment(
         vertices::Vector;
         (cache!) = nothing,
         update_cache = isnothing(cache!),
-        cache_construction_kwargs = default_cache_construction_kwargs(alg, ptn),
+        cache_construction_kwargs = (;),
         cache_update_kwargs = (;)
     )
     if isnothing(cache!)
-        cache! = Ref(cache(alg, ptn; cache_construction_kwargs...))
+        cache! = Ref(initialize_cache(alg, ptn; cache_construction_kwargs...))
     end
 
     if update_cache
diff --git a/src/expect.jl b/src/expect.jl
index 0a831b9e..2ad9e774 100644
--- a/src/expect.jl
+++ b/src/expect.jl
@@ -1,4 +1,3 @@
-using Dictionaries: Dictionary, set!
 using ITensors: Op, contract, op, which_op
 
 default_expect_alg() = "bp"
@@ -31,7 +30,7 @@ function expect(
     )
     ψIψ = QuadraticFormNetwork(ψ)
     if isnothing(cache!)
-        cache! = Ref(cache(alg, ψIψ; cache_construction_kwargs...))
+        cache! = Ref(initialize_cache(alg, ψIψ; cache_construction_kwargs...))
     end
 
     if update_cache
diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl
index 513a3d97..0404acf2 100644
--- a/src/formnetworks/bilinearformnetwork.jl
+++ b/src/formnetworks/bilinearformnetwork.jl
@@ -1,7 +1,7 @@
 using Adapt: adapt
 using DataGraphs: DataGraphs, set_vertex_data!
 using ITensors.NDTensors: datatype, denseblocks
-using ITensors: ITensor, Op, delta, prime, sim
+using ITensors: ITensor, Index, Op, dag, delta, prime, sim
 using NamedGraphs.GraphsExtensions: disjoint_union
 
 default_dual_site_index_map = prime
diff --git a/src/formnetworks/linearformnetwork.jl b/src/formnetworks/linearformnetwork.jl
index 654d858f..9e9fb4c4 100644
--- a/src/formnetworks/linearformnetwork.jl
+++ b/src/formnetworks/linearformnetwork.jl
@@ -1,5 +1,5 @@
 using DataGraphs: DataGraphs, set_vertex_data!
-using ITensors: ITensor, prime
+using ITensors: ITensor, dag, prime
 using NamedGraphs.GraphsExtensions: disjoint_union
 
 default_dual_link_index_map = prime
diff --git a/src/formnetworks/quadraticformnetwork.jl b/src/formnetworks/quadraticformnetwork.jl
index eb4b2674..adba9a3b 100644
--- a/src/formnetworks/quadraticformnetwork.jl
+++ b/src/formnetworks/quadraticformnetwork.jl
@@ -1,4 +1,7 @@
 using DataGraphs: DataGraphs, set_vertex_data!, underlying_graph, vertex_data
+using Dictionaries: Dictionary, set!
+using ITensors: ITensor, commoninds, dag, delta
+using NamedGraphs.PartitionedGraphs: PartitionedGraph, QuotientEdge, quotientedges
 
 default_index_map = prime
 default_inv_index_map = noprime
@@ -85,6 +88,41 @@ function QuadraticFormNetwork(
     return QuadraticFormNetwork(blf, dual_index_map, dual_inv_index_map)
 end
 
+# Build initial BP messages on each quotient edge as `delta(bra, ket)`
+# pairs, one per ket link Index crossing the cut. The bra-side counterpart
+# of each ket Index is computed explicitly via `dual_index_map(fn)`, so
+# the pairing is correct even when multiple link indices share an edge
+# (where `commoninds`-zip ordering between layers is not guaranteed).
+function identity_messages(
+        fn::QuadraticFormNetwork;
+        partitioned_vertices = default_partitioned_vertices(fn)
+    )
+    ptn = PartitionedGraph(fn, partitioned_vertices)
+    messages = Dictionary{QuotientEdge, Vector{ITensor}}()
+    tn = tensornetwork(fn)
+    elt = scalartype(tn)
+    map_idx = dual_index_map(fn)
+    pv = partitioned_vertices
+    ket_s = ket_vertex_suffix(fn)
+    for pe in quotientedges(ptn)
+        src_orig = unique(first.(filter(v -> last(v) == ket_s, pv[parent(src(pe))])))
+        dst_orig = unique(first.(filter(v -> last(v) == ket_s, pv[parent(dst(pe))])))
+        for (from_orig, to_orig, e) in (
+                (src_orig, dst_orig, pe),
+                (dst_orig, src_orig, reverse(pe)),
+            )
+            ms = ITensor[]
+            for v_from in from_orig, v_to in to_orig
+                for k in commoninds(tn[ket_vertex(fn, v_from)], tn[ket_vertex(fn, v_to)])
+                    push!(ms, delta(elt, dag(map_idx(k)), k))
+                end
+            end
+            set!(messages, e, ms)
+        end
+    end
+    return messages
+end
+
 function update(qf::QuadraticFormNetwork, original_state_vertex, ket_state::ITensor)
     state_inds = inds(ket_state)
     bra_state = replaceinds(dag(ket_state), state_inds, dual_index_map(qf).(state_inds))
diff --git a/src/initialize_cache.jl b/src/initialize_cache.jl
new file mode 100644
index 00000000..ca4bf71e
--- /dev/null
+++ b/src/initialize_cache.jl
@@ -0,0 +1,29 @@
+using Dictionaries: Dictionary
+using Graphs: is_tree
+using ITensors.NDTensors: @Algorithm_str, Algorithm
+using NamedGraphs.PartitionedGraphs: PartitionedGraph, quotient_graph
+
+# Build a cache for algorithm `alg` on `tn`. The fallback constructs a
+# plain `BeliefPropagationCache` with no message defaults; the
+# `QuadraticFormNetwork` specialization injects `identity_messages` on
+# loopy quotient graphs (canonical for the structurally ψ-vs-ψ case).
+function initialize_cache(alg::Algorithm"bp", tn; kwargs...)
+    return BeliefPropagationCache(tn; kwargs...)
+end
+
+function initialize_cache(
+        alg::Algorithm"bp",
+        fn::QuadraticFormNetwork;
+        partitioned_vertices = default_partitioned_vertices(fn),
+        messages = nothing
+    )
+    ptn = PartitionedGraph(fn, partitioned_vertices)
+    if isnothing(messages)
+        messages = if is_tree(quotient_graph(ptn))
+            Dictionary()
+        else
+            identity_messages(fn; partitioned_vertices)
+        end
+    end
+    return BeliefPropagationCache(ptn; messages)
+end
diff --git a/src/inner.jl b/src/inner.jl
index 3733de48..02760a7b 100644
--- a/src/inner.jl
+++ b/src/inner.jl
@@ -1,4 +1,4 @@
-using ITensors: inner, scalar
+using ITensors: inner, scalar, sim
 using LinearAlgebra: norm, norm_sqr
 
 default_contract_alg(tns::Tuple) = "bp"
@@ -173,7 +173,9 @@ end
 
 # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl,
 # or `norm_sqr` to match `LinearAlgebra.norm_sqr`
-LinearAlgebra.norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...)
+function LinearAlgebra.norm_sqr(ψ::AbstractITensorNetwork; kwargs...)
+    return scalar(norm_sqr_network(ψ); kwargs...)
+end
 
 function LinearAlgebra.norm(ψ::AbstractITensorNetwork; kwargs...)
     return sqrt(abs(real(norm_sqr(ψ; kwargs...))))
diff --git a/src/normalize.jl b/src/normalize.jl
index 83d77cec..d15b08cb 100644
--- a/src/normalize.jl
+++ b/src/normalize.jl
@@ -15,13 +15,13 @@ function rescale(
         tn::AbstractITensorNetwork,
         args...;
         (cache!) = nothing,
-        cache_construction_kwargs = default_cache_construction_kwargs(alg, tn),
+        cache_construction_kwargs = (;),
         update_cache = isnothing(cache!),
         cache_update_kwargs = (;),
         kwargs...
     )
     if isnothing(cache!)
-        cache! = Ref(cache(alg, tn; cache_construction_kwargs...))
+        cache! = Ref(initialize_cache(alg, tn; cache_construction_kwargs...))
     end
 
     if update_cache
@@ -55,7 +55,7 @@ end
 function LinearAlgebra.normalize(
         alg::Algorithm"exact", tn::AbstractITensorNetwork; kwargs...
     )
-    logn = logscalar(alg, inner_network(tn, tn); kwargs...)
+    logn = logscalar(alg, norm_sqr_network(tn); kwargs...)
     c = inv(exp(logn / (2 * length(vertices(tn)))))
     return map(t -> c * t, tn)
 end
@@ -64,17 +64,14 @@ function LinearAlgebra.normalize(
         alg::Algorithm,
         tn::AbstractITensorNetwork;
         (cache!) = nothing,
-        cache_construction_function = tn ->
-        cache(alg, tn; default_cache_construction_kwargs(alg, tn)...),
         update_cache = isnothing(cache!),
         cache_update_kwargs = (;),
         cache_construction_kwargs = (;)
     )
-    norm_tn = inner_network(tn, tn)
+    norm_tn = norm_sqr_network(tn)
     if isnothing(cache!)
-        cache! = Ref(cache(alg, norm_tn; cache_construction_kwargs...))
+        cache! = Ref(initialize_cache(alg, norm_tn; cache_construction_kwargs...))
     end
-
     vs = collect(vertices(tn))
     verts = vcat([ket_vertex(norm_tn, v) for v in vs], [bra_vertex(norm_tn, v) for v in vs])
     norm_tn = rescale(alg, norm_tn; verts, cache!, update_cache, cache_update_kwargs)
diff --git a/test/test_apply.jl b/test/test_apply.jl
index 2c928b03..21a97086 100644
--- a/test/test_apply.jl
+++ b/test/test_apply.jl
@@ -1,8 +1,8 @@
 using Compat: Compat
 using Graphs: vertices
-using ITensorNetworks:
-    BeliefPropagationCache, apply, environment, norm_sqr_network, siteinds, update
-using ITensors: ITensors, ITensor, inner, op
+using ITensorNetworks: BeliefPropagationCache, apply, environment, initialize_cache,
+    norm_sqr_network, siteinds, update
+using ITensors: ITensors, Algorithm, ITensor, inner, op
 using NamedGraphs.NamedGraphGenerators: named_grid
 using SplitApplyCombine: group
 using StableRNGs: StableRNG
@@ -18,18 +18,29 @@ include("utils.jl")
     ψ = random_tensornetwork(rng, s; link_space = χ)
     v1, v2 = (2, 2), (1, 2)
     ψψ = norm_sqr_network(ψ)
-    # Simple Belief Propagation grouping (one bra+ket per partition) gives
-    # a product environment around `[v1, v2]`, which is what `apply` requires.
-    bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ)))
-    bp_cache = update(bp_cache; maxiter = 20)
-    envsSBP = environment(bp_cache, [(v1, "bra"), (v1, "ket"), (v2, "bra"), (v2, "ket")])
+    # Vertices of `[v1, v2]` across all layers of `ψψ` (bra/ket/operator),
+    # so the environment around them is just the incoming BP messages —
+    # the per-site operator tensors aren't pulled in as central tensors.
+    env_verts(vs) = [
+        (v, suffix) for v in vs for suffix in ("bra", "ket", "operator")
+    ]
+    # Simple Belief Propagation grouping (one bra/ket/operator triple per
+    # partition) gives a product environment around `[v1, v2]`, which is
+    # what `apply` requires.
+    pv_SBP = group(v -> v[1], vertices(ψψ))
+    bp_cache = update(
+        initialize_cache(Algorithm("bp"), ψψ; partitioned_vertices = pv_SBP);
+        maxiter = 20
+    )
+    envsSBP = environment(bp_cache, env_verts((v1, v2)))
     # Column-grouping (one whole column per partition) gives a non-product
     # environment; `apply` should reject it.
-    bp_cache_col = BeliefPropagationCache(ψψ, group(v -> v[1][1], vertices(ψψ)))
-    bp_cache_col = update(bp_cache_col; maxiter = 20)
-    envsGBP = environment(
-        bp_cache_col, [(v1, "bra"), (v1, "ket"), (v2, "bra"), (v2, "ket")]
+    pv_col = group(v -> v[1][1], vertices(ψψ))
+    bp_cache_col = update(
+        initialize_cache(Algorithm("bp"), ψψ; partitioned_vertices = pv_col);
+        maxiter = 20
     )
+    envsGBP = environment(bp_cache_col, env_verts((v1, v2)))
     inner_alg = "exact"
     ngates = 5
     truncerr = 0.0
diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl
index 29aafc68..67d6fd1e 100644
--- a/test/test_belief_propagation.jl
+++ b/test/test_belief_propagation.jl
@@ -1,18 +1,16 @@
 using Compat: Compat
 using Graphs: vertices
-using ITensorNetworks: ITensorNetworks, BeliefPropagationCache, contract,
-    contraction_sequence, environment, inner_network, message, message_diff,
+using ITensorNetworks: ITensorNetworks, BeliefPropagationCache, QuadraticFormNetwork,
+    bra_vertex, contract, contraction_sequence, default_partitioned_vertices, environment,
+    identity_messages, ket_vertex, message, message_diff, norm_sqr_network, operator_vertex,
     partitioned_tensornetwork, scalar, siteinds, tensornetwork, update, update_factor,
-    updated_message, ⊗
+    updated_message
 include("utils.jl")
 using ITensors.NDTensors: array
-using ITensors: ITensors, Algorithm, ITensor, combiner, commoninds, dag, inds, inner, op,
-    prime, random_itensor
+using ITensors: ITensors, ITensor, combiner, dag, inds, inner, prime, random_itensor
 using LinearAlgebra: eigvals, tr
-using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
-using NamedGraphs.PartitionedGraphs: quotientedges
-using NamedGraphs: NamedEdge, NamedGraph
-using SplitApplyCombine: group
+using NamedGraphs.NamedGraphGenerators: named_grid
+using NamedGraphs.PartitionedGraphs: PartitionedGraph, quotientedges
 using StableRNGs: StableRNG
 using TensorOperations: TensorOperations
 using Test: @test, @testset
@@ -27,16 +25,22 @@ using Test: @test, @testset
         χ = 2
         rng = StableRNG(1234)
         ψ = random_tensornetwork(rng, elt, s; link_space = χ)
-        ψψ = ψ ⊗ prime(dag(ψ); sites = [])
-        bpc = BeliefPropagationCache(ψψ, group(v -> first(v), vertices(ψψ)))
+        ψψ = norm_sqr_network(ψ)
+        pv = default_partitioned_vertices(ψψ)
+        ptn = PartitionedGraph(ψψ, pv)
+        bpc = BeliefPropagationCache(
+            ptn; messages = identity_messages(ψψ; partitioned_vertices = pv)
+        )
 
-        #Test updating the tensors in the cache
-        vket, vbra = ((1, 1), 1), ((1, 1), 2)
+        # Test updating the tensors in the cache. QFN bra has both site
+        # and link inds primed (relative to the ket), so the bra-side
+        # tensor is constructed as `dag(prime(new_A))` directly.
+        v = (1, 1)
+        vket = ket_vertex(ψψ, v)
+        vbra = bra_vertex(ψψ, v)
         A = bpc[vket]
         new_A = random_itensor(elt, inds(A))
-        new_A_dag = ITensors.replaceind(
-            dag(prime(new_A)), only(s[first(vket)])', only(s[first(vket)])
-        )
+        new_A_dag = dag(prime(new_A))
         bpc[vket] = new_A
         bpc[vbra] = new_A_dag
         @test bpc[vket] == new_A
@@ -50,28 +54,31 @@ using Test: @test, @testset
             @test eltype(only(message(bpc, pe))) == elt
         end
         #Test updating the underlying tensornetwork in the cache
-        v = first(vertices(ψψ))
+        v_any = first(vertices(ψψ))
         rng = StableRNG(1234)
-        new_tensor = random_itensor(rng, inds(ψψ[v]))
-        bpc_updated = update_factor(bpc, v, new_tensor)
+        new_tensor = random_itensor(rng, inds(ψψ[v_any]))
+        bpc_updated = update_factor(bpc, v_any, new_tensor)
         ψψ_updated = tensornetwork(bpc_updated)
-        @test ψψ_updated[v] == new_tensor
+        @test ψψ_updated[v_any] == new_tensor
 
-        #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD
+        # Two-site RDM at `vs`: ask for the environment around all three
+        # layers (bra/ket/operator) at `vs` so the BP env is just incoming
+        # messages, then contract with only the bra and ket tensors at
+        # `vs` — dropping the per-site identity operator at those
+        # vertices leaves the bra (primed) and ket (unprimed) site inds
+        # open, which is exactly the RDM.
         vs = [(2, 2), (2, 3)]
-
-        # Prime the bra-ket shared site indices on the ket side at the
-        # selected vertices, so the contracted RDM has open primed/unprimed
-        # legs there. Mutates a copy of `ψψ` in place; no graph edits.
-        ψψsplit = copy(ψψ)
-        for v in vs
-            common = commoninds(ψψsplit[(v, 1)], ψψsplit[(v, 2)])
-            ψψsplit[(v, 2)] = prime(ψψsplit[(v, 2)], common)
-        end
-        env_tensors = environment(bpc, [(v, 2) for v in vs])
-        rdm =
-            contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]]))
-
+        env_vs = vcat(
+            [bra_vertex(ψψ, v) for v in vs],
+            [ket_vertex(ψψ, v) for v in vs],
+            [operator_vertex(ψψ, v) for v in vs]
+        )
+        env_tensors = environment(bpc, env_vs)
+        local_tensors = vcat(
+            ITensor[ψψ[bra_vertex(ψψ, v)] for v in vs],
+            ITensor[ψψ[ket_vertex(ψψ, v)] for v in vs]
+        )
+        rdm = contract(vcat(env_tensors, local_tensors))
         rdm = array(
             (rdm * combiner(inds(rdm; plev = 0)...)) * combiner(inds(rdm; plev = 1)...)
         )
@@ -92,3 +99,18 @@ using Test: @test, @testset
         @test iszero(scalar(ψ; alg = "bp"))
     end
 end
+
+# Regression guard for `identity_messages`: a tiny QN-graded loopy
+# network (4-cycle, S=1/2 sites) where the old single-leg `delta(i)`
+# initialization collapsed messages to empty blocks and BP produced
+# NaN. Two-leg `delta(b, k)` keeps the QN sectors aligned, so the
+# QFN-routed `scalar` (which auto-uses `identity_messages` because QFN
+# is structurally ψ-vs-ψ) must come back finite and nonzero.
+@testset "QN-graded loopy BP — identity_messages regression" begin
+    g = named_grid((2, 2); periodic = true)
+    s = siteinds("S=1/2", g; conserve_qns = true)
+    ψ = productstate(v -> isodd(sum(v)) ? "↑" : "↓", s)
+    n2 = scalar(QuadraticFormNetwork(ψ); alg = "bp", cache_update_kwargs = (; maxiter = 25))
+    @test isfinite(n2)
+    @test !iszero(n2)
+end
diff --git a/test/test_forms.jl b/test/test_forms.jl
index 3d71fdae..40923d45 100644
--- a/test/test_forms.jl
+++ b/test/test_forms.jl
@@ -59,15 +59,14 @@ include("utils.jl")
     ∂qf_∂v = only(environment(qf, state_vertices(qf, [v]); alg = "exact"))
     @test (∂qf_∂v) * (qf[ket_vertex(qf, v)] * qf[bra_vertex(qf, v)]) ≈ contract(qf)
 
-    ∂qf_∂v_bp = environment(qf, state_vertices(qf, [v]); alg = "bp", update_cache = false)
-    ∂qf_∂v_bp = contract(∂qf_∂v_bp)
-    ∂qf_∂v_bp /= norm(∂qf_∂v_bp)
-    ∂qf_∂v /= norm(∂qf_∂v)
-    @test ∂qf_∂v_bp != ∂qf_∂v
-
+    # `environment(::AbstractFormNetwork, …; alg = "bp")` builds messages
+    # via the form-network's `identity_messages` path and updates the
+    # cache; on a tree, BP is exact so the BP env matches the exact env
+    # after one sweep.
     ∂qf_∂v_bp = environment(qf, state_vertices(qf, [v]); alg = "bp", update_cache = true)
     ∂qf_∂v_bp = contract(∂qf_∂v_bp)
     ∂qf_∂v_bp /= norm(∂qf_∂v_bp)
+    ∂qf_∂v /= norm(∂qf_∂v)
     @test ∂qf_∂v_bp ≈ ∂qf_∂v
 
     #Test having non-uniform number of site indices per vertex
diff --git a/test/test_normalize.jl b/test/test_normalize.jl
index 9fa7e06b..47a478ac 100644
--- a/test/test_normalize.jl
+++ b/test/test_normalize.jl
@@ -1,9 +1,11 @@
 using Graphs: SimpleGraph, uniform_tree
-using ITensorNetworks: BeliefPropagationCache, QuadraticFormNetwork, edge_scalars, messages,
+using ITensorNetworks: BeliefPropagationCache, QuadraticFormNetwork,
+    default_partitioned_vertices, edge_scalars, identity_messages, messages,
     norm_sqr_network, rescale, scalartype, siteinds, vertex_scalars
 using ITensors: dag, inner, scalar
 using LinearAlgebra: normalize
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
+using NamedGraphs.PartitionedGraphs: PartitionedGraph
 using NamedGraphs: NamedGraph
 using StableRNGs: StableRNG
 using TensorOperations: TensorOperations
@@ -37,7 +39,18 @@ include("utils.jl")
     ψ = normalize(x; alg = "exact")
     @test scalar(norm_sqr_network(ψ); alg = "exact") ≈ 1.0
 
-    ψIψ_bpc = Ref(BeliefPropagationCache(QuadraticFormNetwork(x)))
+    # Loopy BP needs explicit initial messages; the form-network's
+    # `identity_messages` gives a properly-flowing pair-delta starting
+    # point even for QN-graded indices.
+    qfn = QuadraticFormNetwork(x)
+    pv = default_partitioned_vertices(qfn)
+    ptn = PartitionedGraph(qfn, pv)
+    ψIψ_bpc = Ref(
+        BeliefPropagationCache(
+            ptn;
+            messages = identity_messages(qfn; partitioned_vertices = pv)
+        )
+    )
     ψ = normalize(
         x; alg = "bp", (cache!) = ψIψ_bpc, update_cache = true,
         cache_update_kwargs = (; maxiter = 20)
diff --git a/test/utils.jl b/test/utils.jl
index e6492676..ddb54e3f 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -5,24 +5,14 @@
 
 using DataGraphs: underlying_graph, vertex_data
 using Dictionaries: Indices
-using Graphs: AbstractGraph, dst, edges, src, vertices
+using Graphs: AbstractGraph, add_edge!, dst, edges, src, vertices
 using ITensorNetworks: ITensorNetwork, IndsNetwork
 using ITensors.NDTensors: dim
-using ITensors: ITensors, ITensor, Index, QN, dag, hasqns, inds, itensor, onehot
+using ITensors: ITensors, ITensor, Index, QN, dag, hasqns, inds, itensor
 using NamedGraphs.GraphsExtensions: incident_edges
 using NamedGraphs: NamedGraph
 using Random: Random, AbstractRNG
 
-# Pick a length-1 space matching `tn`'s Index style: a dense `1` if no
-# vertex tensor has QN-graded indices, otherwise a single-block `[QN() => 1]`.
-# Used by `_add_edge!` below to thread a trivial Link Index between two
-# tensors without disturbing the QN structure (real `Index` constructors
-# need a space, not just a dimension).
-function _trivial_link_space(tn)
-    any(hasqns ∘ inds, (tn[v] for v in vertices(tn))) || return 1
-    return [QN() => 1]
-end
-
 # --- random_tensornetwork ----------------------------------------------------
 
 # Core: at each vertex of `graph`, place an `itensor(randn(rng, eltype, ...), inds_v)`
@@ -83,27 +73,14 @@ end
 
 # --- productstate -------------------------------------------------------------
 
-# Thread a fresh trivial-QN dim-1 link Index between `src(edge)` and `dst(edge)`
-# by outer-producting `onehot` basis vectors into each endpoint tensor. This is
-# the productstate-specific analog of `ITensorNetworks.add_edge!` — the latter
-# uses QR, which would push site-state QN flux into the link and leave BP's
-# `default_message` (single-index `delta`) with no compatible block. `onehot`
-# defaults to `Float64`, so we typecast it to `elt` to keep `productstate(elt,
-# ...)` element-type-preserving. Reverse-map reconciliation on the second
-# `setindex!` brings the graph edge along for free.
-function _add_edge!(elt::Type, tn, edge)
-    iₑ = Index(_trivial_link_space(tn), "Link")
-    X = ITensors.convert_eltype(elt, onehot(iₑ => 1))
-    tn[src(edge)] = tn[src(edge)] * X
-    tn[dst(edge)] = tn[dst(edge)] * dag(X)
-    return tn
-end
-
 # Build a product-state ITensorNetwork: start from a site-only TN (one tensor
 # per vertex with just the on-site state vector, looked up via
 # `ITensors.state(name, site_index)`), then add each edge from the original
-# IndsNetwork via `_add_edge!`. `state` is anything indexable by vertex (dict,
-# dictionary, array, ...); the `Function` method just materializes a Dict first.
+# IndsNetwork via `Graphs.add_edge!`. The latter threads a fresh link `Index`
+# via QR; QN flux on the link is fine here because BP messages now pair the
+# bra and ket legs explicitly via `identity_messages`. `state` is anything
+# indexable by vertex (dict, dictionary, array, ...); the `Function` method
+# just materializes a Dict first.
 function productstate(elt::Type, state::Function, s::IndsNetwork)
     return productstate(elt, Dict(v => state(v) for v in vertices(s)), s)
 end
@@ -115,7 +92,7 @@ function productstate(elt::Type, state, s::IndsNetwork)
     end
     tn = ITensorNetwork(ts)
     for e in edges(s)
-        _add_edge!(elt, tn, e)
+        add_edge!(tn, e)
     end
     return tn
 end