diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index 1c1895fc..0b7ca2ad 100644
--- a/docs/src/experimental_methods.md
+++ b/docs/src/experimental_methods.md
@@ -29,21 +29,6 @@ Methods which still need to be discussed, modified, or deprecated.
   linkinds(tn::AbstractITensorNetwork)
   ```
 
-* Rewrite every index of a network according to a structural mapping `IndsNetwork => IndsNetwork`
-  (site indices per vertex, link indices per edge). The two `IndsNetwork`s must share the
-  same underlying graph (`abstractitensornetwork.jl`):
-  ```julia
-  replaceinds(tn::AbstractITensorNetwork, is_is′::Pair{<:IndsNetwork, <:IndsNetwork})
-  ```
-
-* "Split" an edge index by applying a map to each copy of it on the adjacent ITensors.
-  By default the `dst(edge)` copy is primed and the `src(edge)` copy is unchanged (`abstractitensornetwork.jl`):
-  ```julia
-  split_index(tn::AbstractITensorNetwork, edges_to_split; 
-              src_ind_map::Function = identity,
-              dst_ind_map::Function = prime)
-  ```
-
 ## TreeTensorNetwork Types
 
 #### OpSum Constructors
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 83eacb51..ed98ed45 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -741,22 +741,6 @@ function combine_linkinds(
     return combine_linkinds(tn, combiners)
 end
 
-function split_index(
-        tn::AbstractITensorNetwork,
-        edges_to_split;
-        src_ind_map::Function = identity,
-        dst_ind_map::Function = prime
-    )
-    tn = copy(tn)
-    for e in edges_to_split
-        inds = commoninds(tn[src(e)], tn[dst(e)])
-        tn[src(e)] = replaceinds(tn[src(e)], inds, src_ind_map(inds))
-        tn[dst(e)] = replaceinds(tn[dst(e)], inds, dst_ind_map(inds))
-    end
-
-    return tn
-end
-
 function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...)
     return LinearFormNetwork(x, y; kwargs...)
 end
diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl
index ca22e679..8776c455 100644
--- a/test/test_belief_propagation.jl
+++ b/test/test_belief_propagation.jl
@@ -2,12 +2,12 @@ using Compat: Compat
 using Graphs: vertices
 using ITensorNetworks: ITensorNetworks, @preserve_graph, BeliefPropagationCache, contract,
     contraction_sequence, environment, inner_network, message, message_diff,
-    partitioned_tensornetwork, scalar, siteinds, split_index, tensornetwork, update,
-    update_factor, updated_message, ⊗
+    partitioned_tensornetwork, scalar, siteinds, tensornetwork, update, update_factor,
+    updated_message, ⊗
 include("utils.jl")
 using ITensors.NDTensors: array
-using ITensors:
-    ITensors, Algorithm, ITensor, combiner, dag, inds, inner, op, prime, random_itensor
+using ITensors: ITensors, Algorithm, ITensor, combiner, commoninds, dag, inds, inner, op,
+    prime, random_itensor
 using LinearAlgebra: eigvals, tr
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
 using NamedGraphs.PartitionedGraphs: quotientedges
@@ -60,7 +60,14 @@ using Test: @test, @testset
         #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD
         vs = [(2, 2), (2, 3)]
 
-        ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs]))
+        # 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]]))