diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 1788564bd..61b69d258 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -96,6 +96,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   and the right virtual leg of `mpo[end]` and contracted the two — which at length 1 is the *same*
   tensor, so it returned `O * O` on twice the physical space instead of `O`.
   ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484))
+- `make_time_mpo` with `TaylorCluster` on a Hamiltonian whose virtual bond dimension varies along
+  the chain. The loopback step labels rows by the left virtual bond and columns by the right one, but deleted the absorbed "three-level" rows and columns using the *right*-bond linear index for both. When the two bonds differ in size this deletes the wrong rows, leaving spurious rows behind and removing valid ones, so the resulting time-evolution MPO was wrong. Rows and columns are now derived from their own bond, and the removal happens in a single pass after all loopback contributions have been accumulated instead of inside the loop that still reads those entries.
+  ([#511](https://github.com/QuantumKitHub/MPSKit.jl/pull/511))
 
 ### Performance
 
diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl
index b129e9d0c..89ad94026 100644
--- a/src/algorithms/timestep/taylorcluster.jl
+++ b/src/algorithms/timestep/taylorcluster.jl
@@ -130,21 +130,57 @@ function _taylor_loopback!(mpo, virtual_sz, linds, ::Val{N}, τ::Number) where {
         V_right = virtual_sz[i + 1]
         linds_right = linds[i + 1]
         cinds_right = CartesianIndices(linds_right)
-        for b in cinds_right[2:end]
-            all(in((1, V_right)), b.I) || continue
 
-            b_lin = linds_right[b]
-            a = count(==(V_right), b.I)
+        threelevel_linds_right = _get_threelevel_linds(V_right, linds_right)
+        for c in threelevel_linds_right
+            c_cart = cinds_right[c]
+            a = count(==(V_right), c_cart.I)
             factor = τ^a * factorial(N - a) / factorial(N)
-            slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin]
-            for I in nonzero_keys(slice)
-                (I[1] == b_lin || I[4] == b_lin) && delete!(slice, I)
-            end
+            slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, c]
         end
     end
+    _remove_threelevels!(mpo, virtual_sz, linds)
     return mpo
 end
 
+function _get_threelevel_linds(virtual_sz, linds)
+    row_threelevels_linds = Int[]
+    for c in CartesianIndices(linds)
+        c_lin = linds[c]
+        cI = c.I
+        all(in((1, virtual_sz)), cI) && !all(==(1), cI) && push!(row_threelevels_linds, c_lin)
+    end
+    return row_threelevels_linds
+end
+
+# removes the rows and columns that are labled by level-labels made up of 1 and "3" only, but not all 1's.
+function _remove_threelevels!(mpo, virtual_sz, linds)
+    for (i, slice) in enumerate(parent(mpo))
+        # three-level rows
+        row_threelevel_linds = _get_threelevel_linds(virtual_sz[i], linds[i])
+
+        # three-level columns
+        col_threelevel_linds = _get_threelevel_linds(virtual_sz[i + 1], linds[i + 1])
+
+        # build the set of all CartesianIndices in a three-level row or column
+        d1, _, _, d4 = size(slice)
+        threelevel_cinds = Set{CartesianIndex{4}}()
+        for r in row_threelevel_linds
+            union!(threelevel_cinds, CartesianIndices((r:r, 1:1, 1:1, 1:d4)))
+        end
+        for c in col_threelevel_linds
+            union!(threelevel_cinds, CartesianIndices((1:d1, 1:1, 1:1, c:c)))
+        end
+
+        # remove the nonzero entries that fall in a three-level row or column
+        # this avoids a lot of misses wrt looping over all nonzero_keys
+        for I in intersect(nonzero_keys(slice), threelevel_cinds)
+            delete!(slice, I)
+        end
+    end
+    return
+end
+
 # Algorithm 2: collapse rows and columns that are equivalent under the
 # permutation symmetry of the Taylor expansion.
 function _taylor_remove_equivalents!(mpo, virtual_sz, linds)