qutip · ytdHuang · Oct 5, 2025 · Oct 5, 2025
diff --git a/docs/src/LS_solvers.md b/docs/src/LS_solvers.md
@@ -1,39 +1,39 @@
 # [LinearSolve solvers](@id LS-solvers)
 
-In this page, we list several recommended solvers provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving [`steadystate`](@ref) and spectrum in hierarchical equations of motion approach.  
+In this page, we list several recommended `alg`orithms provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving [`steadystate`](@ref) and spectrum in hierarchical equations of motion approach.  
 
 Remember to import `LinearSolve.jl`
 
 ```julia
 using LinearSolve
 ```
 
-(click [here](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/) to see the full solver list provided by `LinearSolve.jl`)
+(click [here](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/) to see the full `alg`orithm list provided by `LinearSolve.jl`)
 
-### A generic BICGSTAB implementation from Krylov (Default solver)
+### A generic GMRES implementation from Krylov (Default algorithm)
 
 ```julia
 KrylovJL_GMRES(rtol=1e-12, atol=1e-14)
 ```
 
 ### UMFPACKFactorization
-This solver performs better when there is more structure to the sparsity pattern (depends on the complexity of your system and baths).
+This `alg`orithm performs better when there is more structure to the sparsity pattern (depends on the complexity of your system and baths).
 
 ```julia
 UMFPACKFactorization()
 ```
 
 ### KLUFactorization
-This solver performs better when there is less structure to the sparsity pattern (depends on the complexity of your system and baths).
+This `alg`orithm performs better when there is less structure to the sparsity pattern (depends on the complexity of your system and baths).
 
 ```julia
 KLUFactorization()
 ```
 
 ### Pardiso
-This solver is based on Intel openAPI Math Kernel Library (MKL) Pardiso
+This `alg`orithm is based on Intel openAPI Math Kernel Library (MKL) Pardiso
 !!! note "Note"
-    Using this solver requires adding the package [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl), i.e. `using Pardiso`
+    Using this `alg`orithm requires adding the package [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl), i.e. `using Pardiso`
 
 ```julia
 using Pardiso

diff --git a/docs/src/ODE_solvers.md b/docs/src/ODE_solvers.md
@@ -1,19 +1,19 @@
 # [DifferentialEquations solvers](@id ODE-solvers)
 
-In this page, we list several recommended solvers provided by [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) for solving time evolution in hierarchical equations of motion approach.  
+In this page, we list several recommended `alg`orithms provided by [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) for solving time evolution in hierarchical equations of motion approach.  
 
 Remember to import `OrdinaryDiffEq.jl` (or `DifferentialEquations.jl`)
 
 ```julia
 using OrdinaryDiffEq ## or "using DifferentialEquations" 
 ```
 
-(click [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) to see the full solver list provided by `DifferentialEquations.jl`)
+(click [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) to see the full `alg`orithm list provided by `DifferentialEquations.jl`)
 
 For any extra solver options, we can add it in the function `HEOMsolve` with keyword arguments. These keyword arguments will be directly pass to the solvers in `DifferentialEquations`
 (click [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) to see the documentation for the common solver options)
 
-### DP5 (Default solver)
+### DP5 (Default algorithm)
 Dormand-Prince's 5/4 Runge-Kutta method. (free 4th order interpolant)
 
 ```julia

diff --git a/docs/src/install.md b/docs/src/install.md
@@ -30,10 +30,10 @@ HierarchicalEOM.about()
 In order to get a better experience and take full advantage of `HierarchicalEOM`, we recommend to install the following external packages:
 
 #### [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/)
-`DifferentialEquations` is needed to provide the low-level ODE solvers especially for solving [time evolution](@ref doc-Time-Evolution). For [low dependency usage](https://diffeq.sciml.ai/stable/features/low_dep/), users can use [`OrdinaryDiffEq.jl`](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) instead.
+`DifferentialEquations` is needed to provide the low-level ODE `alg`orithms especially for solving [time evolution](@ref doc-Time-Evolution). For [low dependency usage](https://diffeq.sciml.ai/stable/features/low_dep/), users can use [`OrdinaryDiffEq.jl`](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) instead.
 
 #### [LinearSolve.jl](http://linearsolve.sciml.ai/stable/)
-`LinearSolve` is a unified interface for the linear solving packages of Julia. It interfaces with other packages of the Julia ecosystem to make it easier to test alternative solver packages and pass small types to control algorithm swapping. It is needed to provide the solvers especially for solving [stationary state](@ref doc-Stationary-State) and [spectra](@ref doc-Spectrum) for both bosonic and fermionic systems.
+`LinearSolve` is a unified interface for the linear solving packages of Julia. It interfaces with other packages of the Julia ecosystem to make it easier to test alternative solver packages and pass small types to control algorithm swapping. It is needed to provide the `alg`orithms especially for solving [stationary state](@ref doc-Stationary-State) and [spectra](@ref doc-Spectrum) for both bosonic and fermionic systems.
 
 #### [JLD2.jl](https://juliaio.github.io/JLD2.jl/stable/)
 `JLD2` saves and loads Julia data structures in a format comprising a subset of HDF5. Because the size of matrix in `HierarchicalEOM` is usually super large and leads to long time calculation, we support the functionality for saving and loading the `HierarchicalEOM`-type objects into files by `JLD2 >= 0.4.23`.

diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md
@@ -8,7 +8,7 @@ We briefly summarize how to numerically compute the spectrum associated with the
  - [Power Spectrum](@ref doc-PS)
  - [Density of States](@ref doc-DOS)
 
-`HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [LinearSolve solvers](@ref LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
+`HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many `alg`orithms. It offers quite a few options for the user to tailor the solver to their specific needs. The default `alg`orithm (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [LinearSolve solvers](@ref LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
 
 !!! compat "Extension for CUDA.jl"
     `HierarchicalEOM.jl` provides an extension to support GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) acceleration for solving the spectrum, but this feature requires `Julia 1.9+` and `HierarchicalEOM 1.1+`. See [here](@ref doc-ext-CUDA) for more details.

diff --git a/docs/src/stationary_state.md b/docs/src/stationary_state.md
@@ -12,7 +12,7 @@ The first method is implemented by solving the linear problem
 0=\hat{\mathcal{M}}\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)
 ```
 
-`HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [LinearSolve solvers](@ref LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
+`HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many `alg`orithms. It offers quite a few options for the user to tailor the solver to their specific needs. The default `alg`orithm (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [LinearSolve solvers](@ref LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
 
 ```julia
 # the HEOMLS matrix
@@ -30,7 +30,7 @@ The second method is implemented by solving the ordinary differential equation (
 ```
 until finding a stationary solution.
 
-`HierarchicalEOM.jl` wraps some of the functions in [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/), which is a very rich numerical library for solving the differential equations and provides many ODE solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to the documentation of [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/).
+`HierarchicalEOM.jl` wraps some of the functions in [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/), which is a very rich numerical library for solving the differential equations and provides many ODE `alg`orithms. It offers quite a few options for the user to tailor the solver to their specific needs. The default `alg`orithm (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to the documentation of [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/).
 
 ### Given the initial state as Density Operator (`QuantumObject` type)
 

diff --git a/docs/src/time_evolution.md b/docs/src/time_evolution.md
@@ -1,6 +1,6 @@
 # [Time Evolution](@id doc-Time-Evolution)
 ## Introduction
-`HierarchicalEOM.jl` implements various methods and solvers to simulate the open quantum system dynamics. 
+`HierarchicalEOM.jl` implements various methods and `alg`orithms to simulate the open quantum system dynamics. 
 The [HEOM Liouvillian superoperator (HEOMLS) matrix](@ref doc-HEOMLS-Matrix) ``\hat{\mathcal{M}}`` characterizes the dynamics of the reduce state and in the full extended space of all [auxiliary density operators (ADOs)](@ref doc-ADOs) ``\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)``, namely
 ```math
 \begin{equation}
@@ -76,15 +76,15 @@ end
 ```
 
 ## Ordinary Differential Equation (ODE) Method
-The first method is implemented by solving the ordinary differential equation (ODE). `HierarchicalEOM.jl` wraps some of the functions in [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/), which is a very rich numerical library for solving the differential equations and provides many ODE solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [DifferentialEquations solvers](@ref ODE-solvers) and also the documentation of [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/).
+The first method is implemented by solving the ordinary differential equation (ODE). `HierarchicalEOM.jl` wraps some of the functions in [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/), which is a very rich numerical library for solving the differential equations and provides many ODE `alg`orithms. It offers quite a few options for the user to tailor the solver to their specific needs. The default `alg`orithm (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [DifferentialEquations solvers](@ref ODE-solvers) and also the documentation of [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/).
 
 !!! compat "Extension for CUDA.jl"
     `HierarchicalEOM.jl` provides an extension to support GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) acceleration for [`HEOMsolve`](@ref) (only for ODE method). See [here](@ref doc-ext-CUDA) for more details.
 
 See the docstring of this method:  
 
 ```@docs
-HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, solver::OrdinaryDiffEqAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), verbose::Bool = true, filename::String = "", SOLVEROptions...,) where {T_state<:Union{QuantumObject,ADOs}}
+HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, alg::OrdinaryDiffEqAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), verbose::Bool = true, filename::String = "", kwargs...,) where {T_state<:Union{QuantumObject,ADOs}}
 ```
 
 ```julia

diff --git a/src/density_of_states.jl b/src/density_of_states.jl
@@ -1,7 +1,7 @@
 export DensityOfStates
 
 @doc raw"""
-    DensityOfStates(M, ρ, d_op, ωlist; solver, verbose, filename, SOLVEROptions...)
+    DensityOfStates(M, ρ, d_op, ωlist; alg, verbose, filename, kwargs...)
 Calculate density of states for the fermionic system in frequency domain.
 
 ```math
@@ -13,13 +13,13 @@ Calculate density of states for the fermionic system in frequency domain.
 - `ρ::Union{QuantumObject,ADOs}` :  the system density matrix or the auxiliary density operators.
 - `d_op::QuantumObject` : The annihilation operator (``d`` as shown above) acting on the fermionic system.
 - `ωlist::AbstractVector` : the specific frequency points to solve.
-- `solver::SciMLLinearSolveAlgorithm` : solver in package `LinearSolve.jl`. Default to `KrylovJL_GMRES(rtol=1e-12, atol=1e-14)`.
+- `alg::SciMLLinearSolveAlgorithm` : The solving algorithm in package `LinearSolve.jl`. Default to `KrylovJL_GMRES(rtol=1e-12, atol=1e-14)`.
 - `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`.
 - `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process.
-- `SOLVEROptions` : extra options for solver 
+- `kwargs` : The keyword arguments for `LinearProblem`.
 
 # Notes
-- For more details about `solver` and `SOLVEROptions`, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/)
+- For more details about `alg`, `kwargs`, and `LinearProblem`, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/)
 
 # Returns
 - `dos::AbstractVector` : the list of density of states corresponds to the specified `ωlist`
@@ -29,11 +29,13 @@ Calculate density of states for the fermionic system in frequency domain.
     ρ::Union{QuantumObject,ADOs},
     d_op::QuantumObject,
     ωlist::AbstractVector;
-    solver::SciMLLinearSolveAlgorithm = KrylovJL_GMRES(rtol = 1e-12, atol = 1e-14),
+    alg::SciMLLinearSolveAlgorithm = KrylovJL_GMRES(rtol = 1e-12, atol = 1e-14),
     verbose::Bool = true,
     filename::String = "",
-    SOLVEROptions...,
+    kwargs...,
 )
+    haskey(kwargs, :solver) &&
+        error("The keyword argument `solver` for DensityOfStates has been deprecated, please use `alg` instead.")
 
     # check M
     if M.parity == EVEN
@@ -85,10 +87,10 @@ Calculate density of states for the fermionic system in frequency domain.
         Iω = i * ω * I_total
 
         if prog.counter[] == 0
-            cache_m = init(LinearProblem(M.data.A - Iω, b_m), solver, SOLVEROptions...)
+            cache_m = init(LinearProblem(M.data.A - Iω, b_m), alg, kwargs...)
             sol_m = solve!(cache_m)
 
-            cache_p = init(LinearProblem(M.data.A + Iω, b_p), solver, SOLVEROptions...)
+            cache_p = init(LinearProblem(M.data.A + Iω, b_p), alg, kwargs...)
             sol_p = solve!(cache_p)
         else
             cache_m.A = M.data.A - Iω