From f0a262f2ce5c74613601e5f474c3e801114aca0b Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 15 Jul 2025 23:57:17 +0800
Subject: [PATCH 1/5] fix spelling typo

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

diff --git a/src/correlations.jl b/src/correlations.jl
index 3eae092b5..a8ac92f27 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -19,7 +19,7 @@ end
         C::QuantumObject;
         kwargs...)
 
-Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 """
@@ -96,7 +96,7 @@ end
         reverse::Bool=false,
         kwargs...)
 
-Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 

From 7b57f46e8ff0abd2e35eae1a96061b803fda1255 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Wed, 16 Jul 2025 00:30:06 +0800
Subject: [PATCH 2/5] Store both `tlist` and `times` in time evolution
 solutions

---
 .../src/users_guide/time_evolution/mcsolve.md | 24 ++++-----
 .../src/users_guide/time_evolution/mesolve.md | 11 ++--
 .../src/users_guide/time_evolution/sesolve.md | 13 +++--
 .../users_guide/time_evolution/solution.md    | 24 ++++++---
 src/spectrum.jl                               |  2 +-
 src/time_evolution/brmesolve.jl               |  2 +-
 src/time_evolution/lr_mesolve.jl              | 26 +++++----
 src/time_evolution/mcsolve.jl                 | 11 ++--
 src/time_evolution/mesolve.jl                 |  7 +--
 src/time_evolution/sesolve.jl                 |  7 +--
 src/time_evolution/smesolve.jl                | 11 ++--
 src/time_evolution/ssesolve.jl                | 11 ++--
 src/time_evolution/time_evolution.jl          | 53 +++++++++++--------
 .../time_evolution_dynamical.jl               |  3 +-
 test/core-test/brmesolve.jl                   |  6 +--
 test/core-test/time_evolution.jl              | 26 ++++++---
 16 files changed, 138 insertions(+), 99 deletions(-)

diff --git a/docs/src/users_guide/time_evolution/mcsolve.md b/docs/src/users_guide/time_evolution/mcsolve.md
index e5a3e0309..4f9aaa723 100644
--- a/docs/src/users_guide/time_evolution/mcsolve.md
+++ b/docs/src/users_guide/time_evolution/mcsolve.md
@@ -59,7 +59,7 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 ```
 
 ```@example mcsolve
-times = LinRange(0.0, 10.0, 200)
+tlist = LinRange(0.0, 10.0, 200)
 
 ψ0 = tensor(fock(2, 0), fock(10, 8))
 a  = tensor(qeye(2), destroy(10))
@@ -69,7 +69,7 @@ H = 2 * π * a' * a + 2 * π * σm' * σm + 2 * π * 0.25 * (σm * a' + σm' * a
 c_ops = [sqrt(0.1) * a]
 e_ops = [a' * a, σm' * σm]
 
-sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops)
+sol_500 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -78,8 +78,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, tlist, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, tlist, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -93,9 +93,9 @@ We can also change the number of trajectories (`ntraj`). This can be used to exp
 ```@example mcsolve
 e_ops = [a' * a]
 
-sol_1   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 1)
-sol_10  = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 10)
-sol_100 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 100)
+sol_1   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 1)
+sol_10  = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 10)
+sol_100 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 100)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -104,9 +104,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, tlist, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, tlist, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, tlist, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
@@ -126,8 +126,8 @@ Monte Carlo evolutions often need hundreds of trajectories to obtain sufficient
     See the [documentation of `DifferentialEquations.jl`](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) for more details. Also, see Julia's documentation for more details about multithreading and adding more processes.
 
 ```julia
-sol_serial   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
-sol_parallel = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
+sol_serial   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
+sol_parallel = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
 ```
 
 !!! tip "Parallelization on a Cluster"
diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md
index 9537d3bf8..dc5b3eb13 100644
--- a/docs/src/users_guide/time_evolution/mesolve.md
+++ b/docs/src/users_guide/time_evolution/mesolve.md
@@ -128,15 +128,14 @@ sol = mesolve(H, ψ0, tlist, c_ops, e_ops = [sigmaz(), sigmay()])
 We can therefore plot the expectation values:
 
 ```@example mesolve
-times = sol.times
 expt_z = real(sol.expect[1,:])
 expt_y = real(sol.expect[2,:])
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -210,8 +209,6 @@ tlist = LinRange(0.0, 10.0, 200)
 L = liouvillian(H_a + H_c + H_I, c_ops)
 sol = mesolve(L, ψ0, tlist, e_ops=[σm' * σm, a' * a])
 
-times = sol.times
-
 # expectation value of Number operator
 N_atom   = real(sol.expect[1,:])
 N_cavity = real(sol.expect[2,:])
@@ -219,8 +216,8 @@ N_cavity = real(sol.expect[2,:])
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid)
-lines!(ax, times, N_cavity, label = "cavity photon number", linestyle = :dash)
+lines!(ax, tlist, N_atom, label = "atom excitation probability", linestyle = :solid)
+lines!(ax, tlist, N_cavity, label = "cavity photon number", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index d3ba8a21c..9dee3c005 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -59,8 +59,8 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])
     Here, we call [`sesolve`](@ref) directly instead of pre-defining [`sesolveProblem`](@ref) first (as shown previously).
 
 ```@example sesolve
-times = sol.times
-print(size(times))
+println(size(sol.tlist)) # time points corresponds to stored expectation values 
+println(size(sol.times)) # time points corresponds to stored states
 ```
 
 ```@example sesolve
@@ -77,8 +77,8 @@ expt_y = real(expt[2,:])
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rb)
 
@@ -91,6 +91,9 @@ If the keyword argument `e_ops` is not specified (or given as an empty `Vector`)
 tlist = [0, 10]
 sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
 
+println(size(sol.tlist))
+println(size(sol.times))
+
 sol.states
 ```
 
@@ -104,9 +107,11 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)
 ```
 
 ```@example sesolve
+println(size(sol.tlist))
 sol.expect
 ```
 
 ```@example sesolve
+println(size(sol.times))
 sol.states
 ```
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index b6a694837..f49b9f483 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -12,9 +12,10 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 
 | **Fields (Attributes)** | **Description** |
 |:------------------------|:----------------|
-| `sol.times` | The time list of the evolution. |
-| `sol.states` | The list of result states. |
-| `sol.expect` | The expectation values corresponding to each time point in `sol.times`. |
+| `sol.times` | The list of time points at which the states are stored during the evolution. |
+| `sol.tlist` | The list of time points at which the expectation values are stored during the evolution. |
+| `sol.states` | The list of result states corresponding to each time point in `times`. |
+| `sol.expect` | The expectation values corresponding to each time point in `sol.tlist`. |
 | `sol.alg` | The algorithm which is used during the solving process. |
 | `sol.abstol` | The absolute tolerance which is used during the solving process. |
 | `sol.reltol` | The relative tolerance which is used during the solving process. |
@@ -54,10 +55,10 @@ nothing # hide
 
 Recall that `Julia` uses `Fortran`-style indexing that begins with one (i.e., `[1,:]` represents the 1-st observable, where `:` represents all values corresponding to `tlist`).
 
-Together with the array of times at which these expectation values are calculated:
+Together with the list of time points at which these expectation values are calculated:
 
 ```@example TE-solution
-times = sol.times
+tlist = sol.tlist
 nothing # hide
 ```
 
@@ -67,9 +68,9 @@ we can plot the resulting expectation values:
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = L"t")
-lines!(ax, times, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle")
-lines!(ax, times, expt2, label = L"\langle 1 | \rho(t) | 1 \rangle")
-lines!(ax, times, expt3, label = L"\langle 0 | \rho(t) | 1 \rangle")
+lines!(ax, tlist, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle")
+lines!(ax, tlist, expt2, label = L"\langle 1 | \rho(t) | 1 \rangle")
+lines!(ax, tlist, expt3, label = L"\langle 0 | \rho(t) | 1 \rangle")
 
 ylims!(ax, (-0.5, 1.0))
 axislegend(ax, position = :lb)
@@ -83,6 +84,13 @@ State vectors, or density matrices, are accessed in a similar manner:
 sol.states
 ```
 
+Together with the list of time points at which these states are stored:
+
+```@example TE-solution
+times = sol.times
+nothing # hide
+```
+
 Here, the solution contains only one (final) state. Because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). One can also specify `e_ops` and `saveat` separately.
 
 Some other solvers can have other output.
diff --git a/src/spectrum.jl b/src/spectrum.jl
index b97488575..9a6313ad4 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -313,7 +313,7 @@ end
 Calculate the power spectrum corresponding to a two-time correlation function using fast Fourier transform (FFT).
 
 # Parameters
-- `tlist::AbstractVector`: List of times at which the two-time correlation function is given.
+- `tlist::AbstractVector`: List of time points at which the two-time correlation function is given.
 - `corr::AbstractVector`: List of two-time correlations corresponding to the given time point in `tlist`.
 - `inverse::Bool`: Whether to use the inverse Fourier transform or not. Default to `false`.
 
diff --git a/src/time_evolution/brmesolve.jl b/src/time_evolution/brmesolve.jl
index 05a47a2d3..0fc04d035 100644
--- a/src/time_evolution/brmesolve.jl
+++ b/src/time_evolution/brmesolve.jl
@@ -158,7 +158,7 @@ Solves for the dynamics of a system using the Bloch-Redfield master equation, gi
 
 - `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
 - `ψ0`:  Initial state of the system $|\psi(0)\rangle$. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `a_ops`: Nested list with each element is a `Tuple` of operator-function pairs `(a_op, spectra)`, and the coupling [`Operator`](@ref) `a_op` must be hermitian with corresponding `spectra` being a `Function` of transition energy
 - `c_ops`: List of collapse operators corresponding to Lindblad dissipators
 - `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index bc166a540..15cb4fc02 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -7,10 +7,11 @@ A structure storing the results and some information from solving low-rank maste
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{QuantumObject}`: The list of result states.
-- `expect::Matrix`: The expectation values corresponding to each time point in `times`.
-- `fexpect::Matrix`: The function values at each time point.
+- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
+- `expect::Matrix`: The expectation values corresponding to each time point in `tlist`.
+- `fexpect::Matrix`: The function values corresponding to each time point in `tlist`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -19,7 +20,8 @@ A structure storing the results and some information from solving low-rank maste
 - `B::Vector{QuantumObject}`: The `B` matrix of the low-rank algorithm at each time point.
 """
 struct TimeEvolutionLRSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     RetT<:Enum,
@@ -29,7 +31,8 @@ struct TimeEvolutionLRSol{
     TSZB<:AbstractVector,
     TM<:Vector{<:Integer},
 }
-    times::TT
+    times::TT1
+    tlist::TT2
     states::TS
     expect::TE
     fexpect::TE
@@ -144,19 +147,19 @@ function _periodicsave_func(integrator)
     return u_modified!(integrator, false)
 end
 
-_save_control_lr_mesolve(u, t, integrator) = t in integrator.p.times
+_save_control_lr_mesolve(u, t, integrator) = t in integrator.p.tlist
 
 function _save_affect_lr_mesolve!(integrator)
     ip = integrator.p
     N, M = ip.N, ip.M
-    idx = select(integrator.t, ip.times)
+    idx = select(integrator.t, ip.tlist)
 
     @views z = reshape(integrator.u[1:(N*M)], N, M)
     @views B = reshape(integrator.u[(N*M+1):end], M, M)
     _calculate_expectation!(ip, z, B, idx)
 
     if integrator.p.opt.progress
-        print("\rProgress: $(round(Int, 100*idx/length(ip.times)))%")
+        print("\rProgress: $(round(Int, 100*idx/length(ip.tlist)))%")
         flush(stdout)
     end
     return u_modified!(integrator, false)
@@ -430,7 +433,7 @@ function lr_mesolveProblem(
         e_ops = e_ops,
         f_ops = f_ops,
         opt = opt,
-        times = t_l,
+        tlist = t_l,
         expvals = expvals,
         funvals = funvals,
         Ml = Ml,
@@ -537,7 +540,7 @@ function lr_mesolve(
 end
 
 function lr_mesolve(prob::ODEProblem; kwargs...)
-    sol = solve(prob, prob.p.opt.alg, tstops = prob.p.times)
+    sol = solve(prob, prob.p.opt.alg, tstops = prob.p.tlist)
     prob.p.opt.progress && print("\n")
 
     N = prob.p.N
@@ -550,6 +553,7 @@ function lr_mesolve(prob::ODEProblem; kwargs...)
 
     return TimeEvolutionLRSol(
         sol.t,
+        prob.p.tlist,
         ρt,
         prob.p.expvals,
         prob.p.funvals,
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 7db7048ff..db01499c6 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -89,7 +89,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -195,7 +195,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -254,7 +254,7 @@ function mcsolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
-        prob_mc.times,
+        prob_mc.tlist,
         prob_mc.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
@@ -320,7 +320,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm to use for the ODE solver. Default to `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -411,7 +411,8 @@ function mcsolve(
 
     return TimeEvolutionMCSol(
         ntraj,
-        ens_prob_mc.times,
+        _sol_1.t,
+        ens_prob_mc.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 32305cb1b..e61660c8c 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -34,7 +34,7 @@ where
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -132,7 +132,7 @@ where
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm for the ODE solver. The default value is `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -205,7 +205,8 @@ function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
     end
 
     return TimeEvolutionSol(
-        prob.times,
+        sol.t,
+        prob.tlist,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index c5607f231..5a9bebf3a 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -30,7 +30,7 @@ Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -104,7 +104,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -155,7 +155,8 @@ function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
     ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.dimensions), sol.u)
 
     return TimeEvolutionSol(
-        prob.times,
+        sol.t,
+        prob.tlist,
         ψt,
         _get_expvals(sol, SaveFuncSESolve),
         sol.retcode,
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 710029197..1a782b981 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -49,7 +49,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -194,7 +194,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -269,7 +269,7 @@ function smesolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
-        prob_sme.times,
+        prob_sme.tlist,
         prob_sme.dimensions,
         merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3])),
     )
@@ -321,7 +321,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise).
@@ -424,7 +424,8 @@ function smesolve(
 
     return TimeEvolutionStochasticSol(
         ntraj,
-        ens_prob.times,
+        _sol_1.t,
+        ens_prob.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 6945dc81e..a3d394f34 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -51,7 +51,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: `NullParameters` of parameters to pass to the solver.
@@ -189,7 +189,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and  ``dW_n(t)`` is t
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: `NullParameters` of parameters to pass to the solver.
@@ -262,7 +262,7 @@ function ssesolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
-        prob_sme.times,
+        prob_sme.tlist,
         prob_sme.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
@@ -317,7 +317,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise).
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -418,7 +418,8 @@ function ssesolve(
 
     return TimeEvolutionStochasticSol(
         ntraj,
-        ens_prob.times,
+        _sol_1.t,
+        ens_prob.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 0116ff398..eec44389e 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -14,7 +14,7 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 # Fields (Attributes)
 
 - `prob::AbstractSciMLProblem`: The `ODEProblem` of the time evolution.
-- `times::AbstractVector`: The time list of the evolution.
+- `tlist::AbstractVector`: The time list of the evolution.
 - `dimensions::AbstractDimensions`: The dimensions of the Hilbert space.
 - `kwargs::KWT`: Generic keyword arguments.
 
@@ -23,7 +23,7 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 """
 struct TimeEvolutionProblem{PT<:AbstractSciMLProblem,TT<:AbstractVector,DT<:AbstractDimensions,KWT}
     prob::PT
-    times::TT
+    tlist::TT
     dimensions::DT
     kwargs::KWT
 end
@@ -37,7 +37,7 @@ function Base.getproperty(prob::TimeEvolutionProblem, key::Symbol)
     end
 end
 
-TimeEvolutionProblem(prob, times, dims) = TimeEvolutionProblem(prob, times, dims, nothing)
+TimeEvolutionProblem(prob, tlist, dims) = TimeEvolutionProblem(prob, tlist, dims, nothing)
 
 @doc raw"""
     struct TimeEvolutionSol
@@ -46,16 +46,18 @@ A structure storing the results and some information from solving time evolution
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{QuantumObject}`: The list of result states.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `times`.
+- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
+- `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `tlist`.
 - `retcode`: The return code from the solver.
 - `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.
 """
 struct TimeEvolutionSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     RETT<:Enum,
@@ -63,7 +65,8 @@ struct TimeEvolutionSol{
     AT<:Real,
     RT<:Real,
 }
-    times::TT
+    times::TT1
+    tlist::TT2
     states::TS
     expect::TE
     retcode::RETT
@@ -96,11 +99,12 @@ A structure storing the results and some information from solving quantum trajec
 # Fields (Attributes)
 
 - `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`
+- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
+- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
+- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
+- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `tlist`
 - `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.
@@ -109,7 +113,8 @@ A structure storing the results and some information from solving quantum trajec
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionMCSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     TEA<:Union{AbstractArray,Nothing},
@@ -120,7 +125,8 @@ struct TimeEvolutionMCSol{
     RT<:Real,
 }
     ntraj::Int
-    times::TT
+    times::TT1
+    tlist::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
@@ -158,18 +164,20 @@ A structure storing the results and some information from solving trajectories o
 # Fields (Attributes)
 
 - `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`
+- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
+- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
+- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
+- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `tlist`
 - `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.
 """
 struct TimeEvolutionStochasticSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     TEA<:Union{AbstractArray,Nothing},
@@ -179,7 +187,8 @@ struct TimeEvolutionStochasticSol{
     RT<:Real,
 }
     ntraj::Int
-    times::TT
+    times::TT1
+    tlist::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 2e07047da..1b1ef95db 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -243,7 +243,8 @@ function dfd_mesolve(
     end
 
     return TimeEvolutionSol(
-        dfd_prob.times,
+        sol.t,
+        dfd_prob.tlist,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/test/core-test/brmesolve.jl b/test/core-test/brmesolve.jl
index 68cfad8c1..63b6ba02c 100644
--- a/test/core-test/brmesolve.jl
+++ b/test/core-test/brmesolve.jl
@@ -104,11 +104,11 @@ end
     e_ops = pauli_vectors
     H = δ * 0.5 * sigmax() + ϵ * 0.5 * sigmaz()
     ψ0 = unit(2basis(2, 0) + basis(2, 1))
-    times = LinRange(0, 10, 100)
+    tlist = LinRange(0, 10, 100)
 
     for (me_c_ops, brme_c_ops, brme_a_ops) in arg_sets
-        me = mesolve(H, ψ0, times, me_c_ops, e_ops = e_ops, progress_bar = Val(false))
-        brme = brmesolve(H, ψ0, times, brme_a_ops, brme_c_ops, e_ops = e_ops, progress_bar = Val(false))
+        me = mesolve(H, ψ0, tlist, me_c_ops, e_ops = e_ops, progress_bar = Val(false))
+        brme = brmesolve(H, ψ0, tlist, brme_a_ops, brme_c_ops, e_ops = e_ops, progress_bar = Val(false))
 
         @test all(me.expect .== brme.expect)
     end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 71642227c..72a57afd9 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -84,13 +84,16 @@ end
 
     @test prob.prob.f.f isa MatrixOperator
     @test sum(abs.(sol.expect[1, :] .- amp_rabi .* sin.(Ω_rabi * tlist) .^ 2)) / length(tlist) < 0.1
-    @test length(sol.times) == length(tlist)
+    @test length(sol.times) == 1
+    @test length(sol.tlist) == length(tlist)
     @test length(sol.states) == 1
     @test size(sol.expect) == (length(e_ops), length(tlist))
     @test length(sol2.times) == length(tlist)
+    @test length(sol2.tlist) == length(tlist)
     @test length(sol2.states) == length(tlist)
     @test sol2.expect === nothing
-    @test length(sol3.times) == length(tlist)
+    @test length(sol3.times) == length(saveat)
+    @test length(sol3.tlist) == length(tlist)
     @test length(sol3.states) == length(saveat)
     @test size(sol3.expect) == (length(e_ops), length(tlist))
     @test sol.expect[1, saveat_idxs] ≈ expect(e_ops[1], sol3.states) atol = 1e-6
@@ -167,13 +170,16 @@ end
 
     @test TESetup.prob_me.prob.f.f isa MatrixOperator
     @test isket(sol_me5.states[1])
-    @test length(sol_me.times) == length(tlist)
+    @test length(sol_me.times) == 1
+    @test length(sol_me.tlist) == length(tlist)
     @test length(sol_me.states) == 1
     @test size(sol_me.expect) == (length(e_ops), length(tlist))
     @test length(sol_me2.times) == length(tlist)
+    @test length(sol_me2.tlist) == length(tlist)
     @test length(sol_me2.states) == length(tlist)
     @test sol_me2.expect === nothing
-    @test length(sol_me3.times) == length(tlist)
+    @test length(sol_me3.times) == length(saveat)
+    @test length(sol_me3.tlist) == length(tlist)
     @test length(sol_me3.states) == length(saveat)
     @test size(sol_me3.expect) == (length(e_ops), length(tlist))
     @test sol_me3.expect[1, TESetup.saveat_idxs] ≈ expect(e_ops[1], sol_me3.states) atol = 1e-6
@@ -288,9 +294,11 @@ end
     @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 length(sol_mc.times) == length(tlist)
+    @test length(sol_mc.times) == 1
+    @test length(sol_mc.tlist) == length(tlist)
     @test size(sol_mc.expect) == (length(e_ops), length(tlist))
-    @test length(sol_mc_states.times) == length(tlist)
+    @test length(sol_mc_states.times) == length(saveat)
+    @test length(sol_mc_states.tlist) == length(tlist)
     @test sol_mc_states.expect === nothing
 
     sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
@@ -390,7 +398,8 @@ 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) == 1
+    @test length(sol_sse.tlist) == length(tlist)
     @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)
@@ -524,7 +533,8 @@ end
 
     @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 length(sol_sme.times) == length(tlist)
+    @test length(sol_sme.times) == 1
+    @test length(sol_sme.tlist) == length(tlist)
     @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)

From f69b9aff67710ca55d4e1497d9c43161a1ed6dee Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Wed, 16 Jul 2025 00:41:09 +0800
Subject: [PATCH 3/5] update changelog

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

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9e1d0d804..7aa0469c7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Implement `EnrSpace` and corresponding functionality. ([#500])
 - Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
+- Store both `tlist` and `times` in time evolution solutions. ([#506], [#504])
 
 ## [v0.32.1]
 Release date: 2025-06-24
@@ -258,3 +259,5 @@ Release date: 2024-11-13
 [#494]: https://github.com/qutip/QuantumToolbox.jl/issues/494
 [#500]: https://github.com/qutip/QuantumToolbox.jl/issues/500
 [#501]: https://github.com/qutip/QuantumToolbox.jl/issues/501
+[#504]: https://github.com/qutip/QuantumToolbox.jl/issues/504
+[#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506

From 05b42682c377c355863c94a1a6a3db707f1a4c50 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Wed, 16 Jul 2025 01:03:34 +0800
Subject: [PATCH 4/5] fix typos

---
 docs/src/users_guide/time_evolution/solution.md | 2 +-
 src/time_evolution/lr_mesolve.jl                | 2 +-
 src/time_evolution/time_evolution.jl            | 6 +++---
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index f49b9f483..fe50da434 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -13,7 +13,7 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 | **Fields (Attributes)** | **Description** |
 |:------------------------|:----------------|
 | `sol.times` | The list of time points at which the states are stored during the evolution. |
-| `sol.tlist` | The list of time points at which the expectation values are stored during the evolution. |
+| `sol.tlist` | The list of time points at which the expectation values are calculated during the evolution. |
 | `sol.states` | The list of result states corresponding to each time point in `times`. |
 | `sol.expect` | The expectation values corresponding to each time point in `sol.tlist`. |
 | `sol.alg` | The algorithm which is used during the solving process. |
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index 15cb4fc02..d0ca3b56e 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -8,7 +8,7 @@ A structure storing the results and some information from solving low-rank maste
 # Fields (Attributes)
 
 - `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
 - `expect::Matrix`: The expectation values corresponding to each time point in `tlist`.
 - `fexpect::Matrix`: The function values corresponding to each time point in `tlist`.
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index eec44389e..4187c6f95 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -47,7 +47,7 @@ A structure storing the results and some information from solving time evolution
 # Fields (Attributes)
 
 - `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `tlist`.
 - `retcode`: The return code from the solver.
@@ -100,7 +100,7 @@ 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 states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
 - `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
@@ -165,7 +165,7 @@ 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 states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are stored during the evolution.
+- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
 - `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.

From d0bb36daa81a2e75a5174ba3d2b91d68c2dd6c86 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Mon, 21 Jul 2025 23:54:02 +0800
Subject: [PATCH 5/5] rename field `tlist` as `times`, and `times` as
 `times_states`

---
 CHANGELOG.md                                  |  2 +-
 .../src/users_guide/time_evolution/sesolve.md | 10 ++---
 .../users_guide/time_evolution/solution.md    | 18 ++++----
 src/time_evolution/lr_mesolve.jl              | 24 +++++-----
 src/time_evolution/mcsolve.jl                 |  4 +-
 src/time_evolution/mesolve.jl                 |  2 +-
 src/time_evolution/sesolve.jl                 |  2 +-
 src/time_evolution/smesolve.jl                |  4 +-
 src/time_evolution/ssesolve.jl                |  4 +-
 src/time_evolution/time_evolution.jl          | 44 +++++++++----------
 .../time_evolution_dynamical.jl               |  2 +-
 test/core-test/time_evolution.jl              | 36 +++++++--------
 12 files changed, 76 insertions(+), 76 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8fdc0b4bb..d5c2c198f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Implement `EnrSpace` and corresponding functionality. ([#500])
 - Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
-- Store both `tlist` and `times` in time evolution solutions. ([#506], [#504])
+- Store both `times` and `times_states` in time evolution solutions. ([#506], [#504])
 - Fix errors in `Julia v1.12`. ([#507])
 
 ## [v0.32.1]
diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index 9dee3c005..70372a7f0 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -59,8 +59,8 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])
     Here, we call [`sesolve`](@ref) directly instead of pre-defining [`sesolveProblem`](@ref) first (as shown previously).
 
 ```@example sesolve
-println(size(sol.tlist)) # time points corresponds to stored expectation values 
-println(size(sol.times)) # time points corresponds to stored states
+println(size(sol.times)) # time points corresponds to stored expectation values 
+println(size(sol.times_states)) # time points corresponds to stored states
 ```
 
 ```@example sesolve
@@ -91,8 +91,8 @@ If the keyword argument `e_ops` is not specified (or given as an empty `Vector`)
 tlist = [0, 10]
 sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
 
-println(size(sol.tlist))
 println(size(sol.times))
+println(size(sol.times_states))
 
 sol.states
 ```
@@ -107,11 +107,11 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)
 ```
 
 ```@example sesolve
-println(size(sol.tlist))
+println(size(sol.times))
 sol.expect
 ```
 
 ```@example sesolve
-println(size(sol.times))
+println(size(sol.times_states))
 sol.states
 ```
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index fe50da434..94609f18a 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -12,10 +12,10 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 
 | **Fields (Attributes)** | **Description** |
 |:------------------------|:----------------|
-| `sol.times` | The list of time points at which the states are stored during the evolution. |
-| `sol.tlist` | The list of time points at which the expectation values are calculated during the evolution. |
-| `sol.states` | The list of result states corresponding to each time point in `times`. |
-| `sol.expect` | The expectation values corresponding to each time point in `sol.tlist`. |
+| `sol.times` | The list of time points at which the expectation values are calculated during the evolution. |
+| `sol.times_states` | The list of time points at which the states are stored during the evolution. |
+| `sol.states` | The list of result states corresponding to each time point in `sol.times_states`. |
+| `sol.expect` | The expectation values corresponding to each time point in `sol.times`. |
 | `sol.alg` | The algorithm which is used during the solving process. |
 | `sol.abstol` | The absolute tolerance which is used during the solving process. |
 | `sol.reltol` | The relative tolerance which is used during the solving process. |
@@ -58,7 +58,7 @@ Recall that `Julia` uses `Fortran`-style indexing that begins with one (i.e., `[
 Together with the list of time points at which these expectation values are calculated:
 
 ```@example TE-solution
-tlist = sol.tlist
+times = sol.times
 nothing # hide
 ```
 
@@ -68,9 +68,9 @@ we can plot the resulting expectation values:
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = L"t")
-lines!(ax, tlist, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle")
-lines!(ax, tlist, expt2, label = L"\langle 1 | \rho(t) | 1 \rangle")
-lines!(ax, tlist, expt3, label = L"\langle 0 | \rho(t) | 1 \rangle")
+lines!(ax, times, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle")
+lines!(ax, times, expt2, label = L"\langle 1 | \rho(t) | 1 \rangle")
+lines!(ax, times, expt3, label = L"\langle 0 | \rho(t) | 1 \rangle")
 
 ylims!(ax, (-0.5, 1.0))
 axislegend(ax, position = :lb)
@@ -87,7 +87,7 @@ sol.states
 Together with the list of time points at which these states are stored:
 
 ```@example TE-solution
-times = sol.times
+times = sol.times_states
 nothing # hide
 ```
 
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index d0ca3b56e..1c73010e1 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -7,11 +7,11 @@ A structure storing the results and some information from solving low-rank maste
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
-- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
-- `expect::Matrix`: The expectation values corresponding to each time point in `tlist`.
-- `fexpect::Matrix`: The function values corresponding to each time point in `tlist`.
+- `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{QuantumObject}`: The list of result states corresponding to each time point in `times_states`.
+- `expect::Matrix`: The expectation values corresponding to each time point in `times`.
+- `fexpect::Matrix`: The function values corresponding to each time point in `times`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -32,7 +32,7 @@ struct TimeEvolutionLRSol{
     TM<:Vector{<:Integer},
 }
     times::TT1
-    tlist::TT2
+    times_states::TT2
     states::TS
     expect::TE
     fexpect::TE
@@ -147,19 +147,19 @@ function _periodicsave_func(integrator)
     return u_modified!(integrator, false)
 end
 
-_save_control_lr_mesolve(u, t, integrator) = t in integrator.p.tlist
+_save_control_lr_mesolve(u, t, integrator) = t in integrator.p.times
 
 function _save_affect_lr_mesolve!(integrator)
     ip = integrator.p
     N, M = ip.N, ip.M
-    idx = select(integrator.t, ip.tlist)
+    idx = select(integrator.t, ip.times)
 
     @views z = reshape(integrator.u[1:(N*M)], N, M)
     @views B = reshape(integrator.u[(N*M+1):end], M, M)
     _calculate_expectation!(ip, z, B, idx)
 
     if integrator.p.opt.progress
-        print("\rProgress: $(round(Int, 100*idx/length(ip.tlist)))%")
+        print("\rProgress: $(round(Int, 100*idx/length(ip.times)))%")
         flush(stdout)
     end
     return u_modified!(integrator, false)
@@ -433,7 +433,7 @@ function lr_mesolveProblem(
         e_ops = e_ops,
         f_ops = f_ops,
         opt = opt,
-        tlist = t_l,
+        times = t_l,
         expvals = expvals,
         funvals = funvals,
         Ml = Ml,
@@ -540,7 +540,7 @@ function lr_mesolve(
 end
 
 function lr_mesolve(prob::ODEProblem; kwargs...)
-    sol = solve(prob, prob.p.opt.alg, tstops = prob.p.tlist)
+    sol = solve(prob, prob.p.opt.alg, tstops = prob.p.times)
     prob.p.opt.progress && print("\n")
 
     N = prob.p.N
@@ -552,8 +552,8 @@ function lr_mesolve(prob::ODEProblem; kwargs...)
     ρt = map(x -> Qobj(x[1] * x[2] * x[1]', type = Operator(), dims = prob.p.Hdims), zip(zt, Bt))
 
     return TimeEvolutionLRSol(
+        prob.p.times,
         sol.t,
-        prob.p.tlist,
         ρt,
         prob.p.expvals,
         prob.p.funvals,
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index db01499c6..fe8da0876 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -254,7 +254,7 @@ function mcsolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
-        prob_mc.tlist,
+        prob_mc.times,
         prob_mc.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
@@ -411,8 +411,8 @@ function mcsolve(
 
     return TimeEvolutionMCSol(
         ntraj,
+        ens_prob_mc.times,
         _sol_1.t,
-        ens_prob_mc.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index e61660c8c..85d6ed613 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -205,8 +205,8 @@ function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
     end
 
     return TimeEvolutionSol(
+        prob.times,
         sol.t,
-        prob.tlist,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 5a9bebf3a..ef9d98f9b 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -155,8 +155,8 @@ function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
     ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.dimensions), sol.u)
 
     return TimeEvolutionSol(
+        prob.times,
         sol.t,
-        prob.tlist,
         ψt,
         _get_expvals(sol, SaveFuncSESolve),
         sol.retcode,
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 1a782b981..aa6dbae9b 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -269,7 +269,7 @@ function smesolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
-        prob_sme.tlist,
+        prob_sme.times,
         prob_sme.dimensions,
         merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3])),
     )
@@ -424,8 +424,8 @@ function smesolve(
 
     return TimeEvolutionStochasticSol(
         ntraj,
+        ens_prob.times,
         _sol_1.t,
-        ens_prob.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index a3d394f34..e68987c02 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -262,7 +262,7 @@ function ssesolveEnsembleProblem(
 
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
-        prob_sme.tlist,
+        prob_sme.times,
         prob_sme.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
@@ -418,8 +418,8 @@ function ssesolve(
 
     return TimeEvolutionStochasticSol(
         ntraj,
+        ens_prob.times,
         _sol_1.t,
-        ens_prob.tlist,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 4187c6f95..bfdd6942b 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -14,7 +14,7 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 # Fields (Attributes)
 
 - `prob::AbstractSciMLProblem`: The `ODEProblem` of the time evolution.
-- `tlist::AbstractVector`: The time list of the evolution.
+- `times::AbstractVector`: The time list of the evolution.
 - `dimensions::AbstractDimensions`: The dimensions of the Hilbert space.
 - `kwargs::KWT`: Generic keyword arguments.
 
@@ -23,7 +23,7 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 """
 struct TimeEvolutionProblem{PT<:AbstractSciMLProblem,TT<:AbstractVector,DT<:AbstractDimensions,KWT}
     prob::PT
-    tlist::TT
+    times::TT
     dimensions::DT
     kwargs::KWT
 end
@@ -37,7 +37,7 @@ function Base.getproperty(prob::TimeEvolutionProblem, key::Symbol)
     end
 end
 
-TimeEvolutionProblem(prob, tlist, dims) = TimeEvolutionProblem(prob, tlist, dims, nothing)
+TimeEvolutionProblem(prob, times, dims) = TimeEvolutionProblem(prob, times, dims, nothing)
 
 @doc raw"""
     struct TimeEvolutionSol
@@ -46,10 +46,10 @@ A structure storing the results and some information from solving time evolution
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
-- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times`.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `tlist`.
+- `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{QuantumObject}`: The list of result states corresponding to each time point in `times_states`.
+- `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `times`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -66,7 +66,7 @@ struct TimeEvolutionSol{
     RT<:Real,
 }
     times::TT1
-    tlist::TT2
+    times_states::TT2
     states::TS
     expect::TE
     retcode::RETT
@@ -99,12 +99,12 @@ A structure storing the results and some information from solving quantum trajec
 # Fields (Attributes)
 
 - `ntraj::Int`: Number of trajectories
-- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
-- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
-- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `tlist`
+- `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`
 - `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.
@@ -126,7 +126,7 @@ struct TimeEvolutionMCSol{
 }
     ntraj::Int
     times::TT1
-    tlist::TT2
+    times_states::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
@@ -164,12 +164,12 @@ A structure storing the results and some information from solving trajectories o
 # Fields (Attributes)
 
 - `ntraj::Int`: Number of trajectories
-- `times::AbstractVector`: The list of time points at which the states are stored during the evolution.
-- `tlist::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
-- `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
-- `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `tlist`.
-- `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `tlist`
+- `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`
 - `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.
@@ -188,7 +188,7 @@ struct TimeEvolutionStochasticSol{
 }
     ntraj::Int
     times::TT1
-    tlist::TT2
+    times_states::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 1b1ef95db..98960764e 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -243,8 +243,8 @@ function dfd_mesolve(
     end
 
     return TimeEvolutionSol(
+        dfd_prob.times,
         sol.t,
-        dfd_prob.tlist,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 5e0d386b8..e17db89c2 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -84,16 +84,16 @@ end
 
     @test prob.prob.f.f isa MatrixOperator
     @test sum(abs.(sol.expect[1, :] .- amp_rabi .* sin.(Ω_rabi * tlist) .^ 2)) / length(tlist) < 0.1
-    @test length(sol.times) == 1
-    @test length(sol.tlist) == length(tlist)
+    @test length(sol.times) == length(tlist)
+    @test length(sol.times_states) == 1
     @test length(sol.states) == 1
     @test size(sol.expect) == (length(e_ops), length(tlist))
     @test length(sol2.times) == length(tlist)
-    @test length(sol2.tlist) == length(tlist)
+    @test length(sol2.times_states) == length(tlist)
     @test length(sol2.states) == length(tlist)
     @test sol2.expect === nothing
-    @test length(sol3.times) == length(saveat)
-    @test length(sol3.tlist) == length(tlist)
+    @test length(sol3.times) == length(tlist)
+    @test length(sol3.times_states) == length(saveat)
     @test length(sol3.states) == length(saveat)
     @test size(sol3.expect) == (length(e_ops), length(tlist))
     @test sol.expect[1, saveat_idxs] ≈ expect(e_ops[1], sol3.states) atol = 1e-6
@@ -170,16 +170,16 @@ end
 
     @test TESetup.prob_me.prob.f.f isa MatrixOperator
     @test isket(sol_me5.states[1])
-    @test length(sol_me.times) == 1
-    @test length(sol_me.tlist) == length(tlist)
+    @test length(sol_me.times) == length(tlist)
+    @test length(sol_me.times_states) == 1
     @test length(sol_me.states) == 1
     @test size(sol_me.expect) == (length(e_ops), length(tlist))
     @test length(sol_me2.times) == length(tlist)
-    @test length(sol_me2.tlist) == length(tlist)
+    @test length(sol_me2.times_states) == length(tlist)
     @test length(sol_me2.states) == length(tlist)
     @test sol_me2.expect === nothing
-    @test length(sol_me3.times) == length(saveat)
-    @test length(sol_me3.tlist) == length(tlist)
+    @test length(sol_me3.times) == length(tlist)
+    @test length(sol_me3.times_states) == length(saveat)
     @test length(sol_me3.states) == length(saveat)
     @test size(sol_me3.expect) == (length(e_ops), length(tlist))
     @test sol_me3.expect[1, TESetup.saveat_idxs] ≈ expect(e_ops[1], sol_me3.states) atol = 1e-6
@@ -294,11 +294,11 @@ end
     @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 length(sol_mc.times) == 1
-    @test length(sol_mc.tlist) == length(tlist)
+    @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 length(sol_mc_states.times) == length(saveat)
-    @test length(sol_mc_states.tlist) == length(tlist)
+    @test length(sol_mc_states.times) == length(tlist)
+    @test length(sol_mc_states.times_states) == length(saveat)
     @test sol_mc_states.expect === nothing
 
     sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
@@ -398,8 +398,8 @@ end
     )
 
     @test sum(abs, sol_sse.expect .- sol_me.expect) / length(tlist) < 0.1
-    @test length(sol_sse.times) == 1
-    @test length(sol_sse.tlist) == length(tlist)
+    @test length(sol_sse.times) == length(tlist)
+    @test length(sol_sse.times_states) == 1
     @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)
@@ -533,8 +533,8 @@ end
 
     @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 length(sol_sme.times) == 1
-    @test length(sol_sme.tlist) == length(tlist)
+    @test length(sol_sme.times) == length(tlist)
+    @test length(sol_sme.times_states) == 1
     @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)