diff --git a/CHANGELOG.md b/CHANGELOG.md
index a5535d6f0..02bbec4fe 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Change default solver detection in `eigensolve` when using `sigma` keyword argument (shift-inverse algorithm). If the operator is a `SparseMatrixCSC`, the default solver is `UMFPACKFactorization`, otherwise it is automatically chosen by LinearSolve.jl, depending on the type of the operator. ([#580]) 
+- Add keyword argument `assume_hermitian` to `liouvillian`. This allows users to disable the assumption that the Hamiltonian is Hermitian. ([#581])
 
 ## [v0.38.1]
 Release date: 2025-10-27
@@ -360,3 +361,4 @@ Release date: 2024-11-13
 [#576]: https://github.com/qutip/QuantumToolbox.jl/issues/576
 [#579]: https://github.com/qutip/QuantumToolbox.jl/issues/579
 [#580]: https://github.com/qutip/QuantumToolbox.jl/issues/580
+[#581]: https://github.com/qutip/QuantumToolbox.jl/issues/581
diff --git a/Project.toml b/Project.toml
index faa1148b8..27f6d6db4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -64,7 +64,7 @@ Pkg = "1"
 ProgressMeter = "1.11.0"
 Random = "1"
 SciMLBase = "2.105"
-SciMLOperators = "1.4"
+SciMLOperators = "1.11"
 SparseArrays = "1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 6f6b980d7..4ea7827fe 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -40,11 +40,24 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
 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) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
+function _liouvillian(
+    H::MT,
+    Id::AbstractMatrix,
+    ::Val{assume_hermitian},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator},assume_hermitian}
+    H_spre = _spre(H, Id)
+    H_spost = assume_hermitian ? _spost(H, Id) : _spost(H', Id)
+    return -1im * (H_spre - H_spost)
+end
+function _liouvillian(H::MatrixOperator, Id::AbstractMatrix, assume_hermitian::Val)
+    isconstant(H) ||
+        throw(ArgumentError("The given Hamiltonian for constructing Liouvillian must be constant MatrixOperator."))
+    return MatrixOperator(_liouvillian(H.A, Id, assume_hermitian))
+end
+_liouvillian(H::ScaledOperator, Id::AbstractMatrix, assume_hermitian::Val{true}) =
+    ScaledOperator(H.λ, _liouvillian(H.L, Id, assume_hermitian))
+_liouvillian(H::AddedOperator, Id::AbstractMatrix, assume_hermitian::Val) =
+    AddedOperator(map(op -> _liouvillian(op, Id, assume_hermitian), H.ops))
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
@@ -139,12 +152,25 @@ lindblad_dissipator(O::AbstractQuantumObject{Operator}, Id_cache = I(size(O, 1))
 lindblad_dissipator(O::AbstractQuantumObject{SuperOperator}, Id_cache = nothing) = O
 
 @doc raw"""
-    liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dimensions)))
+    liouvillian(
+        H::AbstractQuantumObject,
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
+        Id_cache=I(prod(H.dimensions));
+        assume_hermitian::Union{Bool,Val} = Val(true),
+    )
+
+Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``.
+
+By default, when the Hamiltonian `H` is assumed to be Hermitian [`assume_hermitian = Val(true)` or `true`], the Liouvillian [`SuperOperator`](@ref) is defined as :
+
+```math
+\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot],
+```
 
-Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
+otherwise [`assume_hermitian = Val(false)` or `false`],
 
 ```math
-\mathcal{L} [\cdot] = -i[\hat{H}, \cdot] + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
+\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}^\dagger\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot],
 ```
 
 where 
@@ -156,27 +182,51 @@ where
 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).
+
+!!! warning "Beware of type-stability!"
+    If you want to keep type stability, it is recommended to use `assume_hermitian = Val(true)` instead of `assume_hermitian = true`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function liouvillian(
     H::AbstractQuantumObject{OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    Id_cache::Diagonal = I(prod(H.dimensions)),
+    Id_cache::Diagonal = I(prod(H.dimensions));
+    assume_hermitian::Union{Bool,Val} = Val(true),
 ) where {OpType<:Union{Operator,SuperOperator}}
-    L = liouvillian(H, Id_cache)
-    if !(c_ops isa Nothing)
-        L += _sum_lindblad_dissipators(c_ops, Id_cache)
+    L = liouvillian(H, Id_cache; assume_hermitian = assume_hermitian)
+    if !isnothing(c_ops)
+        return L + _sum_lindblad_dissipators(c_ops, Id_cache)
     end
     return L
 end
 
-liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal = I(prod(c_ops[1].dims))) =
+liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal = I(prod(c_ops[1].dims)); kwargs...) =
     _sum_lindblad_dissipators(c_ops, Id_cache)
 
-liouvillian(H::Nothing, c_ops::Nothing) = 0
+liouvillian(H::Nothing, c_ops::Nothing; kwargs...) = 0
+
+liouvillian(
+    H::AbstractQuantumObject{Operator},
+    Id_cache::Diagonal = I(prod(H.dimensions));
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = get_typename_wrapper(H)(_liouvillian(H.data, Id_cache, makeVal(assume_hermitian)), SuperOperator(), H.dimensions)
 
-liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.dimensions))) =
-    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator(), H.dimensions)
+liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal; kwargs...) = H
 
-liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal) = H
+_sum_lindblad_dissipators(c_ops::Nothing, Id_cache::Diagonal) = 0
+
+function _sum_lindblad_dissipators(c_ops::AbstractVector, Id_cache::Diagonal)
+    return sum(op -> lindblad_dissipator(op, Id_cache), c_ops; init = 0)
+end
 
-_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
+# Help the compiler to unroll the sum at compile time
+@generated function _sum_lindblad_dissipators(c_ops::Tuple, Id_cache::Diagonal)
+    N = length(c_ops.parameters)
+    if N == 0
+        return :(0)
+    end
+    ex = :(lindblad_dissipator(c_ops[1], Id_cache))
+    for i in 2:N
+        ex = :($ex + lindblad_dissipator(c_ops[$i], Id_cache))
+    end
+    return ex
+end
diff --git a/src/settings.jl b/src/settings.jl
index 668e5f0e7..72ccbabef 100644
--- a/src/settings.jl
+++ b/src/settings.jl
@@ -15,7 +15,7 @@ function Base.show(io::IO, s::Settings)
 end
 
 @doc raw"""
-    QuantumToolbox.settings 
+    QuantumToolbox.settings
 
 Contains all the default global settings of QuantumToolbox.jl.
 
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index e91e6b27e..4a9a9c805 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -1,6 +1,7 @@
 export mesolveProblem, mesolve, mesolve_map
 
-_mesolve_make_L_QobjEvo(H::Union{QuantumObject,Nothing}, c_ops) = QobjEvo(liouvillian(H, c_ops); type = SuperOperator())
+_mesolve_make_L_QobjEvo(H::Union{QuantumObject,Nothing}, c_ops) =
+    QuantumObjectEvolution(liouvillian(H, c_ops), type = SuperOperator())
 _mesolve_make_L_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}, c_ops) = liouvillian(QobjEvo(H), c_ops)
 _mesolve_make_L_QobjEvo(H::Nothing, c_ops::Nothing) = throw(ArgumentError("Both H and
 c_ops are Nothing. You are probably running the wrong function."))
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 305822611..79e392f72 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -571,8 +571,9 @@ function _DSF_mcsolve_Affect!(integrator)
     c_ops0_herm .= map(op -> op' * op, c_ops0)
 
     H_nh = convert(eltype(ψt), 0.5im) * sum(c_ops0_herm)
-    # By doing this, we are assuming that the system is time-independent and f is a MatrixOperator
-    copyto!(integrator.f.f.A, lmul!(-1im, H(op_l2, dsf_params).data - H_nh))
+    # By doing this, we are assuming that the system is time-independent and f is a ScaledOperator
+    # of the form -1im * (H - H_nh)
+    copyto!(integrator.f.f.L, H(op_l2, dsf_params).data - H_nh)
     return u_modified!(integrator, true)
 end
 
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index cb86a74fe..1b596b325 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -3,6 +3,7 @@
     using SparseArrays
     using StaticArraysCore
     using SciMLOperators
+    import SciMLOperators: AddedOperator
 
     # DomainError: incompatible between size of array and type
     @testset "Thrown Errors" begin
@@ -182,15 +183,15 @@
     @testset "Time Dependent Operators and SuperOperators" begin
         N = 10
         a = destroy(N)
-        coef1(p, t) = exp(-1im * p.ω1 * t)
-        coef2(p, t) = sin(p.ω2 * t)
-        coef3(p, t) = sin(p.ω3 * t)
+        coef1(p, t) = exp(-p.γ * t)
+        coef2(p, t) = cos(p.ω1 * t)
+        coef3(p, t) = sin(p.ω2 * t)
         t = rand()
-        p = (ω1 = rand(), ω2 = rand(), ω3 = rand())
+        p = (γ = rand(), ω1 = rand(), ω2 = rand())
 
         # Operator
-        H_td = QobjEvo(((a, coef1), a' * a, (a', coef2)))
-        H_ti = coef1(p, t) * a + a' * a + coef2(p, t) * a'
+        H_td = QobjEvo(((a + a', coef1), a' * a, (-1im * (a - a'), coef2))) # a Hermitian Hamiltonian
+        H_ti = coef1(p, t) * (a + a') + a' * a - 1im * coef2(p, t) * (a - a')
         ψ = rand_ket(N)
         @test H_td(p, t) ≈ H_ti
         @test iscached(H_td) == true
@@ -207,7 +208,7 @@
         X = a * a'
         c_op1 = QobjEvo(a', coef1)
         c_op2 = QobjEvo(((a, coef2), (X, coef3)))
-        c_ops = [c_op1, c_op2]
+        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
@@ -225,6 +226,15 @@
         @test isconstant(L_td) == false
         @test issuper(L_td) == true
 
+        # test number of lazy operators and assume_hermitian = Val(false)
+        L_td_assume_herm = liouvillian(H_td)
+        L_td_assume_not_herm = liouvillian(H_td, assume_hermitian = Val(false))
+        @test L_td_assume_herm.data isa AddedOperator
+        @test L_td_assume_not_herm.data isa AddedOperator
+        @test length(L_td_assume_herm.data.ops) == 3 # 1-time-indep. + 2-time-dep.
+        @test length(L_td_assume_not_herm.data.ops) == 5 # 1-time-indep. + 2 x 2-time-dep.
+        @test L_td_assume_herm(p, t) ≈ L_td_assume_not_herm(p, t) # check the matrix since H_td is itself Hermitian
+
         coef_wrong1(t) = nothing
         coef_wrong2(p, t::ComplexF64) = nothing
         @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
@@ -237,15 +247,15 @@
         @test_throws ArgumentError cache_operator(L_td, ψ)
 
         @testset "Type Inference" begin
-            # we use destroy and create here because they somehow causes type instability before
-            H_td2 = H_td + QobjEvo(destroy(N) + create(N), coef3)
-            c_ops1 = (destroy(N), create(N))
-            c_ops2 = (destroy(N), QobjEvo(create(N), coef1))
+            H_td2 = H_td + QobjEvo(a + a', coef3)
+            c_ops1 = (a, a')
+            c_ops2 = (a, QobjEvo(a', coef1))
 
             @inferred liouvillian(H_td, c_ops1)
             @inferred liouvillian(H_td, c_ops2)
             @inferred liouvillian(H_td2, c_ops1)
             @inferred liouvillian(H_td2, c_ops2)
+            @inferred liouvillian(H_td2, c_ops2, assume_hermitian = Val(false))
         end
     end
 end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index c2100d541..95b462f41 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -62,6 +62,7 @@ end
 
 @testitem "sesolve" setup=[TESetup] begin
     using SciMLOperators
+    import SciMLOperators: ScaledOperator
 
     # Get parameters from TESetup to simplify the code
     H = TESetup.H
@@ -82,7 +83,7 @@ end
     amp_rabi = TESetup.g^2 / Ω_rabi^2
     ##
 
-    @test prob.prob.f.f isa MatrixOperator
+    @test prob.prob.f.f isa ScaledOperator
     @test sum(abs.(sol.expect[1, :] .- amp_rabi .* sin.(Ω_rabi * tlist) .^ 2)) / length(tlist) < 0.1
     @test length(sol.times) == length(tlist)
     @test length(sol.times_states) == 1
@@ -390,6 +391,7 @@ end
 
 @testitem "mcsolve" setup=[TESetup] begin
     using SciMLOperators
+    import SciMLOperators: ScaledOperator
     using Statistics
 
     # Get parameters from TESetup to simplify the code
@@ -432,7 +434,7 @@ end
     expect_mc_states_mean = expect.(Ref(e_ops[1]), average_states(sol_mc_states))
     expect_mc_states_mean2 = expect.(Ref(e_ops[1]), average_states(sol_mc_states2))
 
-    @test prob_mc.prob.f.f isa MatrixOperator
+    @test prob_mc.prob.f.f isa ScaledOperator
     @test sum(abs, sol_mc.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, sol_mc2.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, average_expect(sol_mc3) .- sol_me.expect) / length(tlist) < 0.1
@@ -523,7 +525,7 @@ end
                 progress_bar = Val(false),
                 keep_runs_results = Val(true),
             )
-            @test allocs_tot < n1 * ntraj + 400 # 150 allocations per trajectory + 500 for initialization
+            @test allocs_tot < n1 * ntraj + 600 # 150 allocations per trajectory + 600 for initialization
 
             allocs_tot = @allocations mcsolve(
                 H,
@@ -1046,6 +1048,11 @@ end
     @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)
+
+    # test assume_hermitian = Val(false)
+    L_td_non_herm = liouvillian(H_td2, c_ops, assume_hermitian = Val(false))
+    sol_me_td3 = mesolve(L_td_non_herm, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p)
+    @test sol_me.expect ≈ sol_me_td3.expect atol = 1e-6 * length(tlist)
 end
 
 @testitem "mcsolve, ssesolve and smesolve reproducibility" setup=[TESetup] begin
diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl
index 3f8e601a4..5141e8854 100644
--- a/test/ext-test/cpu/autodiff/autodiff.jl
+++ b/test/ext-test/cpu/autodiff/autodiff.jl
@@ -53,6 +53,7 @@ coef_γ(p, t) = sqrt(p[3])
 H = QobjEvo(a' * a, coef_Δ) + QobjEvo(a + a', coef_F)
 c_ops = [QobjEvo(a, coef_γ)]
 const L = liouvillian(H, c_ops)
+const L_assume_non_herm = liouvillian(H, c_ops, assume_hermitian = Val(false))
 
 function my_f_mesolve(p)
     sol = mesolve(
@@ -67,6 +68,19 @@ function my_f_mesolve(p)
     return real(expect(a' * a, sol.states[end]))
 end
 
+function my_f_mesolve_assume_non_herm(p)
+    sol = mesolve(
+        L_assume_non_herm,
+        ψ0_mesolve,
+        tlist_mesolve,
+        progress_bar = Val(false),
+        params = p,
+        sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()),
+    )
+
+    return real(expect(a' * a, sol.states[end]))
+end
+
 # Analytical solution
 n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
 
@@ -113,6 +127,7 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
 
         my_f_mesolve_direct(params)
         my_f_mesolve(params)
+        my_f_mesolve_assume_non_herm(params)
 
         grad_exact = Zygote.gradient((p) -> n_ss(p[1], p[2], p[3]), params)[1]
 
@@ -122,20 +137,31 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
         end
 
         @testset "Zygote.jl" begin
-            grad_qt = Zygote.gradient(my_f_mesolve, params)[1]
-            @test grad_qt ≈ grad_exact atol=1e-6
+            grad_qt1 = Zygote.gradient(my_f_mesolve, params)[1]
+            grad_qt2 = Zygote.gradient(my_f_mesolve_assume_non_herm, params)[1]
+            @test grad_qt1 ≈ grad_exact atol=1e-6
+            @test grad_qt2 ≈ grad_exact atol=1e-6
         end
 
         @testset "Enzyme.jl" begin
-            dparams = Enzyme.make_zero(params)
+            dparams1 = Enzyme.make_zero(params)
             Enzyme.autodiff(
                 Enzyme.set_runtime_activity(Enzyme.Reverse),
                 my_f_mesolve,
                 Active,
-                Duplicated(params, dparams),
+                Duplicated(params, dparams1),
             )[1]
 
-            @test dparams ≈ grad_exact atol=1e-6
+            dparams2 = Enzyme.make_zero(params)
+            Enzyme.autodiff(
+                Enzyme.set_runtime_activity(Enzyme.Reverse),
+                my_f_mesolve_assume_non_herm,
+                Active,
+                Duplicated(params, dparams2),
+            )[1]
+
+            @test dparams1 ≈ grad_exact atol=1e-6
+            @test dparams2 ≈ grad_exact atol=1e-6
         end
     end
 end