diff --git a/src/correlations.jl b/src/correlations.jl
index 5f3c6d03..38d7cb99 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -20,7 +20,7 @@ ExponentialSeries(; tol = 1e-14, calc_steadystate = false) = ExponentialSeries(t
         A::QuantumObject,
         B::QuantumObject,
         C::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         kwargs...)
 
 Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\expval{\hat{A}(t) \hat{B}(t + \tau) \hat{C}(t)}``
@@ -35,7 +35,7 @@ function correlation_3op_2t(
     A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
     C::QuantumObject{<:AbstractArray{T5},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     kwargs...,
 ) where {
     T1,
@@ -65,7 +65,7 @@ end
         τ_l::AbstractVector,
         A::QuantumObject,
         B::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         reverse::Bool=false,
         kwargs...)
 
@@ -81,7 +81,7 @@ function correlation_2op_2t(
     τ_l::AbstractVector,
     A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     reverse::Bool = false,
     kwargs...,
 ) where {
@@ -108,7 +108,7 @@ end
         τ_l::AbstractVector,
         A::QuantumObject,
         B::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         reverse::Bool=false,
         kwargs...)
 
@@ -122,7 +122,7 @@ function correlation_2op_1t(
     τ_l::AbstractVector,
     A::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T4},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     reverse::Bool = false,
     kwargs...,
 ) where {
@@ -143,7 +143,7 @@ end
         ω_list::AbstractVector,
         A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
         B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         solver::MySolver=ExponentialSeries(),
         kwargs...)
 
@@ -158,7 +158,7 @@ function spectrum(
     ω_list::AbstractVector,
     A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
     B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::MySolver = ExponentialSeries(),
     kwargs...,
 ) where {
diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl
index dee88be2..483b82a5 100644
--- a/src/qobj/eigsolve.jl
+++ b/src/qobj/eigsolve.jl
@@ -323,7 +323,7 @@ end
 
 @doc raw"""
     eigsolve_al(H::QuantumObject,
-        T::Real, c_ops::Union{Nothing,AbstractVector}=nothing;
+        T::Real, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
         H_t::Union{Nothing,Function}=nothing,
         params::NamedTuple=NamedTuple(),
@@ -363,7 +363,7 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 function eigsolve_al(
     H::QuantumObject{MT1,HOpType},
     T::Real,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     H_t::Union{Nothing,Function} = nothing,
     params::NamedTuple = NamedTuple(),
diff --git a/src/qobj/operator_sum.jl b/src/qobj/operator_sum.jl
index 909e27db..5b239821 100644
--- a/src/qobj/operator_sum.jl
+++ b/src/qobj/operator_sum.jl
@@ -7,10 +7,13 @@ A constructor to represent a sum of operators ``\sum_i c_i \hat{O}_i`` with a li
 
 This is very useful when we have to update only the coefficients, without allocating memory by performing the sum of the operators.
 """
-struct OperatorSum{CT<:Vector{<:Number},OT<:AbstractVector} <: AbstractQuantumObject
+struct OperatorSum{CT<:AbstractVector{<:Number},OT<:Union{AbstractVector,Tuple}} <: AbstractQuantumObject
     coefficients::CT
     operators::OT
-    function OperatorSum(coefficients::CT, operators::OT) where {CT<:Vector{<:Number},OT<:AbstractVector}
+    function OperatorSum(
+        coefficients::CT,
+        operators::OT,
+    ) where {CT<:AbstractVector{<:Number},OT<:Union{AbstractVector,Tuple}}
         length(coefficients) == length(operators) ||
             throw(DimensionMismatch("The number of coefficients must be the same as the number of operators."))
         # Check if all the operators have the same dimensions
@@ -22,7 +25,8 @@ struct OperatorSum{CT<:Vector{<:Number},OT<:AbstractVector} <: AbstractQuantumOb
             mapreduce(eltype, promote_type, coefficients),
         )
         coefficients2 = T.(coefficients)
-        return new{Vector{T},OT}(coefficients2, operators)
+        CT2 = typeof(coefficients2)
+        return new{CT2,OT}(coefficients2, operators)
     end
 end
 
diff --git a/src/steadystate.jl b/src/steadystate.jl
index ab60695e..59e6f6ad 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -66,7 +66,7 @@ end
 @doc raw"""
     steadystate(
         H::QuantumObject,
-        c_ops::Union{Nothing,AbstractVector} = nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateSolver = SteadyStateDirectSolver(),
         kwargs...
     )
@@ -75,13 +75,13 @@ Solve the stationary state based on different solvers.
 
 # Parameters
 - `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
-- `c_ops::Union{Nothing,AbstractVector}=nothing`: The list of the collapse operators.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators.
 - `solver::SteadyStateSolver=SteadyStateDirectSolver()`: see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for different solvers.
 - `kwargs...`: The keyword arguments for the solver.
 """
 function steadystate(
     H::QuantumObject{<:AbstractArray,OpType},
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
 ) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
@@ -188,7 +188,7 @@ _steadystate(
         H::QuantumObject,
         ψ0::QuantumObject,
         tspan::Real = Inf,
-        c_ops::Union{Nothing,AbstractVector} = nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateODESolver = SteadyStateODESolver(),
         reltol::Real = 1.0e-8,
         abstol::Real = 1.0e-10,
@@ -213,7 +213,7 @@ or
 - `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
 - `ψ0::QuantumObject`: The initial state of the system.
 - `tspan::Real=Inf`: The final time step for the steady state problem.
-- `c_ops::Union{Nothing,AbstractVector}=nothing`: The list of the collapse operators.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators.
 - `solver::SteadyStateODESolver=SteadyStateODESolver()`: see [`SteadyStateODESolver`](@ref) for more details.
 - `reltol::Real=1.0e-8`: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
 - `abstol::Real=1.0e-10`: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
@@ -223,7 +223,7 @@ function steadystate(
     H::QuantumObject{MT1,HOpType},
     ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
     tspan::Real = Inf,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateODESolver = SteadyStateODESolver(),
     reltol::Real = 1.0e-8,
     abstol::Real = 1.0e-10,
@@ -274,7 +274,7 @@ end
         H_p::QuantumObject{<:AbstractArray,OpType2},
         H_m::QuantumObject{<:AbstractArray,OpType3},
         ωd::Number,
-        c_ops::Union{Nothing,AbstractVector} = nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         n_max::Integer = 2,
         tol::R = 1e-8,
         solver::FSolver = SSFloquetLinearSystem,
@@ -341,7 +341,7 @@ In the case of `SSFloquetEffectiveLiouvillian`, instead, the effective Liouvilli
 - `H_p::QuantumObject`: The Hamiltonian or the Liouvillian of the part of the drive that oscillates as ``e^{i \omega t}``.
 - `H_m::QuantumObject`: The Hamiltonian or the Liouvillian of the part of the drive that oscillates as ``e^{-i \omega t}``.
 - `ωd::Number`: The frequency of the drive.
-- `c_ops::AbstractVector = QuantumObject`: The optional collapse operators.
+- `c_ops::Union{Nothing,AbstractVector} = nothing`: The optional collapse operators.
 - `n_max::Integer = 2`: The number of Fourier components to consider.
 - `tol::R = 1e-8`: The tolerance for the solver.
 - `solver::FSolver = SSFloquetLinearSystem`: The solver to use.
@@ -352,7 +352,7 @@ function steadystate_floquet(
     H_p::QuantumObject{<:AbstractArray,OpType2},
     H_m::QuantumObject{<:AbstractArray,OpType3},
     ωd::Number,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     n_max::Integer = 2,
     tol::R = 1e-8,
     solver::FSolver = SSFloquetLinearSystem(),
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index cfebcbf9..f73c1edc 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -100,9 +100,9 @@ end
     mcsolveProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
         ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
         tlist::AbstractVector,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         jump_callback::TJC=ContinuousLindbladJumpCallback(),
@@ -147,9 +147,9 @@ If the environmental measurements register a quantum jump, the wave function und
 - `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
 - `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -172,9 +172,9 @@ function mcsolveProblem(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray,KetQuantumObject},
     tlist::AbstractVector,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     seeds::Union{Nothing,Vector{Int}} = nothing,
@@ -286,9 +286,9 @@ end
     mcsolveEnsembleProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
         ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
         tlist::AbstractVector,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         jump_callback::TJC=ContinuousLindbladJumpCallback(),
@@ -335,9 +335,9 @@ If the environmental measurements register a quantum jump, the wave function und
 - `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
 - `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -362,9 +362,9 @@ function mcsolveEnsembleProblem(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
     tlist::AbstractVector,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     jump_callback::TJC = ContinuousLindbladJumpCallback(),
@@ -396,9 +396,9 @@ end
     mcsolve(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
         ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
         tlist::AbstractVector,
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         ntraj::Int=1,
@@ -445,9 +445,9 @@ If the environmental measurements register a quantum jump, the wave function und
 - `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `c_ops::Union{Nothing,AbstractVector}`: List of collapse operators ``\{\hat{C}_n\}_n``.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}`: List of collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
 - `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -475,9 +475,9 @@ function mcsolve(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
     tlist::AbstractVector,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     seeds::Union{Nothing,Vector{Int}} = nothing,
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 922f2994..8dc5895d 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -52,9 +52,9 @@ end
     mesolveProblem(H::QuantumObject,
         ψ0::QuantumObject,
         tlist::AbstractVector, 
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         progress_bar::Union{Val,Bool}=Val(true),
@@ -77,9 +77,9 @@ where
 - `H::QuantumObject`: The Hamiltonian ``\hat{H}`` or the Liouvillian of the system.
 - `ψ0::QuantumObject`: The initial state of the system.
 - `tlist::AbstractVector`: The time list of the evolution.
-- `c_ops::Union{Nothing,AbstractVector}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::OrdinaryDiffEqAlgorithm=Tsit5()`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}=nothing`: The list of the operators for which the expectation values are calculated.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the operators for which the expectation values are calculated.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing`: The time-dependent Hamiltonian or Liouvillian.
 - `params::NamedTuple=NamedTuple()`: The parameters of the time evolution.
 - `progress_bar::Union{Val,Bool}=Val(true)`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -101,9 +101,9 @@ function mesolveProblem(
     H::QuantumObject{MT1,HOpType},
     ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
     tlist,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -170,9 +170,9 @@ end
     mesolve(H::QuantumObject,
         ψ0::QuantumObject,
         tlist::AbstractVector, 
-        c_ops::Union{Nothing,AbstractVector}=nothing;
+        c_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         progress_bar::Union{Val,Bool}=Val(true),
@@ -195,9 +195,9 @@ where
 - `H::QuantumObject`: The Hamiltonian ``\hat{H}`` or the Liouvillian of the system.
 - `ψ0::QuantumObject`: The initial state of the system.
 - `tlist::AbstractVector`: The time list of the evolution.
-- `c_ops::Union{Nothing,AbstractVector}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
+- `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Named Tuple of parameters to pass to the solver.
 - `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -219,9 +219,9 @@ function mesolve(
     H::QuantumObject{MT1,HOpType},
     ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
     tlist::AbstractVector,
-    c_ops::Union{Nothing,AbstractVector} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 383ff2a0..42966e12 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -48,7 +48,7 @@ end
         ψ0::QuantumObject,
         tlist::AbstractVector;
         alg::OrdinaryDiffEqAlgorithm=Tsit5()
-        e_ops::Union{Nothing,AbstractVector} = nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         progress_bar::Union{Val,Bool}=Val(true),
@@ -66,7 +66,7 @@ Generates the ODEProblem for the Schrödinger time evolution of a quantum system
 - `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: The time list of the evolution.
 - `alg::OrdinaryDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: The list of operators to be evaluated during the evolution.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: The list of operators to be evaluated during the evolution.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
 - `params::NamedTuple`: The parameters of the system.
 - `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -89,7 +89,7 @@ function sesolveProblem(
     ψ0::QuantumObject{<:AbstractVector{T2},KetQuantumObject},
     tlist::AbstractVector;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -148,7 +148,7 @@ end
         ψ0::QuantumObject,
         tlist::AbstractVector;
         alg::OrdinaryDiffEqAlgorithm=Tsit5(),
-        e_ops::Union{Nothing,AbstractVector} = nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         progress_bar::Union{Val,Bool}=Val(true),
@@ -166,7 +166,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 - `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: List of times at which to save the state of the system.
 - `alg::OrdinaryDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
 - `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -189,7 +189,7 @@ function sesolve(
     ψ0::QuantumObject{<:AbstractVector{T2},KetQuantumObject},
     tlist::AbstractVector;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 9d77861b..a423e66d 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -67,9 +67,9 @@ end
     ssesolveProblem(H::QuantumObject,
         ψ0::QuantumObject,
         tlist::AbstractVector;
-        sc_ops::Union{Nothing,AbstractVector}=nothing;
+        sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing;
         alg::StochasticDiffEqAlgorithm=SRA1()
-        e_ops::Union{Nothing,AbstractVector} = nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         progress_bar::Union{Val,Bool}=Val(true),
@@ -101,9 +101,9 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 - `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: The time list of the evolution.
-- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}=nothing`: The list of operators to be evaluated during the evolution.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of operators to be evaluated during the evolution.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
 - `params::NamedTuple`: The parameters of the system.
 - `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -125,9 +125,9 @@ function ssesolveProblem(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
     tlist::AbstractVector,
-    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::StochasticDiffEqAlgorithm = SRA1(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -209,9 +209,9 @@ end
     ssesolveEnsembleProblem(H::QuantumObject,
         ψ0::QuantumObject,
         tlist::AbstractVector;
-        sc_ops::Union{Nothing,AbstractVector} = nothing;
+        sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         alg::StochasticDiffEqAlgorithm=SRA1()
-        e_ops::Union{Nothing,AbstractVector} = nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         prob_func::Function=_mcsolve_prob_func,
@@ -244,9 +244,9 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 - `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: The time list of the evolution.
-- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}=nothing`: The list of operators to be evaluated during the evolution.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of operators to be evaluated during the evolution.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
 - `params::NamedTuple`: The parameters of the system.
 - `prob_func::Function`: Function to use for generating the SDEProblem.
@@ -269,9 +269,9 @@ function ssesolveEnsembleProblem(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
     tlist::AbstractVector,
-    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::StochasticDiffEqAlgorithm = SRA1(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     prob_func::Function = _ssesolve_prob_func,
@@ -291,7 +291,7 @@ end
         tlist::AbstractVector,
         sc_ops::Union{Nothing, AbstractVector}=nothing;
         alg::StochasticDiffEqAlgorithm=SRA1(),
-        e_ops::Union{Nothing,AbstractVector}=nothing,
+        e_ops::Union{Nothing,AbstractVector,Tuple}=nothing,
         H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
         params::NamedTuple=NamedTuple(),
         ntraj::Int=1,
@@ -329,9 +329,9 @@ Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener i
 - `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
 - `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
 - `tlist::AbstractVector`: List of times at which to save the state of the system.
-- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `sc_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
 - `alg::StochasticDiffEqAlgorithm`: Algorithm to use for the time evolution.
-- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `e_ops::Union{Nothing,AbstractVector,Tuple}`: List of operators for which to calculate expectation values.
 - `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
 - `params::NamedTuple`: Dictionary of parameters to pass to the solver.
 - `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
@@ -358,9 +358,9 @@ function ssesolve(
     H::QuantumObject{MT1,OperatorQuantumObject},
     ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
     tlist::AbstractVector,
-    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::StochasticDiffEqAlgorithm = SRA1(),
-    e_ops::Union{Nothing,AbstractVector} = nothing,
+    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
     params::NamedTuple = NamedTuple(),
     ntraj::Int = 1,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 31e4d2ab..f3878295 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -164,7 +164,7 @@ end
 
 function TimeDependentOperatorSum(
     coefficient_functions,
-    operators::Vector{<:QuantumObject};
+    operators::Union{AbstractVector{<:QuantumObject},Tuple};
     params = nothing,
     init_time = 0.0,
 )
@@ -198,7 +198,7 @@ end
 
 ### LIOUVILLIAN ###
 @doc raw"""
-    liouvillian(H::QuantumObject, c_ops::Union{AbstractVector,Nothing}=nothing, Id_cache=I(prod(H.dims)))
+    liouvillian(H::QuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dims)))
 
 Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
 
@@ -218,7 +218,7 @@ See also [`spre`](@ref), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
 """
 function liouvillian(
     H::QuantumObject{MT1,OpType1},
-    c_ops::Union{AbstractVector,Nothing} = nothing,
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     Id_cache = I(prod(H.dims)),
 ) where {MT1<:AbstractMatrix,OpType1<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
@@ -254,7 +254,7 @@ function liouvillian_floquet(
     Hₚ::QuantumObject{<:AbstractArray{T2},OpType2},
     Hₘ::QuantumObject{<:AbstractArray{T3},OpType3},
     ω::Real,
-    c_ops::Union{AbstractVector,Nothing} = nothing;
+    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     n_max::Int = 3,
     tol::Real = 1e-15,
 ) where {
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index ed850565..1b3df4e7 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -64,7 +64,7 @@
     e_ops = [a_d * a]
     psi0 = fock(N, 3)
     t_l = LinRange(0, 100 * 2π, 1000)
-    H_t_f = TimeDependentOperatorSum((((t, p) -> sin(t)),), [H_t]) # It will be converted to liouvillian internally
+    H_t_f = TimeDependentOperatorSum((((t, p) -> sin(t)),), (H_t,)) # It will be converted to liouvillian internally
     sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, H_t = H_t_f, progress_bar = Val(false))
     ρ_ss1 = steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetLinearSystem())[1]
     ρ_ss2 =
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 9c02da21..7b1f4b42 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -42,6 +42,7 @@
             @inferred sesolve(H, psi0, t_l, e_ops = e_ops, progress_bar = Val(false))
             @inferred sesolve(H, psi0, t_l, progress_bar = Val(false))
             @inferred sesolve(H, psi0, t_l, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
+            @inferred sesolve(H, psi0, t_l, e_ops = (a_d * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple of different types
         end
     end
 
@@ -123,6 +124,7 @@
             @inferred mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolve(H, psi0, t_l, c_ops, progress_bar = Val(false))
             @inferred mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
+            @inferred mesolve(H, psi0, t_l, (a, a'), e_ops = (a_d * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple of different types
         end
 
         @testset "Type Inference mcsolve" begin
@@ -131,6 +133,7 @@
             @inferred mcsolve(H, psi0, t_l, c_ops, ntraj = 500, progress_bar = Val(true))
             @inferred mcsolve(H, psi0, [0, 10], c_ops, ntraj = 500, progress_bar = Val(false))
             @inferred mcsolve(H, Qobj(zeros(Int64, N)), t_l, c_ops, ntraj = 500, progress_bar = Val(false))
+            @inferred mcsolve(H, psi0, t_l, (a, a'), e_ops = (a_d * a, a'), ntraj = 500, progress_bar = Val(false)) # We test the type inference for Tuple of different types
         end
 
         @testset "Type Inference ssesolve" begin