diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index a589600a4..098cfb376 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -37,6 +37,8 @@ import SciMLOperators:
     AbstractSciMLOperator,
     MatrixOperator,
     ScalarOperator,
+    ScaledOperator,
+    AddedOperator,
     IdentityOperator,
     cache_operator,
     update_coefficients!,
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 08f3672a6..2408683f2 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -2,10 +2,10 @@
 Functions for generating (common) quantum super-operators.
 =#
 
-export spre, spost, sprepost, lindblad_dissipator
+export spre, spost, sprepost, liouvillian, lindblad_dissipator
 
 # intrinsic functions for super-operators
-# (keep these because they take AbstractMatrix as input)
+# (keep these because they take AbstractMatrix as input and ensure the output is sparse matrix)
 _spre(A::AbstractMatrix, Id::AbstractMatrix) = kron(Id, sparse(A))
 _spre(A::AbstractSparseMatrix, Id::AbstractMatrix) = kron(Id, A)
 _spost(B::AbstractMatrix, Id::AbstractMatrix) = kron(transpose(sparse(B)), Id)
@@ -14,9 +14,22 @@ _sprepost(A::AbstractMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), spa
 _sprepost(A::AbstractMatrix, B::AbstractSparseMatrix) = kron(transpose(B), sparse(A))
 _sprepost(A::AbstractSparseMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), A)
 _sprepost(A::AbstractSparseMatrix, B::AbstractSparseMatrix) = kron(transpose(B), A)
+_liouvillian(H::AbstractMatrix, Id::AbstractMatrix) = -1im * (_spre(H, Id) - _spost(H, Id))
+
+# (if input is AbstractSciMLOperator)
+_spre(A::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spre(A.A, Id))
+_spre(A::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(A.λ, _spre(A.L, Id))
+_spre(A::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spre(op, Id), +, A.ops)
+_spost(B::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spost(B.A, Id))
+_spost(B::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(B.λ, _spost(B.L, Id))
+_spost(B::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spost(op, Id), +, B.ops)
+_liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
+_liouvillian(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian(H.L, Id))
+_liouvillian(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian(op, Id), +, H.ops)
+# TODO: support `_sprepost`, `sprepost`, and `lindblad_dissipator` for AbstractSciMLOperator (allow c_ops with Vector{QobjEvo})
 
 @doc raw"""
-    spre(A::QuantumObject, Id_cache=I(size(A,1)))
+    spre(A::AbstractQuantumObject, Id_cache=I(size(A,1)))
 
 Returns the [`SuperOperator`](@ref) form of `A` acting on the left of the density matrix operator: ``\mathcal{O} \left(\hat{A}\right) \left[ \hat{\rho} \right] = \hat{A} \hat{\rho}``.
 
@@ -27,14 +40,13 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 ```
 (see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
-The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when
-the same function is applied multiple times with a known Hilbert space dimension.
+The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
 """
-spre(A::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}, Id_cache = I(size(A, 1))) where {T} =
-    QuantumObject(_spre(A.data, Id_cache), SuperOperator, A.dims)
+spre(A::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(A, 1))) where {DT} =
+    get_typename_wrapper(A)(_spre(A.data, Id_cache), SuperOperator, A.dims)
 
 @doc raw"""
-    spost(B::QuantumObject, Id_cache=I(size(B,1)))
+    spost(B::AbstractQuantumObject, Id_cache=I(size(B,1)))
 
 Returns the [`SuperOperator`](@ref) form of `B` acting on the right of the density matrix operator: ``\mathcal{O} \left(\hat{B}\right) \left[ \hat{\rho} \right] = \hat{\rho} \hat{B}``.
 
@@ -45,11 +57,10 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 ```
 (see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
-The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when
-the same function is applied multiple times with a known Hilbert space dimension.
+The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
 """
-spost(B::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}, Id_cache = I(size(B, 1))) where {T} =
-    QuantumObject(_spost(B.data, Id_cache), SuperOperator, B.dims)
+spost(B::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(B, 1))) where {DT} =
+    get_typename_wrapper(B)(_spost(B.data, Id_cache), SuperOperator, B.dims)
 
 @doc raw"""
     sprepost(A::QuantumObject, B::QuantumObject)
@@ -84,21 +95,21 @@ Returns the Lindblad [`SuperOperator`](@ref) defined as
 \hat{O}^\dagger \hat{O} \hat{\rho} - \hat{\rho} \hat{O}^\dagger \hat{O} \right)
 ```
 
-The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when
-the same function is applied multiple times with a known Hilbert space dimension.
+The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
 
-See also [`spre`](@ref) and [`spost`](@ref).
+See also [`spre`](@ref), [`spost`](@ref), and [`sprepost`](@ref).
 """
 function lindblad_dissipator(O::QuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(O, 1))) where {DT}
     Od_O = O' * O
-    return sprepost(O, O') - spre(Od_O, Id_cache) / 2 - spost(Od_O, Id_cache) / 2
+    return sprepost(O, O') - (spre(Od_O, Id_cache) + spost(Od_O, Id_cache)) / 2
 end
+# TODO: suppport collapse operator given as QobjEvo-type
 
 # It is already a SuperOperator
 lindblad_dissipator(O::QuantumObject{DT,SuperOperatorQuantumObject}, Id_cache = nothing) where {DT} = O
 
 @doc raw"""
-    liouvillian(H::QuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dims)))
+    liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dims)))
 
 Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
 
@@ -117,20 +128,18 @@ The optional argument `Id_cache` can be used to pass a precomputed identity matr
 See also [`spre`](@ref), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
 """
 function liouvillian(
-    H::QuantumObject{MT1,OpType1},
+    H::AbstractQuantumObject{DT,OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     Id_cache = I(prod(H.dims)),
-) where {MT1<:AbstractMatrix,OpType1<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
     if !(c_ops isa Nothing)
-        for c_op in c_ops
-            L += lindblad_dissipator(c_op, Id_cache)
-        end
+        L += mapreduce(lindblad_dissipator, +, c_ops)
     end
     return L
 end
 
-liouvillian(H::QuantumObject{<:AbstractMatrix,OperatorQuantumObject}, Id_cache::Diagonal = I(prod(H.dims))) =
-    -1im * (spre(H, Id_cache) - spost(H, Id_cache))
+liouvillian(H::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache::Diagonal = I(prod(H.dims))) where {DT} =
+    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator, H.dims)
 
-liouvillian(H::QuantumObject{<:AbstractMatrix,SuperOperatorQuantumObject}, Id_cache::Diagonal) = H
+liouvillian(H::AbstractQuantumObject{DT,SuperOperatorQuantumObject}, Id_cache::Diagonal) where {DT} = H
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index f191921ae..76714c579 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -42,23 +42,7 @@ function _generate_mesolve_kwargs(e_ops, progress_bar::Val{false}, tlist, kwargs
 end
 
 _mesolve_make_L_QobjEvo(H::QuantumObject, c_ops) = QobjEvo(liouvillian(H, c_ops); type = SuperOperator)
-function _mesolve_make_L_QobjEvo(H::Tuple, c_ops)
-    c_ops isa Nothing && return QobjEvo(H; type = SuperOperator, f = liouvillian)
-    return QobjEvo((H..., mapreduce(op -> lindblad_dissipator(op), +, c_ops)); type = SuperOperator, f = liouvillian)
-end
-_mesolve_make_L_QobjEvo(H::QuantumObjectEvolution{DT,OperatorQuantumObject}, c_ops) where {DT<:AbstractSciMLOperator} =
-    throw(
-        ArgumentError(
-            "This function does not support the data type of time-dependent Operator `H` currently. Try to provide `H` as a time-dependent SuperOperator or Tuple instead.",
-        ),
-    )
-function _mesolve_make_L_QobjEvo(
-    H::QuantumObjectEvolution{DT,SuperOperatorQuantumObject},
-    c_ops,
-) where {DT<:AbstractSciMLOperator}
-    c_ops isa Nothing && return H
-    return H + QobjEvo((mapreduce(op -> lindblad_dissipator(op), +, c_ops)))
-end
+_mesolve_make_L_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}, c_ops) = liouvillian(QobjEvo(H), c_ops)
 
 @doc raw"""
     mesolveProblem(
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 90284b50c..9297cbedd 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -1,6 +1,6 @@
 export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionSSESol
 
-export liouvillian, liouvillian_floquet, liouvillian_generalized
+export liouvillian_floquet, liouvillian_generalized
 
 const DEFAULT_ODE_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep = false, save_end = true)
 const DEFAULT_SDE_SOLVER_OPTIONS = (abstol = 1e-2, reltol = 1e-2, save_everystep = false, save_end = true)
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index b472e0fef..be24e3bfb 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -65,10 +65,15 @@
     end
 
     @testset "Operator and SuperOperator" begin
+        N = 10
+        A = Qobj(rand(ComplexF64, N, N))
+        B = Qobj(rand(ComplexF64, N, N))
+        ρ = rand_dm(N) # random density matrix
+        @test mat2vec(A * ρ * B) ≈ spre(A) * spost(B) * mat2vec(ρ) ≈ sprepost(A, B) * mat2vec(ρ) # we must make sure this equality holds !
+
         a = sprand(ComplexF64, 100, 100, 0.1)
         a2 = Qobj(a)
         a3 = Qobj(a, type = SuperOperator)
-
         @test isket(a2) == false
         @test isbra(a2) == false
         @test isoper(a2) == true
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 4d41df619..4cad02526 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -160,25 +160,37 @@
     #     @test X_warn == X3
     # end
 
-    @testset "Time Dependent Operators" begin
+    @testset "Time Dependent Operators and SuperOperators" begin
         N = 10
         a = destroy(N)
         coef1(p, t) = exp(-1im * p.ω1 * t)
         coef2(p, t) = sin(p.ω2 * t)
+        t = rand()
+        p = (ω1 = rand(), ω2 = rand())
+
+        # Operator
+        H_td = QobjEvo(((a, coef1), a' * a, (a', coef2)))
+        H_ti = coef1(p, t) * a + a' * a + coef2(p, t) * a'
+        ψ = rand_ket(N)
+        @test H_td(p, t) ≈ H_ti
+        @test H_td(ψ, p, t) ≈ H_ti * ψ
+        @test isconstant(a) == true
+        @test isconstant(H_td) == false
+        @test isconstant(QobjEvo(a)) == true
+        @test isoper(H_td) == true
+
+        # SuperOperator
+        c_ops = [sqrt(rand()) * a]
+        L_td = liouvillian(H_td, c_ops)
+        L_td2 = -1im * spre(H_td) + 1im * spost(H_td) + lindblad_dissipator(c_ops[1])
+        ρvec = mat2vec(rand_dm(N))
+        @test L_td(p, t) ≈ L_td2(p, t)
+        @test L_td(ρvec, p, t) ≈ L_td2(ρvec, p, t)
+        @test isconstant(L_td) == false
+        @test issuper(L_td) == true
 
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])
-
-        op1 = QobjEvo(((a, coef1), a' * a, (a', coef2)))
-        op1 = QobjEvo(((a, coef1), a' * a, (a', coef2)))
-
-        p = (ω1 = 1, ω2 = 2)
-        @test op1(p, 0.1) ≈ coef1(p, 0.1) * a + a' * a + coef2(p, 0.1) * a'
-
-        ψ = fock(N, 1)
-        @test op1(ψ, p, 0.1) ≈ (coef1(p, 0.1) * a + a' * a + coef2(p, 0.1) * a') * ψ
-
-        @test isconstant(a) == true
-        @test isconstant(op1) == false
-        @test isconstant(Qobj(a)) == true
+        @test_throws ArgumentError H_td(ρvec, p, t)
+        @test_throws ArgumentError L_td(ψ, p, t)
     end
 end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 666aa38ab..fd7a33f91 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -205,18 +205,9 @@
         # @test sol_sse.expect ≈ sol_sse_td.expect atol = 1e-2 * length(tlist)
 
         H_td2 = QobjEvo(H_td)
-        L_td = QobjEvo(H_td, type = SuperOperator, f = liouvillian)
+        L_td = liouvillian(H_td2)
 
         sol_se_td2 = sesolve(H_td2, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p)
-        @test_throws ArgumentError mesolve(
-            H_td2,
-            ψ0,
-            tlist,
-            c_ops,
-            e_ops = e_ops,
-            progress_bar = Val(false),
-            params = p,
-        )
         sol_me_td2 = mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
         sol_mc_td2 = mcsolve(
             H_td2,
@@ -253,6 +244,7 @@
             @inferred mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, saveat = tlist, progress_bar = Val(false))
             @inferred mesolve(H, ψ0, tlist, (a, a'), e_ops = (a' * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple
             @inferred mesolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
+            @inferred mesolve(H_td2, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
             @inferred mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
         end