diff --git a/src/manybody.jl b/src/manybody.jl
index 4386ad24..dfb7f25a 100644
--- a/src/manybody.jl
+++ b/src/manybody.jl
@@ -167,6 +167,9 @@ number(b::ManyBodyBasis) = number(ComplexF64, b)
     transition([T=ComplexF64,] b::ManyBodyBasis, to, from)
 
 Operator ``|\\mathrm{to}⟩⟨\\mathrm{from}|`` transferring particles between modes.
+
+Note that `to` and `from` can be collections of indices. The resulting operator in this case
+will be equal to ``a^\\dagger_{to_1} a^\\dagger_{to_2} \\ldots a_{from_2} a_{from_1}``.
 """
 function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T}
     Is = Int[]
@@ -175,16 +178,13 @@ function transition(::Type{T}, b::ManyBodyBasis, to, from) where {T}
     buffer = allocate_buffer(b.occupations)
     # <{m}_j| at_to a_from |{m}_i>
     for (i, occ_i) in enumerate(b.occupations)
-        if occ_i[from] == 0
-            continue
-        end
         C = state_transition!(buffer, occ_i, from, to)
         C === nothing && continue
         j = state_index(b.occupations, buffer)
         j === nothing && continue
         push!(Is, j)
         push!(Js, i)
-        push!(Vs, sqrt(occ_i[from]) * sqrt(buffer[to]))
+        push!(Vs, C)
     end
     return SparseOperator(b, sparse(Is, Js, Vs, length(b), length(b)))
 end
@@ -373,7 +373,7 @@ end
 
 Base.@propagate_inbounds function state_transition!(buffer, occ_m, at_indices, a_indices)
     any(==(0), (occ_m[m] for m in at_indices)) && return nothing
-    result = 1.0
+    result = 1
     copyto!(buffer, occ_m)
     for i in at_indices
         result *= buffer[i]
diff --git a/test/test_manybody.jl b/test/test_manybody.jl
index d83fd617..0fb6d056 100644
--- a/test/test_manybody.jl
+++ b/test/test_manybody.jl
@@ -9,6 +9,19 @@ Random.seed!(0)
 D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
 D(x1::StateVector, x2::StateVector) = norm(x2-x1)
 
+# Test `SortedVector`
+ve = [3, 2, 1, 4, 5]
+sv1 = QuantumOpticsBase.SortedVector(ve)
+sv2 = QuantumOpticsBase.SortedVector(ve, Base.Reverse)
+@test length(sv1) == 5
+@test length(sv2) == 5
+@test collect(sv1) == [1, 2, 3, 4, 5]
+@test collect(sv2) == [5, 4, 3, 2, 1]
+@test map(x -> x^2, sv1) == [1, 4, 9, 16, 25]
+@test map(x -> x^2, sv2) == [25, 16, 9, 4, 1]
+@test findfirst(iseven, sv1) == 2
+@test findfirst(iseven, sv2) == 2
+
 # Test state creation
 Nmodes = 5
 b = GenericBasis(Nmodes)
@@ -105,6 +118,10 @@ t32 = transition(b_mb, 3, 2)
 @test 1e-12 > D(t32*basisstate(b_mb, [0, 2, 1, 0]), 2*basisstate(b_mb, [0, 1, 2, 0]))
 @test 1e-12 > D(number(b_mb, 2), transition(b_mb, 2, 2))
 
+t22_33 = transition(b_mb, (2, 2), (3, 3))
+d3 = destroy(b_mb, 3)
+c2 = create(b_mb, 2)
+@test t22_33 ≈ c2 * c2 * d3 * d3
 
 # Single particle operator in second quantization
 b_single = GenericBasis(Nmodes)