diff --git a/CHANGELOG.md b/CHANGELOG.md
index 208ea97e8..8671bd3c8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Improve efficiency of `bloch_redfield_tensor` by avoiding unnecessary conversions. ([#509])
 - Support `SciMLOperators v1.4+`. ([#470])
 - Fix compatibility with `Makie v0.24+`. ([#513])
+- Add `keep_runs_results` option for multi-trajectory solvers to align with `QuTiP`. ([#512])
+  - Breaking changes for multi-trajectory solutions:
+    - the original fields `expect` and `states` now store the results depend on keyword argument `keep_runs_results` (decide whether to store all trajectories results or not).
+    - remove field `average_expect`
+    - remove field `runs_expect`
+  - New statistical analysis functions:
+    - `average_states`
+    - `average_expect`
+    - `std_expect`
 
 ## [v0.33.0]
 Release date: 2025-07-22
@@ -273,4 +282,5 @@ Release date: 2024-11-13
 [#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506
 [#507]: https://github.com/qutip/QuantumToolbox.jl/issues/507
 [#509]: https://github.com/qutip/QuantumToolbox.jl/issues/509
+[#512]: https://github.com/qutip/QuantumToolbox.jl/issues/512
 [#513]: https://github.com/qutip/QuantumToolbox.jl/issues/513
diff --git a/Project.toml b/Project.toml
index 745abe013..81acac6c1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [weakdeps]
@@ -65,5 +66,6 @@ SciMLOperators = "1.4"
 SparseArrays = "1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
+Statistics = "1"
 StochasticDiffEq = "6"
 julia = "1.10"
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index 8c9b8c209..da800fc76 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -197,6 +197,9 @@ TimeEvolutionProblem
 TimeEvolutionSol
 TimeEvolutionMCSol
 TimeEvolutionStochasticSol
+average_states
+average_expect
+std_expect
 sesolveProblem
 mesolveProblem
 mcsolveProblem
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index 94609f18a..e57d30362 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -97,7 +97,47 @@ Some other solvers can have other output.
 
 ## [Multiple trajectories solution](@id doc-TE:Multiple-trajectories-solution)
 
-This part is still under construction, please read the docstrings for the following functions first:
+The solutions are different for solvers which compute multiple trajectories, such as the [`TimeEvolutionMCSol`](@ref) (Monte Carlo) or the [`TimeEvolutionStochasticSol`](@ref) (stochastic methods). The storage of expectation values and states depends on the keyword argument `keep_runs_results`, which determines whether the results of all trajectories are stored in the solution.
 
-- [`TimeEvolutionMCSol`](@ref)
-- [`TimeEvolutionStochasticSol`](@ref)
\ No newline at end of file
+When the keyword argument `keep_runs_results` is passed as `Val(false)` to a multi-trajectory solver, the `states` and `expect` fields store only the average results (averaged over all trajectories). The results can be accessed by the following index-order:
+
+```julia
+sol.states[time_idx]
+sol.expect[e_op,time_idx]
+```
+
+For example:
+
+```@example TE-solution
+tlist = LinRange(0, 1, 11)
+c_ops = (destroy(2),)
+e_ops = (num(2),)
+
+sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(false), progress_bar=Val(false))
+
+size(sol_mc1.expect)
+```
+
+If the keyword argument `keep_runs_results = Val(true)`, the results for each trajectory and each time are stored, and the index-order of the elements in fields `states` and `expect` are:
+
+
+```julia
+sol.states[trajectory,time_idx]
+sol.expect[e_op,trajectory,time_idx]
+```
+
+For example:
+
+```@example TE-solution
+sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=25, keep_runs_results=Val(true), progress_bar=Val(false))
+
+size(sol_mc2.expect)
+```
+
+We also provide the following functions for statistical analysis of multi-trajectory `sol`utions:
+
+| **Functions** | **Description** |
+|:------------|:----------------|
+| [`average_states(sol)`](@ref average_states) | Return the trajectory-averaged result states (as density [`Operator`](@ref)) at each time point. |
+| [`average_expect(sol)`](@ref average_expect) | Return the trajectory-averaged expectation values at each time point. |
+| [`std_expect(sol)`](@ref std_expect) | Return the trajectory-wise standard deviation of the expectation values at each time point. |
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index afdbad0ed..c3d9be754 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -5,6 +5,7 @@ using LinearAlgebra
 import LinearAlgebra: BlasInt, BlasFloat, checksquare
 import LinearAlgebra.LAPACK: hseqr!
 using SparseArrays
+import Statistics: mean, std
 
 # SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
 import SciMLBase:
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index fe8da0876..647be42d4 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -22,7 +22,7 @@ end
 
 function _normalize_state!(u, dims, normalize_states)
     getVal(normalize_states) && normalize!(u)
-    return QuantumObject(u, type = Ket(), dims = dims)
+    return QuantumObject(u, Ket(), dims)
 end
 
 function _mcsolve_make_Heff_QobjEvo(H::QuantumObject, c_ops)
@@ -278,6 +278,7 @@ end
         progress_bar::Union{Val,Bool} = Val(true),
         prob_func::Union{Function, Nothing} = nothing,
         output_func::Union{Tuple,Nothing} = nothing,
+        keep_runs_results::Union{Val,Bool} = Val(false),
         normalize_states::Union{Val,Bool} = Val(true),
         kwargs...,
     )
@@ -332,6 +333,7 @@ If the environmental measurements register a quantum jump, the wave function und
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `prob_func`: Function to use for generating the ODEProblem.
 - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `normalize_states`: Whether to normalize the states. Default to `Val(true)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
 
@@ -363,6 +365,7 @@ function mcsolve(
     progress_bar::Union{Val,Bool} = Val(true),
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
+    keep_runs_results::Union{Val,Bool} = Val(false),
     normalize_states::Union{Val,Bool} = Val(true),
     kwargs...,
 ) where {TJC<:LindbladJumpCallbackType}
@@ -384,7 +387,7 @@ function mcsolve(
         kwargs...,
     )
 
-    return mcsolve(ens_prob_mc, alg, ntraj, ensemblealg, normalize_states)
+    return mcsolve(ens_prob_mc, alg, ntraj, ensemblealg, makeVal(keep_runs_results), normalize_states)
 end
 
 function mcsolve(
@@ -392,6 +395,7 @@ function mcsolve(
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     ntraj::Int = 500,
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
+    keep_runs_results = Val(false),
     normalize_states = Val(true),
 )
     sol = _ensemble_dispatch_solve(ens_prob_mc, alg, ensemblealg, ntraj)
@@ -403,25 +407,24 @@ function mcsolve(
     _expvals_all =
         _expvals_sol_1 isa Nothing ? nothing : map(i -> _get_expvals(sol[:, i], SaveFuncMCSolve), eachindex(sol))
     expvals_all = _expvals_all isa Nothing ? nothing : stack(_expvals_all, dims = 2) # Stack on dimension 2 to align with QuTiP
-    states = map(i -> _normalize_state!.(sol[:, i].u, Ref(dims), normalize_states), eachindex(sol))
+
+    # stack to transform Vector{Vector{QuantumObject}} -> Matrix{QuantumObject}
+    states_all = stack(map(i -> _normalize_state!.(sol[:, i].u, Ref(dims), normalize_states), eachindex(sol)), dims = 1)
+
     col_times = map(i -> _mc_get_jump_callback(sol[:, i]).affect!.col_times, eachindex(sol))
     col_which = map(i -> _mc_get_jump_callback(sol[:, i]).affect!.col_which, eachindex(sol))
 
-    expvals = _expvals_sol_1 isa Nothing ? nothing : dropdims(sum(expvals_all, dims = 2), dims = 2) ./ length(sol)
-
     return TimeEvolutionMCSol(
         ntraj,
         ens_prob_mc.times,
         _sol_1.t,
-        states,
-        expvals,
-        expvals, # This is average_expect
-        expvals_all,
+        _store_multitraj_states(states_all, keep_runs_results),
+        _store_multitraj_expect(expvals_all, keep_runs_results),
         col_times,
         col_which,
         sol.converged,
         _sol_1.alg,
-        NamedTuple(_sol_1.prob.kwargs).abstol,
-        NamedTuple(_sol_1.prob.kwargs).reltol,
+        _sol_1.prob.kwargs[:abstol],
+        _sol_1.prob.kwargs[:reltol],
     )
 end
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 85d6ed613..9dfd22273 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -211,7 +211,7 @@ function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
         sol.alg,
-        NamedTuple(sol.prob.kwargs).abstol,
-        NamedTuple(sol.prob.kwargs).reltol,
+        sol.prob.kwargs[:abstol],
+        sol.prob.kwargs[:reltol],
     )
 end
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index ef9d98f9b..cd635a705 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -161,7 +161,7 @@ function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
         _get_expvals(sol, SaveFuncSESolve),
         sol.retcode,
         sol.alg,
-        NamedTuple(sol.prob.kwargs).abstol,
-        NamedTuple(sol.prob.kwargs).reltol,
+        sol.prob.kwargs[:abstol],
+        sol.prob.kwargs[:reltol],
     )
 end
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index aa6dbae9b..1d28d4e02 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -293,6 +293,7 @@ end
         prob_func::Union{Function, Nothing} = nothing,
         output_func::Union{Tuple,Nothing} = nothing,
         progress_bar::Union{Val,Bool} = Val(true),
+        keep_runs_results::Union{Val,Bool} = Val(false),
         store_measurement::Union{Val,Bool} = Val(false),
         kwargs...,
     )
@@ -333,6 +334,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 - `prob_func`: Function to use for generating the SDEProblem.
 - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `store_measurement`: Whether to store the measurement expectation values. Default is `Val(false)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
 
@@ -365,6 +367,7 @@ function smesolve(
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
     progress_bar::Union{Val,Bool} = Val(true),
+    keep_runs_results::Union{Val,Bool} = Val(false),
     store_measurement::Union{Val,Bool} = Val(false),
     kwargs...,
 ) where {StateOpType<:Union{Ket,Operator,OperatorKet}}
@@ -392,7 +395,7 @@ function smesolve(
         alg = sc_ops_isa_Qobj ? SRIW1() : SRA2()
     end
 
-    return smesolve(ensemble_prob, alg, ntraj, ensemblealg)
+    return smesolve(ensemble_prob, alg, ntraj, ensemblealg, makeVal(keep_runs_results))
 end
 
 function smesolve(
@@ -400,6 +403,7 @@ function smesolve(
     alg::StochasticDiffEqAlgorithm = SRA2(),
     ntraj::Int = 500,
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
+    keep_runs_results = Val(false),
 )
     sol = _ensemble_dispatch_solve(ens_prob, alg, ensemblealg, ntraj)
 
@@ -412,24 +416,22 @@ function smesolve(
         _expvals_sol_1 isa Nothing ? nothing : map(i -> _get_expvals(sol[:, i], SaveFuncMESolve), eachindex(sol))
     expvals_all = _expvals_all isa Nothing ? nothing : stack(_expvals_all, dims = 2) # Stack on dimension 2 to align with QuTiP
 
-    states = map(i -> _smesolve_generate_state.(sol[:, i].u, Ref(dims), ens_prob.kwargs.isoperket), eachindex(sol))
+    # stack to transform Vector{Vector{QuantumObject}} -> Matrix{QuantumObject}
+    states_all = stack(
+        map(i -> _smesolve_generate_state.(sol[:, i].u, Ref(dims), ens_prob.kwargs.isoperket), eachindex(sol)),
+        dims = 1,
+    )
 
     _m_expvals =
         _m_expvals_sol_1 isa Nothing ? nothing : map(i -> _get_m_expvals(sol[:, i], SaveFuncSMESolve), eachindex(sol))
-    m_expvals = _m_expvals isa Nothing ? nothing : stack(_m_expvals, dims = 2)
-
-    expvals =
-        _get_expvals(_sol_1, SaveFuncMESolve) isa Nothing ? nothing :
-        dropdims(sum(expvals_all, dims = 2), dims = 2) ./ length(sol)
+    m_expvals = _m_expvals isa Nothing ? nothing : stack(_m_expvals, dims = 2) # Stack on dimension 2 to align with QuTiP
 
     return TimeEvolutionStochasticSol(
         ntraj,
         ens_prob.times,
         _sol_1.t,
-        states,
-        expvals,
-        expvals, # This is average_expect
-        expvals_all,
+        _store_multitraj_states(states_all, keep_runs_results),
+        _store_multitraj_expect(expvals_all, keep_runs_results),
         m_expvals, # Measurement expectation values
         sol.converged,
         _sol_1.alg,
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index e68987c02..487a3a06a 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -285,6 +285,7 @@ end
         prob_func::Union{Function, Nothing} = nothing,
         output_func::Union{Tuple,Nothing} = nothing,
         progress_bar::Union{Val,Bool} = Val(true),
+        keep_runs_results::Union{Val,Bool} = Val(false),
         store_measurement::Union{Val,Bool} = Val(false),
         kwargs...,
     )
@@ -328,6 +329,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 - `prob_func`: Function to use for generating the SDEProblem.
 - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `store_measurement`: Whether to store the measurement results. Default is `Val(false)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
 
@@ -360,6 +362,7 @@ function ssesolve(
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
     progress_bar::Union{Val,Bool} = Val(true),
+    keep_runs_results::Union{Val,Bool} = Val(false),
     store_measurement::Union{Val,Bool} = Val(false),
     kwargs...,
 )
@@ -386,7 +389,7 @@ function ssesolve(
         alg = sc_ops_isa_Qobj ? SRIW1() : SRA2()
     end
 
-    return ssesolve(ens_prob, alg, ntraj, ensemblealg)
+    return ssesolve(ens_prob, alg, ntraj, ensemblealg, makeVal(keep_runs_results))
 end
 
 function ssesolve(
@@ -394,6 +397,7 @@ function ssesolve(
     alg::StochasticDiffEqAlgorithm = SRA2(),
     ntraj::Int = 500,
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
+    keep_runs_results = Val(false),
 )
     sol = _ensemble_dispatch_solve(ens_prob, alg, ensemblealg, ntraj)
 
@@ -406,24 +410,20 @@ function ssesolve(
     _expvals_all =
         _expvals_sol_1 isa Nothing ? nothing : map(i -> _get_expvals(sol[:, i], SaveFuncSSESolve), eachindex(sol))
     expvals_all = _expvals_all isa Nothing ? nothing : stack(_expvals_all, dims = 2) # Stack on dimension 2 to align with QuTiP
-    states = map(i -> _normalize_state!.(sol[:, i].u, Ref(dims), normalize_states), eachindex(sol))
+
+    # stack to transform Vector{Vector{QuantumObject}} -> Matrix{QuantumObject}
+    states_all = stack(map(i -> _normalize_state!.(sol[:, i].u, Ref(dims), normalize_states), eachindex(sol)), dims = 1)
 
     _m_expvals =
         _m_expvals_sol_1 isa Nothing ? nothing : map(i -> _get_m_expvals(sol[:, i], SaveFuncSSESolve), eachindex(sol))
     m_expvals = _m_expvals isa Nothing ? nothing : stack(_m_expvals, dims = 2)
 
-    expvals =
-        _get_expvals(_sol_1, SaveFuncSSESolve) isa Nothing ? nothing :
-        dropdims(sum(expvals_all, dims = 2), dims = 2) ./ length(sol)
-
     return TimeEvolutionStochasticSol(
         ntraj,
         ens_prob.times,
         _sol_1.t,
-        states,
-        expvals,
-        expvals, # This is average_expect
-        expvals_all,
+        _store_multitraj_states(states_all, keep_runs_results),
+        _store_multitraj_expect(expvals_all, keep_runs_results),
         m_expvals, # Measurement expectation values
         sol.converged,
         _sol_1.alg,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index aa68cfc00..d04c7338b 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -1,4 +1,6 @@
-export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionStochasticSol
+export TimeEvolutionSol
+export TimeEvolutionMultiTrajSol, TimeEvolutionMCSol, TimeEvolutionStochasticSol
+export average_states, average_expect, std_expect
 
 export liouvillian_floquet, liouvillian_generalized
 
@@ -6,6 +8,8 @@ const DEFAULT_ODE_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep
 const DEFAULT_SDE_SOLVER_OPTIONS = (abstol = 1e-3, reltol = 2e-3, save_everystep = false, save_end = true)
 const COL_TIMES_WHICH_INIT_SIZE = 200
 
+abstract type TimeEvolutionMultiTrajSol{Tstates,Texpect} end
+
 @doc raw"""
     struct TimeEvolutionProblem
 
@@ -92,7 +96,7 @@ function Base.show(io::IO, sol::TimeEvolutionSol)
 end
 
 @doc raw"""
-    struct TimeEvolutionMCSol
+    struct TimeEvolutionMCSol <: TimeEvolutionMultiTrajSol
 
 A structure storing the results and some information from solving quantum trajectories of the Monte Carlo wave function time evolution.
 
@@ -101,36 +105,49 @@ A structure storing the results and some information from solving quantum trajec
 - `ntraj::Int`: Number of trajectories
 - `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times_states`.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
-- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
-- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`
+- `states::AbstractVecOrMat`: The list of result states in each trajectory and each time point in `times_states`.
+- `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
 - `col_times::Vector{Vector{Real}}`: The time records of every quantum jump occurred in each trajectory.
 - `col_which::Vector{Vector{Int}}`: The indices of which collapse operator was responsible for each quantum jump in `col_times`.
 - `converged::Bool`: Whether the solution is converged or not.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
 - `reltol::Real`: The relative tolerance which is used during the solving process.
+
+# Notes
+
+When the keyword argument `keep_runs_results` is passed as `Val(false)` to a multi-trajectory solver, the `states` and `expect` fields store only the average results (averaged over all trajectories). The results can be accessed by the following index-order:
+
+- `sol.states[time_idx]`
+- `sol.expect[e_op,time_idx]`
+
+If the keyword argument `keep_runs_results = Val(true)`, the results for each trajectory and each time are stored, and the index-order of the elements in fields `states` and `expect` are:
+
+- `sol.states[trajectory,time_idx]`
+- `sol.expect[e_op,trajectory,time_idx]`
+
+We also provide the following functions for statistical analysis of multi-trajectory solutions:
+
+- [`average_states`](@ref)
+- [`average_expect`](@ref)
+- [`std_expect`](@ref)
 """
 struct TimeEvolutionMCSol{
     TT1<:AbstractVector{<:Real},
     TT2<:AbstractVector{<:Real},
-    TS<:AbstractVector,
-    TE<:Union{AbstractMatrix,Nothing},
-    TEA<:Union{AbstractArray,Nothing},
+    TS<:AbstractVecOrMat,
+    TE<:Union{AbstractArray,Nothing},
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
     AlgT<:OrdinaryDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-}
+} <: TimeEvolutionMultiTrajSol{TS,TE}
     ntraj::Int
     times::TT1
     times_states::TT2
     states::TS
     expect::TE
-    average_expect::TE # Currently just a synonym for `expect`
-    runs_expect::TEA
     col_times::TJT
     col_which::TJW
     converged::Bool
@@ -144,11 +161,11 @@ function Base.show(io::IO, sol::TimeEvolutionMCSol)
     print(io, "(converged: $(sol.converged))\n")
     print(io, "--------------------------------\n")
     print(io, "num_trajectories = $(sol.ntraj)\n")
-    print(io, "num_states = $(length(sol.states[1]))\n")
+    print(io, "num_states = $(size(sol.states, ndims(sol.states)))\n") # get the size of last dimension
     if sol.expect isa Nothing
         print(io, "num_expect = 0\n")
     else
-        print(io, "num_expect = $(size(sol.average_expect, 1))\n")
+        print(io, "num_expect = $(size(sol.expect, 1))\n")
     end
     print(io, "ODE alg.: $(sol.alg)\n")
     print(io, "abstol = $(sol.abstol)\n")
@@ -166,33 +183,46 @@ A structure storing the results and some information from solving trajectories o
 - `ntraj::Int`: Number of trajectories
 - `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times_states`.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
-- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
-- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`
+- `states::AbstractVecOrMat`: The list of result states in each trajectory and each time point in `times_states`.
+- `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
 - `converged::Bool`: Whether the solution is converged or not.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
 - `reltol::Real`: The relative tolerance which is used during the solving process.
+
+# Notes
+
+When the keyword argument `keep_runs_results` is passed as `Val(false)` to a multi-trajectory solver, the `states` and `expect` fields store only the average results (averaged over all trajectories). The results can be accessed by the following index-order:
+
+- `sol.states[time_idx]`
+- `sol.expect[e_op,time_idx]`
+
+If the keyword argument `keep_runs_results = Val(true)`, the results for each trajectory and each time are stored, and the index-order of the elements in fields `states` and `expect` are:
+
+- `sol.states[trajectory,time_idx]`
+- `sol.expect[e_op,trajectory,time_idx]`
+
+We also provide the following functions for statistical analysis of multi-trajectory solutions:
+
+- [`average_states`](@ref)
+- [`average_expect`](@ref)
+- [`std_expect`](@ref)
 """
 struct TimeEvolutionStochasticSol{
     TT1<:AbstractVector{<:Real},
     TT2<:AbstractVector{<:Real},
-    TS<:AbstractVector,
-    TE<:Union{AbstractMatrix,Nothing},
-    TEA<:Union{AbstractArray,Nothing},
+    TS<:AbstractVecOrMat,
+    TE<:Union{AbstractArray,Nothing},
     TEM<:Union{AbstractArray,Nothing},
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-}
+} <: TimeEvolutionMultiTrajSol{TS,TE}
     ntraj::Int
     times::TT1
     times_states::TT2
     states::TS
     expect::TE
-    average_expect::TE # Currently just a synonym for `expect`
-    runs_expect::TEA
     measurement::TEM
     converged::Bool
     alg::AlgT
@@ -205,11 +235,11 @@ function Base.show(io::IO, sol::TimeEvolutionStochasticSol)
     print(io, "(converged: $(sol.converged))\n")
     print(io, "--------------------------------\n")
     print(io, "num_trajectories = $(sol.ntraj)\n")
-    print(io, "num_states = $(length(sol.states[1]))\n")
+    print(io, "num_states = $(size(sol.states, ndims(sol.states)))\n") # get the size of last dimension
     if sol.expect isa Nothing
         print(io, "num_expect = 0\n")
     else
-        print(io, "num_expect = $(size(sol.average_expect, 1))\n")
+        print(io, "num_expect = $(size(sol.expect, 1))\n")
     end
     print(io, "SDE alg.: $(sol.alg)\n")
     print(io, "abstol = $(sol.abstol)\n")
@@ -217,6 +247,62 @@ function Base.show(io::IO, sol::TimeEvolutionStochasticSol)
     return nothing
 end
 
+@doc raw"""
+    average_states(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-averaged result states (as density [`Operator`](@ref)) at each time point.
+"""
+average_states(sol::TimeEvolutionMultiTrajSol{<:Matrix{<:QuantumObject}}) = _average_traj_states(sol.states)
+average_states(sol::TimeEvolutionMultiTrajSol{<:Vector{<:QuantumObject}}) = sol.states  # this case should already be averaged over all trajectories
+
+# TODO: Check if broadcasting division ./ size(states, 1) is type stable
+_average_traj_states(states::Matrix{<:QuantumObject{Ket}}) =
+    map(x -> x / size(states, 1), dropdims(sum(ket2dm, states, dims = 1), dims = 1))
+_average_traj_states(states::Matrix{<:QuantumObject{ObjType}}) where {ObjType<:Union{Operator,OperatorKet}} =
+    map(x -> x / size(states, 1), dropdims(sum(states, dims = 1), dims = 1))
+
+@doc raw"""
+    average_expect(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-averaged expectation values at each time point.
+"""
+average_expect(sol::TimeEvolutionMultiTrajSol{TS,Array{T,3}}) where {TS,T<:Number} = _average_traj_expect(sol.expect)
+average_expect(sol::TimeEvolutionMultiTrajSol{TS,Matrix{T}}) where {TS,T<:Number} = sol.expect  # this case should already be averaged over all trajectories
+average_expect(::TimeEvolutionMultiTrajSol{TS,Nothing}) where {TS} = nothing
+
+_average_traj_expect(expvals::Array{T,3}) where {T<:Number} =
+    dropdims(sum(expvals, dims = 2), dims = 2) ./ size(expvals, 2)
+
+# these are used in multi-trajectory solvers before returning solutions
+_store_multitraj_states(states::Matrix{<:QuantumObject}, keep_runs_results::Val{false}) = _average_traj_states(states)
+_store_multitraj_states(states::Matrix{<:QuantumObject}, keep_runs_results::Val{true}) = states
+_store_multitraj_expect(expvals::Array{T,3}, keep_runs_results::Val{false}) where {T<:Number} =
+    _average_traj_expect(expvals)
+_store_multitraj_expect(expvals::Array{T,3}, keep_runs_results::Val{true}) where {T<:Number} = expvals
+_store_multitraj_expect(expvals::Nothing, keep_runs_results) = nothing
+
+@doc raw"""
+    std_expect(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-wise standard deviation of the expectation values at each time point.
+"""
+function std_expect(sol::TimeEvolutionMultiTrajSol{TS,Array{T,3}}) where {TS,T<:Number}
+    # the following standard deviation (std) is defined as the square-root of variance instead of pseudo-variance
+    # i.e., it is equivalent to (even for complex expectation values):
+    #    dropdims(
+    #        sqrt.(mean(abs2.(sol.expect), dims = 2) .- abs2.(mean(sol.expect, dims = 2))),
+    #        dims = 2
+    #    )
+    # [this should be included in the runtest]
+    return dropdims(std(sol.expect, corrected = false, dims = 2), dims = 2)
+end
+std_expect(::TimeEvolutionMultiTrajSol{TS,Matrix{T}}) where {TS,T<:Number} = throw(
+    ArgumentError(
+        "Can not compute the standard deviation without the expectation values of each trajectory. Try to specify keyword argument `keep_runs_results=Val(true)` to the solver.",
+    ),
+)
+std_expect(::TimeEvolutionMultiTrajSol{TS,Nothing}) where {TS} = nothing
+
 #######################################
 #=
     Callbacks for Monte Carlo quantum trajectories
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 98960764e..305822611 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -249,8 +249,8 @@ function dfd_mesolve(
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
         sol.alg,
-        NamedTuple(sol.prob.kwargs).abstol,
-        NamedTuple(sol.prob.kwargs).reltol,
+        sol.prob.kwargs[:abstol],
+        sol.prob.kwargs[:reltol],
     )
 end
 
diff --git a/test/Project.toml b/test/Project.toml
index 6c91482ab..b783960dd 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,6 +5,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 99bbb198c..d10a7c6d6 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -248,6 +248,7 @@ end
 
 @testitem "mcsolve" setup=[TESetup] begin
     using SciMLOperators
+    using Statistics
 
     # Get parameters from TESetup to simplify the code
     H = TESetup.H
@@ -270,7 +271,9 @@ end
         progress_bar = Val(false),
         jump_callback = DiscreteLindbladJumpCallback(),
     )
-    sol_mc_states = mcsolve(H, ψ0, tlist, c_ops, saveat = saveat, progress_bar = Val(false))
+    sol_mc3 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), keep_runs_results = Val(true))
+    sol_mc_states =
+        mcsolve(H, ψ0, tlist, c_ops, saveat = saveat, progress_bar = Val(false), keep_runs_results = Val(true))
     sol_mc_states2 = mcsolve(
         H,
         ψ0,
@@ -279,26 +282,31 @@ end
         saveat = saveat,
         progress_bar = Val(false),
         jump_callback = DiscreteLindbladJumpCallback(),
+        keep_runs_results = Val(true),
     )
 
-    ρt_mc = [ket2dm.(normalize.(states)) for states in sol_mc_states.states]
-    expect_mc_states = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc)
-    expect_mc_states_mean = sum(expect_mc_states, dims = 2) / size(expect_mc_states, 2)
-
-    ρt_mc2 = [ket2dm.(normalize.(states)) for states in sol_mc_states2.states]
-    expect_mc_states2 = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc2)
-    expect_mc_states_mean2 = sum(expect_mc_states2, dims = 2) / size(expect_mc_states2, 2)
+    # also test function average_states
+    # average the states from all trajectories, and then calculate the expectation value
+    expect_mc_states_mean = expect.(Ref(e_ops[1]), average_states(sol_mc_states))
+    expect_mc_states_mean2 = expect.(Ref(e_ops[1]), average_states(sol_mc_states2))
 
     @test prob_mc.prob.f.f isa MatrixOperator
     @test sum(abs, sol_mc.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, sol_mc2.expect .- sol_me.expect) / length(tlist) < 0.1
-    @test sum(abs, vec(expect_mc_states_mean) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
-    @test sum(abs, vec(expect_mc_states_mean2) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
+    @test sum(abs, average_expect(sol_mc3) .- sol_me.expect) / length(tlist) < 0.1
+    @test sum(abs, expect_mc_states_mean .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
+    @test sum(abs, expect_mc_states_mean2 .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
     @test length(sol_mc.times) == length(tlist)
     @test length(sol_mc.times_states) == 1
     @test size(sol_mc.expect) == (length(e_ops), length(tlist))
+    @test size(sol_mc.states) == (1,)
+    @test length(sol_mc3.times) == length(tlist)
+    @test length(sol_mc3.times_states) == 1
+    @test size(sol_mc3.expect) == (length(e_ops), 500, length(tlist)) # ntraj = 500
+    @test size(sol_mc3.states) == (500, 1) # ntraj = 500
     @test length(sol_mc_states.times) == length(tlist)
     @test length(sol_mc_states.times_states) == length(saveat)
+    @test size(sol_mc_states.states) == (500, length(saveat)) # ntraj = 500
     @test sol_mc_states.expect === nothing
 
     sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
@@ -308,7 +316,7 @@ end
           "(converged: $(sol_mc.converged))\n" *
           "--------------------------------\n" *
           "num_trajectories = $(sol_mc.ntraj)\n" *
-          "num_states = $(length(sol_mc.states[1]))\n" *
+          "num_states = $(size(sol_mc.states, ndims(sol_mc.states)))\n" *
           "num_expect = $(size(sol_mc.expect, 1))\n" *
           "ODE alg.: $(sol_mc.alg)\n" *
           "abstol = $(sol_mc.abstol)\n" *
@@ -318,7 +326,7 @@ end
           "(converged: $(sol_mc_states.converged))\n" *
           "--------------------------------\n" *
           "num_trajectories = $(sol_mc_states.ntraj)\n" *
-          "num_states = $(length(sol_mc_states.states[1]))\n" *
+          "num_states = $(size(sol_mc_states.states, ndims(sol_mc_states.states)))\n" *
           "num_expect = 0\n" *
           "ODE alg.: $(sol_mc_states.alg)\n" *
           "abstol = $(sol_mc_states.abstol)\n" *
@@ -330,17 +338,73 @@ end
     @test_throws ArgumentError mcsolve(H, ψ0, tlist, c_ops, save_idxs = [1, 2], progress_bar = Val(false))
     @test_throws DimensionMismatch mcsolve(H, TESetup.ψ_wrong, tlist, c_ops, progress_bar = Val(false))
 
+    # test average_states, average_expect, and std_expect
+    expvals_all = sol_mc3.expect[:, :, 2:end] # ignore testing initial time point since its standard deviation is a very small value (basically zero)
+    stdvals = std_expect(sol_mc3)
+    @test average_states(sol_mc) == sol_mc.states
+    @test average_expect(sol_mc) == sol_mc.expect
+    @test size(stdvals) == (length(e_ops), length(tlist))
+    @test all(
+        isapprox.(
+            stdvals[:, 2:end], # ignore testing initial time point since its standard deviation is a very small value (basically zero)
+            dropdims(sqrt.(mean(abs2.(expvals_all), dims = 2) .- abs2.(mean(expvals_all, dims = 2))), dims = 2);
+            atol = 1e-6,
+        ),
+    )
+    @test average_expect(sol_mc_states) === nothing
+    @test std_expect(sol_mc_states) === nothing
+    @test_throws ArgumentError std_expect(sol_mc)
+
     @testset "Memory Allocations (mcsolve)" begin
         ntraj = 100
-        allocs_tot = @allocations mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, ntraj = ntraj, progress_bar = Val(false)) # Warm-up
-        allocs_tot = @allocations mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, ntraj = ntraj, progress_bar = Val(false))
-        @test allocs_tot < 120 * ntraj + 400 # 150 allocations per trajectory + 500 for initialization
-
-        allocs_tot =
-            @allocations mcsolve(H, ψ0, tlist, c_ops, ntraj = ntraj, saveat = [tlist[end]], progress_bar = Val(false)) # Warm-up
-        allocs_tot =
-            @allocations mcsolve(H, ψ0, tlist, c_ops, ntraj = ntraj, saveat = [tlist[end]], progress_bar = Val(false))
-        @test allocs_tot < 110 * ntraj + 300 # 100 allocations per trajectory + 300 for initialization
+        for keep_runs_results in (Val(false), Val(true))
+            n1 = QuantumToolbox.getVal(keep_runs_results) ? 120 : 140
+            n2 = QuantumToolbox.getVal(keep_runs_results) ? 110 : 130
+
+            allocs_tot = @allocations mcsolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = Val(true),
+            ) # Warm-up
+            allocs_tot = @allocations mcsolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = Val(true),
+            )
+            @test allocs_tot < n1 * ntraj + 400 # 150 allocations per trajectory + 500 for initialization
+
+            allocs_tot = @allocations mcsolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = Val(true),
+            ) # Warm-up
+            allocs_tot = @allocations mcsolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = Val(true),
+            )
+            @test allocs_tot < n2 * ntraj + 300 # 100 allocations per trajectory + 300 for initialization
+        end
     end
 
     @testset "Type Inference (mcsolve)" begin
@@ -400,6 +464,7 @@ end
     @test sum(abs, sol_sse.expect .- sol_me.expect) / length(tlist) < 0.1
     @test length(sol_sse.times) == length(tlist)
     @test length(sol_sse.times_states) == 1
+    @test size(sol_sse.states) == (1,) # ntraj = 500 but keep_runs_results = Val(false)
     @test size(sol_sse.expect) == (length(e_ops), length(tlist))
     @test isnothing(sol_sse.measurement)
     @test size(sol_sse2.measurement) == (length(c_ops), 20, length(tlist) - 1)
@@ -410,7 +475,7 @@ end
           "(converged: $(sol_sse.converged))\n" *
           "--------------------------------\n" *
           "num_trajectories = $(sol_sse.ntraj)\n" *
-          "num_states = $(length(sol_sse.states[1]))\n" *
+          "num_states = $(size(sol_sse.states, ndims(sol_sse.states)))\n" *
           "num_expect = $(size(sol_sse.expect, 1))\n" *
           "SDE alg.: $(sol_sse.alg)\n" *
           "abstol = $(sol_sse.abstol)\n" *
@@ -422,15 +487,54 @@ end
 
     @testset "Memory Allocations (ssesolve)" begin
         ntraj = 100
-        allocs_tot = @allocations ssesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, ntraj = ntraj, progress_bar = Val(false)) # Warm-up
-        allocs_tot = @allocations ssesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, ntraj = ntraj, progress_bar = Val(false))
-        @test allocs_tot < 1100 * ntraj + 400 # TODO: Fix this high number of allocations
-
-        allocs_tot =
-            @allocations ssesolve(H, ψ0, tlist, c_ops, ntraj = ntraj, saveat = [tlist[end]], progress_bar = Val(false)) # Warm-up
-        allocs_tot =
-            @allocations ssesolve(H, ψ0, tlist, c_ops, ntraj = ntraj, saveat = [tlist[end]], progress_bar = Val(false))
-        @test allocs_tot < 1000 * ntraj + 300 # TODO: Fix this high number of allocations
+        for keep_runs_results in (Val(false), Val(true))
+            n1 = QuantumToolbox.getVal(keep_runs_results) ? 1100 : 1120
+            n2 = QuantumToolbox.getVal(keep_runs_results) ? 1000 : 1020
+
+            allocs_tot = @allocations ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n1 * ntraj + 400 # TODO: Fix this high number of allocations
+
+            allocs_tot = @allocations ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations ssesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n2 * ntraj + 300 # TODO: Fix this high number of allocations
+        end
     end
 
     @testset "Type Inference (ssesolve)" begin
@@ -534,10 +638,11 @@ end
     @test sum(abs, sol_sme3.expect .- sol_me.expect) / length(tlist) < 0.1
     @test length(sol_sme.times) == length(tlist)
     @test length(sol_sme.times_states) == 1
+    @test size(sol_sme.states) == (1,) # ntraj = 500 but keep_runs_results = Val(false)
     @test size(sol_sme.expect) == (length(e_ops), length(tlist))
     @test isnothing(sol_sme.measurement)
     @test size(sol_sme2.measurement) == (length(sc_ops_sme), 20, length(tlist) - 1)
-    @test all([sol_sme4.states[j][i] ≈ vector_to_operator(sol_sme5.states[j][i]) for i in eachindex(saveat), j in 1:10])
+    @test all([sol_sme4.states[i] ≈ vector_to_operator(sol_sme5.states[i]) for i in eachindex(saveat)])
 
     sol_sme_string = sprint((t, s) -> show(t, "text/plain", s), sol_sme)
     @test sol_sme_string ==
@@ -545,7 +650,7 @@ end
           "(converged: $(sol_sme.converged))\n" *
           "--------------------------------\n" *
           "num_trajectories = $(sol_sme.ntraj)\n" *
-          "num_states = $(length(sol_sme.states[1]))\n" *
+          "num_states = $(size(sol_sme.states, ndims(sol_sme.states)))\n" *
           "num_expect = $(size(sol_sme.expect, 1))\n" *
           "SDE alg.: $(sol_sme.alg)\n" *
           "abstol = $(sol_sme.abstol)\n" *
@@ -557,94 +662,109 @@ end
 
     @testset "Memory Allocations (smesolve)" begin
         ntraj = 100
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme,
-            sc_ops_sme,
-            e_ops = e_ops,
-            ntraj = ntraj,
-            progress_bar = Val(false),
-        ) # Warm-up
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme,
-            sc_ops_sme,
-            e_ops = e_ops,
-            ntraj = ntraj,
-            progress_bar = Val(false),
-        )
-        @test allocs_tot < 1100 * ntraj + 2300 # TODO: Fix this high number of allocations
-
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme,
-            sc_ops_sme,
-            ntraj = ntraj,
-            saveat = [tlist[end]],
-            progress_bar = Val(false),
-        ) # Warm-up
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme,
-            sc_ops_sme,
-            ntraj = ntraj,
-            saveat = [tlist[end]],
-            progress_bar = Val(false),
-        )
-        @test allocs_tot < 1000 * ntraj + 1500 # TODO: Fix this high number of allocations
-
-        # Diagonal Noise Case
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme2,
-            sc_ops_sme2,
-            e_ops = e_ops,
-            ntraj = ntraj,
-            progress_bar = Val(false),
-        ) # Warm-up
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme2,
-            sc_ops_sme2,
-            e_ops = e_ops,
-            ntraj = 1,
-            progress_bar = Val(false),
-        )
-        @test allocs_tot < 600 * ntraj + 1400 # TODO: Fix this high number of allocations
-
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme2,
-            sc_ops_sme2,
-            ntraj = ntraj,
-            saveat = [tlist[end]],
-            progress_bar = Val(false),
-        ) # Warm-up
-        allocs_tot = @allocations smesolve(
-            H,
-            ψ0,
-            tlist,
-            c_ops_sme2,
-            sc_ops_sme2,
-            ntraj = 1,
-            saveat = [tlist[end]],
-            progress_bar = Val(false),
-        )
-        @test allocs_tot < 550 * ntraj + 1000 # TODO: Fix this high number of allocations
+        for keep_runs_results in (Val(false), Val(true))
+            n1 = QuantumToolbox.getVal(keep_runs_results) ? 1100 : 1120
+            n2 = QuantumToolbox.getVal(keep_runs_results) ? 1000 : 1020
+            n3 = QuantumToolbox.getVal(keep_runs_results) ? 600 : 620
+            n4 = QuantumToolbox.getVal(keep_runs_results) ? 550 : 570
+
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme,
+                sc_ops_sme,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme,
+                sc_ops_sme,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n1 * ntraj + 2300 # TODO: Fix this high number of allocations
+
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme,
+                sc_ops_sme,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme,
+                sc_ops_sme,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n2 * ntraj + 1500 # TODO: Fix this high number of allocations
+
+            # Diagonal Noise Case
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme2,
+                sc_ops_sme2,
+                e_ops = e_ops,
+                ntraj = ntraj,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme2,
+                sc_ops_sme2,
+                e_ops = e_ops,
+                ntraj = 1,
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n3 * ntraj + 1400 # TODO: Fix this high number of allocations
+
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme2,
+                sc_ops_sme2,
+                ntraj = ntraj,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            ) # Warm-up
+            allocs_tot = @allocations smesolve(
+                H,
+                ψ0,
+                tlist,
+                c_ops_sme2,
+                sc_ops_sme2,
+                ntraj = 1,
+                saveat = [tlist[end]],
+                progress_bar = Val(false),
+                keep_runs_results = keep_runs_results,
+            )
+            @test allocs_tot < n4 * ntraj + 1000 # TODO: Fix this high number of allocations
+        end
     end
 
     @testset "Type Inference (smesolve)" begin
@@ -800,45 +920,114 @@ end
     rng = TESetup.rng
 
     rng = MersenneTwister(1234)
-    sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_mc1 =
+        mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng, keep_runs_results = Val(true))
     rng = MersenneTwister(1234)
-    sol_sse1 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sse1 = ssesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops,
+        ntraj = 50,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
     rng = MersenneTwister(1234)
-    sol_sme1 =
-        smesolve(H, ψ0, tlist, c_ops_sme, sc_ops_sme, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sme1 = smesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops_sme,
+        sc_ops_sme,
+        ntraj = 50,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
 
     rng = MersenneTwister(1234)
-    sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_mc2 =
+        mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng, keep_runs_results = Val(true))
     rng = MersenneTwister(1234)
-    sol_sse2 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sse2 = ssesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops,
+        ntraj = 50,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
     rng = MersenneTwister(1234)
-    sol_sme2 =
-        smesolve(H, ψ0, tlist, c_ops_sme, sc_ops_sme, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sme2 = smesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops_sme,
+        sc_ops_sme,
+        ntraj = 50,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
 
     rng = MersenneTwister(1234)
-    sol_mc3 = mcsolve(H, ψ0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_mc3 = mcsolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops,
+        ntraj = 510,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
     rng = MersenneTwister(1234)
-    sol_sse3 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sse3 = ssesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops,
+        ntraj = 60,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
     rng = MersenneTwister(1234)
-    sol_sme3 =
-        smesolve(H, ψ0, tlist, c_ops_sme, sc_ops_sme, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+    sol_sme3 = smesolve(
+        H,
+        ψ0,
+        tlist,
+        c_ops_sme,
+        sc_ops_sme,
+        ntraj = 60,
+        e_ops = e_ops,
+        progress_bar = Val(false),
+        rng = rng,
+        keep_runs_results = Val(true),
+    )
 
     @test sol_mc1.expect ≈ sol_mc2.expect atol = 1e-10
-    @test sol_mc1.runs_expect ≈ sol_mc2.runs_expect atol = 1e-10
     @test sol_mc1.col_times ≈ sol_mc2.col_times atol = 1e-10
     @test sol_mc1.col_which ≈ sol_mc2.col_which atol = 1e-10
 
-    @test sol_mc1.runs_expect ≈ sol_mc3.runs_expect[:, 1:500, :] atol = 1e-10
+    @test sol_mc1.expect ≈ sol_mc3.expect[:, 1:500, :] atol = 1e-10
 
     @test sol_sse1.expect ≈ sol_sse2.expect atol = 1e-10
-    @test sol_sse1.runs_expect ≈ sol_sse2.runs_expect atol = 1e-10
 
-    @test sol_sse1.runs_expect ≈ sol_sse3.runs_expect[:, 1:50, :] atol = 1e-10
+    @test sol_sse1.expect ≈ sol_sse3.expect[:, 1:50, :] atol = 1e-10
 
     @test sol_sme1.expect ≈ sol_sme2.expect atol = 1e-10
-    @test sol_sme1.runs_expect ≈ sol_sme2.runs_expect atol = 1e-10
 
-    @test sol_sme1.runs_expect ≈ sol_sme3.runs_expect[:, 1:50, :] atol = 1e-10
+    @test sol_sme1.expect ≈ sol_sme3.expect[:, 1:50, :] atol = 1e-10
 end
 
 @testitem "example" begin