diff --git a/Project.toml b/Project.toml
index 937b7894..be0a924d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
-authors = ["Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>"]
-version = "2.3.3"
+authors = ["Yi-Te Huang"]
+version = "2.4.0"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -37,7 +37,7 @@ LinearSolve = "2.4.2 - 2"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "1"
-QuantumToolbox = "0.22 - 0.24"
+QuantumToolbox = "0.25"
 Reexport = "1"
 SciMLBase = "2"
 SciMLOperators = "0.3"
diff --git a/docs/src/ADOs.md b/docs/src/ADOs.md
index ea2b37ce..1535355e 100644
--- a/docs/src/ADOs.md
+++ b/docs/src/ADOs.md
@@ -21,7 +21,7 @@ which is usually obtained after solving the [time evolution](@ref doc-Time-Evolu
 ## Fields
 The fields of the structure [`ADOs`](@ref) are as follows:
  - `data` : the vectorized auxiliary density operators
- - `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+ - `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
  - `N` : the number of auxiliary density operators
 - `parity`: the [parity](@ref doc-Parity) label
 
@@ -31,6 +31,7 @@ One obtain the value of each fields as follows:
 ados::ADOs
 
 ados.data
+ados.dimensions
 ados.dims
 ados.N
 ados.parity
diff --git a/docs/src/heom_matrix/M_Boson.md b/docs/src/heom_matrix/M_Boson.md
index 95a35047..09857df4 100644
--- a/docs/src/heom_matrix/M_Boson.md
+++ b/docs/src/heom_matrix/M_Boson.md
@@ -31,7 +31,7 @@ M_odd  = M_Boson(Hs, tier, Bath, ODD)
 The fields of the structure [`M_Boson`](@ref) are as follows:
  - `data` : the sparse matrix of HEOM Liouvillian superoperator
  - `tier` : the tier (cutoff level) for the bosonic hierarchy
- - `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+ - `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
  - `N` : the number of total [ADOs](@ref doc-ADOs)
  - `sup_dim` : the dimension of system superoperator
  - `parity` : the [parity](@ref doc-Parity) label of the operator which HEOMLS is acting on.
@@ -44,6 +44,7 @@ M::M_Boson
 
 M.data
 M.tier
+M.dimensions
 M.dims
 M.N
 M.sup_dim
diff --git a/docs/src/heom_matrix/M_Boson_Fermion.md b/docs/src/heom_matrix/M_Boson_Fermion.md
index 5f1300df..7c231a80 100644
--- a/docs/src/heom_matrix/M_Boson_Fermion.md
+++ b/docs/src/heom_matrix/M_Boson_Fermion.md
@@ -36,7 +36,7 @@ The fields of the structure [`M_Boson_Fermion`](@ref) are as follows:
  - `data` : the sparse matrix of HEOM Liouvillian superoperator
  - `Btier` : the tier (cutoff level) for bosonic hierarchy
  - `Ftier` : the tier (cutoff level) for fermionic hierarchy
- - `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+ - `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
  - `N` : the number of total [ADOs](@ref doc-ADOs)
  - `sup_dim` : the dimension of system superoperator
  - `parity` : the [parity](@ref doc-Parity) label of the operator which HEOMLS is acting on. 
@@ -51,6 +51,7 @@ M::M_Boson_Fermion
 M.data
 M.Btier
 M.Ftier
+M.dimensions
 M.dims
 M.N
 M.sup_dim
diff --git a/docs/src/heom_matrix/M_Fermion.md b/docs/src/heom_matrix/M_Fermion.md
index 1a858ab8..4e07c4fa 100644
--- a/docs/src/heom_matrix/M_Fermion.md
+++ b/docs/src/heom_matrix/M_Fermion.md
@@ -31,7 +31,7 @@ M_odd  = M_Fermion(Hs, tier, Bath, ODD)
 The fields of the structure [`M_Fermion`](@ref) are as follows:
  - `data` : the sparse matrix of HEOM Liouvillian superoperator
  - `tier` : the tier (cutoff level) for the fermionic hierarchy
- - `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+ - `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
  - `N` : the number of total [ADOs](@ref doc-ADOs)
  - `sup_dim` : the dimension of system superoperator
  - `parity` : the [parity](@ref doc-Parity) label of the operator which HEOMLS is acting on.
@@ -44,6 +44,7 @@ M::M_Fermion
 
 M.data
 M.tier
+M.dimensions
 M.dims
 M.N
 M.sup_dim
diff --git a/docs/src/heom_matrix/schrodinger_eq.md b/docs/src/heom_matrix/schrodinger_eq.md
index 0f7c01e6..63ce31f2 100644
--- a/docs/src/heom_matrix/schrodinger_eq.md
+++ b/docs/src/heom_matrix/schrodinger_eq.md
@@ -31,7 +31,7 @@ M_odd  = M_S(Hs, ODD)
 The fields of the structure [`M_S`](@ref) are as follows:
  - `data` : the sparse matrix of HEOM Liouvillian superoperator
  - `tier` : the tier (cutoff level) for the hierarchy, which equals to `0` in this case
- - `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+ - `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
  - `N` : the number of total [ADOs](@ref doc-ADOs), which equals to `1` (only the reduced density operator) in this case
  - `sup_dim` : the dimension of system superoperator
  - `parity::AbstractParity` : the [parity](@ref doc-Parity) label of the operator which HEOMLS is acting on.
@@ -42,6 +42,7 @@ M::M_S
 
 M.data
 M.tier
+M.dimensions
 M.dims
 M.N
 M.sup_dim
diff --git a/docs/src/libraryAPI.md b/docs/src/libraryAPI.md
index d0c8d7a7..035cfe12 100644
--- a/docs/src/libraryAPI.md
+++ b/docs/src/libraryAPI.md
@@ -62,7 +62,7 @@ ODD
 HEOMSuperOp
 HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix)
 HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs)
-HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int)
+HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int)
 AbstractHEOMLSMatrix
 M_S
 M_S(Hsys::QuantumObject, parity::AbstractParity=EVEN; verbose::Bool=true)
diff --git a/ext/HierarchicalEOM_CUDAExt.jl b/ext/HierarchicalEOM_CUDAExt.jl
index 41183f43..b3b8b805 100644
--- a/ext/HierarchicalEOM_CUDAExt.jl
+++ b/ext/HierarchicalEOM_CUDAExt.jl
@@ -22,17 +22,28 @@ Return a new HEOMLS-matrix-type object with `M.data` is in the type of `CuSparse
 function CuSparseMatrixCSC(M::T) where {T<:AbstractHEOMLSMatrix}
     A_gpu = _convert_to_gpu_matrix(M.data)
     if T <: M_S
-        return M_S(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity)
+        return M_S(A_gpu, M.tier, M.dimensions, M.N, M.sup_dim, M.parity)
     elseif T <: M_Boson
-        return M_Boson(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
+        return M_Boson(A_gpu, M.tier, M.dimensions, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
     elseif T <: M_Fermion
-        return M_Fermion(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
+        return M_Fermion(A_gpu, M.tier, M.dimensions, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
     else
-        return M_Boson_Fermion(A_gpu, M.Btier, M.Ftier, M.dims, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy)
+        return M_Boson_Fermion(
+            A_gpu,
+            M.Btier,
+            M.Ftier,
+            M.dimensions,
+            M.N,
+            M.sup_dim,
+            M.parity,
+            M.Bbath,
+            M.Fbath,
+            M.hierarchy,
+        )
     end
 end
 
-CuSparseMatrixCSC{ComplexF32}(M::HEOMSuperOp) = HEOMSuperOp(_convert_to_gpu_matrix(M.data), M.dims, M.N, M.parity)
+CuSparseMatrixCSC{ComplexF32}(M::HEOMSuperOp) = HEOMSuperOp(_convert_to_gpu_matrix(M.data), M.dimensions, M.N, M.parity)
 
 function _convert_to_gpu_matrix(A::AbstractSparseMatrix)
     if A isa CuSparseMatrixCSC{ComplexF32,Int32}
@@ -53,7 +64,7 @@ _convert_to_gpu_matrix(A::MatrixOperator) = MatrixOperator(_convert_to_gpu_matri
 _convert_to_gpu_matrix(A::ScaledOperator) = ScaledOperator(A.λ, _convert_to_gpu_matrix(A.L))
 _convert_to_gpu_matrix(A::AddedOperator) = AddedOperator(map(op -> _convert_to_gpu_matrix(op), A.ops))
 
-_Tr(M::Type{<:CuSparseMatrixCSC}, dims::SVector, N::Int) = CuSparseVector(_Tr(eltype(M), dims, N))
+_Tr(M::Type{<:CuSparseMatrixCSC}, dimensions::Dimensions, N::Int) = CuSparseVector(_Tr(eltype(M), dimensions, N))
 
 # change the type of `ADOs` to match the type of HEOMLS matrix
 _HandleVectorType(M::Type{<:CuSparseMatrixCSC}, V::SparseVector) = CuArray{_CType(eltype(M))}(V)
diff --git a/src/ADOs.jl b/src/ADOs.jl
index 06cae660..798baa7b 100644
--- a/src/ADOs.jl
+++ b/src/ADOs.jl
@@ -7,10 +7,13 @@ The Auxiliary Density Operators for HEOM model.
 
 # Fields
 - `data` : the vectorized auxiliary density operators
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of auxiliary density operators
 - `parity`: the parity label (`EVEN` or `ODD`).
 
+!!! note "`dims` property"
+    For a given `ados::ADOs`, `ados.dims` or `getproperty(ados, :dims)` returns its `dimensions` in the type of integer-vector.
+
 # Methods
 One can obtain the density matrix for specific index (`idx`) by calling : `ados[idx]`.
 `HierarchicalEOM.jl` also supports the following calls (methods) :
@@ -26,15 +29,17 @@ end
 """
 struct ADOs
     data::SparseVector{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     parity::AbstractParity
-end
 
-# these functions are for forward compatibility
-ADOs(data::SparseVector{ComplexF64,Int64}, dim::Int, N::Int, parity::AbstractParity) = ADOs(data, [dim], N, parity)
-ADOs(data::SparseVector{ComplexF64,Int64}, dims::AbstractVector, N::Int, parity::AbstractParity) =
-    ADOs(data, SVector{length(dims),Int}(dims), N, parity)
+    function ADOs(data::AbstractVector, dims, N::Int, parity::AbstractParity)
+        dimensions = _gen_dimensions(dims)
+        Vsize = size(data, 1)
+        ((Vsize / N) == prod(dimensions)^2) || error("The `dimensions` is not consistent with the ADOs number `N`.")
+        return new(sparsevec(data), dimensions, N, parity)
+    end
+end
 
 @doc raw"""
     ADOs(V, N, parity)
@@ -45,16 +50,7 @@ Generate the object of auxiliary density operators for HEOM model.
 - `N::Int` : the number of auxiliary density operators.
 - `parity::AbstractParity` : the parity label (`EVEN` or `ODD`). Default to `EVEN`.
 """
-function ADOs(V::AbstractVector, N::Int, parity::AbstractParity = EVEN)
-    # check the dimension of V
-    d = size(V, 1)
-    dim = √(d / N)
-    if isinteger(dim)
-        return ADOs(sparsevec(V), SVector{1,Int}(Int(dim)), N, parity)
-    else
-        error("The dimension of vector is not consistent with the ADOs number \"N\".")
-    end
-end
+ADOs(V::AbstractVector, N::Int, parity::AbstractParity = EVEN) = ADOs(V, isqrt(Int(size(V, 1) / N)), N, parity)
 
 @doc raw"""
     ADOs(ρ, N, parity)
@@ -67,11 +63,20 @@ Generate the object of auxiliary density operators for HEOM model.
 """
 function ADOs(ρ::QuantumObject, N::Int = 1, parity::AbstractParity = EVEN)
     _ρ = sparsevec(ket2dm(ρ).data)
-    return ADOs(sparsevec(_ρ.nzind, _ρ.nzval, N * length(_ρ)), ρ.dims, N, parity)
+    return ADOs(sparsevec(_ρ.nzind, _ρ.nzval, N * length(_ρ)), ρ.dimensions, N, parity)
 end
 ADOs(ρ, N::Int = 1, parity::AbstractParity = EVEN) =
     error("HierarchicalEOM doesn't support input `ρ` with type : $(typeof(ρ))")
 
+function Base.getproperty(ados::ADOs, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(ados, :dimensions))
+    else
+        return getfield(ados, key)
+    end
+end
+
 Base.checkbounds(A::ADOs, i::Int) =
     ((i > A.N) || (i < 1)) ? error("Attempt to access $(A.N)-element ADOs at index [$(i)]") : nothing
 
@@ -92,10 +97,10 @@ Base.lastindex(A::ADOs) = length(A)
 function Base.getindex(A::ADOs, i::Int)
     checkbounds(A, i)
 
-    D = prod(A.dims)
+    D = prod(A.dimensions)
     sup_dim = D^2
     back = sup_dim * i
-    return QuantumObject(reshape(A.data[(back-sup_dim+1):back], D, D), Operator, A.dims)
+    return QuantumObject(reshape(A.data[(back-sup_dim+1):back], D, D), Operator, A.dimensions)
 end
 
 function Base.getindex(A::ADOs, r::UnitRange{Int})
@@ -103,11 +108,11 @@ function Base.getindex(A::ADOs, r::UnitRange{Int})
     checkbounds(A, r[end])
 
     result = []
-    D = prod(A.dims)
+    D = prod(A.dimensions)
     sup_dim = D^2
     for i in r
         back = sup_dim * i
-        push!(result, QuantumObject(reshape(A.data[(back-sup_dim+1):back], D, D), Operator, A.dims))
+        push!(result, QuantumObject(reshape(A.data[(back-sup_dim+1):back], D, D), Operator, A.dimensions))
     end
     return result
 end
@@ -115,8 +120,10 @@ Base.getindex(A::ADOs, ::Colon) = getindex(A, 1:lastindex(A))
 
 Base.iterate(A::ADOs, state::Int = 1) = state > length(A) ? nothing : (A[state], state + 1)
 
-Base.show(io::IO, A::ADOs) =
-    print(io, "$(A.N) Auxiliary Density Operators with $(A.parity) and (system) dims = $(A.dims)\n")
+Base.show(io::IO, A::ADOs) = print(
+    io,
+    "$(A.N) Auxiliary Density Operators with $(A.parity) and (system) dims = $(_get_dims_string(A.dimensions))\n",
+)
 Base.show(io::IO, m::MIME"text/plain", A::ADOs) = show(io, A)
 
 @doc raw"""
@@ -130,8 +137,8 @@ Return the density matrix of the reduced state (system) from a given auxiliary d
 - `ρ::QuantumObject` : The density matrix of the reduced state
 """
 function getRho(ados::ADOs)
-    D = prod(ados.dims)
-    return QuantumObject(reshape(ados.data[1:(D^2)], D, D), Operator, ados.dims)
+    D = prod(ados.dimensions)
+    return QuantumObject(reshape(ados.data[1:(D^2)], D, D), Operator, ados.dimensions)
 end
 
 @doc raw"""
@@ -168,9 +175,9 @@ where ``O`` is the operator and ``\rho`` is the reduced density operator in the
 function QuantumToolbox.expect(op, ados::ADOs; take_real::Bool = true)
     if op isa HEOMSuperOp
         _check_sys_dim_and_ADOs_num(op, ados)
-        exp_val = dot(transpose(_Tr(eltype(ados), ados.dims, ados.N)), (SparseMatrixCSC(op) * ados).data)
+        exp_val = dot(transpose(_Tr(eltype(ados), ados.dimensions, ados.N)), (SparseMatrixCSC(op) * ados).data)
     else
-        _op = HandleMatrixType(op, ados.dims, "op (observable)"; type = Operator)
+        _op = HandleMatrixType(op, ados.dimensions, "op (observable)"; type = Operator)
         exp_val = tr(_op.data * getRho(ados).data)
     end
 
@@ -198,7 +205,7 @@ where ``O`` is the operator and ``\rho`` is the reduced density operator in one
 - `exp_val` : The expectation value
 """
 function QuantumToolbox.expect(op, ados_list::Vector{ADOs}; take_real::Bool = true)
-    dims = ados_list[1].dims
+    dimensions = ados_list[1].dimensions
     N = ados_list[1].N
     for i in 2:length(ados_list)
         _check_sys_dim_and_ADOs_num(ados_list[1], ados_list[i])
@@ -208,9 +215,9 @@ function QuantumToolbox.expect(op, ados_list::Vector{ADOs}; take_real::Bool = tr
         _check_sys_dim_and_ADOs_num(op, ados_list[1])
         _op = op
     else
-        _op = HEOMSuperOp(spre(op), EVEN, dims, N)
+        _op = HEOMSuperOp(spre(op), EVEN, dimensions, N)
     end
-    tr_op = transpose(_Tr(eltype(op), dims, N)) * SparseMatrixCSC(_op).data
+    tr_op = transpose(_Tr(eltype(op), dimensions, N)) * SparseMatrixCSC(_op).data
 
     exp_val = [dot(tr_op, ados.data) for ados in ados_list]
 
diff --git a/src/HeomBase.jl b/src/HeomBase.jl
index be5c98aa..c0c51f61 100644
--- a/src/HeomBase.jl
+++ b/src/HeomBase.jl
@@ -2,6 +2,15 @@ export AbstractHEOMLSMatrix
 
 abstract type AbstractHEOMLSMatrix{T} end
 
+function Base.getproperty(M::AbstractHEOMLSMatrix, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(M, :dimensions))
+    else
+        return getfield(M, key)
+    end
+end
+
 @doc raw"""
     (M::AbstractHEOMLSMatrix)(p, t)
 
@@ -32,12 +41,12 @@ _get_SciML_matrix_wrapper(M::AddedOperator) = _get_SciML_matrix_wrapper(M.ops[1]
 _get_SciML_matrix_wrapper(M::AbstractHEOMLSMatrix) = _get_SciML_matrix_wrapper(M.data)
 
 # equal to : sparse(vec(system_identity_matrix))
-function _Tr(T::Type{<:Number}, dims::SVector, N::Int)
-    D = prod(dims)
+function _Tr(T::Type{<:Number}, dimensions::Dimensions, N::Int)
+    D = prod(dimensions)
     return SparseVector(N * D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(T, D))
 end
-_Tr(M::AbstractHEOMLSMatrix) = _Tr(_get_SciML_matrix_wrapper(M), M.dims, M.N)
-_Tr(M::Type{<:SparseMatrixCSC}, dims::SVector, N::Int) = _Tr(eltype(M), dims, N)
+_Tr(M::AbstractHEOMLSMatrix) = _Tr(_get_SciML_matrix_wrapper(M), M.dimensions, M.N)
+_Tr(M::Type{<:SparseMatrixCSC}, dimensions::Dimensions, N::Int) = _Tr(eltype(M), dimensions, N)
 
 function HandleMatrixType(
     M::AbstractQuantumObject,
@@ -55,19 +64,19 @@ function HandleMatrixType(
 end
 function HandleMatrixType(
     M::AbstractQuantumObject,
-    dims::SVector,
+    dimensions::Dimensions,
     MatrixName::String = "";
     type::T = nothing,
 ) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}}
-    if M.dims == dims
+    if M.dimensions == dimensions
         return HandleMatrixType(M, MatrixName; type = type)
     else
-        error("The dims of $(MatrixName) should be: $(dims)")
+        error("The dimensions of $(MatrixName) should be: $(_get_dims_string(dimensions))")
     end
 end
 HandleMatrixType(
     M,
-    dims::SVector,
+    dimensions::Dimensions,
     MatrixName::String = "";
     type::T = nothing,
 ) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}} =
@@ -86,7 +95,7 @@ _HandleVectorType(M::Type{<:SparseMatrixCSC}, V::SparseVector) = Vector{_CType(e
 function _HandleSteadyStateMatrix(M::AbstractHEOMLSMatrix{<:MatrixOperator})
     S = size(M, 1)
     ElType = eltype(M)
-    D = prod(M.dims)
+    D = prod(M.dimensions)
     A = copy(M.data.A)
     A[1, 1:S] .= 0
 
@@ -96,8 +105,8 @@ function _HandleSteadyStateMatrix(M::AbstractHEOMLSMatrix{<:MatrixOperator})
 end
 
 function _check_sys_dim_and_ADOs_num(A, B)
-    if (A.dims != B.dims)
-        error("Inconsistent system dimension (\"dims\").")
+    if (A.dimensions != B.dimensions)
+        error("Inconsistent system dimensions.")
     end
 
     if (A.N != B.N)
diff --git a/src/HierarchicalEOM.jl b/src/HierarchicalEOM.jl
index 6c0a3368..4a4fa09a 100644
--- a/src/HierarchicalEOM.jl
+++ b/src/HierarchicalEOM.jl
@@ -5,7 +5,8 @@ import Reexport: @reexport
 @reexport using QuantumToolbox
 
 # intrinsic QuantumToolbox functions
-import QuantumToolbox: _FType, _CType, _spre, _spost, _sprepost, _liouvillian
+import QuantumToolbox:
+    _FType, _CType, _spre, _spost, _sprepost, _liouvillian, _gen_dimensions, _get_dims_string, dimensions_to_dims
 
 # SciML packages (for OrdinaryDiffEq and LinearSolve)
 import SciMLBase: init, solve, solve!, u_modified!, ODEProblem, FullSpecialize, CallbackSet, NullParameters
diff --git a/src/bath/BosonBath.jl b/src/bath/BosonBath.jl
index a41078d5..f44e819d 100644
--- a/src/bath/BosonBath.jl
+++ b/src/bath/BosonBath.jl
@@ -203,14 +203,17 @@ A bosonic bath for the real part of bath correlation function ``C^{u=\textrm{R}}
 
 # Fields
 - `Comm`  : the super-operator (commutator) for the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` in real part of bath correlation function ``C^{u=\textrm{R}}``.
 - `γ` : the coefficients ``\gamma_i`` in real part of bath correlation function ``C^{u=\textrm{R}}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::bosonReal`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct bosonReal <: AbstractBosonBath
     Comm::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     Nterm::Int
@@ -234,7 +237,7 @@ function bosonReal(op::QuantumObject, η_real::Vector{Ti}, γ_real::Vector{Tj})
         error("The length of \'η_real\' and \'γ_real\' should be the same.")
     end
     Id_cache = I(size(_op, 1))
-    return bosonReal(_spre(_op.data, Id_cache) - _spost(_op.data, Id_cache), _op.dims, η_real, γ_real, N_exp_term)
+    return bosonReal(_spre(_op.data, Id_cache) - _spost(_op.data, Id_cache), _op.dimensions, η_real, γ_real, N_exp_term)
 end
 
 @doc raw"""
@@ -244,15 +247,18 @@ A bosonic bath for the imaginary part of bath correlation function ``C^{u=\textr
 # Fields
 - `Comm`  : the super-operator (commutator) for the coupling operator.
 - `anComm`  : the super-operator (anti-commutator) for the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` in imaginary part of bath correlation function ``C^{u=\textrm{I}}``.
 - `γ` : the coefficients ``\gamma_i`` in imaginary part of bath correlation function ``C^{u=\textrm{I}}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::bosonImag`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct bosonImag <: AbstractBosonBath
     Comm::SparseMatrixCSC{ComplexF64,Int64}
     anComm::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     Nterm::Int
@@ -278,7 +284,7 @@ function bosonImag(op::QuantumObject, η_imag::Vector{Ti}, γ_imag::Vector{Tj})
     Id_cache = I(size(_op, 1))
     spreQ = _spre(_op.data, Id_cache)
     spostQ = _spost(_op.data, Id_cache)
-    return bosonImag(spreQ - spostQ, spreQ + spostQ, _op.dims, η_imag, γ_imag, N_exp_term)
+    return bosonImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_imag, γ_imag, N_exp_term)
 end
 
 @doc raw"""
@@ -288,16 +294,19 @@ A bosonic bath which the real part and imaginary part of the bath correlation fu
 # Fields
 - `Comm`  : the super-operator (commutator) for the coupling operator.
 - `anComm`  : the super-operator (anti-commutator) for the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η_real` : the real part of coefficients ``\eta_i`` in bath correlation function ``\sum_i \eta_i \exp(-\gamma_i t)``.
 - `η_imag` : the imaginary part of coefficients ``\eta_i`` in bath correlation function ``\sum_i \eta_i \exp(-\gamma_i t)``.
 - `γ` : the coefficients ``\gamma_i`` in bath correlation function ``\sum_i \eta_i \exp(-\gamma_i t)``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::bosonRealImag`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct bosonRealImag <: AbstractBosonBath
     Comm::SparseMatrixCSC{ComplexF64,Int64}
     anComm::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η_real::AbstractVector
     η_imag::AbstractVector
     γ::AbstractVector
@@ -330,7 +339,7 @@ function bosonRealImag(
     Id_cache = I(size(_op, 1))
     spreQ = _spre(_op.data, Id_cache)
     spostQ = _spost(_op.data, Id_cache)
-    return bosonRealImag(spreQ - spostQ, spreQ + spostQ, _op.dims, η_real, η_imag, γ, N_exp_term)
+    return bosonRealImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_real, η_imag, γ, N_exp_term)
 end
 
 @doc raw"""
@@ -381,17 +390,20 @@ An bath object which describes the absorption process of the bosonic system by a
 - `spre`   : the super-operator (left side operator multiplication) for the coupling operator.
 - `spost`  : the super-operator (right side operator multiplication) for the coupling operator.
 - `CommD`  : the super-operator (commutator) for the adjoint of the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `γ` : the coefficients ``\gamma_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `η_emit` : the coefficients ``\eta_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::bosonAbsorb`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct bosonAbsorb <: AbstractBosonBath
     spre::SparseMatrixCSC{ComplexF64,Int64}
     spost::SparseMatrixCSC{ComplexF64,Int64}
     CommD::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     η_emit::AbstractVector
@@ -426,7 +438,7 @@ function bosonAbsorb(
         _spre(_op.data, Id_cache),
         _spost(_op.data, Id_cache),
         _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
-        _op.dims,
+        _op.dimensions,
         η_absorb,
         γ_absorb,
         η_emit,
@@ -442,17 +454,20 @@ An bath object which describes the emission process of the bosonic system by a c
 - `spre`   : the super-operator (left side operator multiplication) for the coupling operator.
 - `spost`  : the super-operator (right side operator multiplication) for the coupling operator.
 - `CommD`  : the super-operator (commutator) for the adjoint of the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `γ` : the coefficients ``\gamma_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `η_absorb` : the coefficients ``\eta_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::bosonEmit`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct bosonEmit <: AbstractBosonBath
     spre::SparseMatrixCSC{ComplexF64,Int64}
     spost::SparseMatrixCSC{ComplexF64,Int64}
     CommD::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     η_absorb::AbstractVector
@@ -487,7 +502,7 @@ function bosonEmit(
         _spre(_op.data, Id_cache),
         _spost(_op.data, Id_cache),
         _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
-        _op.dims,
+        _op.dimensions,
         η_emit,
         γ_emit,
         η_absorb,
@@ -495,6 +510,18 @@ function bosonEmit(
     )
 end
 
+function Base.getproperty(
+    b::BType,
+    key::Symbol,
+) where {BType<:Union{bosonReal,bosonImag,bosonRealImag,bosonAbsorb,bosonEmit}}
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(b, :dimensions))
+    else
+        return getfield(b, key)
+    end
+end
+
 @doc raw"""
     correlation_function(bath, tlist)
 Calculate the correlation function ``C(t)`` for a given bosonic bath and time list.
diff --git a/src/bath/FermionBath.jl b/src/bath/FermionBath.jl
index fd2fba6c..426019b1 100644
--- a/src/bath/FermionBath.jl
+++ b/src/bath/FermionBath.jl
@@ -87,18 +87,21 @@ An bath object which describes the absorption process of the fermionic system by
 - `spost`  : the super-operator (right side operator multiplication) for the coupling operator.
 - `spreD`  : the super-operator (left side operator multiplication) for the adjoint of the coupling operator.
 - `spostD` : the super-operator (right side operator multiplication) for the adjoint of the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `γ` : the coefficients ``\gamma_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `η_emit` : the coefficients ``\eta_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::fermionAbsorb`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct fermionAbsorb <: AbstractFermionBath
     spre::SparseMatrixCSC{ComplexF64,Int64}
     spost::SparseMatrixCSC{ComplexF64,Int64}
     spreD::SparseMatrixCSC{ComplexF64,Int64}
     spostD::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     η_emit::AbstractVector
@@ -134,7 +137,7 @@ function fermionAbsorb(
         _spost(_op.data, Id_cache),
         _spre(adjoint(_op).data, Id_cache),
         _spost(adjoint(_op).data, Id_cache),
-        _op.dims,
+        _op.dimensions,
         η_absorb,
         γ_absorb,
         η_emit,
@@ -151,18 +154,21 @@ An bath object which describes the emission process of the fermionic system by a
 - `spost`  : the super-operator (right side operator multiplication) for the coupling operator.
 - `spreD`  : the super-operator (left side operator multiplication) for the adjoint of the coupling operator.
 - `spostD` : the super-operator (right side operator multiplication) for the adjoint of the coupling operator.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `η` : the coefficients ``\eta_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `γ` : the coefficients ``\gamma_i`` of emission bath correlation function ``C^{\nu=-}``.
 - `η_absorb` : the coefficients ``\eta_i`` of absorption bath correlation function ``C^{\nu=+}``.
 - `Nterm` : the number of exponential-expansion term of correlation function
+
+!!! note "`dims` property"
+    For a given `b::fermionEmit`, `b.dims` or `getproperty(b, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct fermionEmit <: AbstractFermionBath
     spre::SparseMatrixCSC{ComplexF64,Int64}
     spost::SparseMatrixCSC{ComplexF64,Int64}
     spreD::SparseMatrixCSC{ComplexF64,Int64}
     spostD::SparseMatrixCSC{ComplexF64,Int64}
-    dims::SVector
+    dimensions::Dimensions
     η::AbstractVector
     γ::AbstractVector
     η_absorb::AbstractVector
@@ -199,7 +205,7 @@ function fermionEmit(
         _spost(_op.data, Id_cache),
         _spre(adjoint(_op).data, Id_cache),
         _spost(adjoint(_op).data, Id_cache),
-        _op.dims,
+        _op.dimensions,
         η_emit,
         γ_emit,
         η_absorb,
@@ -207,6 +213,15 @@ function fermionEmit(
     )
 end
 
+function Base.getproperty(b::BType, key::Symbol) where {BType<:Union{fermionAbsorb,fermionEmit}}
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(b, :dimensions))
+    else
+        return getfield(b, key)
+    end
+end
+
 @doc raw"""
     correlation_function(bath, tlist)
 Calculate the correlation function ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` for a given fermionic bath and time list.
diff --git a/src/density_of_states.jl b/src/density_of_states.jl
index b3515e37..250f75b5 100644
--- a/src/density_of_states.jl
+++ b/src/density_of_states.jl
@@ -54,7 +54,7 @@ Calculate density of states for the fermionic system in frequency domain.
     # Handle d_op
     MType = _get_SciML_matrix_wrapper(M)
     _tr = transpose(_Tr(M))
-    Id_sys = I(prod(d_op.dims))
+    Id_sys = I(prod(d_op.dimensions))
     Id_HEOM = I(M.N)
     d_normal = HEOMSuperOp(spre(d_op, Id_sys), ODD, M; Id_cache = Id_HEOM)
     d_dagger = HEOMSuperOp(spre(d_op', Id_sys), ODD, M; Id_cache = Id_HEOM)
diff --git a/src/deprecated.jl b/src/deprecated.jl
index 0d4230d7..69f4baa1 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -26,7 +26,7 @@ _HEOMSuperOp_deprecated_method_error() = error(
     "Specifying `mul_basis = \"L\", \"R\", or \"LR\"` to `HEOMSuperOp` has been deprecated, try to construct system SuperOperator manually by using `spre`, `spost`, or `sprepost`.",
 )
 
-HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int, mul_basis::AbstractString; Id_cache = I(N)) =
+HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int, mul_basis::AbstractString; Id_cache = I(N)) =
     _HEOMSuperOp_deprecated_method_error()
 HEOMSuperOp(
     op,
@@ -34,6 +34,6 @@ HEOMSuperOp(
     refHEOMLS::AbstractHEOMLSMatrix,
     mul_basis::AbstractString;
     Id_cache = I(refHEOMLS.N),
-) = HEOMSuperOp(op, opParity, refHEOMLS.dims, refHEOMLS.N, mul_basis; Id_cache = Id_cache)
+) = HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N, mul_basis; Id_cache = Id_cache)
 HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString; Id_cache = I(refADOs.N)) =
-    HEOMSuperOp(op, opParity, refADOs.dims, refADOs.N, mul_basis; Id_cache = Id_cache)
+    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N, mul_basis; Id_cache = Id_cache)
diff --git a/src/evolution.jl b/src/evolution.jl
index 90a4a80b..e4c94e04 100644
--- a/src/evolution.jl
+++ b/src/evolution.jl
@@ -95,7 +95,7 @@ function HEOMsolve(
         expvals = Array{ComplexF64}(undef, 0, steps + 1)
         is_empty_e_ops = true
     else
-        Id_sys = I(prod(M.dims))
+        Id_sys = I(prod(M.dimensions))
         Id_HEOM = I(M.N)
         expvals = Array{ComplexF64}(undef, length(e_ops), steps + 1)
         tr_e_ops = _generate_Eops(M, e_ops, Id_sys, Id_HEOM)
@@ -124,9 +124,9 @@ function HEOMsolve(
         if !is_empty_e_ops
             _expect = op -> dot(op, ρvec)
             @. expvals[:, progr.counter[]+1] = _expect(tr_e_ops)
-            n == steps ? ADOs_list[1] = ADOs(ρvec, M.dims, M.N, M.parity) : nothing
+            n == steps ? ADOs_list[1] = ADOs(ρvec, M.dimensions, M.N, M.parity) : nothing
         else
-            ADOs_list[n+1] = ADOs(ρvec, M.dims, M.N, M.parity)
+            ADOs_list[n+1] = ADOs(ρvec, M.dimensions, M.N, M.parity)
         end
 
         ρvec = exp_Mt * ρvec
@@ -214,7 +214,7 @@ function HEOMsolve(
         tr_e_ops = typeof(M.data)[]
         is_empty_e_ops = true
     else
-        Id_sys = I(prod(M.dims))
+        Id_sys = I(prod(M.dimensions))
         Id_HEOM = I(M.N)
         expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
         tr_e_ops = _generate_Eops(M, e_ops, Id_sys, Id_HEOM)
@@ -245,7 +245,7 @@ function HEOMsolve(
         flush(stdout)
     end
     sol = solve(prob, solver)
-    ADOs_list = map(ρvec -> ADOs(Vector{ComplexF64}(ρvec), M.dims, M.N, M.parity), sol.u)
+    ADOs_list = map(ρvec -> ADOs(Vector{ComplexF64}(ρvec), M.dimensions, M.N, M.parity), sol.u)
 
     # save ADOs to file
     if filename != ""
@@ -272,8 +272,8 @@ end
 function _generate_Eops(M::AbstractHEOMLSMatrix, e_ops, Id_sys, Id_HEOM)
     MType = _get_SciML_matrix_wrapper(M)
     tr_e_ops = [
-        transpose(_Tr(M)) * MType(HEOMSuperOp(spre(op, Id_sys), EVEN, M.dims, M.N; Id_cache = Id_HEOM)).data for
-        op in e_ops
+        transpose(_Tr(M)) * MType(HEOMSuperOp(spre(op, Id_sys), EVEN, M.dimensions, M.N; Id_cache = Id_HEOM)).data
+        for op in e_ops
     ]
     return tr_e_ops
 end
@@ -301,7 +301,7 @@ _make_L(M::AbstractHEOMLSMatrix, H_t::Nothing) = M.data
 function _make_L(M::AbstractHEOMLSMatrix, H_t::QuantumObjectEvolution)
     Id_HEOM = I(M.N)
     MType = _get_SciML_matrix_wrapper(M)
-    L_t = HandleMatrixType(liouvillian(H_t), M.dims, "H_t"; type = SuperOperator)
+    L_t = HandleMatrixType(liouvillian(H_t), M.dimensions, "H_t"; type = SuperOperator)
 
     return M.data + _L_t_to_HEOMSuperOp(MType, L_t.data, Id_HEOM)
 end
diff --git a/src/heom_matrices/HierarchyDict.jl b/src/heom_matrices/HierarchyDict.jl
index 29684f37..3b1841e0 100644
--- a/src/heom_matrices/HierarchyDict.jl
+++ b/src/heom_matrices/HierarchyDict.jl
@@ -92,7 +92,7 @@ end
 @noinline function genBathHierarchy(
     B::Vector{T},
     tier::Int,
-    dims::SVector;
+    dimensions::Dimensions;
     threshold::Real = 0.0,
 ) where {T<:AbstractBath}
     Nterm = 0
@@ -101,7 +101,7 @@ end
     if T == BosonBath
         baths = AbstractBosonBath[]
         for (α, b) in enumerate(B)
-            if b.op.dims != dims
+            if b.op.dimensions != dimensions
                 error("The matrix size of the bosonic bath coupling operators are not consistent.")
             end
             push!(baths, b.bath...)
@@ -115,7 +115,7 @@ end
     elseif T == FermionBath
         baths = AbstractFermionBath[]
         for (α, b) in enumerate(B)
-            if b.op.dims != dims
+            if b.op.dimensions != dimensions
                 error("The matrix size of the fermionic bath coupling operators are not consistent.")
             end
             push!(baths, b.bath...)
@@ -177,7 +177,7 @@ end
     fB::Vector{FermionBath},
     tier_b::Int,
     tier_f::Int,
-    dims::SVector;
+    dimensions::Dimensions;
     threshold::Real = 0.0,
 )
     # deal with boson bath
@@ -185,7 +185,7 @@ end
     bosonPtr = Tuple[]
     baths_b = AbstractBosonBath[]
     for (α, b) in enumerate(bB)
-        if b.op.dims != dims
+        if b.op.dimensions != dimensions
             error("The matrix size of the bosonic bath coupling operators are not consistent.")
         end
         push!(baths_b, b.bath...)
@@ -202,7 +202,7 @@ end
     fermionPtr = Tuple[]
     baths_f = AbstractFermionBath[]
     for (α, b) in enumerate(fB)
-        if b.op.dims != dims
+        if b.op.dimensions != dimensions
             error("The matrix size of the fermionic bath coupling operators are not consistent.")
         end
         push!(baths_f, b.bath...)
diff --git a/src/heom_matrices/M_Boson.jl b/src/heom_matrices/M_Boson.jl
index 4d86d643..eb03944a 100644
--- a/src/heom_matrices/M_Boson.jl
+++ b/src/heom_matrices/M_Boson.jl
@@ -7,17 +7,20 @@ HEOM Liouvillian superoperator matrix for bosonic bath
 # Fields
 - `data<:AbstractSciMLOperator` : the matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the bosonic hierarchy
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of total ADOs
 - `sup_dim` : the dimension of system superoperator
 - `parity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
 - `bath::Vector{BosonBath}` : the vector which stores all `BosonBath` objects
 - `hierarchy::HierarchyDict`: the object which contains all dictionaries for boson-bath-ADOs hierarchy.
+
+!!! note "`dims` property"
+    For a given `M::M_Boson`, `M.dims` or `getproperty(M, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct M_Boson{T<:AbstractSciMLOperator} <: AbstractHEOMLSMatrix{T}
     data::T
     tier::Int
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     sup_dim::Int
     parity::AbstractParity
@@ -63,7 +66,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
-    sup_dim = prod(_Hsys.dims)^2
+    sup_dim = prod(_Hsys.dimensions)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
     # the Liouvillian operator for free Hamiltonian term
@@ -74,7 +77,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
         print("Checking the importance value for each ADOs...")
         flush(stdout)
     end
-    Nado, baths, hierarchy = genBathHierarchy(Bath, tier, _Hsys.dims, threshold = threshold)
+    Nado, baths, hierarchy = genBathHierarchy(Bath, tier, _Hsys.dimensions, threshold = threshold)
     idx2nvec = hierarchy.idx2nvec
     nvec2idx = hierarchy.nvec2idx
     if verbose && (threshold > 0.0)
@@ -152,7 +155,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
         println("[DONE]")
         flush(stdout)
     end
-    return M_Boson(L_he, tier, copy(_Hsys.dims), Nado, sup_dim, parity, Bath, hierarchy)
+    return M_Boson(L_he, tier, _Hsys.dimensions, Nado, sup_dim, parity, Bath, hierarchy)
 end
 
 _getBtier(M::M_Boson) = M.tier
diff --git a/src/heom_matrices/M_Boson_Fermion.jl b/src/heom_matrices/M_Boson_Fermion.jl
index 67cea6ee..b88a1771 100644
--- a/src/heom_matrices/M_Boson_Fermion.jl
+++ b/src/heom_matrices/M_Boson_Fermion.jl
@@ -8,19 +8,22 @@ HEOM Liouvillian superoperator matrix for mixtured (bosonic and fermionic) bath
 - `data<:AbstractSciMLOperator` : the matrix of HEOM Liouvillian superoperator
 - `Btier` : the tier (cutoff level) for bosonic hierarchy
 - `Ftier` : the tier (cutoff level) for fermionic hierarchy
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of total ADOs
 - `sup_dim` : the dimension of system superoperator
 - `parity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
 - `Bbath::Vector{BosonBath}` : the vector which stores all `BosonBath` objects
 - `Fbath::Vector{FermionBath}` : the vector which stores all `FermionBath` objects
 - `hierarchy::MixHierarchyDict`: the object which contains all dictionaries for mixed-bath-ADOs hierarchy.
+
+!!! note "`dims` property"
+    For a given `M::M_Boson_Fermion`, `M.dims` or `getproperty(M, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct M_Boson_Fermion{T<:AbstractSciMLOperator} <: AbstractHEOMLSMatrix{T}
     data::T
     Btier::Int
     Ftier::Int
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     sup_dim::Int
     parity::AbstractParity
@@ -100,7 +103,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
-    sup_dim = prod(_Hsys.dims)^2
+    sup_dim = prod(_Hsys.dimensions)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
     # the Liouvillian operator for free Hamiltonian term
@@ -111,7 +114,8 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
         print("Checking the importance value for each ADOs...")
         flush(stdout)
     end
-    Nado, baths_b, baths_f, hierarchy = genBathHierarchy(Bbath, Fbath, Btier, Ftier, _Hsys.dims, threshold = threshold)
+    Nado, baths_b, baths_f, hierarchy =
+        genBathHierarchy(Bbath, Fbath, Btier, Ftier, _Hsys.dimensions, threshold = threshold)
     idx2nvec = hierarchy.idx2nvec
     nvec2idx = hierarchy.nvec2idx
     if verbose && (threshold > 0.0)
@@ -221,7 +225,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
         println("[DONE]")
         flush(stdout)
     end
-    return M_Boson_Fermion(L_he, Btier, Ftier, copy(_Hsys.dims), Nado, sup_dim, parity, Bbath, Fbath, hierarchy)
+    return M_Boson_Fermion(L_he, Btier, Ftier, _Hsys.dimensions, Nado, sup_dim, parity, Bbath, Fbath, hierarchy)
 end
 
 _getBtier(M::M_Boson_Fermion) = M.Btier
diff --git a/src/heom_matrices/M_Fermion.jl b/src/heom_matrices/M_Fermion.jl
index 9bbe8be6..b4eefd2c 100644
--- a/src/heom_matrices/M_Fermion.jl
+++ b/src/heom_matrices/M_Fermion.jl
@@ -7,17 +7,20 @@ HEOM Liouvillian superoperator matrix for fermionic bath
 # Fields
 - `data<:AbstractSciMLOperator` : the matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the fermionic hierarchy
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of total ADOs
 - `sup_dim` : the dimension of system superoperator
 - `parity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
 - `bath::Vector{FermionBath}` : the vector which stores all `FermionBath` objects
 - `hierarchy::HierarchyDict`: the object which contains all dictionaries for fermion-bath-ADOs hierarchy.
+
+!!! note "`dims` property"
+    For a given `M::M_Fermion`, `M.dims` or `getproperty(M, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct M_Fermion{T<:AbstractSciMLOperator} <: AbstractHEOMLSMatrix{T}
     data::T
     tier::Int
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     sup_dim::Int
     parity::AbstractParity
@@ -61,7 +64,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
-    sup_dim = prod(_Hsys.dims)^2
+    sup_dim = prod(_Hsys.dimensions)^2
     I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)
 
     # the Liouvillian operator for free Hamiltonian term
@@ -72,7 +75,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
         print("Checking the importance value for each ADOs...")
         flush(stdout)
     end
-    Nado, baths, hierarchy = genBathHierarchy(Bath, tier, _Hsys.dims, threshold = threshold)
+    Nado, baths, hierarchy = genBathHierarchy(Bath, tier, _Hsys.dimensions, threshold = threshold)
     idx2nvec = hierarchy.idx2nvec
     nvec2idx = hierarchy.nvec2idx
     if verbose && (threshold > 0.0)
@@ -149,7 +152,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
         println("[DONE]")
         flush(stdout)
     end
-    return M_Fermion(L_he, tier, copy(_Hsys.dims), Nado, sup_dim, parity, Bath, hierarchy)
+    return M_Fermion(L_he, tier, _Hsys.dimensions, Nado, sup_dim, parity, Bath, hierarchy)
 end
 
 _getBtier(M::M_Fermion) = 0
diff --git a/src/heom_matrices/M_S.jl b/src/heom_matrices/M_S.jl
index aab502b7..bf99ffdb 100644
--- a/src/heom_matrices/M_S.jl
+++ b/src/heom_matrices/M_S.jl
@@ -12,15 +12,18 @@ where ``[\cdot, \cdot]_-`` stands for commutator.
 # Fields
 - `data<:AbstractSciMLOperator` : the matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the hierarchy, which equals to `0` in this case
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of total ADOs, which equals to `1` (only the reduced density operator) in this case
 - `sup_dim` : the dimension of system superoperator
 - `parity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
+
+!!! note "`dims` property"
+    For a given `M::M_S`, `M.dims` or `getproperty(M, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct M_S{T<:AbstractSciMLOperator} <: AbstractHEOMLSMatrix{T}
     data::T
     tier::Int
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     sup_dim::Int
     parity::AbstractParity
@@ -46,7 +49,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
-    sup_dim = prod(_Hsys.dims)^2
+    sup_dim = prod(_Hsys.dimensions)^2
 
     # the Liouvillian operator for free Hamiltonian
     if verbose
@@ -58,7 +61,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
         println("[DONE]")
         flush(stdout)
     end
-    return M_S(Lsys, 0, copy(_Hsys.dims), 1, sup_dim, parity)
+    return M_S(Lsys, 0, _Hsys.dimensions, 1, sup_dim, parity)
 end
 
 _getBtier(M::M_S) = 0
diff --git a/src/heom_matrices/heom_matrix_base.jl b/src/heom_matrices/heom_matrix_base.jl
index 2285441d..cab51ecd 100644
--- a/src/heom_matrices/heom_matrix_base.jl
+++ b/src/heom_matrices/heom_matrix_base.jl
@@ -8,17 +8,29 @@ General HEOM superoperator matrix.
 
 # Fields
 - `data<:Union{AbstractSparseMatrix}` : the HEOM superoperator matrix
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dimensions` : the dimension list of the coupling operator (should be equal to the system dimensions).
 - `N` : the number of auxiliary density operators
 - `parity`: the parity label (`EVEN` or `ODD`).
+
+!!! note "`dims` property"
+    For a given `M::HEOMSuperOp`, `M.dims` or `getproperty(M, :dims)` returns its `dimensions` in the type of integer-vector.
 """
 struct HEOMSuperOp{DT<:AbstractSparseMatrix}
     data::DT
-    dims::SVector
+    dimensions::Dimensions
     N::Int
     parity::AbstractParity
 end
 
+function Base.getproperty(M::HEOMSuperOp, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(M, :dimensions))
+    else
+        return getfield(M, key)
+    end
+end
+
 @doc raw"""
     HEOMSuperOp(op, opParity, refHEOMLS; Id_cache=I(refHEOMLS.N))
 Construct the HEOM superoperator matrix corresponding to the given system SuperOperator which acts on all `ADOs`.  
@@ -28,10 +40,10 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 # Parameters
 - `op` : The system SuperOperator which will act on all `ADOs`.
 - `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`.
-- `refHEOMLS::AbstractHEOMLSMatrix` : copy the system `dims` and number of `ADOs` (`N`) from this reference HEOMLS matrix
+- `refHEOMLS::AbstractHEOMLSMatrix` : copy the system `dimensions` and number of `ADOs` (`N`) from this reference HEOMLS matrix
 """
 HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix; Id_cache = I(refHEOMLS.N)) =
-    HEOMSuperOp(op, opParity, refHEOMLS.dims, refHEOMLS.N; Id_cache = Id_cache)
+    HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N; Id_cache = Id_cache)
 
 @doc raw"""
     HEOMSuperOp(op, opParity, refADOs; Id_cache=I(refADOs.N))
@@ -42,10 +54,10 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 # Parameters
 - `op` : The system SuperOperator which will act on all `ADOs`.
 - `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`.
-- `refADOs::ADOs` : copy the system `dims` and number of `ADOs` (`N`) from this reference `ADOs`   
+- `refADOs::ADOs` : copy the system `dimensions` and number of `ADOs` (`N`) from this reference `ADOs`   
 """
 HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs; Id_cache = I(refADOs.N)) =
-    HEOMSuperOp(op, opParity, refADOs.dims, refADOs.N; Id_cache = Id_cache)
+    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N; Id_cache = Id_cache)
 
 @doc raw"""
     HEOMSuperOp(op, opParity, dims, N; Id_cache=I(N))
@@ -56,27 +68,22 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 # Parameters
 - `op` : The system SuperOperator which will act on all `ADOs`.
 - `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`.
-- `dims::SVector` : the dimension list of the coupling operator (should be equal to the system dims).
+- `dims` : the dimension list of the coupling operator (should be equal to the system `dimensions`).
 - `N::Int` : the number of `ADOs`.
 """
-function HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int; Id_cache = I(N))
-    sup_op = HandleMatrixType(op, dims, "op (operator)"; type = SuperOperator)
+function HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int; Id_cache = I(N))
+    dimensions = _gen_dimensions(dims)
+    sup_op = HandleMatrixType(op, dimensions, "op (operator)"; type = SuperOperator)
 
-    return HEOMSuperOp(kron(Id_cache, sup_op.data), dims, N, opParity)
+    return HEOMSuperOp(kron(Id_cache, sup_op.data), dimensions, N, opParity)
 end
-HEOMSuperOp(op, opParity::AbstractParity, dims::Int, N::Int; Id_cache = I(N)) =
-    HEOMSuperOp(op, opParity, SVector{1,Int}(dims), N; Id_cache = Id_cache)
-HEOMSuperOp(op, opParity::AbstractParity, dims::Vector{Int}, N::Int; Id_cache = I(N)) =
-    HEOMSuperOp(op, opParity, SVector{length(dims),Int}(dims), N; Id_cache = Id_cache)
-HEOMSuperOp(op, opParity::AbstractParity, dims::Tuple, N::Int; Id_cache = I(N)) =
-    HEOMSuperOp(op, opParity, SVector(dims), N; Id_cache = Id_cache)
 
 function SparseArrays.SparseMatrixCSC{T}(M::HEOMSuperOp) where {T}
     A = M.data
     if typeof(A) == SparseMatrixCSC{T}
         return M
     else
-        return HEOMSuperOp(SparseMatrixCSC{T}(M.data), M.dims, M.N, M.parity)
+        return HEOMSuperOp(SparseMatrixCSC{T}(M.data), M.dimensions, M.N, M.parity)
     end
 end
 SparseArrays.SparseMatrixCSC(M::HEOMSuperOp) = SparseMatrixCSC{ComplexF64}(M)
@@ -124,7 +131,7 @@ function Base.show(io::IO, M::HEOMSuperOp)
     print(
         io,
         "$(M.parity) HEOM superoperator matrix acting on arbitrary-parity-ADOs\n",
-        "system dims = $(M.dims)\n",
+        "system dims = $(_get_dims_string(M.dimensions))\n",
         "number of ADOs N = $(M.N)\n",
         "data =\n",
     )
@@ -147,7 +154,7 @@ function Base.show(io::IO, M::AbstractHEOMLSMatrix)
         io,
         type,
         " type HEOMLS matrix acting on $(M.parity) ADOs\n",
-        "system dims = $(M.dims)\n",
+        "system dims = $(_get_dims_string(M.dimensions))\n",
         "number of ADOs N = $(M.N)\n",
         "data =\n",
     )
@@ -160,23 +167,23 @@ Base.show(io::IO, m::MIME"text/plain", M::AbstractHEOMLSMatrix) = show(io, M)
 function Base.:(*)(Sup::HEOMSuperOp, ados::ADOs)
     _check_sys_dim_and_ADOs_num(Sup, ados)
 
-    return ADOs(Sup.data * ados.data, ados.dims, ados.N, Sup.parity * ados.parity)
+    return ADOs(Sup.data * ados.data, ados.dimensions, ados.N, Sup.parity * ados.parity)
 end
 
 function Base.:(*)(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp)
     _check_sys_dim_and_ADOs_num(Sup1, Sup2)
 
-    return HEOMSuperOp(Sup1.data * Sup2.data, Sup1.dims, Sup1.N, Sup1.parity * Sup2.parity)
+    return HEOMSuperOp(Sup1.data * Sup2.data, Sup1.dimensions, Sup1.N, Sup1.parity * Sup2.parity)
 end
 
-Base.:(*)(n::Number, Sup::HEOMSuperOp) = HEOMSuperOp(n * Sup.data, Sup.dims, Sup.N, Sup.parity)
+Base.:(*)(n::Number, Sup::HEOMSuperOp) = HEOMSuperOp(n * Sup.data, Sup.dimensions, Sup.N, Sup.parity)
 Base.:(*)(Sup::HEOMSuperOp, n::Number) = n * Sup
 
 function Base.:(+)(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp)
     _check_sys_dim_and_ADOs_num(Sup1, Sup2)
     _check_parity(Sup1, Sup2)
 
-    return HEOMSuperOp(Sup1.data + Sup2.data, Sup1.dims, Sup1.N, Sup1.parity)
+    return HEOMSuperOp(Sup1.data + Sup2.data, Sup1.dimensions, Sup1.N, Sup1.parity)
 end
 
 Base.:(+)(M::MatrixOperator, Sup::HEOMSuperOp) = MatrixOperator(M.A + Sup.data)
@@ -193,7 +200,7 @@ function Base.:(-)(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp)
     _check_sys_dim_and_ADOs_num(Sup1, Sup2)
     _check_parity(Sup1, Sup2)
 
-    return HEOMSuperOp(Sup1.data - Sup2.data, Sup1.dims, Sup1.N, Sup1.parity)
+    return HEOMSuperOp(Sup1.data - Sup2.data, Sup1.dimensions, Sup1.N, Sup1.parity)
 end
 
 Base.:(-)(M::MatrixOperator, Sup::HEOMSuperOp) = MatrixOperator(M.A - Sup.data)
@@ -233,17 +240,17 @@ end
 
 function _reset_HEOMLS_data(M::T, new_data::AbstractSciMLOperator) where {T<:AbstractHEOMLSMatrix}
     if T <: M_S
-        return M_S(new_data, M.tier, M.dims, M.N, M.sup_dim, M.parity)
+        return M_S(new_data, M.tier, M.dimensions, M.N, M.sup_dim, M.parity)
     elseif T <: M_Boson
-        return M_Boson(new_data, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
+        return M_Boson(new_data, M.tier, M.dimensions, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
     elseif T <: M_Fermion
-        return M_Fermion(new_data, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
+        return M_Fermion(new_data, M.tier, M.dimensions, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
     else
         return M_Boson_Fermion(
             new_data,
             M.Btier,
             M.Ftier,
-            M.dims,
+            M.dimensions,
             M.N,
             M.sup_dim,
             M.parity,
@@ -274,7 +281,7 @@ Note that if ``V`` is acting on fermionic systems, it should be even-parity to b
 """
 function addBosonDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector{T} = QuantumObject[]) where {T<:QuantumObject}
     if length(jumpOP) > 0
-        Id_cache = I(prod(M.dims))
+        Id_cache = I(prod(M.dimensions))
         L = mapreduce(J -> lindblad_dissipator(J, Id_cache), +, jumpOP)
 
         return M + HEOMSuperOp(L, M.parity, M)
@@ -310,9 +317,9 @@ Note that the parity of the dissipator will be determined by the parity of the g
 """
 function addFermionDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector{T} = QuantumObject[]) where {T<:QuantumObject}
     if length(jumpOP) > 0
-        Id_cache = I(prod(M.dims))
+        Id_cache = I(prod(M.dimensions))
         L_data = mapreduce(J -> _fermion_lindblad_dissipator(J, M.parity, Id_cache), +, jumpOP)
-        L = QuantumObject(L_data, type = SuperOperator, dims = M.dims)
+        L = QuantumObject(L_data, type = SuperOperator, dims = M.dimensions)
 
         return M + HEOMSuperOp(L, M.parity, M)
     else
@@ -322,10 +329,10 @@ end
 addFermionDissipator(M::AbstractHEOMLSMatrix, jumpOP::QuantumObject) = addFermionDissipator(M, [jumpOP])
 
 function _fermion_lindblad_dissipator(
-    J::QuantumObject{DT,OperatorQuantumObject},
+    J::QuantumObject{OperatorQuantumObject},
     parity::AbstractParity,
     Id_cache = I(size(J, 1)),
-) where {DT}
+)
     _J = J.data
     Jd_J = _J' * _J
     return (-1)^(value(parity)) * _sprepost(_J, _J') - (_spre(Jd_J, Id_cache) + _spost(Jd_J, Id_cache)) / 2
@@ -359,8 +366,8 @@ function addTerminator(M::Mtype, Bath::Union{BosonBath,FermionBath}) where {Mtyp
         error("The type of input HEOMLS matrix does not support this functionality.")
     end
 
-    if M.dims != Bath.op.dims
-        error("The system dims between the HEOMLS matrix and Bath coupling operator are not consistent.")
+    if M.dimensions != Bath.op.dimensions
+        error("The system dimensions between the HEOMLS matrix and Bath coupling operator are not consistent.")
     end
 
     if Bath.δ == 0
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 4301adf5..74dfee63 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -40,7 +40,7 @@ function QuantumToolbox.steadystate(
         flush(stdout)
     end
 
-    return ADOs(Vector{ComplexF64}(sol.u), M.dims, M.N, M.parity)
+    return ADOs(Vector{ComplexF64}(sol.u), M.dimensions, M.N, M.parity)
 end
 
 @doc raw"""
@@ -100,7 +100,7 @@ function QuantumToolbox.steadystate(
         flush(stdout)
     end
 
-    return ADOs(Vector{ComplexF64}(sol.u[end]), M.dims, M.N, M.parity)
+    return ADOs(Vector{ComplexF64}(sol.u[end]), M.dimensions, M.N, M.parity)
 end
 
 function _ss_condition(integrator, abstol, reltol, min_t)