generalize the definition of liouvillian

Author: ytdHuang

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. ([#581])
+
 ## [v0.38.1]
 Release date: 2025-10-27
 
@@ -357,3 +359,4 @@ Release date: 2024-11-13
 [#575]: https://github.com/qutip/QuantumToolbox.jl/issues/575
 [#576]: https://github.com/qutip/QuantumToolbox.jl/issues/576
 [#579]: https://github.com/qutip/QuantumToolbox.jl/issues/579
+[#581]: https://github.com/qutip/QuantumToolbox.jl/issues/581

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,4 +180,4 @@ liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.di
 
 liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal) = H
 
-_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
+_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops; init = 0)

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)

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