diff --git a/src/operators/multilinempo.jl b/src/operators/multilinempo.jl
index d9c09aa06..d49413e63 100644
--- a/src/operators/multilinempo.jl
+++ b/src/operators/multilinempo.jl
@@ -49,3 +49,6 @@ for f_space in (:physicalspace, :left_virtualspace, :right_virtualspace)
     @eval $f_space(t::MultilineMPO, I::CartesianIndex{2}) = $f_space(t, Tuple(I)...)
     @eval $f_space(t::MultilineMPO) = map(Base.Fix1($f_space, t), eachindex(t))
 end
+
+TensorKit.leftunit(t::MultilineMPO) = TensorKit.leftunit(t[1]) # same for every line
+TensorKit.rightunit(t::MultilineMPO) = TensorKit.rightunit(t[1])
diff --git a/src/states/multilinemps.jl b/src/states/multilinemps.jl
index a090b67cd..eda3c5c35 100644
--- a/src/states/multilinemps.jl
+++ b/src/states/multilinemps.jl
@@ -98,3 +98,6 @@ for f_space in (:physicalspace, :left_virtualspace, :right_virtualspace)
     @eval $f_space(t::MultilineMPS, I::CartesianIndex{2}) = $f_space(t, Tuple(I)...)
     @eval $f_space(t::MultilineMPS) = map(Base.Fix1($f_space, t), eachindex(t))
 end
+
+TensorKit.leftunit(t::MultilineMPS) = TensorKit.leftunit(t[1]) # same for every line
+TensorKit.rightunit(t::MultilineMPS) = TensorKit.rightunit(t[1])
diff --git a/test/operators/mpo.jl b/test/operators/mpo.jl
index afbd99e09..2f2f900f8 100644
--- a/test/operators/mpo.jl
+++ b/test/operators/mpo.jl
@@ -108,6 +108,8 @@ end
     @test physicalspace(multiH, 1, 1) == P
     @test left_virtualspace(multiH, 1, 1) == left_virtualspace(multiH, 2, 1) == V
     @test right_virtualspace(multiH, CartesianIndex(1, 2)) == V
+    @test leftunit(multiH) == leftunit(H) == unit(sectortype(P))
+    @test rightunit(multiH) == rightunit(H) == unit(sectortype(P))
 end
 
 @testset "Adapt" for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1))
diff --git a/test/states/multilinemps.jl b/test/states/multilinemps.jl
index e9fb8a569..cf5456f6f 100644
--- a/test/states/multilinemps.jl
+++ b/test/states/multilinemps.jl
@@ -31,6 +31,7 @@ using TensorKit: ℙ
     @test physicalspace(ψ) == fill(d, 2, 2)
     @test all(x -> x ≾ D, left_virtualspace(ψ))
     @test all(x -> x ≾ D, right_virtualspace(ψ))
+    @test leftunit(ψ) == unit(sectortype(D))
 
     for i in 1:size(ψ, 1), j in 1:size(ψ, 2)
         @plansor difference[-1 -2; -3] := ψ.AL[i, j][-1 -2; 1] * ψ.C[i, j][1; -3] -