From 205d36fb9947ee7501af7201389f41e63e1f9f71 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 10:36:28 +0800
Subject: [PATCH 01/17] generalize the definition of `liouvillian`

---
 CHANGELOG.md                          | 3 +++
 src/qobj/superoperators.jl            | 9 +++++----
 test/core-test/quantum_objects_evo.jl | 9 ++++-----
 test/core-test/steady_state.jl        | 5 ++---
 4 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3421cd799..a7166319c 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. ([#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
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 6f6b980d7..ebf94b320 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,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)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index cb86a74fe..98cbc64b6 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

From cea1d417082926f244620ba112faa304f34a0611 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 12:24:50 +0800
Subject: [PATCH 02/17] try to fix autodiff issue

---
 src/qobj/superoperators.jl | 14 +++++++++-----
 src/utilities.jl           |  2 ++
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index ebf94b320..28ab2db98 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -40,17 +40,21 @@ 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))
+function _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
+    CType = _complex_float_type(H)
+    CType(-1.0im) * _spre(H, Id) + CType(1.0im) * _spost(H', Id)
+end
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
-_liouvillian(H::ScaledOperator, Id::AbstractMatrix) =
-    -1im * (ScaledOperator(H.λ, _spre(H.L, Id)) - ScaledOperator(conj(H.λ), _spost(H.L', Id)))
+function _liouvillian(H::ScaledOperator, Id::AbstractMatrix)
+    CType = _complex_float_type(H)
+    CType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) + CType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
+end
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
 
 # 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
+    return _sprepost(O, O') + _complex_float_type(H)(-0.5 + 0.0im) * (_spre(Od_O, Id) + _spost(Od_O, Id))
 end
 function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
     _O = O.A
diff --git a/src/utilities.jl b/src/utilities.jl
index 6d647c08b..5bbf41064 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -175,6 +175,7 @@ end
 
 # functions for getting Float or Complex element type
 _float_type(::AbstractArray{T}) where {T<:Number} = _float_type(T)
+_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _float_type(::Type{Int32}) = Float32
 _float_type(::Type{Int64}) = Float64
 _float_type(::Type{Float32}) = Float32
@@ -186,6 +187,7 @@ _float_type(::Type{Complex{Float64}}) = Float64
 _float_type(::Type{Complex{T}}) where {T<:Real} = T
 _float_type(T::Type{<:Real}) = T # Allow other untracked Real types, like ForwardDiff.Dual
 _complex_float_type(::AbstractArray{T}) where {T<:Number} = _complex_float_type(T)
+_complex_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _complex_float_type(::Type{Int32}) = ComplexF32
 _complex_float_type(::Type{Int64}) = ComplexF64
 _complex_float_type(::Type{Float32}) = ComplexF32

From 08c89f75ce75774b08e9e9630138697aa0ef9298 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 12:25:14 +0800
Subject: [PATCH 03/17] format files

---
 src/qobj/superoperators.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 28ab2db98..99728bdb1 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -42,12 +42,13 @@ end
 ## intrinsic liouvillian 
 function _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
     CType = _complex_float_type(H)
-    CType(-1.0im) * _spre(H, Id) + CType(1.0im) * _spost(H', Id)
+    return CType(-1.0im) * _spre(H, Id) + CType(1.0im) * _spost(H', Id)
 end
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
 function _liouvillian(H::ScaledOperator, Id::AbstractMatrix)
     CType = _complex_float_type(H)
-    CType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) + CType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
+    return CType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
+           CType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
 end
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
 

From 09578e682499ec6ee472ce6ee10c1b22da010b65 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 12:28:28 +0800
Subject: [PATCH 04/17] fix typo

---
 src/qobj/superoperators.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 99728bdb1..a62cf85ae 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -55,7 +55,7 @@ _liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _li
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
     Od_O = O' * O
-    return _sprepost(O, O') + _complex_float_type(H)(-0.5 + 0.0im) * (_spre(Od_O, Id) + _spost(Od_O, Id))
+    return _sprepost(O, O') + _complex_float_type(O)(-0.5 + 0.0im) * (_spre(Od_O, Id) + _spost(Od_O, Id))
 end
 function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
     _O = O.A

From 0c001576fadfa2cec2914edb9ba170ee68fe04f5 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 13:56:12 +0800
Subject: [PATCH 05/17] add warning message

---
 src/qobj/superoperators.jl            | 19 +++++++++++--------
 test/core-test/quantum_objects_evo.jl |  4 ++--
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index a62cf85ae..42bc0c42f 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -41,14 +41,14 @@ end
 
 ## intrinsic liouvillian 
 function _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
-    CType = _complex_float_type(H)
-    return CType(-1.0im) * _spre(H, Id) + CType(1.0im) * _spost(H', Id)
+    CFType = _complex_float_type(H)
+    return CFType(-1.0im) * _spre(H, Id) + CFType(1.0im) * _spost(H', Id)
 end
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
 function _liouvillian(H::ScaledOperator, Id::AbstractMatrix)
-    CType = _complex_float_type(H)
-    return CType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
-           CType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
+    CFType = _complex_float_type(H)
+    return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
+           CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
 end
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
 
@@ -180,9 +180,12 @@ liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal =
 
 liouvillian(H::Nothing, c_ops::Nothing) = 0
 
-liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.dimensions))) =
-    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator(), H.dimensions)
+function liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.dimensions)))
+    @warn "The definition of `liouvillian` L for a given Hamiltonian H is changed to L[⋅] = -i( H[⋅] - [⋅]H' ), with ' represents complex conjugation" *
+          "QuantumToolbox.jl no longer expects H to be Hermitian." maxlog = 1
+    return get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator(), H.dimensions)
+end
 
 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; init = 0)
+_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 98cbc64b6..c1f407b80 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -215,7 +215,7 @@
             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)`
+        L_td = @test_logs (:warn,) (:warn,) (:warn,) liouvillian(H_td, c_ops) # two warnings from lazy tensor `lindblad_dissipator(c_op2)`, one from changing liouvillian definition
         ρvec = mat2vec(rand_dm(N))
         @test L_td(p, t) ≈ L_ti
         @test iscached(L_td) == false
@@ -227,7 +227,7 @@
 
         coef_wrong1(t) = nothing
         coef_wrong2(p, t::ComplexF64) = nothing
-        @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
+        @test_logs (:warn,) (:warn,) (:warn,) liouvillian(H_td * H_td) # two warnings from lazy tensor, one from changing liouvillian definition
         @test_throws ArgumentError QobjEvo(a, coef_wrong1)
         @test_throws ArgumentError QobjEvo(a, coef_wrong2)
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])

From 58409f6fb022e79993a9d9f52cd005a4a64fbabb Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 14:24:29 +0800
Subject: [PATCH 06/17] remove warnings

---
 src/qobj/superoperators.jl            | 7 ++-----
 test/core-test/quantum_objects_evo.jl | 4 ++--
 2 files changed, 4 insertions(+), 7 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 42bc0c42f..ec4d479e0 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -180,11 +180,8 @@ liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal =
 
 liouvillian(H::Nothing, c_ops::Nothing) = 0
 
-function liouvillian(H::AbstractQuantumObject{Operator}, Id_cache::Diagonal = I(prod(H.dimensions)))
-    @warn "The definition of `liouvillian` L for a given Hamiltonian H is changed to L[⋅] = -i( H[⋅] - [⋅]H' ), with ' represents complex conjugation" *
-          "QuantumToolbox.jl no longer expects H to be Hermitian." maxlog = 1
-    return get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator(), H.dimensions)
-end
+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) = H
 
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index c1f407b80..98cbc64b6 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -215,7 +215,7 @@
             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,) (:warn,) liouvillian(H_td, c_ops) # two warnings from lazy tensor `lindblad_dissipator(c_op2)`, one from changing liouvillian definition
+        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_ti
         @test iscached(L_td) == false
@@ -227,7 +227,7 @@
 
         coef_wrong1(t) = nothing
         coef_wrong2(p, t::ComplexF64) = nothing
-        @test_logs (:warn,) (:warn,) (:warn,) liouvillian(H_td * H_td) # two warnings from lazy tensor, one from changing liouvillian definition
+        @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
         @test_throws ArgumentError QobjEvo(a, coef_wrong1)
         @test_throws ArgumentError QobjEvo(a, coef_wrong2)
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])

From 341000a1642f8120b182e9d8ead1f1c521e48bc7 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Tue, 28 Oct 2025 17:10:51 +0800
Subject: [PATCH 07/17] add warning message in module initialization

---
 src/QuantumToolbox.jl | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 2da2b5655..cfffd4ecb 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -134,4 +134,13 @@ include("visualization.jl")
 # deprecated functions
 include("deprecated.jl")
 
+# module initialization
+function __init__()
+    # this will be called immediately after the module is loaded (e.g., by import or using) at runtime for the first time.
+    # useful for announcement
+    @warn "The definition of `liouvillian` L for a given Hamiltonian H is changed to L[⋅] = -i( H[⋅] - [⋅]H' ), with ' representing complex conjugation. " *
+          "QuantumToolbox.jl no longer expects H to be Hermitian."
+    return nothing
+end
+
 end

From 12d149a06d626f6ad8288d79c30f0f22f953e683 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Wed, 29 Oct 2025 15:37:00 +0800
Subject: [PATCH 08/17] add keyword argument `assume_hermitian`

---
 CHANGELOG.md                           |  2 +-
 src/QuantumToolbox.jl                  |  9 ---
 src/qobj/superoperators.jl             | 86 ++++++++++++++++++++------
 src/settings.jl                        |  2 +-
 test/core-test/quantum_objects_evo.jl  | 20 ++++--
 test/core-test/time_evolution.jl       |  5 ++
 test/ext-test/cpu/autodiff/autodiff.jl | 36 +++++++++--
 7 files changed, 119 insertions(+), 41 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 063ebc179..02bbec4fe 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,7 +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]) 
-- Generalize the definition of `liouvillian`. It no longer expects the Hamiltonian to be Hermitian. ([#581])
+- 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
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index cfffd4ecb..2da2b5655 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -134,13 +134,4 @@ include("visualization.jl")
 # deprecated functions
 include("deprecated.jl")
 
-# module initialization
-function __init__()
-    # this will be called immediately after the module is loaded (e.g., by import or using) at runtime for the first time.
-    # useful for announcement
-    @warn "The definition of `liouvillian` L for a given Hamiltonian H is changed to L[⋅] = -i( H[⋅] - [⋅]H' ), with ' representing complex conjugation. " *
-          "QuantumToolbox.jl no longer expects H to be Hermitian."
-    return nothing
-end
-
 end
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index ec4d479e0..394419e7a 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -39,18 +39,38 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(B_T, Id)
 end
 
-## intrinsic liouvillian 
-function _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
+## intrinsic liouvillian
+_liouvillian(
+    H::MT,
+    Id::AbstractMatrix,
+    assume_hermitian::Val{true},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_Herm(H, Id)
+_liouvillian(
+    H::MT,
+    Id::AbstractMatrix,
+    assume_hermitian::Val{false},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_NonHerm(H, Id)
+
+## intrinsic liouvillian (assume H is Hermitian)
+_liouvillian_Herm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
+    -1im * (_spre(H, Id) - _spost(H, Id))
+_liouvillian_Herm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_Herm(H.A, Id))
+_liouvillian_Herm(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian_Herm(H.L, Id))
+_liouvillian_Herm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_Herm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
+
+## intrinsic liouvillian (assume H is Non-Hermitian)
+function _liouvillian_NonHerm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
     CFType = _complex_float_type(H)
     return CFType(-1.0im) * _spre(H, Id) + CFType(1.0im) * _spost(H', Id)
 end
-_liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
-function _liouvillian(H::ScaledOperator, Id::AbstractMatrix)
+_liouvillian_NonHerm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_NonHerm(H.A, Id))
+function _liouvillian_NonHerm(H::ScaledOperator, Id::AbstractMatrix)
     CFType = _complex_float_type(H)
     return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
            CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
 end
-_liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
+_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) =
+    mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
@@ -145,12 +165,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``:
+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],
+```
+
+otherwise [`assume_hermitian = Val(false)` or `false`],
 
 ```math
-\mathcal{L} [\cdot] = -i\left(\hat{H}[\cdot] - [\cdot]\hat{H}^\dagger\right) + \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 
@@ -162,27 +195,42 @@ 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)
+    L = liouvillian(H, Id_cache; assume_hermitian = assume_hermitian)
     if !(c_ops isa Nothing)
         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))) =
-    _sum_lindblad_dissipators(c_ops, Id_cache)
-
-liouvillian(H::Nothing, c_ops::Nothing) = 0
-
-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) = H
+liouvillian(
+    H::Nothing,
+    c_ops::Union{AbstractVector,Tuple},
+    Id_cache::Diagonal = I(prod(c_ops[1].dims));
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = _sum_lindblad_dissipators(c_ops, Id_cache)
+
+liouvillian(H::Nothing, c_ops::Nothing; assume_hermitian::Union{Bool,Val} = Val(true)) = 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{SuperOperator},
+    Id_cache::Diagonal;
+    assume_hermitian::Union{Bool,Val} = Val(true),
+) = H
 
 _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
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/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 98cbc64b6..47bac1df6 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -182,15 +182,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
@@ -225,6 +225,13 @@
         @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 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
@@ -245,6 +252,7 @@
             @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..4feba8b2d 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -1046,6 +1046,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

From 0da8418417d6e745ca3500e5512de04ee1c489e9 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Wed, 29 Oct 2025 16:00:56 +0800
Subject: [PATCH 09/17] format files

---
 src/qobj/superoperators.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 394419e7a..a649be7cf 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -69,8 +69,7 @@ function _liouvillian_NonHerm(H::ScaledOperator, Id::AbstractMatrix)
     return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
            CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
 end
-_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) =
-    mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
+_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}

From 3cf00ec71c1a778620ddd4bae40b4bb96a676b66 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Mon, 3 Nov 2025 11:21:16 +0800
Subject: [PATCH 10/17] improve QobjEvo tests

---
 test/core-test/quantum_objects_evo.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 47bac1df6..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
@@ -228,6 +229,8 @@
         # 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

From a16abd3b95d7a0e1d1763624f4b1b019fe00f446 Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <61953577+albertomercurio@users.noreply.github.com>
Date: Wed, 5 Nov 2025 14:11:38 +0100
Subject: [PATCH 11/17] Adapt to new version of SciMLOperators.jl (#585)

---
 src/qobj/superoperators.jl                    | 59 ++++++-------------
 .../time_evolution_dynamical.jl               |  5 +-
 src/utilities.jl                              |  2 -
 test/core-test/time_evolution.jl              |  6 +-
 4 files changed, 26 insertions(+), 46 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index a649be7cf..5b8765d03 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -39,42 +39,29 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(B_T, Id)
 end
 
-## intrinsic liouvillian
-_liouvillian(
+## intrinsic liouvillian 
+function _liouvillian(
     H::MT,
     Id::AbstractMatrix,
-    assume_hermitian::Val{true},
-) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_Herm(H, Id)
-_liouvillian(
-    H::MT,
-    Id::AbstractMatrix,
-    assume_hermitian::Val{false},
-) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_NonHerm(H, Id)
-
-## intrinsic liouvillian (assume H is Hermitian)
-_liouvillian_Herm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
-    -1im * (_spre(H, Id) - _spost(H, Id))
-_liouvillian_Herm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_Herm(H.A, Id))
-_liouvillian_Herm(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian_Herm(H.L, Id))
-_liouvillian_Herm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_Herm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
-
-## intrinsic liouvillian (assume H is Non-Hermitian)
-function _liouvillian_NonHerm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
-    CFType = _complex_float_type(H)
-    return CFType(-1.0im) * _spre(H, Id) + CFType(1.0im) * _spost(H', Id)
+    ::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
-_liouvillian_NonHerm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_NonHerm(H.A, Id))
-function _liouvillian_NonHerm(H::ScaledOperator, Id::AbstractMatrix)
-    CFType = _complex_float_type(H)
-    return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
-           CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
+function _liouvillian(H::MatrixOperator, Id::AbstractMatrix, assume_hermitian::Val)
+    isconstant(H) || error("Liouvillian can only be constructed from constant MatrixOperator.")
+    return MatrixOperator(_liouvillian(H.A, Id, assume_hermitian))
 end
-_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
+_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}}
     Od_O = O' * O
-    return _sprepost(O, O') + _complex_float_type(O)(-0.5 + 0.0im) * (_spre(Od_O, Id) + _spost(Od_O, Id))
+    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
 end
 function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
     _O = O.A
@@ -211,14 +198,10 @@ function liouvillian(
     return L
 end
 
-liouvillian(
-    H::Nothing,
-    c_ops::Union{AbstractVector,Tuple},
-    Id_cache::Diagonal = I(prod(c_ops[1].dims));
-    assume_hermitian::Union{Bool,Val} = Val(true),
-) = _sum_lindblad_dissipators(c_ops, Id_cache)
+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; assume_hermitian::Union{Bool,Val} = Val(true)) = 0
+liouvillian(H::Nothing, c_ops::Nothing; kwargs...) = 0
 
 liouvillian(
     H::AbstractQuantumObject{Operator},
@@ -226,10 +209,6 @@ liouvillian(
     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{SuperOperator},
-    Id_cache::Diagonal;
-    assume_hermitian::Union{Bool,Val} = Val(true),
-) = H
+liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal; kwargs...) = H
 
 _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
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/src/utilities.jl b/src/utilities.jl
index 5bbf41064..6d647c08b 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -175,7 +175,6 @@ end
 
 # functions for getting Float or Complex element type
 _float_type(::AbstractArray{T}) where {T<:Number} = _float_type(T)
-_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _float_type(::Type{Int32}) = Float32
 _float_type(::Type{Int64}) = Float64
 _float_type(::Type{Float32}) = Float32
@@ -187,7 +186,6 @@ _float_type(::Type{Complex{Float64}}) = Float64
 _float_type(::Type{Complex{T}}) where {T<:Real} = T
 _float_type(T::Type{<:Real}) = T # Allow other untracked Real types, like ForwardDiff.Dual
 _complex_float_type(::AbstractArray{T}) where {T<:Number} = _complex_float_type(T)
-_complex_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _complex_float_type(::Type{Int32}) = ComplexF32
 _complex_float_type(::Type{Int64}) = ComplexF64
 _complex_float_type(::Type{Float32}) = ComplexF32
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 4feba8b2d..6d67495e6 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

From 53a07a1000a9f5ed3de33154dd517a5768289561 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Wed, 5 Nov 2025 21:49:27 +0800
Subject: [PATCH 12/17] minor changes

---
 src/qobj/superoperators.jl     | 2 +-
 test/core-test/steady_state.jl | 5 +++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 5b8765d03..1a43175a6 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -50,7 +50,7 @@ function _liouvillian(
     return -1im * (H_spre - H_spost)
 end
 function _liouvillian(H::MatrixOperator, Id::AbstractMatrix, assume_hermitian::Val)
-    isconstant(H) || error("Liouvillian can only be constructed from constant MatrixOperator.")
+    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}) =
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index b95273631..6151af349 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -64,8 +64,9 @@
     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, 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())
+    ρ_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())
 
     @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

From 6d3a8091815fda0f7a6c6ce0a8898838af76c7ae Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Wed, 5 Nov 2025 22:31:20 +0800
Subject: [PATCH 13/17] format files

---
 src/qobj/superoperators.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 1a43175a6..c9b180f81 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -50,7 +50,8 @@ function _liouvillian(
     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."))
+    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}) =

From b6cd1e4d4acfc42ac34c4999ba98851c52993732 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Fri, 14 Nov 2025 11:49:36 +0800
Subject: [PATCH 14/17] update compat for `SciMLOperators`

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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"

From b7f4642131b6e32ed079c1251ac1778f3423095f Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Fri, 14 Nov 2025 10:49:50 +0100
Subject: [PATCH 15/17] Fix type instabilities with `lindblad_dissipator`

---
 src/qobj/superoperators.jl    | 18 +++++++++++++++---
 src/time_evolution/mesolve.jl |  3 ++-
 2 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index c9b180f81..290c8a693 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -193,8 +193,8 @@ function liouvillian(
     assume_hermitian::Union{Bool,Val} = Val(true),
 ) where {OpType<:Union{Operator,SuperOperator}}
     L = liouvillian(H, Id_cache; assume_hermitian = assume_hermitian)
-    if !(c_ops isa Nothing)
-        L += _sum_lindblad_dissipators(c_ops, Id_cache)
+    if !isnothing(c_ops)
+        return L + _sum_lindblad_dissipators(c_ops, Id_cache)
     end
     return L
 end
@@ -212,4 +212,16 @@ liouvillian(
 
 liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal; kwargs...) = H
 
-_sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)
+function _sum_lindblad_dissipators(c_ops::AbstractVector, Id_cache::Diagonal)
+    return sum(op -> lindblad_dissipator(op, Id_cache), c_ops; init = 0)
+end
+
+# 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)
+    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/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."))

From f904ec9be1560f71085768e40f2fe2068ce9dc82 Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Fri, 14 Nov 2025 11:08:51 +0100
Subject: [PATCH 16/17] Add nothing case for c_ops

---
 src/qobj/superoperators.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 290c8a693..2ae35ec11 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -212,6 +212,8 @@ liouvillian(
 
 liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal; kwargs...) = 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

From ea37ce55863ed2e688408bd03e209b6f4646df9b Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Fri, 14 Nov 2025 11:55:56 +0100
Subject: [PATCH 17/17] Add handling for zero length tuple in c_ops and adjust
 allocation test threshold

---
 src/qobj/superoperators.jl       | 3 +++
 test/core-test/time_evolution.jl | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 2ae35ec11..4ea7827fe 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -221,6 +221,9 @@ end
 # 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))
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 6d67495e6..95b462f41 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -525,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,