diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 73da47899..d88b3aa3b 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -47,7 +47,10 @@ function changebonds(
         # Use the nullspaces and SVD decomposition to determine the optimal expansion space
         VL = left_null(ψ.AL[i])
         VR = right_null!(_transpose_tail(ψ.AR[i + 1]; copy = true))
-        intermediate = normalize!(adjoint(VL) * AC2 * adjoint(VR))
+        intermediate = adjoint(VL) * AC2 * adjoint(VR)
+        nrm = norm(intermediate)
+        # skip empty-content bonds; normalizing zero would NaN the SVD
+        nrm > eps(real(scalartype(intermediate)))^(3 / 4) && scale!(intermediate, inv(nrm))
         U, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd)
 
         AL′[i] = VL * U
@@ -72,7 +75,9 @@ function changebonds(ψ::MultilineMPS, H, alg::OptimalExpand, envs = environment
         # Use the nullspaces and SVD decomposition to determine the optimal expansion space
         VL = left_null(ψ.AL[i, j])
         VR = right_null!(_transpose_tail(ψ.AR[i, j + 1]; copy = true))
-        intermediate = normalize!(adjoint(VL) * AC2 * adjoint(VR))
+        intermediate = adjoint(VL) * AC2 * adjoint(VR)
+        nrm = norm(intermediate)
+        nrm > eps(real(scalartype(intermediate)))^(3 / 4) && scale!(intermediate, inv(nrm))
         U, _, Vᴴ = svd_trunc!(intermediate; trunc = alg.trscheme, alg = alg.alg_svd)
 
         AL′[i, j] = VL * U
@@ -99,7 +104,10 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op
 
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * AC2 * adjoint(NR)
-    _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+    nrm = norm(g2)
+    # nothing to expand here; normalizing a zero g2 would NaN the SVD
+    nrm ≤ eps(real(scalartype(g2)))^(3 / 4) && return ψ
+    _, _, Vᴴ = svd_trunc!(scale!(g2, inv(nrm)); trunc = alg.trscheme, alg = alg.alg_svd)
 
     # optimal vectors at site+1
     ar_re = Vᴴ * NR
@@ -126,7 +134,9 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Opt
 
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * AC2 * adjoint(NR)
-    U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+    nrm = norm(g2)
+    nrm ≤ eps(real(scalartype(g2)))^(3 / 4) && return ψ
+    U, _, _ = svd_trunc!(scale!(g2, inv(nrm)); trunc = alg.trscheme, alg = alg.alg_svd)
 
     # optimal vectors at site-1
     Q = NL * U
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index 7445a0e21..2ab840ffe 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -89,7 +89,10 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra
 
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * ac2 * adjoint(NR)
-    _, _, Vᴴ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+    nrm = norm(g2)
+    # nothing to expand here; normalizing a zero g2 would NaN the SVD
+    nrm ≤ eps(real(scalartype(g2)))^(3 / 4) && return ψ
+    _, _, Vᴴ = svd_trunc!(scale!(g2, inv(nrm)); trunc = alg.trscheme, alg = alg.alg_svd)
 
     # optimal vectors at site+1, zero weight at site
     ar_re = Vᴴ * NR
@@ -115,7 +118,9 @@ function changebond!(site::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::Ran
 
     # select the dominant directions in the complement of the current state
     g2 = adjoint(NL) * ac2 * adjoint(NR)
-    U, _, _ = svd_trunc!(normalize!(g2); trunc = alg.trscheme, alg = alg.alg_svd)
+    nrm = norm(g2)
+    nrm ≤ eps(real(scalartype(g2)))^(3 / 4) && return ψ
+    U, _, _ = svd_trunc!(scale!(g2, inv(nrm)); trunc = alg.trscheme, alg = alg.alg_svd)
 
     # optimal vectors at site-1, zero weight at site
     Q = NL * U
diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl
index a2aea1c08..da75e90b4 100644
--- a/test/algorithms/changebonds.jl
+++ b/test/algorithms/changebonds.jl
@@ -96,6 +96,32 @@ end
     @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5))
 end
 
+# Regression: CBE expanders must not crash on bonds with no expansion content. A charged
+# single-particle fermionic (fℤ₂-graded) state carries only one parity sector on the bonds away
+# from the particle, so the projected two-site update `adjoint(NL) * AC2 * adjoint(NR)` is
+# *structurally exactly zero* there. `OptimalExpand`/`RandExpand` used to `normalize!` that zero
+# tensor, producing NaNs that crashed the SVD (`ArgumentError: invalid argument #4 to LAPACK
+# call`); they must instead skip such bonds, leaving the state unchanged.
+@testset "CBE zero-content bonds (graded)" begin
+    L = 6
+    H = kitaev_model(ComplexF64, Trivial; t = 1.0, mu = 0, Delta = 0, L = L)
+    P = physicalspace(H)[1]
+    S = typeof(P)
+    Vtriv = oneunit(P)
+    Vodd = S(c => dim(P, c) for c in sectors(P) if !isone(c))
+    # single particle created on the last site: interior bonds carry only the even sector, so
+    # the expansion content on them is exactly zero
+    state = FiniteMPS(
+        [isometry(ComplexF64, (k <= L ? Vtriv : Vodd) ⊗ P, (k < L ? Vtriv : Vodd)) for k in 1:L]
+    )
+
+    state_oe, _ = changebonds(state, H, OptimalExpand(; trscheme = truncrank(2)))
+    @test abs(dot(state, state_oe)) ≈ 1 atol = 1.0e-8
+
+    state_re = changebonds(state, RandExpand(; trscheme = truncrank(2)))
+    @test abs(dot(state, state_re)) ≈ 1 atol = 1.0e-8
+end
+
 # density-matrix-style MPS: each site carries two physical legs (ket ⊗ bra). The operator-free
 # bond-change algorithms (`RandExpand` expansion, `SvdCut` truncation) must handle the extra
 # physical leg. Operator-based expanders (`OptimalExpand`/`SketchedExpand`) are not covered here