diff --git a/Project.toml b/Project.toml
index f5c609f3c..5f61ffbc2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -37,7 +37,7 @@ MPSKitAdaptExt = "Adapt"
 [compat]
 Accessors = "0.1"
 Adapt = "4"
-BlockTensorKit = "0.3.15"
+BlockTensorKit = "0.3.16"
 Compat = "3.47, 4.10"
 DocStringExtensions = "0.9.3"
 HalfIntegers = "1.6.0"
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index be8bae968..edd1d8d71 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -97,6 +97,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   tensor, so it returned `O * O` on twice the physical space instead of `O`.
   ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484))
 - `leading_boundary` with `IDMRG2` didn't update the left edge of the AC tensor, which could result in space mismatches depending on the truncation scheme. This is corrected for in ([#516](https://github.com/QuantumKitHub/MPSKit.jl/pull/516)).
+- Fix hardcoding of number of physical spaces in the `changebonds` implementations for
+  `FiniteMPS`, enabling its use for systems with composite physical spaces
+  ([#514](https://github.com/QuantumKitHub/MPSKit.jl/pull/514))
 
 ### Performance
 
diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 0a039ddf3..a82060c01 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -179,7 +179,7 @@ function changebond!(
     right_tail = _transpose_tail(right)
     # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new
     # directions), leaving the state unchanged
-    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc_space = (codomain(right_tail)[1] ⊕ _lastspace(Q)') ← domain(right_tail)
     nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
     AL_exp = catdomain(left, Q)
 
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index 484deda5c..5cfff2de8 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -141,7 +141,7 @@ function changebond!(
     Q = NL * U
     # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new directions)
     right_tail = _transpose_tail(right)
-    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc_space = (codomain(right_tail)[1] ⊕ _lastspace(Q)') ← domain(right_tail)
     nc, Qr = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
     AL_exp = catdomain(left, Q)
 
diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl
index a0c43d9c5..4ba2eedeb 100644
--- a/src/algorithms/changebonds/sketchedexpand.jl
+++ b/src/algorithms/changebonds/sketchedexpand.jl
@@ -152,7 +152,7 @@ function changebond!(
     right_tail = _transpose_tail(right)
     # embed `_transpose_tail(right)` into the enlarged codomain (zero weight in the new
     # directions), leaving the state unchanged
-    nc_space = (codomain(right_tail)[1] ⊕ space(Q, 3)') ← domain(right_tail)
+    nc_space = (codomain(right_tail)[1] ⊕ _lastspace(Q)') ← domain(right_tail)
     nc, Qr2 = lq_compact!(absorb!(zerovector!(similar(right_tail, nc_space)), right_tail))
     AL_exp = catdomain(left, Q)
 
diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl
index 2d72fae4d..bbfa57274 100644
--- a/test/algorithms/changebonds.jl
+++ b/test/algorithms/changebonds.jl
@@ -12,6 +12,7 @@ using TensorKit: ℙ
 using Random
 
 spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))]
+maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 2:length(ψ))
 
 
 @testset "MPO $(spacetype(pspace))" for (pspace, Dspace) in spacelist
@@ -78,21 +79,23 @@ end
     state_re = changebonds(
         state, RandExpand(; trunc = truncrank(dim(Dspace) * dim(Dspace)))
     )
+    @test maxbond(state_re) > maxbond(state)
     @test dot(state, state_re) ≈ 1 atol = 1.0e-8
 
     state_oe, _ = changebonds(
         state, H, OptimalExpand(; trunc = truncrank(dim(Dspace) * dim(Dspace)))
     )
+    @test maxbond(state_oe) > maxbond(state)
     @test dot(state, state_oe) ≈ 1 atol = 1.0e-8
 
     state_se, _ = changebonds(
         state, H,
         SketchedExpand(; trunc = truncrank(dim(Dspace) * dim(Dspace)), oversampling = 4)
     )
+    @test maxbond(state_se) > maxbond(state)
     @test dot(state, state_se) ≈ 1 atol = 1.0e-8
 
     state_tr = changebonds(state_oe, SvdCut(; trunc = truncrank(dim(Dspace))))
-
     @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5))
 end
 
@@ -124,15 +127,14 @@ 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
-# because `FiniteMPOHamiltonian`/`FiniteMPO` only support a single physical leg per site.
+# physical leg. Operator-based expanders (`OptimalExpand`/`SketchedExpand`) make use of a
+# one-sided MPO application on the first physical leg.
 @testset "Density-matrix FiniteMPS $(spacetype(pcomp))" for (pcomp, Dspace) in [
         (ℙ^2 ⊗ (ℙ^2)', ℙ^6),
         (Rep[SU₂](1 // 2 => 1) ⊗ Rep[SU₂](1 // 2 => 1)', Rep[SU₂](0 => 4, 1 => 3)),
     ]
     Random.seed!(2468)
     L = 8
-    maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 2:length(ψ))
 
     ψ = FiniteMPS(rand, ComplexF64, fill(pcomp, L), Dspace)
     @test numind(ψ.AC[L ÷ 2]) == 4    # two physical legs + two virtual legs
@@ -146,7 +148,28 @@ end
     # SvdCut truncates the enlarged bond back down, leaving a normalized state
     ψ_tr = changebonds(ψ_re, SvdCut(; trunc = truncrank(dim(Dspace))))
     @test maxbond(ψ_tr) < maxbond(ψ_re)
-    @test abs(dot(ψ_tr, ψ_tr)) ≈ 1 atol = 1.0e-8
+    @test abs(dot(ψ, ψ_tr)) ≈ 1 atol = 1.0e-8
+
+    # use random time-evolution MPO to test operator-based expanders
+    pspace = pcomp[1]
+    nn = rand(ComplexF64, pspace * pspace, pspace * pspace)
+    nn += nn'
+    H = FiniteMPOHamiltonian(fill(pspace, L), (i, i + 1) => nn for i in 1:(L - 1))
+    beta = 0.1
+    O = make_time_mpo(H, beta, TaylorCluster(; N = 2); imaginary_evolution = true)
+
+    ψ_oe, _ = changebonds(
+        ψ, O, OptimalExpand(; trunc = truncrank(dim(Dspace) * 2))
+    )
+    @test maxbond(ψ_oe) > maxbond(ψ)
+    @test dot(ψ, ψ_oe) ≈ 1 atol = 1.0e-8
+
+    ψ_se, _ = changebonds(
+        ψ, O,
+        SketchedExpand(; trunc = truncrank(dim(Dspace) * 2), oversampling = 4)
+    )
+    @test maxbond(ψ_se) > maxbond(ψ)
+    @test dot(ψ, ψ_se) ≈ 1 atol = 1.0e-8
 end
 
 @testset "MultilineMPS $(spacetype(pspace))" for (pspace, Dspace) in spacelist