diff --git a/.gitignore b/.gitignore
index b07886a94..1d017f123 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,7 @@ Manifest.toml
 benchmarks/benchmarks_output.json
 
 .ipynb_checkpoints
-*.ipynb
\ No newline at end of file
+*.ipynb
+
+
+.devcontainer/*
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 84bc5139e..ca788a467 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 - Add error message for bad input in state/operator generating functions ([#603])
+- Extend `sesolve` and `mesolve` to handle `Operator` and `SuperOperator` as initial conditions for propagator calculation. This introduces a `states_type` field to `TimeEvolutionProblem`. ([#606])
 
 ## [v0.39.1]
 Release date: 2025-11-19
@@ -383,3 +384,4 @@ Release date: 2024-11-13
 [#591]: https://github.com/qutip/QuantumToolbox.jl/issues/591
 [#596]: https://github.com/qutip/QuantumToolbox.jl/issues/596
 [#603]: https://github.com/qutip/QuantumToolbox.jl/issues/603
+[#606]: https://github.com/qutip/QuantumToolbox.jl/issues/606
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 1d2cf5689..219bd5854 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -119,10 +119,12 @@ Converts a sparse QuantumObject to a dense QuantumObject.
 to_dense(A::QuantumObject) = QuantumObject(to_dense(A.data), A.type, A.dimensions)
 to_dense(A::MT) where {MT<:AbstractSparseArray} = Array(A)
 to_dense(A::MT) where {MT<:AbstractArray} = A
+to_dense(A::Diagonal) = diagm(A.diag)
 
 to_dense(::Type{T}, A::AbstractSparseArray) where {T<:Number} = Array{T}(A)
 to_dense(::Type{T1}, A::AbstractArray{T2}) where {T1<:Number,T2<:Number} = Array{T1}(A)
 to_dense(::Type{T}, A::AbstractArray{T}) where {T<:Number} = A
+to_dense(::Type{T}, A::Diagonal{T}) where {T<:Number} = diagm(A.diag)
 
 function to_dense(::Type{M}) where {M<:Union{Diagonal,SparseMatrixCSC}}
     T = M
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 8630810c7..7a2c22449 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -261,6 +261,7 @@ function mcsolveEnsembleProblem(
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob_mc.times,
+        prob_mc.states_type,
         prob_mc.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index b9aee6595..e55a655af 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -6,17 +6,20 @@ _mesolve_make_L_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}, c_ops) = liouvil
 _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."))
 
-function _gen_mesolve_solution(sol, times, dimensions, isoperket::Val)
-    if getVal(isoperket)
-        ρt = map(ϕ -> QuantumObject(ϕ, type = OperatorKet(), dims = dimensions), sol.u)
+function _gen_mesolve_solution(
+    sol,
+    prob::TimeEvolutionProblem{ST},
+) where {ST<:Union{Operator,OperatorKet,SuperOperator}}
+    if prob.states_type isa Operator
+        ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = prob.states_type, dims = prob.dimensions), sol.u)
     else
-        ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator(), dims = dimensions), sol.u)
+        ρt = map(ϕ -> QuantumObject(ϕ, type = prob.states_type, dims = prob.dimensions), sol.u)
     end
 
     kwargs = NamedTuple(sol.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility
 
     return TimeEvolutionSol(
-        times,
+        prob.times,
         sol.t,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
@@ -55,7 +58,7 @@ where
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref), [`OperatorKet`](@ref), or [`SuperOperator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -66,6 +69,7 @@ where
 
 # Notes
 
+- The initial state `ψ0` can also be [`SuperOperator`](@ref). This is useful for simulating many density matrices simultaneously or calculating propagator. For example, when `H` is a [`SuperOperator`](@ref), `ψ0` can be given as `qeye_like(H)` (an identity [`SuperOperator`](@ref) matrix).
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolveProblem`](@ref) instead.
@@ -86,7 +90,7 @@ function mesolveProblem(
     progress_bar::Union{Val,Bool} = Val(true),
     inplace::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}}
+) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet,SuperOperator}}
     (isoper(H) && isket(ψ0) && isnothing(c_ops)) && return sesolveProblem(
         H,
         ψ0,
@@ -106,12 +110,18 @@ function mesolveProblem(
     L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
     check_dimensions(L_evo, ψ0)
 
-    T = Base.promote_eltype(L_evo, ψ0)
-    ρ0 = if isoperket(ψ0) # Convert it to dense vector with complex element type
-        to_dense(_complex_float_type(T), copy(ψ0.data))
-    else
-        to_dense(_complex_float_type(T), mat2vec(ket2dm(ψ0).data))
+    # Convert to dense vector with complex element type
+
+    T = _complex_float_type(Base.promote_eltype(L_evo, ψ0))
+    is_operket_or_super = isoperket(ψ0) || issuper(ψ0)
+    ρ0 = is_operket_or_super ? to_dense(T, copy(ψ0.data)) : to_dense(T, mat2vec(ket2dm(ψ0).data))
+    states_type = ψ0.type
+    if !is_operket_or_super
+        states_type = Operator()
     end
+
+
+
     L = cache_operator(L_evo.data, ρ0)
 
     kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
@@ -122,7 +132,7 @@ function mesolveProblem(
 
     prob = ODEProblem{getVal(inplace),FullSpecialize}(L, ρ0, tspan, params; kwargs4...)
 
-    return TimeEvolutionProblem(prob, tlist, L_evo.dimensions, (isoperket = Val(isoperket(ψ0)),))
+    return TimeEvolutionProblem(prob, tlist, states_type, L_evo.dimensions)
 end
 
 @doc raw"""
@@ -154,7 +164,7 @@ where
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref), [`OperatorKet`](@ref), or [`SuperOperator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm for the ODE solver. The default value is `DP5()`.
@@ -165,7 +175,7 @@ where
 - `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
-
+- The initial state `ψ0` can also be [`SuperOperator`](@ref). This is useful for simulating many density matrices simultaneously or calculating propagator. For example, when `H` is a [`SuperOperator`](@ref), `ψ0` can be given as `qeye_like(H)` (an identity [`SuperOperator`](@ref) matrix).
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolve`](@ref) instead.
@@ -188,7 +198,7 @@ function mesolve(
     progress_bar::Union{Val,Bool} = Val(true),
     inplace::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}}
+) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet,SuperOperator}}
     (isoper(H) && isket(ψ0) && isnothing(c_ops)) && return sesolve(
         H,
         ψ0,
@@ -230,7 +240,7 @@ end
 function mesolve(prob::TimeEvolutionProblem, alg::AbstractODEAlgorithm = DP5(); kwargs...)
     sol = solve(prob.prob, alg; kwargs...)
 
-    return _gen_mesolve_solution(sol, prob.times, prob.dimensions, prob.kwargs.isoperket)
+    return _gen_mesolve_solution(sol, prob)
 end
 
 @doc raw"""
@@ -266,7 +276,7 @@ for each combination in the ensemble.
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
+- `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states. It can be either a [`Ket`](@ref), [`Operator`](@ref), [`OperatorKet`](@ref), or [`SuperOperator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm for the ODE solver. The default is `DP5()`.
@@ -278,6 +288,7 @@ for each combination in the ensemble.
 
 # Notes
 
+- The initial state `ψ0` can also be [`SuperOperator`](@ref). This is useful for simulating many density matrices simultaneously or calculating propagator. For example, when `H` is a [`SuperOperator`](@ref), `ψ0` can be given as `qeye_like(H)` (an identity [`SuperOperator`](@ref) matrix).
 - The function returns an array of solutions with dimensions matching the Cartesian product of initial states and parameter sets.
 - If `ψ0` is a vector of `m` states and `params = (p1, p2, ...)` where `p1` has length `n1`, `p2` has length `n2`, etc., the output will be of size `(m, n1, n2, ...)`.
 - If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolve_map`](@ref) instead.
@@ -298,7 +309,7 @@ function mesolve_map(
     params::Union{NullParameters,Tuple} = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}}
+) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet,SuperOperator}}
     (isoper(H) && all(isket, ψ0) && isnothing(c_ops)) && return sesolve_map(
         H,
         ψ0,
@@ -315,7 +326,7 @@ function mesolve_map(
     # Convert to appropriate format based on state type
     ψ0_iter = map(ψ0) do state
         T = _complex_float_type(eltype(state))
-        if isoperket(state)
+        if isoperket(state) || issuper(state)
             to_dense(T, copy(state.data))
         else
             to_dense(T, mat2vec(ket2dm(state).data))
@@ -347,7 +358,7 @@ mesolve_map(
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     kwargs...,
-) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}} =
+) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet,SuperOperator}} =
     mesolve_map(H, [ψ0], tlist, c_ops; kwargs...)
 
 # this method is for advanced usage
@@ -357,14 +368,14 @@ mesolve_map(
 #
 # Return: An array of TimeEvolutionSol objects with the size same as the given iter.
 function mesolve_map(
-    prob::TimeEvolutionProblem{<:ODEProblem},
+    prob::TimeEvolutionProblem{StateOpType,<:AbstractDimensions,<:ODEProblem},
     iter::AbstractArray,
     alg::AbstractODEAlgorithm = DP5(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
     progress_bar::Union{Val,Bool} = Val(true),
-)
+) where {StateOpType<:Union{Ket,Operator,OperatorKet,SuperOperator}}
     # generate ensemble problem
     ntraj = length(iter)
     _prob_func = isnothing(prob_func) ? (prob, i, repeat) -> _se_me_map_prob_func(prob, i, repeat, iter) : prob_func
@@ -380,14 +391,14 @@ function mesolve_map(
     ens_prob = TimeEvolutionProblem(
         EnsembleProblem(prob.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob.times,
+        prob.states_type,
         prob.dimensions,
-        (progr = _output_func[2], channel = _output_func[3], isoperket = prob.kwargs.isoperket),
+        (progr = _output_func[2], channel = _output_func[3]),
     )
 
     sol = _ensemble_dispatch_solve(ens_prob, alg, ensemblealg, ntraj)
 
     # handle solution and make it become an Array of TimeEvolutionSol
-    sol_vec =
-        [_gen_mesolve_solution(sol[:, i], prob.times, prob.dimensions, prob.kwargs.isoperket) for i in eachindex(sol)] # map is type unstable
+    sol_vec = [_gen_mesolve_solution(sol[:, i], prob) for i in eachindex(sol)] # map is type unstable
     return reshape(sol_vec, size(iter))
 end
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index df0b2cd93..8aaa13588 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -2,13 +2,13 @@ export sesolveProblem, sesolve, sesolve_map
 
 _sesolve_make_U_QobjEvo(H) = -1im * QuantumObjectEvolution(H, type = Operator())
 
-function _gen_sesolve_solution(sol, times, dimensions)
-    ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = dimensions), sol.u)
+function _gen_sesolve_solution(sol, prob::TimeEvolutionProblem{ST}) where {ST<:Union{Ket,Operator}}
+    ψt = map(ϕ -> QuantumObject(ϕ, type = prob.states_type, dims = prob.dimensions), sol.u)
 
     kwargs = NamedTuple(sol.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility
 
     return TimeEvolutionSol(
-        times,
+        prob.times,
         sol.t,
         ψt,
         _get_expvals(sol, SaveFuncSESolve),
@@ -40,7 +40,7 @@ Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref) or a [`Operator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -50,6 +50,7 @@ Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 
 # Notes
 
+- The initial state `ψ0` can also be [`Operator`](@ref). This is useful for simulating many states simultaneously or calculating propagator. For example, `ψ0` can be given as `qeye_like(H)` (an identity [`Operator`](@ref) matrix).
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
@@ -61,18 +62,19 @@ Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 """
 function sesolveProblem(
     H::Union{AbstractQuantumObject{Operator},Tuple},
-    ψ0::QuantumObject{Ket},
+    ψ0::QuantumObject{ST},
     tlist::AbstractVector;
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
     inplace::Union{Val,Bool} = Val(true),
     kwargs...,
-)
+) where {ST<:Union{Ket,Operator}}
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
     tlist = _check_tlist(tlist, _float_type(ψ0))
+    states_type = ψ0.type
 
     H_evo = _sesolve_make_U_QobjEvo(H) # Multiply by -i
     isoper(H_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
@@ -90,7 +92,7 @@ function sesolveProblem(
 
     prob = ODEProblem{getVal(inplace),FullSpecialize}(U, ψ0, tspan, params; kwargs4...)
 
-    return TimeEvolutionProblem(prob, tlist, H_evo.dimensions)
+    return TimeEvolutionProblem(prob, tlist, states_type, H_evo.dimensions)
 end
 
 @doc raw"""
@@ -115,7 +117,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref) or a [`Operator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `alg`: The algorithm for the ODE solver. The default is `Vern7(lazy=false)`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -126,6 +128,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 
 # Notes
 
+- The initial state `ψ0` can also be [`Operator`](@ref). This is useful for simulating many states simultaneously or calculating propagator. For example, `ψ0` can be given as `qeye_like(H)` (an identity [`Operator`](@ref) matrix).
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
@@ -138,7 +141,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 """
 function sesolve(
     H::Union{AbstractQuantumObject{Operator},Tuple},
-    ψ0::QuantumObject{Ket},
+    ψ0::QuantumObject{ST},
     tlist::AbstractVector;
     alg::AbstractODEAlgorithm = Vern7(lazy = false),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
@@ -146,7 +149,7 @@ function sesolve(
     progress_bar::Union{Val,Bool} = Val(true),
     inplace::Union{Val,Bool} = Val(true),
     kwargs...,
-)
+) where {ST<:Union{Ket,Operator}}
 
     # Move sensealg argument to solve for Enzyme.jl support.
     # TODO: Remove it when https://github.com/SciML/SciMLSensitivity.jl/issues/1225 is fixed.
@@ -175,7 +178,7 @@ end
 function sesolve(prob::TimeEvolutionProblem, alg::AbstractODEAlgorithm = Vern7(lazy = false); kwargs...)
     sol = solve(prob.prob, alg; kwargs...)
 
-    return _gen_sesolve_solution(sol, prob.times, prob.dimensions)
+    return _gen_sesolve_solution(sol, prob)
 end
 
 @doc raw"""
@@ -204,7 +207,7 @@ for each combination in the ensemble.
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states.
+- `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states. It can be either a [`Ket`](@ref) or [`Operator`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `alg`: The algorithm for the ODE solver. The default is `Vern7(lazy=false)`.
 - `ensemblealg`: Ensemble algorithm to use for parallel computation. Default is `EnsembleThreads()`.
@@ -215,6 +218,7 @@ for each combination in the ensemble.
 
 # Notes
 
+- The initial state `ψ0` can also be [`Operator`](@ref). This is useful for simulating many states simultaneously or calculating propagator. For example, `ψ0` can be given as `qeye_like(H)` (an identity [`Operator`](@ref) matrix).
 - The function returns an array of solutions with dimensions matching the Cartesian product of initial states and parameter sets.
 - If `ψ0` is a vector of `m` states and `params = (p1, p2, ...)` where `p1` has length `n1`, `p2` has length `n2`, etc., the output will be of size `(m, n1, n2, ...)`.
 - See [`sesolve`](@ref) for more details.
@@ -225,7 +229,7 @@ for each combination in the ensemble.
 """
 function sesolve_map(
     H::Union{AbstractQuantumObject{Operator},Tuple},
-    ψ0::AbstractVector{<:QuantumObject{Ket}},
+    ψ0::AbstractVector{<:QuantumObject{ST}},
     tlist::AbstractVector;
     alg::AbstractODEAlgorithm = Vern7(lazy = false),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
@@ -233,9 +237,12 @@ function sesolve_map(
     params::Union{NullParameters,Tuple} = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
     kwargs...,
-)
+) where {ST<:Union{Ket,Operator}}
     # mapping initial states and parameters
-    ψ0_iter = map(get_data, ψ0)
+
+    ψ0 = map(to_dense, ψ0) # Convert all initial states to dense vectors
+
+    ψ0_iter = map(state -> to_dense(_complex_float_type(eltype(state)), copy(state.data)), ψ0)
     if params isa NullParameters
         iter = collect(Iterators.product(ψ0_iter, [params])) |> vec # convert nx1 Matrix into Vector
     else
@@ -255,8 +262,12 @@ function sesolve_map(
 
     return sesolve_map(prob, iter, alg, ensemblealg; progress_bar = progress_bar)
 end
-sesolve_map(H::Union{AbstractQuantumObject{Operator},Tuple}, ψ0::QuantumObject{Ket}, tlist::AbstractVector; kwargs...) =
-    sesolve_map(H, [ψ0], tlist; kwargs...)
+sesolve_map(
+    H::Union{AbstractQuantumObject{Operator},Tuple},
+    ψ0::QuantumObject{ST},
+    tlist::AbstractVector;
+    kwargs...,
+) where {ST<:Union{Ket,Operator}} = sesolve_map(H, [ψ0], tlist; kwargs...)
 
 # this method is for advanced usage
 # User can define their own iterator structure, prob_func and output_func
@@ -265,14 +276,14 @@ sesolve_map(H::Union{AbstractQuantumObject{Operator},Tuple}, ψ0::QuantumObject{
 #
 # Return: An array of TimeEvolutionSol objects with the size same as the given iter.
 function sesolve_map(
-    prob::TimeEvolutionProblem{<:ODEProblem},
+    prob::TimeEvolutionProblem{ST,<:AbstractDimensions,<:ODEProblem},
     iter::AbstractArray,
     alg::AbstractODEAlgorithm = Vern7(lazy = false),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
     progress_bar::Union{Val,Bool} = Val(true),
-)
+) where {ST<:Union{Ket,Operator}}
     # generate ensemble problem
     ntraj = length(iter)
     _prob_func = isnothing(prob_func) ? (prob, i, repeat) -> _se_me_map_prob_func(prob, i, repeat, iter) : prob_func
@@ -288,6 +299,7 @@ function sesolve_map(
     ens_prob = TimeEvolutionProblem(
         EnsembleProblem(prob.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob.times,
+        prob.states_type,
         prob.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
@@ -295,6 +307,6 @@ function sesolve_map(
     sol = _ensemble_dispatch_solve(ens_prob, alg, ensemblealg, ntraj)
 
     # handle solution and make it become an Array of TimeEvolutionSol
-    sol_vec = [_gen_sesolve_solution(sol[:, i], prob.times, prob.dimensions) for i in eachindex(sol)] # map is type unstable
+    sol_vec = [_gen_sesolve_solution(sol[:, i], prob) for i in eachindex(sol)] # map is type unstable
     return reshape(sol_vec, size(iter))
 end
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 1560de45f..2147f1a5f 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -100,11 +100,11 @@ function smesolveProblem(
     check_dimensions(L_evo, ψ0)
     dims = L_evo.dimensions
 
-    T = Base.promote_eltype(L_evo, ψ0)
+    T = _complex_float_type(Base.promote_eltype(L_evo, ψ0))
     ρ0 = if isoperket(ψ0) # Convert it to dense vector with complex element type
-        to_dense(_complex_float_type(T), copy(ψ0.data))
+        to_dense(T, copy(ψ0.data))
     else
-        to_dense(_complex_float_type(T), mat2vec(ket2dm(ψ0).data))
+        to_dense(T, mat2vec(ket2dm(ψ0).data))
     end
 
     sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list))
@@ -146,7 +146,7 @@ function smesolveProblem(
         kwargs4...,
     )
 
-    return TimeEvolutionProblem(prob, tlist, dims, (isoperket = Val(isoperket(ψ0)),))
+    return TimeEvolutionProblem(prob, tlist, ψ0.type, dims, (isoperket = Val(isoperket(ψ0)),))
 end
 
 @doc raw"""
@@ -274,6 +274,7 @@ function smesolveEnsembleProblem(
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
+        prob_sme.states_type,
         prob_sme.dimensions,
         merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3])),
     )
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index df3372858..7a45c625c 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -101,6 +101,7 @@ function ssesolveProblem(
     check_dimensions(H_eff_evo, ψ0)
     dims = H_eff_evo.dimensions
 
+    states_type = ψ0.type
     ψ0 = to_dense(_complex_float_type(ψ0), get_data(ψ0))
 
     sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list))
@@ -142,7 +143,7 @@ function ssesolveProblem(
         kwargs4...,
     )
 
-    return TimeEvolutionProblem(prob, tlist, dims)
+    return TimeEvolutionProblem(prob, tlist, states_type, dims)
 end
 
 @doc raw"""
@@ -268,6 +269,7 @@ function ssesolveEnsembleProblem(
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
+        prob_sme.states_type,
         prob_sme.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index e7c3791d3..7e35c781a 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -19,15 +19,23 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 
 - `prob::AbstractSciMLProblem`: The `ODEProblem` of the time evolution.
 - `times::AbstractVector`: The time list of the evolution.
+- `states_type::QuantumObjectType`: The type of the quantum states during the evolution (e.g., [`Ket`](@ref), [`Operator`](@ref), [`OperatorKet`](@ref), or [`SuperOperator`](@ref)).
 - `dimensions::AbstractDimensions`: The dimensions of the Hilbert space.
 - `kwargs::KWT`: Generic keyword arguments.
 
 !!! note "`dims` property"
     For a given `prob::TimeEvolutionProblem`, `prob.dims` or `getproperty(prob, :dims)` returns its `dimensions` in the type of integer-vector.
 """
-struct TimeEvolutionProblem{PT<:AbstractSciMLProblem,TT<:AbstractVector,DT<:AbstractDimensions,KWT}
+struct TimeEvolutionProblem{
+    ST<:QuantumObjectType,
+    DT<:AbstractDimensions,
+    PT<:AbstractSciMLProblem,
+    TT<:AbstractVector,
+    KWT,
+}
     prob::PT
     times::TT
+    states_type::ST
     dimensions::DT
     kwargs::KWT
 end
@@ -41,7 +49,7 @@ function Base.getproperty(prob::TimeEvolutionProblem, key::Symbol)
     end
 end
 
-TimeEvolutionProblem(prob, times, dims) = TimeEvolutionProblem(prob, times, dims, nothing)
+TimeEvolutionProblem(prob, times, states_type, dims) = TimeEvolutionProblem(prob, times, states_type, dims, nothing)
 
 @doc raw"""
     struct TimeEvolutionSol
diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl
new file mode 100644
index 000000000..b57379e7b
--- /dev/null
+++ b/test/core-test/propagator.jl
@@ -0,0 +1,20 @@
+@testitem "Propagator (by solvers)" begin
+    ϵ0 = 1.0 * 2π
+    Δ = 0.8 * 2π
+    H = (ϵ0/2) * sigmaz() + (Δ/2) * sigmax()
+    L = liouvillian(H)
+    ψ0 = basis(2, 0)
+    ρ0 = mat2vec(ket2dm(ψ0))
+
+    dt = π/5
+    tlist = 0:dt:(2π)
+    ψt = sesolve(H, ψ0, tlist; progress_bar = Val(false)).states[2:end] # ignore the initial state
+    ρt = mesolve(H, ρ0, tlist; progress_bar = Val(false)).states[2:end] # ignore the initial state
+    Prop_se = sesolve(H, qeye_like(H), [0, dt]; progress_bar = Val(false)).states[end]
+    Prop_me = mesolve(L, qeye_like(L), [0, dt]; progress_bar = Val(false)).states[end]
+
+    for n in 1:(length(tlist)-1)
+        @test isapprox(Prop_se^n * ψ0, ψt[n]; atol = 1e-5)
+        @test isapprox(Prop_me^n * ρ0, ρt[n]; atol = 1e-5)
+    end
+end