diff --git a/CHANGELOG.md b/CHANGELOG.md
index 866e0ba20..9d75fd765 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Generalize the definition of `liouvillian`. It no longer expects the Hamiltonian to be Hermitian. ([#541])
+
 ## [v0.36.0]
 Release date: 2025-09-29
 
@@ -321,5 +323,6 @@ Release date: 2024-11-13
 [#536]: https://github.com/qutip/QuantumToolbox.jl/issues/536
 [#537]: https://github.com/qutip/QuantumToolbox.jl/issues/537
 [#539]: https://github.com/qutip/QuantumToolbox.jl/issues/539
+[#541]: https://github.com/qutip/QuantumToolbox.jl/issues/541
 [#544]: https://github.com/qutip/QuantumToolbox.jl/issues/544
 [#546]: https://github.com/qutip/QuantumToolbox.jl/issues/546
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 6bc9a71c3..acf21a9cb 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -41,9 +41,10 @@ end
 
 ## intrinsic liouvillian 
 _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
-    -1im * (_spre(H, Id) - _spost(H, Id))
+    -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::ScaledOperator, Id::AbstractMatrix) =
+    -1im * (ScaledOperator(H.λ, _spre(H.L, Id)) - ScaledOperator(conj(H.λ), _spost(H.L', Id)))
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
 
 # intrinsic lindblad_dissipator
@@ -144,7 +145,7 @@ lindblad_dissipator(O::AbstractQuantumObject{SuperOperator}, Id_cache = nothing)
 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]
+\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}^\dagger\right) + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
 ```
 
 where 
@@ -179,19 +180,4 @@ liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.di
 
 liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal) = H
 
-function _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal)
-    D = 0
-    # sum all the (time-independent) c_ops first
-    c_ops_ti = filter(op -> isa(op, QuantumObject), c_ops)
-    if !isempty(c_ops_ti)
-        D += 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)
-        D += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_td)
-    end
-
-    return D
-end
+_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops; init = 0)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index a417f52ac..9e8479449 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -207,7 +207,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
@@ -237,10 +237,9 @@
         @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)
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index 6151af349..b95273631 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -64,9 +64,8 @@
     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_fourier(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SteadyStateLinearSolver())[1]
-    ρ_ss2 =
-        steadystate_fourier(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian())
+    ρ_ss1 = steadystate_fourier(H, 0.5 * H_t, 0.5 * H_t, 1, c_ops, solver = SteadyStateLinearSolver())[1]
+    ρ_ss2 = steadystate_fourier(H, 0.5 * H_t, 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian())
 
     @test abs(sum(sol_me.expect[1, (end-100):end]) / 101 - expect(e_ops[1], ρ_ss1)) < 1e-3
     @test abs(sum(sol_me.expect[1, (end-100):end]) / 101 - expect(e_ops[1], ρ_ss2)) < 1e-3