diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8e189994c..aec3b9c0c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Introduce new methods of `sesolve_map` and `mesolve_map` for advanced usage. Users can now customize their own `iter`ator structure, `prob_func` and `output_func`. ([#565])
+- Use `ProgressMeter.jl` for progress bar rather than our in-house implementation. ([#569])
 
 ## [v0.37.0]
 Release date: 2025-10-12
@@ -339,3 +340,4 @@ Release date: 2024-11-13
 [#555]: https://github.com/qutip/QuantumToolbox.jl/issues/555
 [#557]: https://github.com/qutip/QuantumToolbox.jl/issues/557
 [#565]: https://github.com/qutip/QuantumToolbox.jl/issues/565
+[#569]: https://github.com/qutip/QuantumToolbox.jl/issues/569
diff --git a/Project.toml b/Project.toml
index b36265928..5331fbdea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumToolbox"
 uuid = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
-authors = ["Alberto Mercurio", "Yi-Te Huang"]
 version = "0.37.0"
+authors = ["Alberto Mercurio", "Yi-Te Huang"]
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -18,6 +18,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
@@ -60,6 +61,7 @@ Makie = "0.24"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqTsit5 = "1"
 Pkg = "1"
+ProgressMeter = "1.11.0"
 Random = "1"
 SciMLBase = "2.105"
 SciMLOperators = "1.4"
diff --git a/docs/src/index.md b/docs/src/index.md
index a8f3baaa8..460e26f1a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -104,6 +104,8 @@ In order to get a better experience and take full advantage of `QuantumToolbox`
   - [`Zygote.jl`](https://github.com/FluxML/Zygote.jl)
   - [`Enzyme.jl`](https://github.com/EnzymeAD/Enzyme.jl)
   - [`ForwardDiff.jl`](https://github.com/JuliaDiff/ForwardDiff.jl)
+- Progress Bars:
+  - [`ProgressMeter.jl`](https://github.com/timholy/ProgressMeter.jl)
 - Packages for other advanced usage:
   - [`StaticArrays.jl`](https://github.com/JuliaArrays/StaticArrays.jl)
   - [`SciMLOperators.jl`](https://github.com/SciML/SciMLOperators.jl)
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index c3d9be754..d43060182 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -62,6 +62,7 @@ import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
 import LaTeXStrings: @L_str
 import Pkg
+import ProgressMeter: Progress, next!
 import Random: AbstractRNG, default_rng, seed!
 import SpecialFunctions: loggamma
 import StaticArraysCore: SVector, MVector
@@ -81,7 +82,6 @@ export cache_operator, iscached, isconstant
 include("settings.jl")
 include("utilities.jl")
 include("versioninfo.jl")
-include("progress_bar.jl")
 include("linear_maps.jl")
 
 # Quantum Object
diff --git a/src/deprecated.jl b/src/deprecated.jl
index b1378f311..74717bc79 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -130,3 +130,11 @@ function MultiSiteOperator(args...)
     )
     return multisite_operator(args...)
 end
+
+export ProgressBar
+function ProgressBar(args...; kwargs...)
+    # Use error instead of depwarn, since ProgressBar and Progress have different arguments
+    return error(
+        "`ProgressBar` is deprecated and will be removed in next major release. Use `Progress` from `ProgressMeter.jl` instead.",
+    )
+end
diff --git a/src/progress_bar.jl b/src/progress_bar.jl
deleted file mode 100644
index 0fe989f52..000000000
--- a/src/progress_bar.jl
+++ /dev/null
@@ -1,85 +0,0 @@
-export ProgressBar, next!
-
-mutable struct ProgressBar{CT,T1<:Integer,T2<:Real,T3,T4<:Real,LT}
-    counter::CT
-    max_counts::T1
-    enable::Bool
-    bar_width::T1
-    start_time::T2
-    previous_time::T3
-    interval::T4
-    lock::LT
-end
-
-function ProgressBar(max_counts::Int; enable::Bool = true, bar_width::Int = 30, interval = 1.0)
-    start_time = time()
-    lock = ReentrantLock()
-    return ProgressBar(
-        Threads.Atomic{Int}(0),
-        max_counts,
-        enable,
-        bar_width,
-        start_time,
-        Threads.Atomic{Float64}(start_time),
-        interval,
-        lock,
-    )
-end
-
-function next!(p::ProgressBar, io::IO = stdout)
-    p.counter[] >= p.max_counts && return
-
-    Threads.atomic_add!(p.counter, 1)
-
-    !p.enable && return
-
-    counter = p.counter[]
-    max_counts = p.max_counts
-    bar_width = p.bar_width
-    start_time = p.start_time
-    previous_time = p.previous_time[]
-    interval = p.interval
-
-    current_time = time()
-
-    ((current_time - previous_time) < interval && counter != p.max_counts) && return
-
-    Threads.atomic_xchg!(p.previous_time, current_time)
-
-    percentage = counter / max_counts
-    percentage_100 = lpad(round(100 * percentage, digits = 1), 5, " ")
-    progress = floor(Int, bar_width * percentage)
-
-    # Calculate the elapsed time in seconds
-    elapsed_time = floor(Int, current_time - start_time)
-    # Convert the elapsed time into a string in hours, minutes and seconds
-    elapsed_time_str = string(
-        elapsed_time ÷ 3600,
-        "h ",
-        lpad((elapsed_time % 3600) ÷ 60, 2, "0"),
-        "m ",
-        lpad(elapsed_time % 60, 2, "0"),
-        "s",
-    )
-
-    # Calculate the estimated time of arrival
-    eta = floor(Int, elapsed_time / counter * (max_counts - counter))
-    # convert eta into a string in hours, minutes and seconds
-    eta_str = string(eta ÷ 3600, "h ", lpad((eta % 3600) ÷ 60, 2, "0"), "m ", lpad(eta % 60, 2, "0"), "s")
-
-    # Construct the progress bar string
-    bar = "[" * repeat("=", progress) * repeat(" ", bar_width - progress) * "]"
-
-    lock(p.lock)
-    try
-        print(io, "\rProgress: $bar $percentage_100% --- Elapsed Time: $elapsed_time_str (ETA: $eta_str)")
-
-        counter == p.max_counts && print(io, "\n")
-
-        flush(io)
-    finally
-        unlock(p.lock)
-    end
-
-    return nothing
-end
diff --git a/src/time_evolution/callback_helpers/callback_helpers.jl b/src/time_evolution/callback_helpers/callback_helpers.jl
index 7dbdc135d..e45de70b4 100644
--- a/src/time_evolution/callback_helpers/callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/callback_helpers.jl
@@ -57,7 +57,10 @@ end
 function _generate_save_callback(e_ops, tlist, progress_bar, method)
     e_ops_data = e_ops isa Nothing ? nothing : _get_e_ops_data(e_ops, method)
 
-    progr = getVal(progress_bar) ? ProgressBar(length(tlist), enable = getVal(progress_bar)) : nothing
+    progr_desc = _get_progress_desc(method)
+    progr =
+        getVal(progress_bar) ?
+        Progress(length(tlist), showspeed = true, enabled = getVal(progress_bar), desc = progr_desc) : nothing
 
     expvals = e_ops isa Nothing ? nothing : Array{ComplexF64}(undef, length(e_ops), length(tlist))
 
@@ -69,7 +72,7 @@ function _generate_stochastic_save_callback(e_ops, sc_ops, tlist, store_measurem
     e_ops_data = e_ops isa Nothing ? nothing : _get_e_ops_data(e_ops, method)
     m_ops_data = _get_m_ops_data(sc_ops, method)
 
-    progr = getVal(progress_bar) ? ProgressBar(length(tlist), enable = getVal(progress_bar)) : nothing
+    progr = getVal(progress_bar) ? Progress(length(tlist), showspeed = true, enabled = getVal(progress_bar)) : nothing
 
     expvals = e_ops isa Nothing ? nothing : Array{ComplexF64}(undef, length(e_ops), length(tlist))
     m_expvals = getVal(store_measurement) ? Array{Float64}(undef, length(sc_ops), length(tlist) - 1) : nothing
diff --git a/src/time_evolution/callback_helpers/mcsolve_callback_helpers.jl b/src/time_evolution/callback_helpers/mcsolve_callback_helpers.jl
index d96b50202..cafbb5550 100644
--- a/src/time_evolution/callback_helpers/mcsolve_callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/mcsolve_callback_helpers.jl
@@ -106,7 +106,7 @@ function _generate_mcsolve_kwargs(ψ0, T, e_ops, tlist, c_ops, jump_callback, rn
     end
 
     if e_ops isa Nothing
-        # We are implicitly saying that we don't have a `ProgressBar`
+        # We are implicitly saying that we don't have a `Progress`
         kwargs2 =
             haskey(kwargs, :callback) ? merge(kwargs, (callback = CallbackSet(cb1, kwargs.callback),)) :
             merge(kwargs, (callback = cb1,))
diff --git a/src/time_evolution/callback_helpers/mesolve_callback_helpers.jl b/src/time_evolution/callback_helpers/mesolve_callback_helpers.jl
index a24229d29..ba98aa6b2 100644
--- a/src/time_evolution/callback_helpers/mesolve_callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/mesolve_callback_helpers.jl
@@ -2,7 +2,7 @@
 Helper functions for the mesolve callbacks.
 =#
 
-struct SaveFuncMESolve{TE,PT<:Union{Nothing,ProgressBar},IT,TEXPV<:Union{Nothing,AbstractMatrix}} <: AbstractSaveFunc
+struct SaveFuncMESolve{TE,PT<:Union{Nothing,Progress},IT,TEXPV<:Union{Nothing,AbstractMatrix}} <: AbstractSaveFunc
     e_ops::TE
     progr::PT
     iter::IT
@@ -14,6 +14,8 @@ end
 
 _get_e_ops_data(e_ops, ::Type{SaveFuncMESolve}) = [_generate_mesolve_e_op(op) for op in e_ops] # Broadcasting generates type instabilities on Julia v1.10
 
+_get_progress_desc(::Type{SaveFuncMESolve}) = "(mesolve) "
+
 ##
 
 # When e_ops is a list of operators
diff --git a/src/time_evolution/callback_helpers/sesolve_callback_helpers.jl b/src/time_evolution/callback_helpers/sesolve_callback_helpers.jl
index 2bbff8bf0..bec73cc30 100644
--- a/src/time_evolution/callback_helpers/sesolve_callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/sesolve_callback_helpers.jl
@@ -2,7 +2,7 @@
 Helper functions for the sesolve callbacks.
 =#
 
-struct SaveFuncSESolve{TE,PT<:Union{Nothing,ProgressBar},IT,TEXPV<:Union{Nothing,AbstractMatrix}} <: AbstractSaveFunc
+struct SaveFuncSESolve{TE,PT<:Union{Nothing,Progress},IT,TEXPV<:Union{Nothing,AbstractMatrix}} <: AbstractSaveFunc
     e_ops::TE
     progr::PT
     iter::IT
@@ -14,6 +14,8 @@ end
 
 _get_e_ops_data(e_ops, ::Type{SaveFuncSESolve}) = get_data.(e_ops)
 
+_get_progress_desc(::Type{SaveFuncSESolve}) = "(sesolve) "
+
 ##
 
 # When e_ops is a list of operators
diff --git a/src/time_evolution/callback_helpers/smesolve_callback_helpers.jl b/src/time_evolution/callback_helpers/smesolve_callback_helpers.jl
index 6fb75a979..677dd8589 100644
--- a/src/time_evolution/callback_helpers/smesolve_callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/smesolve_callback_helpers.jl
@@ -6,7 +6,7 @@ struct SaveFuncSMESolve{
     SM,
     TE,
     TE2,
-    PT<:Union{Nothing,ProgressBar},
+    PT<:Union{Nothing,Progress},
     IT,
     TEXPV<:Union{Nothing,AbstractMatrix},
     TMEXPV<:Union{Nothing,AbstractMatrix},
diff --git a/src/time_evolution/callback_helpers/ssesolve_callback_helpers.jl b/src/time_evolution/callback_helpers/ssesolve_callback_helpers.jl
index e28ccc6ef..437558533 100644
--- a/src/time_evolution/callback_helpers/ssesolve_callback_helpers.jl
+++ b/src/time_evolution/callback_helpers/ssesolve_callback_helpers.jl
@@ -6,7 +6,7 @@ struct SaveFuncSSESolve{
     SM,
     TE,
     TE2,
-    PT<:Union{Nothing,ProgressBar},
+    PT<:Union{Nothing,Progress},
     IT,
     TEXPV<:Union{Nothing,AbstractMatrix},
     TMEXPV<:Union{Nothing,AbstractMatrix},
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index ef195b591..d21055fd6 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -205,7 +205,7 @@ If the environmental measurements register a quantum jump, the wave function und
 - `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `prob_func`: Function to use for generating the ODEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `kwargs`: The keyword arguments for the ODEProblem.
 
 # Notes
@@ -238,7 +238,13 @@ function mcsolveEnsembleProblem(
     _prob_func = isnothing(prob_func) ? _ensemble_dispatch_prob_func(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
+        _ensemble_dispatch_output_func(
+            ensemblealg,
+            progress_bar,
+            ntraj,
+            _mcsolve_output_func;
+            progr_desc = "(mcsolve) ",
+        ) : output_func
 
     prob_mc = mcsolveProblem(
         H,
@@ -332,7 +338,7 @@ If the environmental measurements register a quantum jump, the wave function und
 - `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `prob_func`: Function to use for generating the ODEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `normalize_states`: Whether to normalize the states. Default to `Val(true)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 0c0b659a3..e010533dc 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -351,7 +351,7 @@ mesolve_map(
 # this method is for advanced usage
 # User can define their own iterator structure, prob_func and output_func
 #   - `prob_func`: Function to use for generating the ODEProblem.
-#   - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+#   - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 #
 # Return: An array of TimeEvolutionSol objects with the size same as the given iter.
 function mesolve_map(
@@ -368,7 +368,13 @@ function mesolve_map(
     _prob_func = isnothing(prob_func) ? (prob, i, repeat) -> _se_me_map_prob_func(prob, i, repeat, iter) : prob_func
     _output_func =
         isnothing(output_func) ?
-        _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _standard_output_func) : output_func
+        _ensemble_dispatch_output_func(
+            ensemblealg,
+            progress_bar,
+            ntraj,
+            _standard_output_func;
+            progr_desc = "(mesolve_map) ",
+        ) : output_func
     ens_prob = TimeEvolutionProblem(
         EnsembleProblem(prob.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob.times,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 634b48f50..365d9f331 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -260,7 +260,7 @@ sesolve_map(H::Union{AbstractQuantumObject{Operator},Tuple}, ψ0::QuantumObject{
 # this method is for advanced usage
 # User can define their own iterator structure, prob_func and output_func
 #   - `prob_func`: Function to use for generating the ODEProblem.
-#   - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+#   - `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 #
 # Return: An array of TimeEvolutionSol objects with the size same as the given iter.
 function sesolve_map(
@@ -277,7 +277,13 @@ function sesolve_map(
     _prob_func = isnothing(prob_func) ? (prob, i, repeat) -> _se_me_map_prob_func(prob, i, repeat, iter) : prob_func
     _output_func =
         isnothing(output_func) ?
-        _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _standard_output_func) : output_func
+        _ensemble_dispatch_output_func(
+            ensemblealg,
+            progress_bar,
+            ntraj,
+            _standard_output_func;
+            progr_desc = "(sesolve_map) ",
+        ) : output_func
     ens_prob = TimeEvolutionProblem(
         EnsembleProblem(prob.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob.times,
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 52fedc93b..6779bd70e 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -107,8 +107,6 @@ function smesolveProblem(
         to_dense(_complex_float_type(T), mat2vec(ket2dm(ψ0).data))
     end
 
-    progr = ProgressBar(length(tlist), enable = getVal(progress_bar))
-
     sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list))
 
     K = get_data(L_evo)
@@ -203,7 +201,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 - `ntraj`: Number of trajectories to use. Default is `500`.
 - `ensemblealg`: Ensemble method to use. Default to `EnsembleThreads()`.
 - `prob_func`: Function to use for generating the SDEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `store_measurement`: Whether to store the measurement expectation values. Default is `Val(false)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
@@ -251,7 +249,13 @@ function smesolveEnsembleProblem(
         ) : prob_func
     _output_func =
         output_func isa Nothing ?
-        _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _standard_output_func) : output_func
+        _ensemble_dispatch_output_func(
+            ensemblealg,
+            progress_bar,
+            ntraj,
+            _standard_output_func;
+            progr_desc = "(smesolve) ",
+        ) : output_func
 
     prob_sme = smesolveProblem(
         H,
@@ -332,7 +336,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 - `ntraj`: Number of trajectories to use. Default is `500`.
 - `ensemblealg`: Ensemble method to use. Default to `EnsembleThreads()`.
 - `prob_func`: Function to use for generating the SDEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `store_measurement`: Whether to store the measurement expectation values. Default is `Val(false)`.
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 67b15954d..a8ca0de98 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -103,8 +103,6 @@ function ssesolveProblem(
 
     ψ0 = to_dense(_complex_float_type(ψ0), get_data(ψ0))
 
-    progr = ProgressBar(length(tlist), enable = getVal(progress_bar))
-
     sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list))
 
     # Here the coefficients depend on the state, so this is a non-linear operator, which should be implemented with FunctionOperator instead. However, the nonlinearity is only on the coefficients, and it should be safe.
@@ -198,7 +196,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and  ``dW_n(t)`` is t
 - `ensemblealg`: Ensemble method to use. Default to `EnsembleThreads()`.
 - `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise.
 - `prob_func`: Function to use for generating the SDEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `store_measurement`: Whether to store the measurement results. Default is `Val(false)`.
 - `kwargs`: The keyword arguments for the ODEProblem.
@@ -245,7 +243,13 @@ function ssesolveEnsembleProblem(
         ) : prob_func
     _output_func =
         output_func isa Nothing ?
-        _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _standard_output_func) : output_func
+        _ensemble_dispatch_output_func(
+            ensemblealg,
+            progress_bar,
+            ntraj,
+            _standard_output_func;
+            progr_desc = "(ssesolve) ",
+        ) : output_func
 
     prob_sme = ssesolveProblem(
         H,
@@ -327,7 +331,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 - `ntraj`: Number of trajectories to use. Default is `500`.
 - `ensemblealg`: Ensemble method to use. Default to `EnsembleThreads()`.
 - `prob_func`: Function to use for generating the SDEProblem.
-- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `ProgressBar` object, and the (optional) `RemoteChannel` object.
+- `output_func`: a `Tuple` containing the `Function` to use for generating the output of a single trajectory, the (optional) `Progress` object, and the (optional) `RemoteChannel` object.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
 - `keep_runs_results`: Whether to save the results of each trajectory. Default to `Val(false)`.
 - `store_measurement`: Whether to store the measurement results. Default is `Val(false)`.
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 69b46d5f8..d7c332e27 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -373,10 +373,11 @@ function _ensemble_dispatch_output_func(
     ::ET,
     progress_bar,
     ntraj,
-    output_func,
+    output_func;
+    progr_desc = "Progress: ",
 ) where {ET<:Union{EnsembleSerial,EnsembleThreads}}
     if getVal(progress_bar)
-        progr = ProgressBar(ntraj, enable = getVal(progress_bar))
+        progr = Progress(ntraj, showspeed = true, enabled = getVal(progress_bar), desc = progr_desc)
         f = (sol, i) -> _ensemble_output_func_progress(sol, i, progr, output_func)
         return (f, progr, nothing)
     else
@@ -387,10 +388,11 @@ function _ensemble_dispatch_output_func(
     ::ET,
     progress_bar,
     ntraj,
-    output_func,
+    output_func;
+    progr_desc = "Progress... ",
 ) where {ET<:Union{EnsembleSplitThreads,EnsembleDistributed}}
     if getVal(progress_bar)
-        progr = ProgressBar(ntraj, enable = getVal(progress_bar))
+        progr = Progress(ntraj, showspeed = true, enabled = getVal(progress_bar), desc = progr_desc)
         progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
 
         f = (sol, i) -> _ensemble_output_func_distributed(sol, i, progr_channel, output_func)
diff --git a/test/core-test/progress_bar.jl b/test/core-test/progress_bar.jl
deleted file mode 100644
index 2792a6e0a..000000000
--- a/test/core-test/progress_bar.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-@testitem "Progress Bar" begin
-    bar_width = 30
-    strLength = 67 + bar_width # including "\r" in the beginning of the string
-    prog = ProgressBar(bar_width, enable = true, bar_width = bar_width, interval = 0.2)
-    for p in 1:bar_width
-        sleep(0.3)
-
-        output = sprint((t, s) -> next!(s, t), prog)
-
-        if p < bar_width
-            @test length(output) == strLength
-        else # the last output has an extra "\n" in the end
-            @test length(output) == strLength + 1
-        end
-    end
-
-    prog = ProgressBar(bar_width, enable = true, bar_width = bar_width, interval = 0.2)
-    # No delay, so no output printed
-    for p in 1:bar_width
-        output = sprint((t, s) -> next!(s, t), prog)
-
-        if p < bar_width
-            @test length(output) == 0
-        else
-            @test length(output) == strLength + 1
-        end
-    end
-end