Add keyword argument assume_hermitian to liouvillian

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

src/qobj/superoperators.jl

@@ -39,17 +39,42 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(B_T, Id)
 end
 
-## intrinsic liouvillian 
-_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(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))
+_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_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_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}}
     Od_O = O' * O
-    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
+    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
@@ -139,12 +164,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[\hat{H}, \cdot] + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]
+\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],
 ```
 
 where 
@@ -156,27 +194,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)

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.
 

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

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

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

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

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