From b71a0e3e0707ac08252a373c5079e8ea119b0e68 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Thu, 22 May 2025 21:53:42 +0800
Subject: [PATCH 01/17] align statistical analysis methods for multi-trajectory
 with qutip

---
 docs/src/resources/api.md            |  3 ++
 src/time_evolution/mcsolve.jl        |  4 --
 src/time_evolution/smesolve.jl       |  6 ---
 src/time_evolution/ssesolve.jl       |  6 ---
 src/time_evolution/time_evolution.jl | 79 ++++++++++++++++++++--------
 5 files changed, 59 insertions(+), 39 deletions(-)

diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index 0358ceb1d..ec3124986 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -192,6 +192,9 @@ TimeEvolutionProblem
 TimeEvolutionSol
 TimeEvolutionMCSol
 TimeEvolutionStochasticSol
+average_states
+average_expect
+std_expect
 sesolveProblem
 mesolveProblem
 mcsolveProblem
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 7db7048ff..8ab40cced 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -407,14 +407,10 @@ function mcsolve(
     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,
         states,
-        expvals,
-        expvals, # This is average_expect
         expvals_all,
         col_times,
         col_which,
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 710029197..a9e5f20be 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -418,16 +418,10 @@ function smesolve(
         _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)
-
     return TimeEvolutionStochasticSol(
         ntraj,
         ens_prob.times,
         states,
-        expvals,
-        expvals, # This is average_expect
         expvals_all,
         m_expvals, # Measurement expectation values
         sol.converged,
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 6945dc81e..fd4aa976e 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -412,16 +412,10 @@ function ssesolve(
         _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,
         states,
-        expvals,
-        expvals, # This is average_expect
         expvals_all,
         m_expvals, # Measurement expectation values
         sol.converged,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 30730bc9d..5bbf5aaee 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{Texpect} end
+
 @doc raw"""
     struct TimeEvolutionProblem
 
@@ -89,7 +93,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.
 
@@ -97,34 +101,37 @@ A structure storing the results and some information from solving quantum trajec
 
 - `ntraj::Int`: Number of trajectories
 - `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory.
-- `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::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
+- `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.
+
+# Methods
+
+We also provide the following functions for analyzing statistics from multi-trajectory solutions.
+
+- [`average_states`](@ref)
+- [`average_expect`](@ref)
+- [`std_expect`](@ref)
 """
 struct TimeEvolutionMCSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
-    TE<:Union{AbstractMatrix,Nothing},
-    TEA<:Union{AbstractArray,Nothing},
+    TE<:Union{AbstractArray,Nothing},
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
     AlgT<:OrdinaryDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-}
+} <: TimeEvolutionMultiTrajSol{TE}
     ntraj::Int
     times::TT
     states::TS
     expect::TE
-    average_expect::TE # Currently just a synonym for `expect`
-    runs_expect::TEA
     col_times::TJT
     col_which::TJW
     converged::Bool
@@ -142,7 +149,7 @@ function Base.show(io::IO, sol::TimeEvolutionMCSol)
     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")
@@ -151,7 +158,7 @@ function Base.show(io::IO, sol::TimeEvolutionMCSol)
 end
 
 @doc raw"""
-    struct TimeEvolutionStochasticSol
+    struct TimeEvolutionStochasticSol <: TimeEvolutionMultiTrajSol
 
 A structure storing the results and some information from solving trajectories of the Stochastic time evolution.
 
@@ -159,31 +166,34 @@ A structure storing the results and some information from solving trajectories o
 
 - `ntraj::Int`: Number of trajectories
 - `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory.
-- `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::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
+- `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.
+
+# Methods
+
+We also provide the following functions for analyzing statistics from multi-trajectory solutions.
+
+- [`average_states`](@ref)
+- [`average_expect`](@ref)
+- [`std_expect`](@ref)
 """
 struct TimeEvolutionStochasticSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
-    TE<:Union{AbstractMatrix,Nothing},
-    TEA<:Union{AbstractArray,Nothing},
+    TE<:Union{AbstractArray,Nothing},
     TEM<:Union{AbstractArray,Nothing},
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-}
+} <: TimeEvolutionMultiTrajSol{TE}
     ntraj::Int
     times::TT
     states::TS
     expect::TE
-    average_expect::TE # Currently just a synonym for `expect`
-    runs_expect::TEA
     measurement::TEM
     converged::Bool
     alg::AlgT
@@ -200,7 +210,7 @@ function Base.show(io::IO, sol::TimeEvolutionStochasticSol)
     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")
@@ -208,6 +218,29 @@ function Base.show(io::IO, sol::TimeEvolutionStochasticSol)
     return nothing
 end
 
+@doc raw"""
+    average_states(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-averaged result states at each time point.
+"""
+average_states(sol::TimeEvolutionMultiTrajSol) = mean(sol.states)
+
+@doc raw"""
+    average_expect(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-averaged expectation values at each time point.
+"""
+average_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} = dropdims(mean(sol.runs_expect, dims = 2), dims = 2)
+average_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
+
+@doc raw"""
+    std_expect(sol::TimeEvolutionMultiTrajSol)
+
+Return the trajectory-wise standard deviation of the expectation values at each time point.
+"""
+std_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} = dropdims(std(sol.runs_expect, dims = 2), dims = 2)
+std_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
+
 #######################################
 #=
     Callbacks for Monte Carlo quantum trajectories

From 767725a3480a9de23d1993e975bb769b2e3a7bd3 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 00:02:30 +0800
Subject: [PATCH 02/17] add `Statistics.jl` to dependency

---
 Project.toml          | 2 ++
 src/QuantumToolbox.jl | 1 +
 2 files changed, 3 insertions(+)

diff --git a/Project.toml b/Project.toml
index cdff9028e..c6184bbd0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,6 +23,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]
@@ -63,6 +64,7 @@ SciMLOperators = "0.3, 0.4"
 SparseArrays = "1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
+Statistics = "1"
 StochasticDiffEq = "6"
 Test = "1"
 julia = "1.10"
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 40bfc245f..71f579429 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:

From fed47906d46c2ddf8cd29c4097bf5be952e7052f Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 04:06:27 +0800
Subject: [PATCH 03/17] minor changes

---
 src/time_evolution/time_evolution.jl | 35 +++++++++++++++++++++-------
 1 file changed, 27 insertions(+), 8 deletions(-)

diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 5bbf5aaee..fd9f86d60 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -110,9 +110,14 @@ A structure storing the results and some information from solving quantum trajec
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 
-# Methods
+# Notes
 
-We also provide the following functions for analyzing statistics from multi-trajectory solutions.
+The index-order of the elements in fields `states` and `expect` are:
+
+- `sol.states[trajectory][time]`
+- `sol.expect[e_op,trajectory,time]`
+
+We also provide the following functions for statistical analysis of multi-trajectory solutions.
 
 - [`average_states`](@ref)
 - [`average_expect`](@ref)
@@ -173,9 +178,14 @@ A structure storing the results and some information from solving trajectories o
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 
-# Methods
+# Notes
+
+The index-order of the elements in fields `states` and `expect` are:
 
-We also provide the following functions for analyzing statistics from multi-trajectory solutions.
+- `sol.states[trajectory][time]`
+- `sol.expect[e_op,trajectory,time]`
+
+We also provide the following functions for statistical analysis of multi-trajectory solutions.
 
 - [`average_states`](@ref)
 - [`average_expect`](@ref)
@@ -221,16 +231,17 @@ end
 @doc raw"""
     average_states(sol::TimeEvolutionMultiTrajSol)
 
-Return the trajectory-averaged result states at each time point.
+Return the trajectory-averaged result states (as density [`Operator`](@ref)) at each time point.
 """
-average_states(sol::TimeEvolutionMultiTrajSol) = mean(sol.states)
+average_states(sol::TimeEvolutionMultiTrajSol) = mean(ket2dm, sol.states)
 
 @doc raw"""
     average_expect(sol::TimeEvolutionMultiTrajSol)
 
 Return the trajectory-averaged expectation values at each time point.
 """
-average_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} = dropdims(mean(sol.runs_expect, dims = 2), dims = 2)
+average_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} =
+    dropdims(mean(sol.expect, dims = 2), dims = 2)
 average_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
 
 @doc raw"""
@@ -238,7 +249,15 @@ average_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
 
 Return the trajectory-wise standard deviation of the expectation values at each time point.
 """
-std_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} = dropdims(std(sol.runs_expect, dims = 2), dims = 2)
+function std_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray}
+    # 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
+    #    )
+    return dropdims(std(sol.expect, corrected = false, dims = 2), dims = 2)
+end
 std_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
 
 #######################################

From de132d3333c8e4997135e98a905aea2aaee2d4a4 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 04:06:33 +0800
Subject: [PATCH 04/17] update tests

---
 test/core-test/dynamical-shifted-fock.jl |  6 +--
 test/core-test/time_evolution.jl         | 56 +++++++++++++++---------
 test/runtests.jl                         | 33 +++++++-------
 3 files changed, 56 insertions(+), 39 deletions(-)

diff --git a/test/core-test/dynamical-shifted-fock.jl b/test/core-test/dynamical-shifted-fock.jl
index 663ac6494..690b708fe 100644
--- a/test/core-test/dynamical-shifted-fock.jl
+++ b/test/core-test/dynamical-shifted-fock.jl
@@ -63,7 +63,7 @@
     )
     val_ss = abs2(sol0.expect[1, end])
     @test sum(abs2.(sol0.expect[1, :] .- sol_dsf_me.expect[1, :])) / (val_ss * length(tlist)) < 0.1
-    @test sum(abs2.(sol0.expect[1, :] .- sol_dsf_mc.expect[1, :])) / (val_ss * length(tlist)) < 0.1
+    @test sum(abs2.(sol0.expect[1, :] .- average_expect(sol_dsf_mc)[1, :])) / (val_ss * length(tlist)) < 0.1
 
     # Two cavities case
     F = 2
@@ -143,7 +143,7 @@
 
     val_ss = abs2(sol0.expect[1, end])
     @test sum(abs2.(sol0.expect[1, :] .- sol_dsf_me.expect[1, :])) / (val_ss * length(tlist)) < 0.6
-    @test sum(abs2.(sol0.expect[1, :] .- sol_dsf_mc.expect[1, :])) / (val_ss * length(tlist)) < 0.6
+    @test sum(abs2.(sol0.expect[1, :] .- average_expect(sol_dsf_mc)[1, :])) / (val_ss * length(tlist)) < 0.6
     @test sum(abs2.(sol0.expect[2, :] .- sol_dsf_me.expect[2, :])) / (val_ss * length(tlist)) < 0.6
-    @test sum(abs2.(sol0.expect[2, :] .- sol_dsf_mc.expect[2, :])) / (val_ss * length(tlist)) < 0.6
+    @test sum(abs2.(sol0.expect[2, :] .- average_expect(sol_dsf_mc)[2, :])) / (val_ss * length(tlist)) < 0.6
 end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 5abb7d18d..1c273c892 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -186,6 +186,11 @@
         # Redirect to `sesolve`
         sol_me5 = mesolve(H, ψ0, tlist, progress_bar = Val(false))
 
+        # standard deviation
+        std_mc = dropdims(sqrt.(mean(abs2.(sol_mc.expect), dims = 2) .- abs2.(mean(sol_mc.expect, dims = 2))), dims = 2)
+        std_mc2 =
+            dropdims(sqrt.(mean(abs2.(sol_mc2.expect), dims = 2) .- abs2.(mean(sol_mc2.expect, dims = 2))), dims = 2)
+
         ρ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)
@@ -202,13 +207,17 @@
         @test prob_me.prob.f.f isa MatrixOperator
         @test prob_mc.prob.f.f isa MatrixOperator
         @test isket(sol_me5.states[1])
-        @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 all(isapprox.(average_states(sol_mc), sol_me.states; atol = 0.1))
+        @test all(isapprox.(average_states(sol_mc2), sol_me.states; atol = 0.1))
+        @test all(isapprox.(average_expect(sol_mc), sol_me.expect; atol = 0.1))
+        @test all(isapprox.(average_expect(sol_mc2), sol_me.expect; atol = 0.1))
+        @test all(std_expect(sol_mc) .≈ std_mc)
+        @test all(std_expect(sol_mc2) .≈ std_mc2)
         @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, sol_sse.expect .- sol_me.expect) / length(tlist) < 0.1
-        @test sum(abs, sol_sme.expect .- sol_me.expect) / length(tlist) < 0.1
-        @test sum(abs, sol_sme3.expect .- sol_me.expect) / length(tlist) < 0.1
+        @test sum(abs, average_expect(sol_sse) .- sol_me.expect) / length(tlist) < 0.1
+        @test sum(abs, average_expect(sol_sme) .- sol_me.expect) / length(tlist) < 0.1
+        @test sum(abs, average_expect(sol_sme3) .- sol_me.expect) / length(tlist) < 0.1
         @test length(sol_me.times) == length(tlist)
         @test length(sol_me.states) == 1
         @test size(sol_me.expect) == (length(e_ops), length(tlist))
@@ -221,13 +230,15 @@
         @test sol_me3.expect[1, saveat_idxs] ≈ expect(e_ops[1], sol_me3.states) atol = 1e-6
         @test all([sol_me3.states[i] ≈ vector_to_operator(sol_me4.states[i]) for i in eachindex(saveat)])
         @test length(sol_mc.times) == length(tlist)
-        @test size(sol_mc.expect) == (length(e_ops), length(tlist))
+        @test size(sol_mc.expect) == (length(e_ops), 500, length(tlist))
         @test length(sol_mc_states.times) == length(tlist)
         @test sol_mc_states.expect === nothing
+        @test average_expect(sol_mc_states) === nothing
+        @test std_expect(sol_mc_states) === nothing
         @test length(sol_sse.times) == length(tlist)
-        @test size(sol_sse.expect) == (length(e_ops), length(tlist))
+        @test size(sol_sse.expect) == (length(e_ops), 500, length(tlist))
         @test length(sol_sme.times) == length(tlist)
-        @test size(sol_sme.expect) == (length(e_ops), length(tlist))
+        @test size(sol_sme.expect) == (length(e_ops), 500, length(tlist))
         @test isnothing(sol_sse.measurement)
         @test isnothing(sol_sme.measurement)
         @test size(sol_sse2.measurement) == (length(c_ops), 20, length(tlist) - 1)
@@ -336,7 +347,7 @@
 
         @test sol_se.expect ≈ sol_se_td.expect atol = 1e-6 * length(tlist)
         @test sol_me.expect ≈ sol_me_td.expect atol = 1e-6 * length(tlist)
-        @test sol_mc.expect ≈ sol_mc_td.expect atol = 1e-2 * length(tlist)
+        @test average_expect(sol_mc) ≈ average_expect(sol_mc_td) atol = 1e-2 * length(tlist)
         # @test sol_sse.expect ≈ sol_sse_td.expect atol = 1e-2 * length(tlist)
 
         H_td2 = QobjEvo(H_td)
@@ -350,7 +361,7 @@
 
         @test sol_se.expect ≈ sol_se_td2.expect atol = 1e-6 * length(tlist)
         @test sol_me.expect ≈ sol_me_td2.expect atol = 1e-6 * length(tlist)
-        @test sol_mc.expect ≈ sol_mc_td2.expect atol = 1e-2 * length(tlist)
+        @test average_expect(sol_mc) ≈ average_expect(sol_mc_td2) atol = 1e-2 * length(tlist)
         # @test sol_sse.expect ≈ sol_sse_td2.expect atol = 1e-2 * length(tlist)
 
         @testset "Memory Allocations (mesolve)" begin
@@ -542,7 +553,7 @@
             @inferred mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
         end
 
-        @testset "Type Inference mcsolve" begin
+        @testset "Type Inference mcsolve and statistical analysis functions" begin
             @inferred mcsolveEnsembleProblem(
                 H,
                 ψ0,
@@ -553,10 +564,11 @@
                 progress_bar = Val(false),
                 rng = rng,
             )
-            @inferred mcsolve(H, ψ0, tlist, c_ops, ntraj = 5, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+            sol1 =
+                @inferred mcsolve(H, ψ0, tlist, c_ops, ntraj = 5, e_ops = e_ops, progress_bar = Val(false), rng = rng)
+            sol2 = @inferred mcsolve(H, ψ0_int, tlist, c_ops, ntraj = 5, progress_bar = Val(false), rng = rng)
             @inferred mcsolve(H, ψ0, tlist, c_ops, ntraj = 5, progress_bar = Val(true), rng = rng)
             @inferred mcsolve(H, ψ0, [0, 10], c_ops, ntraj = 5, progress_bar = Val(false), rng = rng)
-            @inferred mcsolve(H, ψ0_int, tlist, c_ops, ntraj = 5, progress_bar = Val(false), rng = rng)
             @inferred mcsolve(
                 H,
                 ψ0,
@@ -578,6 +590,13 @@
                 params = p,
                 rng = rng,
             )
+
+            @inferred average_states(sol1)
+            @inferred average_expect(sol1)
+            @inferred std_expect(sol1)
+            @inferred average_states(sol2)
+            @inferred average_expect(sol2)
+            @inferred std_expect(sol2)
         end
 
         @testset "Type Inference ssesolve" begin
@@ -763,21 +782,18 @@
             )
 
             @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
     end
 
@@ -818,7 +834,7 @@
         t_l = LinRange(0, 20 / γ1, 1000)
         sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = [sp1 * sm1, sp2 * sm2], progress_bar = false) # Here we don't put Val(false) because we want to test the support for Bool type
         sol_mc = mcsolve(H, psi0, t_l, c_ops, e_ops = [sp1 * sm1, sp2 * sm2], progress_bar = Val(false))
-        @test sum(abs.(sol_mc.expect[1:2, :] .- sol_me.expect[1:2, :])) / length(t_l) < 0.1
+        @test sum(abs.(average_expect(sol_mc) .- sol_me.expect)) / length(t_l) < 0.1
         @test expect(sp1 * sm1, sol_me.states[end]) ≈ expect(sigmap() * sigmam(), ptrace(sol_me.states[end], 1))
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7407d82fb..d5f2c57c6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,26 +12,27 @@ const testdir = dirname(@__FILE__)
 
 # Put core tests in alphabetical order
 core_tests = [
-    "block_diagonal_form.jl",
-    "correlations_and_spectrum.jl",
-    "dynamical_fock_dimension_mesolve.jl",
-    "dynamical-shifted-fock.jl",
-    "eigenvalues_and_operators.jl",
-    "entropy_and_metric.jl",
-    "generalized_master_equation.jl",
-    "low_rank_dynamics.jl",
-    "negativity_and_partial_transpose.jl",
-    "progress_bar.jl",
-    "quantum_objects.jl",
-    "quantum_objects_evo.jl",
-    "states_and_operators.jl",
-    "steady_state.jl",
+    # "block_diagonal_form.jl",
+    # "correlations_and_spectrum.jl",
+    # "dynamical_fock_dimension_mesolve.jl",
+    # "dynamical-shifted-fock.jl",
+    # "eigenvalues_and_operators.jl",
+    # "entropy_and_metric.jl",
+    # "generalized_master_equation.jl",
+    # "low_rank_dynamics.jl",
+    # "negativity_and_partial_transpose.jl",
+    # "progress_bar.jl",
+    # "quantum_objects.jl",
+    # "quantum_objects_evo.jl",
+    # "states_and_operators.jl",
+    # "steady_state.jl",
     "time_evolution.jl",
-    "utilities.jl",
-    "wigner.jl",
+    # "utilities.jl",
+    # "wigner.jl",
 ]
 
 if (GROUP == "All") || (GROUP == "Core")
+    import Statistics: mean
     using QuantumToolbox
     import QuantumToolbox: position, momentum
     import Random: MersenneTwister

From ec2b69de194f2ba64c158d4a5a1f04c18b2ad2a0 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 19:43:16 +0800
Subject: [PATCH 05/17] fix errors

---
 src/qobj/functions.jl                | 3 +++
 src/time_evolution/time_evolution.jl | 8 ++++----
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 820c20c30..30d13bd8a 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -16,6 +16,9 @@ ket2dm(ψ::QuantumObject{Ket}) = ψ * ψ'
 
 ket2dm(ρ::QuantumObject{Operator}) = ρ
 
+# for special usage, e.g., average_states
+ket2dm(states::Vector{<:QuantumObject{T}}) where {T<:Union{Ket,Operator}} = map(s -> ket2dm(s), states)
+
 @doc raw"""
     expect(O::Union{AbstractQuantumObject,Vector{AbstractQuantumObject}}, ψ::Union{QuantumObject,Vector{QuantumObject}})
 
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index fd9f86d60..c3116a499 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -132,7 +132,7 @@ struct TimeEvolutionMCSol{
     AlgT<:OrdinaryDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-} <: TimeEvolutionMultiTrajSol{TE}
+} <: TimeEvolutionMultiTrajSol{TS,TE}
     ntraj::Int
     times::TT
     states::TS
@@ -199,7 +199,7 @@ struct TimeEvolutionStochasticSol{
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-} <: TimeEvolutionMultiTrajSol{TE}
+} <: TimeEvolutionMultiTrajSol{TS,TE}
     ntraj::Int
     times::TT
     states::TS
@@ -242,7 +242,7 @@ Return the trajectory-averaged expectation values at each time point.
 """
 average_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray} =
     dropdims(mean(sol.expect, dims = 2), dims = 2)
-average_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
+average_expect(::TimeEvolutionMultiTrajSol{Nothing}) = nothing
 
 @doc raw"""
     std_expect(sol::TimeEvolutionMultiTrajSol)
@@ -258,7 +258,7 @@ function std_expect(sol::TimeEvolutionMultiTrajSol{TE}) where {TE<:AbstractArray
     #    )
     return dropdims(std(sol.expect, corrected = false, dims = 2), dims = 2)
 end
-std_expect(sol::TimeEvolutionMultiTrajSol{Nothing}) = nothing
+std_expect(::TimeEvolutionMultiTrajSol{Nothing}) = nothing
 
 #######################################
 #=

From ef22c1bfb65ce72a73abb5a8e48c606457656644 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 19:43:24 +0800
Subject: [PATCH 06/17] fix runtests

---
 test/core-test/time_evolution.jl | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 1c273c892..8193e067e 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -186,10 +186,14 @@
         # Redirect to `sesolve`
         sol_me5 = mesolve(H, ψ0, tlist, progress_bar = Val(false))
 
-        # standard deviation
-        std_mc = dropdims(sqrt.(mean(abs2.(sol_mc.expect), dims = 2) .- abs2.(mean(sol_mc.expect, dims = 2))), dims = 2)
-        std_mc2 =
-            dropdims(sqrt.(mean(abs2.(sol_mc2.expect), dims = 2) .- abs2.(mean(sol_mc2.expect, dims = 2))), dims = 2)
+        # average states
+        avg_mc_states = average_states(sol_mc_states)
+        avg_mc_states2 = average_states(sol_mc_states2)
+
+        # variance
+        var_mc = dropdims(abs.(mean(abs2.(sol_mc.expect), dims = 2) .- abs2.(mean(sol_mc.expect, dims = 2))), dims = 2)
+        var_mc2 =
+            dropdims(abs.(mean(abs2.(sol_mc2.expect), dims = 2) .- abs2.(mean(sol_mc2.expect, dims = 2))), dims = 2)
 
         ρ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)
@@ -207,12 +211,12 @@
         @test prob_me.prob.f.f isa MatrixOperator
         @test prob_mc.prob.f.f isa MatrixOperator
         @test isket(sol_me5.states[1])
-        @test all(isapprox.(average_states(sol_mc), sol_me.states; atol = 0.1))
-        @test all(isapprox.(average_states(sol_mc2), sol_me.states; atol = 0.1))
-        @test all(isapprox.(average_expect(sol_mc), sol_me.expect; atol = 0.1))
-        @test all(isapprox.(average_expect(sol_mc2), sol_me.expect; atol = 0.1))
-        @test all(std_expect(sol_mc) .≈ std_mc)
-        @test all(std_expect(sol_mc2) .≈ std_mc2)
+        @test all(sum.(abs, get_data.(avg_mc_states .- sol_me2.states[saveat_idxs])) .< 0.1)
+        @test all(sum.(abs, get_data.(avg_mc_states2 .- sol_me2.states[saveat_idxs])) .< 0.1)
+        @test sum(abs, average_expect(sol_mc) .- sol_me.expect) / length(tlist) < 0.1
+        @test sum(abs, average_expect(sol_mc2) .- sol_me.expect) / length(tlist) < 0.1
+        @test sum(abs, std_expect(sol_mc) .^ 2 .- var_mc) < 1e-8
+        @test sum(abs, std_expect(sol_mc2) .^ 2 .- var_mc2) < 1e-8
         @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_sse) .- sol_me.expect) / length(tlist) < 0.1
@@ -591,10 +595,10 @@
                 rng = rng,
             )
 
-            @inferred average_states(sol1)
+            # @inferred average_states(sol1) # TODO: need to change type of sol.states
             @inferred average_expect(sol1)
             @inferred std_expect(sol1)
-            @inferred average_states(sol2)
+            # @inferred average_states(sol2) # TODO: need to change type of sol.states
             @inferred average_expect(sol2)
             @inferred std_expect(sol2)
         end

From 83e937baa45a33a7899e6aa1ece8cc530dded537 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 19:47:29 +0800
Subject: [PATCH 07/17] update changelog

---
 CHANGELOG.md | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 145b17345..230e97fe1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Align statistical analysis methods for multi-trajectory solutions to `QuTiP` [#471]. Remove field `average_expect`, replace field `runs_expect` as `expect`, and introduce the following functions:
+  - `average_states`
+  - `average_expect`
+  - `std_expect`
+
 ## [v0.31.1]
 Release date: 2025-05-16
 
@@ -228,3 +233,4 @@ Release date: 2024-11-13
 [#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
 [#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
 [#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
+[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471

From 0fc8624c531d04c71550ea5e217dc1e9762124cf Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Fri, 23 May 2025 19:51:30 +0800
Subject: [PATCH 08/17] fix typo

---
 src/time_evolution/time_evolution.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index c3116a499..98042bb6b 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -132,7 +132,7 @@ struct TimeEvolutionMCSol{
     AlgT<:OrdinaryDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-} <: TimeEvolutionMultiTrajSol{TS,TE}
+} <: TimeEvolutionMultiTrajSol{TE}
     ntraj::Int
     times::TT
     states::TS
@@ -199,7 +199,7 @@ struct TimeEvolutionStochasticSol{
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
     RT<:Real,
-} <: TimeEvolutionMultiTrajSol{TS,TE}
+} <: TimeEvolutionMultiTrajSol{TE}
     ntraj::Int
     times::TT
     states::TS

From cbdc0fecbbe349f1123df3c1dea623a18160c566 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Fri, 23 May 2025 20:12:20 +0800
Subject: [PATCH 09/17] fix runtest list

---
 test/runtests.jl | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index d5f2c57c6..a4b70f83e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,23 +12,23 @@ const testdir = dirname(@__FILE__)
 
 # Put core tests in alphabetical order
 core_tests = [
-    # "block_diagonal_form.jl",
-    # "correlations_and_spectrum.jl",
-    # "dynamical_fock_dimension_mesolve.jl",
-    # "dynamical-shifted-fock.jl",
-    # "eigenvalues_and_operators.jl",
-    # "entropy_and_metric.jl",
-    # "generalized_master_equation.jl",
-    # "low_rank_dynamics.jl",
-    # "negativity_and_partial_transpose.jl",
-    # "progress_bar.jl",
-    # "quantum_objects.jl",
-    # "quantum_objects_evo.jl",
-    # "states_and_operators.jl",
-    # "steady_state.jl",
+    "block_diagonal_form.jl",
+    "correlations_and_spectrum.jl",
+    "dynamical_fock_dimension_mesolve.jl",
+    "dynamical-shifted-fock.jl",
+    "eigenvalues_and_operators.jl",
+    "entropy_and_metric.jl",
+    "generalized_master_equation.jl",
+    "low_rank_dynamics.jl",
+    "negativity_and_partial_transpose.jl",
+    "progress_bar.jl",
+    "quantum_objects.jl",
+    "quantum_objects_evo.jl",
+    "states_and_operators.jl",
+    "steady_state.jl",
     "time_evolution.jl",
-    # "utilities.jl",
-    # "wigner.jl",
+    "utilities.jl",
+    "wigner.jl",
 ]
 
 if (GROUP == "All") || (GROUP == "Core")

From 4cbc469f58d2e1e9d10a39396905913d6805ec09 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 13:05:24 +0800
Subject: [PATCH 10/17] minor changes

---
 src/time_evolution/mesolve.jl                  | 4 ++--
 src/time_evolution/sesolve.jl                  | 4 ++--
 src/time_evolution/time_evolution_dynamical.jl | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 980148889..ef5b99d56 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -215,7 +215,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 c24970f75..ea68148ce 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -165,7 +165,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/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 2e07047da..1383336c3 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -248,8 +248,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
 

From d789c8736b798622c632ce38f6c82037c312f822 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 13:27:12 +0800
Subject: [PATCH 11/17] add field `rng` for multi-trajectory solutions

---
 src/time_evolution/mcsolve.jl        |  7 ++++---
 src/time_evolution/smesolve.jl       |  3 ++-
 src/time_evolution/ssesolve.jl       |  3 ++-
 src/time_evolution/time_evolution.jl |  7 +++++++
 test/core-test/time_evolution.jl     | 21 +++++++++------------
 5 files changed, 24 insertions(+), 17 deletions(-)

diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 8ab40cced..29db4c3d0 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -256,7 +256,7 @@ function mcsolveEnsembleProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob_mc.times,
         prob_mc.dimensions,
-        (progr = _output_func[2], channel = _output_func[3]),
+        (progr = _output_func[2], channel = _output_func[3], rng = rng),
     )
 
     return ensemble_prob
@@ -412,11 +412,12 @@ function mcsolve(
         ens_prob_mc.times,
         states,
         expvals_all,
+        ens_prob_mc.kwargs.rng,
         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/smesolve.jl b/src/time_evolution/smesolve.jl
index a9e5f20be..0e2096738 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -271,7 +271,7 @@ function smesolveEnsembleProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
         prob_sme.dimensions,
-        merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3])),
+        merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3], rng = rng)),
     )
 
     return ensemble_prob
@@ -423,6 +423,7 @@ function smesolve(
         ens_prob.times,
         states,
         expvals_all,
+        ens_prob.kwargs.rng,
         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 fd4aa976e..99de46848 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -264,7 +264,7 @@ function ssesolveEnsembleProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
         prob_sme.dimensions,
-        (progr = _output_func[2], channel = _output_func[3]),
+        (progr = _output_func[2], channel = _output_func[3], rng = rng),
     )
 
     return ensemble_prob
@@ -417,6 +417,7 @@ function ssesolve(
         ens_prob.times,
         states,
         expvals_all,
+        ens_prob.kwargs.rng,
         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 98042bb6b..6b3dacf2d 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -103,6 +103,7 @@ A structure storing the results and some information from solving quantum trajec
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
+- `rng::AbstractRNG`: Random number generator for reproducibility.
 - `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.
@@ -127,6 +128,7 @@ struct TimeEvolutionMCSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractArray,Nothing},
+    TR<:AbstractRNG,
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
     AlgT<:OrdinaryDiffEqAlgorithm,
@@ -137,6 +139,7 @@ struct TimeEvolutionMCSol{
     times::TT
     states::TS
     expect::TE
+    rng::TR
     col_times::TJT
     col_which::TJW
     converged::Bool
@@ -173,6 +176,8 @@ A structure storing the results and some information from solving trajectories o
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
+- `rng::AbstractRNG`: Random number generator for reproducibility.
+- `measurement::Union{AbstractArray,Nothing}`: Measurements for each trajectories and stochastic collapse operators (`sc_ops`).
 - `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.
@@ -195,6 +200,7 @@ struct TimeEvolutionStochasticSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractArray,Nothing},
+    TR<:AbstractRNG,
     TEM<:Union{AbstractArray,Nothing},
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
@@ -204,6 +210,7 @@ struct TimeEvolutionStochasticSol{
     times::TT
     states::TS
     expect::TE
+    rng::TR
     measurement::TEM
     converged::Bool
     alg::AlgT
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 8193e067e..89b25f839 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -751,11 +751,9 @@
                 rng = rng,
             )
 
-            rng = MersenneTwister(1234)
-            sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
-            sol_sse2 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
+            sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = sol_mc1.rng)
+            sol_sse2 =
+                ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = sol_sse1.rng)
             sol_sme2 = smesolve(
                 H,
                 ψ0,
@@ -765,14 +763,13 @@
                 ntraj = 50,
                 e_ops = e_ops,
                 progress_bar = Val(false),
-                rng = rng,
+                rng = sol_sme1.rng,
             )
 
-            rng = MersenneTwister(1234)
-            sol_mc3 = mcsolve(H, ψ0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
-            sol_sse3 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
+            sol_mc3 =
+                mcsolve(H, ψ0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = sol_mc1.rng)
+            sol_sse3 =
+                ssesolve(H, ψ0, tlist, c_ops, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = sol_sse1.rng)
             sol_sme3 = smesolve(
                 H,
                 ψ0,
@@ -782,7 +779,7 @@
                 ntraj = 60,
                 e_ops = e_ops,
                 progress_bar = Val(false),
-                rng = rng,
+                rng = sol_sme1.rng,
             )
 
             @test sol_mc1.expect ≈ sol_mc2.expect atol = 1e-10

From 8fe048a4526b2ab7f39b8fc49c2be4898ed817ef Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 13:58:11 +0800
Subject: [PATCH 12/17] update documentation

---
 docs/src/users_guide/steadystate.md           |  2 +-
 .../src/users_guide/time_evolution/mcsolve.md | 11 +++++-----
 .../users_guide/time_evolution/solution.md    | 22 ++++++++++++++++---
 .../users_guide/time_evolution/stochastic.md  |  8 +++++--
 .../time_evolution/time_dependent.md          |  2 +-
 5 files changed, 33 insertions(+), 12 deletions(-)

diff --git a/docs/src/users_guide/steadystate.md b/docs/src/users_guide/steadystate.md
index a74353b90..dccf12122 100644
--- a/docs/src/users_guide/steadystate.md
+++ b/docs/src/users_guide/steadystate.md
@@ -94,7 +94,7 @@ tlist = LinRange(0, 50, 100)
 
 # monte-carlo
 sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=100, progress_bar=false)
-exp_mc = real(sol_mc.expect[1, :])
+exp_mc = real(average_expect(sol_mc)[1, :])
 
 # master eq.
 sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=false)
diff --git a/docs/src/users_guide/time_evolution/mcsolve.md b/docs/src/users_guide/time_evolution/mcsolve.md
index e5a3e0309..6e369291a 100644
--- a/docs/src/users_guide/time_evolution/mcsolve.md
+++ b/docs/src/users_guide/time_evolution/mcsolve.md
@@ -70,6 +70,7 @@ c_ops = [sqrt(0.1) * a]
 e_ops = [a' * a, σm' * σm]
 
 sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops)
+expect_avg = average_expect(sol_500)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -78,8 +79,8 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution (500 trajectories)",
 )
-lines!(ax, times, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
-lines!(ax, times, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
+lines!(ax, times, real(expect_avg[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, times, real(expect_avg[2,:]), label = "atom excitation probability", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -104,9 +105,9 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution",
 )
-lines!(ax, times, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
-lines!(ax, times, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
-lines!(ax, times, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
+lines!(ax, times, real(average_expect(sol_1)[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, times, real(average_expect(sol_10)[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, times, real(average_expect(sol_100)[1,:]), label = "100 trajectories", linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index b6a694837..0a2bf9991 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -89,7 +89,23 @@ 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 expectation values and states for all trajectories at each time points will be stored in these solutions. The index-order of the elements in fields `states` and `expect` are:
 
-- [`TimeEvolutionMCSol`](@ref)
-- [`TimeEvolutionStochasticSol`](@ref)
\ No newline at end of file
+```julia
+sol.states[trajectory][time]
+sol.expect[e_op,trajectory,time]
+```
+
+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. |
+
+Multi-trajectory solutions also keep the random number generator (`rng`) to allow recomputing the results:
+
+```julia
+rng = sol.rng # this can be specified as a keyword argument (`rng = rng`) to the solvers
+```
\ No newline at end of file
diff --git a/docs/src/users_guide/time_evolution/stochastic.md b/docs/src/users_guide/time_evolution/stochastic.md
index 852ae089d..d1820b859 100644
--- a/docs/src/users_guide/time_evolution/stochastic.md
+++ b/docs/src/users_guide/time_evolution/stochastic.md
@@ -111,11 +111,13 @@ sse_sol = ssesolve(
 measurement_avg = sum(sse_sol.measurement, dims=2) / size(sse_sol.measurement, 2)
 measurement_avg = dropdims(measurement_avg, dims=2)
 
+expect_avg = average_expect(sse_sol)
+
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time")
 lines!(ax, tlist[1:end-1], real(measurement_avg[1,:]), label = L"J_x", color = :red, linestyle = :solid)
-lines!(ax, tlist, real(sse_sol.expect[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
+lines!(ax, tlist, real(expect_avg[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
@@ -141,6 +143,8 @@ sme_sol = smesolve(
     store_measurement = Val(true),
 )
 
+expect_avg = average_expect(sme_sol)
+
 measurement_avg = sum(sme_sol.measurement, dims=2) / size(sme_sol.measurement, 2)
 measurement_avg = dropdims(measurement_avg, dims=2)
 
@@ -148,7 +152,7 @@ measurement_avg = dropdims(measurement_avg, dims=2)
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time")
 lines!(ax, tlist[1:end-1], real(measurement_avg[1,:]), label = L"J_x", color = :red, linestyle = :solid)
-lines!(ax, tlist, real(sme_sol.expect[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
+lines!(ax, tlist, real(expect_avg[1,:]),  label = L"\langle x \rangle", color = :black, linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
diff --git a/docs/src/users_guide/time_evolution/time_dependent.md b/docs/src/users_guide/time_evolution/time_dependent.md
index 5c5a0b111..f2e9a46c2 100644
--- a/docs/src/users_guide/time_evolution/time_dependent.md
+++ b/docs/src/users_guide/time_evolution/time_dependent.md
@@ -132,7 +132,7 @@ e_ops = [
 
 # solve dynamics 
 exp_me = mesolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, progress_bar = Val(false)).expect
-exp_mc = mcsolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, ntraj = 100, progress_bar = Val(false)).expect
+exp_mc = average_expect(mcsolve(H_t, ψ0, tlist, c_ops; e_ops = e_ops, ntraj = 100, progress_bar = Val(false)))
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))

From e357aa9a94cc917bafc03992dcbe8dfbb891a470 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 13:58:18 +0800
Subject: [PATCH 13/17] update changelog

---
 CHANGELOG.md | 94 +++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 90 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 230e97fe1..f29ea351c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,10 +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)
 
-- Align statistical analysis methods for multi-trajectory solutions to `QuTiP` [#471]. Remove field `average_expect`, replace field `runs_expect` as `expect`, and introduce the following functions:
-  - `average_states`
-  - `average_expect`
-  - `std_expect`
+- Align statistical analysis methods for multi-trajectory solutions to `QuTiP` [#471].
 
 ## [v0.31.1]
 Release date: 2025-05-16
@@ -147,6 +144,95 @@ Release date: 2024-11-13
 - This is just a demonstration about [`Changelog.jl`](https://github.com/JuliaDocs/Changelog.jl). ([#139], [#306])
 
 
+<!-- Links generated by Changelog.jl -->
+
+[v0.21.4]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.21.4
+[v0.21.5]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.21.5
+[v0.22.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.22.0
+[v0.23.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.0
+[v0.23.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.1
+[v0.24.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.24.0
+[v0.25.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.0
+[v0.25.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.1
+[v0.25.2]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.2
+[v0.26.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.26.0
+[v0.27.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.27.0
+[v0.28.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.28.0
+[v0.29.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.29.0
+[v0.29.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.29.1
+[v0.30.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.30.0
+[v0.30.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.30.1
+[v0.31.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.31.0
+[v0.31.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.31.1
+[#86]: https://github.com/qutip/QuantumToolbox.jl/issues/86
+[#139]: https://github.com/qutip/QuantumToolbox.jl/issues/139
+[#271]: https://github.com/qutip/QuantumToolbox.jl/issues/271
+[#292]: https://github.com/qutip/QuantumToolbox.jl/issues/292
+[#306]: https://github.com/qutip/QuantumToolbox.jl/issues/306
+[#309]: https://github.com/qutip/QuantumToolbox.jl/issues/309
+[#311]: https://github.com/qutip/QuantumToolbox.jl/issues/311
+[#315]: https://github.com/qutip/QuantumToolbox.jl/issues/315
+[#318]: https://github.com/qutip/QuantumToolbox.jl/issues/318
+[#324]: https://github.com/qutip/QuantumToolbox.jl/issues/324
+[#330]: https://github.com/qutip/QuantumToolbox.jl/issues/330
+[#335]: https://github.com/qutip/QuantumToolbox.jl/issues/335
+[#338]: https://github.com/qutip/QuantumToolbox.jl/issues/338
+[#339]: https://github.com/qutip/QuantumToolbox.jl/issues/339
+[#342]: https://github.com/qutip/QuantumToolbox.jl/issues/342
+[#346]: https://github.com/qutip/QuantumToolbox.jl/issues/346
+[#347]: https://github.com/qutip/QuantumToolbox.jl/issues/347
+[#349]: https://github.com/qutip/QuantumToolbox.jl/issues/349
+[#350]: https://github.com/qutip/QuantumToolbox.jl/issues/350
+[#353]: https://github.com/qutip/QuantumToolbox.jl/issues/353
+[#360]: https://github.com/qutip/QuantumToolbox.jl/issues/360
+[#370]: https://github.com/qutip/QuantumToolbox.jl/issues/370
+[#371]: https://github.com/qutip/QuantumToolbox.jl/issues/371
+[#374]: https://github.com/qutip/QuantumToolbox.jl/issues/374
+[#375]: https://github.com/qutip/QuantumToolbox.jl/issues/375
+[#376]: https://github.com/qutip/QuantumToolbox.jl/issues/376
+[#378]: https://github.com/qutip/QuantumToolbox.jl/issues/378
+[#380]: https://github.com/qutip/QuantumToolbox.jl/issues/380
+[#383]: https://github.com/qutip/QuantumToolbox.jl/issues/383
+[#386]: https://github.com/qutip/QuantumToolbox.jl/issues/386
+[#388]: https://github.com/qutip/QuantumToolbox.jl/issues/388
+[#389]: https://github.com/qutip/QuantumToolbox.jl/issues/389
+[#392]: https://github.com/qutip/QuantumToolbox.jl/issues/392
+[#393]: https://github.com/qutip/QuantumToolbox.jl/issues/393
+[#394]: https://github.com/qutip/QuantumToolbox.jl/issues/394
+[#395]: https://github.com/qutip/QuantumToolbox.jl/issues/395
+[#396]: https://github.com/qutip/QuantumToolbox.jl/issues/396
+[#398]: https://github.com/qutip/QuantumToolbox.jl/issues/398
+[#402]: https://github.com/qutip/QuantumToolbox.jl/issues/402
+[#403]: https://github.com/qutip/QuantumToolbox.jl/issues/403
+[#404]: https://github.com/qutip/QuantumToolbox.jl/issues/404
+[#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405
+[#408]: https://github.com/qutip/QuantumToolbox.jl/issues/408
+[#410]: https://github.com/qutip/QuantumToolbox.jl/issues/410
+[#411]: https://github.com/qutip/QuantumToolbox.jl/issues/411
+[#413]: https://github.com/qutip/QuantumToolbox.jl/issues/413
+[#414]: https://github.com/qutip/QuantumToolbox.jl/issues/414
+[#416]: https://github.com/qutip/QuantumToolbox.jl/issues/416
+[#418]: https://github.com/qutip/QuantumToolbox.jl/issues/418
+[#419]: https://github.com/qutip/QuantumToolbox.jl/issues/419
+[#420]: https://github.com/qutip/QuantumToolbox.jl/issues/420
+[#421]: https://github.com/qutip/QuantumToolbox.jl/issues/421
+[#423]: https://github.com/qutip/QuantumToolbox.jl/issues/423
+[#428]: https://github.com/qutip/QuantumToolbox.jl/issues/428
+[#430]: https://github.com/qutip/QuantumToolbox.jl/issues/430
+[#436]: https://github.com/qutip/QuantumToolbox.jl/issues/436
+[#437]: https://github.com/qutip/QuantumToolbox.jl/issues/437
+[#438]: https://github.com/qutip/QuantumToolbox.jl/issues/438
+[#440]: https://github.com/qutip/QuantumToolbox.jl/issues/440
+[#443]: https://github.com/qutip/QuantumToolbox.jl/issues/443
+[#448]: https://github.com/qutip/QuantumToolbox.jl/issues/448
+[#450]: https://github.com/qutip/QuantumToolbox.jl/issues/450
+[#453]: https://github.com/qutip/QuantumToolbox.jl/issues/453
+[#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
+[#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
+[#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
+[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471
+
+
 <!-- Links generated by Changelog.jl -->
 
 [v0.21.4]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.21.4

From 823beb3f4f3e04a7e6ef02c309a0415f993cc405 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 14:03:10 +0800
Subject: [PATCH 14/17] fix changelog

---
 CHANGELOG.md | 93 ++--------------------------------------------------
 1 file changed, 2 insertions(+), 91 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f29ea351c..646573ab0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +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)
 
-- Align statistical analysis methods for multi-trajectory solutions to `QuTiP` [#471].
+- Align statistical analysis methods for multi-trajectory solutions to `QuTiP`. ([#471])
 
 ## [v0.31.1]
 Release date: 2025-05-16
@@ -230,93 +230,4 @@ Release date: 2024-11-13
 [#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
 [#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
 [#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
-[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471
-
-
-<!-- Links generated by Changelog.jl -->
-
-[v0.21.4]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.21.4
-[v0.21.5]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.21.5
-[v0.22.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.22.0
-[v0.23.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.0
-[v0.23.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.1
-[v0.24.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.24.0
-[v0.25.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.0
-[v0.25.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.1
-[v0.25.2]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.25.2
-[v0.26.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.26.0
-[v0.27.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.27.0
-[v0.28.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.28.0
-[v0.29.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.29.0
-[v0.29.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.29.1
-[v0.30.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.30.0
-[v0.30.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.30.1
-[v0.31.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.31.0
-[v0.31.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.31.1
-[#86]: https://github.com/qutip/QuantumToolbox.jl/issues/86
-[#139]: https://github.com/qutip/QuantumToolbox.jl/issues/139
-[#271]: https://github.com/qutip/QuantumToolbox.jl/issues/271
-[#292]: https://github.com/qutip/QuantumToolbox.jl/issues/292
-[#306]: https://github.com/qutip/QuantumToolbox.jl/issues/306
-[#309]: https://github.com/qutip/QuantumToolbox.jl/issues/309
-[#311]: https://github.com/qutip/QuantumToolbox.jl/issues/311
-[#315]: https://github.com/qutip/QuantumToolbox.jl/issues/315
-[#318]: https://github.com/qutip/QuantumToolbox.jl/issues/318
-[#324]: https://github.com/qutip/QuantumToolbox.jl/issues/324
-[#330]: https://github.com/qutip/QuantumToolbox.jl/issues/330
-[#335]: https://github.com/qutip/QuantumToolbox.jl/issues/335
-[#338]: https://github.com/qutip/QuantumToolbox.jl/issues/338
-[#339]: https://github.com/qutip/QuantumToolbox.jl/issues/339
-[#342]: https://github.com/qutip/QuantumToolbox.jl/issues/342
-[#346]: https://github.com/qutip/QuantumToolbox.jl/issues/346
-[#347]: https://github.com/qutip/QuantumToolbox.jl/issues/347
-[#349]: https://github.com/qutip/QuantumToolbox.jl/issues/349
-[#350]: https://github.com/qutip/QuantumToolbox.jl/issues/350
-[#353]: https://github.com/qutip/QuantumToolbox.jl/issues/353
-[#360]: https://github.com/qutip/QuantumToolbox.jl/issues/360
-[#370]: https://github.com/qutip/QuantumToolbox.jl/issues/370
-[#371]: https://github.com/qutip/QuantumToolbox.jl/issues/371
-[#374]: https://github.com/qutip/QuantumToolbox.jl/issues/374
-[#375]: https://github.com/qutip/QuantumToolbox.jl/issues/375
-[#376]: https://github.com/qutip/QuantumToolbox.jl/issues/376
-[#378]: https://github.com/qutip/QuantumToolbox.jl/issues/378
-[#380]: https://github.com/qutip/QuantumToolbox.jl/issues/380
-[#383]: https://github.com/qutip/QuantumToolbox.jl/issues/383
-[#386]: https://github.com/qutip/QuantumToolbox.jl/issues/386
-[#388]: https://github.com/qutip/QuantumToolbox.jl/issues/388
-[#389]: https://github.com/qutip/QuantumToolbox.jl/issues/389
-[#392]: https://github.com/qutip/QuantumToolbox.jl/issues/392
-[#393]: https://github.com/qutip/QuantumToolbox.jl/issues/393
-[#394]: https://github.com/qutip/QuantumToolbox.jl/issues/394
-[#395]: https://github.com/qutip/QuantumToolbox.jl/issues/395
-[#396]: https://github.com/qutip/QuantumToolbox.jl/issues/396
-[#398]: https://github.com/qutip/QuantumToolbox.jl/issues/398
-[#402]: https://github.com/qutip/QuantumToolbox.jl/issues/402
-[#403]: https://github.com/qutip/QuantumToolbox.jl/issues/403
-[#404]: https://github.com/qutip/QuantumToolbox.jl/issues/404
-[#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405
-[#408]: https://github.com/qutip/QuantumToolbox.jl/issues/408
-[#410]: https://github.com/qutip/QuantumToolbox.jl/issues/410
-[#411]: https://github.com/qutip/QuantumToolbox.jl/issues/411
-[#413]: https://github.com/qutip/QuantumToolbox.jl/issues/413
-[#414]: https://github.com/qutip/QuantumToolbox.jl/issues/414
-[#416]: https://github.com/qutip/QuantumToolbox.jl/issues/416
-[#418]: https://github.com/qutip/QuantumToolbox.jl/issues/418
-[#419]: https://github.com/qutip/QuantumToolbox.jl/issues/419
-[#420]: https://github.com/qutip/QuantumToolbox.jl/issues/420
-[#421]: https://github.com/qutip/QuantumToolbox.jl/issues/421
-[#423]: https://github.com/qutip/QuantumToolbox.jl/issues/423
-[#428]: https://github.com/qutip/QuantumToolbox.jl/issues/428
-[#430]: https://github.com/qutip/QuantumToolbox.jl/issues/430
-[#436]: https://github.com/qutip/QuantumToolbox.jl/issues/436
-[#437]: https://github.com/qutip/QuantumToolbox.jl/issues/437
-[#438]: https://github.com/qutip/QuantumToolbox.jl/issues/438
-[#440]: https://github.com/qutip/QuantumToolbox.jl/issues/440
-[#443]: https://github.com/qutip/QuantumToolbox.jl/issues/443
-[#448]: https://github.com/qutip/QuantumToolbox.jl/issues/448
-[#450]: https://github.com/qutip/QuantumToolbox.jl/issues/450
-[#453]: https://github.com/qutip/QuantumToolbox.jl/issues/453
-[#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
-[#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
-[#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
-[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471
+[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471
\ No newline at end of file

From e8c16ebe4245bbded47ee1b7b3042b2daa8b80da Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 17:11:12 +0800
Subject: [PATCH 15/17] fix rng reproducibility

---
 docs/src/users_guide/time_evolution/solution.md | 2 +-
 src/time_evolution/mcsolve.jl                   | 7 +++++--
 src/time_evolution/smesolve.jl                  | 6 ++++--
 src/time_evolution/ssesolve.jl                  | 6 ++++--
 src/time_evolution/time_evolution.jl            | 4 ++--
 5 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index 0a2bf9991..5c260a51e 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -104,7 +104,7 @@ We also provide the following functions for statistical analysis of multi-trajec
 | [`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. |
 
-Multi-trajectory solutions also keep the random number generator (`rng`) to allow recomputing the results:
+Multi-trajectory solutions also keep the initial random number generator (`rng`) to allow recomputing the results:
 
 ```julia
 rng = sol.rng # this can be specified as a keyword argument (`rng = rng`) to the solvers
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 29db4c3d0..d9a41bc1a 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -235,7 +235,10 @@ function mcsolveEnsembleProblem(
     output_func::Union{Tuple,Nothing} = nothing,
     kwargs...,
 ) where {TJC<:LindbladJumpCallbackType}
-    _prob_func = isnothing(prob_func) ? _ensemble_dispatch_prob_func(rng, ntraj, tlist, _mcsolve_prob_func) : prob_func
+    copied_rng = copy(rng) # use the copied rng, and keep the initial one to pass it directly to solution later
+
+    _prob_func =
+        isnothing(prob_func) ? _ensemble_dispatch_prob_func(copied_rng, ntraj, tlist, _mcsolve_prob_func) : prob_func
     _output_func =
         output_func isa Nothing ?
         _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _mcsolve_output_func) : output_func
@@ -247,7 +250,7 @@ function mcsolveEnsembleProblem(
         c_ops;
         e_ops = e_ops,
         params = params,
-        rng = rng,
+        rng = copied_rng,
         jump_callback = jump_callback,
         kwargs...,
     )
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 0e2096738..481f68075 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -239,10 +239,12 @@ function smesolveEnsembleProblem(
     store_measurement::Union{Val,Bool} = Val(false),
     kwargs...,
 ) where {StateOpType<:Union{Ket,Operator,OperatorKet}}
+    copied_rng = copy(rng) # use the copied rng, and keep the initial one to pass it directly to solution later
+
     _prob_func =
         isnothing(prob_func) ?
         _ensemble_dispatch_prob_func(
-            rng,
+            copied_rng,
             ntraj,
             tlist,
             _stochastic_prob_func;
@@ -261,7 +263,7 @@ function smesolveEnsembleProblem(
         sc_ops;
         e_ops = e_ops,
         params = params,
-        rng = rng,
+        rng = copied_rng,
         progress_bar = Val(false),
         store_measurement = makeVal(store_measurement),
         kwargs...,
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 99de46848..e3285f859 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -233,10 +233,12 @@ function ssesolveEnsembleProblem(
     store_measurement::Union{Val,Bool} = Val(false),
     kwargs...,
 )
+    copied_rng = copy(rng) # use the copied rng, and keep the initial one to pass it directly to solution later
+
     _prob_func =
         isnothing(prob_func) ?
         _ensemble_dispatch_prob_func(
-            rng,
+            copied_rng,
             ntraj,
             tlist,
             _stochastic_prob_func;
@@ -254,7 +256,7 @@ function ssesolveEnsembleProblem(
         sc_ops;
         e_ops = e_ops,
         params = params,
-        rng = rng,
+        rng = copied_rng,
         progress_bar = Val(false),
         store_measurement = makeVal(store_measurement),
         kwargs...,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 6b3dacf2d..12afbaabb 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -103,7 +103,7 @@ A structure storing the results and some information from solving quantum trajec
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
+- `rng::AbstractRNG`: Initial random number generator for reproducibility.
 - `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.
@@ -176,7 +176,7 @@ A structure storing the results and some information from solving trajectories o
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
-- `rng::AbstractRNG`: Random number generator for reproducibility.
+- `rng::AbstractRNG`: Initial random number generator for reproducibility.
 - `measurement::Union{AbstractArray,Nothing}`: Measurements for each trajectories and stochastic collapse operators (`sc_ops`).
 - `converged::Bool`: Whether the solution is converged or not.
 - `alg`: The algorithm which is used during the solving process.

From f72e9740eee5d77c8b4b7a29e31a09955668ad37 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sat, 24 May 2025 17:11:44 +0800
Subject: [PATCH 16/17] update changelog

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 646573ab0..ab2739880 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -230,4 +230,4 @@ Release date: 2024-11-13
 [#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
 [#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
 [#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
-[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471
\ No newline at end of file
+[#471]: https://github.com/qutip/QuantumToolbox.jl/issues/471

From d49c89c8b6a8045ed1c395dccc5ec582be991b17 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <44385685+ytdHuang@users.noreply.github.com>
Date: Mon, 26 May 2025 10:24:23 +0800
Subject: [PATCH 17/17] Update CHANGELOG.md

---
 CHANGELOG.md | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index ab2739880..fb1fef78c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Align statistical analysis methods for multi-trajectory solutions to `QuTiP`. ([#471])
+  - Breaking Changes:
+    - the origin field `expect` in multi-trajectory solutions now stores the expectation values for all trajectories at each time points.
+    - remove field `average_expect` from multi-trajectory solutions
+    - remove field `runs_expect` from multi-trajectory solutions
+    - add field `rng` (random number generator) to multi-trajectory solutions for reproducibility
+  - New statistical analysis functions:
+    - `average_states`
+    - `average_expect`
+    - `std_expect`
 
 ## [v0.31.1]
 Release date: 2025-05-16