diff --git a/docs/src/api.md b/docs/src/api.md
index 554d1e8b0..e81e0b4e8 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -28,9 +28,10 @@ SuperOperatorQuantumObject
 SuperOperator
 QuantumObject
 QuantumObjectEvolution
-size
-eltype
-length
+Base.size
+Base.eltype
+Base.length
+SciMLOperators.cache_operator
 ```
 
 ## [Qobj boolean functions](@id doc-API:Qobj-boolean-functions)
@@ -46,7 +47,8 @@ LinearAlgebra.ishermitian
 LinearAlgebra.issymmetric
 LinearAlgebra.isposdef
 isunitary
-isconstant
+SciMLOperators.iscached
+SciMLOperators.isconstant
 ```
 
 ## [Qobj arithmetic and attributes](@id doc-API:Qobj-arithmetic-and-attributes)
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 098cfb376..4077ed050 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -3,17 +3,19 @@ module QuantumToolbox
 # Re-export:
 #   1. StaticArraysCore.SVector for the type of dims
 #   2. basic functions in LinearAlgebra and SparseArrays
+#   3. some functions in SciMLOperators
 import Reexport: @reexport
 @reexport import StaticArraysCore: SVector
 @reexport using LinearAlgebra
 @reexport using SparseArrays
+@reexport import SciMLOperators: cache_operator, iscached, isconstant
 
 # other functions in LinearAlgebra
 import LinearAlgebra: BlasReal, BlasInt, BlasFloat, BlasComplex, checksquare
 import LinearAlgebra.BLAS: @blasfunc
 import LinearAlgebra.LAPACK: hseqr!
 
-# SciML packages (for OrdinaryDiffEq and LinearSolve)
+# SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
 import SciMLBase:
     solve,
     solve!,
@@ -34,16 +36,15 @@ import SciMLBase:
     DiscreteCallback
 import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1
 import SciMLOperators:
+    SciMLOperators,
     AbstractSciMLOperator,
     MatrixOperator,
     ScalarOperator,
     ScaledOperator,
     AddedOperator,
     IdentityOperator,
-    cache_operator,
     update_coefficients!,
-    concretize,
-    isconstant
+    concretize
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState
diff --git a/src/qobj/boolean_functions.jl b/src/qobj/boolean_functions.jl
index cccea5321..5262d7e90 100644
--- a/src/qobj/boolean_functions.jl
+++ b/src/qobj/boolean_functions.jl
@@ -3,7 +3,7 @@ All boolean functions for checking the data or type in `QuantumObject`
 =#
 
 export isket, isbra, isoper, isoperbra, isoperket, issuper
-export isunitary, isconstant
+export isunitary
 
 @doc raw"""
     isbra(A)
@@ -91,8 +91,15 @@ 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)
+    SciMLOperators.iscached(A::AbstractQuantumObject)
+
+Test whether the [`AbstractQuantumObject`](@ref) `A` has preallocated caches for inplace evaluations.
+"""
+SciMLOperators.iscached(A::AbstractQuantumObject) = iscached(A.data)
+
+@doc raw"""
+    SciMLOperators.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)
+SciMLOperators.isconstant(A::AbstractQuantumObject) = isconstant(A.data)
diff --git a/src/qobj/quantum_object.jl b/src/qobj/quantum_object.jl
index 3c0121dc3..d31f95ac3 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -3,6 +3,7 @@ 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
+    - SciMLOperators: cache_operator
 =#
 
 export QuantumObject
@@ -167,6 +168,41 @@ SparseArrays.droptol!(A::QuantumObject{<:AbstractSparseArray}, tol::Real) = (dro
 SparseArrays.dropzeros(A::QuantumObject{<:AbstractSparseArray}) = QuantumObject(dropzeros(A.data), A.type, A.dims)
 SparseArrays.dropzeros!(A::QuantumObject{<:AbstractSparseArray}) = (dropzeros!(A.data); return A)
 
+@doc raw"""
+    SciMLOperators.cached_operator(L::AbstractQuantumObject, u)
+
+Allocate caches for [`AbstractQuantumObject`](@ref) `L` for in-place evaluation with `u`-like input vectors.
+
+Here, `u` can be in either the following types:
+- `AbstractVector`
+- [`Ket`](@ref)-type [`QuantumObject`](@ref) (if `L` is an [`Operator`](@ref))
+- [`OperatorKet`](@ref)-type [`QuantumObject`](@ref) (if `L` is a [`SuperOperator`](@ref))
+"""
+SciMLOperators.cache_operator(
+    L::AbstractQuantumObject{DT,OpType},
+    u::AbstractVector,
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    get_typename_wrapper(L)(cache_operator(L.data, sparse_to_dense(similar(u))), L.type, L.dims)
+
+function SciMLOperators.cache_operator(
+    L::AbstractQuantumObject{DT1,OpType},
+    u::QuantumObject{DT2,SType},
+) where {
+    DT1,
+    DT2,
+    OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
+    SType<:Union{KetQuantumObject,OperatorKetQuantumObject},
+}
+    check_dims(L, u)
+
+    if isoper(L) && isoperket(u)
+        throw(ArgumentError("The input state `u` must be a Ket if `L` is an Operator."))
+    elseif issuper(L) && isket(u)
+        throw(ArgumentError("The input state `u` must be an OperatorKet if `L` is a SuperOperator."))
+    end
+    return cache_operator(L, u.data)
+end
+
 # data type conversions
 Base.Vector(A::QuantumObject{<:AbstractVector}) = QuantumObject(Vector(A.data), A.type, A.dims)
 Base.Vector{T}(A::QuantumObject{<:AbstractVector}) where {T<:Number} = QuantumObject(Vector{T}(A.data), A.type, A.dims)
diff --git a/src/qobj/quantum_object_base.jl b/src/qobj/quantum_object_base.jl
index 0190e4c44..ea4f9df29 100644
--- a/src/qobj/quantum_object_base.jl
+++ b/src/qobj/quantum_object_base.jl
@@ -209,8 +209,6 @@ function _check_QuantumObject(type::OperatorBraQuantumObject, dims, m::Int, n::I
     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
index 8786669b5..2b698cfc1 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -158,9 +158,8 @@ 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)
+    op, data = _QobjEvo_generate_data(op_func_list, α)
     dims = op.dims
     if type isa Nothing
         type = op.type
@@ -177,24 +176,21 @@ 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)
+    return QuantumObjectEvolution(_make_SciMLOperator(op, α), 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
+    elseif type != op.type
         throw(
             ArgumentError(
                 "The type of the QuantumObjectEvolution object cannot be changed when using another QuantumObjectEvolution object as input.",
@@ -217,7 +213,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
 - `α`: 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)
+@generated function _QobjEvo_generate_data(op_func_list::Tuple, α)
     op_func_list_types = op_func_list.parameters
     N = length(op_func_list_types)
 
@@ -242,9 +238,9 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
             data_type = op_type.parameters[1]
             dims_expr = (dims_expr..., :($op.dims))
             if i == 1
-                first_op = :(f($op))
+                first_op = :($op)
             end
-            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α, f = f))
+            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α))
         else
             op_type = op_func_type
             (isoper(op_type) || issuper(op_type)) ||
@@ -255,7 +251,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
             if i == 1
                 first_op = :(op_func_list[$i])
             end
-            qobj_expr_const = :($qobj_expr_const + f(op_func_list[$i]))
+            qobj_expr_const = :($qobj_expr_const + op_func_list[$i])
         end
     end
 
@@ -272,20 +268,20 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
     end
 end
 
-function _make_SciMLOperator(op_func::Tuple, α; f::Function = identity)
+function _make_SciMLOperator(op_func::Tuple, α)
     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)
+        return ScalarOperator(zero(T), update_func) * MatrixOperator(op_func[1].data)
     end
-    return ScalarOperator(zero(T), update_func) * MatrixOperator(α * f(op_func[1]).data)
+    return ScalarOperator(zero(T), update_func) * MatrixOperator(α * op_func[1].data)
 end
 
-function _make_SciMLOperator(op::QuantumObject, α; f::Function = identity)
+function _make_SciMLOperator(op::QuantumObject, α)
     if α isa Nothing
-        return MatrixOperator(f(op).data)
+        return MatrixOperator(op.data)
     end
-    return MatrixOperator(α * f(op.data))
+    return MatrixOperator(α * op.data)
 end
 
 @doc raw"""
@@ -375,15 +371,11 @@ function (A::QuantumObjectEvolution)(
     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.",
-            ),
-        )
+        throw(ArgumentError("The input state must be a Ket 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.",
+                "The input state must be an OperatorKet if the QuantumObjectEvolution object is a SuperOperator.",
             ),
         )
     end
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 2408683f2..fbb51b8f3 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -5,7 +5,7 @@ Functions for generating (common) quantum super-operators.
 export spre, spost, sprepost, liouvillian, lindblad_dissipator
 
 # intrinsic functions for super-operators
-# (keep these because they take AbstractMatrix as input and ensure the output is sparse matrix)
+## keep these because they take AbstractMatrix as input and ensure the output is sparse matrix
 _spre(A::AbstractMatrix, Id::AbstractMatrix) = kron(Id, sparse(A))
 _spre(A::AbstractSparseMatrix, Id::AbstractMatrix) = kron(Id, A)
 _spost(B::AbstractMatrix, Id::AbstractMatrix) = kron(transpose(sparse(B)), Id)
@@ -14,19 +14,51 @@ _sprepost(A::AbstractMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), spa
 _sprepost(A::AbstractMatrix, B::AbstractSparseMatrix) = kron(transpose(B), sparse(A))
 _sprepost(A::AbstractSparseMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), A)
 _sprepost(A::AbstractSparseMatrix, B::AbstractSparseMatrix) = kron(transpose(B), A)
-_liouvillian(H::AbstractMatrix, Id::AbstractMatrix) = -1im * (_spre(H, Id) - _spost(H, Id))
+function _sprepost(A, B) # for any other input types
+    Id_cache = I(size(A, 1))
+    return _spre(A, Id_cache) * _spost(B, Id_cache)
+end
 
-# (if input is AbstractSciMLOperator)
+## if input is AbstractSciMLOperator 
+## some of them are optimzed to speed things up
+## the rest of the SciMLOperators will just use lazy tensor (and prompt a warning)
 _spre(A::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spre(A.A, Id))
 _spre(A::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(A.λ, _spre(A.L, Id))
 _spre(A::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spre(op, Id), +, A.ops)
+function _spre(A::AbstractSciMLOperator, Id::AbstractMatrix)
+    _lazy_tensor_warning("spre", A)
+    return kron(Id, A)
+end
+
 _spost(B::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spost(B.A, Id))
 _spost(B::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(B.λ, _spost(B.L, Id))
 _spost(B::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spost(op, Id), +, B.ops)
+function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
+    _lazy_tensor_warning("spost", B)
+    return kron(transpose(B), Id)
+end
+
+## intrinsic liouvillian 
+_liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
+    -1im * (_spre(H, Id) - _spost(H, Id))
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
 _liouvillian(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian(H.L, Id))
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian(op, Id), +, H.ops)
-# TODO: support `_sprepost`, `sprepost`, and `lindblad_dissipator` for AbstractSciMLOperator (allow c_ops with Vector{QobjEvo})
+
+# intrinsic lindblad_dissipator
+function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
+    Od_O = O' * O
+    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
+end
+function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
+    _O = O.A
+    Od_O = _O' * _O
+    return MatrixOperator(_sprepost(_O, _O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2)
+end
+function _lindblad_dissipator(O::ScaledOperator, Id::AbstractMatrix)
+    λc_λ = conj(O.λ) * O.λ
+    return ScaledOperator(λc_λ, _lindblad_dissipator(O.L, Id))
+end
 
 @doc raw"""
     spre(A::AbstractQuantumObject, Id_cache=I(size(A,1)))
@@ -41,6 +73,8 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 (see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
 The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
+
+See also [`spost`](@ref) and [`sprepost`](@ref).
 """
 spre(A::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(A, 1))) where {DT} =
     get_typename_wrapper(A)(_spre(A.data, Id_cache), SuperOperator, A.dims)
@@ -58,12 +92,14 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 (see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
 The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.
+
+See also [`spre`](@ref) and [`sprepost`](@ref).
 """
 spost(B::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(B, 1))) where {DT} =
     get_typename_wrapper(B)(_spost(B.data, Id_cache), SuperOperator, B.dims)
 
 @doc raw"""
-    sprepost(A::QuantumObject, B::QuantumObject)
+    sprepost(A::AbstractQuantumObject, B::AbstractQuantumObject)
 
 Returns the [`SuperOperator`](@ref) form of `A` and `B` acting on the left and right of the density matrix operator, respectively: ``\mathcal{O} \left( \hat{A}, \hat{B} \right) \left[ \hat{\rho} \right] = \hat{A} \hat{\rho} \hat{B}``.
 
@@ -77,16 +113,15 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 See also [`spre`](@ref) and [`spost`](@ref).
 """
 function sprepost(
-    A::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-) where {T1,T2}
+    A::AbstractQuantumObject{DT1,OperatorQuantumObject},
+    B::AbstractQuantumObject{DT2,OperatorQuantumObject},
+) where {DT1,DT2}
     check_dims(A, B)
-
-    return QuantumObject(_sprepost(A.data, B.data), SuperOperator, A.dims)
+    return promote_op_type(A, B)(_sprepost(A.data, B.data), SuperOperator, A.dims)
 end
 
 @doc raw"""
-    lindblad_dissipator(O::QuantumObject, Id_cache=I(size(O,1))
+    lindblad_dissipator(O::AbstractQuantumObject, Id_cache=I(size(O,1))
 
 Returns the Lindblad [`SuperOperator`](@ref) defined as
 
@@ -99,14 +134,11 @@ The optional argument `Id_cache` can be used to pass a precomputed identity matr
 
 See also [`spre`](@ref), [`spost`](@ref), and [`sprepost`](@ref).
 """
-function lindblad_dissipator(O::QuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(O, 1))) where {DT}
-    Od_O = O' * O
-    return sprepost(O, O') - (spre(Od_O, Id_cache) + spost(Od_O, Id_cache)) / 2
-end
-# TODO: suppport collapse operator given as QobjEvo-type
+lindblad_dissipator(O::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(O, 1))) where {DT} =
+    get_typename_wrapper(O)(_lindblad_dissipator(O.data, Id_cache), SuperOperator, O.dims)
 
 # It is already a SuperOperator
-lindblad_dissipator(O::QuantumObject{DT,SuperOperatorQuantumObject}, Id_cache = nothing) where {DT} = O
+lindblad_dissipator(O::AbstractQuantumObject{DT,SuperOperatorQuantumObject}, Id_cache = nothing) where {DT} = O
 
 @doc raw"""
     liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dims)))
@@ -134,7 +166,17 @@ function liouvillian(
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
     if !(c_ops isa Nothing)
-        L += mapreduce(lindblad_dissipator, +, c_ops)
+        # sum all the Qobj first
+        c_ops_ti = filter(op -> isa(op, QuantumObject), c_ops)
+        if !isempty(c_ops_ti)
+            L += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_ti)
+        end
+
+        # sum rest of the QobjEvo together
+        c_ops_td = filter(op -> isa(op, QuantumObjectEvolution), c_ops)
+        if !isempty(c_ops_td)
+            L += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_td)
+        end
     end
     return L
 end
diff --git a/src/qobj/synonyms.jl b/src/qobj/synonyms.jl
index 8155ba3fa..8026ad7d2 100644
--- a/src/qobj/synonyms.jl
+++ b/src/qobj/synonyms.jl
@@ -18,7 +18,7 @@ Note that this functions is same as `QuantumObject(A; type=type, dims=dims)`.
 Qobj(A; kwargs...) = QuantumObject(A; kwargs...)
 
 @doc raw"""
-    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type::Union{Nothing, QuantumObjectType}=nothing, f::Function=identity)
+    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type::Union{Nothing, QuantumObjectType}=nothing)
 
 Generate [`QuantumObjectEvolution`](@ref).
 
@@ -27,7 +27,6 @@ Note that this functions is same as `QuantumObjectEvolution(op_func_list)`. If `
 # 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.
@@ -105,8 +104,7 @@ 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)
+) = QuantumObjectEvolution(op_func_list, α; type = type)
 
 QobjEvo(data::AbstractSciMLOperator, type::QuantumObjectType, dims::Integer) =
     QuantumObjectEvolution(data, type, SVector{1,Int}(dims))
diff --git a/src/utilities.jl b/src/utilities.jl
index 00cf0866c..150751e6a 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -128,6 +128,8 @@ function convert_unit(value::T, unit1::Symbol, unit2::Symbol) where {T<:Real}
     return _FType(T)(value * (_energy_units[unit1] / _energy_units[unit2]))
 end
 
+get_typename_wrapper(A) = Base.typename(typeof(A)).wrapper
+
 _get_dense_similar(A::AbstractArray, args...) = similar(A, args...)
 _get_dense_similar(A::AbstractSparseMatrix, args...) = similar(nonzeros(A), args...)
 
@@ -151,6 +153,9 @@ _non_static_array_warning(argname, arg::AbstractVector{T}) where {T} =
           join(arg, ", ") *
           ")` instead of `$argname = $arg`." maxlog = 1
 
+_lazy_tensor_warning(func_name::String, data::AbstractSciMLOperator) =
+    @warn "The function `$func_name` uses lazy tensor (which can hurt performance) for data type: $(get_typename_wrapper(data))"
+
 # functions for getting Float or Complex element type
 _FType(::AbstractArray{T}) where {T<:Number} = _FType(T)
 _FType(::Type{Int32}) = Float32
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index be24e3bfb..36e720b51 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -80,12 +80,17 @@
         @test issuper(a2) == false
         @test isoperket(a2) == false
         @test isoperbra(a2) == false
+        @test iscached(a2) == true
+        @test isconstant(a2) == true
+        @test isunitary(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 iscached(a3) == true
+        @test isconstant(a3) == true
         @test isunitary(a3) == false
         @test_throws DimensionMismatch Qobj(a, dims = 2)
     end
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 4cad02526..bf9f6116d 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -165,14 +165,18 @@
         a = destroy(N)
         coef1(p, t) = exp(-1im * p.ω1 * t)
         coef2(p, t) = sin(p.ω2 * t)
+        coef3(p, t) = sin(p.ω3 * t)
         t = rand()
-        p = (ω1 = rand(), ω2 = rand())
+        p = (ω1 = rand(), ω2 = rand(), ω3 = rand())
 
         # Operator
         H_td = QobjEvo(((a, coef1), a' * a, (a', coef2)))
         H_ti = coef1(p, t) * a + a' * a + coef2(p, t) * a'
         ψ = rand_ket(N)
         @test H_td(p, t) ≈ H_ti
+        @test iscached(H_td) == true
+        H_td = cache_operator(H_td, ψ)
+        @test iscached(H_td) == true
         @test H_td(ψ, p, t) ≈ H_ti * ψ
         @test isconstant(a) == true
         @test isconstant(H_td) == false
@@ -180,17 +184,36 @@
         @test isoper(H_td) == true
 
         # SuperOperator
-        c_ops = [sqrt(rand()) * a]
-        L_td = liouvillian(H_td, c_ops)
-        L_td2 = -1im * spre(H_td) + 1im * spost(H_td) + lindblad_dissipator(c_ops[1])
+        X = a * a'
+        c_op1 = QobjEvo(((a', coef1),))
+        c_op2 = QobjEvo(((a, coef2), (X, coef3)))
+        c_ops = [c_op1, c_op2]
+        D1_ti = abs2(coef1(p, t)) * lindblad_dissipator(a')
+        D2_ti =
+            abs2(coef2(p, t)) * lindblad_dissipator(a) + # normal dissipator for first  element in c_op2
+            abs2(coef3(p, t)) * lindblad_dissipator(X) + # normal dissipator for second element in c_op2
+            coef2(p, t) * conj(coef3(p, t)) * (spre(a) * spost(X') - 0.5 * spre(X' * a) - 0.5 * spost(X' * a)) + # cross terms
+            conj(coef2(p, t)) * coef3(p, t) * (spre(X) * spost(a') - 0.5 * spre(a' * X) - 0.5 * spost(a' * X))   # cross terms
+        L_ti = liouvillian(H_ti) + D1_ti + D2_ti
+        L_td = @test_logs (:warn,) (:warn,) liouvillian(H_td, c_ops) # warnings from lazy tensor in `lindblad_dissipator(c_op2)`
         ρvec = mat2vec(rand_dm(N))
-        @test L_td(p, t) ≈ L_td2(p, t)
-        @test L_td(ρvec, p, t) ≈ L_td2(ρvec, p, t)
+        @test L_td(p, t) ≈ L_ti
+        @test iscached(L_td) == false
+        L_td = cache_operator(L_td, ρvec)
+        @test iscached(L_td) == true
+        @test L_td(ρvec, p, t) ≈ L_ti * ρvec
         @test isconstant(L_td) == false
         @test issuper(L_td) == true
 
+        @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])
         @test_throws ArgumentError H_td(ρvec, p, t)
+        @test_throws ArgumentError cache_operator(H_td, ρvec)
         @test_throws ArgumentError L_td(ψ, p, t)
+        @test_throws ArgumentError cache_operator(L_td, ψ)
+
+        @testset "Type Inference" begin
+            @inferred liouvillian(H_td, (a, QobjEvo(((a', coef1),))))
+        end
     end
 end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index fd7a33f91..dcfb5f665 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -229,6 +229,8 @@
         # @test sol_sse.expect ≈ sol_sse_td2.expect atol = 1e-2 * length(tlist)
 
         @testset "Type Inference mesolve" begin
+            coef(p, t) = exp(-t)
+            ad_t = QobjEvo(((a', coef),))
             @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(
@@ -242,7 +244,7 @@
             @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, ψ0, tlist, (a, ad_t), e_ops = (a' * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple
             @inferred mesolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
             @inferred mesolve(H_td2, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
             @inferred mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)