diff --git a/src/HeomBase.jl b/src/HeomBase.jl
index 92a1be51..be5c98aa 100644
--- a/src/HeomBase.jl
+++ b/src/HeomBase.jl
@@ -39,28 +39,44 @@ end
 _Tr(M::AbstractHEOMLSMatrix) = _Tr(_get_SciML_matrix_wrapper(M), M.dims, M.N)
 _Tr(M::Type{<:SparseMatrixCSC}, dims::SVector, N::Int) = _Tr(eltype(M), dims, N)
 
-function HandleMatrixType(M::AbstractQuantumObject, MatrixName::String = ""; type::QuantumObjectType = Operator)
-    if M.type == type
-        return M
-    else
+function HandleMatrixType(
+    M::AbstractQuantumObject,
+    MatrixName::String = "";
+    type::T = nothing,
+) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}}
+    if (type isa Nothing)
+        (M.type isa OperatorQuantumObject) ||
+            (M.type isa SuperOperatorQuantumObject) ||
+            error("The matrix $(MatrixName) should be either an Operator or a SuperOperator.")
+    elseif M.type != type
         error("The matrix $(MatrixName) should be an $(type).")
     end
+    return M
 end
 function HandleMatrixType(
     M::AbstractQuantumObject,
     dims::SVector,
     MatrixName::String = "";
-    type::QuantumObjectType = Operator,
-)
+    type::T = nothing,
+) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}}
     if M.dims == dims
         return HandleMatrixType(M, MatrixName; type = type)
     else
         error("The dims of $(MatrixName) should be: $(dims)")
     end
 end
-HandleMatrixType(M, dims::SVector, MatrixName::String = ""; type::QuantumObjectType = Operator) =
+HandleMatrixType(
+    M,
+    dims::SVector,
+    MatrixName::String = "";
+    type::T = nothing,
+) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}} =
     HandleMatrixType(M, MatrixName; type = type)
-HandleMatrixType(M, MatrixName::String = ""; type::QuantumObjectType = Operator) =
+HandleMatrixType(
+    M,
+    MatrixName::String = "";
+    type::T = nothing,
+) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}} =
     error("HierarchicalEOM doesn't support matrix $(MatrixName) with type : $(typeof(M))")
 
 # change the type of `ADOs` to match the type of HEOMLS matrix
diff --git a/src/bath/BosonBath.jl b/src/bath/BosonBath.jl
index 0408f816..2b89c456 100644
--- a/src/bath/BosonBath.jl
+++ b/src/bath/BosonBath.jl
@@ -8,14 +8,14 @@ Base.show(io::IO, B::AbstractBosonBath) =
 Base.show(io::IO, m::MIME"text/plain", B::AbstractBosonBath) = show(io, B)
 
 function _check_bosonic_coupling_operator(op)
-    _op = HandleMatrixType(op, "op (coupling operator)")
+    _op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
     if !ishermitian(_op)
         @warn "The system-bosonic-bath coupling operator \"op\" should be Hermitian Operator."
     end
     return _op
 end
 
-_check_bosonic_RWA_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)")
+_check_bosonic_RWA_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)", type = Operator)
 
 function _combine_same_gamma(η::Vector{Ti}, γ::Vector{Tj}) where {Ti,Tj<:Number}
     if length(η) != length(γ)
@@ -96,7 +96,7 @@ function BosonBath(
     δ::Number = 0.0;
     combine::Bool = true,
 ) where {Ti,Tj<:Number}
-    _op = HandleMatrixType(op, "op (coupling operator)")
+    _op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
     if combine
         ηnew, γnew = _combine_same_gamma(η, γ)
         bRI = bosonRealImag(_op, real.(ηnew), imag.(ηnew), γnew)
@@ -144,7 +144,7 @@ function BosonBath(
     δ::Tm = 0.0;
     combine::Bool = true,
 ) where {Ti,Tj,Tk,Tl,Tm<:Number}
-    _op = HandleMatrixType(op, "op (coupling operator)")
+    _op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
 
     if combine
         ηR, γR = _combine_same_gamma(η_real, γ_real)
@@ -366,7 +366,7 @@ function BosonBathRWA(
     δ::Tm = 0.0,
 ) where {Ti,Tj,Tk,Tl,Tm<:Number}
     _check_gamma_absorb_and_emit(γ_absorb, γ_emit)
-    _op = HandleMatrixType(op, "op (coupling operator)")
+    _op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
 
     bA = bosonAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
     bE = bosonEmit(_op, η_emit, γ_emit, η_absorb)
diff --git a/src/bath/FermionBath.jl b/src/bath/FermionBath.jl
index 53deba42..7ddd2883 100644
--- a/src/bath/FermionBath.jl
+++ b/src/bath/FermionBath.jl
@@ -7,7 +7,7 @@ Base.show(io::IO, B::AbstractFermionBath) =
     print(io, "$(typeof(B))-type bath with $(B.Nterm) exponential-expansion terms\n")
 Base.show(io::IO, m::MIME"text/plain", B::AbstractFermionBath) = show(io, B)
 
-_check_fermionic_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)")
+_check_fermionic_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)", type = Operator)
 
 @doc raw"""
     struct FermionBath <: AbstractBath
@@ -72,7 +72,7 @@ function FermionBath(
 ) where {Ti,Tj,Tk,Tl,Tm<:Number}
     _check_gamma_absorb_and_emit(γ_absorb, γ_emit)
 
-    _op = HandleMatrixType(op, "op (coupling operator)")
+    _op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
     fA = fermionAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
     fE = fermionEmit(_op, η_emit, γ_emit, η_absorb)
     return FermionBath(AbstractFermionBath[fA, fE], _op, fA.Nterm + fE.Nterm, δ)
diff --git a/src/heom_matrices/M_Boson.jl b/src/heom_matrices/M_Boson.jl
index 6a3ab21f..4d86d643 100644
--- a/src/heom_matrices/M_Boson.jl
+++ b/src/heom_matrices/M_Boson.jl
@@ -41,7 +41,7 @@ end
 Generate the boson-type HEOM Liouvillian superoperator matrix
 
 # Parameters
-- `Hsys` : The time-independent system Hamiltonian
+- `Hsys` : The time-independent system Hamiltonian or Liouvillian
 - `tier::Int` : the tier (cutoff level) for the bosonic bath
 - `Bath::Vector{BosonBath}` : objects for different bosonic baths
 - `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
@@ -62,7 +62,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 )
 
     # check for system dimension
-    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
+    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
     sup_dim = prod(_Hsys.dims)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
diff --git a/src/heom_matrices/M_Boson_Fermion.jl b/src/heom_matrices/M_Boson_Fermion.jl
index 8a427905..67cea6ee 100644
--- a/src/heom_matrices/M_Boson_Fermion.jl
+++ b/src/heom_matrices/M_Boson_Fermion.jl
@@ -73,7 +73,7 @@ end
 Generate the boson-fermion-type HEOM Liouvillian superoperator matrix
 
 # Parameters
-- `Hsys` : The time-independent system Hamiltonian
+- `Hsys` : The time-independent system Hamiltonian or Liouvillian
 - `Btier::Int` : the tier (cutoff level) for the bosonic bath
 - `Ftier::Int` : the tier (cutoff level) for the fermionic bath
 - `Bbath::Vector{BosonBath}` : objects for different bosonic baths
@@ -99,7 +99,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 )
 
     # check for system dimension
-    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
+    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
     sup_dim = prod(_Hsys.dims)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
diff --git a/src/heom_matrices/M_Fermion.jl b/src/heom_matrices/M_Fermion.jl
index 93a413e1..9bbe8be6 100644
--- a/src/heom_matrices/M_Fermion.jl
+++ b/src/heom_matrices/M_Fermion.jl
@@ -41,7 +41,7 @@ end
 Generate the fermion-type HEOM Liouvillian superoperator matrix
 
 # Parameters
-- `Hsys` : The time-independent system Hamiltonian
+- `Hsys` : The time-independent system Hamiltonian or Liouvillian
 - `tier::Int` : the tier (cutoff level) for the fermionic bath
 - `Bath::Vector{FermionBath}` : objects for different fermionic baths
 - `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
@@ -60,7 +60,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
 )
 
     # check for system dimension
-    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
+    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
     sup_dim = prod(_Hsys.dims)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
diff --git a/src/heom_matrices/M_S.jl b/src/heom_matrices/M_S.jl
index 0d4fc555..aab502b7 100644
--- a/src/heom_matrices/M_S.jl
+++ b/src/heom_matrices/M_S.jl
@@ -36,7 +36,7 @@ M[\cdot]=-i \left[H_{sys}, \cdot \right]_-,
 where ``[\cdot, \cdot]_-`` stands for commutator.
 
 # Parameters
-- `Hsys` : The time-independent system Hamiltonian
+- `Hsys` : The time-independent system Hamiltonian or Liouvillian
 - `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
 - `verbose::Bool` : To display verbose output during the process or not. Defaults to `true`.
 
@@ -45,7 +45,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 @noinline function M_S(Hsys::QuantumObject, parity::AbstractParity = EVEN; verbose::Bool = true)
 
     # check for system dimension
-    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
+    _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
     sup_dim = prod(_Hsys.dims)^2
 
     # the Liouvillian operator for free Hamiltonian
diff --git a/src/heom_matrices/heom_matrix_base.jl b/src/heom_matrices/heom_matrix_base.jl
index 5183b40e..2285441d 100644
--- a/src/heom_matrices/heom_matrix_base.jl
+++ b/src/heom_matrices/heom_matrix_base.jl
@@ -442,7 +442,7 @@ function bath_sum_γ(nvec, baths::Vector{T}) where {T<:Union{AbstractBosonBath,A
 end
 
 # commutator of system Hamiltonian
-minus_i_L_op(Hsys::QuantumObject, Id = I(size(Hsys, 1))) = _liouvillian(Hsys.data, Id)
+minus_i_L_op(Hsys::QuantumObject, Id = I(size(Hsys, 1))) = liouvillian(Hsys, Id).data
 
 # connect to bosonic (n-1)th-level for "Real & Imag combined operator"
 minus_i_D_op(bath::bosonRealImag, k, n_k) = n_k * (-1.0im * bath.η_real[k] * bath.Comm + bath.η_imag[k] * bath.anComm)
diff --git a/test/M_S.jl b/test/M_S.jl
index 298e7d7c..966456a5 100644
--- a/test/M_S.jl
+++ b/test/M_S.jl
@@ -4,11 +4,13 @@
     t = 10
     Hsys = sigmax()
     L = M_S(Hsys; verbose = false)
+    L_super = M_S(liouvillian(Hsys); verbose = false) # test if input is already Liouvillian
     ψ0 = basis(2, 0)
     @test show(devnull, MIME("text/plain"), L) === nothing
-    @test size(L) == (4, 4)
-    @test L.N == 1
-    @test nnz(L.data.A) == nnz(L(0)) == 8
+    @test size(L) == size(L_super) == (4, 4)
+    @test L.N == L_super.N == 1
+    @test nnz(L.data.A) == nnz(L_super.data.A) == nnz(L(0)) == nnz(L_super(0)) == 8
+    @test L.data == L_super.data
     ados_list = HEOMsolve(L, ψ0, 0:1:t; reltol = 1e-8, abstol = 1e-10, verbose = false).ados
     ados = ados_list[end]
     @test ados.dims == L.dims