diff --git a/src/operators_lazytensor.jl b/src/operators_lazytensor.jl
index 16716cc1..c729d7a5 100644
--- a/src/operators_lazytensor.jl
+++ b/src/operators_lazytensor.jl
@@ -90,7 +90,7 @@ SparseArrays.sparse(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op
 -(a::LazyTensor) = LazyTensor(a, -a.factor)
 
 function *(a::LazyTensor{B1,B2}, b::LazyTensor{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}
-    indices = sortedindices.union(a.indices, b.indices)
+    indices = sort(union(a.indices, b.indices))
     ops = Vector{AbstractOperator}(undef, length(indices))
     for n in 1:length(indices)
         i = indices[n]
diff --git a/src/sortedindices.jl b/src/sortedindices.jl
index 9e9ff2d9..0e8164b2 100644
--- a/src/sortedindices.jl
+++ b/src/sortedindices.jl
@@ -19,43 +19,6 @@ function complement(N::Int, indices::Vector{Int})
     x
 end
 
-"""
-[1, 4, 5], [2, 4, 7] => [1, 2, 4, 5, 7]
-"""
-function union(ind1::Vector{Int}, ind2::Vector{Int})
-    i1 = 1
-    i2 = 1
-    N1 = length(ind1)
-    N2 = length(ind2)
-    xvec = Vector{Int}()
-    while true
-        if i1 > N1
-            for j=i2:N2
-                push!(xvec, ind2[j])
-            end
-            return xvec
-        elseif i2 > N2
-            for j=i1:N1
-                push!(xvec, ind1[j])
-            end
-            return xvec
-        end
-        x1 = ind1[i1]
-        x2 = ind2[i2]
-        if x1 == x2
-            i1 += 1
-            i2 += 1
-            push!(xvec, x1)
-        elseif x1 < x2
-            i1 += 1
-            push!(xvec, x1)
-        else
-            i2 += 1
-            push!(xvec, x2)
-        end
-    end
-end
-
 
 """
 [1, 4, 5], [2, 4, 7] => [1, 5]
@@ -91,46 +54,6 @@ function shiftremove(ind1::Vector{Int}, ind2::Vector{Int})
     x
 end
 
-"""
-[2, 3, 6], [1, 3, 4, 6, 7] => [3, 6]
-"""
-function intersect(ind1::Vector{Int}, ind2::Vector{Int})
-    i1 = 1
-    i2 = 1
-    N1 = length(ind1)
-    N2 = length(ind2)
-    xvec = Vector{Int}()
-    if i1 > N1 || i2 > N2
-        return xvec
-    end
-    x1 = ind1[i1]
-    x2 = ind2[i2]
-    while true
-        if x1 == x2
-            i1 += 1
-            i2 += 1
-            push!(xvec, x1)
-            if i1 > N1 || i2 > N2
-                return xvec
-            end
-            x1 = ind1[i1]
-            x2 = ind2[i2]
-        elseif x1 < x2
-            i1 += 1
-            if i1 > N1
-                return xvec
-            end
-            x1 = ind1[i1]
-        else
-            i2 += 1
-            if i2 > N2
-                return xvec
-            end
-            x2 = ind2[i2]
-        end
-    end
-end
-
 function reducedindices(I_::Vector{Int}, I::Vector{Int})
     N = length(I_)
     x = Vector{Int}(undef, N)
diff --git a/test/test_sortedindices.jl b/test/test_sortedindices.jl
index ac427a9a..cc86e98b 100644
--- a/test/test_sortedindices.jl
+++ b/test/test_sortedindices.jl
@@ -8,16 +8,6 @@ s = QuantumOptics.sortedindices
 
 @test s.complement(6, [1, 4]) == [2, 3, 5, 6]
 
-@test s.union(Int[], Int[]) == Int[]
-@test s.union(Int[], [2, 3]) == [2, 3]
-@test s.union([2, 3], Int[]) == [2, 3]
-@test s.union([1, 4, 5], [2, 4, 7]) == [1, 2, 4, 5, 7]
-@test s.union([2, 4, 7], [1, 4, 5]) == [1, 2, 4, 5, 7]
-@test s.union([1, 4, 5], [2, 4]) == [1, 2, 4, 5]
-@test s.union([2, 4], [1, 4, 5]) == [1, 2, 4, 5]
-@test s.union([1, 4, 5], [2, 4, 5]) == [1, 2, 4, 5]
-@test s.union([2, 4, 5], [1, 4, 5]) == [1, 2, 4, 5]
-
 @test s.remove([1, 4, 5], [2, 4, 7]) == [1, 5]
 @test s.remove([1, 4, 5, 7], [2, 4, 7]) == [1, 5]
 @test s.remove([1, 4, 5, 8], [2, 4, 7]) == [1, 5, 8]
@@ -26,15 +16,6 @@ s = QuantumOptics.sortedindices
 @test s.shiftremove([1, 4, 5, 7], [2, 4, 7]) == [1, 3]
 @test s.shiftremove([1, 4, 5, 8], [2, 4, 7]) == [1, 3, 5]
 
-@test s.intersect([2, 3, 6], [1, 3, 4, 7]) == [3]
-@test s.intersect([1, 3, 4, 7], [2, 3, 6]) == [3]
-@test s.intersect([2, 3, 6], [1, 3]) == [3]
-@test s.intersect([1, 3], [2, 3, 6]) == [3]
-@test s.intersect([2, 3, 6], [1, 3, 6]) == [3, 6]
-@test s.intersect([1, 3, 6], [2, 3, 6]) == [3, 6]
-@test s.intersect(Int[], [2, 3]) == Int[]
-@test s.intersect([2, 3], Int[]) == Int[]
-
 @test s.reducedindices([3, 5], [2, 3, 5, 6]) == [2, 3]
 x = [3, 5]
 s.reducedindices!(x, [2, 3, 5, 6])