diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index e45dcbbf..66b84572 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -71,6 +71,7 @@ jobs:
       - uses: julia-actions/julia-runtest@v1
         env:
           GROUP: ${{ matrix.group }}
+          JULIA_NUM_THREADS: auto
       - uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,ext
diff --git a/docs/src/api.md b/docs/src/api.md
index dae7f72a..554d1e8b 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -13,6 +13,7 @@ Pages = ["api.md"]
 ## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
 
 ```@docs
+AbstractQuantumObject
 BraQuantumObject
 Bra
 KetQuantumObject
@@ -26,10 +27,15 @@ OperatorKet
 SuperOperatorQuantumObject
 SuperOperator
 QuantumObject
-OperatorSum
+QuantumObjectEvolution
 size
 eltype
 length
+```
+
+## [Qobj boolean functions](@id doc-API:Qobj-boolean-functions)
+
+```@docs
 isbra
 isket
 isoper
@@ -40,6 +46,7 @@ LinearAlgebra.ishermitian
 LinearAlgebra.issymmetric
 LinearAlgebra.isposdef
 isunitary
+isconstant
 ```
 
 ## [Qobj arithmetic and attributes](@id doc-API:Qobj-arithmetic-and-attributes)
@@ -154,6 +161,7 @@ lindblad_dissipator
 ## [Synonyms of functions for Qobj](@id doc-API:Synonyms-of-functions-for-Qobj)
 ```@docs
 Qobj
+QobjEvo
 shape
 isherm
 trans
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 9c2920df..549df5b0 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -21,6 +21,7 @@ import SciMLBase:
     reinit!,
     remake,
     u_modified!,
+    ODEFunction,
     ODEProblem,
     SDEProblem,
     EnsembleProblem,
@@ -32,13 +33,21 @@ import SciMLBase:
     ContinuousCallback,
     DiscreteCallback
 import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1
-import SciMLOperators: MatrixOperator
+import SciMLOperators:
+    AbstractSciMLOperator,
+    MatrixOperator,
+    ScalarOperator,
+    IdentityOperator,
+    cache_operator,
+    update_coefficients!,
+    concretize,
+    isconstant
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
 import OrdinaryDiffEqTsit5: Tsit5
-import DiffEqNoiseProcess: RealWienerProcess
+import DiffEqNoiseProcess: RealWienerProcess!
 
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
@@ -62,7 +71,9 @@ include("progress_bar.jl")
 include("linear_maps.jl")
 
 # Quantum Object
+include("qobj/quantum_object_base.jl")
 include("qobj/quantum_object.jl")
+include("qobj/quantum_object_evo.jl")
 include("qobj/boolean_functions.jl")
 include("qobj/arithmetic_and_attributes.jl")
 include("qobj/eigsolve.jl")
@@ -71,7 +82,6 @@ include("qobj/states.jl")
 include("qobj/operators.jl")
 include("qobj/superoperators.jl")
 include("qobj/synonyms.jl")
-include("qobj/operator_sum.jl")
 
 # time evolution
 include("time_evolution/time_evolution.jl")
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 29b872de..96947bb2 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -36,84 +36,80 @@ end
 
 for op in (:(+), :(-), :(*))
     @eval begin
-        function LinearAlgebra.$op(
-            A::QuantumObject{<:AbstractArray{T1},OpType},
-            B::QuantumObject{<:AbstractArray{T2},OpType},
-        ) where {T1,T2,OpType<:QuantumObjectType}
-            A.dims != B.dims &&
-                throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-            return QuantumObject($(op)(A.data, B.data), A.type, A.dims)
+        function LinearAlgebra.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
+            check_dims(A, B)
+            QType = promote_op_type(A, B)
+            return QType($(op)(A.data, B.data), A.type, A.dims)
         end
-        LinearAlgebra.$op(A::QuantumObject{<:AbstractArray{T}}) where {T} = QuantumObject($(op)(A.data), A.type, A.dims)
+        LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dims)
 
-        LinearAlgebra.$op(n::T1, A::QuantumObject{<:AbstractArray{T2}}) where {T1<:Number,T2} =
-            QuantumObject($(op)(n * I, A.data), A.type, A.dims)
-        LinearAlgebra.$op(A::QuantumObject{<:AbstractArray{T1}}, n::T2) where {T1,T2<:Number} =
-            QuantumObject($(op)(A.data, n * I), A.type, A.dims)
+        LinearAlgebra.$op(n::T, A::AbstractQuantumObject) where {T<:Number} =
+            get_typename_wrapper(A)($(op)(n * I, A.data), A.type, A.dims)
+        LinearAlgebra.$op(A::AbstractQuantumObject, n::T) where {T<:Number} =
+            get_typename_wrapper(A)($(op)(A.data, n * I), A.type, A.dims)
     end
 end
 
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,OperatorQuantumObject},
+    B::QuantumObject{DT2,KetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Ket, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},BraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,BraQuantumObject},
+    B::AbstractQuantumObject{DT2,OperatorQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Bra, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},KetQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},BraQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,KetQuantumObject},
+    B::QuantumObject{DT2,BraQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Operator, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},BraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,BraQuantumObject},
+    B::QuantumObject{DT2,KetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
+    B::QuantumObject{DT2,OperatorQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},OperatorBraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorKetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,OperatorBraQuantumObject},
+    B::QuantumObject{DT2,OperatorKetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorKetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
+    B::QuantumObject{DT2,OperatorKetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, OperatorKet, A.dims)
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{<:AbstractArray{T1},OperatorBraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
+    B::AbstractQuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
 ) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, OperatorBra, A.dims)
 end
 
-LinearAlgebra.:(^)(A::QuantumObject{<:AbstractArray{T}}, n::T1) where {T,T1<:Number} =
-    QuantumObject(^(A.data, n), A.type, A.dims)
-LinearAlgebra.:(/)(A::QuantumObject{<:AbstractArray{T}}, n::T1) where {T,T1<:Number} =
-    QuantumObject(/(A.data, n), A.type, A.dims)
+LinearAlgebra.:(^)(A::QuantumObject{DT}, n::T) where {DT,T<:Number} = QuantumObject(^(A.data, n), A.type, A.dims)
+LinearAlgebra.:(/)(A::AbstractQuantumObject{DT}, n::T) where {DT,T<:Number} =
+    get_typename_wrapper(A)(A.data / n, A.type, A.dims)
 
 @doc raw"""
     dot(A::QuantumObject, B::QuantumObject)
@@ -125,93 +121,91 @@ Note that `A` and `B` should be [`Ket`](@ref) or [`OperatorKet`](@ref)
 `A ⋅ B` (where `⋅` can be typed by tab-completing `\cdot` in the REPL) is a synonym for `dot(A, B)`
 """
 function LinearAlgebra.dot(
-    A::QuantumObject{<:AbstractArray{T1},OpType},
-    B::QuantumObject{<:AbstractArray{T2},OpType},
-) where {T1<:Number,T2<:Number,OpType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    A::QuantumObject{DT1,OpType},
+    B::QuantumObject{DT2,OpType},
+) where {DT1,DT2,OpType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
     A.dims != B.dims && throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(A.data, B.data)
 end
 
 @doc raw"""
-    dot(i::QuantumObject, A::QuantumObject j::QuantumObject)
+    dot(i::QuantumObject, A::AbstractQuantumObject j::QuantumObject)
 
-Compute the generalized dot product `dot(i, A*j)` between three [`QuantumObject`](@ref): ``\langle i | \hat{A} | j \rangle``
+Compute the generalized dot product `dot(i, A*j)` between a [`AbstractQuantumObject`](@ref) and two [`QuantumObject`](@ref) (`i` and `j`), namely ``\langle i | \hat{A} | j \rangle``.
 
 Supports the following inputs:
 - `A` is in the type of [`Operator`](@ref), with `i` and `j` are both [`Ket`](@ref).
 - `A` is in the type of [`SuperOperator`](@ref), with `i` and `j` are both [`OperatorKet`](@ref)
 """
 function LinearAlgebra.dot(
-    i::QuantumObject{<:AbstractArray{T1},KetQuantumObject},
-    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-    j::QuantumObject{<:AbstractArray{T3},KetQuantumObject},
-) where {T1<:Number,T2<:Number,T3<:Number}
+    i::QuantumObject{DT1,KetQuantumObject},
+    A::AbstractQuantumObject{DT2,OperatorQuantumObject},
+    j::QuantumObject{DT3,KetQuantumObject},
+) where {DT1,DT2,DT3}
     ((i.dims != A.dims) || (A.dims != j.dims)) &&
         throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 function LinearAlgebra.dot(
-    i::QuantumObject{<:AbstractArray{T1},OperatorKetQuantumObject},
-    A::QuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
-    j::QuantumObject{<:AbstractArray{T3},OperatorKetQuantumObject},
-) where {T1<:Number,T2<:Number,T3<:Number}
+    i::QuantumObject{DT1,OperatorKetQuantumObject},
+    A::AbstractQuantumObject{DT2,SuperOperatorQuantumObject},
+    j::QuantumObject{DT3,OperatorKetQuantumObject},
+) where {DT1,DT2,DT3}
     ((i.dims != A.dims) || (A.dims != j.dims)) &&
         throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 
 @doc raw"""
-    conj(A::QuantumObject)
+    conj(A::AbstractQuantumObject)
 
-Return the element-wise complex conjugation of the [`QuantumObject`](@ref).
+Return the element-wise complex conjugation of the [`AbstractQuantumObject`](@ref).
 """
-Base.conj(A::QuantumObject{<:AbstractArray{T}}) where {T} = QuantumObject(conj(A.data), A.type, A.dims)
+Base.conj(A::AbstractQuantumObject) = get_typename_wrapper(A)(conj(A.data), A.type, A.dims)
 
 @doc raw"""
-    transpose(A::QuantumObject)
+    transpose(A::AbstractQuantumObject)
 
-Lazy matrix transpose of the [`QuantumObject`](@ref).
+Lazy matrix transpose of the [`AbstractQuantumObject`](@ref).
 """
 LinearAlgebra.transpose(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(transpose(A.data), A.type, A.dims)
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    get_typename_wrapper(A)(transpose(A.data), A.type, A.dims)
 
 @doc raw"""
     A'
-    adjoint(A::QuantumObject)
+    adjoint(A::AbstractQuantumObject)
 
-Lazy adjoint (conjugate transposition) of the [`QuantumObject`](@ref)
+Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 
 Note that `A'` is a synonym for `adjoint(A)`
 """
 LinearAlgebra.adjoint(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(adjoint(A.data), A.type, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},KetQuantumObject}) where {T} =
-    QuantumObject(adjoint(A.data), Bra, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},BraQuantumObject}) where {T} =
-    QuantumObject(adjoint(A.data), Ket, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},OperatorKetQuantumObject}) where {T} =
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    get_typename_wrapper(A)(adjoint(A.data), A.type, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,KetQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Bra, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,BraQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Ket, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorKetQuantumObject}) where {DT} =
     QuantumObject(adjoint(A.data), OperatorBra, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},OperatorBraQuantumObject}) where {T} =
+LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorBraQuantumObject}) where {DT} =
     QuantumObject(adjoint(A.data), OperatorKet, A.dims)
 
 @doc raw"""
-    inv(A::QuantumObject)
+    inv(A::AbstractQuantumObject)
 
-Matrix inverse of the [`QuantumObject`](@ref)
+Matrix inverse of the [`AbstractQuantumObject`](@ref). If `A` is a [`QuantumObjectEvolution`](@ref), the inverse is computed at the last computed time.
 """
 LinearAlgebra.inv(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     QuantumObject(sparse(inv(Matrix(A.data))), A.type, A.dims)
 
 LinearAlgebra.Hermitian(
-    A::QuantumObject{<:AbstractArray{T},OpType},
+    A::QuantumObject{DT,OpType},
     uplo::Symbol = :U,
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     QuantumObject(Hermitian(A.data, uplo), A.type, A.dims)
 
 @doc raw"""
@@ -436,8 +430,8 @@ Matrix sine of [`QuantumObject`](@ref), defined as
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
 LinearAlgebra.sin(
-    A::QuantumObject{<:AbstractMatrix{T},ObjType},
-) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) - exp(-1im * A)) / 2im
+    A::QuantumObject{DT,ObjType},
+) where {DT,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) - exp(-1im * A)) / 2im
 
 @doc raw"""
     cos(A::QuantumObject)
@@ -449,8 +443,8 @@ Matrix cosine of [`QuantumObject`](@ref), defined as
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
 LinearAlgebra.cos(
-    A::QuantumObject{<:AbstractMatrix{T},ObjType},
-) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) + exp(-1im * A)) / 2
+    A::QuantumObject{DT,ObjType},
+) where {DT,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) + exp(-1im * A)) / 2
 
 @doc raw"""
     diag(A::QuantumObject, k::Int=0)
@@ -659,11 +653,11 @@ tidyup!(A::AbstractArray{T}, tol::T2 = 1e-14) where {T,T2<:Real} =
     @. A = T(abs(real(A)) > tol) * real(A) + 1im * T(abs(imag(A)) > tol) * imag(A)
 
 @doc raw"""
-    get_data(A::QuantumObject)
+    get_data(A::AbstractQuantumObject)
 
-Returns the data of a QuantumObject.
+Returns the data of a [`AbstractQuantumObject`](@ref).
 """
-get_data(A::QuantumObject) = A.data
+get_data(A::AbstractQuantumObject) = A.data
 
 @doc raw"""
     get_coherence(ψ::QuantumObject)
diff --git a/src/qobj/boolean_functions.jl b/src/qobj/boolean_functions.jl
index 4ec34767..cccea532 100644
--- a/src/qobj/boolean_functions.jl
+++ b/src/qobj/boolean_functions.jl
@@ -3,74 +3,82 @@ All boolean functions for checking the data or type in `QuantumObject`
 =#
 
 export isket, isbra, isoper, isoperbra, isoperket, issuper
-export isunitary
+export isunitary, isconstant
 
 @doc raw"""
-    isbra(A::QuantumObject)
+    isbra(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`BraQuantumObject`](@ref).
+Checks if the [`QuantumObject`](@ref) `A` is a [`BraQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isbra(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} = OpType <: BraQuantumObject
+isbra(A::QuantumObject) = isbra(typeof(A))
+isbra(::Type{QuantumObject{DT,BraQuantumObject,N}}) where {DT,N} = true
+isbra(A) = false # default case
 
 @doc raw"""
-    isket(A::QuantumObject)
+    isket(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`KetQuantumObject`](@ref).
+Checks if the [`QuantumObject`](@ref) `A` is a [`KetQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isket(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} = OpType <: KetQuantumObject
+isket(A::QuantumObject) = isket(typeof(A))
+isket(::Type{QuantumObject{DT,KetQuantumObject,N}}) where {DT,N} = true
+isket(A) = false # default case
 
 @doc raw"""
-    isoper(A::QuantumObject)
+    isoper(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorQuantumObject`](@ref).
+Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`OperatorQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoper(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =
-    OpType <: OperatorQuantumObject
+isoper(A::AbstractQuantumObject) = isoper(typeof(A))
+isoper(::Type{<:AbstractQuantumObject{DT,OperatorQuantumObject,N}}) where {DT,N} = true
+isoper(A) = false # default case
 
 @doc raw"""
-    isoperbra(A::QuantumObject)
+    isoperbra(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorBraQuantumObject`](@ref).
+Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorBraQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoperbra(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =
-    OpType <: OperatorBraQuantumObject
+isoperbra(A::QuantumObject) = isoperbra(typeof(A))
+isoperbra(::Type{QuantumObject{DT,OperatorBraQuantumObject,N}}) where {DT,N} = true
+isoperbra(A) = false # default case
 
 @doc raw"""
-    isoperket(A::QuantumObject)
+    isoperket(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorKetQuantumObject`](@ref).
+Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorKetQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoperket(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =
-    OpType <: OperatorKetQuantumObject
+isoperket(A::QuantumObject) = isoperket(typeof(A))
+isoperket(::Type{QuantumObject{DT,OperatorKetQuantumObject,N}}) where {DT,N} = true
+isoperket(A) = false # default case
 
 @doc raw"""
-    issuper(A::QuantumObject)
+    issuper(A)
 
-Checks if the [`QuantumObject`](@ref) `A` is a [`SuperOperatorQuantumObject`](@ref).
+Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`SuperOperatorQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-issuper(A::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:QuantumObjectType} =
-    OpType <: SuperOperatorQuantumObject
+issuper(A::AbstractQuantumObject) = issuper(typeof(A))
+issuper(::Type{<:AbstractQuantumObject{DT,SuperOperatorQuantumObject,N}}) where {DT,N} = true
+issuper(A) = false # default case
 
 @doc raw"""
-    ishermitian(A::QuantumObject)
+    ishermitian(A::AbstractQuantumObject)
 
-Test whether the [`QuantumObject`](@ref) is Hermitian.
+Test whether the [`AbstractQuantumObject`](@ref) is Hermitian.
 """
-LinearAlgebra.ishermitian(A::QuantumObject{<:AbstractArray{T}}) where {T} = ishermitian(A.data)
+LinearAlgebra.ishermitian(A::AbstractQuantumObject) = ishermitian(A.data)
 
 @doc raw"""
-    issymmetric(A::QuantumObject)
+    issymmetric(A::AbstractQuantumObject)
 
-Test whether the [`QuantumObject`](@ref) is symmetric.
+Test whether the [`AbstractQuantumObject`](@ref) is symmetric.
 """
-LinearAlgebra.issymmetric(A::QuantumObject{<:AbstractArray{T}}) where {T} = issymmetric(A.data)
+LinearAlgebra.issymmetric(A::AbstractQuantumObject) = issymmetric(A.data)
 
 @doc raw"""
-    isposdef(A::QuantumObject)
+    isposdef(A::AbstractQuantumObject)
 
-Test whether the [`QuantumObject`](@ref) is positive definite (and Hermitian) by trying to perform a Cholesky factorization of `A`.
+Test whether the [`AbstractQuantumObject`](@ref) is positive definite (and Hermitian) by trying to perform a Cholesky factorization of `A`.
 """
-LinearAlgebra.isposdef(A::QuantumObject{<:AbstractArray{T}}) where {T} = isposdef(A.data)
+LinearAlgebra.isposdef(A::AbstractQuantumObject) = isposdef(A.data)
 
 @doc raw"""
     isunitary(U::QuantumObject; kwargs...)
@@ -81,3 +89,10 @@ Note that all the keyword arguments will be passed to `Base.isapprox`.
 """
 isunitary(U::QuantumObject{<:AbstractArray{T}}; kwargs...) where {T} =
     isoper(U) ? isapprox(U.data * U.data', I(size(U, 1)); kwargs...) : false
+
+@doc raw"""
+    isconstant(A::AbstractQuantumObject)
+
+Test whether the [`AbstractQuantumObject`](@ref) `A` is constant in time. For a [`QuantumObject`](@ref), this function returns `true`, while for a [`QuantumObjectEvolution`](@ref), this function returns `true` if the operator is contant in time.
+"""
+isconstant(A::AbstractQuantumObject) = isconstant(A.data)
diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl
index 483b82a5..0c972a3f 100644
--- a/src/qobj/eigsolve.jl
+++ b/src/qobj/eigsolve.jl
@@ -3,7 +3,8 @@ Eigen solvers and results for QuantumObject
 =#
 
 export EigsolveResult
-export eigenenergies, eigenstates, eigsolve, eigsolve_al
+export eigenenergies, eigenstates, eigsolve
+export eigsolve_al
 
 @doc raw"""
     struct EigsolveResult{T1<:Vector{<:Number}, T2<:AbstractMatrix{<:Number}, ObjType<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject},N}
@@ -322,33 +323,34 @@ function eigsolve(
 end
 
 @doc raw"""
-    eigsolve_al(H::QuantumObject,
-        T::Real, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        H_t::Union{Nothing,Function}=nothing,
-        params::NamedTuple=NamedTuple(),
-        ρ0::Union{Nothing, AbstractMatrix} = nothing,
-        k::Int=1,
-        krylovdim::Int=min(10, size(H, 1)),
-        maxiter::Int=200,
-        eigstol::Real=1e-6,
-        kwargs...)
+    eigsolve_al(
+        H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
+        T::Real,
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        params::NamedTuple = NamedTuple(),
+        ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
+        k::Int = 1,
+        krylovdim::Int = min(10, size(H, 1)),
+        maxiter::Int = 200,
+        eigstol::Real = 1e-6,
+        kwargs...,
+    )
 
 Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnoldi-Lindblad method.
 
 # Arguments
-- `H`: The Hamiltonian (or directly the Liouvillian) of the system.
-- `T`: The time at which to evaluate the time evolution
+- `H`: The Hamiltonian (or directly the Liouvillian) of the system. It can be a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a tuple of the form supported by [`mesolve`](@ref).
+- `T`: The time at which to evaluate the time evolution.
 - `c_ops`: A vector of collapse operators. Default is `nothing` meaning the system is closed.
-- `alg`: The differential equation solver algorithm
-- `H_t`: A function `H_t(t)` that returns the additional term at time `t`
-- `params`: A dictionary of additional parameters
-- `ρ0`: The initial density matrix. If not specified, a random density matrix is used
-- `k`: The number of eigenvalues to compute
-- `krylovdim`: The dimension of the Krylov subspace
-- `maxiter`: The maximum number of iterations for the eigsolver
-- `eigstol`: The tolerance for the eigsolver
-- `kwargs`: Additional keyword arguments passed to the differential equation solver
+- `alg`: The differential equation solver algorithm. Default is `Tsit5()`.
+- `params`: A `NamedTuple` containing the parameters of the system.
+- `ρ0`: The initial density matrix. If not specified, a random density matrix is used.
+- `k`: The number of eigenvalues to compute.
+- `krylovdim`: The dimension of the Krylov subspace.
+- `maxiter`: The maximum number of iterations for the eigsolver.
+- `eigstol`: The tolerance for the eigsolver.
+- `kwargs`: Additional keyword arguments passed to the differential equation solver.
 
 # Notes
 - For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
@@ -361,11 +363,10 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 - [1] Minganti, F., & Huybrechts, D. (2022). Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of time-independent and Floquet open quantum systems. Quantum, 6, 649.
 """
 function eigsolve_al(
-    H::QuantumObject{MT1,HOpType},
+    H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
     T::Real,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    H_t::Union{Nothing,Function} = nothing,
     params::NamedTuple = NamedTuple(),
     ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
     k::Int = 1,
@@ -373,14 +374,12 @@ function eigsolve_al(
     maxiter::Int = 200,
     eigstol::Real = 1e-6,
     kwargs...,
-) where {MT1<:AbstractMatrix,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
-    L = liouvillian(H, c_ops)
+) where {DT1,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
     prob = mesolveProblem(
-        L,
-        QuantumObject(ρ0, dims = H.dims),
-        [0, T];
-        alg = alg,
-        H_t = H_t,
+        L_evo,
+        QuantumObject(ρ0, type = Operator, dims = H.dims),
+        [zero(T), T];
         params = params,
         progress_bar = Val(false),
         kwargs...,
@@ -389,9 +388,9 @@ function eigsolve_al(
 
     # prog = ProgressUnknown(desc="Applications:", showspeed = true, enabled=progress)
 
-    Lmap = ArnoldiLindbladIntegratorMap(eltype(MT1), size(L), integrator)
+    Lmap = ArnoldiLindbladIntegratorMap(eltype(DT1), size(L_evo), integrator)
 
-    res = _eigsolve(Lmap, mat2vec(ρ0), L.type, L.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
+    res = _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
     # finish!(prog)
 
     vals = similar(res.values)
@@ -399,7 +398,7 @@ function eigsolve_al(
 
     for i in eachindex(res.values)
         vec = view(res.vectors, :, i)
-        vals[i] = dot(vec, L.data, vec)
+        vals[i] = dot(vec, L_evo.data, vec)
         @. vecs[:, i] = vec * exp(-1im * angle(vec[1]))
     end
 
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index de729d08..ef98af4d 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -17,7 +17,7 @@ ket2dm(ψ::QuantumObject{<:AbstractArray{T},KetQuantumObject}) where {T} = ψ *
 ket2dm(ρ::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}) where {T} = ρ
 
 @doc raw"""
-    expect(O::QuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})
+    expect(O::AbstractQuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})
 
 Expectation value of the [`Operator`](@ref) `O` with the state `ψ`. The state can be a [`Ket`](@ref), [`Bra`](@ref) or [`Operator`](@ref).
 
@@ -44,15 +44,15 @@ julia> expect(Hermitian(a' * a), ψ) |> round
 ```
 """
 function expect(
-    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ψ::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-) where {T1,T2}
+    O::AbstractQuantumObject{DT1,OperatorQuantumObject},
+    ψ::QuantumObject{DT2,KetQuantumObject},
+) where {DT1,DT2}
     return dot(ψ.data, O.data, ψ.data)
 end
 function expect(
-    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ψ::QuantumObject{<:AbstractArray{T2},BraQuantumObject},
-) where {T1,T2}
+    O::AbstractQuantumObject{DT1,OperatorQuantumObject},
+    ψ::QuantumObject{DT2,BraQuantumObject},
+) where {DT1,DT2}
     return expect(O, ψ')
 end
 function expect(
@@ -95,11 +95,9 @@ The function returns a real number if `O` is hermitian, and returns a complex nu
 
 Note that `ψ` can also be given as a list of [`QuantumObject`](@ref), it returns a list of expectation values.
 """
-variance(
-    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ψ::QuantumObject{<:AbstractArray{T2}},
-) where {T1,T2} = expect(O^2, ψ) - expect(O, ψ)^2
-variance(O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, ψ::Vector{<:QuantumObject}) where {T1} =
+variance(O::QuantumObject{DT1,OperatorQuantumObject}, ψ::QuantumObject{DT2}) where {DT1,DT2} =
+    expect(O^2, ψ) - expect(O, ψ)^2
+variance(O::QuantumObject{DT1,OperatorQuantumObject}, ψ::Vector{<:QuantumObject}) where {DT1} =
     expect(O^2, ψ) .- expect(O, ψ) .^ 2
 
 @doc raw"""
@@ -149,7 +147,7 @@ function dense_to_sparse(A::VT, tol::Real = 1e-10) where {VT<:AbstractVector}
 end
 
 @doc raw"""
-    kron(A::QuantumObject, B::QuantumObject, ...)
+    kron(A::AbstractQuantumObject, B::AbstractQuantumObject, ...)
 
 Returns the [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product) ``\hat{A} \otimes \hat{B} \otimes \cdots``.
 
@@ -191,13 +189,15 @@ Quantum Object:   type=Operator   dims=[20, 20]   size=(400, 400)   ishermitian=
 ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦
 ```
 """
-LinearAlgebra.kron(
-    A::QuantumObject{<:AbstractArray{T1},OpType},
-    B::QuantumObject{<:AbstractArray{T2},OpType},
-) where {T1,T2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}} =
-    QuantumObject(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
-LinearAlgebra.kron(A::QuantumObject) = A
-function LinearAlgebra.kron(A::Vector{<:QuantumObject})
+function LinearAlgebra.kron(
+    A::AbstractQuantumObject{DT1,OpType},
+    B::AbstractQuantumObject{DT2,OpType},
+) where {DT1,DT2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
+    QType = promote_op_type(A, B)
+    return QType(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
+end
+LinearAlgebra.kron(A::AbstractQuantumObject) = A
+function LinearAlgebra.kron(A::Vector{<:AbstractQuantumObject})
     @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
     return kron(A...)
 end
diff --git a/src/qobj/operator_sum.jl b/src/qobj/operator_sum.jl
deleted file mode 100644
index 5b239821..00000000
--- a/src/qobj/operator_sum.jl
+++ /dev/null
@@ -1,57 +0,0 @@
-export OperatorSum
-
-@doc raw"""
-    struct OperatorSum
-
-A constructor to represent a sum of operators ``\sum_i c_i \hat{O}_i`` with a list of coefficients ``c_i`` and a list of operators ``\hat{O}_i``.
-
-This is very useful when we have to update only the coefficients, without allocating memory by performing the sum of the operators.
-"""
-struct OperatorSum{CT<:AbstractVector{<:Number},OT<:Union{AbstractVector,Tuple}} <: AbstractQuantumObject
-    coefficients::CT
-    operators::OT
-    function OperatorSum(
-        coefficients::CT,
-        operators::OT,
-    ) where {CT<:AbstractVector{<:Number},OT<:Union{AbstractVector,Tuple}}
-        length(coefficients) == length(operators) ||
-            throw(DimensionMismatch("The number of coefficients must be the same as the number of operators."))
-        # Check if all the operators have the same dimensions
-        size_1 = size(operators[1])
-        mapreduce(x -> size(x) == size_1, &, operators) ||
-            throw(DimensionMismatch("All the operators must have the same dimensions."))
-        T = promote_type(
-            mapreduce(x -> eltype(x), promote_type, operators),
-            mapreduce(eltype, promote_type, coefficients),
-        )
-        coefficients2 = T.(coefficients)
-        CT2 = typeof(coefficients2)
-        return new{CT2,OT}(coefficients2, operators)
-    end
-end
-
-Base.size(A::OperatorSum) = size(A.operators[1])
-Base.size(A::OperatorSum, inds...) = size(A.operators[1], inds...)
-Base.length(A::OperatorSum) = length(A.operators[1])
-Base.copy(A::OperatorSum) = OperatorSum(copy(A.coefficients), copy(A.operators))
-Base.deepcopy(A::OperatorSum) = OperatorSum(deepcopy(A.coefficients), deepcopy(A.operators))
-
-function update_coefficients!(A::OperatorSum, coefficients)
-    length(A.coefficients) == length(coefficients) ||
-        throw(DimensionMismatch("The number of coefficients must be the same as the number of operators."))
-    return A.coefficients .= coefficients
-end
-
-@inline function LinearAlgebra.mul!(y::AbstractVector{T}, A::OperatorSum, x::AbstractVector, α, β) where {T}
-    # Note that β is applied only to the first term
-    mul!(y, A.operators[1], x, α * A.coefficients[1], β)
-    @inbounds for i in 2:length(A.operators)
-        A.coefficients[i] == 0 && continue
-        mul!(y, A.operators[i], x, α * A.coefficients[i], 1)
-    end
-    return y
-end
-
-function liouvillian(A::OperatorSum, Id_cache = I(prod(A.operators[1].dims)))
-    return OperatorSum(A.coefficients, liouvillian.(A.operators, Ref(Id_cache)))
-end
diff --git a/src/qobj/quantum_object.jl b/src/qobj/quantum_object.jl
index f55244b2..84e34785 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -1,114 +1,11 @@
 #=
-This file defines:
-    1. the QuantumObject (Qobj) structure
-    2. all the type structures for QuantumObject
-Also support for fundamental functions in Julia standard library:
-    - Base: show, length, size, eltype, getindex, setindex!, isequal, :(==), isapprox, Vector, Matrix
+This file defines the QuantumObject (Qobj) structure.
+It also implements the fundamental functions in Julia standard library:
+    - Base: show, real, imag, Vector, Matrix
     - SparseArrays: sparse, nnz, nonzeros, rowvals, droptol!, dropzeros, dropzeros!, SparseVector, SparseMatrixCSC
 =#
 
-export AbstractQuantumObject, QuantumObject
-export QuantumObjectType,
-    BraQuantumObject,
-    KetQuantumObject,
-    OperatorQuantumObject,
-    OperatorBraQuantumObject,
-    OperatorKetQuantumObject,
-    SuperOperatorQuantumObject
-export Bra, Ket, Operator, OperatorBra, OperatorKet, SuperOperator
-
-abstract type AbstractQuantumObject end
-abstract type QuantumObjectType end
-
-@doc raw"""
-    BraQuantumObject <: QuantumObjectType
-
-Constructor representing a bra state ``\langle\psi|``.
-"""
-struct BraQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::BraQuantumObject) = print(io, "Bra")
-
-@doc raw"""
-    const Bra = BraQuantumObject()
-
-A constant representing the type of [`BraQuantumObject`](@ref): a bra state ``\langle\psi|``
-"""
-const Bra = BraQuantumObject()
-
-@doc raw"""
-    KetQuantumObject <: QuantumObjectType
-
-Constructor representing a ket state ``|\psi\rangle``.
-"""
-struct KetQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::KetQuantumObject) = print(io, "Ket")
-
-@doc raw"""
-    const Ket = KetQuantumObject()
-
-A constant representing the type of [`KetQuantumObject`](@ref): a ket state ``|\psi\rangle``
-"""
-const Ket = KetQuantumObject()
-
-@doc raw"""
-    OperatorQuantumObject <: QuantumObjectType
-
-Constructor representing an operator ``\hat{O}``.
-"""
-struct OperatorQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::OperatorQuantumObject) = print(io, "Operator")
-
-@doc raw"""
-    const Operator = OperatorQuantumObject()
-
-A constant representing the type of [`OperatorQuantumObject`](@ref): an operator ``\hat{O}``
-"""
-const Operator = OperatorQuantumObject()
-
-@doc raw"""
-    SuperOperatorQuantumObject <: QuantumObjectType
-
-Constructor representing a super-operator ``\hat{\mathcal{O}}`` acting on vectorized density operator matrices.
-"""
-struct SuperOperatorQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::SuperOperatorQuantumObject) = print(io, "SuperOperator")
-
-@doc raw"""
-    const SuperOperator = SuperOperatorQuantumObject()
-
-A constant representing the type of [`SuperOperatorQuantumObject`](@ref): a super-operator ``\hat{\mathcal{O}}`` acting on vectorized density operator matrices
-"""
-const SuperOperator = SuperOperatorQuantumObject()
-
-@doc raw"""
-    OperatorBraQuantumObject <: QuantumObjectType
-
-Constructor representing a bra state in the [`SuperOperator`](@ref) formalism ``\langle\langle\rho|``.
-"""
-struct OperatorBraQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::OperatorBraQuantumObject) = print(io, "OperatorBra")
-
-@doc raw"""
-    const OperatorBra = OperatorBraQuantumObject()
-
-A constant representing the type of [`OperatorBraQuantumObject`](@ref): a bra state in the [`SuperOperator`](@ref) formalism ``\langle\langle\rho|``.
-"""
-const OperatorBra = OperatorBraQuantumObject()
-
-@doc raw"""
-    OperatorKetQuantumObject <: QuantumObjectType
-
-Constructor representing a ket state in the [`SuperOperator`](@ref) formalism ``|\rho\rangle\rangle``.
-"""
-struct OperatorKetQuantumObject <: QuantumObjectType end
-Base.show(io::IO, ::OperatorKetQuantumObject) = print(io, "OperatorKet")
-
-@doc raw"""
-    const OperatorKet = OperatorKetQuantumObject()
-
-A constant representing the type of [`OperatorKetQuantumObject`](@ref): a ket state in the [`SuperOperator`](@ref) formalism ``|\rho\rangle\rangle``
-"""
-const OperatorKet = OperatorKetQuantumObject()
+export QuantumObject
 
 @doc raw"""
     struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N}
@@ -135,7 +32,7 @@ julia> a isa QuantumObject
 true
 ```
 """
-struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject
+struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject{MT,ObjType,N}
     data::MT
     type::ObjType
     dims::SVector{N,Int}
@@ -215,57 +112,6 @@ function QuantumObject(
     throw(DomainError(size(A), "The size of the array is not compatible with vector or matrix."))
 end
 
-function _check_dims(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
-    _non_static_array_warning("dims", dims)
-    return (all(>(0), dims) && length(dims) > 0) ||
-           throw(DomainError(dims, "The argument dims must be of non-zero length and contain only positive integers."))
-end
-_check_dims(dims::Any) = throw(
-    ArgumentError(
-        "The argument dims must be a Tuple or a StaticVector of non-zero length and contain only positive integers.",
-    ),
-)
-
-function _check_QuantumObject(type::KetQuantumObject, dims, m::Int, n::Int)
-    (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Ket"))
-    (prod(dims) != m) && throw(DimensionMismatch("Ket with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
-function _check_QuantumObject(type::BraQuantumObject, dims, m::Int, n::Int)
-    (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Bra"))
-    (prod(dims) != n) && throw(DimensionMismatch("Bra with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
-function _check_QuantumObject(type::OperatorQuantumObject, dims, m::Int, n::Int)
-    (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with Operator"))
-    (prod(dims) != m) &&
-        throw(DimensionMismatch("Operator with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
-function _check_QuantumObject(type::SuperOperatorQuantumObject, dims, m::Int, n::Int)
-    (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with SuperOperator"))
-    (prod(dims) != sqrt(m)) &&
-        throw(DimensionMismatch("SuperOperator with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
-function _check_QuantumObject(type::OperatorKetQuantumObject, dims, m::Int, n::Int)
-    (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorKet"))
-    (prod(dims) != sqrt(m)) &&
-        throw(DimensionMismatch("OperatorKet with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
-function _check_QuantumObject(type::OperatorBraQuantumObject, dims, m::Int, n::Int)
-    (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorBra"))
-    (prod(dims) != sqrt(n)) &&
-        throw(DimensionMismatch("OperatorBra with dims = $(dims) does not fit the array size = $((m, n))."))
-    return nothing
-end
-
 function QuantumObject(
     A::QuantumObject{<:AbstractArray{T,N}};
     type::ObjType = A.type,
@@ -294,7 +140,7 @@ function Base.show(
     return show(io, MIME("text/plain"), op_data)
 end
 
-function Base.show(io::IO, QO::QuantumObject{<:AbstractArray{T},OpType}) where {T,OpType<:OperatorQuantumObject}
+function Base.show(io::IO, QO::AbstractQuantumObject)
     op_data = QO.data
     println(
         io,
@@ -310,41 +156,6 @@ function Base.show(io::IO, QO::QuantumObject{<:AbstractArray{T},OpType}) where {
     return show(io, MIME("text/plain"), op_data)
 end
 
-@doc raw"""
-    size(A::QuantumObject)
-    size(A::QuantumObject, idx::Int)
-
-Returns a tuple containing each dimensions of the array in the [`QuantumObject`](@ref).
-
-Optionally, you can specify an index (`idx`) to just get the corresponding dimension of the array.
-"""
-Base.size(A::QuantumObject{<:AbstractArray{T}}) where {T} = size(A.data)
-Base.size(A::QuantumObject{<:AbstractArray{T}}, idx::Int) where {T} = size(A.data, idx)
-
-Base.getindex(A::QuantumObject{<:AbstractArray{T}}, inds...) where {T} = getindex(A.data, inds...)
-Base.setindex!(A::QuantumObject{<:AbstractArray{T}}, val, inds...) where {T} = setindex!(A.data, val, inds...)
-
-@doc raw"""
-    eltype(A::QuantumObject)
-
-Returns the elements type of the matrix or vector corresponding to the [`QuantumObject`](@ref) `A`.
-"""
-Base.eltype(A::QuantumObject) = eltype(A.data)
-
-@doc raw"""
-    length(A::QuantumObject)
-
-Returns the length of the matrix or vector corresponding to the [`QuantumObject`](@ref) `A`.
-"""
-Base.length(A::QuantumObject{<:AbstractArray{T}}) where {T} = length(A.data)
-
-Base.isequal(A::QuantumObject{<:AbstractArray{T}}, B::QuantumObject{<:AbstractArray{T}}) where {T} =
-    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isequal(A.data, B.data)
-Base.isapprox(A::QuantumObject{<:AbstractArray{T}}, B::QuantumObject{<:AbstractArray{T}}; kwargs...) where {T} =
-    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isapprox(A.data, B.data; kwargs...)
-Base.:(==)(A::QuantumObject{<:AbstractArray{T}}, B::QuantumObject{<:AbstractArray{T}}) where {T} =
-    (A.type == B.type) && (A.dims == B.dims) && (A.data == B.data)
-
 Base.real(x::QuantumObject) = QuantumObject(real(x.data), x.type, x.dims)
 Base.imag(x::QuantumObject) = QuantumObject(imag(x.data), x.type, x.dims)
 
@@ -368,7 +179,3 @@ SparseArrays.SparseMatrixCSC(A::QuantumObject{<:AbstractMatrix}) =
     QuantumObject(SparseMatrixCSC(A.data), A.type, A.dims)
 SparseArrays.SparseMatrixCSC{T}(A::QuantumObject{<:SparseMatrixCSC}) where {T<:Number} =
     QuantumObject(SparseMatrixCSC{T}(A.data), A.type, A.dims)
-
-# functions for getting Float or Complex element type
-_FType(::QuantumObject{<:AbstractArray{T}}) where {T<:Number} = _FType(T)
-_CType(::QuantumObject{<:AbstractArray{T}}) where {T<:Number} = _CType(T)
diff --git a/src/qobj/quantum_object_base.jl b/src/qobj/quantum_object_base.jl
new file mode 100644
index 00000000..0190e4c4
--- /dev/null
+++ b/src/qobj/quantum_object_base.jl
@@ -0,0 +1,216 @@
+#=  
+This file defines the AbstractQuantumObject structure, all the type structures for AbstractQuantumObject, and fundamental functions in Julia standard library:  
+    - Base: show, length, size, eltype, getindex, setindex!, isequal, :(==), isapprox  
+=#
+
+export AbstractQuantumObject
+export QuantumObjectType,
+    BraQuantumObject,
+    KetQuantumObject,
+    OperatorQuantumObject,
+    OperatorBraQuantumObject,
+    OperatorKetQuantumObject,
+    SuperOperatorQuantumObject
+export Bra, Ket, Operator, OperatorBra, OperatorKet, SuperOperator
+
+@doc raw"""
+    abstract type AbstractQuantumObject{DataType,ObjType,N}
+
+Abstract type for all quantum objects like [`QuantumObject`](@ref) and [`QuantumObjectEvolution`](@ref).
+
+# Example
+```
+julia> sigmax() isa AbstractQuantumObject
+true
+```
+"""
+abstract type AbstractQuantumObject{DataType,ObjType,N} end
+
+abstract type QuantumObjectType end
+
+@doc raw"""
+    BraQuantumObject <: QuantumObjectType
+
+Constructor representing a bra state ``\langle\psi|``.
+"""
+struct BraQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::BraQuantumObject) = print(io, "Bra")
+
+@doc raw"""
+    const Bra = BraQuantumObject()
+
+A constant representing the type of [`BraQuantumObject`](@ref): a bra state ``\langle\psi|``
+"""
+const Bra = BraQuantumObject()
+
+@doc raw"""
+    KetQuantumObject <: QuantumObjectType
+
+Constructor representing a ket state ``|\psi\rangle``.
+"""
+struct KetQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::KetQuantumObject) = print(io, "Ket")
+
+@doc raw"""
+    const Ket = KetQuantumObject()
+
+A constant representing the type of [`KetQuantumObject`](@ref): a ket state ``|\psi\rangle``
+"""
+const Ket = KetQuantumObject()
+
+@doc raw"""
+    OperatorQuantumObject <: QuantumObjectType
+
+Constructor representing an operator ``\hat{O}``.
+"""
+struct OperatorQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::OperatorQuantumObject) = print(io, "Operator")
+
+@doc raw"""
+    const Operator = OperatorQuantumObject()
+
+A constant representing the type of [`OperatorQuantumObject`](@ref): an operator ``\hat{O}``
+"""
+const Operator = OperatorQuantumObject()
+
+@doc raw"""
+    SuperOperatorQuantumObject <: QuantumObjectType
+
+Constructor representing a super-operator ``\hat{\mathcal{O}}`` acting on vectorized density operator matrices.
+"""
+struct SuperOperatorQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::SuperOperatorQuantumObject) = print(io, "SuperOperator")
+
+@doc raw"""
+    const SuperOperator = SuperOperatorQuantumObject()
+
+A constant representing the type of [`SuperOperatorQuantumObject`](@ref): a super-operator ``\hat{\mathcal{O}}`` acting on vectorized density operator matrices
+"""
+const SuperOperator = SuperOperatorQuantumObject()
+
+@doc raw"""
+    OperatorBraQuantumObject <: QuantumObjectType
+
+Constructor representing a bra state in the [`SuperOperator`](@ref) formalism ``\langle\langle\rho|``.
+"""
+struct OperatorBraQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::OperatorBraQuantumObject) = print(io, "OperatorBra")
+
+@doc raw"""
+    const OperatorBra = OperatorBraQuantumObject()
+
+A constant representing the type of [`OperatorBraQuantumObject`](@ref): a bra state in the [`SuperOperator`](@ref) formalism ``\langle\langle\rho|``.
+"""
+const OperatorBra = OperatorBraQuantumObject()
+
+@doc raw"""
+    OperatorKetQuantumObject <: QuantumObjectType
+
+Constructor representing a ket state in the [`SuperOperator`](@ref) formalism ``|\rho\rangle\rangle``.
+"""
+struct OperatorKetQuantumObject <: QuantumObjectType end
+Base.show(io::IO, ::OperatorKetQuantumObject) = print(io, "OperatorKet")
+
+@doc raw"""
+    const OperatorKet = OperatorKetQuantumObject()
+
+A constant representing the type of [`OperatorKetQuantumObject`](@ref): a ket state in the [`SuperOperator`](@ref) formalism ``|\rho\rangle\rangle``
+"""
+const OperatorKet = OperatorKetQuantumObject()
+
+@doc raw"""
+    size(A::AbstractQuantumObject)
+    size(A::AbstractQuantumObject, idx::Int)
+
+Returns a tuple containing each dimensions of the array in the [`AbstractQuantumObject`](@ref).
+
+Optionally, you can specify an index (`idx`) to just get the corresponding dimension of the array.
+"""
+Base.size(A::AbstractQuantumObject) = size(A.data)
+Base.size(A::AbstractQuantumObject, idx::Int) = size(A.data, idx)
+
+Base.getindex(A::AbstractQuantumObject, inds...) = getindex(A.data, inds...)
+Base.setindex!(A::AbstractQuantumObject, val, inds...) = setindex!(A.data, val, inds...)
+
+@doc raw"""
+    eltype(A::AbstractQuantumObject)
+
+Returns the elements type of the matrix or vector corresponding to the [`AbstractQuantumObject`](@ref) `A`.
+"""
+Base.eltype(A::AbstractQuantumObject) = eltype(A.data)
+
+@doc raw"""
+    length(A::AbstractQuantumObject)
+
+Returns the length of the matrix or vector corresponding to the [`AbstractQuantumObject`](@ref) `A`.
+"""
+Base.length(A::AbstractQuantumObject) = length(A.data)
+
+Base.isequal(A::AbstractQuantumObject, B::AbstractQuantumObject) =
+    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isequal(A.data, B.data)
+Base.isapprox(A::AbstractQuantumObject, B::AbstractQuantumObject; kwargs...) =
+    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isapprox(A.data, B.data; kwargs...)
+Base.:(==)(A::AbstractQuantumObject, B::AbstractQuantumObject) =
+    (A.type == B.type) && (A.dims == B.dims) && (A.data == B.data)
+
+function check_dims(A::AbstractQuantumObject, B::AbstractQuantumObject)
+    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects don't have the same Hilbert dimension."))
+    return nothing
+end
+
+function _check_dims(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
+    _non_static_array_warning("dims", dims)
+    return (all(>(0), dims) && length(dims) > 0) ||
+           throw(DomainError(dims, "The argument dims must be of non-zero length and contain only positive integers."))
+end
+_check_dims(dims::Any) = throw(
+    ArgumentError(
+        "The argument dims must be a Tuple or a StaticVector of non-zero length and contain only positive integers.",
+    ),
+)
+
+function _check_QuantumObject(type::KetQuantumObject, dims, m::Int, n::Int)
+    (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Ket"))
+    (prod(dims) != m) && throw(DimensionMismatch("Ket with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+function _check_QuantumObject(type::BraQuantumObject, dims, m::Int, n::Int)
+    (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Bra"))
+    (prod(dims) != n) && throw(DimensionMismatch("Bra with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+function _check_QuantumObject(type::OperatorQuantumObject, dims, m::Int, n::Int)
+    (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with Operator"))
+    (prod(dims) != m) &&
+        throw(DimensionMismatch("Operator with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+function _check_QuantumObject(type::SuperOperatorQuantumObject, dims, m::Int, n::Int)
+    (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with SuperOperator"))
+    (prod(dims) != sqrt(m)) &&
+        throw(DimensionMismatch("SuperOperator with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+function _check_QuantumObject(type::OperatorKetQuantumObject, dims, m::Int, n::Int)
+    (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorKet"))
+    (prod(dims) != sqrt(m)) &&
+        throw(DimensionMismatch("OperatorKet with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+function _check_QuantumObject(type::OperatorBraQuantumObject, dims, m::Int, n::Int)
+    (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorBra"))
+    (prod(dims) != sqrt(n)) &&
+        throw(DimensionMismatch("OperatorBra with dims = $(dims) does not fit the array size = $((m, n))."))
+    return nothing
+end
+
+get_typename_wrapper(A::AbstractQuantumObject) = Base.typename(typeof(A)).wrapper
+
+# functions for getting Float or Complex element type
+_FType(A::AbstractQuantumObject) = _FType(eltype(A))
+_CType(A::AbstractQuantumObject) = _CType(eltype(A))
diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
new file mode 100644
index 00000000..9d25858e
--- /dev/null
+++ b/src/qobj/quantum_object_evo.jl
@@ -0,0 +1,418 @@
+export QuantumObjectEvolution
+
+@doc raw"""
+    struct QuantumObjectEvolution{DT<:AbstractSciMLOperator,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject
+        data::DT
+        type::ObjType
+        dims::SVector{N,Int}
+    end
+
+Julia struct representing any time-dependent quantum object. The `data` field is a `AbstractSciMLOperator` object that represents the time-dependent quantum object. It can be seen as
+
+```math
+\hat{O}(t) = \sum_{i} c_i(p, t) \hat{O}_i
+```
+
+where ``c_i(p, t)`` is a function that depends on the parameters `p` and time `t`, and ``\hat{O}_i`` are the operators that form the quantum object. The `type` field is the type of the quantum object, and the `dims` field is the dimensions of the quantum object. For more information about `type` and `dims`, see [`QuantumObject`](@ref). For more information about `AbstractSciMLOperator`, see the [SciML](https://docs.sciml.ai/SciMLOperators/stable/) documentation.
+
+# Examples
+This operator can be initialized in the same way as the QuTiP `QobjEvo` object. For example
+```
+julia> a = tensor(destroy(10), qeye(2))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 18 stored entries:
+⎡⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦
+
+julia> σm = tensor(qeye(10), sigmam())
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+⎡⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦
+
+julia> coef1(p, t) = exp(-1im * t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = sin(t)
+coef2 (generic function with 1 method)
+
+julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+```
+
+We can also concretize the operator at a specific time `t`
+```
+julia> op1(0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+```
+
+It also supports parameter-dependent time evolution
+```
+julia> coef1(p, t) = exp(-1im * p.ω1 * t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = sin(p.ω2 * t)
+coef2 (generic function with 1 method)
+
+julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+
+julia> p = (ω1 = 1.0, ω2 = 0.5)
+(ω1 = 1.0, ω2 = 0.5)
+
+julia> op1(p, 0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+```
+"""
+struct QuantumObjectEvolution{
+    DT<:AbstractSciMLOperator,
+    ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    N,
+} <: AbstractQuantumObject{DT,ObjType,N}
+    data::DT
+    type::ObjType
+    dims::SVector{N,Int}
+
+    function QuantumObjectEvolution(
+        data::DT,
+        type::ObjType,
+        dims,
+    ) where {DT<:AbstractSciMLOperator,ObjType<:QuantumObjectType}
+        (type == Operator || type == SuperOperator) ||
+            throw(ArgumentError("The type $type is not supported for QuantumObjectEvolution."))
+
+        _check_dims(dims)
+
+        _size = _get_size(data)
+        _check_QuantumObject(type, dims, _size[1], _size[2])
+
+        N = length(dims)
+
+        return new{DT,ObjType,N}(data, type, SVector{N,Int}(dims))
+    end
+end
+
+function QuantumObjectEvolution(data::AbstractSciMLOperator, type::QuantumObjectType, dims::Integer)
+    return QuantumObjectEvolution(data, type, SVector{1,Int}(dims))
+end
+
+"""
+    QuantumObjectEvolution(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
+
+Generate a [`QuantumObjectEvolution`](@ref) object from a [`SciMLOperator`](https://github.com/SciML/SciMLOperators.jl), in the same way as [`QuantumObject`](@ref) for `AbstractArray` inputs.
+"""
+function QuantumObjectEvolution(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
+    _size = _get_size(data)
+
+    if dims isa Nothing
+        if type == Operator
+            dims = SVector{1,Int}(_size[2])
+        elseif type == SuperOperator
+            dims = SVector{1,Int}(isqrt(_size[2]))
+        end
+    end
+
+    return QuantumObjectEvolution(data, type, dims)
+end
+
+# Make the QuantumObjectEvolution, with the option to pre-multiply by a scalar
+function QuantumObjectEvolution(
+    op_func_list::Tuple,
+    α::Union{Nothing,Number} = nothing;
+    type::Union{Nothing,QuantumObjectType} = nothing,
+    f::Function = identity,
+)
+    op, data = _QobjEvo_generate_data(op_func_list, α; f = f)
+    dims = op.dims
+    if type isa Nothing
+        type = op.type
+    end
+
+    # Preallocate the SciMLOperator cache using a dense vector as a reference
+    v0 = sparse_to_dense(similar(op.data, size(op, 1)))
+    data = cache_operator(data, v0)
+
+    return QuantumObjectEvolution(data, type, dims)
+end
+
+function QuantumObjectEvolution(
+    op::QuantumObject,
+    α::Union{Nothing,Number} = nothing;
+    type::Union{Nothing,QuantumObjectType} = nothing,
+    f::Function = identity,
+)
+    if type isa Nothing
+        type = op.type
+    end
+    return QuantumObjectEvolution(_make_SciMLOperator(op, α, f = f), type, op.dims)
+end
+
+function QuantumObjectEvolution(
+    op::QuantumObjectEvolution,
+    α::Union{Nothing,Number} = nothing;
+    type::Union{Nothing,QuantumObjectType} = nothing,
+    f::Function = identity,
+)
+    f !== identity && throw(ArgumentError("The function `f` is not supported for QuantumObjectEvolution inputs."))
+    if type isa Nothing
+        type = op.type
+    else
+        throw(
+            ArgumentError(
+                "The type of the QuantumObjectEvolution object cannot be changed when using another QuantumObjectEvolution object as input.",
+            ),
+        )
+    end
+    if α isa Nothing
+        return QuantumObjectEvolution(op.data, type, op.dims)
+    end
+    return QuantumObjectEvolution(α * op.data, type, op.dims)
+end
+
+#=
+    _QobjEvo_generate_data(op_func_list::Tuple, α; f::Function=identity)
+
+Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution` object. The `op_func_list` is a tuple of tuples or operators. Each element of the tuple can be a tuple with two elements (operator, function) or an operator. The function is used to generate the time-dependent coefficients for the operators. The `α` parameter is used to pre-multiply the operators by a scalar. The `f` parameter is used to pre-applying a function to the operators before converting them to SciML operators. During the parsing, the dimensions of the operators are checked to be the same, and all the constant operators are summed together.
+
+# Arguments
+- `op_func_list::Tuple`: A tuple of tuples or operators.
+- `α`: A scalar to pre-multiply the operators.
+- `f::Function=identity`: A function to pre-apply to the operators.
+=#
+@generated function _QobjEvo_generate_data(op_func_list::Tuple, α; f::Function = identity)
+    op_func_list_types = op_func_list.parameters
+    N = length(op_func_list_types)
+
+    dims_expr = ()
+    first_op = nothing
+    data_expr = :(0)
+    qobj_expr_const = :(0)
+
+    for i in 1:N
+        op_func_type = op_func_list_types[i]
+        if op_func_type <: Tuple
+            length(op_func_type.parameters) == 2 || throw(ArgumentError("The tuple must have two elements."))
+            op_type = op_func_type.parameters[1]
+            func_type = op_func_type.parameters[2]
+            ((isoper(op_type) || issuper(op_type)) && func_type <: Function) || throw(
+                ArgumentError(
+                    "The first element must be a Operator or SuperOperator, and the second element must be a function.",
+                ),
+            )
+
+            op = :(op_func_list[$i][1])
+            data_type = op_type.parameters[1]
+            dims_expr = (dims_expr..., :($op.dims))
+            if i == 1
+                first_op = :(f($op))
+            end
+            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α, f = f))
+        else
+            op_type = op_func_type
+            (isoper(op_type) || issuper(op_type)) ||
+                throw(ArgumentError("The element must be a Operator or SuperOperator."))
+
+            data_type = op_type.parameters[1]
+            dims_expr = (dims_expr..., :(op_func_list[$i].dims))
+            if i == 1
+                first_op = :(op_func_list[$i])
+            end
+            qobj_expr_const = :($qobj_expr_const + f(op_func_list[$i]))
+        end
+    end
+
+    quote
+        dims = tuple($(dims_expr...))
+
+        length(unique(dims)) == 1 || throw(ArgumentError("The dimensions of the operators must be the same."))
+
+        data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const, α)
+
+        data_expr = data_expr_const + $data_expr
+
+        return $first_op, data_expr
+    end
+end
+
+function _make_SciMLOperator(op_func::Tuple, α; f::Function = identity)
+    T = eltype(op_func[1])
+    update_func = (a, u, p, t) -> op_func[2](p, t)
+    if α isa Nothing
+        return ScalarOperator(zero(T), update_func) * MatrixOperator(f(op_func[1]).data)
+    end
+    return ScalarOperator(zero(T), update_func) * MatrixOperator(α * f(op_func[1]).data)
+end
+
+function _make_SciMLOperator(op::QuantumObject, α; f::Function = identity)
+    if α isa Nothing
+        return MatrixOperator(f(op).data)
+    end
+    return MatrixOperator(α * f(op.data))
+end
+
+@doc raw"""
+    (A::QuantumObjectEvolution)(ψout, ψin, p, t)
+
+Apply the time-dependent [`QuantumObjectEvolution`](@ref) object `A` to the input state `ψin` at time `t` with parameters `p`. The output state is stored in `ψout`. This function mimics the behavior of a `AbstractSciMLOperator` object.
+
+# Arguments
+- `ψout::QuantumObject`: The output state. It must have the same type as `ψin`.
+- `ψin::QuantumObject`: The input state. It must be either a [`KetQuantumObject`](@ref) or a [`OperatorKetQuantumObject`](@ref).
+- `p`: The parameters of the time-dependent coefficients.
+- `t`: The time at which the coefficients are evaluated.
+
+# Returns
+- `ψout::QuantumObject`: The output state.
+
+# Examples
+```
+julia> a = destroy(20)
+Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
+⎡⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⎦
+
+julia> coef1(p, t) = sin(t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = cos(t)
+coef2 (generic function with 1 method)
+
+julia> A = QobjEvo(((a, coef1), (a', coef2)))
+Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+
+julia> ψ1 = fock(20, 3)
+Quantum Object:   type=Ket   dims=[20]   size=(20,)
+20-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+     ⋮
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> ψ2 = zero_ket(20)
+Quantum Object:   type=Ket   dims=[20]   size=(20,)
+20-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+     ⋮
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> A(ψ2, ψ1, nothing, 0.1)
+20-element Vector{ComplexF64}:
+                0.0 + 0.0im
+                0.0 + 0.0im
+ 0.1729165499254989 + 0.0im
+                0.0 + 0.0im
+ 1.9900083305560516 + 0.0im
+                    ⋮
+                0.0 + 0.0im
+                0.0 + 0.0im
+                0.0 + 0.0im
+                0.0 + 0.0im
+```
+"""
+function (A::QuantumObjectEvolution)(
+    ψout::QuantumObject{DT1,QobjType},
+    ψin::QuantumObject{DT2,QobjType},
+    p,
+    t,
+) where {DT1,DT2,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    check_dims(ψout, ψin)
+    check_dims(ψout, A)
+
+    if isoper(A) && isoperket(ψin)
+        throw(
+            ArgumentError(
+                "The input state must be a KetQuantumObject if the QuantumObjectEvolution object is an Operator.",
+            ),
+        )
+    elseif issuper(A) && isket(ψin)
+        throw(
+            ArgumentError(
+                "The input state must be an OperatorKetQuantumObject if the QuantumObjectEvolution object is a SuperOperator.",
+            ),
+        )
+    end
+
+    A.data(ψout.data, ψin.data, p, t)
+
+    return ψout
+end
+
+@doc raw"""
+    (A::QuantumObjectEvolution)(ψ, p, t)
+
+Apply the time-dependent [`QuantumObjectEvolution`](@ref) object `A` to the input state `ψ` at time `t` with parameters `p`. Out-of-place version of [`(A::QuantumObjectEvolution)(ψout, ψin, p, t)`](@ref). The output state is stored in a new [`QuantumObject`](@ref) object. This function mimics the behavior of a `AbstractSciMLOperator` object.
+"""
+function (A::QuantumObjectEvolution)(
+    ψ::QuantumObject{DT,QobjType},
+    p,
+    t,
+) where {DT,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    ψout = QuantumObject(similar(ψ.data), ψ.type, ψ.dims)
+    return A(ψout, ψ, p, t)
+end
+
+@doc raw"""
+    (A::QuantumObjectEvolution)(p, t)
+
+Calculate the time-dependent [`QuantumObjectEvolution`](@ref) object `A` at time `t` with parameters `p`.
+
+# Arguments
+- `p`: The parameters of the time-dependent coefficients.
+- `t`: The time at which the coefficients are evaluated.
+
+# Returns
+- `A::QuantumObject`: The output state.
+"""
+function (A::QuantumObjectEvolution)(p, t)
+    # We put 0 in the place of `u` because the time-dependence doesn't depend on the state
+    update_coefficients!(A.data, 0, p, t)
+    return QuantumObject(concretize(A.data), A.type, A.dims)
+end
+
+(A::QuantumObjectEvolution)(t) = A(nothing, t)
+
+#=
+`promote_type` should be applied on types. Here I define `promote_op_type` because it is applied to operators.
+=#
+promote_op_type(A::QuantumObjectEvolution, B::QuantumObjectEvolution) = get_typename_wrapper(A)
+promote_op_type(A::QuantumObjectEvolution, B::QuantumObject) = get_typename_wrapper(A)
+promote_op_type(A::QuantumObject, B::QuantumObjectEvolution) = get_typename_wrapper(B)
+promote_op_type(A::QuantumObject, B::QuantumObject) = get_typename_wrapper(A)
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index b95c2ee2..08f3672a 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -69,7 +69,7 @@ function sprepost(
     A::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
 ) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    check_dims(A, B)
 
     return QuantumObject(_sprepost(A.data, B.data), SuperOperator, A.dims)
 end
@@ -89,13 +89,48 @@ the same function is applied multiple times with a known Hilbert space dimension
 
 See also [`spre`](@ref) and [`spost`](@ref).
 """
-function lindblad_dissipator(
-    O::QuantumObject{<:AbstractArray{T},OperatorQuantumObject},
-    Id_cache = I(size(O, 1)),
-) where {T}
+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
 end
 
 # It is already a SuperOperator
-lindblad_dissipator(O::QuantumObject{<:AbstractArray{T},SuperOperatorQuantumObject}, Id_cache) where {T} = O
+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)))
+
+Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
+
+```math
+\mathcal{L} [\cdot] = -i[\hat{H}, \cdot] + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
+```
+
+where 
+
+```math
+\mathcal{D}(\hat{C}_n) [\cdot] = \hat{C}_n [\cdot] \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n [\cdot] - \frac{1}{2} [\cdot] \hat{C}_n^\dagger \hat{C}_n
+```
+
+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), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
+"""
+function liouvillian(
+    H::QuantumObject{MT1,OpType1},
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+    Id_cache = I(prod(H.dims)),
+) where {MT1<:AbstractMatrix,OpType1<: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
+    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::QuantumObject{<:AbstractMatrix,SuperOperatorQuantumObject}, Id_cache::Diagonal) = H
diff --git a/src/qobj/synonyms.jl b/src/qobj/synonyms.jl
index 17f59e99..8155ba3f 100644
--- a/src/qobj/synonyms.jl
+++ b/src/qobj/synonyms.jl
@@ -2,7 +2,7 @@
 Synonyms of the functions for QuantumObject
 =#
 
-export Qobj, shape, isherm
+export Qobj, QobjEvo, shape, isherm
 export trans, dag, matrix_element, unit
 export sqrtm, logm, expm, sinm, cosm
 export tensor, ⊗
@@ -13,47 +13,148 @@ export qeye
 
 Generate [`QuantumObject`](@ref)
 
-Note that this functions is same as `QuantumObject(A; type=type, dims=dims)`
+Note that this functions is same as `QuantumObject(A; type=type, dims=dims)`.
 """
 Qobj(A; kwargs...) = QuantumObject(A; kwargs...)
 
 @doc raw"""
-    shape(A::QuantumObject)
+    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type::Union{Nothing, QuantumObjectType}=nothing, f::Function=identity)
 
-Returns a tuple containing each dimensions of the array in the [`QuantumObject`](@ref).
+Generate [`QuantumObjectEvolution`](@ref).
 
-Note that this function is same as `size(A)`
+Note that this functions is same as `QuantumObjectEvolution(op_func_list)`. If `α` is provided, all the operators in `op_func_list` will be pre-multiplied by `α`. The `type` parameter is used to specify the type of the [`QuantumObject`](@ref), either `Operator` or `SuperOperator`. The `f` parameter is used to pre-apply a function to the operators before converting them to SciML operators.
+
+# Arguments
+- `op_func_list::Union{Tuple,AbstractQuantumObject}`: A tuple of tuples or operators.
+- `α::Union{Nothing,Number}=nothing`: A scalar to pre-multiply the operators.
+- `f::Function=identity`: A function to pre-apply to the operators.
+
+!!! warning "Beware of type-stability!"
+    Please note that, unlike QuTiP, this function doesn't support `op_func_list` as `Vector` type. This is related to the type-stability issue. See the Section [The Importance of Type-Stability](@ref doc:Type-Stability) for more details.
+
+# Examples
+This operator can be initialized in the same way as the QuTiP `QobjEvo` object. For example
+```
+julia> a = tensor(destroy(10), qeye(2))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 18 stored entries:
+⎡⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦
+
+julia> σm = tensor(qeye(10), sigmam())
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+⎡⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦
+
+julia> coef1(p, t) = exp(-1im * t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = sin(t)
+coef2 (generic function with 1 method)
+
+julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+```
+
+We can also concretize the operator at a specific time `t`
+```
+julia> op1(0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+```
+
+It also supports parameter-dependent time evolution
+```
+julia> coef1(p, t) = exp(-1im * p.ω1 * t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = sin(p.ω2 * t)
+coef2 (generic function with 1 method)
+
+julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+
+julia> p = (ω1 = 1.0, ω2 = 0.5)
+(ω1 = 1.0, ω2 = 0.5)
+
+julia> op1(p, 0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+```
+"""
+QobjEvo(
+    op_func_list::Union{Tuple,AbstractQuantumObject},
+    α::Union{Nothing,Number} = nothing;
+    type::Union{Nothing,QuantumObjectType} = nothing,
+    f::Function = identity,
+) = QuantumObjectEvolution(op_func_list, α; type = type, f = f)
+
+QobjEvo(data::AbstractSciMLOperator, type::QuantumObjectType, dims::Integer) =
+    QuantumObjectEvolution(data, type, SVector{1,Int}(dims))
+
+"""
+    QobjEvo(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
+
+Synonym of [`QuantumObjectEvolution`](@ref) object from a [`SciMLOperator`](https://github.com/SciML/SciMLOperators.jl).  See the documentation for [`QuantumObjectEvolution`](@ref) for more information.
+"""
+QobjEvo(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing) =
+    QuantumObjectEvolution(data; type = type, dims = dims)
+
+@doc raw"""
+    shape(A::AbstractQuantumObject)
+
+Returns a tuple containing each dimensions of the array in the [`AbstractQuantumObject`](@ref).
+
+Note that this function is same as `size(A)`.
 """
-shape(A::QuantumObject{<:AbstractArray{T}}) where {T} = size(A.data)
+shape(A::AbstractQuantumObject) = size(A.data)
 
 @doc raw"""
-    isherm(A::QuantumObject)
+    isherm(A::AbstractQuantumObject)
 
-Test whether the [`QuantumObject`](@ref) is Hermitian.
+Test whether the [`AbstractQuantumObject`](@ref) is Hermitian.
 
-Note that this functions is same as `ishermitian(A)`
+Note that this functions is same as `ishermitian(A)`.
 """
-isherm(A::QuantumObject{<:AbstractArray{T}}) where {T} = ishermitian(A)
+isherm(A::AbstractQuantumObject) = ishermitian(A)
 
 @doc raw"""
-    trans(A::QuantumObject)
+    trans(A::AbstractQuantumObject)
 
-Lazy matrix transpose of the [`QuantumObject`](@ref).
+Lazy matrix transpose of the [`AbstractQuantumObject`](@ref).
 
-Note that this function is same as `transpose(A)`
+Note that this function is same as `transpose(A)`.
 """
-trans(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = transpose(A)
+trans(A::AbstractQuantumObject{DT,OpType}) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    transpose(A)
 
 @doc raw"""
-    dag(A::QuantumObject)
+    dag(A::AbstractQuantumObject)
 
-Lazy adjoint (conjugate transposition) of the [`QuantumObject`](@ref)
+Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 
-Note that this function is same as `adjoint(A)`
+Note that this function is same as `adjoint(A)`.
 """
-dag(A::QuantumObject{<:AbstractArray{T}}) where {T} = adjoint(A)
+dag(A::AbstractQuantumObject) = adjoint(A)
 
 @doc raw"""
     matrix_element(i::QuantumObject, A::QuantumObject j::QuantumObject)
@@ -110,7 +211,7 @@ sqrtm(A::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}) where {T} = sq
 
 Matrix logarithm of [`QuantumObject`](@ref)
 
-Note that this function is same as `log(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
+Note that this function is same as `log(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref).
 """
 logm(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
@@ -121,7 +222,7 @@ logm(
 
 Matrix exponential of [`QuantumObject`](@ref)
 
-Note that this function is same as `exp(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
+Note that this function is same as `exp(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref).
 """
 expm(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
@@ -134,7 +235,7 @@ Matrix sine of [`QuantumObject`](@ref), defined as
 
 ``\sin \left( \hat{A} \right) = \frac{e^{i \hat{A}} - e^{-i \hat{A}}}{2 i}``
 
-Note that this function is same as `sin(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
+Note that this function is same as `sin(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref).
 """
 sinm(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
@@ -147,7 +248,7 @@ Matrix cosine of [`QuantumObject`](@ref), defined as
 
 ``\cos \left( \hat{A} \right) = \frac{e^{i \hat{A}} + e^{-i \hat{A}}}{2}``
 
-Note that this function is same as `cos(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
+Note that this function is same as `cos(A)` and only supports for [`Operator`](@ref) and [`SuperOperator`](@ref).
 """
 cosm(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
@@ -158,7 +259,7 @@ cosm(
 
 Returns the [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product) ``\hat{A} \otimes \hat{B} \otimes \cdots``.
 
-Note that this function is same as `kron(A, B, ...)`
+Note that this function is same as `kron(A, B, ...)`.
 
 # Examples
 
@@ -191,7 +292,7 @@ tensor(A...) = kron(A...)
 
 Returns the [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product) ``\hat{A} \otimes \hat{B}``.
 
-Note that this function is same as `kron(A, B)`
+Note that this function is same as `kron(A, B)`.
 
 # Examples
 
@@ -240,7 +341,7 @@ Identity operator ``\hat{\mathbb{1}}`` with size `N`.
 
 It is also possible to specify the list of Hilbert dimensions `dims` if different subsystems are present.
 
-Note that this function is same as `eye(N, type=type, dims=dims)`, and `type` can only be either [`Operator`](@ref) or [`SuperOperator`](@ref)
+Note that this function is same as `eye(N, type=type, dims=dims)`, and `type` can only be either [`Operator`](@ref) or [`SuperOperator`](@ref).
 """
 qeye(
     N::Int;
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 59e6f6ad..097611d9 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -65,26 +65,26 @@ end
 
 @doc raw"""
     steadystate(
-        H::QuantumObject,
+        H::QuantumObject{DT,OpType},
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateSolver = SteadyStateDirectSolver(),
-        kwargs...
+        kwargs...,
     )
 
 Solve the stationary state based on different solvers.
 
 # Parameters
-- `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators.
-- `solver::SteadyStateSolver=SteadyStateDirectSolver()`: see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for different solvers.
-- `kwargs...`: The keyword arguments for the solver.
+- `H`: The Hamiltonian or the Liouvillian of the system.
+- `c_ops`: The list of the collapse operators.
+- `solver`: see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for different solvers.
+- `kwargs`: The keyword arguments for the solver.
 """
 function steadystate(
-    H::QuantumObject{<:AbstractArray,OpType},
+    H::QuantumObject{DT,OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
-) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, c_ops)
 
     return _steadystate(L, solver; kwargs...)
@@ -179,20 +179,20 @@ _steadystate(
     kwargs...,
 ) where {T} = throw(
     ArgumentError(
-        "The initial state ψ0 is required for SteadyStateODESolver, use the following call instead: `steadystate(H, ψ0, tspan, c_ops)`.",
+        "The initial state ψ0 is required for SteadyStateODESolver, use the following call instead: `steadystate(H, ψ0, tmax, c_ops)`.",
     ),
 )
 
 @doc raw"""
     steadystate(
-        H::QuantumObject,
-        ψ0::QuantumObject,
-        tspan::Real = Inf,
+        H::QuantumObject{DT1,HOpType},
+        ψ0::QuantumObject{DT2,StateOpType},
+        tmax::Real = Inf,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateODESolver = SteadyStateODESolver(),
         reltol::Real = 1.0e-8,
         abstol::Real = 1.0e-10,
-        kwargs...
+        kwargs...,
     )
 
 Solve the stationary state based on time evolution (ordinary differential equations; `OrdinaryDiffEq.jl`) with a given initial state.
@@ -210,50 +210,45 @@ or
 ```
 
 # Parameters
-- `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
-- `ψ0::QuantumObject`: The initial state of the system.
-- `tspan::Real=Inf`: The final time step for the steady state problem.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators.
-- `solver::SteadyStateODESolver=SteadyStateODESolver()`: see [`SteadyStateODESolver`](@ref) for more details.
-- `reltol::Real=1.0e-8`: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
-- `abstol::Real=1.0e-10`: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
-- `kwargs...`: The keyword arguments for the ODEProblem.
+- `H`: The Hamiltonian or the Liouvillian of the system.
+- `ψ0`: The initial state of the system.
+- `tmax=Inf`: The final time step for the steady state problem.
+- `c_ops=nothing`: The list of the collapse operators.
+- `solver`: see [`SteadyStateODESolver`](@ref) for more details.
+- `reltol=1.0e-8`: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
+- `abstol=1.0e-10`: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
+- `kwargs`: The keyword arguments for the ODEProblem.
 """
 function steadystate(
-    H::QuantumObject{MT1,HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
-    tspan::Real = Inf,
+    H::QuantumObject{DT1,HOpType},
+    ψ0::QuantumObject{DT2,StateOpType},
+    tmax::Real = Inf,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateODESolver = SteadyStateODESolver(),
     reltol::Real = 1.0e-8,
     abstol::Real = 1.0e-10,
     kwargs...,
 ) where {
-    MT1<:AbstractMatrix,
-    T2,
+    DT1,
+    DT2,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    (H.dims != ψ0.dims) && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-
-    N = prod(H.dims)
-    u0 = sparse_to_dense(_CType(ψ0), mat2vec(ket2dm(ψ0).data))
-
-    L = MatrixOperator(liouvillian(H, c_ops).data)
-
     ftype = _FType(ψ0)
-    prob = ODEProblem{true}(L, u0, (ftype(0), ftype(tspan))) # Convert tspan to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
-    sol = solve(
-        prob,
-        solver.alg;
-        callback = TerminateSteadyState(abstol, reltol, _steadystate_ode_condition),
-        reltol = reltol,
-        abstol = abstol,
-        kwargs...,
+    cb = TerminateSteadyState(abstol, reltol, _steadystate_ode_condition)
+    sol = mesolve(
+        H,
+        ψ0,
+        [ftype(0), ftype(tmax)],
+        c_ops,
+        progress_bar = Val(false),
+        save_everystep = false,
+        saveat = ftype[],
+        callback = cb,
     )
 
-    ρss = reshape(sol.u[end], N, N)
-    return QuantumObject(ρss, Operator, H.dims)
+    ρss = sol.states[end]
+    return ρss
 end
 
 function _steadystate_ode_condition(integrator, abstol, reltol, min_t)
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index fb08a28d..31d5b1ec 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -81,7 +81,9 @@ function _mcsolve_prob_func(prob, i, repeat)
         ),
     )
 
-    return remake(prob, p = prm)
+    f = deepcopy(prob.f.f)
+
+    return remake(prob, f = f, p = prm)
 end
 
 # Standard output function
@@ -117,20 +119,33 @@ function _mcsolve_generate_statistics(sol, i, states, expvals_all, jump_times, j
     return jump_which[i] = sol_i.prob.p.jump_which
 end
 
+function _mcsolve_make_Heff_QobjEvo(H::QuantumObject, c_ops)
+    c_ops isa Nothing && return QobjEvo(H)
+    return QobjEvo(H - 1im * mapreduce(op -> op' * op, +, c_ops) / 2)
+end
+function _mcsolve_make_Heff_QobjEvo(H::Tuple, c_ops)
+    c_ops isa Nothing && return QobjEvo(H)
+    return QobjEvo((H..., -1im * mapreduce(op -> op' * op, +, c_ops) / 2))
+end
+function _mcsolve_make_Heff_QobjEvo(H::QuantumObjectEvolution, c_ops)
+    c_ops isa Nothing && return H
+    return H + QobjEvo(mapreduce(op -> op' * op, +, c_ops), -1im / 2)
+end
+
 @doc raw"""
-    mcsolveProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-        ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    mcsolveProblem(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector,
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        rng::AbstractRNG=default_rng(),
-        jump_callback::TJC=ContinuousLindbladJumpCallback(),
-        kwargs...)
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        params::NamedTuple = NamedTuple(),
+        rng::AbstractRNG = default_rng(),
+        jump_callback::TJC = ContinuousLindbladJumpCallback(),
+        kwargs...,
+    )
 
-Generates the ODEProblem for a single trajectory of the Monte Carlo wave function time evolution of an open quantum system.
+Generate the ODEProblem for a single trajectory of the Monte Carlo wave function time evolution of an open quantum system.
 
 Given a system Hamiltonian ``\hat{H}`` and a list of collapse (jump) operators ``\{\hat{C}_n\}_n``, the evolution of the state ``|\psi(t)\rangle`` is governed by the Schrodinger equation:
 
@@ -166,24 +181,21 @@ If the environmental measurements register a quantum jump, the wave function und
 
 # Arguments
 
-- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
-- `jump_callback::LindbladJumpCallbackType`: The Jump Callback type: Discrete or Continuous.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -191,41 +203,37 @@ If the environmental measurements register a quantum jump, the wave function und
 - `prob::ODEProblem`: The ODEProblem for the Monte Carlo wave function time evolution.
 """
 function mcsolveProblem(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray,KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
     jump_callback::TJC = ContinuousLindbladJumpCallback(),
     kwargs...,
-) where {MT1<:AbstractMatrix,TJC<:LindbladJumpCallbackType}
-    H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-
+) where {DT1,DT2,TJC<:LindbladJumpCallbackType}
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
     c_ops isa Nothing &&
         throw(ArgumentError("The list of collapse operators must be provided. Use sesolveProblem instead."))
 
-    t_l = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
+    tlist = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
 
-    H_eff = H - 1im * mapreduce(op -> op' * op, +, c_ops) / 2
+    H_eff_evo = _mcsolve_make_Heff_QobjEvo(H, c_ops)
 
     if e_ops isa Nothing
-        expvals = Array{ComplexF64}(undef, 0, length(t_l))
+        expvals = Array{ComplexF64}(undef, 0, length(tlist))
         is_empty_e_ops_mc = true
-        e_ops2 = MT1[]
+        e_ops_data = ()
     else
-        expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
-        e_ops2 = get_data.(e_ops)
+        expvals = Array{ComplexF64}(undef, length(e_ops), length(tlist))
+        e_ops_data = get_data.(e_ops)
         is_empty_e_ops_mc = isempty(e_ops)
     end
 
-    saveat = is_empty_e_ops_mc ? t_l : [t_l[end]]
+    saveat = is_empty_e_ops_mc ? tlist : [tlist[end]]
     # We disable the progress bar of the sesolveProblem because we use a global progress bar for all the trajectories
     default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., saveat = saveat, progress_bar = Val(false))
     kwargs2 = merge(default_values, kwargs)
@@ -243,9 +251,9 @@ function mcsolveProblem(
 
     params2 = (
         expvals = expvals,
-        e_ops_mc = e_ops2,
+        e_ops_mc = e_ops_data,
         is_empty_e_ops_mc = is_empty_e_ops_mc,
-        progr_mc = ProgressBar(length(t_l), enable = false),
+        progr_mc = ProgressBar(length(tlist), enable = false),
         traj_rng = rng,
         c_ops = c_ops_data,
         c_ops_herm = c_ops_herm_data,
@@ -259,39 +267,35 @@ function mcsolveProblem(
         params...,
     )
 
-    return mcsolveProblem(H_eff, ψ0, t_l, alg, H_t, params2, jump_callback; kwargs2...)
+    return mcsolveProblem(H_eff_evo, ψ0, tlist, params2, jump_callback; kwargs2...)
 end
 
 function mcsolveProblem(
-    H_eff::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-    t_l::AbstractVector,
-    alg::OrdinaryDiffEqAlgorithm,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum},
+    H_eff_evo::QuantumObjectEvolution{DT1,OperatorQuantumObject},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
+    tlist::AbstractVector,
     params::NamedTuple,
     jump_callback::DiscreteLindbladJumpCallback;
     kwargs...,
-) where {T1,T2}
+) where {DT1,DT2}
     cb1 = DiscreteCallback(LindbladJumpDiscreteCondition, LindbladJumpAffect!, save_positions = (false, false))
-    cb2 = PresetTimeCallback(t_l, _save_func_mcsolve, save_positions = (false, false))
+    cb2 = PresetTimeCallback(tlist, _save_func_mcsolve, save_positions = (false, false))
     kwargs2 = (; kwargs...)
     kwargs2 =
         haskey(kwargs2, :callback) ? merge(kwargs2, (callback = CallbackSet(cb1, cb2, kwargs2.callback),)) :
         merge(kwargs2, (callback = CallbackSet(cb1, cb2),))
 
-    return sesolveProblem(H_eff, ψ0, t_l; alg = alg, H_t = H_t, params = params, kwargs2...)
+    return sesolveProblem(H_eff_evo, ψ0, tlist; params = params, kwargs2...)
 end
 
 function mcsolveProblem(
-    H_eff::QuantumObject{<:AbstractArray,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray,KetQuantumObject},
-    t_l::AbstractVector,
-    alg::OrdinaryDiffEqAlgorithm,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum},
+    H_eff_evo::QuantumObjectEvolution{DT1,OperatorQuantumObject},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
+    tlist::AbstractVector,
     params::NamedTuple,
     jump_callback::ContinuousLindbladJumpCallback;
     kwargs...,
-)
+) where {DT1,DT2}
     cb1 = ContinuousCallback(
         LindbladJumpContinuousCondition,
         LindbladJumpAffect!,
@@ -299,33 +303,34 @@ function mcsolveProblem(
         interp_points = jump_callback.interp_points,
         save_positions = (false, false),
     )
-    cb2 = PresetTimeCallback(t_l, _save_func_mcsolve, save_positions = (false, false))
+    cb2 = PresetTimeCallback(tlist, _save_func_mcsolve, save_positions = (false, false))
     kwargs2 = (; kwargs...)
     kwargs2 =
         haskey(kwargs2, :callback) ? merge(kwargs2, (callback = CallbackSet(cb1, cb2, kwargs2.callback),)) :
         merge(kwargs2, (callback = CallbackSet(cb1, cb2),))
 
-    return sesolveProblem(H_eff, ψ0, t_l; alg = alg, H_t = H_t, params = params, kwargs2...)
+    return sesolveProblem(H_eff_evo, ψ0, tlist; params = params, kwargs2...)
 end
 
 @doc raw"""
-    mcsolveEnsembleProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-        ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    mcsolveEnsembleProblem(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector,
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        ntraj::Int=1,
-        ensemble_method=EnsembleThreads(),
-        jump_callback::TJC=ContinuousLindbladJumpCallback(),
-        prob_func::Function=_mcsolve_prob_func,
-        output_func::Function=_mcsolve_output_func,
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
-
-Generates the `EnsembleProblem` of `ODEProblem`s for the ensemble of trajectories of the Monte Carlo wave function time evolution of an open quantum system.
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        params::NamedTuple = NamedTuple(),
+        rng::AbstractRNG = default_rng(),
+        ntraj::Int = 1,
+        ensemble_method = EnsembleThreads(),
+        jump_callback::TJC = ContinuousLindbladJumpCallback(),
+        prob_func::Function = _mcsolve_prob_func,
+        output_func::Function = _mcsolve_dispatch_output_func(ensemble_method),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
+
+Generate the `EnsembleProblem` of `ODEProblem`s for the ensemble of trajectories of the Monte Carlo wave function time evolution of an open quantum system.
 
 Given a system Hamiltonian ``\hat{H}`` and a list of collapse (jump) operators ``\{\hat{C}_n\}_n``, the evolution of the state ``|\psi(t)\rangle`` is governed by the Schrodinger equation:
 
@@ -361,29 +366,26 @@ If the environmental measurements register a quantum jump, the wave function und
 
 # Arguments
 
-- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
-- `ntraj::Int`: Number of trajectories to use.
-- `ensemble_method`: Ensemble method to use.
-- `jump_callback::LindbladJumpCallbackType`: The Jump Callback type: Discrete or Continuous.
-- `prob_func::Function`: Function to use for generating the ODEProblem.
-- `output_func::Function`: Function to use for generating the output of a single trajectory.
-- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `ntraj`: Number of trajectories to use.
+- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`.
+- `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
+- `prob_func`: Function to use for generating the ODEProblem.
+- `output_func`: Function to use for generating the output of a single trajectory.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -391,13 +393,11 @@ If the environmental measurements register a quantum jump, the wave function und
 - `prob::EnsembleProblem with ODEProblem`: The Ensemble ODEProblem for the Monte Carlo wave function time evolution.
 """
 function mcsolveEnsembleProblem(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
     ntraj::Int = 1,
@@ -407,7 +407,7 @@ function mcsolveEnsembleProblem(
     output_func::Function = _mcsolve_dispatch_output_func(ensemble_method),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2,TJC<:LindbladJumpCallbackType}
+) where {DT1,DT2,TJC<:LindbladJumpCallbackType}
     progr = ProgressBar(ntraj, enable = getVal(progress_bar))
     if ensemble_method isa EnsembleDistributed
         progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
@@ -427,9 +427,7 @@ function mcsolveEnsembleProblem(
             ψ0,
             tlist,
             c_ops;
-            alg = alg,
             e_ops = e_ops,
-            H_t = H_t,
             params = merge(params, (global_rng = rng, seeds = seeds)),
             rng = rng,
             jump_callback = jump_callback,
@@ -448,13 +446,13 @@ function mcsolveEnsembleProblem(
 end
 
 @doc raw"""
-    mcsolve(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-        ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    mcsolve(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         alg::OrdinaryDiffEqAlgorithm = Tsit5(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
         params::NamedTuple = NamedTuple(),
         rng::AbstractRNG = default_rng(),
         ntraj::Int = 1,
@@ -502,22 +500,21 @@ If the environmental measurements register a quantum jump, the wave function und
 
 # Arguments
 
-- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
-- `ntraj::Int`: Number of trajectories to use.
-- `ensemble_method`: Ensemble method to use.
-- `jump_callback::LindbladJumpCallbackType`: The Jump Callback type: Discrete or Continuous.
-- `prob_func::Function`: Function to use for generating the ODEProblem.
-- `output_func::Function`: Function to use for generating the output of a single trajectory.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
-- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `alg`: The algorithm to use for the ODE solver. Default to `Tsit5()`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `ntraj`: Number of trajectories to use.
+- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`.
+- `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
+- `prob_func`: Function to use for generating the ODEProblem.
+- `output_func`: Function to use for generating the output of a single trajectory.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
@@ -530,16 +527,15 @@ If the environmental measurements register a quantum jump, the wave function und
 
 # Returns
 
-- `sol::TimeEvolutionMCSol`: The solution of the time evolution. See also [`TimeEvolutionMCSol`](@ref)
+- `sol::TimeEvolutionMCSol`: The solution of the time evolution. See also [`TimeEvolutionMCSol`](@ref).
 """
 function mcsolve(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
     ntraj::Int = 1,
@@ -549,7 +545,7 @@ function mcsolve(
     output_func::Function = _mcsolve_dispatch_output_func(ensemble_method),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2,TJC<:LindbladJumpCallbackType}
+) where {DT1,DT2,TJC<:LindbladJumpCallbackType}
     ens_prob_mc = mcsolveEnsembleProblem(
         H,
         ψ0,
@@ -557,7 +553,6 @@ function mcsolve(
         c_ops;
         alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         rng = rng,
         ntraj = ntraj,
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 994bdb08..21b8170a 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -19,17 +19,10 @@ function _save_func_mesolve(integrator)
     return u_modified!(integrator, false)
 end
 
-mesolve_ti_dudt!(du, u, p, t) = mul!(du, p.L, u)
-function mesolve_td_dudt!(du, u, p, t)
-    mul!(du, p.L, u)
-    L_t = p.H_t(t, p)
-    return mul!(du, L_t, u, 1, 1)
-end
-
 _generate_mesolve_e_op(op) = mat2vec(adjoint(get_data(op)))
 
-function _generate_mesolve_kwargs_with_callback(t_l, kwargs)
-    cb1 = PresetTimeCallback(t_l, _save_func_mesolve, save_positions = (false, false))
+function _generate_mesolve_kwargs_with_callback(tlist, kwargs)
+    cb1 = PresetTimeCallback(tlist, _save_func_mesolve, save_positions = (false, false))
     kwargs2 =
         haskey(kwargs, :callback) ? merge(kwargs, (callback = CallbackSet(kwargs.callback, cb1),)) :
         merge(kwargs, (callback = cb1,))
@@ -37,30 +30,42 @@ function _generate_mesolve_kwargs_with_callback(t_l, kwargs)
     return kwargs2
 end
 
-function _generate_mesolve_kwargs(e_ops, progress_bar::Val{true}, t_l, kwargs)
-    return _generate_mesolve_kwargs_with_callback(t_l, kwargs)
+function _generate_mesolve_kwargs(e_ops, progress_bar::Val{true}, tlist, kwargs)
+    return _generate_mesolve_kwargs_with_callback(tlist, kwargs)
 end
 
-function _generate_mesolve_kwargs(e_ops, progress_bar::Val{false}, t_l, kwargs)
+function _generate_mesolve_kwargs(e_ops, progress_bar::Val{false}, tlist, kwargs)
     if e_ops isa Nothing
         return kwargs
     end
-    return _generate_mesolve_kwargs_with_callback(t_l, kwargs)
+    return _generate_mesolve_kwargs_with_callback(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)
+    return QobjEvo((H..., mapreduce(op -> lindblad_dissipator(op), +, c_ops)); type = SuperOperator, f = liouvillian)
+end
+# TODO: Add support for Operator type QobEvo
+function _mesolve_make_L_QobjEvo(H::QuantumObjectEvolution, c_ops)
+    issuper(H) || throw(ArgumentError("The time-dependent Hamiltonian must be a SuperOperator."))
+    c_ops isa Nothing && return H
+    return H + QobjEvo((mapreduce(op -> lindblad_dissipator(op), +, c_ops)))
 end
 
 @doc raw"""
-    mesolveProblem(H::QuantumObject,
-        ψ0::QuantumObject,
-        tlist::AbstractVector, 
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
-
-Generates the ODEProblem for the master equation time evolution of an open quantum system:
+    mesolveProblem(
+        H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
+        ψ0::QuantumObject{DT2,StateOpType},
+        tlist,
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        params::NamedTuple = NamedTuple(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
+
+Generate the ODEProblem for the master equation time evolution of an open quantum system:
 
 ```math
 \frac{\partial \hat{\rho}(t)}{\partial t} = -i[\hat{H}, \hat{\rho}(t)] + \sum_n \mathcal{D}(\hat{C}_n) [\hat{\rho}(t)]
@@ -74,23 +79,20 @@ where
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian ``\hat{H}`` or the Liouvillian of the system.
-- `ψ0::QuantumObject`: The initial state of the system.
-- `tlist::AbstractVector`: The time list of the evolution.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::OrdinaryDiffEqAlgorithm=Tsit5()`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the operators for which the expectation values are calculated.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing`: The time-dependent Hamiltonian or Liouvillian.
-- `params::NamedTuple=NamedTuple()`: The parameters of the time evolution.
-- `progress_bar::Union{Val,Bool}=Val(true)`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
-- `kwargs...`: The keyword arguments for the ODEProblem.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -98,84 +100,74 @@ where
 - `prob::ODEProblem`: The ODEProblem for the master equation time evolution.
 """
 function mesolveProblem(
-    H::QuantumObject{MT1,HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
+    ψ0::QuantumObject{DT2,StateOpType},
     tlist,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
 ) where {
-    MT1<:AbstractMatrix,
-    T2,
+    DT1,
+    DT2,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
-    is_time_dependent = !(H_t isa Nothing)
+    tlist = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
 
-    ρ0 = sparse_to_dense(_CType(ψ0), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
+    L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
+    check_dims(L_evo, ψ0)
 
-    t_l = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
+    ρ0 = sparse_to_dense(_CType(ψ0), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
+    L = L_evo.data
 
-    L = liouvillian(H, c_ops).data
-    progr = ProgressBar(length(t_l), enable = getVal(progress_bar))
+    progr = ProgressBar(length(tlist), enable = getVal(progress_bar))
 
     if e_ops isa Nothing
-        expvals = Array{ComplexF64}(undef, 0, length(t_l))
-        e_ops2 = mat2vec(MT1)[]
+        expvals = Array{ComplexF64}(undef, 0, length(tlist))
+        e_ops_data = ()
         is_empty_e_ops = true
     else
-        expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
-        e_ops2 = [_generate_mesolve_e_op(op) for op in e_ops]
+        expvals = Array{ComplexF64}(undef, length(e_ops), length(tlist))
+        e_ops_data = [_generate_mesolve_e_op(op) for op in e_ops]
         is_empty_e_ops = isempty(e_ops)
     end
 
-    if H_t isa TimeDependentOperatorSum
-        H_t = liouvillian(H_t)
-    end
-
     p = (
-        L = L,
-        progr = progr,
-        Hdims = H.dims,
-        e_ops = e_ops2,
+        e_ops = e_ops_data,
         expvals = expvals,
-        H_t = H_t,
-        times = t_l,
+        progr = progr,
+        times = tlist,
+        Hdims = L_evo.dims,
         is_empty_e_ops = is_empty_e_ops,
         params...,
     )
 
-    saveat = is_empty_e_ops ? t_l : [t_l[end]]
+    saveat = is_empty_e_ops ? tlist : [tlist[end]]
     default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., saveat = saveat)
     kwargs2 = merge(default_values, kwargs)
-    kwargs3 = _generate_mesolve_kwargs(e_ops, makeVal(progress_bar), t_l, kwargs2)
-
-    dudt! = is_time_dependent ? mesolve_td_dudt! : mesolve_ti_dudt!
+    kwargs3 = _generate_mesolve_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2)
 
-    tspan = (t_l[1], t_l[end])
-    return ODEProblem{true,FullSpecialize}(dudt!, ρ0, tspan, p; kwargs3...)
+    tspan = (tlist[1], tlist[end])
+    return ODEProblem{true,FullSpecialize}(L, ρ0, tspan, p; kwargs3...)
 end
 
 @doc raw"""
-    mesolve(H::QuantumObject,
-        ψ0::QuantumObject,
-        tlist::AbstractVector, 
-        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
+    mesolve(
+        H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
+        ψ0::QuantumObject{DT2,StateOpType},
+        tlist::AbstractVector,
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        params::NamedTuple = NamedTuple(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
 
 Time evolution of an open quantum system using Lindblad master equation:
 
@@ -191,16 +183,15 @@ where
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian ``\hat{H}`` or the Liouvillian of the system.
-- `ψ0::QuantumObject`: The initial state of the system.
-- `tlist::AbstractVector`: The time list of the evolution.
-- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Named Tuple of parameters to pass to the solver.
-- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `alg`: The algorithm for the ODE solver. The default value is `Tsit5()`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
@@ -215,19 +206,18 @@ where
 - `sol::TimeEvolutionSol`: The solution of the time evolution. See also [`TimeEvolutionSol`](@ref)
 """
 function mesolve(
-    H::QuantumObject{MT1,HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
+    H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
+    ψ0::QuantumObject{DT2,StateOpType},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
 ) where {
-    MT1<:AbstractMatrix,
-    T2,
+    DT1,
+    DT2,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
@@ -238,7 +228,6 @@ function mesolve(
         c_ops;
         alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         progress_bar = progress_bar,
         kwargs...,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index f2699776..87d9535b 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -16,15 +16,8 @@ function _save_func_sesolve(integrator)
     return u_modified!(integrator, false)
 end
 
-sesolve_ti_dudt!(du, u, p, t) = mul!(du, p.U, u)
-function sesolve_td_dudt!(du, u, p, t)
-    mul!(du, p.U, u)
-    H_t = p.H_t(t, p)
-    return mul!(du, H_t, u, -1im, 1)
-end
-
-function _generate_sesolve_kwargs_with_callback(t_l, kwargs)
-    cb1 = PresetTimeCallback(t_l, _save_func_sesolve, save_positions = (false, false))
+function _generate_sesolve_kwargs_with_callback(tlist, kwargs)
+    cb1 = PresetTimeCallback(tlist, _save_func_sesolve, save_positions = (false, false))
     kwargs2 =
         haskey(kwargs, :callback) ? merge(kwargs, (callback = CallbackSet(kwargs.callback, cb1),)) :
         merge(kwargs, (callback = cb1,))
@@ -32,29 +25,29 @@ function _generate_sesolve_kwargs_with_callback(t_l, kwargs)
     return kwargs2
 end
 
-function _generate_sesolve_kwargs(e_ops, progress_bar::Val{true}, t_l, kwargs)
-    return _generate_sesolve_kwargs_with_callback(t_l, kwargs)
+function _generate_sesolve_kwargs(e_ops, progress_bar::Val{true}, tlist, kwargs)
+    return _generate_sesolve_kwargs_with_callback(tlist, kwargs)
 end
 
-function _generate_sesolve_kwargs(e_ops, progress_bar::Val{false}, t_l, kwargs)
+function _generate_sesolve_kwargs(e_ops, progress_bar::Val{false}, tlist, kwargs)
     if e_ops isa Nothing
         return kwargs
     end
-    return _generate_sesolve_kwargs_with_callback(t_l, kwargs)
+    return _generate_sesolve_kwargs_with_callback(tlist, kwargs)
 end
 
 @doc raw"""
-    sesolveProblem(H::QuantumObject,
-        ψ0::QuantumObject,
+    sesolveProblem(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5()
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
+        params::NamedTuple = NamedTuple(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
 
-Generates the ODEProblem for the Schrödinger time evolution of a quantum system:
+Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 
 ```math
 \frac{\partial}{\partial t} |\psi(t)\rangle = -i \hat{H} |\psi(t)\rangle
@@ -62,22 +55,19 @@ Generates the ODEProblem for the Schrödinger time evolution of a quantum system
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: The time list of the evolution.
-- `alg::OrdinaryDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: The list of operators to be evaluated during the evolution.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
-- `params::NamedTuple`: The parameters of the system.
-- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
-- `kwargs...`: The keyword arguments passed to the `ODEProblem` constructor.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -85,73 +75,68 @@ Generates the ODEProblem for the Schrödinger time evolution of a quantum system
 - `prob`: The `ODEProblem` for the Schrödinger time evolution of the system.
 """
 function sesolveProblem(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractVector{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2}
-    H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-
+) where {DT1,DT2}
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
-    is_time_dependent = !(H_t isa Nothing)
+    tlist = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
 
-    ϕ0 = sparse_to_dense(_CType(ψ0), get_data(ψ0)) # Convert it to dense vector with complex element type
+    H_evo = QobjEvo(H, -1im) # pre-multiply by -i
+    isoper(H_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
+    check_dims(H_evo, ψ0)
 
-    t_l = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
+    ψ0 = sparse_to_dense(_CType(ψ0), get_data(ψ0)) # Convert it to dense vector with complex element type
+    U = H_evo.data
 
-    U = -1im * get_data(H)
-    progr = ProgressBar(length(t_l), enable = getVal(progress_bar))
+    progr = ProgressBar(length(tlist), enable = getVal(progress_bar))
 
     if e_ops isa Nothing
-        expvals = Array{ComplexF64}(undef, 0, length(t_l))
-        e_ops2 = MT1[]
+        expvals = Array{ComplexF64}(undef, 0, length(tlist))
+        e_ops_data = ()
         is_empty_e_ops = true
     else
-        expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
-        e_ops2 = get_data.(e_ops)
+        expvals = Array{ComplexF64}(undef, length(e_ops), length(tlist))
+        e_ops_data = get_data.(e_ops)
         is_empty_e_ops = isempty(e_ops)
     end
 
     p = (
-        U = U,
-        e_ops = e_ops2,
+        e_ops = e_ops_data,
         expvals = expvals,
         progr = progr,
-        Hdims = H.dims,
-        H_t = H_t,
-        times = t_l,
+        times = tlist,
+        Hdims = H_evo.dims,
         is_empty_e_ops = is_empty_e_ops,
         params...,
     )
 
-    saveat = is_empty_e_ops ? t_l : [t_l[end]]
+    saveat = is_empty_e_ops ? tlist : [tlist[end]]
     default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., saveat = saveat)
     kwargs2 = merge(default_values, kwargs)
-    kwargs3 = _generate_sesolve_kwargs(e_ops, makeVal(progress_bar), t_l, kwargs2)
-
-    dudt! = is_time_dependent ? sesolve_td_dudt! : sesolve_ti_dudt!
+    kwargs3 = _generate_sesolve_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2)
 
-    tspan = (t_l[1], t_l[end])
-    return ODEProblem{true,FullSpecialize}(dudt!, ϕ0, tspan, p; kwargs3...)
+    tspan = (tlist[1], tlist[end])
+    return ODEProblem{true,FullSpecialize}(U, ψ0, tspan, p; kwargs3...)
 end
 
 @doc raw"""
-    sesolve(H::QuantumObject,
-        ψ0::QuantumObject,
+    sesolve(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector;
-        alg::OrdinaryDiffEqAlgorithm=Tsit5(),
+        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
+        params::NamedTuple = NamedTuple(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
 
 Time evolution of a closed quantum system using the Schrödinger equation:
 
@@ -161,15 +146,14 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
-- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
@@ -184,27 +168,16 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 - `sol::TimeEvolutionSol`: The solution of the time evolution. See also [`TimeEvolutionSol`](@ref)
 """
 function sesolve(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractVector{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2}
-    prob = sesolveProblem(
-        H,
-        ψ0,
-        tlist;
-        alg = alg,
-        e_ops = e_ops,
-        H_t = H_t,
-        params = params,
-        progress_bar = progress_bar,
-        kwargs...,
-    )
+) where {DT1,DT2}
+    prob = sesolveProblem(H, ψ0, tlist; e_ops = e_ops, params = params, progress_bar = progress_bar, kwargs...)
 
     return sesolve(prob, alg)
 end
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 6482f433..a691ff0e 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -1,32 +1,48 @@
 export ssesolveProblem, ssesolveEnsembleProblem, ssesolve
 
-#TODO: Check if works in GPU
-function _ssesolve_update_coefficients!(ψ, coefficients, sc_ops)
-    _get_en = op -> real(dot(ψ, op, ψ)) #this is en/2: <Sn + Sn'>/2 = Re<Sn>
-    @. coefficients[2:end-1] = _get_en(sc_ops) #coefficients of the OperatorSum: Σ Sn * en/2
-    coefficients[end] = -sum(x -> x^2, coefficients[2:end-1]) / 2 #this last coefficient is -Σen^2/8
-    return nothing
+#=
+    struct DiffusionOperator
+
+A struct to represent the diffusion operator. This is used to perform the diffusion process on N different Wiener processes.
+=#
+struct DiffusionOperator{T,OT<:Tuple{Vararg{AbstractSciMLOperator}}} <: AbstractSciMLOperator{T}
+    ops::OT
+    function DiffusionOperator(ops::OT) where {OT}
+        T = mapreduce(eltype, promote_type, ops)
+        return new{T,OT}(ops)
+    end
 end
 
-function ssesolve_drift!(du, u, p, t)
-    _ssesolve_update_coefficients!(u, p.K.coefficients, p.sc_ops)
-
-    mul!(du, p.K, u)
+@generated function update_coefficients!(L::DiffusionOperator, u, p, t)
+    ops_types = L.parameters[2].parameters
+    N = length(ops_types)
+    quote
+        Base.@nexprs $N i -> begin
+            update_coefficients!(L.ops[i], u, p, t)
+        end
 
-    return nothing
+        nothing
+    end
 end
 
-function ssesolve_diffusion!(du, u, p, t)
-    @inbounds en = @view(p.K.coefficients[2:end-1])
-
-    # du:(H,W). du_reshaped:(H*W,). 
-    # H:Hilbert space dimension, W: number of sc_ops
-    du_reshaped = reshape(du, :)
-    mul!(du_reshaped, p.D, u) #du[:,i] = D[i] * u
-
-    du .-= u .* reshape(en, 1, :) #du[:,i] -= en[i] * u
+@generated function LinearAlgebra.mul!(v::AbstractVecOrMat, L::DiffusionOperator, u::AbstractVecOrMat)
+    ops_types = L.parameters[2].parameters
+    N = length(ops_types)
+    quote
+        M = length(u)
+        S = size(v)
+        (S[1] == M && S[2] == $N) || throw(DimensionMismatch("The size of the output vector is incorrect."))
+        Base.@nexprs $N i -> begin
+            mul!(@view(v[:, i]), L.ops[i], u)
+        end
+        v
+    end
+end
 
-    return nothing
+# TODO: Implement the three-argument dot function for SciMLOperators.jl
+# Currently, we are assuming a time-independent MatrixOperator
+function _ssesolve_update_coeff(u, p, t, op)
+    return real(dot(u, op.A, u)) #this is en/2: <Sn + Sn'>/2 = Re<Sn>
 end
 
 function _ssesolve_prob_func(prob, i, repeat)
@@ -37,25 +53,15 @@ function _ssesolve_prob_func(prob, i, repeat)
     traj_rng = typeof(global_rng)()
     seed!(traj_rng, seed)
 
-    noise = RealWienerProcess(
+    noise = RealWienerProcess!(
         prob.tspan[1],
-        zeros(length(internal_params.sc_ops)),
-        zeros(length(internal_params.sc_ops)),
+        zeros(internal_params.n_sc_ops),
+        zeros(internal_params.n_sc_ops),
         save_everystep = false,
         rng = traj_rng,
     )
 
-    noise_rate_prototype = similar(prob.u0, length(prob.u0), length(internal_params.sc_ops))
-
-    prm = merge(
-        internal_params,
-        (
-            expvals = similar(internal_params.expvals),
-            progr = ProgressBar(size(internal_params.expvals, 2), enable = false),
-        ),
-    )
-
-    return remake(prob, p = prm, noise = noise, noise_rate_prototype = noise_rate_prototype, seed = seed)
+    return remake(prob, noise = noise, seed = seed)
 end
 
 # Standard output function
@@ -86,58 +92,62 @@ function _ssesolve_generate_statistics!(sol, i, states, expvals_all)
     return nothing
 end
 
+_ScalarOperator_e(op, f = +) = ScalarOperator(one(eltype(op)), (a, u, p, t) -> f(_ssesolve_update_coeff(u, p, t, op)))
+
+_ScalarOperator_e2_2(op, f = +) =
+    ScalarOperator(one(eltype(op)), (a, u, p, t) -> f(_ssesolve_update_coeff(u, p, t, op)^2 / 2))
+
 @doc raw"""
-    ssesolveProblem(H::QuantumObject,
-        ψ0::QuantumObject,
-        tlist::AbstractVector;
-        sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
-        alg::StochasticDiffEqAlgorithm=SRA1()
+    ssesolveProblem(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
+        tlist::AbstractVector,
+        sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        rng::AbstractRNG=default_rng(),
-        kwargs...)
+        params::NamedTuple = NamedTuple(),
+        rng::AbstractRNG = default_rng(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
 
-Generates the SDEProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
+Generate the SDEProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
     
-    ```math
-    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
-    ```
+```math
+d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+```
 
 where 
     
 ```math
-    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
-    ```
-    ```math
-    M_n = C_n - \frac{e_n}{2},
-    ```
-    and
-    ```math
-    e_n = \langle C_n + C_n^\dagger \rangle.
-    ```
+K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+```
+```math
+M_n = C_n - \frac{e_n}{2},
+```
+and
+```math
+e_n = \langle C_n + C_n^\dagger \rangle.
+```
 
-Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
+Above, `C_n` is the `n`-th collapse operator and `dW_j(t)` is the real Wiener increment associated to `C_n`.
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: The time list of the evolution.
-- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of operators to be evaluated during the evolution.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
-- `params::NamedTuple`: The parameters of the system.
-- `rng::AbstractRNG`: The random number generator for reproducibility.
-- `kwargs...`: The keyword arguments passed to the `SDEProblem` constructor.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `sc_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -145,146 +155,139 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 - `prob`: The `SDEProblem` for the Stochastic Schrödinger time evolution of the system.
 """
 function ssesolveProblem(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::StochasticDiffEqAlgorithm = SRA1(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
+    progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2}
-    H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-
+) where {DT1,DT2}
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
     sc_ops isa Nothing &&
         throw(ArgumentError("The list of collapse operators must be provided. Use sesolveProblem instead."))
 
-    !(H_t isa Nothing) && throw(ArgumentError("Time-dependent Hamiltonians are not currently supported in ssesolve."))
+    tlist = convert(Vector{Float64}, tlist) # Convert it into Float64 to avoid type instabilities for StochasticDiffEq.jl
 
-    t_l = convert(Vector{Float64}, tlist) # Convert it into Float64 to avoid type instabilities for StochasticDiffEq.jl
+    H_eff_evo = _mcsolve_make_Heff_QobjEvo(H, sc_ops)
+    isoper(H_eff_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
+    check_dims(H_eff_evo, ψ0)
+    dims = H_eff_evo.dims
 
-    ϕ0 = get_data(ψ0)
+    ψ0 = sparse_to_dense(_CType(ψ0), get_data(ψ0))
 
-    H_eff = get_data(H - T2(0.5im) * mapreduce(op -> op' * op, +, sc_ops))
-    sc_ops2 = get_data.(sc_ops)
-
-    coefficients = [1.0, fill(0.0, length(sc_ops) + 1)...]
-    operators = [-1im * H_eff, sc_ops2..., MT1(I(prod(H.dims)))]
-    K = OperatorSum(coefficients, operators)
-    _ssesolve_update_coefficients!(ϕ0, K.coefficients, sc_ops2)
-
-    D = reduce(vcat, sc_ops2)
+    progr = ProgressBar(length(tlist), enable = getVal(progress_bar))
 
     if e_ops isa Nothing
-        expvals = Array{ComplexF64}(undef, 0, length(t_l))
-        e_ops2 = MT1[]
+        expvals = Array{ComplexF64}(undef, 0, length(tlist))
+        e_ops_data = ()
         is_empty_e_ops = true
     else
-        expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
-        e_ops2 = get_data.(e_ops)
+        expvals = Array{ComplexF64}(undef, length(e_ops), length(tlist))
+        e_ops_data = get_data.(e_ops)
         is_empty_e_ops = isempty(e_ops)
     end
 
+    sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops))
+
+    # Here the coefficients depend on the state, so this is a non-linear operator, which should be implemented with FunctionOperator instead. However, the nonlinearity is only on the coefficients, and it should be safe.
+    K_l = sum(
+        op -> _ScalarOperator_e(op, +) * op + _ScalarOperator_e2_2(op, -) * IdentityOperator(prod(dims)),
+        sc_ops_evo_data,
+    )
+
+    K = -1im * get_data(H_eff_evo) + K_l
+
+    D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
+    D = DiffusionOperator(D_l)
+
     p = (
-        K = K,
-        D = D,
-        e_ops = e_ops2,
-        sc_ops = sc_ops2,
+        e_ops = e_ops_data,
         expvals = expvals,
-        Hdims = H.dims,
-        H_t = H_t,
-        times = t_l,
+        progr = progr,
+        times = tlist,
+        Hdims = dims,
         is_empty_e_ops = is_empty_e_ops,
+        n_sc_ops = length(sc_ops),
         params...,
     )
 
-    saveat = is_empty_e_ops ? t_l : [t_l[end]]
+    saveat = is_empty_e_ops ? tlist : [tlist[end]]
     default_values = (DEFAULT_SDE_SOLVER_OPTIONS..., saveat = saveat)
     kwargs2 = merge(default_values, kwargs)
-    kwargs3 = _generate_sesolve_kwargs(e_ops, Val(false), t_l, kwargs2)
-
-    tspan = (t_l[1], t_l[end])
-    noise = RealWienerProcess(t_l[1], zeros(length(sc_ops)), zeros(length(sc_ops)), save_everystep = false, rng = rng)
-    noise_rate_prototype = similar(ϕ0, length(ϕ0), length(sc_ops))
-    return SDEProblem{true}(
-        ssesolve_drift!,
-        ssesolve_diffusion!,
-        ϕ0,
-        tspan,
-        p;
-        noise_rate_prototype = noise_rate_prototype,
-        noise = noise,
-        kwargs3...,
-    )
+    kwargs3 = _generate_sesolve_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2)
+
+    tspan = (tlist[1], tlist[end])
+    noise =
+        RealWienerProcess!(tlist[1], zeros(length(sc_ops)), zeros(length(sc_ops)), save_everystep = false, rng = rng)
+    noise_rate_prototype = similar(ψ0, length(ψ0), length(sc_ops))
+    return SDEProblem{true}(K, D, ψ0, tspan, p; noise_rate_prototype = noise_rate_prototype, noise = noise, kwargs3...)
 end
 
 @doc raw"""
-    ssesolveEnsembleProblem(H::QuantumObject,
-        ψ0::QuantumObject,
-        tlist::AbstractVector;
+    ssesolveEnsembleProblem(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
+        tlist::AbstractVector,
         sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-        alg::StochasticDiffEqAlgorithm=SRA1()
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        rng::AbstractRNG=default_rng(),
-        ntraj::Int=1,
-        ensemble_method=EnsembleThreads(),
-        prob_func::Function=_mcsolve_prob_func,
-        output_func::Function=_ssesolve_dispatch_output_func(ensemble_method),
-        progress_bar::Union{Val,Bool}=Val(true),
-        kwargs...)
-
-Generates the SDE EnsembleProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
+        params::NamedTuple = NamedTuple(),
+        rng::AbstractRNG = default_rng(),
+        ntraj::Int = 1,
+        ensemble_method = EnsembleThreads(),
+        prob_func::Function = _ssesolve_prob_func,
+        output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
+        progress_bar::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
+
+Generate the SDE EnsembleProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
     
-    ```math
-    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
-    ```
+```math
+d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+```
 
 where 
     
 ```math
-    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
-    ```
-    ```math
-    M_n = C_n - \frac{e_n}{2},
-    ```
-    and
-    ```math
-    e_n = \langle C_n + C_n^\dagger \rangle.
-    ```
+K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+```
+```math
+M_n = C_n - \frac{e_n}{2},
+```
+and
+```math
+e_n = \langle C_n + C_n^\dagger \rangle.
+```
 
 Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
 
 # Arguments
 
-- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: The time list of the evolution.
-- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of operators to be evaluated during the evolution.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
-- `params::NamedTuple`: The parameters of the system.
-- `rng::AbstractRNG`: The random number generator for reproducibility.
-- `ntraj::Int`: Number of trajectories to use.
-- `ensemble_method`: Ensemble method to use.
-- `prob_func::Function`: Function to use for generating the SDEProblem.
-- `output_func::Function`: Function to use for generating the output of a single trajectory.
-- `progress_bar::Union{Val,Bool}`: Whether to show a progress bar.
-- `kwargs...`: The keyword arguments passed to the `SDEProblem` constructor.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `sc_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `ntraj`: Number of trajectories to use.
+- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`.
+- `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
+- `prob_func`: Function to use for generating the ODEProblem.
+- `output_func`: Function to use for generating the output of a single trajectory.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
-- For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
 # Returns
@@ -292,13 +295,11 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 - `prob::EnsembleProblem with SDEProblem`: The Ensemble SDEProblem for the Stochastic Shrödinger time evolution.
 """
 function ssesolveEnsembleProblem(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::StochasticDiffEqAlgorithm = SRA1(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
     ntraj::Int = 1,
@@ -307,7 +308,7 @@ function ssesolveEnsembleProblem(
     output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2}
+) where {DT1,DT2}
     progr = ProgressBar(ntraj, enable = getVal(progress_bar))
     if ensemble_method isa EnsembleDistributed
         progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
@@ -327,14 +328,15 @@ function ssesolveEnsembleProblem(
             ψ0,
             tlist,
             sc_ops;
-            alg = alg,
             e_ops = e_ops,
-            H_t = H_t,
             params = merge(params, (global_rng = rng, seeds = seeds)),
             rng = rng,
+            progress_bar = Val(false),
             kwargs...,
         )
 
+        # safetycopy is set to true because the K and D functions cannot be currently deepcopied.
+        # the memory overhead shouldn't be too large, compared to the safetycopy=false case.
         ensemble_prob = EnsembleProblem(prob_sse, prob_func = prob_func, output_func = output_func, safetycopy = true)
 
         return ensemble_prob
@@ -347,67 +349,66 @@ function ssesolveEnsembleProblem(
 end
 
 @doc raw"""
-    ssesolve(H::QuantumObject,
-        ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    ssesolve(
+        H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+        ψ0::QuantumObject{DT2,KetQuantumObject},
         tlist::AbstractVector,
-        sc_ops::Union{Nothing, AbstractVector}=nothing;
-        alg::StochasticDiffEqAlgorithm=SRA1(),
-        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
-        params::NamedTuple=NamedTuple(),
-        rng::AbstractRNG=default_rng(),
-        ntraj::Int=1,
-        ensemble_method=EnsembleThreads(),
-        prob_func::Function=_ssesolve_prob_func,
-        output_func::Function=_ssesolve_dispatch_output_func(ensemble_method),
+        sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
+        alg::StochasticDiffEqAlgorithm = SRA1(),
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        params::NamedTuple = NamedTuple(),
+        rng::AbstractRNG = default_rng(),
+        ntraj::Int = 1,
+        ensemble_method = EnsembleThreads(),
+        prob_func::Function = _ssesolve_prob_func,
+        output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
         progress_bar::Union{Val,Bool} = Val(true),
-        kwargs...)
+        kwargs...,
+    )
 
 
 Stochastic Schrödinger equation evolution of a quantum system given the system Hamiltonian ``\hat{H}`` and a list of stochadtic collapse (jump) operators ``\{\hat{C}_n\}_n``.
 The stochastic evolution of the state ``|\psi(t)\rangle`` is defined by:
     
-    ```math
-    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
-    ```
+```math
+d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+```
 
 where 
     
 ```math
-    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
-    ```
-    ```math
-    M_n = C_n - \frac{e_n}{2},
-    ```
-    and
-    ```math
-    e_n = \langle C_n + C_n^\dagger \rangle.
-    ```
+K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+```
+```math
+M_n = C_n - \frac{e_n}{2},
+```
+and
+```math
+e_n = \langle C_n + C_n^\dagger \rangle.
+```
 
-Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
+Above, `C_n` is the `n`-th collapse operator and `dW_j(t)` is the real Wiener increment associated to `C_n`.
 
 
 # Arguments
 
-- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
-- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
-- `alg::StochasticDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
-- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
-- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
-- `ntraj::Int`: Number of trajectories to use.
-- `ensemble_method`: Ensemble method to use.
-- `prob_func::Function`: Function to use for generating the SDEProblem.
-- `output_func::Function`: Function to use for generating the output of a single trajectory.
-- `progress_bar::Union{Val,Bool}`: Whether to show a progress bar.
-- `kwargs...`: Additional keyword arguments to pass to the solver.
+- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `sc_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
+- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRA1()`.
+- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
+- `params`: `NamedTuple` of parameters to pass to the solver.
+- `rng`: Random number generator for reproducibility.
+- `ntraj`: Number of trajectories to use.
+- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`.
+- `prob_func`: Function to use for generating the ODEProblem.
+- `output_func`: Function to use for generating the output of a single trajectory.
+- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
 
-- `ensemble_method` can be one of `EnsembleThreads()`, `EnsembleSerial()`, `EnsembleDistributed()`
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
@@ -416,16 +417,15 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 
 # Returns
 
-- `sol::TimeEvolutionSSESol`: The solution of the time evolution. See also [`TimeEvolutionSSESol`](@ref)
+- `sol::TimeEvolutionSSESol`: The solution of the time evolution. See also [`TimeEvolutionSSESol`](@ref).
 """
 function ssesolve(
-    H::QuantumObject{MT1,OperatorQuantumObject},
-    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
+    ψ0::QuantumObject{DT2,KetQuantumObject},
     tlist::AbstractVector,
     sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::StochasticDiffEqAlgorithm = SRA1(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     rng::AbstractRNG = default_rng(),
     ntraj::Int = 1,
@@ -434,21 +434,13 @@ function ssesolve(
     output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {MT1<:AbstractMatrix,T2}
-    progr = ProgressBar(ntraj, enable = getVal(progress_bar))
-    progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
-    @async while take!(progr_channel)
-        next!(progr)
-    end
-
+) where {DT1,DT2}
     ens_prob = ssesolveEnsembleProblem(
         H,
         ψ0,
         tlist,
         sc_ops;
-        alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         rng = rng,
         ntraj = ntraj,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 971d53c6..90284b50 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -1,4 +1,3 @@
-export TimeDependentOperatorSum
 export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionSSESol
 
 export liouvillian, liouvillian_floquet, liouvillian_generalized
@@ -21,14 +20,22 @@ A structure storing the results and some information from solving time evolution
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
-struct TimeEvolutionSol{TT<:Vector{<:Real},TS<:AbstractVector,TE<:Matrix{ComplexF64}}
+struct TimeEvolutionSol{
+    TT<:AbstractVector{<:Real},
+    TS<:AbstractVector,
+    TE<:Matrix,
+    RETT<:Enum,
+    AlgT<:OrdinaryDiffEqAlgorithm,
+    AT<:Real,
+    RT<:Real,
+}
     times::TT
     states::TS
     expect::TE
-    retcode::Enum
-    alg::OrdinaryDiffEqAlgorithm
-    abstol::Real
-    reltol::Real
+    retcode::RETT
+    alg::AlgT
+    abstol::AT
+    reltol::RT
 end
 
 function Base.show(io::IO, sol::TimeEvolutionSol)
@@ -63,12 +70,15 @@ A structure storing the results and some information from solving quantum trajec
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionMCSol{
-    TT<:Vector{<:Real},
+    TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     TEA<:Array{ComplexF64,3},
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
+    AlgT<:OrdinaryDiffEqAlgorithm,
+    AT<:Real,
+    RT<:Real,
 }
     ntraj::Int
     times::TT
@@ -78,9 +88,9 @@ struct TimeEvolutionMCSol{
     jump_times::TJT
     jump_which::TJW
     converged::Bool
-    alg::OrdinaryDiffEqAlgorithm
-    abstol::Real
-    reltol::Real
+    alg::AlgT
+    abstol::AT
+    reltol::RT
 end
 
 function Base.show(io::IO, sol::TimeEvolutionMCSol)
@@ -114,12 +124,13 @@ A structure storing the results and some information from solving trajectories o
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionSSESol{
-    TT<:Vector{<:Real},
+    TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     TEA<:Array{ComplexF64,3},
-    T1<:Real,
-    T2<:Real,
+    AlgT<:StochasticDiffEqAlgorithm,
+    AT<:Real,
+    RT<:Real,
 }
     ntraj::Int
     times::TT
@@ -127,9 +138,9 @@ struct TimeEvolutionSSESol{
     expect::TE
     expect_all::TEA
     converged::Bool
-    alg::StochasticDiffEqAlgorithm
-    abstol::T1
-    reltol::T2
+    alg::AlgT
+    abstol::AT
+    reltol::RT
 end
 
 function Base.show(io::IO, sol::TimeEvolutionSSESol)
@@ -155,86 +166,8 @@ struct DiscreteLindbladJumpCallback <: LindbladJumpCallbackType end
 
 ContinuousLindbladJumpCallback(; interp_points::Int = 10) = ContinuousLindbladJumpCallback(interp_points)
 
-## Time-dependent sum of operators
-
-struct TimeDependentOperatorSum{CFT,OST<:OperatorSum}
-    coefficient_functions::CFT
-    operator_sum::OST
-end
-
-function TimeDependentOperatorSum(
-    coefficient_functions,
-    operators::Union{AbstractVector{<:QuantumObject},Tuple};
-    params = nothing,
-    init_time = 0.0,
-)
-    # promote the type of the coefficients and the operators. Remember that the coefficient_functions si a vector of functions and the operators is a vector of QuantumObjects
-    coefficients = [f(init_time, params) for f in coefficient_functions]
-    operator_sum = OperatorSum(coefficients, operators)
-    return TimeDependentOperatorSum(coefficient_functions, operator_sum)
-end
-
-Base.size(A::TimeDependentOperatorSum) = size(A.operator_sum)
-Base.size(A::TimeDependentOperatorSum, inds...) = size(A.operator_sum, inds...)
-Base.length(A::TimeDependentOperatorSum) = length(A.operator_sum)
-
-function update_coefficients!(A::TimeDependentOperatorSum, t, params)
-    @inbounds @simd for i in 1:length(A.coefficient_functions)
-        A.operator_sum.coefficients[i] = A.coefficient_functions[i](t, params)
-    end
-end
-
-(A::TimeDependentOperatorSum)(t, params) = (update_coefficients!(A, t, params); A)
-
-@inline function LinearAlgebra.mul!(y::AbstractVector, A::TimeDependentOperatorSum, x::AbstractVector, α, β)
-    return mul!(y, A.operator_sum, x, α, β)
-end
-
-function liouvillian(A::TimeDependentOperatorSum, Id_cache = I(prod(A.operator_sum.operators[1].dims)))
-    return TimeDependentOperatorSum(A.coefficient_functions, liouvillian(A.operator_sum, Id_cache))
-end
-
 #######################################
 
-### LIOUVILLIAN ###
-@doc raw"""
-    liouvillian(H::QuantumObject, 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``:
-
-```math
-\mathcal{L} [\cdot] = -i[\hat{H}, \cdot] + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
-```
-
-where 
-
-```math
-\mathcal{D}(\hat{C}_n) [\cdot] = \hat{C}_n [\cdot] \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n [\cdot] - \frac{1}{2} [\cdot] \hat{C}_n^\dagger \hat{C}_n
-```
-
-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), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
-"""
-function liouvillian(
-    H::QuantumObject{MT1,OpType1},
-    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    Id_cache = I(prod(H.dims)),
-) where {MT1<:AbstractMatrix,OpType1<: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
-    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::QuantumObject{<:AbstractMatrix,SuperOperatorQuantumObject}, Id_cache::Diagonal) = H
-
 function liouvillian_floquet(
     L₀::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
     Lₚ::QuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 4ffe3187..8830de40 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -129,25 +129,25 @@ function _DFDIncreaseReduceAffect!(integrator)
 
     resize!(integrator, size(L, 1))
     copyto!(integrator.u, mat2vec(ρt))
-    integrator.p = merge(internal_params, (L = L, e_ops = e_ops2, dfd_ρt_cache = similar(integrator.u)))
+    # By doing this, we are assuming that the system is time-independent and f is a MatrixOperator
+    integrator.f = ODEFunction{true,FullSpecialize}(MatrixOperator(L))
+    integrator.p = merge(internal_params, (e_ops = e_ops2, dfd_ρt_cache = similar(integrator.u)))
 
     return nothing
 end
 
 function dfd_mesolveProblem(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T1},StateOpType},
-    t_l::AbstractVector,
+    ψ0::QuantumObject{DT1,StateOpType},
+    tlist::AbstractVector,
     c_ops::Function,
     maxdims::Vector{T2},
     dfd_params::NamedTuple = NamedTuple();
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (dim_list) -> Vector{Vector{T1}}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (dim_list) -> Vector{Vector{DT1}}([]),
     params::NamedTuple = NamedTuple(),
     tol_list::Vector{<:Number} = fill(1e-8, length(maxdims)),
     kwargs...,
-) where {T1,T2<:Integer,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
+) where {DT1,T2<:Integer,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
     length(ψ0.dims) != length(maxdims) &&
         throw(DimensionMismatch("'dim_list' and 'maxdims' do not have the same dimension."))
 
@@ -158,8 +158,8 @@ function dfd_mesolveProblem(
 
     dim_list_evo_times = [0.0]
     dim_list_evo = [dim_list]
-    reduce_list = MVector(ntuple(i -> false, length(dim_list)))
-    increase_list = MVector(ntuple(i -> false, length(dim_list)))
+    reduce_list = MVector(ntuple(i -> false, Val(length(dim_list))))
+    increase_list = MVector(ntuple(i -> false, Val(length(dim_list))))
     pillow_list = _dfd_set_pillow.(dim_list)
 
     params2 = merge(
@@ -187,16 +187,15 @@ function dfd_mesolveProblem(
         haskey(kwargs2, :callback) ? merge(kwargs2, (callback = CallbackSet(cb_dfd, kwargs2.callback),)) :
         merge(kwargs2, (callback = cb_dfd,))
 
-    return mesolveProblem(H₀, ψ0, t_l, c_ops₀; e_ops = e_ops₀, alg = alg, H_t = H_t, params = params2, kwargs2...)
+    return mesolveProblem(H₀, ψ0, tlist, c_ops₀; e_ops = e_ops₀, params = params2, kwargs2...)
 end
 
 @doc raw"""
     dfd_mesolve(H::Function, ψ0::QuantumObject,
-        t_l::AbstractVector, c_ops::Function, maxdims::AbstractVector,
+        tlist::AbstractVector, c_ops::Function, maxdims::AbstractVector,
         dfd_params::NamedTuple=NamedTuple();
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
         e_ops::Function=(dim_list) -> Vector{Vector{T1}}([]),
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         tol_list::Vector{<:Number}=fill(1e-8, length(maxdims)),
         kwargs...)
@@ -210,13 +209,12 @@ Time evolution of an open quantum system using master equation, dynamically chan
 function dfd_mesolve(
     H::Function,
     ψ0::QuantumObject{<:AbstractArray{T1},StateOpType},
-    t_l::AbstractVector,
+    tlist::AbstractVector,
     c_ops::Function,
     maxdims::Vector{T2},
     dfd_params::NamedTuple = NamedTuple();
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     e_ops::Function = (dim_list) -> Vector{Vector{T1}}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     tol_list::Vector{<:Number} = fill(1e-8, length(maxdims)),
     kwargs...,
@@ -224,13 +222,11 @@ function dfd_mesolve(
     dfd_prob = dfd_mesolveProblem(
         H,
         ψ0,
-        t_l,
+        tlist,
         c_ops,
         maxdims,
         dfd_params;
-        alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         tol_list = tol_list,
         kwargs...,
@@ -279,7 +275,7 @@ end
 
 function _DSF_mesolve_Affect!(integrator)
     internal_params = integrator.p
-    op_l = internal_params.op_l
+    op_list = internal_params.op_list
     op_l_vec = internal_params.op_l_vec
     αt_list = internal_params.αt_list
     δα_list = internal_params.δα_list
@@ -293,11 +289,10 @@ function _DSF_mesolve_Affect!(integrator)
     dsf_identity = internal_params.dsf_identity
     dsf_displace_cache_full = internal_params.dsf_displace_cache_full
 
-    op_l_length = length(op_l)
-    fill!(dsf_displace_cache_full.coefficients, 0)
+    op_l_length = length(op_list)
 
-    for i in eachindex(op_l)
-        # op = op_l[i]
+    for i in eachindex(op_list)
+        # op = op_list[i]
         op_vec = op_l_vec[i]
         αt = αt_list[i]
         δα = δα_list[i]
@@ -318,12 +313,17 @@ function _DSF_mesolve_Affect!(integrator)
             # arnoldi!(expv_cache, Aᵢ, dsf_cache)
             # expv!(integrator.u, expv_cache, one(αt), dsf_cache)
 
-            dsf_displace_cache_full.coefficients[i] = Δα
-            dsf_displace_cache_full.coefficients[i+op_l_length] = -conj(Δα)
-            dsf_displace_cache_full.coefficients[i+2*op_l_length] = conj(Δα)
-            dsf_displace_cache_full.coefficients[i+3*op_l_length] = -Δα
+            dsf_displace_cache_full.ops[i].λ.val = Δα
+            dsf_displace_cache_full.ops[i+op_l_length].λ.val = -conj(Δα)
+            dsf_displace_cache_full.ops[i+2*op_l_length].λ.val = conj(Δα)
+            dsf_displace_cache_full.ops[i+3*op_l_length].λ.val = -Δα
 
             αt_list[i] += Δα
+        else
+            dsf_displace_cache_full.ops[i].λ.val = 0
+            dsf_displace_cache_full.ops[i+op_l_length].λ.val = 0
+            dsf_displace_cache_full.ops[i+2*op_l_length].λ.val = 0
+            dsf_displace_cache_full.ops[i+3*op_l_length].λ.val = 0
         end
     end
 
@@ -331,52 +331,47 @@ function _DSF_mesolve_Affect!(integrator)
     arnoldi!(expv_cache, dsf_displace_cache_full, dsf_cache)
     expv!(integrator.u, expv_cache, 1, dsf_cache)
 
-    op_l2 = op_l .+ αt_list
+    op_l2 = op_list .+ αt_list
     e_ops2 = e_ops(op_l2, dsf_params)
     _mat2vec_data = op -> mat2vec(get_data(op)')
     @. e_ops_vec = _mat2vec_data(e_ops2)
-    return copyto!(internal_params.L, liouvillian(H(op_l2, dsf_params), c_ops(op_l2, dsf_params), dsf_identity).data)
+    # By doing this, we are assuming that the system is time-independent and f is a MatrixOperator
+    copyto!(integrator.f.f.A, liouvillian(H(op_l2, dsf_params), c_ops(op_l2, dsf_params), dsf_identity).data)
+    return u_modified!(integrator, true)
 end
 
 function dsf_mesolveProblem(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T},StateOpType},
-    t_l::AbstractVector,
+    ψ0::QuantumObject{<:AbstractVector{T},StateOpType},
+    tlist::AbstractVector,
     c_ops::Function,
-    op_list::Vector{TOl},
+    op_list::Union{AbstractVector,Tuple},
     α0_l::Vector{<:Number} = zeros(length(op_list)),
     dsf_params::NamedTuple = NamedTuple();
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (op_list, p) -> Vector{TOl}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (op_list, p) -> (),
     params::NamedTuple = NamedTuple(),
     δα_list::Vector{<:Real} = fill(0.2, length(op_list)),
     krylov_dim::Int = max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
     kwargs...,
-) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},TOl}
-    op_l = op_list
-    H₀ = H(op_l .+ α0_l, dsf_params)
-    c_ops₀ = c_ops(op_l .+ α0_l, dsf_params)
-    e_ops₀ = e_ops(op_l .+ α0_l, dsf_params)
+) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
+    op_list = deepcopy(op_list)
+    H₀ = H(op_list .+ α0_l, dsf_params)
+    c_ops₀ = c_ops(op_list .+ α0_l, dsf_params)
+    e_ops₀ = e_ops(op_list .+ α0_l, dsf_params)
 
     αt_list = convert(Vector{T}, α0_l)
-    op_l_vec = map(op -> mat2vec(get_data(op)'), op_l)
+    op_l_vec = map(op -> mat2vec(get_data(op)'), op_list)
     # Create the Krylov subspace with kron(H₀.data, H₀.data) just for initialize
     expv_cache = arnoldi(kron(H₀.data, H₀.data), mat2vec(ket2dm(ψ0).data), krylov_dim)
     dsf_identity = I(prod(H₀.dims))
-    dsf_displace_cache_left = map(op -> Qobj(kron(op.data, dsf_identity)), op_l)
-    dsf_displace_cache_left_dag = map(op -> Qobj(kron(sparse(op.data'), dsf_identity)), op_l)
-    dsf_displace_cache_right = map(op -> Qobj(kron(dsf_identity, op.data)), op_l)
-    dsf_displace_cache_right_dag = map(op -> Qobj(kron(dsf_identity, sparse(op.data'))), op_l)
-    dsf_displace_cache_full = OperatorSum(
-        zeros(length(op_l) * 4),
-        vcat(
-            dsf_displace_cache_left,
-            dsf_displace_cache_left_dag,
-            dsf_displace_cache_right,
-            dsf_displace_cache_right_dag,
-        ),
-    )
+    dsf_displace_cache_left = sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(op.data, dsf_identity)), op_list)
+    dsf_displace_cache_left_dag =
+        sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(sparse(op.data'), dsf_identity)), op_list)
+    dsf_displace_cache_right = sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(dsf_identity, op.data)), op_list)
+    dsf_displace_cache_right_dag =
+        sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(dsf_identity, sparse(op.data'))), op_list)
+    dsf_displace_cache_full =
+        dsf_displace_cache_left + dsf_displace_cache_left_dag + dsf_displace_cache_right + dsf_displace_cache_right_dag
 
     params2 = params
     params2 = merge(
@@ -385,7 +380,7 @@ function dsf_mesolveProblem(
             H_fun = H,
             c_ops_fun = c_ops,
             e_ops_fun = e_ops,
-            op_l = op_l,
+            op_list = op_list,
             op_l_vec = op_l_vec,
             αt_list = αt_list,
             δα_list = δα_list,
@@ -403,19 +398,18 @@ function dsf_mesolveProblem(
         haskey(kwargs2, :callback) ? merge(kwargs2, (callback = CallbackSet(cb_dsf, kwargs2.callback),)) :
         merge(kwargs2, (callback = cb_dsf,))
 
-    return mesolveProblem(H₀, ψ0, t_l, c_ops₀; e_ops = e_ops₀, alg = alg, H_t = H_t, params = params2, kwargs2...)
+    return mesolveProblem(H₀, ψ0, tlist, c_ops₀; e_ops = e_ops₀, params = params2, kwargs2...)
 end
 
 @doc raw"""
     dsf_mesolve(H::Function,
         ψ0::QuantumObject,
-        t_l::AbstractVector, c_ops::Function,
+        tlist::AbstractVector, c_ops::Function,
         op_list::Vector{TOl},
         α0_l::Vector{<:Number}=zeros(length(op_list)),
         dsf_params::NamedTuple=NamedTuple();
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
         e_ops::Function=(op_list,p) -> Vector{TOl}([]),
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         δα_list::Vector{<:Number}=fill(0.2, length(op_list)),
         krylov_dim::Int=max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
@@ -429,31 +423,28 @@ Time evolution of an open quantum system using master equation and the Dynamical
 """
 function dsf_mesolve(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T},StateOpType},
-    t_l::AbstractVector,
+    ψ0::QuantumObject{<:AbstractVector{T},StateOpType},
+    tlist::AbstractVector,
     c_ops::Function,
-    op_list::Vector{TOl},
+    op_list::Union{AbstractVector,Tuple},
     α0_l::Vector{<:Number} = zeros(length(op_list)),
     dsf_params::NamedTuple = NamedTuple();
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (op_list, p) -> Vector{TOl}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (op_list, p) -> (),
     params::NamedTuple = NamedTuple(),
     δα_list::Vector{<:Real} = fill(0.2, length(op_list)),
     krylov_dim::Int = max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
     kwargs...,
-) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},TOl}
+) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
     dsf_prob = dsf_mesolveProblem(
         H,
         ψ0,
-        t_l,
+        tlist,
         c_ops,
         op_list,
         α0_l,
         dsf_params;
-        alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         δα_list = δα_list,
         krylov_dim = krylov_dim,
@@ -465,31 +456,29 @@ end
 
 function dsf_mesolve(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T},StateOpType},
-    t_l::AbstractVector,
-    op_list::Vector{TOl},
+    ψ0::QuantumObject{<:AbstractVector{T},StateOpType},
+    tlist::AbstractVector,
+    op_list::Union{AbstractVector,Tuple},
     α0_l::Vector{<:Number} = zeros(length(op_list)),
     dsf_params::NamedTuple = NamedTuple();
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (op_list, p) -> Vector{TOl}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (op_list, p) -> (),
     params::NamedTuple = NamedTuple(),
     δα_list::Vector{<:Real} = fill(0.2, length(op_list)),
     krylov_dim::Int = max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
     kwargs...,
-) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},TOl}
-    c_ops = op_list -> Vector{TOl}([])
+) where {T,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
+    c_ops = op_list -> ()
     return dsf_mesolve(
         H,
         ψ0,
-        t_l,
+        tlist,
         c_ops,
         op_list,
         α0_l,
         dsf_params;
         alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         δα_list = δα_list,
         krylov_dim = krylov_dim,
@@ -501,14 +490,14 @@ end
 
 function _DSF_mcsolve_Condition(u, t, integrator)
     internal_params = integrator.p
-    op_l = internal_params.op_l
+    op_list = internal_params.op_list
     δα_list = internal_params.δα_list
 
     ψt = u
 
     condition = false
-    @inbounds for i in eachindex(op_l)
-        op = op_l[i]
+    @inbounds for i in eachindex(op_list)
+        op = op_list[i]
         δα = δα_list[i]
         Δα = dot(ψt, op.data, ψt) / dot(ψt, ψt)
         if δα < abs(Δα)
@@ -520,7 +509,7 @@ end
 
 function _DSF_mcsolve_Affect!(integrator)
     internal_params = integrator.p
-    op_l = internal_params.op_l
+    op_list = internal_params.op_list
     αt_list = internal_params.αt_list
     δα_list = internal_params.δα_list
     H = internal_params.H_fun
@@ -537,11 +526,10 @@ function _DSF_mcsolve_Affect!(integrator)
     copyto!(ψt, integrator.u)
     normalize!(ψt)
 
-    op_l_length = length(op_l)
-    fill!(dsf_displace_cache_full.coefficients, 0)
+    op_l_length = length(op_list)
 
-    for i in eachindex(op_l)
-        op = op_l[i]
+    for i in eachindex(op_list)
+        op = op_list[i]
         αt = αt_list[i]
         δα = δα_list[i]
         Δα = dot(ψt, op.data, ψt)
@@ -556,10 +544,13 @@ function _DSF_mcsolve_Affect!(integrator)
             # arnoldi!(expv_cache, Aᵢ, dsf_cache)
             # expv!(integrator.u, expv_cache, one(αt), dsf_cache)
 
-            dsf_displace_cache_full.coefficients[i] = conj(Δα)
-            dsf_displace_cache_full.coefficients[i+op_l_length] = -Δα
+            dsf_displace_cache_full.ops[i].λ.val = conj(Δα)
+            dsf_displace_cache_full.ops[i+op_l_length].λ.val = -Δα
 
             αt_list[i] += Δα
+        else
+            dsf_displace_cache_full.ops[i].λ.val = 0
+            dsf_displace_cache_full.ops[i+op_l_length].λ.val = 0
         end
     end
 
@@ -567,13 +558,15 @@ function _DSF_mcsolve_Affect!(integrator)
     arnoldi!(expv_cache, dsf_displace_cache_full, dsf_cache)
     expv!(integrator.u, expv_cache, 1, dsf_cache)
 
-    op_l2 = op_l .+ αt_list
+    op_l2 = op_list .+ αt_list
     e_ops2 = e_ops(op_l2, dsf_params)
     c_ops2 = c_ops(op_l2, dsf_params)
     @. e_ops0 = get_data(e_ops2)
     @. c_ops0 = get_data(c_ops2)
-    H_eff = H(op_l2, dsf_params).data - lmul!(convert(eltype(ψt), 0.5im), mapreduce(op -> op' * op, +, c_ops0))
-    return mul!(internal_params.U, -1im, H_eff)
+    H_nh = lmul!(convert(eltype(ψt), 0.5im), mapreduce(op -> op' * op, +, c_ops0))
+    # By doing this, we are assuming that the system is time-independent and f is a ScaledOperator
+    copyto!(integrator.f.f.L.A, H(op_l2, dsf_params).data - H_nh)
+    return u_modified!(integrator, true)
 end
 
 function _dsf_mcsolve_prob_func(prob, i, repeat)
@@ -582,9 +575,8 @@ function _dsf_mcsolve_prob_func(prob, i, repeat)
     prm = merge(
         internal_params,
         (
-            U = copy(internal_params.U),
-            e_ops_mc = copy(internal_params.e_ops_mc),
-            c_ops = copy(internal_params.c_ops),
+            e_ops_mc = deepcopy(internal_params.e_ops_mc),
+            c_ops = deepcopy(internal_params.c_ops),
             expvals = similar(internal_params.expvals),
             cache_mc = similar(internal_params.cache_mc),
             weights_mc = similar(internal_params.weights_mc),
@@ -598,27 +590,24 @@ function _dsf_mcsolve_prob_func(prob, i, repeat)
             dsf_cache1 = similar(internal_params.dsf_cache1),
             dsf_cache2 = similar(internal_params.dsf_cache2),
             expv_cache = copy(internal_params.expv_cache),
-            dsf_displace_cache_full = OperatorSum(
-                copy(internal_params.dsf_displace_cache_full.coefficients),
-                internal_params.dsf_displace_cache_full.operators,
-            ),
+            dsf_displace_cache_full = deepcopy(internal_params.dsf_displace_cache_full), # This brutally copies also the MatrixOperators, and it is inefficient.
         ),
     )
 
-    return remake(prob, p = prm)
+    f = deepcopy(prob.f.f)
+
+    return remake(prob, f = f, p = prm)
 end
 
 function dsf_mcsolveEnsembleProblem(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T},KetQuantumObject},
-    t_l::AbstractVector,
+    ψ0::QuantumObject{<:AbstractVector{T},KetQuantumObject},
+    tlist::AbstractVector,
     c_ops::Function,
-    op_list::Vector{TOl},
+    op_list::Union{AbstractVector,Tuple},
     α0_l::Vector{<:Number} = zeros(length(op_list)),
     dsf_params::NamedTuple = NamedTuple();
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (op_list, p) -> Vector{TOl}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (op_list, p) -> (),
     params::NamedTuple = NamedTuple(),
     ntraj::Int = 1,
     ensemble_method = EnsembleThreads(),
@@ -627,18 +616,18 @@ function dsf_mcsolveEnsembleProblem(
     krylov_dim::Int = min(5, cld(length(ψ0.data), 3)),
     progress_bar::Union{Bool,Val} = Val(true),
     kwargs...,
-) where {T,TOl,TJC<:LindbladJumpCallbackType}
-    op_l = op_list
-    H₀ = H(op_l .+ α0_l, dsf_params)
-    c_ops₀ = c_ops(op_l .+ α0_l, dsf_params)
-    e_ops₀ = e_ops(op_l .+ α0_l, dsf_params)
+) where {T,TJC<:LindbladJumpCallbackType}
+    op_list = deepcopy(op_list)
+    H₀ = H(op_list .+ α0_l, dsf_params)
+    c_ops₀ = c_ops(op_list .+ α0_l, dsf_params)
+    e_ops₀ = e_ops(op_list .+ α0_l, dsf_params)
 
     αt_list = convert(Vector{T}, α0_l)
     expv_cache = arnoldi(H₀.data, ψ0.data, krylov_dim)
 
-    dsf_displace_cache = map(op -> Qobj(op.data), op_l)
-    dsf_displace_cache_dag = map(op -> Qobj(sparse(op.data')), op_l)
-    dsf_displace_cache_full = OperatorSum(zeros(length(op_l) * 2), vcat(dsf_displace_cache, dsf_displace_cache_dag))
+    dsf_displace_cache = sum(op -> ScalarOperator(one(T)) * MatrixOperator(op.data), op_list)
+    dsf_displace_cache_dag = sum(op -> ScalarOperator(one(T)) * MatrixOperator(sparse(op.data')), op_list)
+    dsf_displace_cache_full = dsf_displace_cache + dsf_displace_cache_dag
 
     params2 = merge(
         params,
@@ -646,7 +635,7 @@ function dsf_mcsolveEnsembleProblem(
             H_fun = H,
             c_ops_fun = c_ops,
             e_ops_fun = e_ops,
-            op_l = op_l,
+            op_list = op_list,
             αt_list = αt_list,
             δα_list = δα_list,
             dsf_cache1 = similar(ψ0.data),
@@ -666,11 +655,9 @@ function dsf_mcsolveEnsembleProblem(
     return mcsolveEnsembleProblem(
         H₀,
         ψ0,
-        t_l,
+        tlist,
         c_ops₀;
         e_ops = e_ops₀,
-        alg = alg,
-        H_t = H_t,
         params = params2,
         ntraj = ntraj,
         ensemble_method = ensemble_method,
@@ -684,13 +671,12 @@ end
 @doc raw"""
     dsf_mcsolve(H::Function,
         ψ0::QuantumObject,
-        t_l::AbstractVector, c_ops::Function,
+        tlist::AbstractVector, c_ops::Function,
         op_list::Vector{TOl},
         α0_l::Vector{<:Number}=zeros(length(op_list)),
         dsf_params::NamedTuple=NamedTuple();
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
         e_ops::Function=(op_list,p) -> Vector{TOl}([]),
-        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         δα_list::Vector{<:Real}=fill(0.2, length(op_list)),
         ntraj::Int=1,
@@ -708,15 +694,14 @@ Time evolution of a quantum system using the Monte Carlo wave function method an
 """
 function dsf_mcsolve(
     H::Function,
-    ψ0::QuantumObject{<:AbstractArray{T},KetQuantumObject},
-    t_l::AbstractVector,
+    ψ0::QuantumObject{<:AbstractVector{T},KetQuantumObject},
+    tlist::AbstractVector,
     c_ops::Function,
-    op_list::Vector{TOl},
+    op_list::Union{AbstractVector,Tuple},
     α0_l::Vector{<:Number} = zeros(length(op_list)),
     dsf_params::NamedTuple = NamedTuple();
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Function = (op_list, p) -> Vector{TOl}([]),
-    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    e_ops::Function = (op_list, p) -> (),
     params::NamedTuple = NamedTuple(),
     δα_list::Vector{<:Real} = fill(0.2, length(op_list)),
     ntraj::Int = 1,
@@ -725,18 +710,17 @@ function dsf_mcsolve(
     krylov_dim::Int = min(5, cld(length(ψ0.data), 3)),
     progress_bar::Union{Bool,Val} = Val(true),
     kwargs...,
-) where {T,TOl,TJC<:LindbladJumpCallbackType}
+) where {T,TJC<:LindbladJumpCallbackType}
     ens_prob_mc = dsf_mcsolveEnsembleProblem(
         H,
         ψ0,
-        t_l,
+        tlist,
         c_ops,
         op_list,
         α0_l,
         dsf_params;
         alg = alg,
         e_ops = e_ops,
-        H_t = H_t,
         params = params,
         ntraj = ntraj,
         ensemble_method = ensemble_method,
diff --git a/src/utilities.jl b/src/utilities.jl
index 67f62697..00cf0866 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -141,6 +141,7 @@ getVal(x) = x # getVal for any other type
 
 _get_size(A::AbstractMatrix) = size(A)
 _get_size(A::AbstractVector) = (length(A), 1)
+_get_size(A::AbstractSciMLOperator) = size(A)
 
 _non_static_array_warning(argname, arg::Tuple{}) =
     throw(ArgumentError("The argument $argname must be a Tuple or a StaticVector of non-zero length."))
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index 8ef46b78..b472e0fe 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -120,6 +120,24 @@
         @test_throws DimensionMismatch Qobj(ρ_bra.data, type = OperatorBra, dims = 4)
     end
 
+    @testset "Checks on non-QuantumObjects" begin
+        x = 1
+        @test isket(x) == false
+        @test isbra(x) == false
+        @test isoper(x) == false
+        @test issuper(x) == false
+        @test isoperket(x) == false
+        @test isoperbra(x) == false
+
+        x = rand(ComplexF64, 2)
+        @test isket(x) == false
+        @test isbra(x) == false
+        @test isoper(x) == false
+        @test issuper(x) == false
+        @test isoperket(x) == false
+        @test isoperbra(x) == false
+    end
+
     @testset "arithmetic" begin
         a = sprand(ComplexF64, 100, 100, 0.1)
         a2 = Qobj(a)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
new file mode 100644
index 00000000..933590e6
--- /dev/null
+++ b/test/core-test/quantum_objects_evo.jl
@@ -0,0 +1,180 @@
+@testset "Quantum Objects Evolution" verbose = true begin
+    # DomainError: incompatible between size of array and type
+    @testset "Thrown Errors" begin
+        a = MatrixOperator(rand(ComplexF64, 3, 2))
+        for t in [Operator, SuperOperator]
+            @test_throws DomainError QobjEvo(a, type = t)
+        end
+
+        a = MatrixOperator(rand(ComplexF64, 3, 2))
+        for t in (Ket, Bra, OperatorKet, OperatorBra)
+            @test_throws ArgumentError QobjEvo(a, type = t)
+        end
+
+        a = QobjEvo(destroy(20))
+        @test_throws ArgumentError QobjEvo(a, type = SuperOperator)
+
+        a = MatrixOperator(rand(ComplexF64, 5, 5))
+        @test_throws DimensionMismatch QobjEvo(a, type = SuperOperator)
+
+        ψ = fock(10, 3)
+        @test_throws TypeError QobjEvo(ψ)
+    end
+
+    # unsupported type of dims
+    @testset "unsupported dims" begin
+        a = MatrixOperator(rand(2, 2))
+        @test_throws ArgumentError QobjEvo(a, dims = 2.0)
+        @test_throws ArgumentError QobjEvo(a, dims = 2.0 + 0.0im)
+        @test_throws DomainError QobjEvo(a, dims = 0)
+        @test_throws DomainError QobjEvo(a, dims = (2, -2))
+        @test_logs (
+            :warn,
+            "The argument dims should be a Tuple or a StaticVector for better performance. Try to use `dims = (2, 2)` or `dims = SVector(2, 2)` instead of `dims = [2, 2]`.",
+        ) QobjEvo(MatrixOperator(rand(4, 4)), dims = [2, 2])
+    end
+
+    @testset "Operator and SuperOperator" begin
+        a = MatrixOperator(sprand(ComplexF64, 100, 100, 0.1))
+        a2 = QobjEvo(a)
+        a3 = QobjEvo(a, type = SuperOperator)
+
+        @test isket(a2) == false
+        @test isbra(a2) == false
+        @test isoper(a2) == true
+        @test issuper(a2) == false
+        @test isoperket(a2) == false
+        @test isoperbra(a2) == false
+        @test isket(a3) == false
+        @test isbra(a3) == false
+        @test isoper(a3) == false
+        @test issuper(a3) == true
+        @test isoperket(a3) == false
+        @test isoperbra(a3) == false
+        @test_throws DimensionMismatch QobjEvo(a, dims = 2)
+    end
+
+    @testset "Promote Operators Type" begin
+        a = destroy(20)
+        A = QobjEvo(a)
+        @test QuantumToolbox.promote_op_type(a, A) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(A, a) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(A, A) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(a, a) == QuantumObject
+    end
+
+    @testset "arithmetic" begin
+        a = MatrixOperator(sprand(ComplexF64, 100, 100, 0.1))
+        a2 = QobjEvo(a)
+        a3 = QobjEvo(a, type = SuperOperator)
+
+        @test +a2 == a2
+        @test -(-a2) == a2
+        @test a2 + 2 == 2 + a2
+        @test (a2 + 2).data == a2.data + 2 * I
+        @test a2 * 2 == 2 * a2
+
+        @test trans(trans(a2)) == a2
+        @test trans(a2).data == transpose(a2.data)
+        # @test adjoint(a2) ≈ trans(conj(a2)) # Currently doesn't work
+        @test adjoint(adjoint(a2)) == a2
+        @test adjoint(a2).data == adjoint(a2.data)
+
+        N = 10
+        # We use MatrixOperator instead of directly using a Qobj to increase coverage
+        a = QobjEvo(MatrixOperator(sprand(ComplexF64, N, N, 5 / N)), Operator, N)
+        a_d = a'
+        X = a + a_d
+        # Y = 1im * (a - a_d) # Currently doesn't work. Fix in SciMLOperators.jl
+        Z = a + trans(a)
+        @test isherm(X) == true
+        # @test isherm(Y) == true
+        # @test issymmetric(Y) == false
+        @test issymmetric(Z) == true
+    end
+
+    # TODO: Implement a new show method for QuantumObjectEvolution
+    # @testset "REPL show" begin
+    #     N = 10
+    #     a = QobjEvo(destroy(N))
+
+    #     opstring = sprint((t, s) -> show(t, "text/plain", s), a)
+    #     datastring = sprint((t, s) -> show(t, "text/plain", s), a.data)
+    #     a_dims = a.dims
+    #     a_size = size(a)
+    #     a_isherm = isherm(a)
+    #     @test opstring ==
+    #           "Quantum Object:   type=Operator   dims=$a_dims   size=$a_size   ishermitian=$a_isherm\n$datastring"
+
+    #     a = spre(a)
+    #     opstring = sprint((t, s) -> show(t, "text/plain", s), a)
+    #     datastring = sprint((t, s) -> show(t, "text/plain", s), a.data)
+    #     a_dims = a.dims
+    #     a_size = size(a)
+    #     a_isherm = isherm(a)
+    #     @test opstring == "Quantum Object:   type=SuperOperator   dims=$a_dims   size=$a_size\n$datastring"
+    # end
+
+    @testset "Type Inference (QuantumObject)" begin
+        for T in [ComplexF32, ComplexF64]
+            N = 4
+            a = MatrixOperator(rand(T, N, N))
+            @inferred QobjEvo(a)
+            for type in [Operator, SuperOperator]
+                @inferred QobjEvo(a, type = type)
+            end
+        end
+
+        @testset "Math Operation" begin
+            a = QobjEvo(destroy(20))
+            σx = QobjEvo(sigmax())
+            @inferred a + a
+            @inferred a + a'
+            # @inferred a + 2 # TODO fix in SciMLOperators.jl
+            @inferred 2 * a
+            @inferred a / 2
+            @inferred a * a
+            @inferred a * a'
+
+            # TODO: kron is currently not supported
+            # @inferred kron(a)
+            # @inferred kron(a, σx)
+            # @inferred kron(a, eye(2))
+        end
+    end
+
+    # TODO: tensor is currently not supported
+    # @testset "tensor" begin
+    #     σx = sigmax()
+    #     X3 = kron(σx, σx, σx)
+    #     @test tensor(σx) == kron(σx)
+    #     @test tensor(fill(σx, 3)...) == X3
+    #     X_warn = @test_logs (
+    #         :warn,
+    #         "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead.",
+    #     ) tensor(fill(σx, 3))
+    #     @test X_warn == X3
+    # end
+
+    @testset "Time Dependent Operators" begin
+        N = 10
+        a = destroy(N)
+        coef1(p, t) = exp(-1im * p.ω1 * t)
+        coef2(p, t) = sin(p.ω2 * t)
+
+        @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
+    end
+end
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index 1b3df4e7..507f1227 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -64,8 +64,11 @@
     e_ops = [a_d * a]
     psi0 = fock(N, 3)
     t_l = LinRange(0, 100 * 2π, 1000)
-    H_t_f = TimeDependentOperatorSum((((t, p) -> sin(t)),), (H_t,)) # It will be converted to liouvillian internally
-    sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, H_t = H_t_f, progress_bar = Val(false))
+
+    coeff(p, t) = sin(t)
+    H_td = (H, (H_t, coeff))
+
+    sol_me = mesolve(H_td, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
     ρ_ss1 = steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetLinearSystem())[1]
     ρ_ss2 =
         steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian())
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 50d3306e..666aa38a 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -59,18 +59,44 @@
         sol_me2 = mesolve(H, psi0, t_l, c_ops, progress_bar = Val(false))
         sol_me3 = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
         sol_mc = mcsolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
+        sol_mc2 = mcsolve(
+            H,
+            psi0,
+            t_l,
+            c_ops,
+            ntraj = 500,
+            e_ops = e_ops,
+            progress_bar = Val(false),
+            jump_callback = DiscreteLindbladJumpCallback(),
+        )
         sol_mc_states = mcsolve(H, psi0, t_l, c_ops, ntraj = 500, saveat = t_l, progress_bar = Val(false))
+        sol_mc_states2 = mcsolve(
+            H,
+            psi0,
+            t_l,
+            c_ops,
+            ntraj = 500,
+            saveat = t_l,
+            progress_bar = Val(false),
+            jump_callback = DiscreteLindbladJumpCallback(),
+        )
         sol_sse = ssesolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
 
         ρt_mc = [ket2dm.(normalize.(states)) for states in sol_mc_states.states]
         expect_mc_states = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc)
         expect_mc_states_mean = sum(expect_mc_states, dims = 2) / size(expect_mc_states, 2)
 
+        ρt_mc2 = [ket2dm.(normalize.(states)) for states in sol_mc_states2.states]
+        expect_mc_states2 = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc2)
+        expect_mc_states_mean2 = sum(expect_mc_states2, dims = 2) / size(expect_mc_states2, 2)
+
         sol_me_string = sprint((t, s) -> show(t, "text/plain", s), sol_me)
         sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
         sol_sse_string = sprint((t, s) -> show(t, "text/plain", s), sol_sse)
         @test sum(abs.(sol_mc.expect .- sol_me.expect)) / length(t_l) < 0.1
+        @test sum(abs.(sol_mc2.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test sum(abs.(vec(expect_mc_states_mean) .- vec(sol_me.expect))) / length(t_l) < 0.1
+        @test sum(abs.(vec(expect_mc_states_mean2) .- vec(sol_me.expect))) / length(t_l) < 0.1
         @test sum(abs.(sol_sse.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test length(sol_me.times) == length(t_l)
         @test length(sol_me.states) == 1
@@ -117,37 +143,219 @@
               "abstol = $(sol_sse.abstol)\n" *
               "reltol = $(sol_sse.reltol)\n"
 
+        # Time-Dependent Hamiltonian
+        # ssesolve is slow to be run on CI. It is not removed from the test because it may be useful for testing in more powerful machines.
+
+        N = 10
+        a = tensor(destroy(N), qeye(2))
+        σm = tensor(qeye(N), sigmam())
+        σz = tensor(qeye(N), sigmaz())
+        ω = 1.0
+        ωd = 1.02
+        Δ = ω - ωd
+        F = 0.05
+        g = 0.1
+        γ = 0.1
+        nth = 0.001
+
+        # Time Evolution in the drive frame
+
+        H = Δ * a' * a + Δ * σz / 2 + g * (a' * σm + a * σm') + F * (a + a')
+        c_ops = [sqrt(γ * (1 + nth)) * a, sqrt(γ * nth) * a', sqrt(γ * (1 + nth)) * σm, sqrt(γ * nth) * σm']
+        e_ops = [a' * a, σz]
+
+        ψ0 = tensor(basis(N, 0), basis(2, 1))
+        tlist = range(0, 2 / γ, 1000)
+
+        rng = MersenneTwister(12)
+
+        sol_se = sesolve(H, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false))
+        sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
+        sol_mc = mcsolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+        # sol_sse = ssesolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+
+        # Time Evolution in the lab frame
+
+        H = ω * a' * a + ω * σz / 2 + g * (a' * σm + a * σm')
+
+        coef1(p, t) = p.F * exp(1im * p.ωd * t)
+        coef2(p, t) = p.F * exp(-1im * p.ωd * t)
+
+        H_td = (H, (a, coef1), (a', coef2))
+        p = (F = F, ωd = ωd)
+
+        sol_se_td = sesolve(H_td, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p)
+        sol_me_td = mesolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
+        sol_mc_td = mcsolve(
+            H_td,
+            ψ0,
+            tlist,
+            c_ops,
+            ntraj = 500,
+            e_ops = e_ops,
+            progress_bar = Val(false),
+            params = p,
+            rng = rng,
+        )
+        # sol_sse_td = ssesolve(H_td, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng)
+
+        @test sol_se.expect ≈ sol_se_td.expect atol = 1e-6 * length(tlist)
+        @test sol_me.expect ≈ sol_me_td.expect atol = 1e-6 * length(tlist)
+        @test sol_mc.expect ≈ sol_mc_td.expect atol = 1e-2 * length(tlist)
+        # @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)
+
+        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,
+            ψ0,
+            tlist,
+            c_ops,
+            ntraj = 500,
+            e_ops = e_ops,
+            progress_bar = Val(false),
+            params = p,
+            rng = rng,
+        )
+        # sol_sse_td2 =
+        # ssesolve(H_td2, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng)
+
+        @test sol_se.expect ≈ sol_se_td2.expect atol = 1e-6 * length(tlist)
+        @test sol_me.expect ≈ sol_me_td2.expect atol = 1e-6 * length(tlist)
+        @test sol_mc.expect ≈ sol_mc_td2.expect atol = 1e-2 * length(tlist)
+        # @test sol_sse.expect ≈ sol_sse_td2.expect atol = 1e-2 * length(tlist)
+
         @testset "Type Inference mesolve" begin
-            @inferred mesolveProblem(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mesolveProblem(H, psi0, [0, 10], c_ops, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mesolveProblem(H, Qobj(zeros(Int64, N)), t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mesolve(H, psi0, t_l, c_ops, progress_bar = Val(false))
-            @inferred mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
-            @inferred mesolve(H, psi0, t_l, (a, a'), e_ops = (a_d * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple of different types
+            @inferred mesolveProblem(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
+            @inferred mesolveProblem(H, ψ0, [0, 10], c_ops, e_ops = e_ops, progress_bar = Val(false))
+            @inferred mesolveProblem(
+                H,
+                tensor(Qobj(zeros(Int64, N)), Qobj([0, 1])),
+                tlist,
+                c_ops,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+            )
+            @inferred mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
+            @inferred mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
+            @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(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
         end
 
         @testset "Type Inference mcsolve" begin
-            @inferred mcsolveEnsembleProblem(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mcsolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
-            @inferred mcsolve(H, psi0, t_l, c_ops, ntraj = 500, progress_bar = Val(true))
-            @inferred mcsolve(H, psi0, [0, 10], c_ops, ntraj = 500, progress_bar = Val(false))
-            @inferred mcsolve(H, Qobj(zeros(Int64, N)), t_l, c_ops, ntraj = 500, progress_bar = Val(false))
-            @inferred mcsolve(H, psi0, t_l, (a, a'), e_ops = (a_d * a, a'), ntraj = 500, progress_bar = Val(false)) # We test the type inference for Tuple of different types
+            @inferred mcsolveEnsembleProblem(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = 5,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+                rng = rng,
+            )
+            @inferred mcsolve(H, ψ0, tlist, c_ops, ntraj = 5, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+            @inferred mcsolve(H, ψ0, tlist, c_ops, ntraj = 5, progress_bar = Val(true), rng = rng)
+            @inferred mcsolve(H, ψ0, [0, 10], c_ops, ntraj = 5, progress_bar = Val(false), rng = rng)
+            @inferred mcsolve(
+                H,
+                tensor(Qobj(zeros(Int64, N)), Qobj([0, 1])),
+                tlist,
+                c_ops,
+                ntraj = 5,
+                progress_bar = Val(false),
+                rng = rng,
+            )
+            @inferred mcsolve(
+                H,
+                ψ0,
+                tlist,
+                (a, a'),
+                e_ops = (a' * a, a'),
+                ntraj = 5,
+                progress_bar = Val(false),
+                rng = rng,
+            ) # We test the type inference for Tuple of different types
+            @inferred mcsolve(
+                H_td,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = 5,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+                params = p,
+                rng = rng,
+            )
         end
 
         @testset "Type Inference ssesolve" begin
+            c_ops_tuple = Tuple(c_ops) # To avoid type instability, we must have a Tuple instead of a Vector
             @inferred ssesolveEnsembleProblem(
                 H,
-                psi0,
-                t_l,
-                c_ops,
-                ntraj = 500,
+                ψ0,
+                tlist,
+                c_ops_tuple,
+                ntraj = 5,
                 e_ops = e_ops,
                 progress_bar = Val(false),
+                rng = rng,
+            )
+            @inferred ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_tuple,
+                ntraj = 5,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+                rng = rng,
+            )
+            @inferred ssesolve(H, ψ0, tlist, c_ops_tuple, ntraj = 5, progress_bar = Val(true), rng = rng)
+            @inferred ssesolve(H, ψ0, [0, 10], c_ops_tuple, ntraj = 5, progress_bar = Val(false), rng = rng)
+            @inferred ssesolve(
+                H,
+                tensor(Qobj(zeros(Int64, N)), Qobj([0, 1])),
+                tlist,
+                c_ops_tuple,
+                ntraj = 5,
+                progress_bar = Val(false),
+                rng = rng,
+            )
+            @inferred ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_tuple,
+                ntraj = 5,
+                e_ops = (a' * a, a'),
+                progress_bar = Val(false),
+                rng = rng,
+            ) # We test the type inference for Tuple of different types
+            @inferred ssesolve(
+                H_td,
+                ψ0,
+                tlist,
+                c_ops_tuple,
+                ntraj = 5,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+                params = p,
+                rng = rng,
             )
-            @inferred ssesolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
-            @inferred ssesolve(H, psi0, t_l, c_ops, ntraj = 500, progress_bar = Val(true))
         end
 
         @testset "mcsolve and ssesolve reproducibility" begin
@@ -170,17 +378,14 @@
             tlist = range(0, 20 / γ, 1000)
 
             rng = MersenneTwister(1234)
-            sleep(0.1) # If we don't sleep, we get an error (why?)
             sol_mc1 = mcsolve(H, psi0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng)
             sol_sse1 = ssesolve(H, psi0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
 
             rng = MersenneTwister(1234)
-            sleep(0.1)
             sol_mc2 = mcsolve(H, psi0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng)
             sol_sse2 = ssesolve(H, psi0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
 
             rng = MersenneTwister(1234)
-            sleep(0.1)
             sol_mc3 = mcsolve(H, psi0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = rng)
 
             @test sol_mc1.expect ≈ sol_mc2.expect atol = 1e-10
@@ -199,15 +404,16 @@
         N = 10
         a = destroy(N)
         H = a' * a
+        c_ops = [sqrt(0.1) * a]
         psi0 = basis(N, 3)
         t_l = LinRange(0, 100, 1000)
         psi_wrong = basis(N - 1, 3)
         @test_throws DimensionMismatch sesolve(H, psi_wrong, t_l)
-        @test_throws DimensionMismatch mesolve(H, psi_wrong, t_l)
-        @test_throws DimensionMismatch mcsolve(H, psi_wrong, t_l)
+        @test_throws DimensionMismatch mesolve(H, psi_wrong, t_l, c_ops)
+        @test_throws DimensionMismatch mcsolve(H, psi_wrong, t_l, c_ops)
         @test_throws ArgumentError sesolve(H, psi0, t_l, save_idxs = [1, 2])
-        @test_throws ArgumentError mesolve(H, psi0, t_l, save_idxs = [1, 2])
-        @test_throws ArgumentError mcsolve(H, psi0, t_l, save_idxs = [1, 2])
+        @test_throws ArgumentError mesolve(H, psi0, t_l, c_ops, save_idxs = [1, 2])
+        @test_throws ArgumentError mcsolve(H, psi0, t_l, c_ops, save_idxs = [1, 2])
     end
 
     @testset "example" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 566104eb..57407fb0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,7 @@ using Pkg
 using QuantumToolbox
 using QuantumToolbox: position, momentum
 using Random
+using SciMLOperators
 
 const GROUP = get(ENV, "GROUP", "All")
 
@@ -21,6 +22,7 @@ core_tests = [
     "permutation.jl",
     "progress_bar.jl",
     "quantum_objects.jl",
+    "quantum_objects_evo.jl",
     "states_and_operators.jl",
     "steady_state.jl",
     "time_evolution.jl",