diff --git a/src/algorithms/excitation/exci_transfer_system.jl b/src/algorithms/excitation/exci_transfer_system.jl
index bb0b90095..c2b4b217d 100644
--- a/src/algorithms/excitation/exci_transfer_system.jl
+++ b/src/algorithms/excitation/exci_transfer_system.jl
@@ -1,47 +1,3 @@
-function left_excitation_transfer_system(
-        lBs, H, exci; mom = exci.momentum,
-        solver = Defaults.linearsolver
-    )
-    len = length(H)
-    found = zero.(lBs)
-    odim = length(lBs)
-
-    for i in 1:odim
-        # this operation can in principle be even further optimized for larger unit cells
-        # as we only require the terms that end at level i.
-        # this would require to check the finite state machine, and discard non-connected
-        # terms.
-        H_partial = map(site -> H.data[site, 1:i, 1:i], 1:len)
-        T = TransferMatrix(exci.right_gs.AR, H_partial, exci.left_gs.AL)
-        start = scale!(last(found[1:i] * T), cis(-mom * len))
-        if exci.trivial && isid(H, i)
-            @plansor start[-1 -2; -3 -4] -= start[1 4; -3 2] * r_RL(exci.right_gs)[2; 3] *
-                τ[3 4; 5 1] * l_RL(exci.right_gs)[-1; 6] * τ[5 6; -4 -2]
-        end
-
-        found[i] = add!(start, lBs[i])
-
-        if reduce(&, contains.(H.data, i, i))
-            if isid(H, i)
-                tm = TransferMatrix(exci.right_gs.AR, exci.left_gs.AL)
-                if exci.trivial
-                    tm = regularize(tm, l_RL(exci.right_gs), r_RL(exci.right_gs))
-                end
-            else
-                tm = TransferMatrix(
-                    exci.right_gs.AR, getindex.(H.data, i, i), exci.left_gs.AL
-                )
-            end
-
-            found[i], convhist = linsolve(
-                flip(tm), found[i], found[i], solver, 1, -cis(-mom * len)
-            )
-            convhist.converged == 0 &&
-                @warn "GBL$i failed to converge: normres = $(convhist.normres)"
-        end
-    end
-    return found
-end
 function left_excitation_transfer_system(
         GBL, H::InfiniteMPOHamiltonian, exci;
         mom = exci.momentum, solver = Defaults.linearsolver
@@ -50,6 +6,11 @@ function left_excitation_transfer_system(
     found = zerovector(GBL)
     odim = length(GBL)
 
+    if istrivial(exci)
+        ρ_left = l_RL(exci.right_gs)
+        ρ_right = r_RL(exci.right_gs)
+    end
+
     for i in 1:odim
         # this operation can in principle be even further optimized for larger unit cells
         # as we only require the terms that end at level i.
@@ -58,9 +19,8 @@ function left_excitation_transfer_system(
         H_partial = map(h -> getindex(h, 1:i, 1, 1, 1:i), parent(H))
         T = TransferMatrix(exci.right_gs.AR, H_partial, exci.left_gs.AL)
         start = scale!(last(found[1:i] * T), cis(-mom * len))
-        if exci.trivial && isidentitylevel(H, i)
-            @plansor start[-1 -2; -3 -4] -= start[1 4; -3 2] * r_RL(exci.right_gs)[2; 3] *
-                τ[3 4; 5 1] * l_RL(exci.right_gs)[-1; 6] * τ[5 6; -4 -2]
+        if istrivial(exci) && isidentitylevel(H, i)
+            regularize!(start, ρ_right, ρ_left)
         end
 
         found[i] = add!(start, GBL[i])
@@ -68,8 +28,8 @@ function left_excitation_transfer_system(
         if !isemptylevel(H, i)
             if isidentitylevel(H, i)
                 T = TransferMatrix(exci.right_gs.AR, exci.left_gs.AL)
-                if exci.trivial
-                    T = regularize(T, l_RL(exci.right_gs), r_RL(exci.right_gs))
+                if istrivial(exci)
+                    T = regularize(T, ρ_left, ρ_right)
                 end
             else
                 T = TransferMatrix(
@@ -86,49 +46,7 @@ function left_excitation_transfer_system(
     end
     return found
 end
-function right_excitation_transfer_system(
-        rBs, H, exci; mom = exci.momentum, solver = Defaults.linearsolver
-    )
-    len = length(H)
-    found = zero.(rBs)
-    odim = length(rBs)
-
-    for i in odim:-1:1
-        # this operation can in principle be even further optimized for larger unit cells
-        # as we only require the terms that end at level i.
-        # this would require to check the finite state machine, and discard non-connected
-        # terms.
-        H_partial = map(site -> H.data[site, i:odim, i:odim], 1:len)
-        T = TransferMatrix(exci.left_gs.AL, H_partial, exci.right_gs.AR)
-        start = scale!(first(T * found[i:odim]), cis(mom * len))
-        if exci.trivial && isid(H, i)
-            @plansor start[-1 -2; -3 -4] -= τ[6 2; 3 4] * start[3 4; -3 5] *
-                l_LR(exci.right_gs)[5; 2] * r_LR(exci.right_gs)[-1; 1] * τ[-2 -4; 1 6]
-        end
-
-        found[i] = add!(start, rBs[i])
 
-        if reduce(&, contains.(H.data, i, i))
-            if isid(H, i)
-                tm = TransferMatrix(exci.left_gs.AL, exci.right_gs.AR)
-                if exci.trivial
-                    tm = regularize(tm, l_LR(exci.left_gs), r_LR(exci.right_gs))
-                end
-            else
-                tm = TransferMatrix(
-                    exci.left_gs.AL, getindex.(H.data, i, i), exci.right_gs.AR
-                )
-            end
-
-            found[i], convhist = linsolve(
-                tm, found[i], found[i], solver, 1, -cis(mom * len)
-            )
-            convhist.converged < 1 &&
-                @warn "GBR$i failed to converge: normres = $(convhist.normres)"
-        end
-    end
-    return found
-end
 function right_excitation_transfer_system(
         GBR, H::InfiniteMPOHamiltonian, exci;
         mom = exci.momentum,
@@ -138,6 +56,11 @@ function right_excitation_transfer_system(
     found = zerovector(GBR)
     odim = length(GBR)
 
+    if istrivial(exci)
+        ρ_left = l_LR(exci.right_gs)
+        ρ_right = r_LR(exci.right_gs)
+    end
+
     for i in odim:-1:1
         # this operation can in principle be even further optimized for larger unit cells
         # as we only require the terms that end at level i.
@@ -146,9 +69,8 @@ function right_excitation_transfer_system(
         H_partial = map(h -> h[i:end, 1, 1, i:end], parent(H))
         T = TransferMatrix(exci.left_gs.AL, H_partial, exci.right_gs.AR)
         start = scale!(first(T * found[i:odim]), cis(mom * len))
-        if exci.trivial && isidentitylevel(H, i)
-            @plansor start[-1 -2; -3 -4] -= τ[6 2; 3 4] * start[3 4; -3 5] *
-                l_LR(exci.right_gs)[5; 2] * r_LR(exci.right_gs)[-1; 1] * τ[-2 -4; 1 6]
+        if istrivial(exci) && isidentitylevel(H, i)
+            regularize!(start, ρ_left, ρ_right)
         end
 
         found[i] = add!(start, GBR[i])
@@ -156,8 +78,8 @@ function right_excitation_transfer_system(
         if !isemptylevel(H, i)
             if isidentitylevel(H, i)
                 tm = TransferMatrix(exci.left_gs.AL, exci.right_gs.AR)
-                if exci.trivial
-                    tm = regularize(tm, l_LR(exci.left_gs), r_LR(exci.right_gs))
+                if istrivial(exci)
+                    tm = regularize(tm, ρ_left, ρ_right)
                 end
             else
                 tm = TransferMatrix(
@@ -166,8 +88,7 @@ function right_excitation_transfer_system(
             end
 
             found[i], convhist = linsolve(
-                tm, found[i], found[i], solver, 1,
-                -cis(mom * len)
+                tm, found[i], found[i], solver, 1, -cis(mom * len)
             )
             convhist.converged < 1 &&
                 @warn "GBR$i failed to converge: normres = $(convhist.normres)"
diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl
index a2fc2759c..0eb400258 100644
--- a/src/algorithms/excitation/quasiparticleexcitation.jl
+++ b/src/algorithms/excitation/quasiparticleexcitation.jl
@@ -57,7 +57,7 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renv
 end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs; num = 1, kwargs...)
     # Infer `renvs` in function body as it depends on `solver`.
-    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
+    renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
     return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...)
 end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP; num = 1, kwargs...)
@@ -136,7 +136,7 @@ end
 function excitations(
         H, alg::QuasiparticleAnsatz, ϕ₀::FiniteQP,
         lenvs = environments(ϕ₀.left_gs, H),
-        renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H);
+        renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H);
         num = 1
     )
     E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs)
@@ -223,7 +223,7 @@ function excitations(
         kwargs...
     )
     # Infer `renvs` in function body as it depends on `solver`.
-    renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
+    renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...)
     return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...)
 end
 function excitations(
@@ -372,7 +372,7 @@ function effective_excitation_renormalization_energy(H, ϕ, lenvs, renvs)
         E[i] = contract_mpo_expval(
             ψ_left.AC[i], leftenv(lenvs, i, ψ_left), H[i], rightenv(lenvs, i, ψ_left)
         )
-        if !ϕ.trivial
+        if istopological(ϕ)
             E[i] += contract_mpo_expval(
                 ψ_right.AC[i], leftenv(renvs, i, ψ_right), H[i], rightenv(renvs, i, ψ_right)
             )
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 9bf7c180e..150f276fb 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -225,7 +225,7 @@ end
 
 function variance(state::InfiniteQP, H::InfiniteMPOHamiltonian, envs = environments(state, H))
     # I remember there being an issue here @gertian?
-    state.trivial ||
+    istopological(state) &&
         throw(ArgumentError("variance of domain wall excitations is not implemented"))
     gs = state.left_gs
 
diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl
index acc9fa9f3..4aacf0dae 100644
--- a/src/environments/qp_envs.jl
+++ b/src/environments/qp_envs.jl
@@ -33,7 +33,7 @@ function environments(exci::Union{InfiniteQP, MultilineQP}, H; kwargs...)
     return environments(exci, H, lenvs; kwargs...)
 end
 function environments(exci::Union{InfiniteQP, MultilineQP}, H, lenvs; kwargs...)
-    renvs = exci.trivial ? lenvs : environments(exci.right_gs, H; kwargs...)
+    renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H; kwargs...)
     return environments(exci, H, lenvs, renvs; kwargs...)
 end
 
@@ -46,7 +46,7 @@ function environments(qp::MultilineQP, operator::MultilineMPO, lenvs, renvs; kwa
 end
 
 function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; kwargs...)
-    ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1))
+    ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1)) .+ 1
     solver = environment_alg(exci, H, exci; kwargs...)
 
     AL = exci.left_gs.AL
@@ -62,11 +62,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
         lBs[pos + 1] += leftenv(lenvs, pos, exci.left_gs) *
             TransferMatrix(exci[pos], H[pos], AL[pos]) / cis(exci.momentum)
 
-        if exci.trivial # regularization of trivial excitations
+        if istrivial(exci) && !isempty(ids) # regularization of trivial excitations
+            ρ_left = l_RL(exci.left_gs, pos + 1)
+            ρ_right = r_RL(exci.left_gs, pos)
             for i in ids
-                @plansor lBs[pos + 1][i][-1 -2; -3 -4] -= lBs[pos + 1][i][1 4; -3 2] *
-                    r_RL(exci.left_gs, pos)[2; 3] * τ[3 4; 5 1] *
-                    l_RL(exci.left_gs, pos + 1)[ -1; 6 ] * τ[5 6; -4 -2]
+                regularize!(lBs[pos + 1][i], ρ_right, ρ_left)
             end
         end
     end
@@ -78,12 +78,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
         rBs[pos - 1] += TransferMatrix(exci[pos], H[pos], AR[pos]) *
             rightenv(renvs, pos, exci.right_gs) * cis(exci.momentum)
 
-        if exci.trivial
+        if istrivial(exci) && !isempty(ids)
+            ρ_left = l_LR(exci.left_gs, pos)
+            ρ_right = r_LR(exci.left_gs, pos - 1)
             for i in ids
-                ρ_left = l_LR(exci.left_gs, pos)
-                ρ_right = r_LR(exci.left_gs, pos - 1)
-                @plansor rBs[pos - 1][i][-1 -2; -3 -4] -= τ[6 4; 1 3] *
-                    rBs[pos - 1][i][1 3; -3 2] * ρ_left[2; 4] * ρ_right[-1; 5] * τ[-2 -4; 5 6]
+                regularize!(rBs[pos - 1][i], ρ_left, ρ_right)
             end
         end
     end
@@ -101,12 +100,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
     for i in 1:(length(exci) - 1)
         lB_cur = lB_cur * TransferMatrix(AR[i], H[i], AL[i]) / cis(exci.momentum)
 
-        if exci.trivial
+        if istrivial(exci) && !isempty(ids)
+            ρ_left = l_RL(exci.left_gs, i + 1)
+            ρ_right = r_RL(exci.left_gs, i)
             for k in ids
-                ρ_left = l_RL(exci.left_gs, i + 1)
-                ρ_right = r_RL(exci.left_gs, i)
-                @plansor lB_cur[k][-1 -2; -3 -4] -= lB_cur[k][1 4; -3 2] *
-                    ρ_right[2; 3] * τ[3 4; 5 1] * ρ_left[-1; 6] * τ[5 6; -4 -2]
+                regularize!(lB_cur[k], ρ_right, ρ_left)
             end
         end
 
@@ -117,12 +115,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs;
     for i in length(exci):-1:2
         rB_cur = TransferMatrix(AL[i], H[i], AR[i]) * rB_cur * cis(exci.momentum)
 
-        if exci.trivial
+        if istrivial(exci) && !isempty(ids)
+            ρ_left = l_LR(exci.left_gs, i)
+            ρ_right = r_LR(exci.left_gs, i - 1)
             for k in ids
-                ρ_left = l_LR(exci.left_gs, i)
-                ρ_right = r_LR(exci.left_gs, i - 1)
-                @plansor rB_cur[k][-1 -2; -3 -4] -= τ[6 4; 1 3] *
-                    rB_cur[k][1 3; -3 2] * ρ_left[2; 4] * ρ_right[-1; 5] * τ[-2 -4; 5 6]
+                regularize!(rB_cur[k], ρ_left, ρ_right)
             end
         end
 
@@ -135,7 +132,7 @@ end
 function environments(
         exci::FiniteQP, H::FiniteMPOHamiltonian,
         lenvs = environments(exci.left_gs, H),
-        renvs = exci.trivial ? lenvs : environments(exci.right_gs, H);
+        renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H);
         kwargs...
     )
     AL = exci.left_gs.AL
@@ -164,7 +161,7 @@ function environments(
 end
 
 function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...)
-    exci.trivial ||
+    istopological(exci) &&
         @warn "there is a phase ambiguity in topologically nontrivial statmech excitations"
     solver = environment_alg(exci, O, exci; kwargs...)
 
@@ -206,7 +203,7 @@ function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...)
     T_RL = TransferMatrix(right_gs.AR, O, left_gs.AL)
     T_LR = TransferMatrix(left_gs.AL, O, right_gs.AR)
 
-    if exci.trivial
+    if istrivial(exci)
         @plansor rvec[-1 -2; -3] := rightenv(lenvs, 0, left_gs)[-1 -2; 1] *
             conj(left_gs.C[0][-3; 1])
         @plansor lvec[-1 -2; -3] := leftenv(lenvs, 1, left_gs)[-1 -2; 1] *
diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl
index bd1f827b2..cf476f914 100644
--- a/src/states/quasiparticle_state.jl
+++ b/src/states/quasiparticle_state.jl
@@ -133,7 +133,7 @@ function Base.convert(
 
     tm = TransferMatrix(input.left_gs.AL, input.right_gs.AR)
 
-    if input.trivial
+    if istrivial(input)
         tm = regularize(tm, l_LR(input.right_gs), r_LR(input.right_gs))
     end
 
@@ -180,7 +180,7 @@ function Base.convert(
     end
 
     tm = TransferMatrix(input.right_gs.AR, input.left_gs.AL)
-    if input.trivial
+    if istrivial(input)
         tm = regularize(tm, l_RL(input.right_gs), r_RL(input.right_gs))
     end
 
@@ -222,9 +222,12 @@ left_virtualspace(state::QP) = map(Base.Fix1(left_virtualspace, state), eachsite
 right_virtualspace(state::QP, i::Int) = right_virtualspace(state.right_gs, i)
 right_virtualspace(state::QP) = map(Base.Fix1(right_virtualspace, state), eachsite(state))
 auxiliaryspace(state::QP) = space(state.Xs[1], 2)
-
+auxiliarysector(state::QP) = only(sectors(auxiliaryspace(state)))
 eachsite(state::QP) = eachsite(state.left_gs)
 
+istopological(qp::QP) = qp.left_gs !== qp.right_gs
+istrivial(qp::QP) = !istopological(qp) && isone(auxiliarysector(qp))
+
 Base.copy(a::QP) = copy!(similar(a), a)
 Base.copyto!(a::QP, b::QP) = copy!(a, b)
 function Base.copy!(a::T, b::T) where {T <: QP}
@@ -233,11 +236,12 @@ function Base.copy!(a::T, b::T) where {T <: QP}
     end
     return a
 end
-function Base.getproperty(v::QP, s::Symbol)
+Base.@constprop :aggressive function Base.getproperty(qp::QP, s::Symbol)
     if s == :trivial
-        return v.left_gs === v.right_gs
+        Base.depwarn("`qp.trivial` is deprecated in favor of `istrivial` and `istopological`", :trivial)
+        return !istopological(qp)
     else
-        return getfield(v, s)
+        return getfield(qp, s)
     end
 end
 
@@ -402,7 +406,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S <: FiniteMPS}
     return FiniteMPS(Ls + Rs + Bs; normalize = false)
 end
 
-function Base.getproperty(exci::MultilineQP, s::Symbol)
+Base.@constprop :aggressive function Base.getproperty(exci::MultilineQP, s::Symbol)
     if s == :momentum
         return first(exci.data).momentum
     elseif s == :left_gs
@@ -416,6 +420,10 @@ function Base.getproperty(exci::MultilineQP, s::Symbol)
     end
 end
 
+# These should really all be the same, so it might make sense to take the first instead
+istrivial(exci::MultilineQP) = all(istrivial, exci.data)
+istopological(exci::MultilineQP) = all(istopological, exci.data)
+
 # VectorInterface
 # ---------------
 
diff --git a/src/transfermatrix/transfermatrix.jl b/src/transfermatrix/transfermatrix.jl
index 88030e91d..31a33ab06 100644
--- a/src/transfermatrix/transfermatrix.jl
+++ b/src/transfermatrix/transfermatrix.jl
@@ -87,6 +87,6 @@ function regularize!(v::MPOTensor, lvec::MPSTensor, rvec::MPSTensor)
 end
 
 function regularize!(v::MPOTensor, lvec::MPSBondTensor, rvec::MPSBondTensor)
-    return @plansor v[-1 -2; -3 -4] -= τ[6 2; 3 4] * v[3 4; -3 5] * lvec[5; 2] * rvec[-1; 1] *
-        τ[-2 -4; 1 6]
+    λ = @plansor lvec[2; 1] * removeunit(removeunit(v, 3), 2)[1; 2]
+    return add!(v, insertleftunit(insertrightunit(rvec, 1; dual = isdual(space(v, 2))), 3), -λ)
 end
diff --git a/test/algorithms.jl b/test/algorithms.jl
index 469506e39..b009c121b 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -475,6 +475,21 @@ module TestAlgorithms
             @test energies[1] ≈ 0.41047925 atol = 1.0e-4
             @test variance(ϕs[1], H) < 1.0e-8
         end
+        @testset "infinite sector convention" begin
+            g = 4
+            H = repeat(transverse_field_ising(Z2Irrep; g = g, L = Inf), 2)
+            V = Z2Space(0 => 24, 1 => 24)
+            ψ = InfiniteMPS(physicalspace(H), [V, V])
+            ψ, envs = find_groundstate(ψ, H, VUMPS(; tol = 1.0e-10, maxiter = 400))
+
+            # testing istrivial and istopological
+            momentum = 0
+            exc0, qp0 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(0))
+            exc1, qp1 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(1))
+            @test isapprox(first(exc1), abs(2 * (g - 1)); atol = 1.0e-6) # charged excitation lower in energy
+            @test first(exc0) > 2 * first(exc1)
+            @test variance(qp1[1], H) < 1.0e-8
+        end
         @testset "infinite (mpo)" begin
             H = repeat(sixvertex(), 2)
             ψ = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^10, ℂ^10])
diff --git a/test/setup.jl b/test/setup.jl
index 09a18a83f..a6b07bbc0 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -94,6 +94,14 @@ function S_x(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) w
         throw(ArgumentError("spin $spin not supported"))
     end
 end
+function S_x(::Type{Z2Irrep}, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number}
+    spin == 1 // 2 || throw(ArgumentError("spin $spin not supported"))
+    pspace = Z2Space(0 => 1, 1 => 1)
+    X = zeros(T, pspace, pspace)
+    block(X, Z2Irrep(0)) .= one(T) # no times 2
+    block(X, Z2Irrep(1)) .= -one(T)
+    return X
+end
 function S_y(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number}
     return if spin == 1 // 2
         TensorMap(T[0 -im; im 0], ℂ^2 ← ℂ^2)
@@ -121,17 +129,29 @@ end
 function S_zz(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number}
     return S_z(Trivial, T; spin) ⊗ S_z(Trivial, T; spin)
 end
+function S_zz(::Type{Z2Irrep}, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number}
+    spin == 1 // 2 || throw(ArgumentError("spin $spin not supported"))
+    P = Z2Space(0 => 1, 1 => 1)
+    ZZ = zeros(ComplexF64, P ⊗ P ← P ⊗ P)
+    flip_charge(charge::Z2Irrep) = only(charge ⊗ Z2Irrep(1))
+    for (s, f) in fusiontrees(ZZ)
+        if s.uncoupled == map(flip_charge, f.uncoupled)
+            ZZ[s, f] .= 1 # no divide by 4
+        end
+    end
+    return ZZ
+end
 
-function transverse_field_ising(::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number}
-    X = S_x(Trivial, T; spin = 1 // 2)
-    ZZ = S_zz(Trivial, T; spin = 1 // 2)
+function transverse_field_ising(sym::Type{<:Sector} = Trivial, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number}
+    X = S_x(sym, T; spin = 1 // 2)
+    ZZ = S_zz(sym, T; spin = 1 // 2)
 
     if L == Inf
-        lattice = PeriodicArray([ℂ^2])
+        lattice = PeriodicArray([space(X, 1)])
         H₁ = InfiniteMPOHamiltonian(lattice, i => -g * X for i in 1:1)
         H₂ = InfiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:1)
     else
-        lattice = fill(ℂ^2, L)
+        lattice = fill(space(X, 1), L)
         H₁ = FiniteMPOHamiltonian(lattice, i => -g * X for i in 1:L)
         H₂ = FiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:(L - 1))
     end